Optogenetic Stimulation Shifts the Excitability of Cerebral Cortex from Type I to Type II: Oscillation Onset and Wave Propagation

Constant optogenetic stimulation targeting both pyramidal cells and inhibitory interneurons has recently been shown to elicit propagating waves of gamma-band (40–80 Hz) oscillations in the local field potential of non-human primate motor cortex. The oscillations emerge with non-zero frequency and small amplitude—the hallmark of a type II excitable medium—yet they also propagate far beyond the stimulation site in the manner of a type I excitable medium. How can neural tissue exhibit both type I and type II excitability? We investigated the apparent contradiction by modeling the cortex as a Wilson-Cowan neural field in which optogenetic stimulation was represented by an external current source. In the absence of any external current, the model operated as a type I excitable medium that supported propagating waves of gamma oscillations similar to those observed in vivo. Applying an external current to the population of inhibitory neurons transformed the model into a type II excitable medium. The findings suggest that cortical tissue normally operates as a type I excitable medium but it is locally transformed into a type II medium by optogenetic stimulation which predominantly targets inhibitory neurons. The proposed mechanism accounts for the graded emergence of gamma oscillations at the stimulation site while retaining propagating waves of gamma oscillations in the non-stimulated tissue. It also predicts that gamma waves can be emitted on every second cycle of a 100 Hz oscillation. That prediction was subsequently confirmed by re-analysis of the neurophysiological data. The model thus offers a theoretical account of how optogenetic stimulation alters the excitability of cortical neural fields.


Introduction
Lu and colleagues [1] recently transduced small regions of primary motor (M1) and ventral premotor (PMv) cortices of macaque monkeys using red-shifted opsin C1V1(T/T). They found that constant optical stimulation of the targeted tissue induced intrinsic gamma-band (40-80 Hz) oscillations in the local field potential ( Fig 1A). The gamma oscillations were manifest in 4x4 mm 2 microelectrode recordings as patterns of concentric rings and spiral waves that propagated into the surrounding tissue well beyond the stimulation site ( Fig 1B). When the optogenetic stimulation was slowly ramped from zero, the oscillations emerged abruptly at  [1]. A: Gamma oscillations in the local field potentials at five recording sites on the microelectrode array for subject T. The oscillation phase has a spatial gradient that indicates wave propagation. Black dots indicate the peaks of one gamma cycle across neighboring electrodes; B: Optogenetic stimulation induces expanding waves, as summarized in the phase-triggered average of gamma (40-110 Hz) spatial field potential, based on the phase of the optogeneticallyinduced 50 Hz gamma oscillation. The dot indicates the point where the fiber optic light source was surgically inserted. The tip of the optical fiber was likely slanted to the right of this point, corresponding to the origin of the waves. C: Trial-averaged spectrogram of the local field potential when the optical stimulation was ramped from 0 mW to 6 mW over 4 seconds. The mean power within each frequency band for the 500 ms preceding stimulation was subtracted (in dB) from the power during stimulation to enhance visualization of the optogentically-induced changes. D: Power of the oscillations in the local field potential during the ramp protocol. a non-zero frequency ( Fig 1C) with low amplitude (Fig 1D). Lu et al recognized that these two characteristics are consistent with a dynamical system undergoing a supercritical Hopf bifurcation [1].
The supercritical Hopf bifurcation is the hallmark of type II neural excitability [2]. It encapsulates the dynamical properties of neurons that can fire arbitrarily small spikes but have a relatively fixed firing rate [3]. Type I neurons, on the other hand, are characterized by fixed amplitude spikes. Those dynamics can arise from a subcritical Hopf bifurcation or a saddlenode bifurcation on an invariant circle (SNIC) [2,4]. The SNIC bifurcation allows arbitrarily small firing rates whereas the subcritical Hopf bifurcation has relatively fixed firing rates. Izhikevich [4] characterizes type I neurons as integrators and type II neurons as resonators. In the present study, we apply the classifications of type I and II excitability in single neurons to the excitability of populations of neurons in spatially extended neural media. Those same classifications determine how well an excitable medium can sustain propagating waves of activity. For the case of a type I medium, small disturbances to the resting state produce large amplitude responses that readily propagate over long distances. Whereas small disturbances in type II media typically elicit only small responses that tend not to propagate. There is an exception however-type II media can become highly excitable when the time course of the recovery variable is substantially slower than that of excitation [5]. In such cases, the dynamics are that of a relaxation oscillator [6] which is capable of supporting propagating waves because of its explosive response to small inputs.
Lu et al's [1] observation that optogenetically-induced gamma waves propagate far beyond the site of stimulation suggests that the cortical tissue normally operates as a type I excitable medium. However this seems to contradict the finding that the optogenetically stimulated tissue operates as a type II excitable medium as revealed by the ramped stimulus protocol. This apparent contradiction offers a glimpse into the effect of optogenetic stimulation on the excitability of neural tissue. We explored the theoretical implications by simulating the propagation of gamma waves in a continuum neural field model of cortex. The model comprised of recurrently connected populations of excitatory and inhibitory neurons which were driven by an external current source that represented the ionic currents induced by optogenetic stimulation. We sought to determine (i) the conditions under which a type II excitable medium could sustain propagating waves of gamma oscillations and (ii) whether the act of optogenetic stimulation could transform a type I excitable medium into type II that produces graded oscillations.

