Apart from the aforementioned research on neural codes, in the realm of neural signal processing, it becomes natural to link the brain signals arising from EEG and MEG to an underlying oscillatory process which connects to the mechanistic underpinnings of brain circuitry. Utilizing these ideas, a large body of literature exists in the domain of EEG related applications like Brain Computed Interfaces (BCIs). Often in these studies motor imagery EEG signals are recorded, classified and the results of classification are used to drive a machine like the wheelchair (Leeb et al., 2007a,b). The dependence on the stationarity of signals is very important for current methods, including optimal spatial filtering (Ramoser et al., 2000) to solve these class of problems posing difficulty in reliable processing of EEG. The stochastic and non-linear nature of EEG signal thus poses critical challenges in its processing such as feature extraction and further classification (Pfurtscheller and Neuper, 2001). As of now, there exists no benchmark method to decipher this problem of EEG processing. We believe that a better understanding of the oscillatory neural network models, mimicking the underlying neural process, could pave way to a novel class of algorithms for processing EEG signals.

Although the objective of the proposed model is to shed light on the oscillatory neural code, we would also like to briefly cite literature on time series data mining and time series representations. Time series data mining is apparently a challenging one because of the unique characteristic features of the time series data such as presence of noise, and nonlinear relation of the data elements (Wilson, 2017). A problem that often arises in time series data processing is to form an optimal representation of the data either by reducing or approximating it, but making sure that the approximated version of the data still carries the local/global features of the original version. For the ease and efficient use of the data, the main challenge is to choose an optimal representation of the same. Time series data representation is a well-studied area where methods such as Discrete Fourier Transform (DFT) (Faloutsos et al., 1994), Discrete Wavelet Transform (DWT) (Percival and Walden, 2006), time series Piecewise approximation (Keogh et al., 2001a,b) have been proposed. Due to the current trends in the use of "big data" processing, other novel methods such as transformation of the time series data to discrete variables or symbols has become popular (Lin et al., 2007). The main idea behind this type of methodology is to transform time series data to a sequential data of symbols by initially discretizing the time series using methods like Piecewise Aggregate Approximation (PAA) (Keogh et al., 2001b). This can be treated as a way to reduce the number of points in the time series data and this is followed by converting the approximated numerical data to corresponding symbols using popular algorithms like SAX (Lin et al., 2007). The advantage of converting the time series to symbolic sequences is that, once the transformation is made, standard pattern matching algorithms can be applied to the sequences for further processing. The aforementioned methods are successful in data mining area, but carry little information on the neural processing of time series data. This is not a flaw of the aforementioned methods because they are not intended to provide any neural perspective on time series data processing. However, the real brain is adept at time series processing since most of the sensory inputs coming from different sensory modalities such as vision, proprioception, auditory, vestibular, tactile, and olfactory stimulus are dynamic in nature. Hence, the objective of this study is to propose a computational model that implements the autoencoding of time series data using biologically plausible neural principles. The very next sub section named as "background" gives the impetus behind the proposed modeling architecture.

#### Background

In response to the aforementioned general motivation, we now present a network of neural oscillators that serves as an oscillatory autoencoder. The reason why we choose the autoencoder architecture is due to the function it serves i.e., encoding the high dimensional input to a low dimensional abstract representation and further decoding it back to the original input signal. From a neural perspective this can be broadly viewed as different stages of neural information transfer. The first stage starts with the encoding of high dimensional sensory stimulus coming from multiple sensory modalities to a more compatible abstract representation in the subcortical structures. For example, visual information fetched by ∼125 million retinal photoreceptors converge to ∼1 million neurons of the lateral geniculate nucleus in the thalamus (Hubel, 1995). This is one of the instances (among many) of huge dimensionality reduction that takes place in the real brain. The decoder can be viewed as the stage in which the information is transferred from the subcortical structures to other cortical structures with more number neurons i.e., transfer of information from lower dimension to higher dimension (Guillery and Sherman, 2002). Standard autoencoder networks use static neurons that have limitations in capturing the temporal features of the input in a naturalistic fashion. The proposed model uses the dynamics of oscillatory system such as phase synchronization, frequency tuning, and also uses the signal processing concepts such as frequency modulation (FM) and multiplexing (MUX) to shed light on the possible information transfer mechanisms in the brain.

We brief out here the methods that are adopted to accomplish the aforementioned objective (this is explained in detail in the following methods section). In this model, a set of bandlimited signals are frequency modulated by a layer of neural oscillators, multiplexed by a layer of rate-coded neurons, and subsequently demultiplexed and demodulated by oscillatory neurons. The network is a hybrid model consisting of two kinds of oscillator models (Kuramoto and Hopf oscillators) and ratecoded dynamic neurons. The signals obtained at the output of the MUX stage may be considered as a reduced-dimensional representation of the input signals. Finally we test the model on actual EEG signals (real world data). The paper is outlined as follows. Section II presents the methods and the model equations, followed by the results in Section III and finally the discussion in Section IV.

#### METHODS

Here we propose the architecture of an autoencoder using oscillatory neurons. The motivation for an oscillatory autoencoder is explained above in the introduction section. The model architecture described here consists of Encoder and Decoder modules as shown in **Figure 1.** The encoder process the input signals and makes a lower dimensional representation of the same. The decoder module reconstructs back the original input signal from this abstract representation.

The encoder receives inputs as an array of N band limited signals, s1(t),.., sN(t). These signals are frequency modulated and multiplexed by the encoder. The multiplexed signals are demultiplexed and demodulated by the decoder. Both the encoder and the decoder are networks of oscillators. The networks are hybrids of Hopf and Kuramoto oscillators (Kuramoto, 1984; Righetti et al., 2006). The motivation for choosing two different phase oscillators is described in the decoder section. The encoder and decoder modules are modeled as follows.

# A. Encoder

The encoder has two stages viz. Frequency Modulation (FM) stage and MUX stage.

#### FM Stage

FM stage has N phase oscillators each with different intrinsic frequencies. N is equal to the dimension of the input. Each of the input signals is connected to one of these oscillators. Input is encoded by the phase dynamics as given in (1). This phase dynamics is equivalent to FM (Haykin et al., 1989) and hence the name FM stage.

$$\stackrel{\bullet}{\theta\_i} = \omega\_i^E + s\_i(t) \tag{1}$$

θi is the phase of the ith oscillator in the encoder layer. ω E i is the intrinsic angular frequency of the ith oscillator in the encoder layer. (Note: The superscript E stands for Encoder layer).

#### MUX Stage

A classical MUX in electronics literature ensures harmonious transfer of information between the sender and receiver by acting like a multiple switch (Omotayo, 1985). Hence, a MUX usually has n number of input lines and 1 output line. However, in the proposed model we do not use this strict definition of MUX instead we take the idea of compressing the n input signals to m dimensions where m<n. This is what is exactly achieved through the hidden layers of a traditional autoencoder. The reason why we named it MUX is to bring about a direct comparison of neural information transfer to the radio FM communication principles.

The MUX stage is implemented by a neural network architecture known as Lateral Anti-Hebbian Network (LAHN) that has Hebbian (excitatory) afferent and anti-Hebbian (inhibitory) lateral connections (Földiak, 1990). The dynamics of a neuron in LAHN is given by Equation (2). Hebbian learning applied to the afferent weight connections (Equation 5) brings the afferent weight vector close to the input data ensuring

feature selection by that particular neuron. The anti-Hebbian rule applied to the lateral connections induces competition among the LAHN neurons. Hence each LAHN neuron learns different features from the input data. This network was shown to extract optimal features from the input data by converging transformation weight vectors to the subspace of the principal components of the input data (Földiak, 1990). Since this network maximizes the variance of the output (Földiak, 1990), it extracts optimal features from the input data. The information required for the unsupervised learning of LAHN neuron is available locally at its synaptic connections (Equations 4, 5) and this makes the network biologically plausible.

This LAHN layer acts as the hidden layer for the oscillatory autoencoder. The low-dimensional representations constructed by the hidden layer of a traditional autoencoder are constructed by this MUX stage in the proposed model. Hence, the number of inputs going to the MUX layer is same as the number of encoder oscillators in the FM stage and the number of outputs from the MUX should be essentially lesser in number than its input to achieve a dimensionality reduction.

The dynamics of a neuron in the MUX stage is given in (2) and (3).

$$Y\_i(t) = \sum\_{j=1}^{N} q\_{ij} O\_j(t) + \sum\_{k=1}^{n} \omega\_{ik} Y\_k(t-1) \tag{2}$$

$$O\_{\dot{\jmath}} = \stackrel{\prime}{\sin} (\theta\_{\dot{\jmath}}) \tag{3}$$

Yi is the output of ith neuron, q and w are the afferent and lateral weight connections of MUX respectively, N is the dimension of the input, n is the total number of neurons in the LAHN. Oj is the state of the jth input oscillator. In MUX, lateral weights are updated using anti-Hebbian learning and afferent weights are updated using Hebbian learning (Földiak, 1990) as given in (4) and (5).

$$
\Delta \boldsymbol{w}\_{ik} = -\eta\_L \boldsymbol{Y}\_i(t) \boldsymbol{Y}\_k(t-1) \tag{4}
$$

$$
\Delta q\_{\rm ij} = \eta\_F [O\_\flat(t) Y\_i(t) - q\_{\rm ij} Y\_i^2(t)] \tag{5}
$$

ηL and ηF are the learning rates for lateral and feed forward weights respectively. MUX with n nodes trained using (4) and (5) mixes the input FM signals with a minimal overlap in their frequency spectrums which further decreases the reconstruction error.

#### B. Decoder

The decoder has three stages such as Frequency Tracking (**FT**) performed by adaptive Hopf oscillators, Demodulation (**DM**) using Kuramoto oscillators, and final smoothening of signal by low-pass filtering (**LPF**) using leaky integrator neurons stages respectively. Each section is explained in detail below.

#### FT Stage

Initially the responses of MUX are passed onwards to the FT stage. The purpose of this stage is to tease out the individual frequencies which are mixed by the MUX stage. This frequency tracking is achieved by using Hopf oscillators with adaptive frequency dynamics. Hopf oscillators were successfully implemented as an adaptive frequency system that updates its intrinsic frequency in an iterative way until it converges to one of the frequencies of the input data (Righetti et al., 2006).This system of Hopf oscillators was previously shown by Righetti et al. (2006) to learn the frequency components of its input signals. This was achieved by adding a frequency adaptation variable to the classical two variable Hopf oscillator dynamics (Righetti et al., 2006).This was shown in phase oscillators having unit circle phase plane limit cycles i.e., using Hopf oscillators. They have further explained similar frequency adaptation dynamics for relaxation oscillators too. However, in this model we are using harmonic phase oscillators for the frequency tracking stage as explained below.

Here, we wanted to achieve the aforementioned phenomena of tracking the frequency of input data. The adaptive frequency Hopf oscillators act like band-pass filters and filter out different frequency bands from the mixed input signal. The adaptive frequency dynamics is accomplished using the following equations:

$$\stackrel{\bullet}{r\_i} = r\_i(\mu - r\_i^2) \tag{6}$$

$$\stackrel{\bullet}{\phi}\_i = \omega\_i^D - \frac{\varepsilon}{r\_i} Y \sin(\phi\_i) \tag{7}$$

$$
\stackrel{\bullet}{\phi}\_i^D = -\varepsilon Y \sin(\phi\_i) \tag{8}
$$

r, ϕ and ωD are the radius, phase and angular frequency variables of a Hopf oscillator respectively (Note: the superscript D stands for Decoder module). µ is the parameter that controls the radius of the limit cycle. For µ =1, it produces a unit circle limit cycle. ε is the coupling factor between the MUX and the Hopf oscillators (Righetti et al., 2006). Because of linearity of the MUX, ε can be computed directly using (9).

$$
\kappa \, = P^+ \tag{9}
$$

$$P = (I - W)^{-1}Q\tag{10}$$

P is the transformation matrix of the MUX and P+ is the pseudo inverse of matrix P. I is the identity matrix. W and Q are the lateral and afferent weight matrices of LAHN respectively. P can be derived by virtue of the linearity of LAHN as given in (2).

#### DM Stage

The purpose of the DM stage is to extract the low-frequency, band limited message signals from the outputs of the FT stage. The DM consists of a layer of Kuramoto oscillators. This shift from Hopf oscillator (in FT stage) to Kuramoto oscillator (in DM stage) is to implement the process of phase synchronization. Kuramoto oscillatory dynamics have been previously implemented to achieve phase synchronization (Kuramoto, 1984). This synchronization in the phase of two oscillators is essential for extracting the message from the FM signal (Haykin et al., 1989) (see Supplementary Material). Each Hopf oscillator in the FT stage is coupled in a one-to-one fashion to a Kuramoto oscillator in the DM stage. The pairs of oscillators (the Kuramoto oscillators of DM and the Hopf oscillators of FT stage) are coupled through their respective phase variables as shown in (11) and (12).

$$\stackrel{\bullet}{\mathcal{V}}\_i = \omega\_i^E + KD\_i \tag{11}$$

$$D\_i = \sin(\phi\_i - \varphi\_i) \tag{12}$$

