Resonance with subthreshold oscillatory drive organizes activity and optimizes learning in neural networks.

Network oscillations across and within brain areas are critical for learning and performance of memory tasks. While a large amount of work has focused on the generation of neural oscillations, their effect on neuronal populations' spiking activity and information encoding is less known. Here, we use computational modeling to demonstrate that a shift in resonance responses can interact with oscillating input to ensure that networks of neurons properly encode new information represented in external inputs to the weights of recurrent synaptic connections. Using a neuronal network model, we find that due to an input current-dependent shift in their resonance response, individual neurons in a network will arrange their phases of firing to represent varying strengths of their respective inputs. As networks encode information, neurons fire more synchronously, and this effect limits the extent to which further "learning" (in the form of changes in synaptic strength) can occur. We also demonstrate that sequential patterns of neuronal firing can be accurately stored in the network; these sequences are later reproduced without external input (in the context of subthreshold oscillations) in both the forward and reverse directions (as has been observed following learning in vivo). To test whether a similar mechanism could act in vivo, we show that periodic stimulation of hippocampal neurons coordinates network activity and functional connectivity in a frequency-dependent manner. We conclude that resonance with subthreshold oscillations provides a plausible network-level mechanism to accurately encode and retrieve information without overstrengthening connections between neurons.

Oscillations in local field potential (LFP) largely reflect coherent postsynaptic potentials among neurons (1). These rhythms are behaviorally relevant, and their features are highly predictive of cognitive processes, and plasticity, in underlying neural networks (1–5). Network oscillations have long been thought to promote plasticity by precisely timing firing between pairs of neurons [i.e., driving spike–timing-dependent plasticity (STDP)] (6, 7). However, it is still unclear whether specific network rhythms are critical for specific neural computations and, if so, why this is the case.

Neurons display complex behavior in response to oscillatory input. Many neuronal subtypes show enhanced membrane voltage responses to periodic subthreshold inputs within narrow frequency bands (8–10). Critically, the frequency at which neurons resonate can shift in response to depolarizing or hyperpolarizing inputs (11–13). Thus, in addition to simply integrating inputs to generate an action potential, neurons are biophysically suited to perform time-dependent computations, including input filtration, based on their periodicity.

The theta (4 to 10 Hz) rhythm is a prominent oscillation present in mammalian brain networks (3). Within the hippocampus, theta plays a central role in the function of place cells, which encode spatial and contextual information (14, 15). Place cells show several interesting features associated with theta-resonant firing. First, their firing phase varies over time, relative to hippocampal theta—a phenomenon called phase precession, which occurs as animals move through their environment (6, 16, 17). Second, sequences of place cell activation occurring during spatial exploration are replayed during subsequent theta oscillations, and surprisingly, these replay events can occur in either the forward or reverse direction (18–21). While the idea that theta (and other hippocampal oscillations) plays a role in hippocampal function is widely accepted, the underlying mechanisms for phase precession, forward replay, and reverse replay—and the link of these features to memory formation—are still largely unknown.

Networks of neurons that display resonance shifts (i.e., the firing response to subthreshold oscillating input changes as a neuron is depolarized) show enhanced pattern formation and separation when rhythmic inputs are present (22, 23). Here we show that resonating networks have a firing pattern that is highly beneficial for both encoding and retrieving patterns of external inputs. Using conductance-based model neurons that display resonance with subthreshold oscillatory input, we show that networks will organize the firing of neurons around an oscillation in a manner that represents an external input. When synapses are able to evolve via a STDP rule, an input will be reliably encoded within the synaptic weights of a network. This leads to the subsequent reproduction of the input-induced firing pattern in the absence of the external pattern, for both static and temporally dynamic inputs. We also show that resonance with subthreshold oscillations provides a network-level mechanism both for theta phase precession and for forward and reverse replay, which reliably happens across any resonant frequency. Finally, we find that subthreshold periodic input induces stable, highly organized functional connectivity over the theta band, in both simulated and in vivo networks. This work demonstrates that resonance with subthreshold oscillations organizes neuronal firing phase with respect to network rhythms and thereby facilitates the encoding and retrieval of information.

Results

We investigated how resonance with subthreshold oscillations affects pattern and sequence learning, using modeled networks of neurons that receive three types of input (Fig. 1). First, each neuron in the network receives a unique level of external, direct current (DC) indicated by the color map. Second, the entire network receives uniform oscillating input (with modifiable frequency and magnitude). Third, individual neurons receive the summed presynaptic input from other neurons in the network. The weights of individual synapses evolve via STDP across the learning phase of simulations.

