#### Abstract

Stochastic processes that exhibit cross-frequency coupling (CFC) are introduced. The ability of these processes to model observed CFC in neural recordings is investigated by comparison with published spectra. One of the proposed models, based on multiplying a pulsatile function of a low-frequency oscillation () with an unobserved and high-frequency component, yields a process with a spectrum that is consistent with observation. Other models, such as those employing a biphasic pulsatile function of a low-frequency oscillation, are demonstrated to be less suitable. We introduce the full stochastic process time series model as a summation of three component weak-sense stationary (WSS) processes, namely, , , and , with a noise process. The process is constructed as a product of a latent and unobserved high-frequency process with a function of the lagged, low-frequency oscillatory component (). After demonstrating that the model process is WSS, an appropriate method of simulation is introduced based upon the WSS property. This work may be of interest to researchers seeking to connect inhibitory and excitatory dynamics directly to observation in a model that accounts for known temporal dependence or to researchers seeking to examine what can occur in a multiplicative time-domain CFC mechanism.

#### 1. Introduction

Cross frequency coupling (CFC) is a statistical relation between the phase or amplitude of a low frequency and the phase or amplitude of a high frequency. In this work, focus is placed upon phase-amplitude CFC, which can be thought of as the correlation of the amplitude of a relatively high-frequency oscillation () with the phase of a lower frequency oscillation ().

The rationale behind this proposal is based on the experimental observations that (i) relatively slow frequency oscillations tend to be coordinated over large regions of neural tissue, unlike higher frequency oscillations [1, 2], and (ii) oscillatory activity reflects changes in the excitability of neural tissue [3]. Hence, a lower frequency oscillation may provide time intervals in which high-frequency activity may occur and may consequently coordinate high-frequency oscillations that are spatially disparate.

Recently, the theoretical interest in phase-amplitude CFC has been reinforced. Recent observations of the phenomenon have occurred in varied species [4–6], brain regions [7–9], and states of vigilance [10]. Experimental studies have shown that the nature of phase-amplitude CFC can be altered by predictive cues and attentional demands; the phase of the low-frequency oscillation can be reset such that a stimulus of attentional interest arrives at the phase of maximal excitability [11]. Furthermore, phase-amplitude CFC has been implicated in learning and memory [6, 12], and the dynamics of phase-amplitude CFC have been shown to change over the course of a cognitive task [8].

In this work, in a fashion akin to that classically employed for parametrically modeling the spectra of time-series, stochastic time-series models are constructed which exhibit cross frequency coupling (see, e.g., [13–15]). Through simulation it is shown that these models can produce time-series that exhibit CFC similar to that exhibited in neural recordings. Some mathematical properties of these models are given, and some consequences are discussed.

#### 2. Methodology

Stochastic processes exhibiting cross frequency coupling are introduced. The ability for these processes to model observed cross frequency coupling in neural recordings is investigated and mathematical properties of the new models are given. Investigations are conducted through the use of mathematical analysis and simulation.

##### 2.1. Model Specification

The stochastic, or random, processes introduced in this work model CFC phenomena in the following way. For each recorded measurement a random variable is introduced; the measurement is modeled as a realization of this random variable. The collection of random variables comprises a discrete-time random process. The observed time-series is modeled as the corresponding collection of realizations of each of the random variables in the random, or stochastic, process. This is the standard setup in classical time-series analysis [13] and it includes the independent and identically distributed (IID) random sample as a special case.

Each of the random processes is characterized by all of the possible moments between the random variables; in this work consideration will be restricted to random processes whose joint distributions are Gaussian. In this situation, specification of the first two joint moments completely specifies the model.

Without restriction upon the time dependence of pair-wise correlation, the number of unique pair-wise correlations increases quadratically with every new observation (in a single trial). This is a more challenging regime to perform inference than is typically considered, as in this case the number of unknowns is growing rapidly with increasing observation length. Contemporary work deals with this issue by recording many trials, or by using models with other restrictions.