γ i is the phase variable of ith Kuramoto oscillator. It has the same intrinsic frequency, ω E i , as that of the encoder oscillators (Equation 1) and K is a positive coupling factor (Kuramoto, 1984). This stage is crucial since phase synchronization occurs at this stage and the synchronization dynamics further decodes the low frequency message signal embedded in the output of the Hopf Oscillator (see Supplementary Material).

#### LPF Stage

Di shown in (12) is the output of the decoder which is further passed through a leaky integrator to smoothen the outputs, i.e., low pass filtering (LPF stage). Leaky integrator acts as a low pass filter which further smoothens out the decoded signal, and eliminates any high frequency components present. Dynamics of leaky integrator is given in (13).

$$\frac{d\stackrel{\triangle}{s}\_i}{dt} = -A\stackrel{\triangle}{s\_i} + D\_i(t) \tag{13}$$

sI is the state of ith leaky integrator which is the reconstructed version of the input signal si(t); A is the leakage factor which is a positive constant.

Hence the proposed model is a hybrid one consisting of oscillatory layers sandwiching a rate coded layer. Hopf oscillators are used in the model for FM. A layer of linear neurons with lateral connections is used for frequency multiplexing which essentially mixes the FM signals. Hopf oscillators with adaptive frequency are used to track the carrier frequencies of the FM signals. Finally, Kuramoto oscillators are used to demodulate the FM signal and extract the message signal. Parameter values used for the simulation is given in **Table 1**.

#### RESULTS

We now test the model described in the previous section on an array of synthetic signals and also on real world EEG signals.

#### A. Simulation of the Model on Synthetic Signals

The synthetic signals used for the simulation are of the general form s(t) = A1sin(ω1t) + A2sin(ω2t).

Specifically, we consider 4 signals shown in (14), (15), (16) and (17) (**Figure 2)**.

$$s\_1(t) = \sin(10\pi t) + 0.5\sin(12\pi t)\tag{14}$$

$$s\_2(t) = \sin(20\pi t) + 0.5\sin(28\pi t)\tag{15}$$

$$s\_{\mathbb{B}}(t) = \sin(50\pi t) + 0.5\sin(56\pi t) \tag{16}$$

$$s\_4(t) = \sin(70\pi t) + 0.5\sin(80\pi t)\tag{17}$$

The initial intrinsic angular frequencies of the FM oscillators are taken as ω E = [200 Hz, 350 Hz, 850 Hz, 1000 Hz]. The input signals, as given by Equations (14)–(17), are used to modulate the encoder oscillators as per Equation (1). Let the resultant frequency modulated signals be O1, O2, O3, O4 respectively as given by Equation (3). **Figures 2A–D** shows the waveforms of the input signals (for a short duration) as given by Equations. (14)–(17). **Figures 3A–D** shows the corresponding frequency spectra. All the frequency spectra are obtained using the Fourier Transformation on the input signals. **Figures 3A–D** clearly show that the input signals are modulated to the higher frequency regime corresponding to the respective carrier waves.

The modulated signals are passed through a MUX which has two neurons (n = 2). The outputs of MUX neurons (MUX composite signal) are Y1 and Y2 as per Equation (2). The spectra of Y1 and Y2 are depicted in **Figures 3E,F**. It is evident from **Figures 3E,F** that the MUX selectively picks and mixes the frequency components of the input signals in such a way that their frequency spectra have minimum overlap. In **Figure 3E**, one neuron of the MUX was more biased to frequency spectra of O1 and O4. In **Figure 3F**, the second neuron of the MUX was more biased to frequency spectra of O2 and O3. The tendency of the hidden layer in

autoencoder network. The waveforms (A–D) follow the equations (14) – (17) respectively.

a traditional autoencoder to decorrelate the input signals, is manifesting in the present context as a tendency to remix the input signals so that there is minimal overlap in the spectrum (Földiak, 1990).

The FT stage has four Hopf oscillators, which are intended to track the four modulating frequencies. Tracking the frequency is similar to tuning the intrinsic oscillations to that particular channel frequency to fetch the information passed through that respective channel. **Figure 4** depicts the frequency adaptation of Hopf oscillators at the FT stage. The intrinsic frequencies of Hopf oscillators are initialized randomly and during the course of time their frequencies get entrained to a specific modulator frequency. Through this adaptation, oscillators are able to select a specific channel of information from a mixture of MUX signals.

**Figure 5** shows the FFT of the four Hopf oscillators' responses. It is evident from the spectrum that each Hopf oscillator is able to pick individual channel that carries the message signal and hence implements the demodulation of the frequency modulated signals. This is an interesting phenomenon which is also observed in the real brain where two cortical regions get entrained to a similar LFP frequency for information transfer or feature binding (Singer and Gray, 1995; Fell and Axmacher, 2011). Synchronization phenomena also circumvent the need for any training between the cortical structures to learn the transmitted information. That is, simply by tuning to a common frequency, two neural structures can communicate over a temporary channel, without any retraining of connections. This is discussed further in detail in the discussion section.

**Figures 6A–D** shows the original and the reconstructed signals (shown for a short duration). The demodulated signals are of lower amplitude and phase shifted compared to the original signal. This can be further corrected using proper amplification and lag shift operation on the output signals. To quantify the accuracy of reconstructed signal, we compute the reconstruction error. This gives an idea on how good the system is with regard to its function as an autoencoder. It is not advisable to directly compare the input and the raw output signal because of the phase shift present in the output signal. To this end, we first corrected the phase shift in the output signal, by computing cross correlation between the input and the output signal. Next, we computed the lag corresponding to the maximum correlation value and circular shifted the output signal using the previously found maximum lag value to correct the phase shift. The percentage (%) reconstruction error is then computed as the deviation of the Pearson correlation coefficient between the input and the phase corrected output signal from unity.

$$\%error = [1 - corr(\chi, \chi)] \times 100$$

where x is the input signal and y is the amplitude and phase corrected reconstructed output signal. The percentage (%) reconstruction error with respect to the number of nodes in the LAHN network shows a decreasing trend indicating a better recovery of signal with increasing number of nodes in the LAHN layer (**Figure 6E**). The reconstruction error is computed after

lateral weights were initialized to zero. η*L* and η*F* were taken as 10−4 .

phase correction of the output signals as explained above. This result shows that choosing an optimal number of neurons in the hidden layer (based on the reconstruction error), it is possible to form a more efficient abstract representation of the input signal.

# B. Simulation of the Model on Real World Signals (EEG Signals)

This section explains the simulation results of the oscillatory autoencoder model on real world data (i.e., data obtained through empirical ways). For this, we considered empirically recorded EEG signals obtained from BCI Competition 2008 - Graz data set B (Leeb et al., 2008). The dataset essentially consists of two class motor imagery EEG signals recorded from three channels (C3, Cz, and C4) (Leeb et al., 2008) at a sampling rate of 250 Hz. **Figure 7** shows the 1 s duration EEG signals from the aforementioned channels. These EEG signals were recorded during a motor imagery task. For further information on the experimental protocol the readers may refer (Leeb et al., 2008).

These three EEG signals form the input to the model, which is further used to modulate the frequency of the phase oscillators with intrinsic frequencies (500, 600, and 750 Hz). **Figures 8A–C** shows the frequency spectrum of the frequency modulated signals (EEG FM signals). These signals were further forward passed to the LAHN layer (with two nodes) to get the low dimensional representation of the same and to perform MUX operation. **Figures 8D,E** shows the frequency spectrum of the MUX-LAHN signals. It is vivid from the figure that each LAHN neuron captures the frequency information of the EEG FM signals and hence forms a low dimensional representation of the raw EEG signals.

Composite MUX signals are further forward passed to the adaptive Hopf oscillators where each Hopf oscillator tunes its intrinsic frequency to each channel frequency. Adaptive Hopf oscillators thus separate the signals from the composite MUX signal (as shown in **Figures 9A–B**) and this is evident from the frequency spectrum of each Hopf oscillator (**Figures 9C–E**).

The adaptive Hopf oscillator outputs are further passed to the demodulator Kuramoto oscillators for phase locking and extracts out the embedded EEG signals. **Figure 10** shows the reconstructed EEG signal from three channels along with the original signal. The reconstructed EEG signals from the oscillatory autoencoder are smoother than the original EEG signals. This smoothing could be due to the large time scales that govern learning in LAHN and adaptive Hopf oscillatory stage. The subtle changes in the reconstructed signal are due to the lower dimensional representation of the LAHN hidden layer. However, the hidden layer serves as a reliable low dimensional representation of the EEG signals which is further delineated in the following discussion section.

# Comparison of the Model Result With the Benchmark Method for Dimensionality Reduction

In the case of aforementioned EEG result, apart from computing the % reconstruction error, we compare the obtained values with a benchmark dimensionality reduction and reconstruction method to check the goodness of the proposed model. To accomplish this we performed standard Principal Component Analysis (PCA) on the input EEG data to reduce its dimension and further reconstructed back the signal to compute the % reconstruction error. After computing the reconstruction error for each signal, an average reconstruction error is computed to compare it with that of the proposed oscillatory network model. The average reconstruction error of PCA is obtained as 5.68% (computed using MATLAB custom code) considering the first two principal components (because LAHN layer in the model has 2 neurons). The average reconstruction error of the oscillatory autoencoder model is obtained as 5.31%. The average reconstruction error of the proposed oscillatory network model is slightly lower than the standard PCA method by 0.31%. Apart from the decrease in the reconstruction error, the neural attributes of the proposed model and also the theory that the model embodies on the mechanisms of neural information transfer in the brain enhances the significance of the proposed model.

#### DISCUSSION

#### Summary of the Work

We propose here an oscillatory autoencoder that reconstructs the input signal using a well defined encoder and decoder

using the principles of FM, MUX, adaptive frequency dynamics, and phase synchronization. We simulated the model using synthetic (linear combination of sinusoids) and also real world EEG signals, thus showing the robustness of the model. The proposed study gives a proof of principle for the potentiality of the oscillatory neural networks in non-trivial applications where oscillations are seldom used such as autoencoder (problem addressed in this paper), feature extraction, clustering, classification etc where mostly rate coded networks are used. The criticality of oscillations in neurobiology, as mentioned below, is the motivation of this work.

FIGURE 8 | Frequency spectrum of EEG FM signals and EEG MUX signal: (A–C) show the frequency spectrum of the EEG frequency modulated signals. (D,E) show the frequency spectrum of the composite MUX signal obtained from LAHN. It is vivid from the figures that the MUX signals cover the spectral information of the EEG FM signals and form a low dimensional representation of the same.

reconstructed EEG signals of the respective channels by the oscillatory autoencoder model. The % of reconstruction error for each reconstruction is 5.8, 6.23, and 3.9% respectively. The absence of noise in the reconstructed signal and also the dimensionality reduction can influence the reconstruction error.

#### Criticality of the Oscillations

Although in computational neuroscience literature, oscillatory neurons are not as common as rate-coded or spiking neuron models, oscillations figure prominently in experimental neurobiology. There exists a large corpus of experimental literature that correlates animal behavior with the aspects of neural oscillations (Buzsáki, 2002; Lisman and Buzsáki, 2008; Adhikari et al., 2010; Fell and Axmacher, 2011). Instances can be found from experimental neurobiology wherein all the major components of neural information processing viz., communication, representation and learning are implemented by neural oscillations. Colgin et al. (2009) reported that CA1 region of hippocampus communicates with Medial Entorhinal Cortex (MEC) via fast gamma synchronization (65–140 Hz) and with CA3 region via slow gamma synchronization (25–50 Hz) (Colgin et al., 2009). That is, by changing the frequency of the signal, it is possible to select the route by which communication takes place. Spatially distributed neurons can encode for several individual features of an object by synchronizing the neural discharges of the features, a phenomenon known as feature binding (Singer and Gray, 1995). For instance, the presentation of an optimally oriented bar gives rise to synchronized spiking of neurons, which are spatially distributed, in the area 17 of the visual cortex (Gray and Singer, 1987). Synchronization in the neural discharge is mirrored in the phase of the corresponding oscillatory LFP activity too. Hence there is a high correlation between the spike timing with the phase of the LFP oscillations. In case of feature binding, synchronization may not sometimes be evident from the spiking activity of the neurons, but the LFP activity shows robust phase synchronization (Alonso and Garcia-Austt, 1987; Buzsáki et al., 1992). Thus understanding the system dynamics in terms of oscillations becomes crucial. In the perspective of learning, a volley of high frequency pre-synaptic pulses with simultaneous depolarization at the postsynaptic side leads to Long Term Potentiation (LTP) (Bliss and Lømo, 1973; Lüscher and Malenka, 2012). These high frequency spikes can be correlated with the corresponding LFP oscillations. Hence the same LTP defined in terms of the spikes can be redefined using oscillatory LFP (Chauvette et al., 2012). Oscillations also have a pivotal role in cognition in both normal and pathological conditions. For example, the disconnectivity hypothesis of schizophrenia relates the disease symptoms to the dysfunction in the communication between different brain regions (Williams and Boksa, 2010). Gamma rhythm has been reported to have a role in the information transfer between the brain regions (Gray et al., 1989). In the early onset schizophrenic patients there is a reduction in the power of the gamma oscillation in the Prefrontal Cortex (PFC) a reason accounted for impaired working memory (Haenschel et al., 2009). Longer time scale oscillations like circadian rhythms are also known to play a critical role in major psychological disorders like bipolar disorder, depression, addiction (McClung, 2007; Alloy et al., 2015). Thus, from circadian to high gamma rhythms, oscillator models can be used to describe brain dynamics over a wide range of frequencies.