Fig. 1.

Input-dependent resonance shift allows for selectively activating subsets of neurons. (A) Model neurons receive three types of input. External input is DC, which varies in magnitude with neuron identity, represented by the color mapped arrow. All neurons receive an identical oscillating input, represented by the sine wave. Additionally neurons receive the synaptic inputs from neighboring neurons according to network connectivity and synaptic weights. (B) The input-dependent resonance shift manifests as a broadening of the resonance curve with increasing excitation of the neurons. (C) Broadening of the resonance curve also occurs for changes in synaptic weights, which provides for selective activation of subsets of neurons based on synaptic coupling. Dashed lines show the frequencies corresponding to the raster plots in D–F, which show the divergent activation for frequencies between 12 and 16 Hz. Error bars, SEM.

Input-Dependent Resonance Shift Allows for Selective Activation of Subsets of Neurons.

The neurons in the model display input-dependent resonance shifts (Fig. 1). A neuron will respond to a wider range of oscillation frequencies if it receives a larger DC input. There are two main regimes apparent in the resonance profile: (i) a 1:1 regime where the neuron fires one spike per cycle at low frequencies and (ii) a 1:2 regime where the neuron fires every other cycle at high input frequencies. For an oscillation of 0 Hz (i.e., in the absence of any oscillation), an additional DC current is added to the DC input so that neurons receive the same total input magnitude as when an oscillation is present. This case does not lead to neuronal spiking.

The broadening of the resonance response occurs within networks as well (Fig. 1). To show this, we formed three clusters within a network with varying intracluster coupling (0.2, 1.0, and 1.4 ), while keeping intercluster coupling constant. This leads to groups with high (green), moderate (light blue), and low (dark blue) synaptic input. The raster plots in Fig. 1 show network activity at 12, 14, and 16 Hz and demonstrate how increasing the frequency of the oscillation provides for selective activation of clusters with stronger coupling.

Networks Learn Patterns of External Input and Reproduce the Reverse.

To investigate the basis of learning through synaptic plasticity in this model, we had networks encode a pattern of external input (a set of DC inputs with varied magnitude across the network) to connections (Fig. 2). We monitored the phase at which neurons fired relative to the oscillations, as a function of their input magnitude. The simulations were split into five phases: before the input pattern (red in Fig. 2), during patterned input (yellow), after pattern learning has saturated (green), and two subsequent replay periods (replay periods 1 and 2; with and without prior patterned DC input). During the period before the input pattern and the replay periods, all neurons received the same moderate DC input and STDP was disabled. The first replay period shows the effect of learning the input pattern, and the second shows the effect of playing the stored pattern back (i.e., no input pattern is present) with active STDP.

Fig. 2.

Resonating networks learn by mapping input patterns to synaptic weights. (A) Raster plots show the relationship between the phase of firing and the external input to the neuron. Black lines show the trace of the oscillating input, and the color of the rasters shows the DC input to the given neuron. Neurons are sorted by their input rank. Subpanels in A correspond to before DC input distribution is applied (pre-), with DC input distribution (input), after learning has saturated (post-), after learning/no DC distribution (replay 1), and after a second period of learning with no DC distribution (replay 2). (B) The relationship between firing phase and DC input varies between negatively, positively, and not correlated for different epochs of the simulation. Data are averaged over 10 cycles of the oscillation. Error bars = SEM. (C) Transitioning from the input pattern-depending firing phases to synchronous firing is gradual. Lines trace the firing phase of 12 neurons with varying input magnitudes across time. The horizontal bars above indicate when the external input and learning are present (white, input but no learning; black, learning and input; gray, no input and no learning).

The raster plots in Fig. 2 show the evolution of firing phase across each period of the simulation. The color indicates the magnitude of input current a neuron receives, and neurons are sorted by this value, with highly activated neurons having a higher input rank. Before any input, neurons fire randomly over a narrow band of phases (Fig. 2, Far Left). The input pattern leads to organized firing with highly activated neurons firing at earlier phases (Fig. 2, Inner Left),with the neurons receiving larger current firing earlier on the oscillatory cycle and neurons that receive smaller DC input following, with the range of firing phases being determined by the spread of activating input (Fig. S1). This variable phase locking is a well-known phenomenon observed during synchronization of weakly forced oscillators where there is a small detuning of mutual frequencies of the drive and the oscillator (for example, see ref. 24). The neurons in resonance behave as oscillators, and their specific frequency depends on the properties (height and width) of their resonance curve, the shape of which is in turn contingent on the magnitude of DC input (Fig. 1).