Here, as is often customary, the introduced models are weak-sense stationary (WSS) random processes [13–15]. The weak-sense stationary property implies that the correlation between any pair of random variables in the process does not depend upon absolute time, but, rather, only upon the difference in the times associated with the pair. It also implies that the mean of each of the random variables in the process is equal; thus the process mean is also independent of time. Gaussian weak-sense stationary random processes are amongst the simplest of time-series models and the number of pair-wise correlations for these processes is of the same order as the number of measurements.

In this work, the random processes exhibiting CFC are constructed from component WSS processes modeling -rhythm, -rhythm, background activity, and sensor noise. The mean of these processes is taken to be zero, consistent with randomly observed neural phenomena (nonevoked). Based upon observed spectra, the autocorrelation of the , , and noise components is specified in the Fourier domain. Based upon the discrete-time analog of the Wiener-Khintchine theorem, the autocorrelation of each of these component processes is obtained by inverse discrete-time Fourier transforming the specified spectra [13].

###### 2.1.1. The Component Processes

Let be the integer-valued time-index of the length WSS zero-mean random process with autocovariance sequence :Here denotes the expected value of the random variable . Similarly, specify WSS zero-mean processes for the rhythm and noise, , components of the model. That is, Further, specify the and components as uncorrelated. Because both and are also zero-mean, it follows that is equal to zero. These components are jointly Gaussian, uncorrelated, and hence they are independent. The components and are linked to model CFC. This linking and its consequence are discussed in Section 2.1.2. It remains to specify the autocovariance sequences , , and . As described, this is accomplished by specificying their respective spectra, , , and , and using the example relation obtained by applying the discrete-time analog to the Wiener-Khintchine theorem [13]: Here is the Nyquist frequency, equal to (in Hz), specified in terms of the sample period (in ). Figures 1 and 2 depict the specified model autocovariance sequences and spectra for the and components (resp.). The component is further detailed in Section 2.1.2.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

###### 2.1.2. Specification of Exhibiting Cross Frequency Coupling

In this work, CFC is specified by equating to a function , relating to an unobserved, latent random process . Let be a zero-mean, WSS, random process possessing a spectrum , such that is near zero outside of the frequency interval associated with a -rhythm. Candidate frequency intervals include low- (30 Hz–50 Hz) and high- (100 Hz–140 Hz). (These ranges can differ depending on the lab, species, and neural area. In [4], the low- frequency interval is specified to be 30 Hz–50 Hz, and the high- frequency interval is specified to be 80 Hz–150 Hz.) Consider where is a vector of random variables belonging to the -process. The th element of the random vector is equal to . The function can be chosen to model various types of CFC. Together the process and the function control the dependence of upon . The relation (4) is general and can be used to model many types of CFC. For example, when exhibits a phase shifted and sinusoidal variance (see Figure 3). When exhibits a punctate interval of increased variance centered upon the prefered phase , (for the equal to case, with , , , , and , see Figure 3). (To be pulsatile the Fourier transform of must have peaks spaced by . Here this is accomplished by summing taken to higher powers). When is squared in (6), and biphasic coupling is modeled (see Figure 3, , and as specified in Figure 3). Note that doubling the power of in (7) doubles the harmonic frequency spacing present in the specified in (6). In the Appendix it is shown that the inverse tranform of a sequence with four harmonics in the frequency domain has a period equal to the inverse of the harmonic spacing.

###### 2.1.3. The Complete Time-Series Model

The model random process generating the observed time-series is specified as the sum of the , , and components. Consider

##### 2.2. Properties of the Model

Mathematical properties resulting from the specification (8) are explored by calculation.

###### 2.2.1. The Computation of

By definition the mean, , of the component is

###### 2.2.2. The Computation of

By definition, the process mean, , evaluated at time-index is Since all three components , , and are mean-zero processes. Thus, the model process is also mean-zero.

