Sensory Stream Adaptation in Chaotic Networks

Implicit expectations induced by predictable stimuli sequences affect neuronal response to upcoming stimuli at both single cell and neural population levels. Temporally regular sensory streams also phase entrain ongoing low frequency brain oscillations but how and why this happens is unknown. Here we investigate how random recurrent neural networks without plasticity respond to stimuli streams containing oddballs. We found the neuronal correlates of sensory stream adaptation emerge if networks generate chaotic oscillations which can be phase entrained by stimulus streams. The resultant activity patterns are close to critical and support history dependent response on long timescales. Because critical network entrainment is a slow process stimulus response adapts gradually over multiple repetitions. Repeated stimuli generate suppressed responses but oddball responses are large and distinct. Oscillatory mismatch responses persist in population activity for long periods after stimulus offset while individual cell mismatch responses are strongly phasic. These effects are weakened in temporally irregular sensory streams. Thus we show that network phase entrainment provides a biologically plausible mechanism for neural oddball detection. Our results do not depend on specific network characteristics, are consistent with experimental studies and may be relevant for multiple pathologies demonstrating altered mismatch processing such as schizophrenia and depression.

In predictably structured sensory streams like music, speech and visual input during motion the brain doesn't simply respond to external events but anticipates them [1][2][3][4] . Studies of temporal expectation, also known as implicit timing 5 , demonstrate the importance of predictability of when stimuli will occur. Task performance is enhanced if stimuli occur at expected times according to an established rhythm [6][7][8] . The sensory streams used are also found to phase entrain low-frequency brain oscillations [9][10][11][12] . Phase entrainment increases with the temporal regularity of the sensory stream and correlates with behavioural performance. It is not limited to sensory cortices 13 ; activity patterns in larger scale cortical 14,15 and subcortical networks 16 are also modulated by temporal expectation.
Related studies show that implicit expectations of upcoming stimulus itself can also strongly affect neural response. Stimulus specific adapatation (SSA) occurs in single cells 17 and in multi-unit activity 18 . Here the response to a stimulus is in general suppressed after multiple presentations while responses to oddball stimuli may be enhanced 19 . SSA occurs in both excitatory and inhibitory cells 20,21 and in multiple areas from midbrain to cortex 22 . Type expectation also occurs at the neuronal population level as oscillatory mismatch responses (MMR) after stimulus offset in EEG/MEG [23][24][25][26][27] , ECoG 19 , and LFP 28 . Most importantly adaptation can occur to whole stimulus sequences on multiple time scales simultaneously 18,29 while the effects are stronger when streaming is temporally regular 30 . Mismatch responses are altered in pathologies, notably schizophrenia, traumatic brain injury and Alzheimers disease, and can be used as diagnostic tools [31][32][33][34][35] .
Whether population and single cell modulations originate from the same underlying mechanism is controversial [36][37][38][39] . SSA is thought to depend on single cell adaptation but studies have also implicated network mechanisms 21,40,41 including recurrent inhibition 42,43 . MMR has been suggested to involve synaptic adaptation 44 as well as memory 45 and predictive coding 46 based mechanisms. Intriguing recent investigations have highlighted a possible interaction of SSA and MMR through recurrent network inhibition 20,47 .
Neural information processing may be embedded in recurrent network dynamics [48][49][50][51][52][53] which are also thought to underlie ongoing brain oscillations [54][55][56] . Recent modeling work 57 has demonstrated entrainment of low frequency stable network oscillations by transcranial magnetic stimulation. However low frequency brain oscillations may also be chaotic [58][59][60][61][62] . We find that the neuronal correlates of sensory stream adaptation emerge when and sensory input by controlling network stability 59 . As a simple proxy for the LFP/EEG we use the total absolute size of the average input current to a rate unit 57,70,83,84 , which we term the 'neuroelectric activity' 57 .
The sensory streams we investigate are made up of several stimuli applied in sequence. A stimulus represents a static visual noise patch or an auditory stimulus composed of random frequencies, for example. Each stimulus is modeled as a fixed set of excitatory driving inputs (see Methods) one of which is applied to each of the units in the network for the duration of the stimulus. The driving inputs which make up any given stimulus are random, drawn independently from a single given distribution once and then fixed for the duration of the simulation. Therefore all stimuli are identical statistically (see Methods).
Most of the sensory streams we investigate include only two sensory stimuli, denoted A and B. They are always applied to the network for brief 50 ms periods, as is common in sensory streaming experiments. During inter-trial-intervals (ITIs) between stimuli the network is driven by a different input we denote 'background' input X. Input X can represent external sensory input during ITI periods or alternatively purely internally generated excitatory driving or a combination of the two. In order to demonstrate that our results do not require that sensory stimuli provide stronger levels of excitation than the background input we also draw input X from the same distribution as the sensory stimuli A and B. The ITIs in a sensory stream can be of fixed or variable length but the mean length investigated is usually 800 ms. This generates sensory streams with stimulus presentation frequencies typically used in experiments.
Simulations (except those for demonstration of dynamics and for calculations of Lyapunov exponents (see below)) include noise. The noise is not necessary for the results here. It has been included to show that our results are resilient to its presence (see Methods). In the first part of this paper we illustrate how chaos entrainment works to generate reliable, distinct and history dependent stimulus response in this model. In the second part we show how these properties can explain several experimental results in sensory stream adaptation.
Sensory streams entrain chaotic network activity. First we illustrate how sensory streams entrain chaos in a particular exemplar 500 unit inhibitory network simulation. Figure 1(A,B) show the spontaneous activity generated by the network when it is constantly driven by the background input X. (In the following we refer to the activity endogenously generated by the network under the constant driving of the background input X as the 'autonomous' activity.) The random looking oscillations in the mean neuroelectric activity ( Fig. 1(A)) and in the individual units activity ( Fig. 1(B)) suggest that the autonomous activity is chaotic.
We use principal component analysis to visualize the network chaotic attractor. The grey line in Fig. 2(A) shows the trajectory of the autonomous activity projected in the space of its first two principal components. One particular trajectory is shown. However there are multiple other trajectories winding around the attractor, tangled up with the one shown. Which trajectory the system follows depends on the choice of initial condition. Each trajectory gradually fills the attractor as time progresses but trajectories also diverge from each other on shorter timescales. Due to this mixing property initially nearby points eventually become spread across the whole attractor.
This trajectory divergence is quantified by the maximal Lyapunov exponent. It measures the rate at which nearby dynamical trajectories diverge or converge. When the maximal Lyapunov exponent is positive network activity is chaotic and unstable to small perturbations, like noise, which can knock the system from one trajectory onto another. On the other hand in networks with stable dynamical attractors, like limit cycles and fixed points, small perturbations to trajectories decay and the maximal Lyapunov exponent is negative or zero. In this network simulation example (Figs 1(A) and 2(A, grey line)) the maximal Lyapunov exponent of the autonomous activity, denoted λ U , (see Methods) is positive, λ U = 0.0012, confirming the presence of chaos. Figure 1(C,D) shows what happens when the same network simulation is driven by a particular sensory stream, denoted the AB sensory stream. Here the two 50 ms stimuli, A and B, are presented in alternating sequence. The ITI following stimulus A is always 415 ms while that after B is always 1185 ms. Evidently the irregular variations in activity in the autonomous network simulation ( Fig. 1(A,B)) become locked to the sensory stream in the driven simulation. The solid lines in Fig. 2(A) show the sensory stream driven trajectory projected onto the principal components of the autonomous activity. For visualization the trajectory has been divided into different sections coloured according to the epoch in the sensory stream. The trajectory appears to circulate in the vicinity of the autonomous attractor and is periodic with period of the sensory driving (here and in the following we use the terminology 'periodic sensory stream' to denote a sensory stream which repeats exactly at fixed time intervals). This driven network trajectory is fully stable as confirmed by calculation of the driven maximal Lyapunov exponent, denoted λ D (see Methods) which is negative, λ D = −0.0013, for this simulation in this particular sensory stream. This example therefore demonstrates the complete stabilization of autonomous rate chaos by the sensory stream.
Entrainment of chaos can also be partial. Figure 2(B,C) show examples from the same network simulation under the alternating AB sensory stream with the same stimuli A and B except that the ITI lengths have been changed slightly, while the mean ITI remains fixed at 800 ms. In Fig. 2(B) the driven trajectory is again fully stable, λ D < 0, but it is now 'mode locked' with double the period of the driving, a behaviour also shown by stable network oscillators 57 . In Fig. 2(C) chaos is only partially stabilized, λ U > λ D > 0. Although the driven trajectory never quite repeats it is strongly restricted within a 'bundle' of much lower dimension than the chaotic attractor (grey). This is a general characteristic of chaotic synchronization 75,76,[85][86][87] and, if the sensory stream is periodic, phase synchronization 63 . Here phase synchronization is evident because any point in the stream occurs roughly at the same fixed phase around the trajectory on each cycle.
Stabilized chaotic trajectories are strongly history dependent. In sensory stream adaptation experiments repeat presentations of a given stimulus must generate different distinct responses which are however reliable when the stream itself is repeated 51,52,74,77,88-91 . Here we illustrate why (fully and partially) stabilized chaos provides the ideal substrate for this 92,93 . Figure 2(D) shows the trajectory of the same network driven by the same stimuli A and B but now composed in a different sensory stream where the sequence AAABB is presented repeatedly with fixed 800 ms ITIs between all stimulus presentations. The trajectory is again stable, λ D < 0, but stimulus responses are totally different to those under the AB stream, Fig. 2(A). The response to any given stimulus is strongly history dependent and adapted to the periodic sensory stream. Three distinct stimulus A responses and two distinct stimulus B responses can be seen. Such history dependent neural response is known in various task settings and brain regions. 94 provides a recent example.
The distinctiveness of stimulus responses depends on how chaotic the autonomous network dynamics is, given by λ U . Points on trajectories which are preceded by more similar sensory streams tend to be closer in the state space. This is visible in Fig. 2(D) where the stimulus response endpoints (green circles) for the second and third presentations of A are closer together than they are to the first stimulus A endpoint. However as the strength and dimension of autonomous chaos increases ITI period trajectories become more divergent and the driven trajectory becomes spread over a larger region. Therefore the responses generated by individual stimulus presentations reflect longer and longer sequences of preceding stimuli.
On the other hand the reliability of a stimulus response depends on two factors. First, if chaos is not completely suppressed, Fig. 2 Fig. 2(B), repeat presentations of the same stimulus can generate slightly different responses, even though the dynamics is completely deterministic. This type of reliability is described by λ D and in general decreases as it increases. It is evident from these examples that the effects of unsupressed chaos on response reliability can be very weak. This is because the 'widths' of the trajectory 'bundles' are small compared to the total variation around the trajectory phase. In fact partially stabilized chaos is particularly relevant for sensory stream adaptation since it combines fairly high reliability with high distinctiveness (as illustrated below.) A much stronger effect on stimulus response reliability originates directly from the sensory stream itself. In fully and partially stabilized chaotic networks, λ D < λ U , response reliabilities for given stimuli decrease as the multiplicity of sensory streams which precede individual presentations of that stimulus increases. For example in the AAABB sensory stream, Fig. 2(D), stimuli of type A generate much less reliable responses than stimuli of type B.
Chaos stabilization only weakly depends on stream predictability. Since stabilized chaos facilitates sensory stream adapted response it is important to investigate how chaos stabilization depends on network properties as well as on the properties of the sensory stream.
In Fig. 3 driven Lyapunov exponents λ D are plotted versus autonomous Lyapunov exponents λ U for many inhibitory network simulations generated with various network connection probabilities (see Methods.) Fig. 3

