Fig. 1 (b) is for noisy perturbations of a ﬁxed point. Bifurcations of the ﬁxed point itself (under changes in parameters other than the drive) affect entrainment. Numerically we ﬁnd that entrainment requires the system to lie near a linear stability boundary, in accord with evi-dence that the brain operates near marginal stability ( Breakspear et al., 2006; Robinson et al., 1997; Stam, 2005; Stam et al., 1999 ), and also near a nonlinear stability boundary, as seen in Fig. 3 (e). Entrain-ment like that in Fig. 1 also requires that no other attractor coexist near the ﬁxed point; otherwise the spectrum would change marked-ly. For a bistable limit cycle associated with seizure dynamics, such as near the subcritical 10 Hz bifurcation in the model ( Breakspear et al., 2006 ), a sinusoidal drive near resonance perturbs the system away from the stable state onto a large-amplitude attractor, even for rela-tively weak stimuli ( Kim et al., 2009 ); a typical case in Fig. 5 shows a dramatic increase in power at all frequencies for drives near fα≈9 Hz and 2 fα, with a waveform similar to the 10 Hz seizure in Robinson et al. (2002) and Breakspear et al. (2006) ; such behavior is characteristic of photosensitive epilepsy ( Kim et al., 2009; Parra et al., 2005 ). The seizure frequency f≈10 Hz is little affected by the drive, which hardly perturbs the large-amplitude attractor. Numerical exploration reveals that entrainment of fαis favored for parameters yielding decreases (2-to 10-fold relative to typical alert resting values) in mean ﬁring rates ϕe,ϕi,a n dϕs,w i t h ϕrincreased (pos-sibly substantially) or unchanged. Such decreases are qualitatively con-sistent with decreased metabolic load in these structures during entrainment ( Parkes et al., 2004 ) (assuming that background activity is analogous to the aperiodic stimulation used there for comparison), suggesting that subnetworks activated by stimuli likely differ from those underlying resting activity, agreeing with recent experiments (Birca et al., 2006; Kerr et al., 2008 ). A relative increase in thalamic inhi-bition is plausibly a form of automatic gain control ( Schwartz and Simoncelli, 2001 ), attenuating strong sensory stimuli before they are re-layed to the cortex. Discussion We have studied driven nonlinear brain dynamics in a physiologically-based model of the cortex and thalamus, have compared its predictions with experiments from the literature, and have made a series of newpredictions that can be tested in future experimental work. The main re-sults are as follows: (i) Our theory reproduces the key features of entrainment of the alpha rhythm to periodic stimuli, including entrainment tosubharmonics of the drive, showing extensive agreement with experiment ( Herrmann, 2001 ). Further good agreement is shown by the presence of nonlinearly-generated harmonics of the drive. The comparison in Fig. 1 demonstrates that the model agrees with experiment over wide ranges of drive and response frequencies. (ii) The model predicts additional nonlinear dynamics to be found in future experiments. We predict bistability between different entrained states having the same drive amplitude and drive frequency, which can be tested experimentally by looking for hysteresis when slowly ramping the drive frequency up and down. Spontaneous transitions between bistable states during a long recording session might also contribute to intertrial var-iability, a possibility that could be tested experimentally. (iii) At large drive amplitudes we predict period doubling of phase-locked states leading to quasiperiodic and chaotic dynamics. The 1:1 state bifurcates to a 1: Nregime for minimal drive am-plitude when f D≈Nfα, consistent with experiment ( Herrmann, 2001 ). (iv) We predict nonlinear sum and difference frequency generation in cases where background activity is sharply peaked but not entrained by the drive. This has been observed in studies using two stimulus frequencies simultaneously ( Regan, 1989 ). (v) The model predicts that when a stable limit cycle coexists with the stable ﬁxed point corresponding to the resting state, peri-odic stimuli near the alpha frequency and its harmonics can drive the system into a seizure state. This plausibly explains seizures induced by high amplitude ﬂashing light ( Kim et al., 2009; Parra et al., 2005 ), and is open to further test. These ﬁndings provide new veri ﬁcations of the model, and of neu-ralﬁeld theory more generally, complementary to previous studies. Since the same model has previously successfully described many lin-ear and nonlinear dynamics ( Breakspear et al., 2006; Robinson et al., 1997, 2002, 2004 ), this paper uni ﬁes these phenomena with nonlinear SSVEPs within the same framework using compatible pa-rameters. Moreover, this work opens new interdisciplinary research avenues by predicting a rich variety of nonlinear dynamics that should be experimentally veri ﬁable in human subjects. Quantitative modeling is essential for distinguishing between potential physiolog-ical mechanisms in such studies, and a model that also covers multi-ple other phenomena is essential to enable uni ﬁcation and avoid ad hoc interpretations. It is worth pointing out that the data sets of Herrmann (2001) and Rager and Singer (1998) were acquired under quite different experi-mental conditions, yet our model reproduces the main features of both. Herrmann (2001) recorded scalp EEGs from awake human sub-jects, while Rager and Singer (1998) recorded intracortical multi-unit activity (MUA) and local ﬁeld potentials (LFPs) from anesthetized cats. While EEGs necessarily involve ﬁltering through the scalp, this is predominantly linear ( David et al., 2006; Nunez, 1995; Robinson et al., 2001a ) and so does not change the observed frequencies. In any case, EEGs and cortical LFPs are closely related because the same brain structures (mainly pyramidal cells) are ultimately respon-sible for both ( Nunez, 1995; Steriade et al., 1997 ), and indeed the re-sults of Herrmann (2001) and Rager and Singer (1998) are in good agreement. Two differences are that the cat LFP data do not show ev-idence of the subharmonics and gamma-band resonances seen in the human data; however, preliminary investigations indicate that our model can also reproduce this situation for plausible parameter values, though the exploration of appropriate parameters for the anesthetized cat is beyond the scope of this paper, and any such 50f (Hz) f D(Hz)45 40353025201510 50 10 20 30 40 50 Fig. 5. Spectrum vs. f Dfor parameters with a bistable seizure limit cycle attractor for the parameters used by Breakspear et al. (2006). Dark shading corresponds to high values. 1953 J. A. Roberts, P. A. Robinson / Neuro Image 62 (2012) 1947-1955