Abstract
Gamma rhythm (30–70 Hz), thought to represent the interactions between excitatory and inhibitory populations, can be induced by presenting achromatic gratings in the primary visual cortex (V1) and is sensitive to stimulus properties such as size and contrast. In addition, gamma occurs in short bursts and shows a “frequency falloff” effect where its peak frequency is high after stimulus onset and slowly decreases to a steady state. Recently, these size-contrast properties and temporal characteristics were replicated in a self-oscillating Wilson–Cowan (WC) model operating as an inhibition stabilized network (ISN), stimulated by Ornstein–Uhlenbeck (OU) type inputs. In particular, frequency falloff was explained by delayed and slowly accumulated inputs arriving at local inhibitory populations. We hypothesized that if the stimulus is preceded by another higher contrast stimulus, frequency falloff could be abolished or reversed, since the excessive inhibition will now take more time to dissipate. We presented gratings at different contrasts consecutively to two female monkeys while recording gamma using microelectrode arrays in V1 and confirmed this prediction. Further, this model also replicated a characteristic pattern of gamma frequency modulation to counter-phasing stimuli as reported previously. These phenomena were also replicated by an ISN model subject to slow adaptation in feedforward excitatory input. Thus, ISN model with delayed surround input or adapted feedforward input replicates gamma frequency responses to time-varying contrasts.
Significance Statement
Gamma rhythms represent interactions between excitatory and inhibitory populations during visual stimulation. Gamma power and center frequency varies depending on stimulus features, and onset of stimulus produces a “frequency fall” trend where onset frequency is higher and subsequently plateaus to a lower value. In an earlier work, we argued, using a noisy rate model of V1, that a delayed onset of inhibition-drive from the surround populations produced the gamma “frequency falloff.” We tested a key prediction of this hypothesis that the frequency falloff can be abolished or reversed if the stimulus is preceded by a higher contrast stimulus and confirmed the same by recording from primate primary visual cortex while presenting multiple stimuli consecutively at varying contrasts.
Introduction
Gamma rhythm (30–70 Hz) can be induced in the primary visual cortex (V1) by presentation of suitable achromatic gratings or plain hues and represents interactions between excitatory and recurrently connected inhibitory populations. Such gamma generation has been demonstrated in network models of V1 (Buzsáki and Wang, 2012; Chariker et al., 2018; Zachariou et al., 2021), suggesting gamma rhythm characteristics are tightly linked to the recurrent circuit within V1. Gamma properties tend to reflect possible changes in cortical circuitry due to aging (Murty et al., 2020) and mild cognitive impairment (MCI; Murty et al., 2021).
Gamma rhythm exhibits differences in power and peak frequency depending on the properties of stimulus presented. For instance, gamma frequency is higher during presentation of higher contrast stimuli (Ray and Maunsell, 2010), while gamma power increases and frequency decreases when larger-size stimuli are presented (Gieselmann and Thiele, 2008; Ray and Maunsell, 2011; Peter et al., 2019). These properties of V1 gamma had been replicated in a Wilson–Cowan (WC)-type model, which we refer to as the Jadi–Sejnowski (JS) model, that generates gamma as limit cycles produced by interaction between an excitatory and an inhibitory population (Jadi and Sejnowski, 2014; Shirhatti et al., 2022). Shirhatti and colleagues have previously shown that the existence of bifurcation in the JS model emulates the attenuation of gamma rhythm induced by a full-screen achromatic grating when a small discontinuity is introduced near the receptive field (RF) center (Shirhatti et al., 2022). Furthermore, gamma rhythms exhibit specific temporal characteristics such as nonsinusoidal waveform shape (Krishnakumaran et al., 2022) and bursty nature (Burns et al., 2011; Xing et al., 2012; Chandran et al., 2017), which have also been replicated in the JS model (Krishnakumaran and Ray, 2023). Gamma bursts with realistic burst length distributions were replicated in the JS model by supplying noisy inputs of Ornstein–Uhlenbeck (OU) type, obtained by low-pass filtering of Poisson distributed inputs. Furthermore, gamma rhythms exhibit a “frequency falloff,” where the gamma frequency is high near stimulus onset, and slowly reduces to a lower steady-state value, which is visible over the initial 300–400 ms in time–frequency (TF) spectra [for examples, see TF spectra in Xing et al. (2012), their Fig. 1; Peter et al. (2019), their Fig. 5, for LFP in macaques; Murty et al. (2020), their Fig. 4A, for EEG in humans; and Perry et al. (2020), their Fig. 3, for MEG in humans; the falloff is more visible in macaque LFP owing to better signal-to-noise ratio of LFP]. Replicating the frequency falloff in the JS model required the cutoff of filter be lower for the input drive to the inhibitory population than for the excitatory population (Krishnakumaran and Ray, 2023). This suggested a slower accumulation of inhibition, analogous to the temporal integration of inhibitory drive coming from surround populations over unmyelinated lateral connections arriving at different delays (Angelucci et al., 2017), and the differential synaptic dynamics of excitation and inhibition.
As inhibition is slower in accumulation and decay in our model, we predicted that a prior buildup of inhibition, caused by presenting a higher contrast before a lower contrast, would either reverse or suppress the trend of the frequency falloff for the lower contrast stimulus. In the current study, we explore this effect of contrast adaptation on the frequency falloff of gamma rhythm by presenting full-screen gratings of different contrasts consecutively with no interstimulus interval. Using time–frequency spectral analysis, we traced the frequency of gamma bursts at different times from stimulus onset and compared the frequency transient for a stimulus, under unadapted and adapted conditions. Further, we tested whether the noisy JS model replicates a characteristic pattern of gamma frequency modulation in response to slow time-varying contrast stimuli as reported previously (Ray and Maunsell, 2010).
Although the delayed inhibition model described above makes a specific prediction about the frequency falloff effect, observing the same in real data does not rule out other mechanisms which could also explain the same effect. In particular, the experimental paradigm used in this study has been used extensively to study the effect of adaptation, which could arise from multiple mechanisms including one where the feedforward input gets adapted [see Kohn (2007) for details; mechanisms are discussed in more detail later]. We therefore also tested a noisy JS model with a slow subtractive adaptation applied to the excitatory input drive to emulate feedforward input adaptation and tested its responses to the time-varying contrast stimuli.
Materials and Methods
Ethics statement
All the experiments on monkeys were conducted in adherence to guidelines approved by the Institutional Animal Ethics Committee of the Indian Institute of Science-Bangalore and the Committee for the Purpose of Control and Supervision of Experiments on Animals (CPCSEA).
Animal preparation and training
Two adult female monkeys (Macaca radiata; M1: 13 years, ∼3.3 kg and M2: 15 years, ∼5.6 kg) were used. Each monkey was surgically implanted with a titanium headpost over the anterior/frontal region of the skull under general anesthesia. Once recovered, the monkey was trained on a visual passive fixation task. After the monkey learnt to maintain fixation for an adequate duration, a second surgery was performed under general anesthesia to insert a microelectrode array (Utah array, 81 active platinum microelectrodes in M1 and 48 electrodes in M2, 1 mm long each, 400 µm interelectrode distance; Blackrock Microsystems) in the primary visual cortex area V1 (right hemisphere, centered ∼10–12 mm rostral from the occipital ridge and ∼10–12 mm lateral from the midline, with location varying slightly in the two monkeys). In M2, a second microelectrode array was placed in V4 in the same surgery, whose signals were not utilized in the current study. The RFs of the recorded V1 neurons were located in the lower left quadrant of the visual space with respect to fixation (at an eccentricity of ∼3–4.5° in M1 and ∼1.5–3° in M2). Following a period of postsurgery care and monitored recovery of at least 2 weeks, the monkeys performed the experimental task regularly while microelectrode data were recorded.
Data acquisition
Raw signals from microelectrodes were recorded using the 128-channel Cerebus neural signal processor (Blackrock Microsystems). The signals were filtered online between 0.3 and 500 Hz (Butterworth filters; first-order analog and fourth-order digital, respectively) to get the LFP data, which were then recorded at a sampling rate of 2 kHz and a resolution of 16 bits. This LFP timeseries was used directly in our analyses.
Experimental setup and behavior
During the experiment, the monkey's head was held still by the headpost as it sat in a monkey chair and viewed a monitor (BenQ XL2411, LCD, 1,280 × 720 resolution, 100 Hz refresh rate) placed ∼50 cm from its eyes. A Faraday enclosure housed the subject and the display setup with a dedicated ground for isolation from external electrical noise. The monitor was calibrated and gamma-corrected using i1Display Pro (X-Rite PANTONE). Mean luminance was set to 60 cd/m2 at the monitor surface and gamma was set to unity for each of the three primaries.
The monkeys were subjected to a passive fixation task, in which they were rewarded with juice for fixating at a small dot of 0.05–0.10° radius at the center of the screen throughout the trial (4 s duration; fixation spot was displayed throughout). Each trial began with fixation period of 1 s duration, followed by consecutive presentations of an achromatic grating at three (sometimes two for M1) different contrasts (see below, Stimuli, on details on the grating), 1 s each with an interstimulus interval of 0 s, as depicted in Figure 1. The order of different contrasts in each trial was randomly chosen. Juice was awarded only if fixation was maintained within 2° from the fixation point. Trials in which fixation was broken were not considered in our analyses. Eye position data was recorded using the ETL-200 Primate Eye Tracking System (ISCAN) which gave horizontal and vertical coordinates. During the task, monitoring of eye position data, control of task flow, generation, and randomization of stimuli were performed using a custom software on macOS. Although fixation was required to be within 2° from the fixation spot, the monkeys were able to hold their fixation well within these limits during the task with a standard deviation of less than 0.4° across all sessions for both monkeys. These small eye movements are unlikely to affect our results since full-screen stimuli were used.
Summary of experimental paradigm and data analysis. Illustration of stimulus presentation during an example trial of passive fixation task. A full-screen achromatic grating is presented in three contrasts, and a blank gray screen is presented during the fixation period and as 0% contrast stimulus. A single trial of the passive fixation task consists of a fixation period followed by three consecutive contrast presentations without interstimulus intervals. Maintaining fixation at the center of the screen throughout the fixation period and the stimulus presentations results in a reward, whereas breaking it causes the trial to end without reward.
Stimuli
The set of stimuli comprised full-screen static achromatic grating of different contrasts (Fig. 1). The achromatic grating was at an orientation of 90° for both monkeys and had a spatial frequency of 4 cpd for M1 and 2 cpd for M2. These grating parameters were chosen to induce strong fast gamma and reduced slow gamma (Murty et al., 2018). These full-screen stimuli subtended a visual angle of ∼56° in the horizontal direction and ∼33° in the vertical direction. The gratings were presented at three different contrasts in M1 (25, 50 and 100%) and four different contrasts in M2 (0, 25, 50 and 100%; 0% contrast produces the same blank gray screen as during fixation period). An additional 0% contrast was included in M2 to consider unadapted conditions occurring later in a trial whereas in M1, unadapted conditions for a stimulus would refer to trials where a given contrast was the first to be displayed after fixation period. The results did not differ between unadapted conditions positioned at fixation period and later in trial in M2, and so the trials are pooled together for analyses. Adapted condition refers to the presentation of a stimulus, immediately preceded by a stimulus with nonzero contrast.
Electrode selection
As with our previous reports (Krishnakumaran et al., 2022; Krishnakumaran and Ray, 2023), we considered only those electrodes for analysis that had robust estimates of the RFs which were estimated by presenting small sinusoidal gratings at different eccentricities from the central fixation point. These criteria were determined by a RF mapping protocol that was run across multiple days (Dubey and Ray, 2019; Krishnakumaran and Ray, 2023). Further, high artifact electrodes were identified using the impedances measured on the days of recording and discarded. This procedure yielded 77 and 32 electrodes from M1 and M2, respectively.
Data analysis
We discarded stimulus presentations with excessive artifacts for each session (<5.9% of presentations of each stimulus in M1 and M2). The resultant number of useable trials for each contrast pair is given in Figure 3. For the 50% contrast condition in Figure 2, the number of useable trials were 93 trials in M1 and 87 trials in M2 for each unadapted presentations of 50% contrast (i.e., 0–50%) and 86 trials in M1 and 35 trials in M2 for the 100% contrast adapted (i.e., 100–50%). The trials in the unadapted case in M2 were more in number than the adapted case since this set included all trials where the nonzero contrast stimulus was presented immediately after the fixation period and other trials where it was presented later in the stimulus sequence but preceded by a 0% contrast stimulus. Furthermore, the software used to generate our trials as pseudorandomly generated sequences of contrasts used a blocked design where different contrasts (25, 50 and 100% for M1 and 0%, 25, 50 and 100% for M2) were chosen pseudorandomly presented without replacement and the next block started only once all contrasts were shown. Therefore, two consecutive stimuli with the same contrast could occur only when one block finished with a particular contrast and the next block started in the same trial with the same contrast, which happened rarely. This gave rise to the very low number of trials with the same contrast presented consecutively (or equivalently a single contrast presented for 2 s), as shown in the diagonal-occupying plots of Figure 3.
Effect of contrast adaptation on gamma in time–frequency spectrum. A, Change in TF spectrum of LFP from M1 from unadapted trials where 50% contrast was preceded by blank gray screen. The 0 s mark corresponds to the onset of 50% contrast. B, Same as A but from trials with 50% contrast following a 100% contrast. C, Peak frequency of gamma band for each condition. Solid line traces the median peak frequency of gamma across electrodes. The width of the line at each timepoint represents the standard error (SE), computed using bootstrapping. Black and blue traces indicate gamma-band peak-frequency traces corresponding to change in TF spectra in A and B, respectively. Horizontal blue bars at the top edge mark the timepoints where the difference between adapted and unadapted peak frequencies is statistically significant (p < 0.01). D–F, LFP spectra and gamma frequency analysis from M2.
Change in time–frequency spectra for different consecutively presented contrast pairs. A, Change in TF spectra from M1. Each column corresponds to a different “adaptor” contrast. Each row corresponds to a different “adapted” contrast that follows the adaptor. The sequence of contrasts is indicated in text above the respective plots. The plots corresponding to trials where the same contrast was presented for >2 s continuously occupy a diagonal of the figure. Note that these same contrast trials are much lower compared with other conditions, due to which the TF spectra are very patchy. For each contrast pair, the number of good trials used for TF analysis are indicated in text within the respective plot. Colorbars for each row are given at the right of the respective row. Note that colorbars have different scales across rows for better visualization of gamma frequency variation. B, Change in TF spectra from M2.
Time–frequency analysis of LFP and gamma peak-frequency estimation
Time–frequency (TF) spectra were computed using multitape method (MT) using the Chronux toolbox (Bokil et al., 2010), with moving windows of 0.4 s and three tapers. This gave us a frequency resolution of 2.5 Hz, which was needed to trace the gamma peak-frequency falloff. To compute the electrode averaged change in TF spectral power from baseline, for each stimulus, the TF spectrum from each trial corresponding to the presentation of the stimulus from each electrode were averaged. Then the baseline power at each frequency was computed by averaging the power at the respective frequency over the baseline period interval between −750 and 0 ms from stimulus onset. This baseline for each frequency was subtracted (on a log scale) from the averaged TF spectrum at every timepoint at the frequency. The peak-frequency timeseries of gamma was computed for each stimulus as the frequency within 20–80 Hz range with the maximum power in the electrode averaged change in TF spectrum at each timepoint. We used peak frequencies instead of mean power within a band because due to the “1/f” nature of the power spectral density in real data, raw power within a band is dominated by the lower frequency range where absolute power is higher.
Experimental design and statistical analysis
To estimate gamma peak-frequency trace at each timepoint, the peak frequencies estimated from all the electrodes at the respective timepoint were pooled and the median of the resulting distribution was computed and plotted in Figure 2. The standard error (SE) was estimated by bootstrapping over N iterations (where N is the number of electrodes). This involved random sampling with replacement of the data N times and estimating their median in each iteration, which resulted in N medians, whose standard deviation was taken as the SE.
To test any difference in temporal trend in gamma frequency between adapted and unadapted conditions, the peak frequencies were subject to a paired tailed sign test. First, the difference in peak frequencies recorded in a single electrode between adapted and unadapted conditions was taken at each timepoint. This resulted in a difference in frequency timeseries for each electrode. This trace was binned into nonoverlapping consecutive time windows of 10 ms starting from 0 ms. In each time window, the difference in frequency values between adapted and unadapted condition across all electrodes were pooled and subject to a right-tailed sign test. A p value <0.01 obtained for a time window rejects the null hypothesis that the median difference in frequency between adapted and unadapted conditions is 0 within the time window, suggesting an alteration in gamma frequency trace in the adapted condition, and such time ranges are marked in Figure 2C–F by horizontal bars at the top.
Jadi–Sejnowski model
Jadi and Sejnowski (2014) used a simple rate model consisting of an excitatory and an inhibitory population with sigmoidal activation, operating as an ISN, and constrained the input drives to the populations to reproduce the increase in power and decrease in peak frequency of gamma with increasing stimulus size as earlier studies have observed in V1 (Gieselmann and Thiele, 2008; Ray and Maunsell, 2011).
The model defines the population firing rates rE and rI of excitatory and inhibitory populations, respectively, as follows:
Increasing the contrast of a large stimulus increases both iE and iI and hence the system moves along the values indicated by magenta symbols in Figure 6A. For low contrast (Fig. 6B, left), the system only has stable or unstable fixed points and hence is out of the oscillatory regime. As contrast is further increased (Fig. 6B, middle), the system moves to the boundary of the superlinear regime (Fig. 6A, white contour), and the model exhibits a supercritical Hopf bifurcation. Across this bifurcation point, the model's local stability about the fixed point changes from stable to unstable (Fig. 6B, right). Hence when the system is at values (rE, rI) farther from the fixed point, the trajectory is attracted toward the fixed point, whereas at the neighborhood of the fixed point, the trajectory is repelled away from the fixed point, resulting in a self-sustained oscillatory trajectory [limit cycle; Jadi and Sejnowski (2014), their Fig. 3]. When the model operated close to a supercritical Hopf bifurcation, the amplitude and frequency of oscillations in firing rates could be closely approximated by linearization of the model [see Jadi and Sejnowski (2014) for details].
The parameter combination used in this study is given in Table 1. The connectivity weights were retained as that of the original model by Jadi and Sejnowski, (2014) but cE and cI [θE and θI in Jadi and Sejnowski (2014)] were changed to ensure that our baseline period input combination (iE, iI) = (0, 0) does not settle to multiple steady states but yields a fixed value of firing rates close to 0. Further, the time constants τE and τI were scaled by a factor of 0.77 or 0.66 for different input models (see below) from the original parameters to ensure that the model generates a wide range of limit cycle frequencies spanning the gamma range (listed in Tables 2, 3).
Parameter values used in JS model
Parameter values for noisy inputs with delayed inhibition
Parameter values for noisy inputs with adaptive excitation
Noisy input to JS model with delayed inhibition
In our study, we subject the model to noisy input drives iE and iI, which were both Ornstein–Uhlenbeck (OU) noise (as in our previous study in Krishnakumaran and Ray, 2023). The OU noise formulation of iE and iI was a first-order low-pass filtered white Gaussian noise:
The above formulation could be considered as an additive white Gaussian (AWG) noise passed through a first-order low-pass filter with a cutoff frequency at θP Hz for a population P. Further, at steady state, the variance would be
Noisy input to JS model with adaptive excitation
Additionally, we subject the JS model with OU noisy inputs with subtractive adaptation applied specifically to excitation. The unadaptive excitation and inhibition input drives are both OU type with cutoff at 16 Hz. But an adaptation signal, derived as a low-pass filtered version of the unadapted excitatory input, is subtracted from the unadapted excitatory input to obtain the net excitatory input applied to the model. The resultant inputs iE and iI to the JS model in Equation 1 is formulated as follows:
Simulation of time-varying contrast experiments
First, we simulated the step-contrast adaptation experiment, conducted in the current study. The contrast variable c(t) could range from 0 to 1 indicating the minimum and maximum inputs. For the contrast adaptation simulation, c(t) was varied between c0%, c50%, and c100%, indicating 0, 50, and 100% contrasts, respectively. The unadapted case was emulated by switching c(t) from c0% to c50%, whereas the adapted condition was from c100% to c50%. These values are specified in Tables 2 and 3 for the different input models. The experiment was simulated by presenting three contrasts values in sequence, namely, c(t) = c50%, c100%, and c50%, during intervals 0–1, 1–2, and 2–3 s, respectively. The first and second presentations of c(t) = c50% were considered as the unadapted and the adapted conditions, respectively.
In the final simulation experiment, contrast of input was varied as a rectified sine profile as in Ray and Maunsell (2010), their Figure 3. The contrast variation was simulated as follows:
Spectral analysis of model output
Since the model was used with stochastic inputs, each input combination was simulated for 50 repeats. Each repeat spanned from −0.547 to 3.548 s and was simulated by the forward Euler method with a step size of 2 × 10−5 s. An LFP proxy was computed from rE and rI as −rE−rI to emulate the recording depth expected to be closer to apical dendrite distribution of L23 [as in our previous studies, (Krishnakumaran et al., 2022; Krishnakumaran and Ray, 2023); see Krishnakumaran et al. (2022) for more discussion, and Mazzoni et al. (2015) for a computational study of LFP proxies in a morphological network model].
The TF spectrum of the LFP proxy for the contrast adaptation experiment was computed using MT in the same way as for the data, treating all 50 iterations of each stimulus condition as analogous to different trials from a single electrode.
However, for the rectified-sinusoidal contrast, the TF spectra were computed using matching pursuit (MP; Mallat and Zhang, 1993). MP was used so that our simulated TF spectra can be compared against the MP TF spectra of LFP data in Ray and Maunsell (2010), their Figure 3. To reduce computational load during MP, the LFP proxy signal recorded was low-pass filtered at a cutoff of 200 Hz using a fourth-order Butterworth filter and downsampled by 50, yielding 4,096 datapoints. The resultant signal was subject to MP analysis with a dyadic dictionary, which iteratively fits gabors to the LFP traces and computes the TF spectrum of the traces by adding together the Wigner–Ville distribution of all these gabors [see Ray et al. (2008) and Chandran et al. (2016) for details of the algorithm and the merits of MP analysis over other methods].
Results
LFP data was recorded from V1 in two awake macaques M1 and M2 in response to a full-screen achromatic grating presented in different contrasts consecutively without interstimulus interval (as shown in Fig. 1) from 77 and 32 electrodes, respectively.
Gamma frequency transient vanishes after contrast adaptation
Figure 2 compares the change in TF spectra of LFP corresponding to a 50% contrast when preceded by a blank gray screen (unadapted) versus by a 100% contrast (adapted). Figure 2A–D shows the change in TF spectral power from baseline corresponding to the unadapted condition from M1 and M2, respectively. Shortly after onset of 50% contrast grating at 0 s, there was an increase in power at all frequencies corresponding to the onset response, followed by gamma rhythm starting at a high frequency and gradually decreasing in frequency to a steady-state value over the initial 500 ms. This is a classic “frequency falloff” effect observed in many studies and modeled using slowly arriving inhibitory inputs in our previous study (Krishnakumaran and Ray, 2023). Figure 2B–E shows the change in TF spectra for the adapted condition with 100% contrast presented for 1 s before the 50% contrast. We hypothesized that in this case, both excitation and inhibition would be higher before the onset of the 50% contrast, and since the inhibition would take longer to dissipate, the frequency falloff should be abolished or reversed. Consistent with this hypothesis, we indeed observed that the gamma band for the adapted condition did not exhibit a frequency falloff in M1 (Fig. 2B) and was even mildly reversed in M2 (Fig. 2E). The peak-frequency trace of gamma was computed individually for each selected electrode and shown as median ± SE across electrodes in Figure 2C–F. Using a paired sign test to verify this suppression of gamma frequency transient at each timepoint, we found p < 0.01 for the entire duration in M1 and an initial ∼500 ms after the onset response in M2 after the onset response as indicated by the shaded areas at the top.
In Figure 3, we present the change in TF spectra from baseline for all contrast pairs for both monkeys. In the unadapted cases for all contrasts (left column), frequency falloff could be observed. For both 25% (top row) and 50% (middle row), increasing the contrast of the preceding (adaptor) stimulus caused a reduction and subsequently a slight reversal of the frequency falloff effect in both monkeys (columns from left to right).
Noisy JS model with delayed inhibition model replicates contrast adaptation transients
In Figure 4, we replicated the contrast adaptation experiment in the noisy JS model with delayed inhibition (Eqs. 1, 2) with 50 iterations. In Figure 4A, mean input drive to excitation iE and inhibition iI, respectively, were estimated using Equation 2, for the unadapted (left) and adapted (right) conditions. Due to the difference in θ between excitatory and inhibition inputs, the excitatory drive iE had almost instantaneous onset/offset response, following the contrast variations, whereas the inhibitory drive iI was slower and smoother with delayed response to stimulus onset. In Figure 4B, we show an example LFP proxy obtained from the noisy JS in one iteration, for each condition. After the lower contrast stimulus onset (marked as 0 s), in the unadapted condition (left), gamma rhythms were generated in bursts throughout the contrast presentation. Additionally, we found an initial evoked response in the baseline of the LFP proxy itself. In the unadapted condition, the nonzero stimulus onset evoked a negative deflection in the LFP proxy which disappears over time to settle at a steady-state value (higher baseline). However, in the adapted condition (right), the initial baseline deflection with respect to the later steady-state was in the opposite direction to that of the unadapted condition. In Figure 4C, the TF spectra averaged across 50 iterations of each condition are shown along with the peak frequency traces. We found an initial decrease in gamma frequency in the unadapted condition. However, in the adapted condition, the gamma frequency during the initial transient period was lower than the unadapted case. These differences in transients after adaptation were due to the slow-accumulating inhibitory inputs iI as adaptation with higher contrast gave rise to a higher initial drive to inhibitory populations after onset of the 50% contrast than in the unadapted case (Fig. 4A). In the unadapted condition, the slow increase of inhibitory input drive gave rise to the frequency falloff, whereas in the adapted condition, the slow decay of inhibitory drive produced the opposite trend.
Gamma response in noisy JS model with delayed inhibition to contrast adaptation. A, Mean input drives to excitatory population (IE, red trace) and inhibitory population (II, blue trace) estimated across 50 iterations for unadapted (left) and higher contrast-adapted (right) conditions. B, LFP proxy output from an example iteration from the model in each condition. C, Mean LFP time–frequency spectrum across iterations for each condition.
Gamma frequency profile in delayed inhibition model for smoothly varying contrast
In Figure 5, we present the output of the noisy JS with delayed surround model (Eqs. 1, 2) for rectified sinusoidally time-varying contrasts to compare against the macaque LFP data from Ray and Maunsell (2010), as replicated here in Figure 5E. It shows time-varying gamma frequency with asymmetric rise and fall trends at higher frequencies, with sudden rises and gradual falls in frequency around a contrast peak. In their experiment, an achromatic grating was presented in time-varying contrast profile, where the profile varied as a rectified sine function of time, emulated in Equation 4, with frequencies f = 0.625, 1.25, and 2.5 Hz. Figure 5A shows the mean excitatory and inhibitory drives emulating this experiment over 50 iterations. The excitatory drive traced the contrast variation closely, whereas the inhibitory drive rose and fell slowly resulting in delayed onset at lower frequencies and an accumulation across successive contrast crests at 2.5 Hz (Fig. 5B). Figure 5C shows example LFP proxies for the corresponding cases. For 0.625 Hz, although iE increased early, the LFP proxy did not show any activity until it reached a higher value. Figure 5D shows the TF spectra averaged across 50 iterations of each case. Gamma band increased in frequency slowly after the onset response for 0.625 Hz and showed asymmetric frequency rise and fall around the contrast peak time for higher frequencies for 1.25 and 2.5 Hz, as seen in Ray and Maunsell (2010), their Figure 3. Thus, this model explained the asymmetry in oscillatory dynamics during the rising versus falling phases, as observed in data.
Gamma response in noisy JS model with delayed inhibition to time-varying contrasts. A, Trajectory of mean inputs (iE, iI) across 50 iterations for rectified sine contrast modulation of 0.625 Hz (left), 1.25 Hz (middle), and 2.5 Hz (right) conditions. B, Mean input drives to excitatory population (red trace) and inhibitory population (blue trace) for corresponding conditions. C, LFP proxy output from an example iteration in each condition. D, Mean LFP time–frequency spectrum across iterations in each condition. E, Mean change in TF spectra of LFP recorded in macaque V1 during the presentation of the respective stimuli, replotted from Ray and Maunsell (2010), their Figure 3D.
There are multiple reasons that contribute to this asymmetry. First, due to slow accumulation/dissipation of inhibition, the model operates in different regimes during the rising versus falling phases, as shown in Figure 5A,B. In addition, in our experiment, we operated the model between nonoscillating and oscillating inputs for lower and higher contrasts, respectively. Figure 6A shows the gamma peak power and frequency generated at different input combinations in the JS model when the inputs are kept constant throughout the simulation. The inputs within the white contour correspond to the superlinear regime within which size and contrast effects on gamma are replicated by the model. Input combinations traced by iE and iI during the rising and falling phases in the simulation of 1.25 Hz sinusoidal contrast (Fig. 5A,B, middle column ) are indicated in magenta and black markers. The input trajectory starts outside the superlinear regime but enters the white contour in Figure 6A as contrast increases. Figure 6B depicts the phase diagrams and vector fields, for example, inputs emulating increasing contrasts (from left to right) during the rising phase of sinusoidally varying contrast with frequency 1.25 Hz (corresponding to the middle column of Fig. 5), to illustrate the change of dynamics across this transition, while Figure 6C depicts the same for the falling phase. In the rising phase of the stimulus (Fig. 6B), the contrast is initially very low [left column–input combination (1)] and there are two stable fixed points separated by an unstable fixed point in between. The system starts at a lower activity state (stable fixed point where rE and rI close to 0) during the initial part of the rising phase and remains there as the threshold of inputs required to cause the system to “jump” across the unstable fixed point is huge. As the contrast increases, the lower activity fixed point approaches the unstable fixed point [reducing the “jump” threshold to zero; Fig. 6B, middle panel–input combination (2)] until the two merge to form a saddle point. When the contrast increases beyond this “Saddle-node bifurcation” point, only one fixed point exists [Fig. 6B, right column; Fig. 6C, left column–input combinations (3–4)]. In our model, since the lower stable fixed point is close to 0, this results in the delayed onset of nonoscillatory activity yielding zero LFP proxy in the initial part in Figure 5C.
Dynamics in the JS model. A, Gamma peak power (left) and frequency (right) generated by the noiseless JS model (Eq. 1) with different combinations of constant inputs iE and iI. The white contour encircles input combinations within the “superlinear regime” identified by Jadi and Sejnowski (2014) as replicating the size and contrast effects on gamma power and frequency. The numbered magenta diamond markers indicate the input combinations (iE, iI) sampled at different timepoints during the rising phase of contrast in Figure 5A,B (middle column) derived from the delayed inhibition model (Eq. 2) subject to rectified sinusoidally varying contrast (Eq. 4) at f = 1.25 Hz. The numbered black circles mark example input combinations sampled during the falling phase of contrast. B, Phase diagrams of the JS model at example input combinations, corresponding to the rising phase samples (magenta diamonds in A; combinations numbered 1–3 in A are arranged from left to right). Red trace corresponds to the E nullcline which comprises pairs of (rE, rI) that result in
Meanwhile, the third fixed point [rightmost fixed point in left and middle panels in Fig. 6B–input combination (2)] has spiral trajectories attracting toward it. During the rising phase, this fixed point changes from stable spiral to locally unstable repelling the trajectories in its vicinity whereas the field farther from it is attractive, giving rise to a limit cycle at suitably high size and contrast (within the superlinear regime enclosed within the white contour in Fig. 6A). This transition from stable spiral to limit cycle by destabilization of the fixed point is the “Hopf bifurcation point.” As the system is driven across this bifurcation [from middle to right panel in Fig. 6B–input combination (3)], the amplitude of limit cycles increases from zero at the Hopf bifurcation point and increases. This gives rise to the delayed onset of oscillations in the rising phase where the delay is proportional to the time taken by the inputs to cross the bifurcation point. The frequency varies depending on the direction of the trajectory (Fig. 6A).
In the falling phase of the contrast (Fig. 6C), as the contrast reduces to lower values below the saddle-node bifurcation point, a saddle point emerging by a second tangential intersection of nullclines is observed in the middle panel in Figure 6C [input combination (5)]. The trajectory remains oscillatory due to the limit cycle around the unstable fixed point, remaining after the Hopf bifurcation [as observed in the orbiting vector fields encircling the neighborhood of the rightmost fixed point in the left and middle panels of Fig. 6C–input combinations (4–5)].
As the inputs fall further, a stable–unstable pair of fixed points emerge by separation of the saddle node, restoring the low-activity stable point. Simultaneously, the amplitude of limit cycles become larger (Fig. 6A) until its trajectory passes through the stable–unstable pair. At this stage, the previously limit cycle generating unstable fixed point behaves as an unstable spiral [Fig. 6C, right panel–input combination (6)], repelling the trajectories around it to slowly converge onto the zero-activity stable state. However, at this stage, especially when the newly separated stable–unstable pair is close by, a suitably strong noisy input fluctuation (say, a fluctuation in inhibitory input which shifts the I nullcline leftward) could cause the current state of the trajectory to “jump” across the unstable fixed point, resulting in a trajectory that loops around the rightmost unstable spiral point and returns slowly to the stable zero state, which we indeed observed in a few rare iterations. Hence for the same low contrasts, the system exhibits oscillations in the falling phase of contrast and zero activity in the rising phase.
Aside from the above bifurcations, we observed in the simulated trajectories of Figure 5, the rightmost panel in Figure 6C [input combination (6)] suggests that a third bifurcation may be possible for a slightly different parameter set or input parameters, called the “Saddle-node on invariant circle” (SNIC) bifurcation. In the middle panel of Figure 6C, instead of the saddle point emerging outside of the limit cycle as observed here, it may form on the trajectory of the limit cycle [for instance, see Pérez-Cervera et al. (2020), their Fig. 19], before transitioning into the state illustrated in the right panel in Figure 6C. The trajectory through the saddle point forms a loop over the unstable fixed point and in turn, the limit cycle generating fixed point now becomes an unstable spiral which becomes infinitely slower as it nears the saddle point. Hence, in the presence of noisy inputs perturbing the system about this saddle point, continuous oscillations with very low frequencies would be observed.
Thus, apart from the differences in iE and iI values during rising and falling phases as shown in Figure 5A,B, the emergence of bifurcations during contrast changes, the speed at which these dynamics unfold (which depends on the temporal frequency of the stimulus), location of the system relative to the fixed points all contribute to the overall dynamics of gamma oscillations.
Noisy JS model with slow-adaptive excitation replicates time-varying contrast effects
In Figure 7, we repeated the experiment in the noisy JS model with subtractively adaptive Excitation (Eqs. 1, 3) with 50 iterations. Figure 7A shows the mean input drives iE and iI computed using Equation 3, for the unadapted (left) and adapted (right) conditions. The slowly adapting excitatory drive iE had almost instantaneous onset in the unadapted case but decayed over time to a steady value. In the adapted case, the onset of the second lower contrast stimulus caused iE to momentarily decrease to a lower value and rise gradually to steady state. In Figure 7B, we show an example LFP proxy obtained in a single iteration, for each condition, exhibiting bursty gamma throughout stimulus presentation and initial evoked response in the baseline depending on the preceding contrast, as observed in the delayed surround model as well. In Figure 7C, the TF spectra averaged across 50 iterations of each condition are shown along with the peak frequency traces. We observe gamma frequency falloff as desired in the unadapted condition and its reversal in the adapted condition. The frequency transients in this model arise from the slowly varying adaptation signal which results in a transient reduction of iE in the unadapted condition and a transient increase in the adapted condition (Fig. 7A). The JS model responds to decreasing excitatory input drive with decrease in gamma frequency in the unadapted condition and vice versa in the adapted condition (analogous to the stimulus contrast effect on gamma frequency).
Gamma response in noisy JS model with slow adapted excitation to contrast adaptation. A, Mean input drives to excitatory population (IE, red trace) and inhibitory population (II, blue trace) estimated across 50 iterations for unadapted (left) and higher contrast-adapted (right) conditions. B, LFP proxy output from an example iteration from the model in each condition. C, Mean LFP time–frequency spectrum across iterations for each condition.
In Figure 8, we present the output of the noisy JS with adaptive Excitation (Eqs. 1, 3) for rectified sinusoidally time-varying contrasts as specified in Equation 4, with frequencies f = 0.625, 1.25, and 2.5 Hz to compare against the experimentally observed TF spectra (Fig. 5E). Figure 8, A and B, shows the trajectories and timeseries traces of mean excitatory and inhibitory drives emulating this experiment over 50 iterations. The inhibitory drive traced the contrast variation closely, whereas the adapted excitatory drive consistently lead inhibition with crests in excitation preceding those in inhibition slightly, exhibiting a similar trajectory shape and IE-leads-II phase relationship as in Figure 5A,B. Figure 8C shows example LFP proxies for the corresponding cases, while Figure 8D shows the TF spectra averaged across 50 iterations of each case, exhibiting the asymmetric frequency rise and fall around the contrast peak time, as seen in Ray and Maunsell (2010), their Figure 3 (replotted in Fig. 5E) and in the delayed surround model. Here, the asymmetry arises because of iE values being different during the rising and falling phases of contrast due to adaptation (as shown in Fig. 8A,B) and the differences in dynamics related to the emergence of bifurcation points, etc. as discussed earlier. Thus, a model with adaptation in feedforward inputs could also explain our results provided the excitatory input adaptation exhibits slow dynamics equivalent to the time constant of
Gamma response in noisy JS model with adapted excitation to rectified sinusoidally varying contrasts. A, Trajectory of mean inputs (iE, iI) across 50 iterations for rectified sine contrast modulation of 0.625 Hz (left), 1.25 Hz (middle), and 2.5 Hz (right) conditions. B, Mean input drives to excitatory population (red trace) and inhibitory population (blue trace) for corresponding conditions. C, LFP proxy output from an example iteration in each condition. D, Mean LFP time–frequency spectrum across iterations in each condition.
Discussion
We show that adaptation with a higher contrast grating abolishes or reverses the transient reduction in gamma frequency after onset of a grating. We show that this difference in gamma frequency transient after adaptation could be replicated in the noisy JS model by virtue of its delayed inhibitory input drive, which also replicates the variation in gamma onset and frequency modulation to stimuli with time-varying contrasts. Another model in which only excitatory inputs are adapted could also explain these results.
Visual adaptation
Typical experiments for studying V1 circuitry involve sustained visual stimulation with constant, or “regularly” varying stimuli, optimized for V1 receptive fields, whereas realistic stimulation and the corresponding responses are less parametrized but have been less tangible to inference so far. Stimulus adaptation studies have shown that neuronal responses are dependent on the sequence of stimulus presentations, at every stage of the visual pathway from retinal ganglia to V1 [see Kohn (2007) for a review] and exhibit correlates to perceptual aftereffects (Jin et al., 2005). Prior-dependent response to subsequent stimuli is ubiquitous in the sensory cortex, resulting in theories, such as predictive-coding or energy efficient-coding strategies (Webster, 2015; Weber et al., 2019), treating adaptation as a signature of the underlying perceptual-coding principle.
In V1, adaptation to prior exposure to achromatic stimuli (adaptors) reportedly result in temporary suppression of firing rate responses to the test stimuli with different recovery timescales depending on the duration of exposure to the adaptor [Patterson et al. (2013), their Figs. 7, 8; Fig. 8 shows adaptation recovery time constants of the order of 1 s for a 4 s exposure to adaptor]. Further recordings from areas 17 and 18 in anesthetized cats exhibit different extents of contrast adaptation in response to 2 s exposure to ramping contrasts with different frequencies (Crowder et al., 2008). Within this time, temporary changes to orientation tuning are observed during subsequent stimulus presentations in the neurons recorded (Patterson et al., 2013; Ghodrati et al., 2019) whereas adaptation in LGN only causes gain changes and not tuning changes (Dhruv and Carandini, 2014). The extent of these orientation changes postadaptation themselves vary in a graded fashion depending on the orientation of the adaptor stimulus, suggesting a significant contribution from within V1 [Wissig and Kohn, 2012; Patterson et al. (2013), their Fig. 9]. Further, Westerberg et al. (2019) studied adaptation effects in laminar recordings from V1 and reported a profound change in activity within the supragranular layers, where intracortical connections are predominant, and suggest that the observed adaptation-driven decrease in synaptic activity occurs mainly outside the input layer where feedforward inputs arrive. On the other hand, contributions to surround suppression coming from lateral unmyelinated connections are orientation selective (Angelucci et al., 2017), which is emulated in inhibitory drive dynamics in our model (see below, Delayed surround drive to inhibition, for elaboration). The feedforward inputs to V1 from LGN are bound to exhibit temporary suppressive dynamics due to contrast adaptation themselves (Camp et al., 2009; Raghavan et al., 2023), as well as its own integration time with contrast (Mante et al., 2005). To emulate this, in our adaptive excitation model, we explicitly introduce these adaptation effects in the LGN-driven input drive iE which also introduces the desired frequency falloff and adaptation effects. In both our input models, a slow timescale for inputs (time constant of 1/(2π1) s = 159 ms) was required to match the data. Adaptation in V1 has been observed to last from brief to long durations, depending on the duration of exposure to the adaptor [Patterson et al. (2013), their Fig. 7]. However, this dependence of adaptation timescales on adaptor duration has not been explicitly incorporated yet in either the delayed inhibition or the adaptive excitation models formulated in this study, both of which employ simple first-order linear filters with fixed, slow time constants to emulate the respective phenomena.
Delayed surround drive to inhibition
The noisy JS model replicates gamma frequency trends by setting θI << θE in the input model, as the lower value of θI delays the change in surround drive to inhibition. Delayed inhibition over millisecond timescales is observed in cortical activity and emerges intrinsically in PING models [Krishnakumaran and Ray (2023); see Discussion for review]. Typically, a local neuronal population as in the noisy JS model is thought to be strongly “local” in interconnectivity (within a cortical column), but they may send lateral inputs to similarly tuned surround populations. Such lateral inputs are more convergent onto inhibitory neurons, mediating surround suppression (Shushruth et al., 2012). Further, these surround inputs arriving through unmyelinated axons are expected to have a considerable latency depending on distance between source and target populations. Bair and colleagues (Bair et al., 2003) analyzed surround suppression in macaque V1, produced by far and near surround stimuli, and found two distinct suppression components with different timescales—a slower component thought to arise through lateral connections at latencies up to ∼20 ms and a fast component arriving from feedback from higher areas with much less delay [Bair et al. (2003), their Fig. 7A; see Angelucci et al. (2017) for detailed review of different surround inputs]. Therefore, while the incoming excitation to a particular area responds with lower latency to the feedforward projection from the immediately upstream area (LGN in case of V1 neurons), the inhibition which is activated by the local and widespread surround network with distance-dependent latencies (Girard et al., 2001) will be both delayed and smoothed out. Stimulus-selective surround suppression indeed exhibits delayed onset (∼10 ms after evoked response onset) in LFP recordings (Angelucci et al., 2017; Wang et al., 2020). Such slow dynamics of inhibition also yield transients in nonoscillatory activity. In recent studies (Zhou et al., 2019, 2023), a delayed normalization (DN) model has been fitted to data from intracranial EEG under different stimulus conditions. The DN model has a different structure than the noisy JS model but also incorporates relative delays in excitation and surround-driven inhibition using different time constants. In these studies, fitting data to DN model yielded a surround time constant of over 2× the excitatory time constant, arguing for delayed effect of surround suppression contributing to the transient trends in nonoscillatory activity.
Inhibition network in V1
Interneurons in V1 are thought to be driven mostly by intracortical excitatory populations than by feedforward inputs to V1. This excitation-driven inhibition is evidenced by different surround modulation and transient properties (Ozeki et al., 2009; Shushruth et al., 2012) typically observed in recordings from cats and macaques, and serves as the principle behind all variations of PING models (Buzsáki and Wang, 2012) and ISN models (Tsodyks et al., 1997; Jadi and Sejnowski, 2014) such as the JS model used in this study. Most network models simulated gamma rhythms generated by interactions between homogeneous populations of excitatory and inhibitory neurons (Chariker et al., 2018; Zachariou et al., 2021). However, recent studies in mice have shown distinct interactions of different interneuronal subtypes (Adesnik et al., 2012; Veit et al., 2017). Network models incorporating these details may help study their roles in cortical processing (Wagatsuma et al., 2023). Two major interneuronal networks involve parvalbumin (PV)-expressing neurons and somatostatin (SOM) and vasoactive intestinal peptide (VIP)-expressing neurons. The PV and SOM→VIP neuronal networks differ in spatial connectivity ranges (Adesnik et al., 2012), with SOM neuron activities showing longer-range spatial integration than PV. Further, while PV neurons synapse onto a pyramidal cell near the soma, SOM neurons’ synapses are found at more distal regions. The synapse location difference has been demonstrated, by morphological modeling (Headley et al., 2024), to affect the nature and extent of the interneurons’ participation during rhythmic activity of different frequencies. Further the wider-ranged SOM network could be implicated in emergence of a second, slow gamma of lower frequency range, with different orientation selectivity [Murty et al. (2018), their Fig. 1] and a wider spatial coherence profile than previously observed fast gamma, during larger-size stimulus presentations (Murty et al., 2018). While our current model does not generate slow gamma, it has been modeled either using a second inhibitory population with different dynamics (Keeley et al., 2017) or simulating an altered effective connectivity caused by larger stimuli activating farther networks interacting with the local network by wider horizontal connections (Han et al., 2021).
Notes on the JS model
JS model was chosen for its self-oscillating behavior across a supercritical Hopf bifurcation, exhibiting V1 gamma-like trends in power and frequency. The model represents a local population of strongly interconnected excitatory and inhibitory neurons with similar stimulus selectivity. The strong recurrent excitatory causes the model to operate in ISN regime (Tsodyks et al., 1997; Jadi and Sejnowski, 2014). It is expected that such a local population be uniformly interconnected so that a mean-field approximation could be made under uniform stimulus conditions. Larger stimuli evoke activity in surrounding populations, which also drive the local population along with the LGN feedforward input from the stimulus center. This surround input is expected to be stronger for the wide-connecting inhibitory neurons than the strongly locally recurrent excitatory neurons. Hence, JS model emulates larger-size stimuli using higher values of input drive to inhibition II (Jadi and Sejnowski, 2014). The bifurcation property and its sensitivity to recurrent-connectivity parameters also allow it to emulate the attenuation of gamma power by different discontinuities in the full-screen grating (Shirhatti et al., 2022). Shirhatti and colleagues modeled discontinuity effects as a change in recurrent connectivity among the activated local population. This could result from either some neurons becoming inactive due to the discontinuity or from the segment of discontinuity recruiting a different set of neurons at different spatial locations which have less interconnectivity than the population activated by a uniform grating. The resultant decrease in recurrent connectivity was shown to shrink the bifurcation boundary enclosing the oscillation-generating inputs. Hence a certain input would be out of the oscillatory regime in the discontinuous condition. Further, the nonlinearity of the model produces nonsinusoidal gamma waveforms. In Krishnakumaran et al. (2022), we found that the waveform shape found in LFP could be replicated by a subset of oscillatory inputs by the JS model.
Another WC-type model, we refer to as the Jia–Xing–Kohn (JXK) model, had been demonstrated to generate bursty gamma rhythms in the form of quasicycles when subject to temporally uncorrelated Poisson noise as inputs and exhibits size and contrast dependencies (Jia et al., 2013) but fails to generate consistent nonsinusoidal waveforms (Krishnakumaran et al., 2022). The JXK model accounts for wider cortical excitation by larger stimuli using a third global excitation population whose influence is larger for larger stimuli. The JS model, in contrast, specifies only local excitatory and inhibitory populations. The activation of surrounding populations during presentation of larger stimuli is emulated by increased surround suppression drive to the inhibitory population (Jadi and Sejnowski, 2014). In the delayed inhibition model, the timescale of local inhibition is determined by τI in Equation 1, whereas that of the global surround suppression arriving through unmyelinated connections is determined by the time constant
As inputs in the original JS model were constant over time, the model would give rise to sustained gamma oscillations as opposed to the bursty gamma observed in LFP. Furthermore, response of the model to stimulus change was rapid and lacked the transients observed in data. Gamma bursts and transient trends in frequency have been replicated by driving the model with OU-type noise in Krishnakumaran and Ray (2023) and in this study. While the core model and its interconnectivity has not been changed, the rich temporal properties of gamma rhythms emerge naturally when realistic input dynamics are incorporated.
Still, the model's specific nonlinearity poses a hindrance to modeling the realistic gradual decrease in activity expected when contrast is slowly increased, as observed in Figure 5C for 0.625 Hz stimuli. Also, between nonoscillating and oscillatory values of IE, the model produces widely deviating trajectories (some exhibiting oscillations with declining frequencies while others with abrupt reduction to zero-activity state) owing to multiple fixed points, depending on noisy fluctuations in the inputs, contributing to the irregular low-power bursts of oscillation in the TF spectrum in Figure 5D for 2.5 Hz in the interval around 0.2–0.3 ms instead of a smooth sustained decline in frequency. Further, since the low-activity fixed point in Figure 6B,C is at zero at any input combination it emerges in, the activity of the system is zero at low contrasts and sharply rises to a different value after a bifurcation point. A smoother transition may be brought about if the nullcline shapes were modified to allow a slowly-shifting stable fixed point instead. But this would indeed have side effects on the dynamical properties of the model. Overcoming the shortcomings of sharp decline to zero-activity state and deviating trajectories may require abandoning the existing, simple nonlinearity of the model to accommodate realistic behavior outside the oscillation-generating input domain. However, such a data-driven alteration of the model's nonlinearity has not yet been attempted.
Footnotes
We thank Surya S. Prakash, Ankan Biswas, and Divya Gulati for their assistance in data collection. We thank Prof. Adam Kohn for valuable discussions and insights. This work was supported by Wellcome Trust/DBT India Alliance (Senior fellowship IA/S/18/2/504003 to S.R.) and a grant from Pratiksha Trusts.
The authors declare no competing financial interests.
- Correspondence should be addressed to Supratim Ray at sray{at}iisc.ac.in.