As the pattern is learned, the overall phase shifts, but neurons return to firing at a uniform phase, independent of their DC input (Fig. 2, Center). This convergence is due to the universal learning rule, which mimics STDP (25), where the synapse is being strengthened (or weakened) when the presynaptic neuron fires within a narrow window before (or after) the postsynaptic neuron. As long as the neuronal pair fires in an ordered sequence, the corresponding synapse gets systematically potentiated or (weakened), leading to increased synaptic input to the neuron having lower DC input. When synaptic input offsets the difference in DC input between the two neurons, the neurons fire simultaneously, resulting in the termination of synaptic potentiation (depression). For this process to be effective, the time length of the EPSP has to be on the order of , where is the oscillation frequency. For theta frequencies, this constitute a time constant of 100 to 300 ms, roughly corresponding to an activation time constant of NMDA receptors (26). However, if the reactivation happens at higher resonant frequencies, as shown in the next section, this activation time constant can be significantly smaller.

When learning is suspended and the external input pattern is removed and all neurons receive the same intermediate DC input, the network shows the reverse pattern of activation (Fig. 2, Inner Right), as now the relative patterns of cellular input are dominated by synaptic currents. After a second period of learning (but with a uniform external input), the network returns to firing at a uniform phase, effectively erasing the stored pattern (Fig. 2, Far Right). The above relationships are summarized in Fig. 2 as we plot relative phase of neuronal spiking as a function of their DC input magnitude for each phase described above (red, before input pattern; gold, input pattern; green, after learning saturates; blue, replay of stored pattern; violet, replay after erasure). Fig. 2 depicts the time course of the evolution of firing phase for 11 neurons having different DC input values. The bars below indicate the timeline when input and learning are present (white, input but no learning; black, learning and input; gray, no input and no learning).

The precise firing phase versus input relationship is dependent on total input to neurons being subthreshold; superthreshold input disrupts this relationship and impedes subsequent learning (Fig. S2). On the other hand, the sign of the current in oscillatory drive does not affect the observed results. Namely, if an oscillation is purely hyperpolarizing, the same pattern of phase organization is observed (Fig. S3). The critical components to this learning and replay mechanism are resonance at the single neuron level and the presence of a subthreshold oscillation (Table S1). The LFP is a complex oscillation with a waveform that superimposes multiple frequencies. For example, sharp-wave ripples are composed of a high-frequency ripple riding on top of a lower frequency sharp wave (27). We tested the robustness of this input learning mechanism to a complicated waveform combining 6 Hz and 120 Hz oscillations (Fig. S4). The input versus phase relationship and pattern reversal after learning were both reproduced with this waveform.

Stored Patterns Can Be Replayed for Any Resonant Frequency.

To demonstrate the generality of the pattern storage and replay mechanism, we introduce a second conductance-based neuronal model based on classic Hodgkin–Huxley (HH) dynamics (28). This model neuron displays spiking resonance in response to subthreshold oscillating input in the gamma band between 40 and 90 Hz (Fig. S5), which is well above the resonance band of the previous model. For ease, the neuronal models will be referred to as Ks for the neuron that resonates in the theta band and HH for the gamma-resonating neuron. In Fig. 3 we show that patterns stored during resonance at one frequency (theta band in our case) will be replayed at a higher frequency (gamma band) in a similar reverse firing order. Here, synaptic weights (and corresponding adjacency matrix) taken from a network of Ks neurons after a pattern was learned are used to connect a network of HH neurons (after adjustment for differences in excitability; see Fig. S6). HH networks replay a pattern on a similar phase range as a Ks network, though with higher variability for late firing neurons. This result further indicates that the described resonance mechanism is very robust to frequency modulation that may occur, for example, during cycles between active behavior and rest.

Fig. 3.

Replay of stored pattern occurs independent of neuronal model and frequency band. (A) Input-induced pattern of firing phase for a network of Ks neurons driven with a 6 Hz oscillation. (B) Reversal of pattern during replay after learning for a network of Ks neurons driven by a 6 Hz oscillation. (C) Reversed pattern replayed by a network of HH neurons at 60 Hz. All raster plots include spike from 10 cycles of the oscillation, and the color of a neuron’s raster indicates the magnitude of DC input it gets in a pattern. (D) Firing phase versus DC input relations for the three above cases (red Ks neuron before learning, green Ks neuron replay, blue HH neuron replay). Error bars = SEM.