Models
The cortex was modeled as a continuum neural field of recurrently connected populations of excitatory and inhibitory neurons following the methods of Wilson and Cowan [7,8]. The neural field represents the spatiotemporal intensity function for neural firing at spatial position x firing a spike at time t. The excitatory and inhibitory populations were treated as separate but interconnected neural fields. The equations governing their firing rates were defined as where U e (x, t) and U i (x, t) represent the mean firing rates of excitatory and inhibitory populations respectively. The local field potential (LFP) was defined as a weighted sum of the local mean firing rates, with the contribution of excitatory cells weighted four times that of inhibitory cells. This weighting reflects the higher prevalence of excitatory cells as well as their greater contribution to the electric field due to the arrangement of their dipoles. The firing rates were related to synaptic input by the sigmoidal function, where u(x, t) represents the local synaptic input at spatial position x at time t. The local synaptic input was computed as a weighted sum of excitatory and inhibitory spiking activity in the immediate vicinity plus external currents J e (x, t) and J i (x, t) that represented the additional synaptic currents induced by optogenetic stimulation. Spatial summation is denoted by the convolution operator, where K(x) is the spatial density profile of the lateral neural projections. That spatial profile was assumed to be Gaussian with distance, where σ is the spatial spread and k ei is the weight associated with the specific connection type indicated by the subscript (inhibitory-to-excitatory in this case). The connection weights k ee , k ei , k ie , k ii differed by connection type and connections emanating from excitatory populations were assumed to have twice the spatial spread of those emanating from inhibitory populations. See Table 1 for specific parameter values. Parameters b e and b i represent the firing thresholds of the excitatory and inhibitory neurons. Parameters τ e and τ i are the time constants of excitation and inhibition. The connections weights and the relative time course of excitation and inhibition both impact the excitability of the neural dynamics which, in turn, effects the capacity of the neural field to support propagating waves. We identified the parameters under which steady optogenetic stimulation replicated the emergence of gamma oscillations in the local field potential. We then explored how far those oscillations propagated away from the stimulation site. See [9] and [10] for reviews of the Wilson-Cowan model. See [11] for a review of continuum neural fields in general.