## Relation of Neural Information Transfer to Radio Communication Principles

The proposed work reinforces the hypothesis of information transfer between the brain regions to FM radio principles proposed by Hoppensteadt and Izhikevich (1998), that cortical areas communicate each other by making sure that their oscillations satisfy a resonant condition (Hoppensteadt and Izhikevich, 1998). They hypothesized that cortical oscillations are frequency modulated (FM) and, when the frequencies of two cortical areas match, they communicate by phase modulation. Thus, cortical communication is proposed to operate on the lines similar to FM radio. Although their paper proposed that signals can be frequency modulated and demodulated, it does not present how these concepts could actually be exploited to perform autoencoding, i.e., the input messages getting frequency modulated, multiplexed, demultiplexed and frequency demodulated. The proposed oscillatory autoencoder model realizes this concept by invoking the adaptive frequency and phase synchronization dynamics which take care of the frequency tuning to the incoming FM signal and hence offers a neurally plausible mechanism for the signal transmission and reconstruction (autoencoding) in the brain. This is achieved by the use of Hopf and Kuramoto dynamics. Both Kuramoto and Hopf oscillators have been previously used as models of neural oscillations in many instances (Cumin and Unsworth, 2007; Righetti et al., 2009). The Kuramoto model has been used to explain neuronal synchronization in large connected networks (Cumin and Unsworth, 2007) and especially building generative models of cortical oscillations (Breakspear et al., 2010). On the other hand, adaptive Hopf oscillators have been used for the generation of rhythmic output patterns such as central pattern generators involved in locomotion (Ijspeert et al., 2005). However, we have not come across any literature exploiting the phase synchronization properties of Kuramoto oscillators and adaptive frequency aspects of Hopf oscillator to model frequency multiplexing and demultiplexing. One of the interesting achievements of the proposed model is to show that Kuramoto—Hopf oscillator combination could act as a neural phase locked loop (PLL) which can be used to decode information from a given cortical region.

# Possible Applications of the Model

Autoencoder networks are usually constructed out of rate coded neurons, though in the recent times autoencoder networks with spiking neurons have also been proposed (Burbank, 2015). In its simplest form, a rate-coded autoencoder is a feedforward network with a single hidden layer and is trained such that the target output is the same as the input; the hidden layer has fewer neurons than the input or the output layer. Then the hidden layer learns to represent the input using fewer dimensions and therefore achieves dimensionality reduction of the input space (Hinton and Salakhutdinov, 2006). A similar reduction is achieved in the proposed oscillatory model since the hidden layer, LAHN, is of lower dimension than the input layer. The connection between Hebbian learning rule and Principal Component Analysis (PCA) is not a new idea since Oja has previously shown how a linear neuron adapting its synaptic weight connections using Hebbian learning rule can converge to the first principal component of the input data (Oja, 1989). This was further extended by Sanger using an asymmetric Generalized Hebbian Algorithm (GHA) learning rule that makes the network to learn the first n principal components (Sanger, 1989) instead of just one principal component. Hebbian/anti-Hebbian network also comes under the category of subspace learning network. This type of network, reduces the input data dimension by learning the principal subspace of the input data (Földiak, 1990; Hu et al., 2015; Pehlevan et al., 2015). Other neural networks in this line are subspace network, Rubner's network (Rubner and Tavan, 1989; Rubner and Schulten, 1990) etc. Although these networks were initially modeled to explain the computations behind the processing of streaming sensory inputs, the synaptic plasticity rules based on the local activity of the neurons neuronal activity were postulated rather than derived from a cost function (Földiak, 1990). This gap was further bridged by computing the local learning rules from a principled cost function (Hu et al., 2015; Pehlevan et al., 2015). Changing the non-linearity of the neuronal activation function explained the potentiality of these networks in extracting the higher order moments of the input data and hence qualified them as the neural architectures for Independent Component Analysis (ICA) (Oja, 1997). The aforementioned studies prove the criticality of this type of network in various applications that include subspace learning, source separation problem, dimensionality reduction etc.

This dimensionality reduction has further implications especially in EEG processing. The model reconstructs the original EEG signals from their lower dimensional LAHN representations. This means that these LAHN signals can serve as the reliable representations especially for high channel EEG signals. These representations could potentially be useful in BCI related processing such as classification of EEG signals, feature clustering, movement signature detection etc. The EEG signals used for the model simulation are two class motor imagery signals which are of particular interest in BCI application. Hence the proposed model not only provides a biologically plausible explanation for the information transfer in the brain but also shows its possible potential application in BCI related EEG processing. Another important feature that makes the current model suitable for EEG processing is its ability to average out the noise present in the input signal. As shown in the results section, the input EEG signals have high frequency ripples in its original form which is further averaged out to produce a smooth reconstructed signal as the output from the model (**Figure 10**). This could possibly be due to the large temporal scale Hebbian learning that happens in the hidden LAHN layer which could thus average out the noise present in the input.

#### Future Extensions of the Proposed Model

A possible extension of the current model could be to add additional circuitry that will enable routing of the signal from the ith input channel to the jth output channel. It must be possible to choose the input/output channels to be coupled through another layer that projects to the current LAHN layer that performs multiplexing of FM signals. In such an extended model, the LAHN layer and the additional circuitry for route selection can be compared to the functions of the thalamus with respect to cortico-thalamic information processing. Hence the proposed model serves as the proof of principle for the potentiality of the oscillatory neural networks in information transfer i.e., encoding and decoding of real world signals using the principles of modulation, MUX, frequency adaptation and phase synchronization and also shows its possible potential role in EEG related applications. Another direction that the proposed model could possibly take is to pick the brain components from the EEG signals. By brain component we mean the sources inside the brain responsible for the generation of the EEG signal. Current approaches like Independent Component Analysis (ICA) require manual selection of components which has a source inside the brain for further analysis. We envisage that using a hierarchical network of LAHN, we could possibly isolate the brain components better due to its inherent ability to filter out noise. As a future work, we envisage that the current model could be possibly used (may be by invoking minor changes) to study the EEG related phenomena like mu band Event Related Desynchronization (ERD), Visual Evoked Event Related Potential (ERP) etc which can possibly shed light on the neural principles behind the occurrence of these phenomena.

#### CONCLUSION

We propose a hybrid oscillatory network model that performs the function of an autoencoder. Using this network, we are able to encode the information onto oscillations, reduce the dimensionality of information and effectively decode them using a neural phase locked loop. The model was successfully applied to both synthetic as well as real world EEG signals. Hence the proposed model shows an oscillatory neural framework in describing information transfer in the brain. By reconstructing the EEG signals from its abstract representations in the hidden layer we have shown the model's ability in better feature extraction of EEG signal which is a critical part in EEG processing. Finally, we conclude that exploring the universality of oscillator networks would open avenues for developing an entirely new class of neural network models that describe brain function in terms of oscillatory properties— amplitude, frequency, and phase. The whole motivation of this work was to show a proof of principle for the potentiality of the oscillatory networks in other domains where usually rate coded or spiking neurons were used. In the future, we plan to apply the model to a wider variety of real world time series signals.

#### DATA AVAILABILITY STATEMENT

All the simulations are done in MATLAB R2016a and the code is made available in Model DB repository. URL is http://senselab.med.yale.edu/ModelDB/showModel.cshtml? model=243595

#### AUTHOR CONTRIBUTIONS

All authors contributed equally to the work. KS performed designing, coding, analysis of the model, and manuscript preparation. VM performed analysis of the model, and

#### REFERENCES

manuscript preparation. VC performed designing the model, analysis of the model and manuscript preparation.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom. 2018.00052/full#supplementary-material

multidimensional scaling of streaming data. Neural Comput. 27, 1461–1495. doi: 10.1162/NECO\_a\_00745

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Soman, Muralidharan and Chakravarthy. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Sensorimotor Synchronization With Auditory and Visual Modalities: Behavioral and Neural Differences

Daniel C. Comstock 1 , Michael J. Hove2 and Ramesh Balasubramaniam1 \*

*1 Cognitive and Information Sciences, University of California, Merced, Merced, CA, United States, 2 Department of Psychological Science, Fitchburg State University, Fitchburg, MA, United States*

It has long been known that the auditory system is better suited to guide temporally precise behaviors like sensorimotor synchronization (SMS) than the visual system. Although this phenomenon has been studied for many years, the underlying neural and computational mechanisms remain unclear. Growing consensus suggests the existence of multiple, interacting, context-dependent systems, and that reduced precision in visuomotor timing might be due to the way experimental tasks have been conceived. Indeed, the appropriateness of the stimulus for a given task greatly influences timing performance. In this review, we examine timing differences for sensorimotor synchronization and error correction with auditory and visual sequences, to inspect the underlying neural mechanisms that contribute to modality differences in timing. The disparity between auditory and visual timing likely relates to differences in the processing specialization between auditory and visual modalities (temporal vs. spatial). We propose this difference could offer potential explanation for the differing temporal abilities between modalities. We also offer suggestions as to how these sensory systems interface with motor and timing systems.

Keywords: sensorimotor synchronization, timing, rhythm, visual perception, auditory perception

#### INTRODUCTION

Many behavioral studies have examined human timing ability in tasks of sensorimotor synchronization (SMS) where subjects synchronize their movements to an external rhythm.

Comparisons between auditory metronomes and visual flashing metronomes reveal that movement synchronization is less variable and can occur at faster rates with auditory metronomes (Chen et al., 2002; Repp, 2003; Repp and Penel, 2004; Lorås et al., 2012). However, visuo-motor synchronization greatly improves when synchronizing with a moving periodic visual metronome (Hove et al., 2010). Adding a changing velocity profile to the moving visual metronome further reduces variability in SMS tapping (Hove et al., 2013a; Iversen et al., 2015), and Gan et al. (2015) suggests that a more realistic velocity profile can bring visual SMS to be as temporally precise as auditory SMS, at moderate but not fast tempi. While most studies of SMS look at finger tapping, others have included synchronized circle drawing, gait, dancing, and eye movements in the context of modality-specific timing effects (e.g., Repp and Su, 2013).

#### Edited by:

*Daya Shankar Gupta, Camden County College, United States*

#### Reviewed by:

*Edward W. Large, University of Connecticut, United States David Ian Anderson, San Francisco State University, United States Todd Troyer, University of Texas at San Antonio, United States*

#### \*Correspondence:

*Ramesh Balasubramaniam<EMAIL_ADDRESS>
Received: *05 March 2018* Accepted: *19 June 2018* Published: *18 July 2018*

#### Citation:

*Comstock DC, Hove MJ and Balasubramaniam R (2018) Sensorimotor Synchronization With Auditory and Visual Modalities: Behavioral and Neural Differences. Front. Comput. Neurosci. 12:53. doi: 10.3389/fncom.2018.00053*

Studies on auditory and visual interference also suggest auditory timing is more prominent. When concurrent auditory metronomes and visual flashing metronomes are presented out-of-phase, the auditory sequences interfere with visuomotor timing, but not vice versa (Repp and Penel, 2002, 2004). The interference effect is considerably reduced with moving visual metronomes and is tied to training and experience as the auditory dominance is stronger in musicians and weaker in video gamers (Hove et al., 2013a). Similarly, auditory cues can improve visual temporal discrimination (Morein-Zamir et al., 2003; Parise and Spence, 2008). This effect only holds for the temporal domain however, as the visual system dominates when auditory and visual stimuli conflict in the spatial domain; spatial dominance in the visual modality is apparent in the well-known "ventriloquist effect" (Vroomen et al., 2001).

## ROLE OF ERROR CORRECTION IN TIMING

Error correction is a crucial component of any SMS task. By inducing perturbations and errors in SMS, we can gain insight into the underlying timing mechanisms. A common method to induce errors in a SMS task is to occasionally perturb an otherwise isochronous metronome (Repp, 2000, 2001a,b; Praamstra et al., 2003; Repp and Keller, 2004; Jang et al., 2016; Jantzen et al., 2018). Error correction in SMS can be broken down into two distinct mechanisms: a phasecorrection mechanism for correcting errors in relative phase, and a period-correction mechanism that corrects changes to the internal timekeeper period (Repp, 2001b; Repp and Keller, 2004). Period corrections require conscious awareness of the error as it involves a conscious updating of the internal rhythm; while a phase correction can happen even with errors too small for conscious awareness and does not involve updating the central timekeeper period and so is considered a more peripheral process than period correction (Repp, 2001b, 2005). An error corrected under the phase-correction mechanism is typically a gradual adjustment that occurs over several beats, while an error corrected under the period-correction mechanism will be evidenced by a pronounced correction, usually followed by a more gradual phase-correction-like pattern after the initial large correction (Repp, 2001b).