Pattern Learning Saturates Naturally in Resonating Networks.

The results described above indicate that neuronal firing phases rapidly converge during learning and that this process minimizes the firing phase difference between neurons. This behavior should result in two interesting phenomena: (i) synaptic strengths will stop changing when the phases converge, and (ii) input differences between neurons will map onto their synaptic weights. To test these effects, we presented an input DC pattern to network for a long time period and tracked the time course of synaptic change. If the learning rate (the magnitude of synaptic change corresponding to ) allows, both the maximum (Fig. 4) and mean (Fig. 4) synaptic weight will saturate before the end of the simulation. Regardless of learning rate, there is a large increase in synaptic change followed by a gradual decline to no change in synapse strength (Fig. 4). The time of peak synaptic change is delayed for slower learning rates. Note that the input pattern is the same for all conditions in Fig. 4 ). Both the final mean synapse strength (Fig. 4 black) and time it takes to saturate (Fig. 4 red) depend on the range of currents in the external pattern. The time to saturation is the time it takes for the mean change in synaptic strength to fall permanently below .

Fig. 4.

Learning saturates naturally after input pattern is completely mapped to synapses. Saturation of learning reliably occurs given that the learning rate is high enough for the given time. Both maximum (A) and mean (B) synaptic weight saturate. Line color indicates network learning rate. (C) The majority of synaptic change occurs early during the learning period and then gradually decreases to zeros. (D) Final mean synapse strength and time until learning saturates depends on the spread of the input distribution. Error bars = SEM.

Saturation of learning occurs when the input pattern is fully mapped to the synaptic weights in the network, a phenomenon quantified in Fig. 5. The mapping of the input pattern is reversed in the synaptic weights. Highly activated neurons, which fire at an earlier phase, strengthen outward connections (black trace) while weakening inputs (red trace). Neurons given lower external inputs do the opposite, strengthening inputs and weakening outputs. This leads to the external input pattern and the synaptic input pattern being complimentary, leading to all neurons receiving the same net input.

Fig. 5.

Input pattern maps to both synaptic inputs and outputs. After learning, input strength (black) is anticorrelated with input magnitude of a neuron in the pattern, and output strength (red) is correlated. Error bars = SEM.

Overall, neurons with the lowest DC current within the input pattern strengthen inputs more than the rest of the network, while highly activated neurons do the opposite. The new pattern of synaptic connectivity is complementary to the input pattern, which leads to all neurons firing at the same phase. Synchronous firing terminates learning, because as spike–time differences between neurons approach zero, there is no net synaptic change (simplified in our model as zero synaptic change for ms). When the external input is removed, the complementary synaptic input distribution lead to a reversal in firing order from the input pattern (see Fig. 2 in previous section).

Resonance with Subthreshold Oscillations Facilitates Sequence Learning and Replay.

Next we investigated whether we can use the observed resonance shifts to store sequential neuronal activation to model the phenomenon of sequential replay following experience (14). Sequences were generated by delivering a slowly varying current to sequentially activate subsets of neurons (Fig. 6, solid lines), with each group resonating with the oscillating current in turn. This current is to model the preferential activation of subpopulations of place cells as an animal traverses a series of spatial locations. The asymmetry in its shape is to model the forward approach of the animal to a given location. It also provides temporal input relationships between neurons, to strengthen connections between neurons activated in a prior location and those activated in the current location. During the course of sequence presentation, groups of neurons display dynamic phase relationships (Fig. 6 ), where neurons that are highly activated fire earlier. For a single group, during the rising phase of activation, the firing phase will move earlier for each cycle of the oscillation (i.e., firing phase precession is observed). Between groups, those that are at peak activation will fire at earlier phases than less activated groups. These phenomena result from the relationship between activation and firing phase (Fig. 6) and the result of the input-dependent resonance shift (Fig. 1). The activation sequences were presented to the network 10 times, during which synapses were allowed to evolve using the same STDP-based learning rule as before.

Fig. 6.

Sequential activation of network subgroups leads to phase precession through time. (A) Networks were sequentially activated by a slowly varying depolarizing current delivered to subsets of neurons (solid lines; color indicates group). Spiking activity of each group is represented by the raster plots of different colors. Neurons are sorted on the y axis based on their lattice location, such that neurons closer on the y axis are more likely to be connected. (B and C) At the transition between the activation of two groups, phase order changes so that the neurons receiving the highest activation always fire at an earlier phase, leading to a phase precession through time. This is shown (B) in relation to the activation of the groups and (C) with reference to the oscillation. (D) The activation versus phase relationship shows that neurons fire earlier and with less variability with more depolarizing input. Error bars = SEM.