(A)
shows λ D for streams where stimuli A and B are presented in alternating sequence separated by fixed and equal ITIs, for various ITI lengths. Chaotic network activity is generally stabilized somewhat by the sensory driving since λ D < λ U when λ U > 0. The magnitude of stabilization, λ U − λ D , is largest for networks with λ U >≈ 0 and decreases to zero as networks become more strongly chaotic. The higher the stimulus presentation frequency the greater the stabilization, as also found in 95,96 . We find some stabilization occurs for ITI periods as short as 100 ms and as long as several seconds. These are natural timescales often employed in sensory streaming studies 6,9,10,18,19,25,[98][99][100][101] . Figure 3(B) compares λ D for the alternating AB sensory stream with fixed 800 ms ITI and a maximally irregular sensory stream where ITI intervals are exponentially distributed with mean 800 ms and where stimuli A and B are applied in scrambled sequence. Although for any particular simulation λ D can depend on the stream (as demonstrated by the example in Fig. 2), on average both periodic and random streams stabilize chaos roughly equally ( Fig. 3(B) lines), (although there is a weak effect of randomness). The somewhat surprising finding that even completely random streams suppress chaos as much as periodic ones is important because it means stimuli in random streams also have the capacity to show stream history dependent response. If this were not the case stimulus response would become fully unreliable suddenly at some level of stream randomness. The point is illustrated in Fig. 2(E,F) where the trajectories under the maximally irregular and a weakly irregular sensory stream are compared. Stimulus responses and ITI period trajectories are much more reliable when activity is less irregular. Since both trajectories are almost equally stable relative response reliabilities are determined by the regularity of the sensory stream. There have been multiple investigations of chaos control in neural networks, for example 53,77,92,95,97,102 . Here we don't address it in detail. However simulations show large negative deflections in λ D at stimulus onset and offset (see Supplemental Fig. S8 and Supplemental Text). This suggests that the strong dependence of stabilization on stimulus presentation frequency and weak dependence on randomness may be very roughly explained by a simple mechanism. When the autonomous Lyapunov exponent, λ U , generated when the network is driven by the background input X, is small there are only a few unstable directions which are confined within a low dimensional attractor manifold within the much higher dimensional space. During stimulus presentations the input X generated attractor is replaced by a different stimulus dependent one. This causes the activity state to move away from the X attractor region towards the new attractor region and then rapidly return after stimulus offset. Both these transient attractor switching pieces of the driven trajectory are stable. If the stabilizing effect of these kicks is sufficient to overcome the divergence on the X attractor then the driven trajectory will become stable overall, regardless of the temporal regularity or type of sensory stimulus presentations.
In the following we denote the approximate regime where chaos is fully or partially stabilized, 0 < λ U < 0.003, for sensory streams with mean ITI around 800 ms, typically used in experiments the 'suppressed chaos regime' . Here we focus on the role of network spontaneous activity in sensory stream adaptation, rather than how this activity arises from network structure. However for completeness in Supplemental Fig. S1 we investigate this regime in more detail. It is quite broad in biologically realistic regimes of network parameter space, spanning recurrent connection probabilities from about 0.1 to 0.25, while networks including excitatory units are inhibition dominated. The point where chaos is just suppressed in fully inhibitory networks has connection probability around 0.17 and mean synpatic conductance around 0.47 nS/cm 2 when average firing rates are 1 Hz (see Methods and Supplemental Methods.)