Results
We began by analyzing the excitability of an isolated pair of excitatory-inhibitory populations. The spatial coupling kernels K(x) in Eqs (1 and 2) were replaced with scalar connection weights (k ee = 15, k ei = 15, k ie = 15, k ii = 7). The connection weights and threshold parameters (b e = 4, b i = 4) were chosen so that the sigmoidal nullcline of the inhibitory population U i intersected the cubic-shaped nullcline of the excitatory population U e near the left knee of the cubic (Fig 2A). This particular configuration is known to undergo a supercritical Hopf bifurcation when an external current is applied to the excitatory population [12,13]. As anticipated, the model produced intrinsic oscillations in the simulated LFP when a steady external current (J e = 2) was applied to the excitatory cell ( Fig 2B). The time constants of excitation (τ e = 2 ms) and inhibition (τ i = 4 ms) were chosen so that the frequency of this oscillation always fell within the 40-80 Hz gamma band ( Fig 2C). The gamma oscillations emerged at zero amplitude and grew monotonically as the external current was increased ( Fig  2F). In all, the isolated pair of excitatory-inhibitory populations replicated the characteristics of gamma oscillations observed by [1] during ramped optogenetic stimulation. The next step was to investigate whether those gamma oscillations would propagate in a spatially extended medium.
The present type II model is only weakly excitable because small perturbations do not induce large responses in U e (Fig 2D and 2E). Nonetheless, excitability can be enhanced by increasing the time course of inhibition relative to excitation [5]. The same model with τ e = 1 ms instead of τ e = 2 ms readily evokes large responses to small perturbations (Fig 2G and 2H). The degree of excitability is evident in the bifurcation structure of U e versus J e . The amplitude of the oscillation rises gradually in the case of the weakly excitable system ( Fig 2F) and explosively in the case of the highly excitable system (Fig 2I). In the next section, we investigate how the degree of excitability of a type II spatial medium governs its ability to sustain propagating waves of gamma oscillations.

Wave propagation in a type II spatial medium
We modeled a 4 mm long strip of cortex as a type II neural medium using a chain of reciprocally coupled excitatory-inhibitory populations. It represented neural tissue spanning the width of the microelectrode array used by Lu et al [1]. The populations were evenly spaced at intervals of dx = 0.01 mm and coupled with a Gaussian spatial density profile Eqs (5 and 6). The Gaussian spread parameters were σ e = 0.2 mm for excitatory cells and σ i = 0.1 mm for inhibitory cells. The axons of the excitatory cells thus reached further than those of the inhibitory cells. Optogenetic stimulation was approximated by a focal current source with a square spatial profile (0.4 mm wide) that was centered on the midpoint of the chain (x = 0). We surveyed the distance that the waves propagated from the stimulation site for a range of excitatory time constants τ e while fixing τ i = 4 ms to preserve the frequency of the gamma oscillations as much as possible. The amplitude of the external current was also fixed (J e = 1.1). The propagation distance was designated as that point x where the maximal value of U e (x, t) fell below 0.05. Absorbing boundary conditions were imposed at both ends of the chain to prevent waves from reflecting back into the medium.
We found that gamma oscillations failed to propagate at all for τ i /τ e < 4. Waves that propagated a finite distance ( Fig 3A) were observed for values of τ i /τ e between 4 and approximately 12.5. The propagation distance grew rapidly as τ i /τ e approached 12.5 and the waves appeared to propagate indefinitely (Fig 3B and 3C) for τ i /τ e ≳ 12.5 ( Fig 3D). The bifurcation diagram (Fig 3E) reveals how the oscillation amplitude rises almost instantaneously for the case of τ i /τ e = 12.5. Such explosive growth is contrary to the slow rise in gamma power that is observed in the optogenetic ramp data (Fig 1D). We conclude that, while a type II neural medium could produce sustained traveling waves with a large difference in the timescales of excitation and inhibition, such differences were biologically unreasonable and the resulting medium could not support graded oscillations.