After this learning phase, the sequence can be reproduced in both the forward (Fig. 7) and reverse directions (Fig. 7). Both types of replay occur under different dynamical conditions. Reverse replay occurs when the whole network is depolarized to resonate with the oscillating input, but all neurons are activated to the same extent (i.e., each neuron receives the same DC input). This is due to the fact that neuronal groups in the end of the sequence receive larger overall input than groups activated at the beginning due to asymmetry in connection strengths. This results in an earlier phase of activation when the network resonates with the oscillatory current. In contrast, forward sequential replay occurs when the network is driven by external noise, in the absence of an oscillation. The neurons that fire early in the sequence subsequently depolarize neurons at the adjacent location, making them more prone to fire. Summary data are shown in Fig. 7 for reverse replay firing phase among the five groups. During reverse replay, groups activated earlier in the sequence reliably fire at a later phase of the oscillation (red trace). Without any learning (i.e., without STDP), groups generally fire at the same phase of the oscillation (black trace). During forward replay, the feed-forwardness of the intergroup connections dominate. The original firing order of the groups is reproduced, and early groups fire before late groups (Fig. 7, red trace).

Fig. 7.

Learned sequences can be replayed in both the forward and reverse directions. After a sequence is learned, (A) reverse replay occurs when a network is driven by an oscillation and (B) forward when the network is driven by noise. Raster plots follow the same organization as Fig. 6. (C and D) The firing relationships between groups is stable across cycles and different from groups without learning. (E) Synaptic connections between groups encode the sequence direction between groups, while weakening the reverse direction. (F) The directionality of intergroup connections emerges gradually after repeated sequence presentations, with increasing variations between groups. Error bars = SEM.

Sequential learning leads to connections being strengthened in the same direction of the sequence (feedforward) and weakens connections in the reverse (feedback). Mean synaptic weights between groups show strengthened connections in the direction of the sequence and weakened connections in reverse (Fig. 7). This is quantified for the entire network by the direction index, which is , where is the mean synaptic weight of connections between groups (Fig. 7). Critically, the feed-forwardness of the connections varies between groups and increases with every sequence presentation.

Functional Network Structure Emerges in the Theta Band.

We next sought to compare the behavior of simulated networks with experimentally observed pattern formation in the in vivo networks. Information representation and subsequent encoding using STDP-type learning rules require stable spike–time relationships. In our model, resonance with periodic input leads to stable spike–timing phase relationships. To quantify this effect, we measure functional network connectivity and stability of the observed functional relationships using metrics that were developed and validated in our laboratory (29) and compare it to results of the same analysis on experimental data.

In networks driven by oscillatory input (a 0.3 amplitude sine wave with a 0.3 DC offset) and background noise, oscillatory input leads to highly organized functional network structure between 4 and 10 Hz (Fig. 8). We quantified functional connectivity in three ways: spike–LFP coherence, mean average minimum difference (AMD) z score, and functional network stability (29). Spike–LFP coherence, which represents the reliability of the time of spikes within the LFP oscillation across the entire network, shows a noise-dependent resonance effect for stimulation between 3 and 13 Hz (Fig. 8). AMD z score and functional network stability are related measures that are based on the pairwise relationships between spike times of neurons across the network. The average significance (z score) of AMD measures between neurons shows a narrow resonance effect between 4 and 10 Hz, with a peak effect at 6 Hz, which depends on the level of background noise (Fig. 8). Functional network stability, which captures how similar AMD z scores are across time and reports the stability of spike–time relationships across pairs of neurons, displays a similarly narrow resonance effect between 4 and 10 Hz but maintains a near maximal value throughout this band (Fig. 8). We compare these results to the ones obtained when we optogenetically stimulated in vivo hippocampal networks (30). Rhythmic stimulation of parvalbumin-expressing (PV+) interneurons in PV::ChR2 transgenic mice was used to ensure that principle cells within the network were received subthreshold periodic inhibitory stimulation. Rhythmic optogenetic stimulation of PV+ interneurons leads to significant increases in both spike–LFP coherence and functional network stability for frequencies between 4 and 10 Hz among the principle cells within the network (Fig. 8). This suggests that in vivo CA1 hippocampal network stably organizes its firing activity within resonant frequency band of principal cells, while such organization is not observed when oscillatory drive is outside of this range.

Fig. 8.