Reliable and distinct deviancy level dependent response in roving sensory streams.
Multiple studies show that a given stimulus produces different neuronal responses depending on whether the stimulus appears as implicitly expected according to the previous streaming 7, 18,19,[23][24][25][26][27][28][29][30]39,99,101,[103][104][105][106] . To investigate this we employ a 'roving' sensory stream, (Fig. 4), often used in experimental studies. Again two stimuli, A and B, are used but each one is repeated several, 'n' , times before switching to the other. The first stimulus after a switch is the most 'deviant' of that stimulus type, denoted 'D1' , the second after a switch is the second most deviant, denoted 'D2' , and so on. The final stimulus before a switch, 'Dn' , is the 'standard' , also denoted 'S' . (We also sometimes use 'D(n-1)' and 'D(n-2)' as the standard (see below.)) Roving sensory streams switch between the two 'single stimulus sensory streams' which we denote AA and BB.
We first illustrate the population activity. Figure 5(A) shows differential peri-stimulus time histograms (dPSTH) of the neuroelectric activity for the network simulation investigated above (Figs 1,2) under the n = 4 roving sensory stream with ITI 800 ms. dPSTH are generated by subtracting the activity generated by the standard S (D4) presentation of a given stimulus type, A or B, from the activity generated by deviant presentations of the same stimulus type (see Methods). Each stimulus type is shown separately. D1 stimuli show strong differential responses at stimulus presentation and modulations which persist throughout the ITI period.
Differential modulations are highly significant, as can be seen from the small error bars. This is because individual presentations of a given stimulus at a given deviance level are always preceded by identical sensory streams in this roving paradigm and because chaos is stabilized, λ D < 0 < λ U . Given stimuli therefore generate distinct and reliable response at different deviancy levels. This is readily appreciated from the sensory stream driven activity projected in the principal components of the autonomous activity, Fig. 5(B).
Oscillatory differential ITI period modulations, (Figs 6(A,B) and 7(A)) which are qualitatively as observed in population level studies 19,23-27 occur because the low dimensional chaotic dynamics is dominated by low frequency oscillatory components. As illustrated in Fig. 5(B) D1 deviants produce responses which are relatively close to, but shifted away from the standards. This results in deviant and standard ITI period trajectories which rotate in a roughly concentric way generating the differential oscillations in the neuroelectric activity dPSTH, Fig. 5