###### 2.2.3. The Computation of

By construction , , and are WSS random processes. To establish as a weak-sense stationary process, it is necessary to compute the autocovariance function . Begin with the definition for the general autocovariance function evaluated at the time indices and . One has From an application of the Isserlis formula [14, 16], the following relation results:which is independent of .

###### 2.2.4. The Spectrum, , of

The discrete-time Fourier transform of is equal to . From an application of the discrete-time convolution theorem, the following expression relating and to is obtained:

###### 2.2.5. Cross Frequency Coupling and the Frequency of

To explore the effect of center frequency upon cross frequency coupling, in simulation, two center frequencies are used. These frequencies are 6 Hz (Figure 4) and 15 Hz (Figure 5). The effect of making this change can be seen by comparing the spectra plotted in Figure 4 with those plotted in Figure 5.

###### 2.2.6. The Effect of Bandwidth on Pulse Duration

Since is approximately nonzero only on a 40 Hz interval, the absolute-autocovariance, , is approximately nonzero for lags of seconds. Thus, pulses will tend to have about a 25-millisecond duration.

###### 2.2.7. Simulated Time-Series: Sample Path Generation

To obtain a length sample-path, that is, a truncated realization, of any one of the random processes , , or described in this work, a draw from IID standard Gaussian random variables is made (mean-zero, unit variance). These Gaussian random variables can be collected into an -dimensional vector . The difference between and the desired random vector , with covariance matrix , is since Thus, if is known, then a realization of can be computed by multiplying a realization of by . It remains to specify . For a weak-sense stationary process, the covariance between random variables depends only upon their separation in time-index. Thus, is Toeplitz with constant diagonals and is completely specified by the autocovariance function. For example, for the process , This method of computing simulated realizations is used for all processes in this work except for , and hence, for . In the former case, a realization of is computed from realizations of and . Similarly, a realization of is computed by summing realizations of , , and .

##### 2.3. Model Assessment

The appropriateness of as a model of actual recordings exhibiting CFC is assessed, somewhat crudely, by comparing visually the sample paths of generated with the three different functions relating and to and by comparing the spectrum of and to that exhibited [5] in actual neural recordings.

#### 3. Results

##### 3.1. Simulated Time-Series Exhibiting CFC

Figure 3 depicts sample paths of constructed, using (5), (6), and (7). Here (5) corresponds to the label “Sinusoidal ,” (6) corresponds to the label “Pulsatile ,” and (7) corresponds to the label “Biphasic .” These time-series approximate high-pass filtered neural recordings exhibiting CFC.

##### 3.2. The Autocovariance Sequence

Having established that is mean-zero, the autocovariance of is equal to Due to independence and the fact that is equal to zero, the following is considered: but Similarly, is also zero. Then

##### 3.3. Is Weak-Sense Stationary (WSS)

The three components , , and are each mean-zero, and hence is also mean-zero. Since the autocovariance is completely specified in (20) in terms of , , and , each of which does not depend upon absolute time, also does not depend upon absolute time. Since the first two moments of the process are independent of absolute time, is WSS.

##### 3.4. The Spectrum of

Applying (3) to (20) results in the following relation:

##### 3.5. Time-Series Models of Cross Frequency Coupling (CFC)

Using the three , coupling functions (sinusoidal, pulsatile, and biphasic pulsatile) used to specify and the simulation method discussed in Section 2.2.7, 5000 realizations of are simulated. From each of these realizations, an estimate of is computed. These estimates are averaged to obtain the three curves depicted in Figure 4. The convolution between and produces prominent side-lobes in the spectra of for the sinusoidal and biphasic cases. For the pulsatile case the spectrum of appears similar to those exhibited by actual neural recordings [5].

##### 3.6. for a Frequency-Shifted