Resonating networks have organized functional structure over a narrow frequency band. Theta band resonance leads to a highly organized functional network structure. In simulated networks, spike–LFP coherence (A), mean AMD z score (B), and functional network stability (C) all dramatically increase between 4 and 10 Hz. This effect is robust to noise, which is indicated by line color. (D) In vivo optogenetic stimulation of hippocampal PV+ neurons lead to similar increases in spike–LFP coherence and functional network stability at these frequencies. Error bars = SEM.

Discussion

We demonstrated in a biophysical model that shifting resonances facilitate learning of static and sequential patterns in neural networks. Our model combines subthreshold activation of neurons by stable and oscillating currents, which leads to firing in a narrow frequency band. The firing rate resonance of our model neurons displays an input-dependent broadening that allows for selective activation of subsets of neurons within a network. The resonance effect also leads to detailed mapping of a firing phase versus input relationship beneficial for the encoding of patterns into synaptic weights and for the autonomous termination of learning. The resonant effect at the single neuron level leads to the emergence of highly organized spike–time relationships at the theta band, which we have also shown in in vivo experiments.

The input-dependent broadening of the resonance curve in firing rate (Fig. 1) allows for selective activation of subsets of neurons within a network with increasing input frequency as has been demonstrated in other computational models, indicating this is a general property of neural networks with resonance (23). This provides a mechanism for networks to change representations by shifting the pattern of input strengths or, alternatively, by modulation of the oscillatory input frequency. Such a mechanism would operate similarly for both externally generated (i.e., sensory input) and internal (i.e., stored representations within synapses) inputs.

The mechanism we describe here can simultaneously promote both forward and reverse replay of recently learned sequences in neural networks, consistent with prior reports of replayed patterns in both directions, across even short intervals of in vivo recording (31). The reverse firing phase relationship and learning saturation seen in our external pattern simulations together provide a plausible mechanism for the generation of reverse replay events in vivo (Fig. 9). This mechanism relies on the fact that neurons with high input fire at early phases of oscillatory drive when in resonance. Before any synaptic change occurs, the firing phase is governed by the distribution of the external inputs the neurons receive. As learning progresses, neurons with the lowest external input strengthen their synaptic inputs more than the rest of the population, while highly activated neurons do the opposite, as shown in Fig. 5. The emerging pattern of synaptic connectivity is complementary to the input pattern, which leads to all neurons firing at the same phase (i.e., in synchrony). Synchronous firing leads to no net synaptic change and thus terminates learning. As the complementary input pattern is now represented within synaptic weights, in the absence of external input, neurons fire in the reverse order.

Fig. 9.

The model proposes a mechanism for the generation of reverse replay. Reverse replay due to how an input pattern imposes a phase procession of neuron firing due with respect to the oscillation. As the network learns the pattern, inputs to weakly excited neurons are strengthened while those to highly excited neurons are weakened. When the pattern is removed, inputs from synaptic connections dominate, and the reverse mapping of synaptic weights leads to reverse reactivation.

The mechanism for the replay of the reversed pattern is not dependent on the encoding frequency. Fig. 3 shows that a pattern can be encoded in one frequency band (6 Hz) and replayed at another (60 Hz), provided the neurons within the network can resonate at both frequencies. Such a mechanism could explain why sequential place cell activation during exploration (usually in the context of theta oscillations) can lead to subsequent replay events occurring in the context of higher frequency oscillations, such as a sharp-wave ripples (19). Here we use separate models to generate spiking responses to inputs of varying frequencies, but a neuron with resonances in both bands would behave in a similar manner. The importance of this frequency generality is that the encoding and replay of patterns in neural firing often occur when different frequencies are dominating the LFP. For example, sequences of place cell activation, and any synaptic encoding, occur when theta is the most dominant frequency band in the LFP, but instances of replay occur during sharp-wave ripples where gamma (40 to 100 Hz) is most prominent (14, 19).

Learning through STDP requires either saturation or compensatory plasticity mechanisms to counteract the inherent positive feedback effects on firing rate, leading to network instability. Previous implementations of STDP have used boundaries on synaptic weights, dynamic asymmetries between potentiation and depression, or renormalization of synaptic weights to preserve firing rates (reviewed in ref. 32). Our model proposes an alternative mode for preventing instability (Fig. 4). As the input pattern is encoded into synaptic weights and the firing phase distribution becomes more uniform, changes in synaptic weights decrease and stop due to features of the STDP curve around , which is a reasonable fit to experimental data (25). While many plasticity mechanisms exist both at the cellular and network level, the current mechanism provides an elegant solution to the question of when neural networks terminate learning of input patterns.