(A).
Responses adapt slowly in the suppressed chaos regime. Experiments show that differential responses adapt gradually over surprisingly large numbers of stimulus presentations (e.g 19,25 .). Evidently as the number of repeat presentations of a given stimulus, D2, D3, …, increases, the activity adapts to the trajectories Furthermore due to the bundling of trajectories described above we find a similar gradual relaxation process occurs in the regime where chaos is only partially suppressed. Figure 6(A) shows neuroelectric activity dPSTH from a different purely inhibitory network simulation driven by the n = 12 roving sensory stream while Fig. 7(A) shows an example from a network including both excitatory and inhibitory units and in which only half of the units are directly driven by the sensory stream. In both cases driven activity is only partially stabilized, λ D > 0, and adaptation occurs slowly over five or six stimulus presentations. These observations are summarized in Fig. 8(A) which shows the absolute size of differential modulations, dPSTH , (|..| denotes absolute value and 〈..〉 denotes average) averaged across the whole 850 ms stimulus and ITI period for various levels of deviancy in inhibitory networks driven by the n = 12 roving sensory stream. Only in the suppressed chaos regime do deviants generate large differential responses, peaking close to where chaos is just stabilized, λ D = 0. The strongest dependence on repetition number is also shown in this regime.
The above example indicates the following explanation for these results. As λ U increases responses at different levels of deviancy become more distinct and differential fluctuations (signal) larger. When λ D increases past zero unreliability due to unsuppressed chaos (noise) increases. Because the dPSTH are calculated from standardized neuroelectric activity time series, dPSTH reflect the relative size of signal to noise which decreases as the chaotic noise increases. Differential dPSTH depend on deviancy level throughout the suppressed chaos regime because stimulus response continues to reflect the preceding sensory stream even when chaos is only partially suppressed.
Differential modulations decrease in magnitude when streaming is irregular. Differential responses are also known to depend on the temporal regularity of the adapting sensory stream 30 . Deviants in irregular sensory streams generate differential responses of smaller magnitude than their counterparts in regular streams. Here we investigate this using an irregularly timed roving sensory stream (see Fig. 4). Differential responses to deviants under the irregularly timed stream (Fig. 6(B)) have similar form to their counterparts under the regularly timed stream (Fig. 6(A)) but are weaker and the error bars larger. This is confirmed in Fig. 8(B) which shows the magnitude of differential modulations, dPSTH , versus λ U for inhibitory network simulations under the irregularly timed n = 12 roving sensory stream. The reduced size of the D1 differential response compared to Fig. 8(A) again occurs because the relative size of the signal component is reduced compared to noise. However in contrast to the case of temporally regular streaming described above where unsuppressed chaos was the sole source of unreliability, here unreliability is also caused by variability in the preceding streaming even when chaos is fully suppressed, λ D < 0. In particular the increased variability of the S response (compare D7 in Fig. 8(A) and (B)) prevents adaptation after the third repetition, D3, and reduces the size of the D1 differential response. These observations are in qualitative agreement with the experimental studies (see Discussion.)

