Wilson and Cowan, 1973 ), neural ﬁeld models ( Baker and Cowan, 2009; Bressloff et al., 2002; Coombes, 2005; Hutt and Atay, 2005; Jirsa, 2009; Jirsa and Haken, 1996; Pinotsis and F riston, 2011; Pinotsis et al., 2012; Robinson et al., 1997; Taylor and Baier, 2011; Wright and Liley, 1996 ) and their limiting case of neural mass models ( Cosandier-Rimélé et al., 2008; David and Friston, 2003; David et al., 2006; Goodfellow et al., 2011; Jansen and Rit, 1995; Lopes da Silva et al., 2003; Moran et al.,2009; Wendling et al., 2000 ) have successfully de s c r i b e dab r o a dr a n g e of healthy and pathological brain dynamics at the spatial and temporal scales probed by EEGs and magnetoencephalograms. Neural ﬁeld models have been notably successful in describing spatiotemporal dynamics in vi-sual cortex ( Baker and Cowan, 2009; Bre ssloff et al., 2002; Coombes, 2005; Ermentrout and Cowan, 1979 ), in particular explaining pattern for-mation during hallucinations. Visual evoked responses to individually-ﬂashed stimuli have also been studied using neural mass models ( David et al., 2006; Jansen and Rit, 1995 ), and ( Spiegler et al., 2011 ) very recently applied the ( Jansen and Rit, 1995 ) model to study entrainment and other nonlinear features of photic-driven activity. Both classes of model have explained a range of seizure phenomena, relating this pathological activ-ity to healthy dynamics generated by the same models ( Breakspear et al., 2006; Cosandier-Rimélé et al., 2008; Goodfellow et al., 2011; Kim et al., 2009; Lopes da Silva et al., 2003; Nevado-Holgado et al., 2012; Roberts and Robinson, 2008; Robinson et al., 2002; Taylor and Baier, 2011; Wendling et al., 2000 ). In this paper we explain experimentally-observed nonlinear inter-actions between brain activity and ﬂicker by incorporating periodic driving into a corticothalamic model ( Roberts, 2010; Robinson et al., 1997, 2002 ). The results reproduce multiple features of published data and constrain the physiologically allowable states. We further predict various other nonlinear phenomena whose existence can be tested in future experiments, including period doubling, bistable phase-locking, hysteresis, wave-wave interactions, and chaos, and show how to treat entrainment and seizures in the same framework. This extends neural ﬁeld modeling into a new regime and enables new experimental tests of the theory. We focus on EEG responses to sinusoidally modulated light, which are termed steady-state visual evoked potentials (SSVEPs) ( Nunez, 1995; Regan, 1989 ), but our approach is general and could be applied to any periodic stimulus. The SSVEP is dominated by synchronous corti-cal activity concentrated at the drive frequency f Dand its harmonics. Ex-periments have measured the response over a wide range of f Din cats (Rager and Singer, 1998 ) and humans ( Herrmann, 2001 ), showing clear spectral peaks, as seen in Fig. 1 (a) for spatially uniform square wave modulated stimuli. Here, as in the original ﬁgure in Herrmann (2001), the spectrum for each f Dis multiplied by fto enhance peaks at higher frequencies, and the peaks appear as discrete dots due to being sampled on a coarse grid. These peaks involve both responses at f D (which can be linear or nonlinear) and inherently nonlinear effects. The clearest nonlinear effect is the generation of N:1 harmonics (Nresponse oscillations phase-locked to each drive oscillation). Since the square wave input signal used here contains only odd har-monics, the even harmonics in the output (in particular the prominent 2:1 harmonic) must be generated nonlinearly; power at the odd har-monic frequencies likely contains both linear and nonlinear compo-nents. The position of this nonlinearity in the visual pathway is not yet well established, and could involve effects in retinal, thalamic, and/or cortical neuronal populations. Suppression of background alpha activity over a range of f Dis another signiﬁcant nonlinear effect seen in Fig. 1 (a), where the dominant re-sponse tracks f Dorf D/2, while the dominant frequency is unchanged out-side this range. (These oscillations are network-level collective modes, whose frequencies are not individual neural ﬁring rates in general. ) This is an example of entrainment, characterized by a reduction of ongoing ac-tivity in favor of activity phase-locked to the drive, at frequencies harmon-ically related to the drive frequency. Fig. 1 (a) thus shows entrainment to f D(1:1 phase locking) and subharmonic entrainment to f D/2 (1:2 locking). Quantitative modeling, such as that used in the present work, will be critical in analyzing the data to distinguish between the above possibili-ties. In this context, it is particularly important to use established models, rather than introducing a new one ad hoc for each new experiment. Established models enable new phenomena to be incorporated into wider, uni ﬁed frameworks, thereby relating them to other measures and situations, and permitting paramet ers to be constrained a priori, rath-er than being treated as free. The wide use of the Hodgkin-Huxley model is perhaps the best known example of this approach in neuroscience. This paper is organized as follows. The Methods section outlines the nonlinear corticothalamic model. In the Results section the periodically driven dynamics of this model system are shown to reproduce the main features of the experimental SSVEP spectra in Fig. 1 ( a )a sw e l la s predicting rich nonlinear dynamics whose existence and properties can be explored in future experiments. Methods The Corticothalamic model section presents the nonlinear cor-ticothalamic model, the dynamics of which will be analyzed using the time series methods described in the Time series analysis section. Fig. 1. Spectral response for f D=1-50 Hz. Responses to the Nth harmonic of f Dare la-beled N:1; subharmonic entrainment is labeled 1:2. Darker shading denotes higher power, and spectra have been multiplied by fto enhance power at high frequencies as in the original Fig. 5 of Herrmann (2001). (a) Adult human EEG. Figure adapted from Fig. 5 of Herrmann (2001) after cropping, changing aspect ratio to 1:1, and put-ting independent variable on the horizontal axis. (b) Model with ( νee,νei,νes,νse,νsr, νsn,νre,νrs)=(1. 7,-2. 8, 0. 19, 0. 84,-0. 62, 1. 0, 0. 28, 3. 0) m V s, ( α,β,αsr,βsr)=(80, 800, 10, 60)s-1,bϕn0>=18 s-1,σn=1s-1,Φn=2s-1; other parameters as in Breakspear et al. (2006). Spectra are smoothed with a 0. 4 Hz moving window to display sharp peaks more clearly. 1948 J. A. Roberts, P. A. Robinson / Neuro Image 62 (2012) 1947-1955