We have shown experimentally that predictions of our model agree with observed, network-wide pattern formation in hippocampal networks when channelrhodopsin-expressing PV+ interneurons are rhythmically stimulated (30). Within the hippocampus, functional network structure emerges and stabilizes during stimulation in the theta band (4 to 10 Hz) but not outside of it. Using several methods of measuring functional connectivity within networks, we found a robust resonance effect in the formation of stable network structure (Fig. 8). This effect is due to the organizing of the firing of the network around the phase of the oscillatory input. The fact that this effect is reproducible in various neuronal models (22) and also in vivo suggests that it may be a general feature of activity organized in neural networks, to optimize encoding of input patterns.

The input-dependent organization of network activity facilitated by resonance provides a network-level substrate for sequential learning (Figs. 6 and 7). When subsets of neurons have overlapping activation curves, the relationship between input and firing phase creates spike–time differences that are optimized for encoding the sequence order. One requirement for this result is that the activation of neurons needs to be skewed in time—in other words, repolarization occurs more rapidly than depolarization (Fig. 6). This ensures that connections strengthened by a balanced STDP regime are feedforward with respect to the sequence order, while feedback connections are weakened. Within the context of hippocampal place cell sequences, there is some evidence for this required skewness in activation (7, 33), though in an experience-dependent manner (34). Replay is the most direct readout of sequential learning. In the hippocampus, replay of place cell sequences occurs both in the forward and reverse direction (7, 18, 19, 21). These replay modes are represented in different proportions across behavioral states, with forward replay being more prevalent during sleep (21, 35). In our model, forward replay occurs when a network is driven by noise (i.e., randomly activated), and reverse occurs when the network is reactivated by oscillating input (Fig. 7 ).

Hippocampal place cells show a theta phase precession in their firing, as an animal approaches a location, neurons that code for a nearby place will fire in the troughs of the theta oscillation while those that code for a far place fire near the peak (6). This phenomenon has also been shown in the entorhinal cortex (16) and in the ventral striatum (17). In our model, neurons in resonance with an oscillating rhythm show a similar firing versus phase relationship.

Beyond the context of place cells, our model demonstrates how a network can translate information between the two main modes of neural coding rate (36) and phase (37–40) coding. Both rate coding, where stimuli are represented by the firing rate of neurons, and phase coding, where information is represented in the time differences between spikes, are observed in nervous systems. Rate coding is simple and reliable, however it is limited in its capacity for dynamic pattern separation (41). Our results provide a mechanism for the translation between these two coding schemes and allow for networks to switch through neuromodulation (23). Whether the mechanism described here mediates information encoding in the brain remains an open question. However, our present data suggest that such a mechanism has explanatory value for many of the observed in vivo phenomena surrounding learning.

Materials and Methods

Neuronal Network Model.

We use a network model that is composed of (or for the data in Figs. 6 and 7) excitatory neurons. Neuronal dynamics were based on a conductance-based model (Ks model) and governed by the current balance equation:The gating variables , , and were of the form . The slow potassium conductance, whose maximum value is , is largely responsible for the resonance displayed by this neuron model, and its value was set to 1.5 . Additional details of the neuronal dynamics can be found in ref. 42. Ks model neurons display a depolarization-dependent spiking resonance to subthreshold inputs in the 4 to 20 range.

Additionally, we used a second conductance-based neuronal model using the HH (28) model and parameters that resonated between 40 and 90 Hz to produce the data in Fig. 3. Membrane potential dynamics were governed by the current balance equation:The gating variables , , and evolved according to , where , , and other parameters are taken from ref. 28.

For both neuronal models, was split into two components. The first is (except for Fig. S3), which is identical for each neuron in the network. Cosine was chosen so when was set to zero (i.e., no oscillation), all neurons would receive the same peak current as DC. The second component was either , which is unique for each neuron, or in the case of data in Figs. 6 and 7 , which is a slowly varying activation current defined by the modified Gaussian function:where is the group to which a neuron is assigned (one of five groups), is the time of maximum activation of that group, ms is the width of the activation function, and is the skewness parameter. This leads to an activation time course that slowly grows to 227 and then rapidly decays to zero (Fig. 6).