Time-frequency responses.
Multiple studies show that sensory streams phase entrain brain oscillations at the streaming frequency, and that this is weakened by temporal irregularity, for example [9][10][11][12][13][14][15][16]107 . Standard and deviant presentations of stimuli are also known to show characteristic responses 31,33,107,108 and differences 19,35,109,110 in their neuroelectric activity spectrograms. Stimulus presentations sometimes generate bursts of coherence in the theta and alpha bands associated with phase resetting of endogenous oscillations 35,108 . Here we investigate these issues in the n = 12 roving sensory stream. We employ two manipulations to remove the evoked stimulus response and reveal the underlying 'induced' activity.
The first 'weaker' manipulation is sometimes employed in experimental studies. From each time point τ ms after the onset of a given stimulus the mean value for that time point τ averaged across all presentations of that stimulus is subtracted (see Methods). We term this the 'reduced' activity. Supplemental Fig.S3(A) shows inter-trial phase coherence (ITC) across trials of a given deviancy level D of this reduced activity in a narrow band around the streaming frequency 1/0.85 Hz under the regularly and irregularly timed n = 12 roving streams versus the autonomous Lyapunov exponent λ U . Phase coherence peaks around the point where chaos is just stabilized for all levels of deviancy and is weaker when streaming is irregular. However it has a non-monotonic dependence on deviancy level D. This is because phase coherence is mainly due to the stimulus response which, due to response adaptation, is not fully removed by this manipulation. Indeed phase coherence at first decreases with increasing deviancy and then increases again, as expected if the D1 response is larger than the subtracted average and the D12 response smaller than the subtracted average. Phase coherence is however very low in stable networks demonstrating their absence of adaptation. , show strong bursts of excess coherence and power locked to stimulus presentation, while these effects are weaker in more strongly chaotic and in stable networks (note different colour scales.) Such responses, which last several hundred milliseconds after stimulus offset and can extend up to around 80 Hz, are sometimes seen in experimental studies 19,31,33,35,[107][108][109][110] .
Of more interest in the current study is a second 'stronger' manipulation. Here from each time point τ ms after the onset of a given stimulus of given deviancy level D we subtract the mean value for that time point τ averaged across all presentations of that stimulus of that deviancy level D. We term this the 'induced' activity. It reveals the underlying chaotic 'noise fluctuations' . Figs. 6(C,D) illustrate this manipulation applied to the network simulation investigated in Fig. 6(A). PSTH of induced phase at the streaming frequency, Fig. 6(C), is highly similar across stimulus B presentations of different levels of deviancy, but less so for stimulus A. Phase for stimulus B appears to have twice the frequency of the sensory stream indicating higher mode locking for this particular network and stimulus. The phase PSTH for the two stimuli peak at different time points after stimulus offset, a common experimental observation. The time frequency plot of ITC, Fig. 6(D), averaged across all stimulus presentations (except the D1 and D12 stimuli, see Methods) for stimulus B shows a strong burst of ITC locked to stimulus onset around the theta to alpha frequency band (3-15 Hz), also often observed experimentally 19,31,33,35,[107][108][109][110] .
In Fig. 9(A) we investigate how ITC in the induced neuroelectric activity at the streaming frequency depends on the autonomous Lyapunov exponent λ U . In common with the reduced activity, Supplemental Fig. S3(A), this quantity peaks around the point where chaos is just stabilized and is weakened in irregular sensory streams. In contrast to the reduced activity it does not strongly depend on the level of deviancy. The fact that ITC in the induced activity does not strongly depend on deviancy does not mean that stimuli of different deviancy levels necessarily have similar phase PSTH profiles. Since phase PSTH can depend on the stimulus, as shown in Fig. 6(C), we expect it may take several stimulus presentations for the shift in phase locking to occur after the stimulus is switched in the roving sensory stream. To investigate this we calculate the distance between the phase distribution for deviancy level D and that for the standard D11 stimulus. This quantity is shown in Fig. 9(B) versus deviancy level D averaged across both stimuli and all network simulations within certain prescribed ranges of autonomous Lyapunov exponent λ U . It shows a gradual decay only in the suppressed chaos regime 0 < λ U < 2, confirming the slow adaptation to a new phase locked state in the induced activity after a stimulus switch. This behaviour is also found in the reduced activity, Supplemental Fig. S3(B). Time frequency ITC Fig. 9(C-F) averaged across both stimuli and all network simulations within certain prescribed ranges of autonomous Lyapunov exponent λ U show a strong burst of theta to alpha coherence at stimulus onset, especially in the suppressed chaos regime, Fig. 9(D).
Several studies show that low frequency brain oscillations can become phase entrained by sensory streams even when the streams are composed of multiple different stimuli 9,30,[99][100][101]104 . This is possible in the current model because even streams composed of different stimuli can stabilize chaos by attractor switching as described above. To demonstrate this we apply the model to a simplified version of the task recently investigated by 9 . Participants discriminated the orientation of target visual stimuli embedded within streams of 'random' visual stimuli termed 'distractors' (see Supplemental Fig. S6(A)). The authors found phase entrainment was strong at the streaming frequency and increased with the temporal regularity of the sensory stream.
We investigate phase coherence of the 'induced' neuroelectric activity. To generate this the stimulus conditioned PSTH is subtracted from the neuroelectric activity for each different stimulus type separately, including targets and distractors, in the same way as described above (see Methods). We find phase coherence at the streaming frequency peaks (Supplemental Fig. S6(B)) close to the point where chaos is just stabilized at this ITI of 400 ms (compare Fig. 3(A)) when distractors are identical. A strong peak is also shown when distractors are random but similar to each other (clustered.) The peak is weakened by temporal irregularity. However the peak is absent when stimuli are totally random. Indeed while totally random stimuli streams can stabilize chaos (so that the trajectory of the driven network activity is independent of the initial network state at the onset of streaming) they do not phase synchronize the activity. The distractor stimuli used in experimental studies like 9 should be considered 'clustered' since although different they have many features in common. Deviants generate larger instantaneous population responses than standards. We often find in inhibition dominated networks that stimulus presentations generate large negative instantaneous deflections in the neuroelectric activity. This can seen from the PSTH in the network simulation with both excitatory and inhibitory units, Fig. 7(B), as well as in the purely inhibitory network, Fig. 1(C). Here we investigate whether the magnitude of these instantaneous responses adapts with repetition number as often found in studies 18,28,47 . This effect is complicated by the possible presence of mode-locking. In networks where this occurs responses to odd numbered deviants (D1, D3, D5, …) may be larger than responses to even numbered deviants (D2, D4,…) for example. We also often find that D1 deviants generate instantaneous population responses of larger magnitude than standards for one stimulus type but not the other stimulus type in two stimuli roving sensory streams. The opposing character of the differential responses for the two stimulus types in the purely inhibitory network, Fig. 1(C), is apparent in the 50 ms stimulus presentation epoch in Fig. 5(A).
We measure the instantaneous stimulus response size, R Z D to stimulus type Z = A,B of deviancy level D as the difference between the mean neuroelectric activity during the 50 ms stimulus presentation epoch and the 50 ms immediately preceding the stimulus presentation (see Methods). The differential response size for D stimuli is then defined as where again |..| denotes absolute value and S = D12 is the standard stimulus. The stimulus averaged differential response, , is shown in Fig. 8(C, solid lines) versus λ U for inhibitory network simulations, while the opposing nature of the differential response sizes is illustrated by the quantity √ ± (dR dR ) Fig. 8(C, dashed line). Here √(x) ± denotes sign (x)√(|x|). Throughout the suppressed chaos regime this latter quantity is large and negative for D1 stimuli. This effect is much weaker for later presentations after D1 (not shown). Although D1 differential responses sizes are often opposing in the two stimulus types, deviants generate larger magnitude responses in the neuroelectric activity than standards on average over both stimuli throughout the suppressed chaos regime, an effect which relaxes slowly over about four stimulus presentations, Fig. 8(C, solid lines). In fact both these effects can also be found in a simple phenomenological model where activity does not fully relax to an equilibrium state between stimulus presentations (see Supplemental). To show these effects the phenomenological model relaxation timescale parameter must be of the order of the time between stimulus presentations. In contrast here long timescales are created when network chaos is (partially) stabilized while network model timescale parameters do not exceed 50 ms and remain biologically reasonable. Thus the experimental observations of large excess instantaneous deflections in LFP 18,28,47 for deviant stimuli compared to standards are qualitatively reproduced (see Discussion) only in the suppressed chaos regime.

