Corticothalamic model Neural ﬁeld theories have successfully reproduced many aspects of the dynamics of populations of neurons by averaging over small spatial and temporal scales to give a description at scales ≳0. 1 mm; see Deco et al. (2008) for a recent review. We use a recent corticothalamic model (Breakspear et al., 2006; Robinson et al., 1997, 2002, 2004 )t oa n a l y z e nonlinear SSVEPs ( Roberts, 2010; Robinson et al., 2002, 2008 ). We fol-low the exposition and notation of Robinson et al. (2002) here, except where otherwise noted. The neuronal populations included are (see Fig. 2 ) excitatory ( e) and inhibitory ( i) cortical neurons, and the speci ﬁc relay ( s) and reticular ( r) thalamic neurons. The relay in the visual sys-tem is the lateral geniculate nucleus, which has reciprocal connections with both the cortex and the reticular nucleus as shown in Fig. 2 (Sherman and Guillery, 1996; Steriade et al., 1997 ). External input from the brainstem and sensory systems ( n) drives the system via s. U n-structured stimuli presented to the entire visual ﬁeld drive a large area of cortex, so we concentrate on spatially uniform activity. A nonlinear sigmoidal function S(Va), given by ( Freeman, 1975; Wilson and Cowan, 1973 ) SVað Þ ¼Qmax 1þexp-Va-θ ð Þ =σ′/C2/C3 ; ð1Þ relates the mean ﬁring rate Qa=S(Va) of each neuronal population a to its mean cell body potential Va(relative to resting), where Qmaxis the maximum ﬁring rate, θis the mean threshold voltage and σ′π=ﬃﬃﬃ 3p is its population standard deviation. Here Va=∑b Vab, with Dab Vabr;tð Þ ¼ νabϕbr;t-τab ð Þ ; ð2Þ Dab¼1 αabβabd2 dt2þ1 αabþ1 βab/C18/C19d dtþ1; ð3Þ where Vabis the contribution due to input ϕbfrom neurons in popu-lation barriving after mean axonal delay τab,Dabgives the soma po-tential response allowing for synaptodendritic dynamics and soma capacitance, connection strength νab=Nabsab,Nabis the mean num-ber of synapses per neuron of type afrom neurons of type b,sabis the strength of response to a unit input from neurons of type b, and αabandβabare the inverse decay and rise times of the soma response to input b, respectively. For the system in Fig. 2, the only nonzero de-laysτabareτes=τis=τse=τre=t0/2 where t0is the total cor-ticothalamic loop delay. Average afferent ﬁring rates ϕapropagate as ﬁelds in the continu-um limit, approximately governed by a damped wave equation with source Qa(Jirsa and Haken, 1996; Nunez, 1995; Robinson et al., 1997 ), giving 1 γ2 a∂2 ∂t2þ2 γa∂ ∂tþ1-r2 a∇2"# ϕar;tð Þ ¼ SVar;tð Þ1/2/C138 ; ð4Þwhere γa=va/ra,vais the axonal conduction speed, and rais the char-acteristic axon length. We assume that the activity is spatially uni-form over the cortical area of interest, consistent with a large-scale mode driven by ﬂicker stimuli; thus the ∇2term is zero and for nota-tional simplicity we omit the rargument hereafter. We set γa=∞for a=i,r,sowing to the short ranges of these axons ( Robinson et al., 2002 ). We retain long-range excitatory intracortical connectivity via nonzero reandﬁniteγe. Note that, although we neglect position de-pendence of activity, we have not reduced our model to a neural mass model, since this would also require shrinking the cortical man-ifold to be zero-dimensional by setting re=0 and γe=∞. We assume random intracortical connectivity, which implies νib=νeb forb=e,i,sand hence Vi=Ve, as derived in detail previously ( Robinson et al., 1998 ). The inhibitory population is modeled explicitly ( νei,νii,νie, andνisare all nonzero); random connectivity implies that its mean-ﬁeld dynamics are closely related to (but not identical to) those of the excit-atory population with which it is spatially commingled (the excitatory neurons project much further, giving rise to the parameter γe). For sim-plicity we assume that all the synaptodendritic time constants are equal except for the inhibitory intrathalamic srconnection, which is mediated by GABA Bdynamics slower than the other channels ( Steriade et al., 1997 ). Hence we set αab=α,βab=β,Dab=Dαfor all abexcept for the srconnection; αandβare thus to be interpreted as effective values. We do not explicitly model the bursting mode of thalamic neurons, which tends to only be relevant to frequencies lower than those typical-ly observed in SSVEPs ( Steriade et al., 1997 ). We set parameters at or near previously published values (see ﬁgure captions) ( Breakspear et al., 2006 ), which enables uni ﬁcation of the ﬁndings here with those of previous studies. Here the parameters place the system near a stability boundary associated with the intrathalamic loop ( Robinson et al., 2002 ). All model parameters corre-spond to physiological quantities and have been rigorously constrained across numerous independent experimental measures ( Robinson et al., 2004 ). Thus for the connectivity in Fig. 2 and the above assumptions, the model equations are 1 γ2 ed2 dt2þ2 γed dtþ1"# ϕetð Þ¼SVetð Þ1/2/C138 ; ð5Þ DαVetð Þ¼ νeeϕetð Þ þνeiϕitð Þ þνesϕst-t0=2 ð Þ ; ð6Þ DαVrtð Þ¼ νreϕet-t0=2 ð Þ þ νrsϕs; ð7Þ Vstð Þ¼ Vsetð Þ þ Vsrtð Þ þ Vsntð Þ; ð8Þ DαVsetð Þ¼νseϕet-t0=2 ð Þ ; ð9Þ Dsr Vsrtð Þ¼ νsrϕrtð Þ; ð10Þ DαVsntð Þ¼ νsnϕntð Þ; ð11Þ with ϕb(t)=S[Vb(t)] for b=i,r,s. External stimuli ϕndrive the brain via the relay nuclei. Resting EEG spectra are successfully generated when driving the model with white noise ( Robinson et al., 1997, 2001a, 2001b, 2002, 2004 ), which assumes that the collective input from all sensory pathways is so complex that it does not favor any particular frequency in the range of interest ( Engel et al., 2001; Lopes da Silva et al., 1974; Stam et al., 1999 ). Impulsive stimuli yield evoked response potentials (Kerr et al., 2008; Robinson et al., 2001a, 2004 ). Sinusoidally modulat-edﬂicker stimuli generate SSVEPs ( Robinson et al., 2002, 2008 ) and in the linear regime the model yields cortical phase velocities for SSVEPs consistent with experiment ( Robinson et al., 2008 ). Here we concen-trate on the nonlinear regime, where the drive is strong.,rs sn se reesrelay nuclei reticular nucleus cortex ee,eisr Fig. 2. Connectivities of cortical (excitatory e, inhibitory i) and thalamic (relay s, retic-ular r) neural populations; arrowheads abdenote where ﬁeldsϕbproject, with connec-tion strength νabfor input to afrom b;ϕndenotes sensory input. 1949 J. A. Roberts, P. A. Robinson / Neuro Image 62 (2012) 1947-1955