While error correction has been well documented in auditory SMS, relatively little work has investigated error correction in visual SMS. In a recent study comparing error correction for auditory and flashing visual sequences, we observed error corrections for perturbations in the auditory condition that were modulated by the direction of the perturbations, but no such modulation was found for perturbations in the visual condition (Comstock and Balasubramaniam, 2017a). This suggests the visual system may not engage in the same SMS timing mechanisms as the auditory system. Additional evidence for a discrepancy in error correction for auditory and visual sequences can be gleaned from the autocorrelation structure of adjacent taps: unlike auditory SMS, tapping with visual flashes does not produce a negative lag1 autocorrelation that can indicate of the presence of a robust central timekeeping and error-correction mechanism (Hove and Keller, 2010). However, visuomotor synchronization with moving and apparent-motion metronomes do produce a negative lag1 autocorrelation, suggesting that a moving visual metronome may engage error correction (Hove and Keller, 2010; Hove et al., 2010); note that negative lag1 autocorrelation does not necessarily stem from error correction and can arise from other timing factors (e.g., Wing and Kristofferson, 1973). It remains unclear if error correction will occur with perturbations in moving visual metronomes or with larger phase perturbations in a flashing visual metronome.

#### UNDERLYING PHYSIOLOGY OF THE AUDITORY AND VISUAL TIMING SYSTEM

#### Brain Networks Involved in Timing Activity

Investigating the neural underpinnings in auditory and visual timing is a massive undertaking due to the many different timing subprocesses and tasks, including: SMS, interval timing, rhythm perception, timing recall, time perception, etc.. Excellent reviews of the brain mechanisms involved in various timing activities include: a review of neural activity in music production (Zatorre et al., 2007); a review of neural activity involved in time perception (Wiener et al., 2010); and an overview of neural activation in SMS as part of a larger review of SMS (Repp and Su, 2013). This body of work consistently demonstrates that temporal processing across tasks and sensory modalities relies heavily on the motor system. This motor network includes the supplemental motor area (SMA), primary motor cortex, lateral premotor cortex, anterior cingulate, basal ganglia, and cerebellum (Repp and Su, 2013). Auditory rhythm perception activates the motor system and is closely linked to movement (Janata et al., 2012; Iversen and Balasubramaniam, 2016; Ross et al., 2016a,b). The SMA is also strongly implicated in motor timing (Coull et al., 2016; Merchant and Yarrow, 2016), and along with the pre-SMA could be a hub of motor timing (Schwartze et al., 2012). Subcortical regions are especially active during sub-second time perception (Wiener et al., 2010), subsecond interval timing (Repp and Su, 2013), and rhythm timing (Grahn and Rowe, 2009; Wiener et al., 2010; Coull et al., 2011; Teki et al., 2011; Hove et al., 2013b). There is evidence of a dorsal auditory stream connecting the auditory cortex to the motor cortex through the posterior parietal cortex that plays a role in rhythm perception (Patel and Iversen, 2014; Ross et al., 2018). Interestingly this dorsal stream is also implicated in visual and tactile rhythm perception (Araneda et al., 2017; Rauschecker, 2017), adding to the idea of a common timing system tied to the motor system. Further evidence of the common timing system is found in a study of auditory and visual synchronization that dissociated modality and tapping stability – putamen activation was highest when synchronizing to auditory beeps, moderate with a frequency-modulated siren and with a moving visual metronome, and lowest with a flashing visual metronome, closely paralleling behavioral performance (Hove et al., 2013b).

While visual SMS activates many of the same motor regions as auditory SMS (Hove et al., 2013b; Araneda et al., 2017), some activations are specific to the visual system. The visual cortex shows activity related to interval timing that follows the expected scalar property, such that size of timing errors measured in the visual cortex scale in proportion to size of the interval being timed as predicted by Weber's law (Shuler, 2016). Additionally, Zhou et al. (2014) found evidence that visual feature processing in the early visual cortex can contribute to duration perception, furthering the notion that at least some timing information is processed independently within the visual cortex. Additionally, in visual rhythm perception, the visual cortex plays a role predicting rhythmic onsets (Comstock and Balasubramaniam, 2017b, 2018). The additional activations with visual timing tasks, taken together with behavioral results, suggest the timing accuracy in visual processing may be compared to the auditory system due to the additional computational demands of processing the higher complexity of visual spatial information along with temporal information.

#### Role of Cortical Oscillations in Timing Encoding and Spreading Information Across the Brain

In addition to looking at the networks and regions involved in temporal processing, a growing body of work shows the role of cortical oscillations in encoding timing across multiple frequency bands. Cortical oscillations play a role in connecting regions across the brain, with higher frequencies utilized for localized interaction and lower frequencies for longer range interaction (Sarnthein et al., 1998; Von Stein and Sarnthein, 2000). This pattern of oscillations is used to connect and calibrate disparate timing systems in the brain (Gupta and Chen, 2016). Oscillations relating to timing appear to arise from multiple context-specifc timing systems in the brain (Wiener and Kanai, 2016). The question is then how these functionally and anatomically disparate systems integrate and interact. It appears that oscillations from different timing systems are coordinated within the striatum (Matell and Meck, 2004; Gu et al., 2015).

Beta band activity (∼20 Hz) is tied to the motor system and several studies indicate beta's role in predicting timing of auditory rhythms (Fujioka et al., 2009, 2012, 2015). Additionally, beta activity reflects top-down imposition of metrical structure on auditory rhythms (Iversen et al., 2009). Recently, beta activity has also been linked to timing predictions within the visual system in response to visual rhythms (Comstock and Balasubramaniam, 2017b).

With rhythm perception, evidence shows that internal oscillations arise to match the fundamental frequency of the rhythm, and frequency of the meter (Nozaradan et al., 2011), as well as to the frequency of imagined rhythms (Okawa et al., 2017). These findings align with the Neural Resonance Theory that posits neural rhythms synchronize to auditory rhythms, and these neural rhythms can influence attention, expectancy, and motor planning (Large and Snyder, 2009). As of yet, it is unclear if this same neural resonance to meter would arise with visual stimuli.

# Neural Underpinnings of Error Correction

The neural correlates of error correction reveal more evidence for multiple interacting and overlapping timing mechanisms. Error detection of timing perturbations in auditory SMS tasks modulates the P1, N1, and N2 auditory ERP components depending on both the size and direction of the perturbation (Praamstra et al., 2003; Jang et al., 2016). Jantzen et al. (2018) also found a theta response stemming from the Pre-SMA and anterior cingulate for error detection, an increase in theta coupling between the SMA and the motor cortex for late perturbations. In visual error detection, the visual P1 component is reduced in latency only for large late perturbations (Comstock and Balasubramaniam, 2017a). Each of these instances show cortical activation specific to a type of perturbation, although these effects are generally limited to larger perturbations.

Smaller perturbations that elicit a phase-correction response are believed to be driven primarily by subcortical mechanisms. Applying repetitive TMS to downregulate motor and premotor cortices produced no effect on phase correction (Doumas et al., 2005), whereas phase-correction was impaired by repetitive TMS to the cerebellum (Bijsterbosch et al., 2011). This fits with the suggestion that phase-correction is primarily subcortical based on evidence from how rapidly the movement trajectory changes after a perturbation (Hove et al., 2014). A possible network that exhibits the rapid timing required for the phase-correction response is a cortico-striatal circuit connecting the cerebellum to the SMA-striatal network via the thalamus (Kotz et al., 2016).

The data on the neural underpinnings of error correction suggest multiple timing systems, each with specific roles, yet able to coordinate for rapid response. Commensurate with this idea is work suggesting the basal ganglia integrates various timing systems through oscillation comparators (Matell and Meck, 2004; Gu et al., 2015). The limited data on visual error correction, however, leave open how well this network can interface with the visual timing systems.

### EVIDENCE THE AUDITORY SYSTEM HAS PRIVILEGED ACCESS TO TIMING SYSTEMS

Considering the auditory system's timing advantage along with the prominence of the motor system in timing processing, we suggest that the auditory system's advantage in timing stems from its stronger coupling to the motor system. Auditory timing compared to visual timing tasks often yield more activation in motor structures, such as the SMA and premotor cortex (Jäncke et al., 2000). Even when visual SMS tasks employed the modality-appropriate moving visual metronomes, audiomotor synchronization with auditory beeps yielded greater activation in the putamen (Hove et al., 2013b). Likewise, priming a visual rhythm with a similar auditory rhythm resulted in increased putamen activation compared to a visual rhythm alone, while a visual rhythm yielded no priming effect on an auditory rhythm (Grahn et al., 2011). The finding that the increased visual synchronization ability provided by a bouncing ball does not transfer to purely perceptual rhythm perception provides further evidence of the role of motor coupling in timing tasks (Silva and Castro, 2016). Additionally, the privileged link between auditory and motor systems can be seen in Parkinson's disease, a disorder that impairs movement due to cell loss within the basal ganglia (Davie, 2008). For example, Parkinsonian gait can improve when cued by an external rhythm, and these interventions are more effective when synchronizing with auditory metronomes than with flashing visual metronomes (Rochester et al., 2005; Arias and Cudeiro, 2008).

Visual timing activities recruit timing centers within the visual system that, based on behavioral results, are less precise compared to the auditory timing system. In Jäncke et al. (2000), visual timing tasks resulted in increased activity in the right superior cerebellum, vermis, and right inferior parietal lobe compared to auditory timing tasks. Visual timing tasks also recruit areas MT, V5, and the superior parietal lobe, tying into the dorsal visual stream (Jantzen et al., 2005), and visual rhythm perception induces increased beta activity at event onsets arising from the visual cortex (Comstock and Balasubramaniam, 2017b). It is unclear if these timing activations in the visual system are the result of compensating for a weaker connection to the motor timing system. It may be that the temporal processing in the visual system is additional processing of visual information required to interface with the motor system.

While differences in coupling strength to the motor system are crucial for modality timing differences, other factors are likely. To that end, it is clear that the visual system is able to pick out high speed temporal information, for example, V1 will phase lock its input/output to up to a 100 Hz visual flashing stimuli (Williams et al., 2004). This suggest that entrainment is not easily transferred to the systems involved in time/rhythm perception, especially at the time frame usually involved in rhythm perception, indicating that the issue may be one of translation. A likely place for that translation would be within the dorsal pathway, which has been found to have neurons with high temporal resolution in macaques, with higher temporal resolution in the auditory dorsal stream (Rauschecker, 2017). If there is a higher temporal resolution of the auditory dorsal stream than in the visual dorsal stream, then it may give explanation as to why the visual system cannot synchronize at the higher frequencies achieved by the auditory system. Of course, it cannot be ruled out that the difference in temporal resolution is due to different levels of timing precision available to the dorsal stream. Reduced timing precision in the visual stream may be caused by increased necessary processing due to richer sensory input of the visual system compared to the auditory system. Indeed, greater processing requirements and longer processing time may help to account for the inability of the visual system to allow for synchronization at the higher tempos allowed by the auditory system.

#### ROLE OF THE VESTIBULAR-TACTILE-SOMATOSENSORY SYSTEM

Another link between auditory and motor systems is that auditory rhythm perception may be tied to the vestibular-tactilesomatosensory (VTS) system, which is important for movement and dance, and therefore closely tied to the motor system and attuned to timing (Todd and Lee, 2015). In addition to its ties for movement, the VTS system is clearly tied to the auditory system with regards to rhythm perception (Phillips-Silver and Trainor, 2005, 2007, 2008; Trainor et al., 2009), and through common neural activation (Araneda et al., 2017). These ties between the auditory and VTS system may be an additional factor in the dominance of the auditory system in the temporal domain.

Since VTS rhythms are ubiquitous in fetal life through the mother's gait, heart rate, breathing, etc., and since these networks are tied into auditory rhythm systems, it is likely that the VTS system is heavily tied into the timing systems used in auditory rhythm perception and in motor rhythm production (Provasi et al., 2014). This is further strengthened by the fact that movement and rhythms are linked and proprioception (part of the VTS system) plays a large role in perception of rhythms that is tied into auditory rhythm perception and production (Trainor et al., 2009). Interactions between the VTS system with visual rhythm perception remains mostly unexplored at this point however, so it is unclear how much this system plays a supramodal role in the timing involved in rhythm perception/production, or if it is only tied to the auditory and motor rhythm timing systems. Further research in this area is needed to answer these questions.

## EVOLUTIONARY ORIGINS OF SENSORIMOTOR SYNCHRONIZATION

In an evolutionary context, it makes sense that auditory and motor systems would be tightly interconnected. First, rhythms in language are critical for both perception and production and may be a driver of SMS ability (Patel, 2006). Beyond language, matching movement to sound is a necessary result of human evolution that allows for the social and cultural inclination of humanity via music (Hagen and Bryant, 2003; Brown and Jordania, 2013). Dance is also tightly connected with music and culture and can provide a further explanatory account of human SMS capability and the connection between the motor and auditory systems (Fitch, 2016; Iversen, 2016; Laland et al., 2016; Ravignani and Cook, 2016).

Beyond humans, common adaptations appear to increase SMS ability in several non-human species capable of some level of audio-motor entrainment such as parrots (Patel et al., 2009), bonobos (Large and Gray, 2015), and sea-lions (Cook et al., 2013). Although some animals can exhibit rhythmic capabilities, some remarkably well like Ronan the sea-lion (Rouse et al., 2016), they are in some ways limited compared to humans (Patel and Iversen, 2014; Merker et al., 2015). Even though there are animals that can entrain to auditory rhythms, only humans appear to be naturally inclined to do so (Wilson and Cook, 2016). Finally, there is some evidence that non-human primates are able to synchronize their movements to predictable visual stimuli (Takeya et al., 2017), yet there has been much less research on visual SMS compared to auditory SMS in non-humans.

