For computational simplicity, the integrate-and-fire (IF) neuron model is used. The equation governing the IF neuron’s membrane potential is:

where τ_m (33 ms) is the membrane time constant, V_R (−70 mV) is the resting membrane potential, g_e is the synaptic input variable, E_e (0 mV) is the reversal potential for an excitatory synapse, R_m (200 MOhm) is the membrane resistance for the cell, and I is the applied DC current. Note that this equation for the IF neuron is a CUBA (current-based) model, following Vogels and Abbott (2005). When the membrane potential reaches the threshold value V_th (−54 mV), a spike is generated, and the membrane potential is reset to the resting membrane potential V_R.

Input connections to the IF neuron are mediated by exponential synapses, without delays. When a presynaptic spike occurs, the g_e parameter is incremented by an amount w (g_e → g_e + w). Otherwise, the g_e variable follows the equation:

where τ_e (5 ms) is the decay time constant for the exponential synapse. In these simulations, w is set to induce a subthreshold membrane potential oscillation, and a linearly spaced range of DC currents is then found to span the 1:1 phase-locking regime. It is important to note that the effect described in this work is extremely robust to the parameters of the output model neuron, and different parameter regimes have been tested to ensure this.

The 5000 presynaptic inputs to the IF neuron were modeled as inhomogeneous Poisson processes whose rate parameters oscillate in time, with a peak rate of 10 Hz. The input oscillation frequency is set to 20 Hz. To ensure that the lack of refractory period in the Poisson process inputs does not affect the results reported here, additional simulations have been performed with inputs that follow an inhomogeneous Gamma process.

The STDP rule is implemented in accordance with Song et al. (2000), with all-to-all spike pairings. Synaptic weights are restricted to range between 0 and w_max. For a given spike pairing with temporal difference s, the synapse between the pre- and post-synaptic neurons is modified according to (w → w + f(s) w_max). Thus, the value of the STDP kernel f(s) for a given temporal difference between a pre- and post-synaptic spike determines the percentage of w_max added to the synaptic weight (Song et al., 2000). Because we study a plastic system with intrinsic stability in the coupling between STDP and oscillations, it is not necessary to include a stabilizing mechanism in the STDP rule, provided that the learning rate A₊ is a reasonable value. Thus, we focus on an STDP rule with hard bounds (Song et al., 2000), so that we do not introduce extra stability arising from an STDP rule with soft bounds (van Rossum et al., 2000), though the combination of these two mechanisms could be interesting for future work.