Wave propagation in a type I spatial medium
Type I excitablilty is associated with either a saddle-node bifurcation on an invariant circle (SNIC) or a subcritical Hopf bifurcation [2,4,12,13]. In our model, the supercritical Hopf regime (type II excitability) is readily transformed to a SNIC (type I excitability) by shifting the inhibitory (S-shaped) nullcline rightwards until it intersects the middle branch of the excitatory (cubic) nullcline ( Fig 4A). This was done by increasing the threshold of inhibition from b i = 4 to b i = 8. The SNIC yielded large responses to small perturbations from the rest state ( Fig  4B) and its high excitability was also evident in the instantaneous onset of large oscillations in the bifurcation diagram ( Fig 4C) even with modest time constants (τ i = 4 ms, τ e = 2 ms). The SNIC is therefore an ideal dynamical regime for sustaining propagating waves in spatial media. More importantly, the SNIC can be transformed back into the supercritical Hopf regime simply by injecting an external current into the inhibitory cell (J i = 4).
We propose that optogenetic stimulation likewise transforms the excitability of cortical tissue from type I to type II. Graded gamma oscillations would thus emerge at the stimulation site via a supercritical Hopf bifurcation while the non-stimulated tissue would continue to support propagating waves via the highly excitable dynamics of the SNIC regime. We tested these predictions using the same spatial model as before but in this case we varied J e while holding J i = 4 fixed. This scenario represents the gradual recruitment of excitatory cells by optogenetic stimulation coinciding with the instantaneous recruitment of inhibitory cells. As always, J e and J i were both set to zero in the un-stimulated spatial region (|x|>0.25). Low-amplitude gamma oscillations still emerged at the stimulation site, as expected, although they did not emit propagating waves when the injection current was low (Fig 5A). Waves were emitted at higher stimulation stimulation currents but not necessarily on every cycle of the gamma oscillation. For example, waves were emitted on every third gamma cycle for the case of J e = 2.2 ( Fig 5B) and every second gamma cycle for the case of J e = 3 ( Fig 5C). Importantly, the waves propagated indefinitely in the medium whenever they were emitted, as is expected of a type I excitable medium. These findings suggest that an excitable neural medium can operate in either the type I or type II regimes depending upon the influence of an external current source.
Intriguingly, one-to-one emission of waves on every gamma cycle was never observed for the model with τ i /τ e = 2. Although it could be observed with τ i /τ e !3.1, we considered this scenario unphysiological as it precluded the slow ramp in oscillation amplitude observed experimentally (Fig 1D), giving rise instead to a rapid rise in amplitude as shown in Fig 3E. In light of this observation, we returned to the experimental data to investigate whether the traveling

Comparison to the neurophysiological data
We re-analyzed the neurophysiological data from [1] to allow direct comparison with our simulations in one spatial dimension. Non-human primate recordings and optogenetic stimulation were implemented as described in [1] under the approval of the Institutional Animal Care and Use Committee (IACUC). We observed a second peak in the LFP power spectral density at *100 Hz confined to electrodes within the region of direct optogenetic stimulation, suggesting that the *50 Hz traveling waves may indeed originate from a 2:1 resonance with local, higher-frequency gamma oscillations. For visualization, broad gamma-band (40-110 Hz) signals in the multi-electrode array were averaged according to radial distance from the site of the optogenetic source. Fig 6A shows a typical example of gamma waves being emitted from neural tissue under steady optogenetic stimulation at 6 mW. Color indicates the amplitude of the local field potential after band limiting to 40-100 Hz. Traveling gamma waves emerged within 1 mm of the optogenetic source and propagated into the surrounding tissue at 18.9 cm/s on average (SD 4.76). Remarkably, the neurophysiological data itself exhibits wave emission on every second cycle of the gamma oscillation. It is most clearly visible in the phase-averaged data ( Fig 6B) where each wavefront in the *50 Hz oscillation is emitted on the second cycle of the *100 Hz oscillation at the source. In Fig 6B the LFP time-points were binned according to the phase of 48 Hz oscillations at the center of stimulation, and 45-100 Hz LFP was then averaged within each phase bin. Until now, we had only seen this behavior in simulations. Fig 6C  shows the simulated results where the space constants (σ e = 0.6 mm, σ i = 0.3 mm) and time constants (τ e = 3.6 ms, τ i = 1.8 ms) of the model have been re-scaled to match the wave speed of the neurophysiological data. We find that these simulations show a remarkable agreement to the phase-averaged neurophysiological data, indicating that the 2:1 resonance is a surprising but physiologically realistic prediction of the neural field model.