## GENERAL SYNTHESIS AND FUTURE DIRECTIONS

In looking at how the brain processes timing information, it is clear that many context sensitive mechanisms interact and coordinate to provide optimal timing output. Much of this interaction appears to happen within the motor system and likely involves the subcortical systems to coordinate the various mechanisms. Current research suggests that oscillations play a key role coordinating the interactions among various timing circuits. However, it is not clear if the various timing systems compute measures of time in the same way. When considering that auditory and visual systems take in very different kinds of information and use it in different ways, i.e., auditory has a stronger temporal precision, and visual has a strong spatial bias, it seems likely that the timing mechanisms themselves may greatly differ.

Consider the difference between extracting timing information between a moving visual rhythm and an auditory rhythm. Moving visual stimuli contain more information than auditory stimuli, such that while entraining to auditory stimuli, prediction of the onset of the next event involves encoding the interval between two events and utilizing that information to predict the onset of the next event. With a moving visual rhythmic stimulus, that interval information is present, but so is information on position/velocity/acceleration. This means predictions of the onset of the next event can be made as part of a continuous process. The fact that even with this information, visual SMS is at best equal to auditory SMS except at fast speeds, begs the question as to why visual SMS is less capable. One possible explanation for this is that the visual system has to encode much more information, and further, encoding that information into a form that is usable by the motor network may require extra processing. This may explain the timing activity found within the visual cortex during visual SMS. Even when there is a simple flashing metronome, there is a measure of timing activity originating from the visual cortex. Considering the reduced temporal ability with visual flashing metronomes, it suggests there may be a translation issue in harnessing a system not optimized to temporal processing the way the auditory system has been, resulting in a weaker connection to the motor timing network.

Different timing systems likely employ varying mechanisms and computational principles that are appropriate to the time scale, cellular properties, and general needs of the system. Existing computational models that capture a range of these phenomena across levels include: pacemaker accumulator models, multiple oscillator models, memory trace models, random process models, ramping activity models, delay line models, and state space trajectory-based models (Addyman et al., 2016; Hass and Durstewitz, 2016). Such models help illustrate the variety of ways to process timing information within a neural network. Evidence also suggests that cells with specific timing mechanisms exist in the basal ganglia and cerebellum (Lusk et al., 2016), yet other areas with multiple functional properties also process timing, such as in the prefrontal cortex (Hyman et al., 2012) and hippocampus (MacDonald et al., 2011). The areas that have multiple functions, as in the hippocampus and prefrontal cortex, will then likely have different computational approach than more specialized timing structures.

Given that there are multiple ways to process timing, and that many forms of cognition require some form of temporal processing, it would be surprising to find that timing mechanisms are not ubiquitous in the brain. This raises an important question. If many different timing mechanisms are available for a given task, and only one output (through action), how do neural systems arrive at the best timing information to use? A strong candidate explanation for this would implicate a mechanism that helps integration through an optimal Bayesian process (Hass and Durstewitz, 2016). Evidence from multimodal sensory integration suggests that when timing information is presented from multiple modalities, the modalities are combined and weighted based on reliability in Bayesian optimal solution (Ernst and Banks, 2002). Since most timing related activity requires motor output, we would expect that the source of timing to be utilized would be determined before, or as that timing information becomes available to the motor system. This seems to make the case that the striatal cells operating as a comparator may be the seat of the Bayesian process to determine the optimal timing source for motor timing.

Since there is some disparity in the amount of work on auditory and visual SMS error correction, there is a need to further study the error correction capabilities within visual SMS. It is currently unknown if visual error correction can be as fast as auditory error correction when dealing modality appropriate stimuli, such as a moving visual sequence or bouncing ball. Another major area of needed work is in understanding the mechanisms by which the Bayesian optimal timing source is chosen in cases where multiple sources are available. If timing mechanisms are as ubiquitous in the brain as evidence suggests, then there may be a variety of ways these mechanisms interface with the motor timing system to produce a single output. Further imaging and computational work is required to understanding this mechanism.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

This work was partially supported by a grant from the National Science Foundation BCS-1460633.

#### REFERENCES

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Comstock, Hove and Balasubramaniam. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Effect of Stimulus Contrast and Visual Attention on Spike-Gamma Phase Relationship in Macaque Primary Visual Cortex

Aritra Das and Supratim Ray\*

*Centre for Neuroscience, Indian Institute of Science, Bangalore, India*

Brain signals often show rhythmic activity in the so-called gamma range (30–80 Hz), whose magnitude and center frequency are modulated by properties of the visual stimulus such as size and contrast, as well as by cognitive processes such as attention. How gamma rhythm can potentially influence cortical processing remains unclear; previous studies have proposed a scheme called phase coding, in which the intensity of the incoming stimulus is coded in the position of the spike relative to the rhythm. Using chronically implanted microelectrode arrays in the primary visual cortex (area V1) of macaques engaged in an attention task while presenting stimuli of varying contrasts, we tested whether the phase of the gamma rhythm relative to spikes varied as a function of stimulus contrast and attentional state. A previous study had found no evidence of gamma phase coding for either contrast or attention in V1, but in that study spikes and local field potential (LFP) were recorded from the same electrode, due to which spike-gamma phase estimation could have been biased. Further, the filtering operation to obtain LFP could also have biased the gamma phase. By analyzing spikes and LFP from different electrodes, we found a weak but significant effect of attention, but not stimulus contrast, on gamma phase relative to spikes. The results remained consistent even after correcting the filter induced lags, although the absolute magnitude of gamma phase shifted by up to ∼15◦ . Although we found a significant effect of attention, we argue that a small magnitude of phase shift as well as the dependence of phase angles on gamma power and center frequency limits a potential role of gamma in phase coding in V1.

Keywords: attention, spike-field coherence, spike-gamma phase, contrast, area V1, stLFP

# INTRODUCTION

Gamma oscillations are rhythmic fluctuations in a frequency range between 30 and 80 Hz in brain signals (Buzsaki, 2006; Buzsáki et al., 2013), which have been consistently linked with high-level cognitive processes such as attention (Fries et al., 2001; Gregoriou et al., 2009), perception (Rodriguez et al., 1999) and feature binding (Singer, 1999). In recordings from the primary visual cortex (area V1), gamma is also known to be highly dependent on the properties of visual stimulus, such as size (Gieselmann and Thiele, 2008; Ray and Maunsell, 2011a; Jia et al., 2013),

#### Edited by:

*Daya Shankar Gupta, Camden County College, United States*

#### Reviewed by:

*Pedro E. Maldonado, Universidad de Chile, Chile Ko Sakai, University of Tsukuba, Japan*

> \*Correspondence: *Supratim Ray<EMAIL_ADDRESS>
Received: *02 April 2018* Accepted: *20 July 2018* Published: *14 August 2018*

#### Citation:

*Das A and Ray S (2018) Effect of Stimulus Contrast and Visual Attention on Spike-Gamma Phase Relationship in Macaque Primary Visual Cortex. Front. Comput. Neurosci. 12:66. doi: 10.3389/fncom.2018.00066* orientation (Berens et al., 2008; Jia et al., 2011), and contrast (Ray and Maunsell, 2010; Jia et al., 2013). Although several hypotheses about how gamma rhythm could influence neural processing have been proposed, such as binding by synchrony (Singer, 1999) and communication-through-coherence (Fries, 2015), whether gamma plays a functional role remains unclear (Ray and Maunsell, 2015).

Here we test a specific hypothesis called phase coding (PC), originally proposed in the context of theta rhythm in the hippocampus (O'Keefe and Recce, 1993; Buzsáki and Chrobak, 1995), in which information is coded in the position of the spike relative to the rhythm. In the context of gamma rhythm (Fries et al., 2007), which is thought to be associated with an inhibitory network of interneurons (Bartos et al., 2007; Cardin et al., 2009; Sohal et al., 2009), this hypothesis posits that the rhythmic network inhibition interacts with excitatory input to pyramidal cells such that the more excited cells (which can overcome the inhibition earlier) fire earlier in the gamma cycle. Thus, stimulus intensity can be coded in the gamma phase relative to the spike. However, whether gamma PC occurs is controversial, with evidence both in favor and against the hypothesis. We have earlier shown that in macaque secondary somatosensory cortex, the phase of gamma rhythm does not vary with stimulus intensity (Ray et al., 2008). In V1, one study showed some evidence of PC with different orientations (which they took as a proxy for stimulus intensity), at least for sites that had weak gamma power and weak gamma-spike phase locking (Vinck et al., 2010). Other studies in V1 showed no evidence of PC when the stimulus contrast (a more direct index of stimulus intensity as compared to orientation) was varied (Chalk et al., 2010; Ray and Maunsell, 2010). Importantly, Chalk and colleagues further showed that even attention, which increases the effective contrast of the stimulus (Carrasco et al., 2004), does not cause a shift in spike-gamma phase. They also showed that in V1, attention causes a reduction in gamma power and spike-gamma coupling (Chalk et al., 2010), exactly opposite of what has been shown in higher cortical areas such as V4 (Fries et al., 2001).

While Chalk et al. (2010) failed to provide evidence in favor of PC, one limitation of their study was that spikes and local field potential (LFP) were collected from the same electrode, which can potentially bias the spike-gamma phase relationship because of the presence of spike-related transients (Ray, 2015). Specifically, spikes are associated with a "transient" in the LFP recorded from the same electrode, which could be due to synaptic activity that leads to the spike as well as low-frequency component of the action potential (spike "bleedthrough"; see Ray, 2015, for details). The remaining studies either removed this transient using signal processing techniques such as Matching Pursuit (Ray et al., 2008), or used spikes and LFP from different electrodes (Ray and Maunsell, 2010; Vinck et al., 2010) that reduces the bias (see section Discussion for more details on this), but none of these reports studied the effect of attention on spike-gamma phase. To test whether stimulus contrast can be coded in the phase of the gamma rhythm, we here trained monkeys to do a demanding attention task while presenting stimuli that varied in contrast to study the effect of both contrast and attention on spike-gamma phase, while recording from chronically implanted microelectrode arrays such that spike-gamma phase could be estimated using spikes and LFPs recorded from different electrodes. Note that since attention is thought to increase the effective stimulus contrast (Carrasco et al., 2004), testing whether gamma phase varies with attentional state is also a test for gamma PC for contrast. We further studied the effect of the online causal filter used to obtain the LFP, which introduces a delay in the LFP and has been shown to influence spike-LFP relationships (Okun, 2017), but has not been accounted for in previous studies.

# MATERIALS AND METHODS

Experimental procedures have been described in detail in earlier studies (Ray and Maunsell, 2010, 2011b; Shirhatti et al., 2016); we provide a brief description here.

## Ethics Statement

The animal protocols reported in this study were approved by the Institutional Animal Care and Use Committee of Harvard Medical School.

 Once they learned the behavioral task, a microelectrode array (Blackrock Microsystems, 96 active electrodes) was implanted in the right V1 cortex (about 15 mm anterior from the occipital ridge and 15 mm lateral from the midline). The microelectrodes were 1 mm long and 400µm apart from each other, with impedance between 0.3 and 1 M at 1 kHz. Although histological analysis had not been performed to identify the exact location of the microelectrode tips, they are expected to be in cortical layer 2/3 or 4 based on the approximate thickness of V1 (2 mm; Hubel and Wiesel, 1977). Electrical signals were recorded using commercial hardware and software (Blackrock Microsystems), referenced to a wire placed on the dura near the microelectrode grid. Raw electrical signals were filtered between 0.3 Hz (Butterworth filter, 1st order, analog) and 500 Hz (Butterworth, 4th order, digital) and digitized at 2 kHz (16-bit resolution) to get the LFP. Multi-units were extracted by filtering the raw signal between 250 Hz (Butterworth filter, fourth order, digital) and 7,500 Hz (Butterworth filter, third order, analog) followed by an amplitude threshold (set at ∼6.25 and ∼4.25 of the signal SD for the two monkeys). To improve the quality of unit isolation, multi-units were subsequently sorted offline (Offline Sorter, Plexon Inc.). The receptive fields, obtained by flashing small Gabor stimuli on a rectangular grid that encompassed the receptive fields of all the electrodes in the array, were located in the lower left quadrant of the visual space at an eccentricity of about 3–5◦ . As in previous studies, only electrodes for which stable estimates of the receptive fields could be obtained (27 and 66 electrodes for the two monkeys), were used for subsequent analysis.

#### Behavioral Task Paradigm

The monkeys were required to maintain their gaze within 1◦ of a small central dot (0.05◦ -0.10◦ diameter) during the task while two achromatic odd-symmetric static Gabor stimuli were synchronously flashed for 400 ms with a mean inter-stimulus period of 600 ms. One of the two Gabor stimuli was centered on the receptive field of one of the recorded sites (new location for every session) while the second stimulus was located at an equal eccentricity on the opposite side of the central fixation point. The monkeys were cued to pay attention to one of the two stimulus locations in different blocks of trials by presenting two instruction trials (not included in the analysis) at the start of the block, in which there was only a single stimulus. The contrasts of the attended and unattended Gabor stimuli were equal on each presentation and could take any of the eight possible values: 0, 1.6, 3.1, 6.2, 12.5, 25, 50, and 100%, chosen pseudorandomly. At an unstipulated time drawn from an exponential distribution (mean 2,000 ms, range 1,000–7,000 ms for Monkey 1; mean 3,000 ms, range 1,000–7,000 ms for Monkey 2), the orientation of the stimulus at the cued location changed by 90◦ . An exponential distribution was used to minimize expectation of target appearance and to keep the attentional state uniform during a trial since the hazard function is flat for an exponentially distributed target onset time. The monkeys were rewarded with a drop of juice for making a saccade to the location of the altered stimulus within 500 ms of orientation change. To account for saccade latency and to minimize guessing, monkeys were rewarded only for saccades beginning at least 100 ms after the orientation change. Trials were terminated at 7,000 ms if the target had not appeared, in which case the monkeys were rewarded for maintaining fixation throughout that trial. These catch trials were excluded from analysis (for more details, see Ray and Maunsell, 2010 and **Figure 1**).

The Gabor stimuli used for this task were both static with SD of 0.5◦ , spatial frequency of 4 cycles per degree, with one of the Gabor stimuli located at the center of the receptive field of one of the recorded sites (new recording site for each session), at its preferred orientation. Data from the two monkeys were collected in 10 and 17 recording sessions, respectively. Only correct trials were used for analysis. For each of the correct trials, only the second stimulus up to the last stimulus before target onset were used for analysis. We only used stimulus contrasts for which salient gamma oscillations were observed (25, 50, and 100% contrasts). For each contrast and attention condition, on average we obtained 79 ± 4 (range 55–101) stimulus repeats for Monkey 1 and 74 ± 5 (range 47–120) for Monkey 2.

### Electrodes and Electrode Pair Selection

Electrodes with receptive field centers within 0.2◦ of the stimulus center in each of the recording sessions were used for analysis, yielding 63 electrodes (23 unique; many electrodes were selected in several recording sessions) for Monkey 1 and 89 electrodes (53 unique) for Monkey 2. These are referred to as "LFP" electrodes. For spike-field coherence (SFC), spike-triggered LFP (stLFP) and spike-gamma phase histograms, we selected a subset of the LFP electrodes from which at least 20 spikes could be recorded in the analysis interval (150–400 ms after stimulus onset) and the

circle for clarity; not visible to the monkey) and the other stimulus appeared at a location of equal eccentricity in the opposite hemifield. The monkey was cued to covertly attend to one of the two locations in different blocks of trials (indicated by black dotted circle, not visible to the monkey). At an unsignaled time, during one of the stimulus presentations, the orientation of the cued stimulus was changed by 90◦ . The monkey was rewarded with a drop of juice for making a saccade to the location of orientation change. If there was no change during a trial (catch trial), the monkey was rewarded for maintaining fixation throughout that trial.