Single units also show adaptation.
Stimulus response adaptation is also know to occur in multi-unit 18,28 and single unit 17,111 activity in both excitatory and inhibitory cells 20,21 . Figure 7(C-F) shows examples of single unit PSTH for standards and deviants from the network simulation including excitatory and inhibitory units under the n = 12 roving sensory stream. As can be seen the dependence on repetition number can be quite complex, both in the initial stimulus response and in the later modulations during the ITI. Stimulus presentations can generate phasic facilitating or suppressing responses to one or other stimulus type or both, Fig. 7(C,D). Stimulus induced modulations of tonically active units are also seen, Fig. 7(E,F). Stimulus induced modulations occur in both excitatory Fig. 7(C) and inhibitory units Fig. 7(D-F) and in units not directly driven by the sensory stream Fig. 7(C-E) as well as those that are, Fig. 7(F). Notice that this network simulation is almost completely stablized and the D11 and D12 responses are very similar.
These phasic responses result from the stimulus presentation induced 'kicks' described above. When the stimuli are presented the trajectory transiently moves away from the chaotic attractor region and then returns. During this period some single units silent during ITI periods are briefly activated. These units can suppress the ones tonically active during the ITI period resulting in negative deflections like those shown in Fig. 7(E,F) which may also appear in the population neuroelectric activity, Fig. 7(B). According to the above results we expect deviants to produce larger single unit responses than standards for one stimulus type but not the other, while deviant responses will be on average larger. On the other hand due to the presence of phasic activations in some units and suppression of tonic activity in others this may not show up if firing rates during stimulus presentation are simply averaged across all units.
To confirm these observations we calculate stimulus specific adaptation indices (SI) for single units, as in 17,111 , (see Methods.) SI vary from 1 to −1, depending on whether a unit is more active during the D1 or S presentation of a stimulus. The black line in Fig. 8(D) shows SI averaged across all the active units (see Methods) for many inhibitory network simulations under the n = 4 roving sensory stream versus λ U . Both stimulus types, A and B, are combined. When calculated in this way SI are close to zero throughout the range of λ U and do not show a strong dependence on it. However when the SI average is restricted to include only units which show positive phasic responses to either the D1 or S stimulus presentations (see Methods), or both, Fig. 8(D, red), SI is significantly positive but only in the suppressed chaos regime. When the calculation includes only active units which do not show positive phasic responses to either D1 or S stimulus presentations, Fig. 8(D, green), the opposite is true. Thus phasic units are more strongly activated on average by deviants than standards while tonic units are on average more strongly suppressed by deviants than standards.

Discussion
We investigated how stimulus streams interact with spontaneously generated dynamics in neural networks without any explicit learning mechanisms. Sensory stream adaptation may be a general property common to many different brain areas 13,22 embedded in multiple possible neural substrates on different spatial and temporal scales during various behavioural task settings. To illustrate how a simple mechanism can unify a range of observations we have necessarily sacrificed model specificity for model generality. The model is the simplest one may imagine; a random network of identical units, each one described by a single dynamical variable with a single timescale driven by completely random stimuli. Despite this its emergent behaviour is so complex that it should be applied to each experimental setting individually to understand the related results. It can therefore provide deep insight into the origin of many seemingly contradictory experimental results. It will be particularly relevant in music, language and birdsong however repetitive sensory patterns may be generated by the motion of the animal itself. In particular chaotic network oscillations may occur as part of central pattern generators.
Remarkably we showed that stimulus sequences produce realistic repetition effects on neural response even without synaptic modification. Deviant stimuli generate long lived differential oscillatory modulations such as mismatch negativity and repetition positivity in population activity 19,[23][24][25][26][27] . They are automatic, independent of the behavioural task and of attention 27 . Quantitative comparisons in a model of this generality are not possible since the form of the modulations depends on the task and brain area under study and also varies across individuals. We found good qualitative agreement with studies (e.g 19,25 .) when chaos is (partially) stabilized by stimulus streams. For greater agreement with experiment in future work the delay between the time of the physical stimulus and the time the signal reaches cortex, not included in the present model, should be taken into account, as part of more detailed modeling of the pre-cortical sensory processing. We showed partial stabilization allowed stimuli to generate distinct and reliable deviancy-level dependent responses in periodic roving streams.
We also found that deviants generate strong differential responses in temporally irregular streams 30,112 but the response magnitude was weaker, in striking agreement with 30 . We showed this was because PSTH variability reflects stream irregularity in (partially) stabilized networks. Actually the effects of different possible forms of sensory stream variability on this model are complex. In the future this model will be profitably applied case by case to understand neural response in sensory streaming experiments with interacting factors. For example when stimulus contrast and stimulus novelty interact (e.g 38 .), or where temporal and type expectations interact [99][100][101] .
We also investigated the stimulus response specifically. In inhibitory networks D1 deviants on average generated larger population responses in neuroelectric activity than standards in one stimulus type. This tendency was only weakly shown for later repetitions. The average neuroelectric response size over both stimulus types was also larger for deviants than standards as is often observed in LFP studies 18,28,47 . Single units showed complex repetition dependent modulations both in their initial stimulus induced phasic activations and for longer periods during the ITI which are intriguingly similar to those found in recent studies 20,21,47 and thought to result from recurrent network inhibition. On average single units showed stronger phasic responses than standards in striking agreement with studies at the multi-unit 18,28 and single unit 17,111 level.
In roving sensory streams we found adaptation to the new single stimulus stream after a deviant can happen very slowly over multiple trials, in agreement with multiple studies (e.g 25 .). In contrast to models based on synaptic mechanisms [113][114][115] where adaptation timescales reflect synaptic time constants here slow adaptation occurs because stabilization of chaos brings the driven dynamics close to critical λ D ≈ 0. The resulting adaptation timescales are around 20 times longer than the 50 ms synaptic time constant. There are several observations which provide convincing support for this which can't be easily reproduced with purely synaptic mechanisms. First it accounts for the dependence of adaptation on temporal regularity of stimulus presentation 30 . Second it accounts for the generation of multiple adaptation timescales and adaptation to whole stimulus sequences 29,111 . Third it accounts for the fact that deviant stimuli may produce stronger responses than unadapted frequency matched controls 19,28,40,41 . However we do not expect network adaptation to operate in isolation from synaptic effects. Rather we hope recognition of the important fact that networks can adapt in a biologically plausible way will provide exciting insight into recent studies of how recurrent local inhibition and synaptic modifications interact 20,21,42,43,47 Our results are in good agreement with experiments that demonstrate that low frequency brain oscillations become phase locked to stimuli presentations in streaming perceptual discrimination tasks [9][10][11]14,15,106,[116][117][118] . Phase locking occurs across a broad range of cortical areas 14 , while locking outside the streaming frequency band is mostly limited to primary sensory areas. Stimuli presentations are associated with increased ITC and power predominantly in the theta to alpha band 31,33,107,108 , while deviants also generate differential responses in roving paradigms 19,35,109,110 , in qualitative agreement with this model. Ongoing brain oscillations are often thought to be emergent properties of network mediated population dynamics [54][55][56]119,120 close to a phase transition 121, 122 . Several experimental studies have convincingly demonstrated they may be chaotic 123,124 . Interestingly low dimensional chaos may be stronger during slow-wave sleep than REM or wakefulness 58,125,126 , consistent with our hypothesis of chaos suppression by sensory driving. Intriguingly several studies have also found direct evidence for unstable periodic orbits in neural dynamics [127][128][129] , which are related to phase entrainment of chaos, while other studies directly demonstrate phase synchronization 130 .
Recurrent networks provide a natural substrate for various forms of neural computation and sequence learning [48][49][50][51][52]131,132 . SORN networks 133 utilize plasticity to modify network structure and are therefore based on completely different principles to the current model. Like other studies 64,65,89,[134][135][136] in various task settings we find optimal performance in a critical state or at the edge of chaos. However we also show that consistency with experiment does not require this. We found stream history dependence in a wide regime where chaos is partially suppressed. This is surprising given the inherent instability of chaos 79 but other studies find that significantly chaotic dynamics is consistent with neural information processing 77 . In the partially stabilized regime driven trajectories remain 'bundled' together so that the separation between them tends to a constant finite level. This is characteristic of chaotic synchronization 75,76,85 and phase synchronization 63 . Phase stabilized chaotic trajectories combine divergence (high signal) with reliability (low internal noise) and therefore facilitate stream adapted response. Furthermore we found that the amount of stabilization, while strongly dependent on stimulus presentation frequency 95 , was only weakly dependent on the temporal or type predictability of stimuli in the stream. This allowed stimulus response reliabilites to be strongly determined by the variability of the preceding sensory stream. Finally we demonstrated this property emerged because stimuli move the network state away from its low dimensional chaotic attractor region. Thus we expect the results presented here to generally apply to chaotic networks in parameter regions where they can be stabilized, irrespective of some details of network physiology and structure.