What is the relationship between J e and J i ?
It remains to ask how the inhibitory current (J i ) might realistically vary with excitatory current (J e ) during optogenetic stimulation. We reasoned that the optogenetically-induced currents must originate from J e = 0 and J i = 0 and increase smoothly with stimulation. Furthermore, the inhibitory current must saturate at J i = 4 for the type I regime to be transformed to a type II regime. To this end we assumed a sigmoidal relationship, where β is an unknown slope parameter. The curve of (J e , J i ) points (Eq 7) must pass through a Hopf bifurcation where the dynamics shift from fixed points to oscillations. We used numerical continuation to map the critical values of J e and J i where those Hopf bifurcation points occur (Fig 7A). We chose the slope parameter β = 3 so that the curve of (J e , J i ) points intersected the line of Hopf bifurcation points near J e = 1. More pertinently, that choice allows the (J e , J i ) curve to closely graze the Hopf bifurcation points in the vicinity of the intersection point. Doing so facilitates the gradual rise in the oscillation amplitude as the critical point is traversed. That slow growth in the oscillation is evident in the bifurcation diagram ( Fig 7B) where J e is varied while J i is governed by Eq (7). It is also observed in the simulated ramp protocol (Fig 7C) where the excitatory current was slowly increased from J e = 0 to J e = 3 over a period of t = 4 seconds to mimic the ramp protocol by [1]. We argue that this slow increase in the amplitude of the gamma oscillation is what accounts for the slow rise in the power of the gamma-band oscillations reported by [1]. The slight delay in the onset of the oscillations in the simulated ramp (Fig 7C) is due to critical slowing in the vicinity of the Hopf bifurcation. We expect that the neural tissue would likewise exhibit critical slowing. This prediction could be tested by comparing the time course of optogeneticallyinduced gamma oscillations in the ramp-down protocol versus the existing ramp-up protocol.