signal to noise ratio of the isolation (Kelly et al., 2007) was greater than 2. This generated 23 (12 unique) and 39 (27 unique) "spike" electrodes for Monkeys 1 and 2, respectively. For each session, we took all combinations of spike and LFP electrodes with receptive fields within 0.2◦ of the stimulus center, yielding 23 (12 unique) and 39 (27 unique) "same" spike-LFP pairs (**Figure 3**), and 163 (120 unique) and 170 (147 unique) pairs of "different" spike-LFP electrodes for Monkeys 1 and 2, respectively (**Figures 4**, **5**).

# DATA ANALYSIS

All data were analyzed using custom codes written in MATLAB (The MathWorks, RRID:SCR\_001622). Individual data analysis methods are briefly summarized below.

# Change in Power Spectral Density (PSD) Plots (Figure 2)

Stimulus-induced responses were first obtained by subtracting the mean LFP across all stimulus repeats for each condition (i.e., the event-related potential) from individual single trial time series data. Subsequent analyses were performed on these stimulus-induced responses. Power spectral densities (PSDs) for different stimulus and attention conditions were computed using the multi-taper method with 5 tapers using the Chronux toolbox (Bokil et al., 2010); http://chronux.org/, (RRID:SCR\_005547).

triangles and circles respectively. (B) Change in alpha power (in decibels) between attend-in and attend-out condition for 25% (blue), 50% (green), and 100% (red) contrasts. Error bars represent standard error of mean across 63 electrodes from Monkey 1 (top row) and 89 electrodes from Monkey 2 (bottom row). Significant differences (*p* < 0.05, *t*-test) are indicated by "\*". (C) Same as (B), but for gamma power. (D) Change in gamma peak frequency (Hz) between attend-in and

attend-out condition for contrasts 25% (blue), 50% (green), and 100% (red). Error bars represent standard error of mean across 63 electrodes from Monkey 1 (top row) and 89 electrodes from Monkey 2 (bottom row).

The analysis period was selected between 150 and 400 ms after stimulus onset to avoid stimulus onset related transients and compared against a "baseline period" between −300 and −50 ms of stimulus onset. To ensure that the change in power from baseline was not affected due to differences in the baseline power for different attention conditions, change in PSDs were plotted with respect to the baseline response of "attend-out" (attention directed outside the receptive field) condition for each stimulus contrast value:

$$
\Delta PSD\_i = 10 \left( \log\_{10} \left( \text{ST} \right)\_i - \log\_{10} \left( BL\_{Att\ Out} \right)\_i \right)\_i
$$

Here i represents the contrast condition (25, 50, or 100%), 1PSDi represents the change in PSD in decibels, (ST)i denotes the PSD in the stimulus epoch and (BLAtt Out)i denotes the baseline PSD for attend-out condition.

For the change in alpha power shown in **Figure 2B**, we first averaged the power between 8 and 12 Hz (triangles in **Figure 2A**; note that because we used an analysis interval of 250 ms, we had a frequency resolution of 4 Hz) and subsequently took the difference between the "attend-in" (attention directed inside the receptive field) and attend-out power on a log scale:

$$
\DeltaPower\_{\overline{i}} = 10 \left( \log\_{10} \left( ST\_{Att \, \, \text{Ir}} \right)\_{i} - \log\_{10} \left( ST\_{Att \, \, \text{Out}} \right)\_{i} \right)\_{i}
$$

Here (ST) i denotes the alpha power in the stimulus epoch for the ith contrast condition. For gamma power, the same procedure was used with three frequency bins centered around the peak gamma frequency (shown in circles in **Figure 2A**). For computing peak gamma frequency, we choose the frequency bin for which 1PSDi attained its maximum value between 30 and 60 Hz.

#### Coherency Analysis (Figures 3–5)

The coherency between two signals x and y is computed using the following equation:

$$\text{Coherence}\_{\text{xy}}\left(f\right) = \frac{\text{S}\_{\text{xy}}(f)}{\sqrt{\text{S}\_{\text{xx}}(f)\text{S}\_{\text{yy}}(f)}}$$

Where Sxy (f) denotes the cross-spectrum between the signals x and y and Sxx (f) and Syy (f) denote the auto spectra of each signal. The coherency values were computed using the multitaper method implemented in Chronux toolbox using five tapers. All the coherence analyses were performed using the sorted multiunit dataset. For spike-field coherence, the spike time series was converted to a binary time series (at 0.5 ms resolution) with a "1" at each time position containing a spike and "0" otherwise (500 data points for the stimulus period). The results were similar for three tapers. All the circular statistical analyses were performed using an open source circular statistics toolbox in MATLAB (CircStat; Berens, 2009). Spike-triggered LFP (stLFP) were computed by taking a ±25 ms segment of the LFP around each spike in the stimulus period and subsequently taking the average of those segments.

# Removing Filtering Effect (Figure 5)

obtained from Watson-Williams test to compare the mean phases.

Any causal filter necessarily introduces a delay in the signal, which may be dependent on the frequency of the signal. Butterworth filters have a linear relationship between phase delay and frequency, such that the group delay (which roughly translates to how much each frequency component of the signal shifts in time due to the filtering process) is almost constant over a wide frequency range. We removed the low-pass filtering effect by dividing the Fourier Transform of the LFP by the Fourier Transform of the low-pass LFP filter (4th order Butterworth filter with a low-pass cutoff at 500 Hz; constructed in MATLAB using the command "butter") and subsequently taking the inverse Fourier Transform (Okun, 2017). The correction was only done between 0 and 500 Hz because the power of the LFP (as well as the filter) was very less beyond 500 Hz. The group delay of this Butterworth filter was ∼0.8 ms over almost the entire frequency range of interest (including the gamma range), such that the stLFP constructed from the corrected LFP signal had a similar shape as the uncorrected stLFP but was shifted leftward by ∼0.8 ms (**Figure 5B** vs. **Figure 4B**).

Note that in addition to this low-pass filter, three other filtering operations also need to be accounted for. The data acquisition system had two analog hardware filters: a high pass filter at 0.3 Hz (first order, Butterworth) and a low-pass filter at 7,500 Hz (third order, Butterworth). In addition, to obtain spike data, the signal was high-pass filtered at 250 Hz (fourth order, Butterworth, digital). However, all three filters had negligible group delay (<0.1 ms) between 500–5,000 Hz, suggesting that these filtering operations did not change the position of the spike appreciably. Similarly, the high-pass filter had a large group delay at very low frequencies, as shown by Okun (2017), but it was negligible in the gamma range. Therefore, these three filters did not have an appreciable effect on the spike-gamma phase estimation.

#### RESULTS

## Spatial Attention Reduces Alpha and Gamma Power and Increases Peak Gamma Frequency in Area V1

We first analyzed changes in alpha and gamma power and gamma peak frequency for attend-in versus attend-out conditions to test whether our results were in agreement with previous attention studies in macaque primary visual cortex (Chalk et al., 2010). **Figure 2A** shows the average change in PSD during the stimulus

FIGURE 4 | Relationship between Spikes and LFPs recorded from different electrodes, as a function of Contrast and Attention. Same as Figure 3, but the analysis is performed on 163 and 170 pairs of different spike-LFP electrodes (the receptive fields of both were within 0.2◦ ) for Monkeys 1 and 2, respectively.

period (150–400 ms after stimulus onset) from the pre-stimulus baseline (−300 to −50 ms), for attend-out (black trace) and attend-in condition (color trace) for 25% (blue), 50% (green), and 100% (red) contrasts. Gamma peak frequency increased with increasing contrast (traces for different contrasts are overlaid in the rightmost plot for comparison), as reported previously (Ray and Maunsell, 2010). To account for potential differences in the baseline activity due to attention, all changes were computed with respect to the baseline activity of the unattended condition (see section Materials and Methods). Consistent with previous studies, we found a strong suppression of alpha power due to attention in both monkeys, which could be observed in the baseline PSD as well (dark gray trace), confirming that the monkeys were indeed attending to the stimuli.

To test these results quantitatively, we first performed a threeway ANOVA test with factors of monkey (2 levels: Monkey 1 and 2), attention (attend-out, attend-in) and contrast (25, 50, and 100%). Alpha power was averaged over 8 and 12 Hz (inverted triangles in **Figure 2A**), while gamma power was averaged in an eight Hz band around the peak frequency for each contrast (40, 48, and 56 Hz for Monkey 1, and 40, 44, and 56 Hz for Monkey 2; **Figure 2A**; see section Materials and Methods for details). The factor monkey was significant for alpha power (Fmonkey = 243.23, p = 9.3 × 10−49), gamma power (Fmonkey = 417.45, p = 1.3 × 10−76), and peak gamma frequency (Fmonkey = 4.52, p = 0.004). Thus, we performed a two-way ANOVA with factors attention and contrast separately for the two monkeys for alpha power, gamma power and peak gamma frequency. The effect of contrast on alpha power was not significant for Monkey 1 (Fcontrast = 0.65, p = 0.522) but significant effect for Monkey 2 (Fcontrast = 7.06, p = 9 × 10−4 ). The effect of attention on alpha power was significant for both monkeys (Fattention = 24.88, p = 9.4 × 10−7 for Monkey 1 and Fattention = 10.17, p = 1.5 × 10−3 for Monkey 2). However, there was no significant interaction between the two factors on alpha power for either monkey (Fcontrast × attention = 0.87, p = 0.42 for Monkey 1 and Fcontrast × attention = 0.91, p = 0.40 for Monkey 2). For gamma power, the effect of contrast was significant for both monkeys (Fcontrast = 11.9, p = 9.8 × 10−6 for Monkey 1 and Fcontrast = 37.78, p = 4.6 × 10−16 for Monkey 2) but the effect of attention was not significant (Fattention = 0.3, p = 0.58 for Monkey 1 and Fattention = 1.72, p = 0.19 for Monkey 2). Again, there was no interaction between the factors (Fcontrast × attention = 0.05, p = 0.95 for Monkey 1 and Fcontrast × attention = 0.17, p = 0.85 for Monkey 2). For peak gamma frequency, there was significant effect of both contrast and attention in both monkeys (Fcontrast = 415.46, p = 1.6 × 10−95 for Monkey 1 and Fcontrast = 973.36, p = 7.7 × 10−178 for Monkey 2; Fattention = 5.07, p = 0.025 for Monkey 1 and Fattention = 33.31, p = 1.4 × 10−8 for Monkey 2). The interaction of the two factors was significant only in Monkey 2 (Fcontrast × attention = 1.11, p = 0.33 for Monkey 1 and Fcontrast × attention = 3.37, p = 0.04 for Monkey 2). Similar results were obtained in a 2-factor ANOVA performed on the data pooled over the two monkeys.

