To study interaction between ﬂicker and background activity we model ϕnas ϕntð Þ¼ ϕnoise n tð Þ þϕD ntð Þ; ð12Þ where ϕnnoiseis Gaussian white noise with mean ϕn(0)and standard de-viation σn, and ϕn Dis a periodic drive with zero mean. We use two forms for ϕn Din the Results section. For comparison with Fig. 1 (a) we use a square wave drive given by ϕD ntð Þ¼Φnsgn cos 2 πf Dt ð Þ 1/2/C138 ; ð13Þ where Φnis the drive amplitude and sgnindicates the sign function. This square wave consists of a sum of sinusoids, such that the drive contains power at all odd multiples of f D; system responses at even multiples are necessarily nonlinear. To study the simpler case of a sin-gle sinusoid we use the form ϕD ntð Þ¼ Φncos 2 πf Dt ð Þ ; ð14Þ for which all higher harmonics in the model output are necessarily generated nonlinearly. Our use of the square wave drive Eq. (13) matches the photic drive used in Herrmann (2001), but neglects possible nonlinearities in the retina that might have modi ﬁed the latter drive before it ar-rived at the corticothalamic system. We thus use a square wave so as to introduce the fewest additional assumptions, but note that it would be straightforward to implement other more complicated forms for the input. For example, the model could be explicitly coupled to models of retinal dynamics, which would enable it to compare and evaluate putative nonlinearities in the earliest parts of the visual stream. However, these topics are beyond the scope of this paper. To solve the model equations numerically we use a standard ex-plicit forward step integration scheme with a timestep of 0. 1 ms (Press et al., 1992 ), of the type used by Robinson et al. (2002), for ex-ample. We have checked that the results are numerically stable and that the timestep is suf ﬁciently short that they do not change signif-icantly when it is further reduced. Time series analysis To characterize the frequency content of the model's response to a range of stimulus parameters, we calculate power spectra from model time series. We estimate spectra from windowed Fourier transforms via Welch's method (as implemented in MATLAB R2008b), using 10 s windows overlapping by 50%. To display sharp peaks more clearly, we additionally smooth spectra with a 0. 4 Hz moving window. To elucidate the role of nonlinear frequency coupling in generating peaks in the model EEG spectra we calculate bicoherence from model time series ( Kim and Powers, 1979 ). Bicoherence is a nonlinear coher-ence measure that is sensitive to phase locking between Fourier compo-nents at different frequencies, and is derived from the bispectrum, a higher-order spectrum that generalizes the usual power spectrum (Kim and Powers, 1979 ). The bicoherence b2(ω1,ω2) between frequen-ciesω1andω2of a signal with Fourier transform X(ω)i sg i v e nb y( Kim and Powers, 1979 ) b2ω1;ω2 ð Þ ¼Xω1ð Þ Xω2ð Þ X/C3ω1þω2 ð Þ jj2 Xω1ð Þ Xω2ð Þ jj2Xω1þω2 ð Þjj2; ð15Þ where X*is the complex conjugate of X. W ee s t i m a t et h e X(ω) terms by averaging FFTs of windowed time series (16384-point Hanning win-dow, mean subtracted from each). To characterize nonlinear dynamics at large drive amplitudes we estimate the largest Lyapunov exponent (LLE) ( Parlitz, 1998 ) frommodel time series using the software package Open TSTOOL v1. 2 (Merkwirth et al., 2009 ). A positive LLE implies that nearby trajecto-ries diverge exponentially and is indicative of chaos; periodic orbits have an LLE of zero. We use the method of time-delay embedding on long (80 s) time series for ϕesampled at 1000 Hz, with the mini-mum embedding dimension chosen using the method of Cao (1997) with three nearest neighbors and 1000 reference points, and thedelay chosen using the ﬁrst minimum of the auto-mutual information function (see Open TSTOOL documentation for further details). A ben-eﬁt of using this method for LLE estimation is that it is also applicable to experimental time series (i. e., does not require access to the full underlying dynamical system). To analyze phase-locking regimes as a function of drive amplitude and frequency, we calculate for each point in parameter space the ratio of the drive's period to that of time series simulated in the ab-sence of noise. We use the sequence of turning points in ϕ eto calcu-late both the time series period (thus identifying 1: Nphase-locking regimes) and the number of oscillations per period (identifying N:1 regimes, where the drive and ϕehave the same period). Results We begin in the Results section by driving the system with a square wave photic input to compare the model output to the data of Herrmann (2001) in Fig. 1 (a). We analyze the driven dynamics in more detail using sinusoidal input in the Sinusoidal drive section. Square wave drive Fig. 1 (b) shows the model response spectrum vs. f Dfor a square wave drive given by Eq. (13), with each spectrum multiplied by ffor comparison with Fig. 1 (a). Key features are peaks at f D(1:1) and its har-monics ( N:1), in good agreement with Fig. 1 (a) (Herrmann, 2001; Rager and Singer, 1998 ). Most signi ﬁcantly, entrainment to f Dis evidenced by suppression of the background alpha activity fα≈11 Hz in favor of f Dfor f D≲15 Hz, consistent with Fig. 1 (a) where no strong alpha peaks are seen off the 1:1 diagonal. The alpha peak is suppressed without any changes to cortical or thalamic parameters; only the drive frequency changes. Further striking agreement is shown by 1:2 subharmonic en-trainment for f D≈15-24 Hz, where peaks follow f=f D/2. Suppression of background activity is necessarily nonlinear; in the linear regime the drive and background superpose without interacting. The suppres-sion observed experimentally cannot be attributed to averaging over trials, or else the alpha peak would not be observed for any ﬁxed drive frequency contradicting Fig. 1 (a). The extent to which the alpha activity is suppressed by the drive in the experimental data is not as clear as in the model; peaks in Fig. 1 (b) are narrower than in Fig. 1 (a) due to the FFT used here yielding better resolution than the autoregressive spec-trum in Fig. 1 (a). The experimental data also appears “lumpier ”than the model data due to being sampled on only a coarse grid of frequen-cies; similarly coarse sampling in the model would yield an even more similar appearance. Additional responses in the model at f=3f D/2 and f=5f D/2 are not seen in Fig. 1 (a), likely due to Fig. 1 (a) being an average over ten subjects, thus further blurring peaks via variability between subjects. Sinusoidal drive In this subsection we drive the system with a single sinusoid to avoid ambiguity between spectral peaks that would be present in a nonsinusoidal drive and peaks generated by nonlinear processes. Fig. 3 (a) shows that a sinusoidal drive given by Eq. (14) reproduces the main features of the EEG response to a square wave drive, demon-strating that entrainment and harmonic generation are not speci ﬁct o square wave input. Here and in Fig. 3 (b) additional features are re-vealed by plotting a wider dynamic range than in Fig. 1 (a) (here we1950 J. A. Roberts, P. A. Robinson / Neuro Image 62 (2012) 1947-1955