Discussion
Lu et al [1] correctly recognized that the onset of optogenetically-induced gamma oscillations in their experimental setup and protocols is consistent with a supercritical Hopf bifurcation. The supercritical Hopf is the hallmark of a type II excitable system. Yet they also observed wave propagation which is typically associated with type I excitability. We used a Wilson-Cowan [7,8] neural field model to identify those conditions under which a type II excitable medium might support propagating waves like those observed in vivo. We found that the model could only do so if the time course of inhibition was substantially slower than that of excitation (τ i /τ e ≳ 12). Under such conditions, the oscillation amplitude grows explosively as the stimulation is increased. However such explosive growth is inconsistent with the slow rise in the gamma power observed by [1] in their stimulus ramp protocol. We conclude that a type II excitable medium with the dynamical properties reported by [1] is not capable of supporting propagating waves.
We instead propose that the neural tissue typically operates as a type I excitable medium and that it can be locally transformed into type II by the very act of optogenetic stimulation activating the inhibitory interneurons. Type I neural medium are exemplified by arbitrarily small firing frequencies and large responses to small perturbations-characteristics which are ideal for sustaining indefinite propagating waves. Our analysis shows that a Wilson-Cowan neural field with type I excitability due to a SNIC bifurcation can be readily transformed into a type II excitable system by strongly stimulating the inhibitory cells. Simulations confirm that under these conditions the model reproduces the graded increase in gamma oscillations at the stimulation site while simultaneously supporting wave propagation in the regions beyond the stimulation site. The model thus accounts for the two major observations of [1] and resolves the apparent contradiction of type I and type II excitability.
The question remains as to how optogenetic stimulation might differentially activate inhibitory and excitatory neurons as our model suggests. Analysis of the phase plane shows that the inhibitory neurons must be recruited early and strongly in order to shift the nullcline of U i from the left hand branch of the nullcline of U e to the right hand branch. For simplicity, we have presented this in our model as a sigmoidal relationship (Eq 7) between the external currents J e and J i that represent the ionic currents induced by optogenetic stimulation. It is known that the optogenetic construct (CaMKII alpha promoter and AAv5 viral vector) used by Lu et al [1] expresses primarily in layer 5 pyramidal (excitatory) neurons and to a lesser extent in parvalbumin-positive inhibitory interneurons [14]. However it is not clear from that work whether the optogenetic construct activates inhibitory neurons earlier than excitatory neurons as our model suggests. Further studies are required to investigate the rate at which excitatory and inhibitory neurons are recruited by optogenetic stimulation.
It is reasonable to also ask whether the neural medium can be transformed into type I excitability by shifting the dynamics from the supercritical to the subcritical Hopf regime rather than to the SNIC as we suggest. Neural fields with subcritical Hopf dynamics have previously been shown to emit solitary and N-pulse traveling waves when subjected to constant stimulation [15][16][17]. The bistability associated with the subcritical Hopf bifurcation allows those models to support co-existing resting state and wave pulse solutions. The resting state loses stability when a constant injection current is applied to the field, leaving only the stable wave solution. When that stimulation is applied focally, it induces a local oscillation that emits wave pulses which propagate throughout the resting medium [16]. As with the present model, the wave pulses are emitted with a range of n:m mode locking regimes when the stimulation is weak [16]. In our model and theirs, the mode locking is due to the time course of the recovery variable which can block wave propagation if it remains high from a previous oscillation cycle. However the bistable models lack the supercitical Hopf dynamics needed to replicate the graded rise in the gamma oscillation observed by Lu and colleagues [1]. Moreover, it is not clear how the dynamics in those models could be transformed from subcritical to supercritical Hopf in a biologically plausible fashion. Especially since the existence of the subcritical Hopf regime seems to depend on the high gain limit of the Heaviside firing rate function, which is presumably a fixed property of the biology. Thus we regard the present model as a more compelling account of the neurophysiological data.
Interestingly, wave emission in the model of Folias and Bressloff [16] is only observed when the stimulation is ramped from high to low. In that scenario, the initial strong stimulation produces a stable stationary wave solution which is transformed into a pulsating 'breather' when the stimulation falls below a critical level. The pulsations become more pronounced with subsequent reductions in stimulation until the breather eventually emits traveling pulses. Yet the breather fails to materialize when the stimulation is instead ramped from low to high because of hysteresis in the bistable dynamics. Importantly, that hysteresis presents an opportunity to test the competing models against the neurophysiological data. Since our model is monostable, it predicts that waves will be emitted by optogenetic stimulation irrespective of the direction of the stimulus ramping protocol. Whereas the bistable model predicts that waves will only be emitted when the optogenetic stimulus is ramped downwards. The differences ought to distinguishable in the neurophysiological data-notwithstanding the effects of critical slowing mentioned earlier.
Unfortunately, Lu and colleagues [1] did not conduct their ramp protocol in both directions, so a new experiment must be conducted to test the hypothesis. Nonetheless, our re-analysis of their data has already confirmed some predictions of the present model. Notably, the 2:1 mode-locking of the wave emission correctly predicted the existence of *100 Hz oscillations at the stimulation site with the *50 Hz gamma oscillations being restricted to the peripheral tissue. Although we acknowledge that the same prediction also follows from the model of Folias and Bressloff [16]. The same behaviors have yet to be fully investigated in spiking neuron models. Conductance-based models of recurrently connected inhibitory neurons with type I excitability have been shown to elicit gamma-band oscillations in the macroscopic network behavior [18]. That oscillation is due to the synchronization of individual neurons which do not necessarily fire on every cycle. Moreover, it emerges from the asynchronous state through a supercritical Hopf bifurcation. Hence the dynamics of the population oscillation in the spiking neuron model are consistent with that of our neural field model. Whether those oscillations can propagate in a spatially extended spiking neural network has yet to be investigated.
In summary, the present model suggests that optogenetic stimulation can locally transform the excitability of cortical tissue from type I to type II by recruiting inhibitory interneurons prior to recruiting excitatory neurons. In doing so, it accounts for the seemingly contradictory observations of traveling waves and and supercritical Hopf bifurcation dynamics in the neurophysiological data. Further neurophysiological studies are required to determine whether optogenetic stimulation does indeed differentially recruit inhibitory and excitatory neurons as we propose. In addition, our model makes the testable prediction that gamma waves are induced at the same critical level of optical stimulation, irrespective of whether the stimulation is ramped upwards or downwards. This prediction allows our model to be empirically distinguished from the bistable model of Folias and Bressloff [16] which predicts that wave-emitting breathers only arise when the optogenetic stimulation is ramped downwards.

Author Contributions
Conceptualization: SH BE.