These results were further confirmed using pairwise t-tests. In almost all conditions, alpha power significantly reduced with attention (**Figure 2B**; Monkey 1: t-test, N = 63, p = 3.9 × 10−7 , 1.1 × 10−13, 3.5 × 10−5 for 25, 50, and 100% contrasts respectively; p = 8.9 × 10−22 for all the contrast conditions combined; Monkey 2: t-test, N = 89, p = 0.057, 2.4 × 10−14 , 3.1 × 10−13 and 8.6 × 10−20 for 25, 50, 100% contrasts and the combined condition, respectively). The reduction in gamma power was significant only for the 50% contrast condition for Monkey 1 (t-test, N = 63, p = 0.41, 0.013, 0.86, and 0.052 for 25, 50, 100% and combined contrast conditions, respectively), and for 25 and 100% contrasts for Monkey 2 (t-test, N = 89, p = 3.4 × 10−8 , 0.72, 0.014 and 1.5 × 10−6 for 25, 50, 100% and combined contrast conditions, respectively). Similarly, the increase in gamma peak frequency was significant only for 25% contrast condition for Monkey 1 (t-test, N = 63, p = 0.007, 0.15, 0.13, and 7.8 × 10−4 for 25, 50, 100% and combined contrast conditions, respectively), and all contrasts for Monkey 2 (t-test, N = 89, p = 1.7 × 10−9 , 6.7 × 10−11, 0.013 and 6.7 × 10−19 for 25, 50, 100% and combined contrast conditions, respectively).

The weak effect of attention on gamma power is not surprising for two reasons. First, because we recorded from a chronically implanted microelectrode array, the stimuli were optimized only for a single site in each session and therefore were non-optimal for most electrodes, unlike the study by Chalk et al. (2010) where stimuli were better optimized. Second, since the stimuli were only presented for 400 ms (to minimize attentional fluctuations within the stimulus duration) and the analysis duration was only 250 ms, the frequency resolution was 4 Hz, which made it difficult to correctly estimate peak frequency shifts that are typically only 2–3 Hz in V1 (Ray and Maunsell, 2010; Bosman et al., 2012). Note that the second limitation can be partially overcome by using Matching Pursuit (Chandran et al., 2016), which allowed us to better characterize the gamma peak frequency shifts in a previous study (see Supplementary Figure 2 of Ray and Maunsell, 2010); we have used multi-taper analysis here because the spike-field coherence (SFC), which was also used to get spike-gamma phase, was obtained using the same technique. In general, the effects of attention on V1 were consistent with the findings of Chalk et al. (2010), and were almost always significant when the results were pooled across contrasts.

### Effect of Contrast and Attention on SFC, stLFP, and Spike-Gamma Phase Computed Using Spike and LFP Recorded From the Same Electrode

Next, we analyzed how attention modulated SFC, stLFP and spike-gamma phase when spikes and LFP were recorded from the same electrode (23 and 39 sites for the two monkeys; see Materials and Methods for details), as was the case in the study by Chalk et al. (2010). The magnitude of the SFC (**Figure 3A**) showed clear peaks in the gamma frequency range, and the peak gamma frequency shifted with an increase in contrast. Consistent with Chalk et al. (2010) and the results obtained using power (**Figure 2**), we found a reduction in SFC magnitude and an increase in peak gamma frequency with attention in almost all conditions. The stLFP plots (**Figure 3B**) showed the presence of a prominent rhythm around the time of the spike, especially for Monkey 2, whose trough was shifted 3–4 ms away from zero. These results were reflected in the spike-gamma phase histograms (**Figure 3C**), obtained by taking the average angle of the SFC across the three frequency bins around the peak gamma frequency (as highlighted in **Figure 3A**). Following the convention used by Chalk and colleagues, phase angles were defined such that trough of the gamma rhythm was at 180◦ and rightward shift of the trough increased the phase angle. The mean phase angles were ∼210◦ for Monkey 1 and ∼235◦ for Monkey 2 and were not significantly different across attention conditions (circular mean phases and the associated p-values obtained from Watson-Williams test are shown in the legend). Even when pooled across contrasts, the mean phases between attend-in and attend-out conditions were not significantly different (Watson-Williams multi-sample test, p = 0.28 and 0.33 for Monkeys 1 and 2). Similarly, the mean phases at different contrasts were not significantly different from each other in either attend-out or attend-in conditions (Watson-Williams multi-sample test, p = 0.89 (attend-out) and p = 0.9 (attend-in) for Monkey 1; p = 0.66 (attend-out) and p = 0.55 (attend-in) for Monkey 2). These results are consistent with Chalk et al. (2010), who obtained a median phase of ∼-0.65π, which translates to ∼243◦ . An offset of ∼30◦ -50◦ from the trough (180◦ ) is also consistent with the findings of Vinck et al. (2010) and Ray and Maunsell (2010), although in these two studies the convention was chosen such that rightward shift of the trough led to a reduction of phase angle below 180◦ (such that the phase angles were between ∼130◦ and ∼150◦ ).

Although our results are consistent with previous studies, there are two serious flaws in these results, which can be clearly observed in the stLFP plots (**Figure 3B**). First, there is a large spike-related transient (sharp negative dip near time zero), which biases the estimation of the gamma phase. Specifically, this transient can be decomposed into a series of sinusoids with their troughs aligned to the trough of the transient, effectively "pulling" the phase of any true phase-locked rhythm toward 180◦ (for a detailed discussion, see Ray, 2015). This can be observed in the two monkeys: the estimated spike-gamma phase is closer to 180◦ for Monkey 1 compared to Monkey 2 (∼210◦ vs. ∼235◦ ), simply because the relative magnitude of the transient compared to the gamma rhythm is larger for Monkey 1. The second flaw is that the spike-related transient, which should be around the time of the spike itself, is shifted toward positive values. We address both these concerns below.

## Effect of Using Different Electrodes for Spikes and LFP on SFC, stLFP, and Spike-Gamma Phase

One popular method to reduce the spike-related transient is to take spikes and LFP from different electrodes (Ray and Maunsell, 2010; Vinck et al., 2010; Ray, 2015). We, therefore, repeated the analysis on 163 and 170 "different" spike-LFP pairs for the two monkeys, such that the receptive fields of both were located within 0.2◦ of stimulus center (see section Materials and Methods for details). Mean SFC showed similar results as before, with clear peaks in the gamma frequency range and an increase in peak gamma frequency with increasing contrast, and a slight reduction in SFC magnitude and an increase in peak frequency with attention in some cases. Spike-related transient, which was prominent in **Figure 3B**, was now much reduced, better revealing the true gamma rhythm in the stLFP (**Figure 4B**) whose trough was 3–4 ms after the spike in both the monkeys. Mean spike-gamma phases were now ∼235◦ and ∼245◦ for the two monkeys (**Figure 4C**; note that the shift in mean phase between **Figures 3**, **4** is much larger for Monkey 1 because the spike transient was relatively much larger for that monkey). Interestingly, for both monkeys and for all contrast conditions, attention appeared to shift the mean gamma-spike phase away from 180◦ . Although this phase difference did not reach significance for many contrast levels (circular means and p-values obtained using Watson-Williams multi-sample test are shown in the bottom of **Figure 4C**), the phase differences were highly significant when combined across contrasts (Watson-Williams multi-sample test, p = 7 × 10−5 and 5 × 10−3 for Monkeys 1and 2), albeit the actual magnitude of the difference was small (∼19◦ and ∼6 ◦ ). The mean phases at different contrasts were not significantly different from each other in either the attend-out or the attend-in condition [Watson-Williams multi-sample test, p= 0.36 (attend-out), p = 0.11 (attend-in) for Monkey 1 and p = 0.7 (attend-out), p = 0.34 (attend-in) for Monkey 2].

## Effect of Removing the Filtering Artifact on SFC, stLFP, and Spike-Gamma Phase Relation

The rightward shift of the spike-related transient away from zero (**Figures 3B**, **4B**) is simply due to the effect of the filtering operation to obtain the LFP. We, therefore, removed this filtering effect (see Materials and Methods for details) and reanalyzed SFC, stLFP and spike-gamma phase for "different" pair condition (**Figure 5**; for the "same" electrode condition, this operation caused the trough of the spike-transient to shift near zero; data not shown). While this operation did not change any of the results shown in **Figure 4**, the mean phases decreased by ∼10◦ at both 25 and 50% and ∼13◦ at 100% contrast (for the same shift in time, the shift in degrees depends on the frequency of the rhythm; **Figures 4B,C** are overlaid as dashed-dot traces on the corresponding panels in **Figure 5** to show the outcome of filtering-effect removal). Otherwise, like **Figure 4**, the effect of attention on spike-gamma phase remained significant when phases were pooled across contrast conditions (Watson-Williams multi-sample test, p = 6 × 10−5 for Monkey 1 and 8 × 10−3 for Monkey 2). Similarly, the mean phases at different contrasts were not significantly different from each other in either the attend-out or the attend-in condition [Watson-Williams multisample test, p = 0.58 (attend-out), p = 0.07 (attend-in) for Monkey 1 and p = 0.69 (attend-out), p = 0.73 (attend-in) for Monkey 2].

#### DISCUSSION

We investigated whether increasing stimulus contrast or allocating more attention to a stimulus (which increases its effective contrast) shifts the position of the spike relative to the phase of the gamma rhythm, as posited by the PC hypothesis. We highlighted two issues that can bias the phase estimation: the presence of the spike-related transient and the effect of filtering to obtain the LFP. After accounting for these issues, we found no effect of stimulus contrast and a weak but significant effect of attention on spike-gamma phase. Although these results are consistent with the PC hypothesis in the context of attention, we discuss three issues that severely limit the efficacy of gamma PC in V1.

#### Issue 1: Magnitude of Gamma PC in V1

For a rhythm occurring at 50 Hz (time period of 20 ms), the interval between the peak and the subsequent trough (the interval over which the inhibition fades away) is 10 ms, which is the maximum range over which PC can operate. It is clear from our results, as well as prior reports, that even if PC occurs, it only uses a small sub-interval within this interval. Since spikes occur away from the trough of the rhythm with increasing stimulus intensity under PC, the delay of the trough from the spike at 100% contrast sets the dynamic range of this coding scheme. In our data, spikes occurred at ∼230◦ at 100% contrast, similar to the value reported by Chalk et al. (2010) (∼240◦ ) and Vinck et al. (2010) (∼137◦ for preferred orientation, which translates to ∼223◦ as per our convention). For a rhythm at ∼50 Hz, a shift of ∼50◦ translates to only ∼3 ms out of the available ∼10 ms for coding. Further, even when contrast was reduced to 25%, there was no discernable change in the trough position. The only study that did show any evidence of phase coding (Vinck et al., 2010) showed a shift of ∼20◦ between the best and worst orientation, which translates to only ∼1 ms shift (in addition, see other issues with their results below). In our data, the shift in phase due to attention is even lesser, especially for Monkey 2 (in addition, see Issue 3 below). It can be argued that the phase could shift down to 180◦ for very low contrasts (providing a dynamic range of ∼3 ms), but it is well known that gamma rhythm itself is weak or absent at very low contrasts (Henrie and Shapley, 2005; Jia et al., 2013) and also peaks at a lower frequency (see Issue 3). Thus, if we consider the range of contrasts for which gamma is reliable, the magnitude of PC (i.e., the range over which the spike varies with respect to the rhythm) appears to be very small in V1. In this context, our filtering correction becomes significant, since even though the group delay is only ∼0.8 ms, it still decreases the dynamics range by a further ∼20–25%.

## Issue 2: Effect of Changing Gamma Amplitude

As shown in **Figure 3**, spikes are associated with a transient in the LFP recorded from the same electrode, which biases the spike-LFP phase analysis. Because the amplitude of an extracellular action potential generally decreases rapidly as a microelectrode is moved away from the neuron (Gold et al., 2006; Schomburg et al., 2012), the spatial spread of a spike is thought to be very local (for example, Xing and colleagues used a range between 30 and 100µm for single units; Xing et al., 2009). Therefore, one way to reduce the spike transient is to take the LFP from a neighboring electrode that is separated from the spike electrode by at least a few hundred microns (for a representative case, see Vinck et al., 2010). There are, however, two issues with this approach. First, although taking spikes and LFPs from different electrodes drastically reduces the spike-related transient, it does not completely eliminate it (Ray, 2015). For example, as shown in Figures 2A,E of Ray and Maunsell (2011b) where stLFPs were constructed using spikes and LFP electrodes separated by different distances for the same two monkeys as used in this study, the spike-transient could be seen up to electrode pairs separated by ∼400µm for Monkey 1 and ∼0.4–1.6 mm for Monkey 2, albeit the magnitude of the spike-transient was much smaller than when stLFP was constructed from the same electrode (d = 0 condition in those plots). This happens because neurons near the LFP electrode are often correlated with the neuron being recorded from the spike electrode, and those neurons produce a transient in the LFP electrode that are locked to the spikes on the spike electrode. The second issue is that this procedure implicitly assumes that gamma oscillations recorded from two nearby electrodes are similar, but the spatial spread of LFP itself is a topic of debate. While some studies have shown that the spatial spread of LFP could be large (up to a few mm; Kajikawa and Schroeder, 2011), others have shown that it could be only a few hundred microns (Katzner et al., 2009; Xing et al., 2009; Dubey and Ray, 2016). A modeling study showed that the spread could depend on the degree of correlation in the neural population (Lindén et al., 2011). Consequently, there might be differences in the gamma recorded from neighboring microelectrodes. For example, we have shown that when a Gabor stimulus is presented, two microelectrodes separated by as little as 0.2◦ can exhibit significantly different center frequencies (Ray and Maunsell, 2010). Therefore, some studies have used other techniques to remove the spike-transient, such as Matching Pursuit (Ray et al., 2008) or a Bayesian Framework (Zanos et al., 2011). All these methods substantially reduce the spike-transient, although it is unlikely that they completely eliminate it.