The simulation presented in Section 3.5 is repeated for a modified . Specifically, the frequency about which the component is centered is shifted from 6 Hz to 15 Hz. The resulting spectra associated with the sinusoidal, pulsatile, and biphasic - coupling are depicted in Figure 5. Sample paths of the associated are shown in Figure 6.

#### 4. Discussion

The stochastic models considered in this work are weak-sense stationary. Phenomena without the WSS property, such as stimulus transients and time varying changes in autocovariance, are not captured by these models. While the topic is left for future study, one expects a convolution similar to that producing to appear in the nonstationary analogs of the time-series models considered in this work.

In [17] the probability density function for a phase-amplitude coupling estimator is provided under the assumption that phase-amplitude CFC is absent. The parametric regression based phase-amplitude coupling estimators [18, 19] are based upon conditional expectations which may not completely reflect the statistical properties of recordings that exhibit phase-amplitude coupling. The time-series models explicitly proposed in this work account for realistic temporal dependency and exhibit phase-amplitude CFC. As in [17], sampling properties of proposed estimators can be assessed based on a time-series model.

Owing to the convolution integral present in (13) and (21), the time-series introduced in this work exhibit, in general, multiple peaks about the 120 Hz center frequency of . This effect is minimized for CFC developed through the product of the latent process with the pulsatile function specified in (6). Of the considered models, it is this latter model which most closely matches observed CFC phenomena. The pulsatile function can be loosely thought of as imposing a window over each cycle of the slow oscillation during which the latent and relatively high-frequency process is amplified. In the neural context, one might infer the existence of a similar window defining intervals of time when neural activity is more easily generated, consistent with the notion that oscillatory activity represents cyclic changes in excitability [3].

The excitability described in [3] is considered to be a function of the relative levels of synaptic inhibition and excitation. Pulsatile functions that more closely match the shape and time-course of synaptic currents may more closely match the observed features of phase-amplitude CFC in neural data. By comparing the simulated data resulting from the shape and time course of a simulated pulsatile function to actual neural data, insight may be gleaned regarding the nature of synaptic currents underlying phase-amplitude CFC. This modeling process, beginning with a stochastic time-series model of observed phenomena, may also provide insight into the other types of CFC, namely, amplitude-amplitude and phase-phase cross frequency coupling.

The process of time-series model simulation and the comparison of simulation results to results obtained from neural data can be facilitated with the use of statistical methods. The quantities in the proposed time-series models can be estimated, confidence intervals can be provided, model selection can be performed, and across-condition comparisons and hypothesis tests can be developed. Work refining these models and developing associated statistical methodology is underway.

The time-series models introduced in this work are, to the best of the authors’ knowledge, the first reported time-series models capable of generating phase-amplitude cross frequency coupling. They may provide a basis for a more accurate and statistically principled assessment of cross frequency coupling and underlying synaptic activity that accounts for the temporal dependencies present in actual recordings.

Finally, the method of simulating sample-paths of weak-sense stationary processes introduced in Section 2.2.7 seems to be novel to the neuroscience community. Its use may facilitate increased accuracy when numerically investigating the behavior of statistics computed from neural recordings.

#### 5. Conclusion

Full appreciation of the effect of different time-domain manifestations of cross frequency coupling may provide insight into neural processes. In particular, prominent spectral features are produced by somewhat mild changes to the nature of a time-series model exhibiting phase-amplitude cross frequency coupling. These changes may provide the basis for refined statistical methodology capable of shedding light into the nature of the inhibitory and excitatory synaptic currents supporting phase-amplitude cross frequency coupling.

#### Appendix

#### DTFT of a Periodic Sequence

Let be a sequence with discrete-time Fourier transform (DTFT) . Then Further, let be nonzero only for , . Then Because for each , is periodic with period , is also periodic with period .

#### Conflict of Interests

The authors declare that they have no competing interests.

#### Acknowledgment

The authors are supported by NSF Grant DMS-1042134.