These simulations were performed using the Brian simulator (Goodman and Brette, 2008, 2009), and a Python script to reproduce part of Figure 4 will be made available on ModelDB (https://senselab.med.yale.edu/modeldb/).

We first assume an STDP rule with linear, all-to-all spike pairings. Following previous theoretical work on STDP (Song et al., 2000), we formulate the STDP rule as a piecewise exponential function:

where s is defined as the difference in pre- and post-synaptic spike times (t_post − t_pre), τ₊ and τ₋ are the STDP time constants for potentiation and de-potentiation, A₊ is the amount of potentiation for an optimally potentiating PRE-POST pairing, and A₋ is the same for a de-potentiating POST-PRE pairing. The ratio of de-potentiation to potentiation (A₋τ₋/A₊τ₊) will be termed the STDP ratio. For the case of identical time constants considered here, we simply write (A₋/A₊). The value of the function f(s) determines the percent change for a synaptic weight on a given pairing, relative to w_max.

The rates of the presynaptic spike trains are approximated by a continuous sinusoid (in simulations, below, we model inputs as inhomogeneous Poisson processes):

where v is the angular frequency, r is the peak firing rate, and c is a parameter controlling the depth of modulation (max-minmax; Skaggs et al., 1996) of the oscillation:

The output spike train is formulated as a series of delta functions:

where φ is the phase offset of the output spikes. Here, the phase distribution of the output neuron is a delta function, reflecting the assumption of the 1:1 phase-locking regime; however, results will be similar for any moderately peaked, unimodal phase distribution.

To calculate the correlation function between the input oscillation and the output spike train S(t), one considers the integral:

where s is the temporal difference between pre- and post-synaptic spikes (t_post − t_pre). This integral simplifies to:

To calculate the expected weight change over time (dw/dt) as a function of output spike phase, we take the integral over the correlation function between pre- and post-synaptic spike trains multiplied by the STDP kernel:

From Eq. 12, one can clearly see that the location of the zeros depends on the angular frequency (v), the STDP time constants (τ₊ and τ₋), the depth of modulation parameter c, and the strength of potentiation and de-potentiation in the STDP rule (A₊ and A₋). One can think of the A₊ terms as phase advancing and the A₋ terms as phase delaying, for the weight change versus spike phase curve. The location of the zero crossing of this solution depends on the STDP ratio (A₋/A₊), and not on the absolute magnitude of A₊ or A₋, because increasing the magnitude of these terms while keeping the ratio constant corresponds to scaling of the function, and thus will not change the location of the zero crossings. Similarly, the firing rate of the inputs (r) does not affect the location of the zeros.

The result of the calculation (dw/dt as a function of output phase) is plotted in Figure 2A, with the parameter c set equal to 1 for full depth of modulation in the input sinusoid. The oscillation frequency f is 20 Hz, the STDP time constants are equal and set to 20 ms, A₊ is 0.01, and the STDP ratio is 1.05 (Bi and Poo, 1998; Song et al., 2000). Note that the form of this solution matches the behavior expected from the heuristic argument depicted in Figures 1C,D – net potentiation for the late-phase spikes, net de-potentiation for early phase spikes. It is also important to note that the extrema of this analytic expression fall near 90° and 270°, which are the points at which the input firing rate is changing fastest in time.

The analytic solution has two zero crossings – one stable, the other unstable. As spike phase is inversely related to input current in the 1:1 phase-locking regime for a model neuron with a positive phase response curve (PRC), it is straightforward to see that the early fixed point will be stable, with the later unstable. For an STDP ratio of 1.05, the location of this stable fixed point falls just after 180°. By changing the STDP ratio, however, stable phase-locking can be achieved throughout a wide range of the oscillation cycle. Thus, by increasing the STDP ratio (i.e., increasing the relative amount of de-potentiation), the point of stable phase-locking will occur later in the oscillation cycle; conversely, by decreasing the STDP ratio (i.e., increasing the relative amount of potentiation), the point of stable phase-locking will occur earlier in the oscillation cycle. The change in the shape of the analytic solution with different values of the STDP ratio is illustrated in Figure 2B, and the location of the fixed point as a function of the STDP ratio is plotted in Figure 2C.

The effect of the parameter c is detailed in Figure 3, in which the depth of modulation is varied from 100 (c = 1) to 40% (c = 4). In Figure 3A, the expected weight change is plotted as a function of output spike phase for an STDP ratio of 1.05, as in Figure 2A. In Figure 3B, the location of the stable fixed point is plotted as a function of the STDP ratio, for the three values of c. For input sinusoids with a smaller depth of modulation, a change in the STDP ratio has a greater effect on the location of the fixed point. The range, however, is similar for all cases.

To corroborate the analytic results, which approximate the input spike trains as a continuous sinusoid, with biophysically realistic input and output spiking processes, we performed computational simulations of an integrate-and-fire (IF) neuron and oscillating inhomogeneous Poisson process inputs. In these simulations, the IF neuron receives 5000 periodically modulated Poisson process spike inputs (see red inputs, Figure 1A). The number of inputs is chosen to be biologically realistic (Amaral et al., 1990), but also computationally feasible. Inputs are connected to the IF neuron by exponential, current-based synapses with a 5-ms decay time constant. Notably, as the analytic solution demonstrates, our STDP-based phase coupling mechanism is robust to the output neuron model. The only requirement is that the output model neuron has a positive PRC; thus, the IF neuron was chosen for computational simplicity.

In the control simulation without STDP (Figure 4A), a range of DC inputs are applied to the output neuron, so that the output will phase-lock to a range of phases of the input oscillation cycle. The mean phase of the output spike as a function of oscillation cycle over 50 trials for each DC input case is plotted in Figure 4A, with different colors for each DC input. It is important to note that after 1 s of simulation, stable phase-locking is achieved, and the initial transient has been integrated out, so that each neuron has reached an individual stable phase-locking point. This set of neurons with the various pre-STDP steady-state phase points can be interpreted as a population of cells, each with slightly different membrane parameters, such as membrane resistance or mean synaptic input.

The synapses connecting the oscillating inputs to the output model neuron express STDP with linear, all-to-all spike pairings. In these simulations, all synapses are initialized to the same weight. For the first 2 s of the simulation, STDP is turned off, during the initial 1 s transient. In this way, each neuron starts out at its individual stable phase. When STDP is turned on, our putative population of output neurons quickly converges to the unique stable spike phase predicted by the analytic solution in Figure 2B (prediction marked by black arrows, Figures 4B–D). The simulations in Figures 4B–D are done with STDP ratios increasingly weighted toward de-potentiation (1.05, 1.50, and 1.70), and correspond to the ratios for the analytic solution plotted in Figure 2B. For an STDP ratio of 1.05, the analytic solution predicts stable phase-locking at 185°; for STDP ratios of 1.50 and 1.70, stable phase-locking is predicted at 220° and 235°, respectively. Thus, by changing the STDP ratio, phase-locking can be achieved throughout the oscillation cycle.

To verify that a population of neurons will converge from outside the 1:1 phase-locking regime to our theoretically predicted stable fixed phase point, we made an additional test with an 800-neuron population receiving input from 10,000 excitatory cells. Each neuron is connected to an input with a 10% connection probability, giving an average of 1,000 input synapses per output neuron. The inputs oscillate, as above, but in this simulation, no extra DC input is added to the neurons. The simulation is run as follows. First, the simulation runs for 5 s, to integrate out any transients. Second, the simulation runs for an additional 5 s, to record the phase distribution before STDP (blue distribution, Figure 5). Third, STDP is applied during 30 s of simulation. Fourth, the simulation runs for a final 5 s, to record the converged phase distribution after STDP (green distribution, Figure 5). In comparing the recorded mean phase after STDP with the theoretically predicted value (cf. small vertical lines, Figure 5), we observe that the output phase distribution leaves the two spikes per cycle phase-locking regime and converges to within 1° of the theoretical value after 30 s of plasticity. The weight distribution of the inputs is depicted in the inset of Figure 5; as shown here, it is important that the learning rate (A₊) is relatively small, so that the weight distribution does not saturate to the hard bounds. The stability of the weight and phase distributions has been checked in longer simulations. Additionally, we have checked with simulations of neurons initially firing less than 1 spike per cycle that the mean of the phase distribution also converges to the theoretically predicted stable fixed point (data not shown).