A small spike-transient is unlikely to influence the estimation of gamma phase when the rhythm itself is very strong but may shift the phase toward 180◦ when the rhythm is weak. For example, even when spikes and LFPs were recorded from separate electrodes (**Figure 4**), the mean phases for Monkey 1 were about ∼10◦ less than Monkey 2, who had a much stronger gamma rhythm than Monkey 1. A visual inspection of the stLFP (**Figure 4B**) reveals a small spike-transient like structure in Monkey 1, which could have contributed to the reduction in spike-gamma phase as compared to Monkey 2. Importantly, in cases where the magnitude of gamma itself varies across conditions, an apparent shift in spike-gamma phase could just be due to a differential contribution of the spike-transient which "pulls" the phase toward 180◦ . For example, Vinck et al. (2010) showed that gamma PC was stronger when gamma power and gamma phase locking was very weak (see their Figure 6). Because they did not show the stLFPs, it is unclear whether the apparent phase shift they documented was because of a genuine leftward shift of the gamma trough or the presence of a spike-transient whose contribution was larger when the gamma rhythm itself was weak.

## Issue 3: Effect of Changing Gamma Peak Frequency

PC hypothesis makes sense when the rhythm has a stable frequency. However, the center frequency of gamma rhythm varies systematically with changes in a variety of stimulus parameters, such as size (Gieselmann and Thiele, 2008; Ray and Maunsell, 2011a; Jia et al., 2013), contrast (Ray and Maunsell, 2010; Bosman et al., 2012; Jia et al., 2013), and drift rates (Gray and Viana Di Prisco, 1997; Friedman-Hill et al., 2000). For example, although we show that the spike-gamma phase angles do not vary with stimulus contrast, note that these angles are computed for different gamma frequencies, making it harder to interpret and compare these phase values. Vinck et al. (2010) used different orientations for comparison, but gamma center frequencies can vary even for different orientations, although the trends are not always consistent (see Figure 2D of Jia et al., 2013 and Figure 1 of Murty et al., 2018). For the same delay between the spike and gamma trough, the effective phase angle is greater when the rhythm is faster. For example, in our data, the stLFP troughs appear to coincide between the attend-in and attendout cases in almost all conditions (**Figure 5B**). However, since attention slightly increases the gamma frequency, the effective phase lag in degrees could be larger, which could explain the small but consistent increase in phase angles.

We note, however, that we computed phase over a 250 ms window (similar results were obtained for 200 ms window), which cover more than 10 cycles of the rhythm. During natural vision, we make 3–4 saccades every second (even during fixation, we make several micro-saccades per second), and such eye movements can change or reset the phase of LFP oscillations (Bosman et al., 2009; Ito et al., 2011). It is possible that PC occurs within a single or a few cycles of gamma rhythm, for which gamma need not even have stable frequency over time. It is also possible that PC occurs differently during natural viewing as opposed to a paradigm where animals are trained to fixate for long durations. For example, Ito and colleagues showed that in freely viewing monkeys, fixation-related spike synchronization occurred at an early phase of the rate response after fixationonset, and the first spikes after the onset of a fixation were locked to a specific epoch of the LFP modulation (Ito et al., 2011). Other studies have also shown that gamma rhythm tends to appear in short bursts over a few cycles (Xing et al., 2012; Lundqvist et al., 2016; Chandran Ks et al., 2017), and therefore PC could theoretically occur over shorter duration than what was considered here. Comparable recordings from monkeys during natural viewing conditions as well as advanced signal processing techniques are required to test this hypothesis.

## Weak Effect of Attention in V1

The effect of attention was weaker in our data than the findings of Chalk et al. (2010), possibly due to the use of sub-optimal stimuli for many sites, fewer sites, and a shorter analysis window. However, it is unlikely that our results would change drastically if these limitations could be overcome. First, the effect of attention on gamma in V1 is in general weak (Chalk et al., 2010; Buffalo et al., 2011). Second, although the reduction in gamma power and SFC with attention were small, we obtained a pronounced reduction in alpha power in all cases. Similarly, the increase in gamma peak frequency (1–3 Hz in our data) was comparable to a previous study by Bosman and colleagues, who reported an increase in gamma peak frequency of 2–3 Hz (Bosman et al., 2012). Third, although the analysis window was shorter than previous studies, which yielded a poor frequency resolution, the stLFP plots were computed in the timedomain itself and therefore did not suffer from the poor frequency resolution, but even these did not show a substantial rightward shift as is expected from the PC hypothesis. Fourth, while we had fewer recording sites that may have yielded less statistical power for power analysis (**Figure 2**), we had a substantial number of pairs (163 and 170 for the two monkeys), so the main result regarding the PC hypothesis (**Figure 5**) did not suffer from the lack of statistical power. Finally, although the effect of attention was weak, contrast had a strong effect on gamma power and frequency, but the PC hypothesis for contrast did not yield a significant result.

In summary, although we did find a weak effect of attention on spike-gamma phase relationship, based on the variety of issues that we have discussed, gamma PC is at best expected to play a minor role in the coding the stimulus contrast in V1.

# AUTHOR CONTRIBUTIONS

AD and SR conceived the idea of research. SR collected data; AD and SR analyzed data. AD and SR wrote the paper.

# FUNDING

This work was supported by Wellcome Trust/DBT India Alliance (Grant 500145/Z/09/Z; intermediate Fellowship to SR), Tata Trusts, Department of Biotechnology-Indian Institute of Science (DBT-IISc) Partnership Programme and a Junior Research Fellowship awarded by IISc from Ministry of Human Resource Development, Government of India (to AD).

# ACKNOWLEDGMENTS

We thank Dr. John Maunsell for his help in experimental design and data collection and Steven Sleboda and Vivian Imamura for technical support.

# REFERENCES

potentials in the rat hippocampus. J. Neurosci. 32, 11798–11811. doi: 10.1523/JNEUROSCI.0656-12.2012

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Das and Ray. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Gender-Dependent Changes in Time Production Following Quadrato Motor Training in Dyslexic and Normal Readers

Tal Dotan Ben-Soussan1 \* and Joseph Glicksohn2,3

1Research Institute for Neuroscience, Education and Didactics, Patrizio Paoletti Foundation for Development and Communication, Assisi, Italy, 2Department of Criminology, Bar-Ilan University, Ramat Gan, Israel, 3The Leslie and Susan Gonda (Goldschmied) Multidisciplinary Brain Research Center, Bar-Ilan University, Ramat Gan, Israel

Time estimation is an important component of the ability to organize and plan sequences of actions as well as cognitive functions, both of which are known to be altered in dyslexia. While attention deficits are accompanied by short Time Productions (TPs), expert meditators have been reported to produce longer durations, and this seems to be related to their increased attentional resources. In the current study, we examined the effects of a month of Quadrato Motor Training (QMT), which is a structured sensorimotor training program that involves sequencing of motor responses based on verbal commands, on TP using a pre-post design. QMT has previously been found to enhance attention and EEG oscillatory activity, especially within the alpha range. For the current study, 29 adult Hebrew readers were recruited, of whom 10 dyslexic participants performed the QMT. The normal readers were randomly assigned to QMT (n = 9) or Verbal Training (VT, identical cognitive training with no overt motor component, and only verbal response, n = 10). Our results demonstrate that in contrast to the controls, longer TP in females was found following 1 month of intensive QMT in the dyslexic group, while the opposite trend occurred in control females. We suggest that this longer TP in the female dyslexics is related to their enhanced attention resulting from QMT. The current findings suggest that the combination of motor and mindful training, embedded in QMT, has a differential effect depending on gender and whether one is dyslexic or not. These results have implications for educational and contemplative neuroscience, emphasizing the connection between specifically-structured motor training, time estimation and attention.

Keywords: time production, quadrato motor training, dyslexia, gender difference, time and motion studies

# INTRODUCTION

Timing deficits in dyslexia include those concerned with time estimation (Nicolson et al., 1995; Ramus et al., 2003; Hölzel et al., 2011), rhythm tapping (Wolff et al., 1990; Wolff, 2002), detecting complex timing patterns (Kujala et al., 2007), rapid temporal processing (Tallal et al., 1993), auditory temporal sensitivity (Witton et al., 1998) and visual motion detection

#### Edited by:

Daya Shankar Gupta, Camden County College, United States

#### Reviewed by:

Marc Wittmann, Institut für Grenzgebiete der Psychologie und Psychohygiene, Germany Matthew S. Matell, Villanova University, United States

#### \*Correspondence:

Tal Dotan Ben-Soussan research@ fondazionepatriziopaoletti.org

Received: 26 March 2018 Accepted: 08 August 2018 Published: 29 August 2018

#### Citation:

Ben-Soussan TD and Glicksohn J (2018) Gender-Dependent Changes in Time Production Following Quadrato Motor Training in Dyslexic and Normal Readers. Front. Comput. Neurosci. 12:71. doi: 10.3389/fncom.2018.00071 (Talcott et al., 2000). Consequently, dyslexia-related timing deficits have at different times been hypothesized as underlying dyslectics' visual and auditory perception problems, motor coordination problems and fluency and automatization problems (Nicolson et al., 2001), all of which have been proposed as adversely affecting the development of language and literacy skills (Overy et al., 2003). The full extent of the timing deficits is yet to be established, but suggests the need for further investigation.

A critical factor here seems to be the cerebellum, which is implicated in both timing functions (Ivry and Hazeltine, 1992; Ivry, 2000, p. 187; Ivry et al., 2002; Rubia, 2006) and in dyslexia (Reynolds et al., 2003; Kujala et al., 2007; Reynolds and Nicolson, 2007; Ben-Soussan et al., 2014a). Thus, cerebellar deficits are thought to affect articulation and working memory, due to deficits in timing which interfere with automatization of learning (Thomas and Karmiloff-Smith, 2002; Overy et al., 2003; Ram-Tsur et al., 2013).

In fact, cerebellar oscillatory function and its role in motor acquisition and timing have long been acknowledged (Andres et al., 1999; Swinnen, 2002; De Zeeuw et al., 2011). Given that impaired motor skills are often observed in dyslexics, some researchers have attributed dyslexics' cognitive and motor deficiencies to abnormal development and functioning of the cerebellum (Nicolson et al., 1995, 2001). These findings have led to the claim that the role of the cerebellum is not limited to regulating the timing, rate, force, rhythm and accuracy of movements, but also to the speed, capacity, consistency and appropriateness of cognitive processes (Schmahmann, 2004; Hölzel et al., 2011; Buckner, 2013).

Ben-Soussan et al. (2015) have recently presented a general model tying cerebellar function to cognitive improvement, by means of a particular form of motor training, which might be viewed as meditation-in-action, namely Quadrato Motor Training (QMT). QMT has been found to increase creativity, reflectivity and spatial cognition (Ben-Soussan et al., 2013, 2014b, 2015), as well as to increase neuronal synchronization and connectivity, especially within the alpha (8–12 Hz) range (Lasaponara et al., 2017). In addition, a month of daily QMT was found to improve reading and increase cerebellar alpha oscillations in dyslexic adults (Ben-Soussan et al., 2014a). Recently, QMT was further found to increase fractional anisotropy (FA) in tracts related to sensorimotor and cognitive functions and mindfulness, including the corticospinal tracts, anterior thalamic radiations and uncinate fasciculi, as well as in the left inferior frontooccipital, superior and inferior longitudinal fasciculi (Piervincenzi et al., 2017), reflecting better white matter integrity as a result of greater intravoxel coherence of fiber orientation, axon density and diameter and/or myelination (Beaulieu et al., 1996; Sen and Basser, 2005; Caminiti et al., 2013).

Let us contrast a hypothesized QMT-based improvement in the functioning of dyslexic adults to that found for another movement-based form of training, which has also reported an improvement in reading fluency for dyslectic participants. Reynolds et al. (2003) study was conducted on dyslectic pupils, reporting improvements in balance, dexterity and reading. As they subsequently reported in a follow-up study (Reynolds and Nicolson, 2007), both dyslexic and non-dyslexic children benefited from the training, while alternative hypotheses raised by critics of the original study (e.g., based on potential artifacts, such as a Hawthorne effect) could be ruled out. In these studies, training was comprised of a home-based exercise program. In the study we report below, our QMT is also a home-based motorexercise program. We do, however, also employ a home-based verbal-exercise program, as a suitable control.