Synaptic input was modeled as a double exponential conductance pulse with the dynamics:The decay constants, and , were set to 250.0 and 0.3 ms respectively. The synaptic delay constant, , was set to 0.08 ms and , where is the time of the last spike of the presynaptic neuron . is a synaptic multiple used to account for differences in the input resistance of the two neuronal models; it is set to 1.0 for the Ks model and 10.0 for the data in Fig. 3. The behavior of the HH model is robust to a range of values (Fig. S6). Total synaptic current to a neuron was defined as , where is 0 mV. Networks had a ring lattice structure and a connectivity rate of 6%. The connectivity scheme was small world and achieved through the Watts–Strogatz method with a rewiring probability of 0.2 (43).

Synapses evolved according to an additive STDP rule, where the weight change of a synapse between a presynaptic neuron and a postsynaptic neuron is defined by:Here , where is the time of the last spike fired by a given neuron. is the time constant of the effect of a spike decays and is set to 10 ms. is a symmetrical region around for which there is no synaptic change and is set to 1.5 ms. was the learning rate and was set to 20 for all simulations except in Fig. 4. Synapses were bounded in the region and initialized at 0.2 nS.

All numerical simulations were performed at a time step of 0.05 ms for the Ks model and 0.01 ms for the HH model using a fourth-order Runge–Kutta algorithm. All summary data take data from five realizations of the model, except for data in Fig. 3, which showed average SEM firing phase over 10 periods in one simulation.

Stimulation and Recording of Hippocampal Networks.

All procedures were approved by the University of Michigan Institutional Animal Care and Use Committee. Pvalb-IRES-CRE mice ((B6;129P2-Pvalbtm1(cre)Arbr/J; Jackson) were crossed to B6;129S-Gt(ROSA)26Sortm32(CAG-OP4*H134R/EYFP)Hze/J mice (Jackson) to generate PV::ChR2 mice, which expressed channelrhodopsin (ChR2) in PV-expressing (PV+) interneurons. By rhythmically activating these neurons in the hippocampus with 473 nm light, principle cells within the network were received subthreshold periodic inhibitory stimulation. For all recordings, PV::ChR2 mice ages 2 to 5 mo (n = 4) were anesthetized with isoflurane and chlorprothixene (1 mg/kg intraperitoneal injection). Mice were head-fixed, and a 1 mm × 1 mm matrix multielectrode [250 electrode spacing; Frederick Haer Co. (FHC)] was slowly advanced into CA1 until stable recordings (with consistent spike waveforms continuously present for at least 30 min before baseline recording) were obtained. An optical fiber was placed adjacent to the recording array for delivery of 473 nm laser light (CrystaLaser). Power output at the fiber tip was estimated at 3 to 10 mW for all experiments. CA1 neurons were recorded over a 15-min baseline period, after which PV+ interneurons were stimulated over multiple successive 15-min periods with a range of frequencies (2 to 18 Hz, 40 ms pulses). The various stimulation frequencies were presented in a random interleaved manner, during which neuronal activity continued to be recorded. Only those neurons recorded throughout the entire experiment were included in analyses of optogenetically induced spike-field coherence and network stability changes. For in vivo data, 80 and 68 neurons, respectively, met inclusion criteria for coherence and stability analysis. This dataset also appeared in ref. 30.

Functional Network Structure.

Functional network structure was calculated for both simulated and recorded networks in a similar manner. The first measure was spike wave coherence, which was calculated as the range of the spike-triggered average of the LFP over a window of 50 ms normalized by the peak amplitude of the LFP. In simulated networks, the LFP was the sum of all synaptic currents. This value ranges between 0, when spikes occur randomly in the LFP oscillation, and 1, when spikes always occur at the same time. In simulated networks, the LFP was the sum of synaptic currents.

The second measure of functional network structure was the stability of functional connections through time (30, 44). The basis of functional connectivity was the average temporal proximity of spikes between neurons and given by for the -th to -th neurons. Here is the time difference between the -th spike fired by neuron and the nearest spike fired by neuron . To determine whether neurons and are functionally connected, is compared with the null value given the firing rate of neuron and random firing of neuron by the score . The null distribution of MD is dependent on the interspike intervals (ISIs) of neuron . For an ISI of length , the first two moments of MD are and . We will find an ISI of length L within a spike train of length T with a probability of . Thus, of all of the intervals in the spike train of neuron , the expected value is . The expected SD is , where . To measure the stability of inferred functional connections, spiking data were separated into nonoverlapping time windows for which values were aggregated into matrices . Between adjacent time windows, cosine similarity, defined by , was used to quantify the change in functional network structure as a value between 0 (randomized) and 1 (no change). The stability of the functional network was quantified as the average similarity between adjacent time windows. Time windows were 2 s for simulated data and 1 min for recorded data.