Methods
Network simulations. We use a rate network model, where f i (t) are the firing rates for the N firing rate units i. I i (t) is the input current to unit i given by, This equation models the variations of synaptic neurotransmitters which control cell firing rates. The synaptic timescale is set as τ g = 50 ms. The term γ − θ + I t I ( ( ) ) , ( − θ I t I ( ) for I(t) − I θ > 0 and zero otherwise) is the dependence of firing rate on input current close to saddle-node on invariant circle bifurcation to firing from quiescence as occurs in 'Type I' cell models. γ = 0.09 is a fixed parameter. V Z,K = 60.0,5.0 mV are fixed parameters. I θ = 4.51 μ A/cm 2 is the current at firing threshold. These firing rate units may be used to represent single cells or alternatively local populations of cells, as in 69 for example. k ij is the fixed network connection matrix. We investigate two network models (i) a 500 unit purely inhibitory network and (ii) a 1200 unit network with 500 inhibitory units and 700 excitatory units. Most simulations for statistical calculations are performed on the 500 unit network. However we also demonstrate that the larger network has the capacity to show similar behaviour (see Supplemental Methods for more details.) g i Z is the driving excitation for unit i during sensory stimulus Z. Simulations generally include two sensory stimuli, denoted A and B. These stimuli are always applied to the network for 50 ms time periods separated by longer inter-trial-intervals (ITIs). During ITIs a different excitatory input, denoted X, is applied in place of the sensory stimuli. The input g i Z is drawn randomly from the same given distribution for each unit i and each stimulus Z in a simulation, whether it be A, B or X, and fixed for the duration of the simulation. The parameters are set so that the average input current to each unit i is just above firing threshold I θ (see Supplemental Methods for more details.) All units i in the network are continuously driven above firing threshold by the external driving excitation by all stimuli Z. Only network inhibition can cause them to become quiescent for periods.
For generality we also investigate an excitatory-inhibitory network simulation where only half of the excitatory units and only half of the inhibitory units are directly driven by the sensory stimuli while the other half receive constant unvarying excitatory driving input. In this case the sensory stimuli A, B and background input X are determined as above but then g i A and g i B for half of the units i are set to g i X . As a proxy for the LFP/EEG we use the absolute synaptic input current per unit where |x| denotes the absolute value of x 57,70 . In the simulation where only half of the units are directly driven by the sensory stream we divide this into four separate contributions, originating from the excitatory and inhibitory units directly and not directly driven by external sensory input.
In Eq. 1 η i (t) are iid Gaussian distributed noise terms with 〈η i (t)〉 = 0 and 〈η i (t 1 )η j (t 2 )〉 = δ ij δ(t 1 − t 2 ). The parameter α controls the strength of this stochastic term and is always set to α = 0.1 so that fluctuations have the expected size 10% of the rates f i , except in Lyapunov exponent calculations and the time series examples where α is set to zero. The noise represents fluctuations arising from multiple sources including incoherent spike arrival times which have been averaged out in the present rate dynamics network model. Noise has no dynamical role in the results presented here; it is included purely to demonstrate their robustness.
All simulations of the noisy rate network were performed with a stochastic weak second order Runge-Kutta integrator 71   The mean absolute differential dPSTH shown in Fig. 8(A,B) for each network simulation, stimulus Z and repetition D are given by where |x| denotes the absolute value of x. dPSTH are calculated from standardized mean activity time series under the n = 12 roving sensory stream. The sum runs over all time lags τ from stimulus onset to the onset of the following stimulus. Figure 8(C) shows results based on instantaneous neuroelectric activity responses R Z,D under the n = 12 roving sensory stream.
where the sum over τ runs over only the 50 ms stimulus presentation epoch and −τ runs over the 50 ms time epoch immediately preceding stimulus presentation in 5 ms increments. Then differential instantaneous response is given by dR Z = |R Z,1 | − |R Z,12 |. Figure 8 To calculate time-frequency plots a Hanning window of length 640 ms is moved along the time series in 5 ms steps, t. From the complex coefficients r(t, f) (e iθ (t, f)), log power, 2log(r(t, f)), and complex phase e iθ (t, f) are computed and averaged to create stimulus onset locked power PSTH and inter-trial coherence (ITC). The absolute value of the mean complex phase is taken for the ITC which varies between zero and unity. PSTH are created in the same way as described above, for both stimuli A and B separately and all deviancy levels D separately. The time-frequency plot shown in Fig. 6(D) is the ITC averaged across stimulus B presentations of all deviancy levels except the first and last. These are omitted to avoid any end-effects originating from the stimulus switching in the 640 ms time window. The time-frequency ITC plots shown in Fig. 9(C-F) are created the same way except both stimuli and all network simulations within given ranges of λ U are averaged. To create excess ITC and excess log power PSTH plots in Supplemental Fig. S4 and Supplemental Fig. S5 time-frequency spectrograms calculated from D11 stimuli of given type are subtracted from D1 stimuli of the same type. Then both stimuli and all network simulations within given ranges of λ U are averaged. Again the D12 stimulus is not used to avoid possible effects originating from the stimulus switch.
To investigate phase coherence specifically at the streaming frequency in Figs 6(C) and 9(A,B), and Supplemental Fig. S3 a more frequency specific method is used. First the reduced or induced neuroelectric is divided into sections. Each section starts 480 ms after the offset of a D12 presentation of a given stimulus, say A, and ends 520 ms after the offset of the subsequent D12 presentation of the other stimulus type, say B. Thus each section includes presentations of only one stimulus type. Since the sampling period is 5 ms each section includes exactly 2048 grid points. In the irregular sensory stream 2048 grid point sections starting 480 ms after the offset of a D12 presentation are used. If a 2048 grid point section extends past the point of first presentation of the other stimulus type it is zero-padded from that point. Each section is band passed in a narrow band around the streaming frequency, (1/0.85)±0.3 Hz and the phase time series, θ(t), extracted using Hilbert transform. To calculate stimulus locked phase coherence, the distribution of phase, p(θ|τ, Z, D), τ ms after the onset of stimulus type Z of deviancy level D is obtained from multiple sections. For phase coherence we use the modulation index 137 which is defined as 1 minus the normalized entropy of this phase distribution, where N is the number of bins i. This quantity is averaged across the period 350 < τ < 650 ms and over both stimuli Z to obtain the phase coherence for each deviancy D. For the irregular stream we only include ITIs longer than 600 ms in this calculation. To calculate the difference between phase coherence for different deviancy levels in Fig. 9(B), Supplemental Fig. S3 where P is the cumulative distribution, ∑ runs over the whole distribution and |..| denotes absolute value. This is averaged across 350 < τ < 650 ms and over both stimuli Z to obtain KS(D). Phase coherence at the streaming frequency in Supplemental Fig. S6 is calculated in a similar way. First 'induced' activity is obtained from the neuroelectric activity by subtracting the stimulus conditioned PSTH H Z (τ), where stimulus Z refers to the targets A and B and each of the six distractors, a total of eight stimuli. Next sections of length 16384 grid points (with time period of 16384 × 5 ms) are extracted, band passed in a narrow band around the streaming frequency (1/0.45) ± 0.4 Hz and the phase time series generated by Hilbert transform. Next the distribution of phase, p(θ|τ, Z), τ ms after the onset of target stimulus type Z = A,B is obtained from multiple target stimulus presentations. Finally the modulation index is calculated from the distribution and averaged across the period 300 < τ < 440 ms after stimulus onset and over both targets A,B. Irregular sensory streams are treated exactly the same way because each target is always followed by a 400 ms ITI.
Single unit phasic analysis. Each unit was defined as active for a particular epoch if its mean firing rate exceeded 10 −5 Hz during the epoch, for example during the final 400 ms of the ITI period following a D1 stimulus of type A. This minimum firing rate could be varied without change in results however some minimum firing rate is necessary for the calculation of stimulus specific adaptation indices SI 17 . These are defined by = − + SI SER Z D SER Z S SER Z D SER Z S ( ( , 1) ( , ))/( ( , 1) ( , )) i Z i i i i for a unit i and stimulus type Z = A,B where 'stimulus epoch rates' SER i (Z,D1) and SER i (Z,S) are the firing rates of unit i during the 50 ms stimulus presentation epoch for deviant D1 and standard S presentations of stimulus type Z. For a unit to have an SI i Z it must be active in either stimulus presentation, i.e. SER i (Z, D1) > 10 −5 or SER i (Z, S) > 10 −5 .
The late epoch rates, LER i (Z, D1) and LER i (Z, S) for cell i and stimulus type Z are defined as the mean firing rate during the final 400 ms of the ITI period following the D1 or S presentations of stimulus type Z. The associated stimulus responsivities, SR i (Z, D1) and SR i (Z, S), are defined by SER i (Z, D1)/LER i (Z, D1) and SER i (Z, S)/LER i (Z, S). Supplemental Fig. S7 shows SR i (Z, D1) versus LER i (Z, D1) and SR i (Z, S) versus LER i (Z, S) for all 500 units in both D1 and S presentations of both stimulus types Z, (a total of 2000 points), for the partially stabilized network simulation shown in Fig. 8(A) under the n = 4 sensory stream. Single units are divided mainly into two distinct clusters. The cluster with SR close to unity and LER around 20 Hz contains the units tonically (or intermittently) activated by backgound input X during ITI periods. The cluster with LER around 10 −3 Hz and SR around 10 −3 contains the units phasically activated by D1 presentations (black plusses) and S presentations (red crosses).
A unit i is denoted 'phasically responsive' for a particular stimulus type Z and presentation type, D1 or S, if the associated SR i > 10 and LER i > 10 −5 . These phasically responsive units are the ones in the upper right quadrant demarcated by the pink dashed lines in Supplemental Fig. S7.