Stochastic Synchronization in Purkinje Cells with Feedforward Inhibition Could Be Studied with Equivalent Phase-Response Curves

Simple-spike synchrony between Purkinje cells projecting to a common neuron in the deep cerebellar nucleus is emerging as an important factor in the encoding of output information from cerebellar cortex. A phenomenon known as stochastic synchronization happens when uncoupled oscillators synchronize due to correlated inputs. Stochastic synchronization is a viable mechanism through which simple-spike synchrony could be generated, but it has received scarce attention, perhaps because the presence of feedforward inhibition in the input to Purkinje cells makes insights difficult. This paper presents a method to account for feedforward inhibition so the usual mathematical approaches to stochastic synchronization can be applied. The method consists in finding a single Phase Response Curve, called the equivalent PRC, that accounts for the effects of both excitatory inputs and delayed feedforward inhibition from molecular layer interneurons. The results suggest that a theory of stochastic synchronization for the case of feedforward inhibition may not be necessary, since this case can be approximately reduced to the case of inputs characterized by a single PRC. Moreover, feedforward inhibition could in many situations increase the level of synchrony experienced by Purkinje cells.

understood in terms of its structure [1]. There is agreement that the cerebellum may play a role in a variety of cognitive functions, in addition to its involvement in motor control [2]. Figure 1 shows the basic anatomical organization of the cerebellar cortex. See [3] or [1] for reviews. In order to designate specific types of neurons and axons in the cerebellum the abbreviations of Fig. 1 will be used. Synapses from a source neuron/axon type toward a target neuron type will be denoted by the abbreviation of the source and target connected by a dash; e.g. PF-PC denotes the synapse between parallel fibers and Purkinje cells.
Perhaps the most influential set of ideas regarding how the cerebellum works is the Marr-Albus model [6,7], which has led to a variety of models in which Purkinje cells act like a perceptron whose learning signal comes from the climbing fibers (e.g. [8][9][10][11][12][13][14]). The discovery of conjunctive Long-Term Depression (LTD) in the PF-PC synapses [15][16][17] has added plausibility to these models. This type of LTD happens when elevated simple-spike activity occurs concurrently with climbing fiber inputs, weakening the synapses that caused the simple spikes, as hypothesized in the Albus model. Over time, however, it has become increasingly clear that conjunctive LTD alone may not explain learning in the cerebellum.
In most Marr-Albus models conjunctive LTD is the sole form of learning underlying cerebellar function. However, several studies suggest that conjunctive LTD is an incomplete explanation. First, it has been shown that cerebellar motor learning can take place in the absence of PF-PC LTD [18,19]. Second, the correlation of Purkinje cell firing and muscle EMG can show both positive or negative correlations, with positive correlations being the more prevalent [20,21]. If the role of Purkinje cells was to gate motor commands through just inhibition, negative correlations should be the most common. Third, synaptic inhibition of Purkinje cells, whose complex spikeelicited plasticity acts to counteract PF-PC conjunctive LTD, seems to play a role in motor learning [22], which is ignored by models that rely exclusively on PF-PC LTD. Moreover, other studies suggest that the timing of Purkinje cells' spikes is important, not only their firing rate. Tottering mutant (tg) mice have virtually the same firing rate as that of wild types during spontaneous activity and in response to optokinetic stimulation; nevertheless, tg mutants show abnormal compensatory eye movements and severe ataxia [23].
When explaining how the timing of PC simple spikes affect their DCN targets, and how different types of CF-mediated plasticity affect cerebellar output, it may be important to pay attention to synchrony among Purkinje cells innervating the same DCN cell. This synchrony can modulate the response of the target DCN cell [24][25][26][27]. It has been observed that there exists simple-spike synchrony among PCs separated by several hundred micrometers, and that this synchrony seems to depend on afferent input [28][29][30][31][32]. This synchrony does not seem to be fully explained by firing rate comodulation or PC recurrent collaterals. Firing rate modulation may be insufficient in this case, because there are cases where the modulation in synchrony is unrelated to the modulation in firing rate [29,30]. Purkinje cell recurrent collaterals tend to generate oscillations whose coherence decays with distance [33], which is inconsistent with the distances across which synchrony is found; moreover, it is unclear how sensory inputs could modify the functional coupling of Purkinje cells [31] if this coupling depended on a fast oscillatory regime. It is thus appropriate to study stochastic synchronization as a candidate mechanism to explain how Purkinje cells can activate synchronously.
The phenomenon of stochastic synchrony happens when several uncoupled oscillators synchronize their phases when receiving correlated inputs [34,35]. The intuition behind this is that if the oscillators become entrained to the inputs then they will respond similarly, thus acquiring similar phases. An interesting aspect of stochastic synchronization is that the degree of synchrony can be controlled by the way that the oscillators respond to inputs, which opens the possibility of its modulation by plasticity mechanisms. If synchrony plays a role in shaping the response of cerebellar cortex, it seems feasible that there are plasticity mechanisms capable of creating synchrony.
One possible reason why stochastic synchrony has not been largely considered in the case of Purkinje cells is the complication arising from the feedforward inhibition in the parallel fibers. As shown in Fig. 1, PFs stimulate MLIs, which in turn stimulate PCs. This inhibition has been observed as IPSPs arising shortly after EPSPs [36], and seems to be fundamental in understanding the response of PCs [37,38]. Although considerable advances have been made in understanding stochastic synchrony [34,35,[39][40][41][42][43][44][45][46][47][48][49][50], no study has explored how feedforward inhibition affects this process.
Exploring stochastic synchronization of Purkinje cells requires to represent their activity in terms of their phase. A neuron that fires periodically can be understood as a dynamical system whose trajectory in phase space follows an asymptotically stable limit cycle. Such a dynamical system can be described by a single variable called its phase; the phase describes how far the current state is along the limit-cycle trajectory. Perturbations to the system (such as synaptic inputs in the case of a neuron) can be described by how they shift the phase of the system when they are received [51,52]. The PRC (Phase Response Curve or Phase Resetting Curve) of the system plots the shift in phase that an input produces as a function of the system's phase when the input is received. PRCs are a standard tool when understanding the behavior of coupled oscillators, and have been extensively used to describe networks of neurons [53,54]. Also, as expected, PRCs are also a standard tool in analytical studies of stochastic synchronization.
This paper presents a prospective method to understand stochastic synchronization when feedforward inhibition is present by using an equivalent PRC. To model feedforward inhibition, the effect of feedforward excitation coming from the PFs to the PCs is represented with an excitatory PRC, and the effect of feedforward inhibition coming to the PCs from the MLIs is represented with an inhibitory PRC. For each excitatory input, a corresponding inhibitory input arrives after a fixed propagation delay. The equivalent PRC lumps the effect of both excitatory and delayed inhibitory PRCs so we only have one type of inputs, and yet the response of the oscillator to a spike train is similar to that when using both excitation and delayed feedforward inhibition. Studies on stochastic synchronization have not dealt with systems where there is feedforward inhibition, but the equivalent PRC allows their insights to be applied in this case.
The equivalent PRC is defined in two different ways. The first definition obtains an equivalent PRC such that the oscillator with feedforward inhibition, and the oscillator with the equivalent PRC have the same statistical distribution of their phases (this is made precise in the text). The equivalent PRC from the second definition is such that an oscillator using it and the oscillator with feedforward inhibition will have on average similar phases at particular points in time when they receive the same inputs. Monte Carlo simulations were performed to verify that oscillators using the equivalent PRCs respond similarly to an oscillator with excitatory and delayed inhibitory inputs. The results show that the oscillators with feedforward inhibition and those with equivalent PRCs produce output spike trains with a significant level of coherence, so if stochastic synchronization happens among oscillators with feedforward inhibition, it will most likely happen in oscillators with the corresponding equivalent PRCs. Furthermore, in most scenarios the effect of delayed feedforward inhibition would be a synchronizing one.

Models
The equivalent PRCs mentioned above will be developed in the next three subsections. An equivalent PRC can only emulate the effects of excitation and delayed inhibition in an approximate manner, and there are several ways to define it. This paper presents two different definitions for the equivalent PRC. The first subsection of this section presents introductory material, and the next two subsections develop the two different definitions of the equivalent PRC.

Reduction to a Phase Equation, and the Stationary Phase PDF
This subsection briefly outlines some basic results from [51,52] using notation based on [39]. The results show how N dynamical systems oscillating in a limit cycle and receiving impulsive inputs can be represented with N phase variables. PRCs, the phase transition function, and the phase evolution equation are introduced for individual oscillators, which allows one to find an equation for their stationary phase Probability Density Function (PDF). If we measure the phase of the oscillator at some random point in time, the phase PDF can provide the probability that the sampled phase is in a particular interval.
Consider N oscillators whose dynamics can be expressed aṡ for i = 1, . . . , N, where the vector X i (t) denotes the state of oscillator i at time t, F is the function describing the dynamics of each oscillator, and I i (t) represents external random inputs consisting of impulsive displacements in phase space. We assume that the oscillators have an asymptotically stable limit cycle X 0 (t). The impulsive inputs are given by where t i n represents the time of the nth input to the ith oscillator, and e i n provides the direction of the shift in phase space caused by the corresponding input, so that at time t i n the oscillator i receives an immediate shift in phase space from point X i to point X i + e i n . Since our oscillator represents a neuron, the value of e i n should be determined by the synapse that receives the input. Furthermore, it is assumed that the inputs will not take the system outside the basin of attraction of X 0 (t). The inputs that will be considered in this paper behave as Poisson random point processes. If an input has a mean firing rate r, then its interimpulse interval T has an exponential distribution: where P (T ) is the probability density function for T . Although Poisson point processes constitute a reasonable null hypothesis regarding the statistical structure of inputs to Purkinje cells, this hypothesis may not always be valid [55]. We define a phase variable θ along the limit cycle so that θ(t) = θ(X 0 (t)), θ has a constant angular velocity ω, and its range is [0, 1). This means that in the absence of external inputs we will haveθ In the first part of the results section we work with systems where all the oscillators have the same angular frequency ω. The phase of points not directly on X 0 is defined through the use of isochrons, which are the set of points that asymptotically converge to a particular trajectory X * 0 in the periodic orbit. If the points on an isochron converge to the trajectory X * 0 , then their phase at time t is θ(X * 0 (t)). For all functions in this paper whose arguments include a phase value, I assume that this phase value is taken modulo 1; e.g. phases 0, 1 are the same phase, and the same can be said of phases −0.1, 0.9, and 2.9. When an oscillator has phase θ at time t, and an input shifts its state at that moment by an amount e then this state moves to a new isochron, with the consequent shift in phase denoted by G(θ, e). We define the phase transition function as with its output taken modulo 1. If the phase of an oscillator at time t n right before its nth input is θ n , then the phase right after the input can be written as φ n ≡ F (θ n , e n ) = θ n + G(θ n , e n ). The evolution equation of the phase is an iterative equation describing the phase of the oscillator at the time when the (n + 1)th input arrives: where T n = t n+1 − t n . The dynamics of the phase in continuous time are described by the equationθ Given that the input times t n and the input effects e n are random variables, so are the phases θ n . The Probability Density Function (PDF) of the phase θ at time step n is denoted ρ(θ, n). This PDF can be described by the following generalized Frobenius-Perron equation [56]: The term W (θ − φ) is a density function for the input interimpulse interval T n expressed as a periodized change of phase with magnitude θ − φ (see Eq. (9)). Q(e) is the probability density function for e. Intuitively, the two innermost integrals produce the expected phase after the nth input with ψ being the starting phase; this expected phase, represented by the integration variable φ is taken in the outermost integral in order to calculate the probability that during the interspike interval the phase transitions from φ to θ . The transition kernel W (θ) can be explicitly obtained in the case of Poisson inputs, considering that its arguments are taken to be modulo 1: with A = r/ω. In the limit of a large number of transitions the PDFs ρ(θ, n) will reach a stationary state ρ(θ) obeying: I refer to this equation as the phase PDF equation.

The Equivalent PRC as a Function of the Phase PDF
The first idea to obtain an equivalent PRC, denoted in this subsection as Δ, is to take the phase PDF ρ(θ) produced by the excitatory and inhibitory inputs, and define Δ so its phase PDF matches ρ(θ). In practice, the formulas developed in this section allow one to take a numerically or experimentally measured phase PDF ρ(θ), and obtain a single PRC Δ so that with a single Poisson input, Δ produces a phase PDF similar to ρ(θ). The formulas from this section are tested by numerically obtaining the phase PDF created by an oscillator with feedforward inhibition, and comparing it with the phase PDF created by an oscillator using the equivalent PRC. I start by assuming a single oscillator and a single Poisson process that produces excitatory inputs. Feedforward inhibition is modeled by assuming that for each excitatory input at time t n there will be a corresponding inhibitory input at time t n + d, where d represents the feedforward delay. All excitatory inputs will produce a shift in phase space e exc , whereas inhibitory inputs produce a shift e inh . The excitatory and inhibitory PRCs are defined, respectively, as Δ exc (θ ) = G(θ, e exc ), Δ inh (θ ) = G(θ, e inh ), where the function G maps shifts in phase space to shifts in phase. An oscillator using instead the equivalent PRC will present a shift in phase Δ(θ ) at the times when the excitatory inputs arrive.
I assume a perturbation from the system where the PRC is zero for all phases and the phase PDF is uniform. Furthermore, the perturbation is small enough so that the phase transition function F is still invertible. Let ρ(θ) = 1 + ερ 1 (θ ), where ε is a scalar that determines the size of the perturbation. Expand the equivalent PRC as Δ(θ ) = εΔ 1 (θ ) + ε 2 Δ 2 (θ ) + O(ε 3 ) (using big Oh notation), so that as ε goes to zero we recover the unperturbed system. Before substituting these terms, the phase PDF equation (10) is simplified in three steps: 1. The middle integral disappears. There is only one type of input, so the PDF Q is a delta function that can be integrated out. 2. We perform the innermost integral using the basic formula for performing change of variables with Dirac δ functions. If g(x) is a real function with a root at x 0 the formula is δ(g(x)) = δ(x−x 0 ) |g (x 0 )| . 3. We substitute the transition kernel W for its expression in Eq. (9).
This yields Notice that the argument (θ − φ) of W is still taken modulo 1, so we need separate integrals for the cases when this argument is positive and negative: We now assume . Substituting these two previous expressions in the identity F (F −1 (φ)) = φ we obtain where the argument has been omitted in some functions for brevity of notation. Using this expression for F −1 we find where we use the fact that since F is invertible, F > 0. Using long division to eliminate the quotient and substituting into (12) yields The terms in this equation can be grouped according to the power of ε that they contain. The zeroth-order terms yield an identity corresponding to the unperturbed case. The equation corresponding to the first power of ε is This equation can be differentiated with respect to θ , and in the resulting equation we can solve for Δ 1 to obtain Integrating this provides an expression to evaluate Δ 1 : where C 1 is an integration constant. Notice that equation (14) only provides constraints for the derivative of Δ 1 , so it cannot be used to determine the value of C 1 , reflecting the fact that oscillators with different frequencies can have the same stationary phase PDF. Since the quantity εC 1 /A becomes a constant term in Δ, it adds an extra amount of advance or retardation to the phase whenever an input is received, and we can adjust its value so that the mean firing rate of the oscillator with two PRCs matches that of the oscillator with the equivalent PRC.
The equation corresponding to ε 2 in Eq. (13) is As in the previous case we can differentiate with respect to θ and solve for Δ 2 , which gives We can use the integration by parts formula ρ 1 ρ 1 = ρ 1 ρ 1 − (ρ 1 ) 2 to find the antiderivative of this expression, which is where C 2 is an integration constant. Alternatively, we can avoid having another integration constant by finding the definite integral of Eq. (17) from 0 to θ , obtaining Equations (15) and (19) provide the relationship between a phase PDF, and the PRC that lead to it. Notice that these results are consistent with those of [43], where a perturbative expansion is used to go from the PRC to the phase PDF, instead of going from the phase PDF to the PRC, as done here.
Monte Carlo simulations were performed in order to verify that an oscillator using the Δ approximation from the formulas above would respond similarly to an oscillator using Δ exc and Δ inh when provided with the same input, which consisted of a Poisson spike train with a frequency r = 600 Hz. This high rate comes from the assumption that the input comes from many homogeneous synapses receiving spike trains at a lower rate. Notice that the superposition of Poisson spike trains is also a Poisson spike train, whose rate is the sum of the rates for the superposed trains. The case of heterogeneous synapses will be treated further ahead. It is assumed that when the oscillator's phase transitions from 1 to 0 a spike is emitted; phase is not allowed to go from 0 to 1 due to inhibition. The shape of Δ exc was generated as a scaled version of the function 1 − cos(2πx), and Δ inh is a scaled version of cos(2πx) − 1, corresponding to type I PRCs [57]. The general shape of Δ exc is justified by the measurements that have been taken of Purkinje cells' PRCs at high frequencies [58], indicating that this is a positive, unimodal, type I PRC. The shape of Δ inh is justified by noticing that stimulation of MLIs tends to consistently induce an increase in the latency of the next spike in Purkinje cells [36,59]. Although the specific shapes of Δ exc and Δ inh used here may be somewhat different from that in real cells, the formulas in this paper can nevertheless be applied when the excitatory and inhibitory PRCs have any smooth shape. The phase PDF required in Eqs. (15) and (19) were obtained by sampling the phase whenever an input arrived. The phase histogram corresponding to these samples was normalized and smoothed with the Savitzky-Golay method. The derivative of the phase PDF was obtained using a 2 point rule with the smoothed phase PDF histogram. The constant C 1 of Eq. (15) was adjusted for each simulation up to two decimal places so that the firing rates of the oscillators with one and two PRCs would match for an input rate of 600 Hz.
The derivation of the formulas for Δ requires that the amplitudes a exc of Δ exc and a inh of Δ inh be moderate. The performance of the formulas gradually deteriorates with larger amplitudes. The values of a exc and a inh were chosen to be near the limit where the agreement obtained from using Δ is still acceptable. The feedforward delay d is taken to be 5 ms, which is somewhat larger than usual estimations of 1 or 2 ms. This long delay is used to test the limitations of the approaches used here to obtain an equivalent PRC, particularly the ones to be presented in the next subsection. Considering that in the absence of inputs the oscillators have a frequency of 50 Hz, a delay of 5 ms constitutes one quarter of the period, which is a considerable interval. Figure 2 presents the results of substituting the two PRCs Δ exc and Δ inh by the equivalent PRC Δ. Simulations were performed for three cases, each corresponding to a column in Fig. 2. The first case is when a exc < a inh , the second case when a exc = a inh , and the third case is for a exc > a inh . Each of these three cases produces a different shape for the equivalent PRC (panel b). Panel a shows that there is a good match between the stationary phase PDFs, as would be expected since the formulas for Δ were derived with this result in mind. The PDF curves come from 400 seconds of simulation, which permitted on average 240000 sample points.
It seems reasonable that the output spikes of the oscillators with one and two PRCs would come at similar times, since they have similar rates and phase PDFs. This is what is shown in panel c of Fig. 2, which shows the cross-correlograms of the output spike trains for both oscillators. The height of each bin in the cross-correlograms can be interpreted as the fraction of the spikes in the oscillator with two PRCs that coincide with a spike of the oscillator with one PRC. For comparison purposes sim-

Fig. 2
Comparison of oscillators with two and one PRCs when matching phase PDFs. Simulations for three different amplitude combinations of Δ exc and Δ inh . For all rows, the figures on the left correspond to a exc = 1/30, a inh = 2/30, figures on the middle correspond to a exc = 2/30, a inh = 2/30, and figures on the right correspond to a exc = 2/30, a inh = 1/30. The input rate for all simulations in a-d was r = 600 Hz. a: Stationary phase PDFs for oscillators with feedforward inhibition (red) and with their equivalent PRC (blue). b: PRCs used in the simulation. Δ exc (black), Δ inh (red), Δ (blue). c: Cross-correlograms of the oscillators' output spike trains with one and two PRCs. The vertical axis corresponds to normalized spike count per time bin, and the horizontal axis to time shift. Normalization was done by dividing the bin spike counts between the total number of spikes in the output of the oscillator with two PRCs. d: Cross-correlograms of the output spike train of the oscillator with two PRCs and a periodic spike train with the same mean frequency. Normalization was done as in c. e: Firing rate response of the oscillators with two (red) and one (blue) PRCs to five different levels of input rates. Input rates range from 240 Hz to 1200 Hz in constant increments of 240 Hz ilar cross-correlograms were produced in panel d, between the spike trains of the oscillators with two PRCs and a spike train with the same frequency and constant interspike intervals. It should be remembered that the parameter C 1 was adjusted so that both oscillators would have the same frequency, but from panel d it is apparent that the peaks in the top cross-correlograms are not just a product of periodicity in the signals.
Finally, I compared the response of both oscillators when the input firing rates were changed through the simulation. As can be seen in panel e, slight inaccuracies in the calculation of the C 1 parameter were amplified for larger firing rates. Moreover, in the case of balanced excitatory and inhibitory amplitudes the output firing rate tended to increase for larger input firing rates. This effect is amplified for larger PRC amplitudes, reflecting the fact that if the phase shift of Δ exc is large, the oscillator will spike before the inhibition arrives.
So far it has been assumed that all inputs of the same type (excitatory or inhibitory) will produce the same effect on the oscillator. On the other hand, a real neuron tends to have heterogeneous synapses. One way to represent this is to have separate excitatory and inhibitory PRCs for each input, representing different synapses. If an oscillator has N syn different inputs, we will have PRCs Δ i exc , Δ i inh , for i = 1, . . . , N syn . In this case, for each excitatory/inhibitory pair we may create an equivalent PRC Δ i . One approach to create Δ i is to consider the stationary phase PDF ρ i and output firing rate that would be produced if the inputs with Δ i exc , Δ i inh were considered in isolation. The formulas above could then be used to create Δ i . Instead of doing this, in the next subsection I develop a different way of obtaining Δ i , based on a more direct calculation of the phase-shifting effects produced by feedforward inhibition.

Obtaining an Equivalent PRC Using the Expected Inhibition
For the purposes of this paper, the ideal result of an equivalent PRC would be to have identical output spike trains for the oscillator with feedforward inhibition, and for the oscillator using the equivalent PRC. This would cause all the effects of stochastic synchronization to be the same for both oscillators. The method presented above to obtain Δ is based on creating an oscillator with a single PRC that has the same phase PDF as the one with two PRCs. However, having the same phase PDF and firing rate as the oscillator with feedforward inhibition does not guarantee that the output spike trains of the two oscillators will be the same. Indeed, given a phase PDF and a firing rate, there are many spike trains satisfying these two conditions. The equivalent PRC Δ defined above is therefore not necessarily optimal, and this raises the question on whether we can find an alternative definition for an equivalent PRC.
To obtain an equivalent PRC we could take a more direct approach, again trying to ensure that the phases of the oscillators with one and two PRCs are often the same. To develop this idea, let us consider two oscillators O ei and O eq . O ei experiences delayed feedforward inhibition, with excitatory and inhibitory PRCs Δ exc and Δ inh . On the other hand O eq has only an equivalent PRC Δ eq , and receives its inputs at the times when O ei receives excitatory spikes. At time t we can denote the phase of O ei by φ(t), and the phase of O eq by φ eq (t). Let us assume that at time t 0 an excitatory input is received by O ei , and the corresponding inhibitory input is received at time t 0 + d. The phase of O ei right before the inhibitory input at time t 0 + d can be denoted as φ(t 0 + d), and the phase right after the inhibitory input as φ . We could define Δ eq with the aim that the input at time t 0 creates a phase shift in φ eq such that φ eq (t 0 + d) = φ + (t 0 + d).
A brief consideration of this goal shows that it is impossible to achieve in an exact manner when the inputs are random. In the period between t 0 and t 0 + d there may be several inputs shifting the phase φ, and the magnitude of the inhibitory shift Δ inh (φ(t 0 + d)) depends on what the phase is in that moment. In other words, Δ inh (φ(t 0 + d)) is a random variable, and Δ eq must shift the phase φ eq (t 0 ) without knowing its value. Another important consideration about this goal is that the equivalent PRC that we are trying to obtain will generate a phase PDF that should be somewhat different from the one arising in the oscillator with two PRCs. This comes from the fact that at time t 0 the phase shifts from Δ eq and Δ exc will be different, so the time spent in the phases near φ(t 0 ) should be different.
The approach taken in this subsection to obtain Δ eq is based on the following idea: where is the expected value of the inhibitory shift that will occur at time t 0 + d. Obtaining this expected value is not trivial, so this will be done by a combination of analytical and numerical methods in two stages. The first stage provides a very simple version of Δ eq , which assumes no inputs between times t 0 and t 0 + d. The second stage goes through the trouble of finding an approximation to E[Δ inh (φ(t 0 + d))] given Poisson inputs, yielding a better version of Δ eq , which approximates Eq. (20). As it turns out, our work is not ready here; this is only the second of five Δ eq versions in this subsection. The third version of Δ eq is required because the equivalent PRC in Eq. (20) is insufficient to ensure that on average φ + (t 0 + d) = φ eq (t 0 + d). Although Eq. (20) ensures that the phase shifts resulting from the inputs at t 0 and at t 0 + d in O ei are on average the same as the phase shift at t 0 in O eq , we still have to consider the phase shifts from inputs between t 0 and t 0 + d. During this period the inputs will shift the phase of O ei according to Δ exc and Δ inh ; on the other hand O eq will shift its phase according to Δ eq at the time of excitatory inputs, resulting in different values for φ + (t 0 + d) and φ eq (t 0 + d). An iterative numerical method will thus be introduced to solve this deficiency, resulting in the third version of Δ eq .
The fourth version of Δ eq will be the equivalent of the second version (based on Eq. (20)), but for the case of heterogeneous inputs (e.g. many different inputs to the oscillator, each with its own excitatory and inhibitory PRCs). The fifth version of Δ eq will be the analog of the third version for the case of heterogeneous inputs.
We now begin obtaining the five versions of Δ eq .

First Δ eq Version
Let us first consider the case of an oscillator with a single Poisson input and feedforward inhibition. As before, it should be remembered that all phases are interpreted to be modulo 1. If an excitatory input arrives at time t n when the phase is φ n , then the phase will immediately experience a shift Δ exc (φ n ). At time t n + d the phase will be (φ n + Δ exc (φ n ) + ωd), where ω is the angular frequency of the oscillator. At that moment the inhibitory input will arrive, causing a phase shift Δ inh (φ n + Δ exc (φ n ) + ωd).
A simple way to define Δ eq for this case is This constitutes the first of the five versions for the equivalent PRC. One difficulty that quickly becomes apparent with it, is that in the time interval between t n and t n + d there will usually be other inputs arriving at the oscillator, so that the phase when the inhibitory input arrives will generally not be (φ + Δ exc (φ) + ωd). Indeed, this approach only produces reasonable results when inputs are unlikely to arrive between t n and t n + d, which could happen when the value of d is very small.

Second Δ eq Version
One way to improve our equivalent PRC is to substitute Δ inh (φ + Δ exc (φ) + ωd) by the expected value of the inhibitory shift given the phase when the excitatory shift happened, as was done in Eq. (20). Let θ(t) be the function that gives the phase of the oscillator at time t. Assume that an excitatory input arrives at time t 0 , when the phase is φ 0 , meaning φ 0 = θ(t 0 ). Furthermore, assume that between the times t 0 and t 0 + d there arrive k e excitatory inputs at the times t e j , j = 1, . . . , k e ; and k i inhibitory inputs at the times t i m , m = 1, . . . , k i . For these particular initial phase and inputs define the phase deviation as D is a random variable that tells us how much the phase will change due to inputs during the interval between t 0 and t 0 + d. For notational convenience let us define a ≡ φ + Δ exc (φ), and b ≡ φ + Δ exc (φ) + ωd. Our goal is to calculate the expected value of Δ inh (φ + Δ exc (φ) + ωd + D), which is denoted by E(Δ inh (b + D)|φ). This notation indicates the expected value of the inhibitory shift given that the excitatory shift happened when the phase was φ. The equivalent PRC can then be defined as: This is the second version of the equivalent PRC in this subsection. The following paragraphs deal with finding a practical way to calculate E(Δ inh (b + D)|φ), culminating with Eq. (27), which can be used in conjunction to Eq. (24) to obtain an approximation for this expected value.
To calculate E(Δ inh (b + D)|φ), we can start by calculating this expected value when we know exactly how many excitatory and inhibitory inputs arrived during the delay period. Assume that the inputs are independent Poisson point processes, with rate r e for the excitatory ones, and rate r i for the inhibitory ones. Using the PDF for the Poisson distribution we can obtain where E(Δ inh (b + D)|φ, k e , k i ) is the expected value of the inhibitory shift given that there were k e excitatory and k i inhibitory inputs during the delay period, making no assumptions about the order in which they arrived. The assumption of independence between excitatory and inhibitory inputs is based on the fact that we are restricted to a time interval of length d, during which none of the inhibitory inputs is the result of feedforward inhibition from one of the excitatory inputs. Notice that the first two factors decay exponentially, so in practice it is only necessary to use a moderate number of terms. The strategy to obtain E(Δ inh (b + D)|φ, k e , k i ) is to first find the PDF of D, so we can then find the expected value through integration. This calculation can become involved, so I will first focus on the simpler case when there is only a single excitatory input and no inhibitory inputs during the delay interval. Under these circumstances, if the initial excitatory stimulus arrived at phase φ, the PDF of D is denoted by p(D|φ, k e = 1, k i = 0). What follows is some formal reasoning to arrive at an expression for E(Δ inh (b + D)|a, k e = 1, k i = 0). The reader may just go directly to Eq. (25), which is intuitive enough.
Let B denote the Borel sets in the interval [0, a exc ], where a exc is the largest value on the range of Δ exc , and define a function I e : B → B that maps each set A ∈ B to its preimage under Δ exc . Given the fact that there was only a single input coming from the Poisson process in the phase interval [a, b], I make the assumption that the input could have arrived with equal probability at any phase between a and b.

This implies that for an interval H in [0, a exc ] the probability of D ∈ H is given by λ(I e (H ) ∩ [a, b])/λ([a, b])
, where λ is the standard Lebesgue measure for the real numbers. Put into words, the probability that the phase advance D produced by the sole excitatory input is in the set H is the probability of the phases that are mapped into H by Δ exc .
Notice that λ([a, b]) = ωd. If we define a measure μ(H ) = λ(I e (H ) ∩ [a, b])/ωd, then the PDF of D will be the Radon-Nikodym derivative of μ with respect to λ. A practical way to calculate this PDF starts by partitioning the interval [0, a exc ] into subintervals were Δ exc is invertible or constant, which should be possible for any reasonable PRC. If Δ exc is invertible in the interval [α, β], and x = Δ exc (α), y = Δ exc (β), then dλ.
The change of variables formula shows that ds.
If Δ exc = 0 in [α, β] in sets other than those of zero measure we can write P D ∈ [x, y]|φ, k e = 1, which means that |Δ exc (Δ −1 exc (s))| −1 /ωd is the PDF of D. If we have an interval where Δ exc is equal to a constant c, the integral becomes undefined due to a zero in the denominator. In this case we have P (D = c|φ, k e = 1, k i = 0) = μ(c), and finding the expected value of the inhibition is trivial. For invertible intervals we can now write the expected value of the inhibition as Using a change of variables this becomes the more intuitive formula Notice that after the change of variables the integral can handle the case where Δ exc is constant. Considering that the value of an integral on the interval [a, b] is the sum of the integral on the subintervals [α, β] where Δ exc is invertible or constant, we can now write In a similar manner it can be shown that The complexity of these equations increases once we have more than one input, and once we have both excitatory and inhibitory inputs, because the order in which they arrive is important. In this case the expected value for the inhibition comes from averaging over all the possible phases when the first and second stimuli could have arrived, and over the possible orders for the arrival of stimuli. To illustrate this, let us look at the formula for k e = 1, k i = 1 Intuitively, the integration variable ψ 1 stands for the phase when the first input arrived, and ψ 2 for the phase when the second input arrived. The first pair of nested integrals are for the case when the excitatory input happened first, and the second ones for the case when the inhibitory input was the first to arrive. The innermost integrals obtain the average inhibition given that the first stimulus arrived at phase ψ 1 , and in the limit when ψ 1 → b, they converge to the Δ inh expression with ψ 1 substituted by b, and ψ 2 substituted by b + Δ exc/inh (b).
In order to write the formula for the case with arbitrary values for k e and k i some preliminary definitions are required. Notice that if we have k e excitatory and k i inhibitory inputs, there are C k e +k i k e different input sequences according to whether the j th input was excitatory or inhibitory. Let us index those sequences and denote them by σ i . In other words, we create C The general formula can now be written as For this equation I have used the convention of writing the differential sign next to its corresponding integration sign. Although Eq. (26) expresses the expected inhibition values that we want to calculate, its complexity makes it virtually useless for practical purposes. Fortunately, a simple assumption can simplify this expression. Assume that for each input sequence, the inputs happen at regular time intervals (the time periods between any two inputs are equal). It is simple to calculate the expected value of the inhibition for this case. If we have K = k e + k i inputs, define γ = ωd/(K + 1), and for i = 1, . . . , C k e +k i k e let θ i 0 = a, We then have where RT stands for "Regular Times," meaning that the inputs arrive at regular time intervals. For a smooth function Δ inh and relatively small values of the feedforward delay d we will have

Third Δ eq Version
Even if we now can obtain a good approximation to the expected phase shift caused by feedforward inhibition, the equivalent PRC from Eq. (23) may still not achieve the goal of reaching, on average, the same phase as the oscillator with two PRCs after the feedforward delay. To make this explicit, assume that the function θ(t) provides the phase of the oscillator with feedforward inhibition at time t, just as it is described for Eq. (22), and let θ eq (t) give the phase of an oscillator using the corresponding equivalent PRC from Eq. (23) when the input times are the same. Using an equivalent PRC instead of Δ exc and Δ inh causes the phase deviation of Eq. (22) to become Δ eq θ eq t e j .
In general, D = D eq during the delay period; we can calculate the expected phase difference that this will cause right after the feedforward delay. If an initial excitatory input arrives at time t when the phase is φ, the expected value of the phase for the oscillator with two PRCs at time t + d right after the feedforward inhibition is On the other hand, the expected value of the phase for the oscillator with one PRC at Subtracting the previous two expressions gives us the expected value of the phase difference between the two oscillators at time t + d given that there was an excitatory input at time t when the phase was φ, denoted by ξ(φ): If we are capable of calculating E(D|φ) and E(D eq |φ), then we can use ξ(φ) in order to create an equivalent PRC that produces a smaller value of ξ(φ). A simple algorithm for doing this is as follows. Let M be an integer, and a small real number.
Define Δ (0) eq as the equivalent PRC from Eq. (23), and D 0 eq as its corresponding phase deviation.

end for
At each step in this algorithm the functions ξ (i) and Δ (i) eq are calculated for all the values of φ, so it is similar to gradient descent performed for a whole function. The resulting PRC Δ (M) eq constitutes the third version of an equivalent PRC we have obtained, and can already provide very good results for oscillators with a single input, or with homogeneous excitatory and inhibitory PRCs.

Fourth Δ eq Version
The fourth version of an equivalent PRC that I will obtain extends the second version to the case when there are heterogeneous inputs. The reason why the approach taken so far to obtain Δ eq may fail when we consider several types of inputs, each with its own excitatory and inhibitory PRCs, is that when Eq. (27) is obtained the phase is assumed to advance at a steady rate between inputs (the phase would advance an amount γ between inputs). If there is only one type of input this is justified, since the oscillator has a constant angular frequency. When we have different types of inputs we consider each one separately, so even if inputs of one type arrive at regular intervals in time, the phase will be shifted between consecutive times by inputs of other types. This will become explicit in the following calculation.
Consider an oscillator with feedforward inhibition that receives N syn different inputs. We consider that for each input there are two "synapses," one excitatory and one inhibitory, characterized by the PRCs Δ i exc , and Δ i inh for the ith input. We need to obtain N syn equivalent PRCs, with Δ i eq being used to replace Δ i exc and Δ i inh . The goal is therefore to obtain a version of Eq. (27) that works for each ith input individually. In order to model how the oscillator's phase changes between repetitions of the ith input I will use the concept of variable phase velocity, which will be explained next.
Let us say an oscillator has a non-constant phase PDF ρ(φ). We can think that this oscillator's phase has a constant rate of change dφ/dt = ω, but its inputs reshape the phase PDF so it becomes ρ(φ). Alternatively, we could image that the oscillator receives no inputs, but instead has a phase whose velocity dφ/dt = ω(t) is changing over time so that ρ(φ) is produced. The idea to be introduced here is to use this oscillator with no inputs and variable phase velocity in order to model the oscillator with constant phase velocity and random inputs.
Let T denote the period that we will assign to our oscillator. Assume that at time t = 0 the oscillator has phase 0, and denote the time it will take to reach phase φ ∈ [0, 1] by τ (φ). It is easy to show that when τ (φ) = T φ 0 ρ(ψ) dψ the phase PDF of the oscillator is ρ(φ). Since ρ > 0, τ (φ) is a monotone increasing function, which implies that it will be invertible on [0, 1]. Let P denote this inverse. The function P provides the phase as a function of the time since phase 0 was crossed. We consider the argument of P to be modulo T , so it is always in the interval [0, T ]. We are now ready to create a version of Eq. (27) for the case of heterogeneous inputs.
Let ρ(φ) be the phase of the oscillator with feedforward inhibition and heterogeneous PRCs. Assume that during the feedforward delay period the synapses for input j receive k j e excitatory inputs, and k j i inhibitory inputs.
where the * symbol next to the expected value denotes the fact that we used the variable phase velocity approach. The definition of the phase deviation D in this equation reflects the case of heterogeneous synapses: Obtaining the equivalent PRC proceeds as before.
where r j e and r j i are the excitatory and inhibitory firing rates of the j th input. Since we are considering the case of feedforward inhibition, we will have r j e = r j i . The fourth version of the equivalent PRC is

Fifth Δ eq Version
As would be expected, the fourth version of Δ eq has the same limitations as the second version, and those limitations can be surmounted using the same approach that led from the second to the third equivalent PRC. The phase deviation when using the equivalent synapses from the fourth PRC is k j e l=1 Δ j eq θ eq t e l,j .
If at time t an excitatory input is received by the j th synapse, the expected phase for the oscillator with feedforward inhibition and heterogeneous synapses at time t + d is The notation E(D|φ, j ) indicates the assumption that a shift of magnitude Δ j exc (φ) happens at time t. In the case of the oscillator with the equivalent PRCs of Eq. (29), the expected phase at time t + d is Similarly to Eq. (28), the difference in phase at time t + d is Let the functions Δ j,(0) eq come from Eq. (29), for j = 1, . . . , N syn . The procedure to obtain the fifth version of the equivalent PRC is as follows:

end for end for
The value E(D (i−1) eq |φ, j ) is the expected phase deviation calculated using the equivalent PRCs Δ j,(i−1) eq . This calculation can be time consuming depending on the method used, so the algorithm can be modified by dividing the N syn equivalent PRCs into N batch batches. Then for each batch we use the same ξ j,(i) j in order to update the PRCs. Fig. 3 Comparison of the equivalent PRCs for homogeneous inputs. Result of a simulation using four types of oscillators. The first one has both excitatory and inhibitory PRCs of amplitudes a exc = 2/30 and a inh = 2/30, respectively, with a feedforward delay of 5 ms. The second one has the equivalent PRC Δ from the previous subsection (Eqs. (15) and (19)), labeled "phase PDF Eq. PRC." The third one has the second equivalent PRC from this subsection (Eq. PRC 2). The fourth one has the third equivalent PRC from this subsection (Eq. PRC 3). The three oscillators with the equivalent PRCs only get stimulated at the time of the excitatory inputs, which have a rate of 600 Hz. The other oscillator experiences feedforward inhibition. a: Shape of the PRCs. The dashed black and red lines show Δ exc and Δ inh , respectively. The blue line corresponds to the PRC Δ from the previous subsection. The green line corresponds to Eq. PRC 2, and the magenta line to Eq. PRC3. b: Phase PDFs for the four oscillators. The red line (partially occluded by the blue line) shows the phase PDF for the oscillator with feedforward inhibition. The blue line corresponds to the oscillator with Eq. PRC Δ from the previous subsection. The green line is from the oscillator with Eq. PRC 2, and the magenta line for the oscillator with Eq. PRC 3. c: Average firing rate for the four oscillators. d-f: Phase difference density (see text) for the three oscillators with equivalent PRCs. g-i: Cross-correlograms between the output spike trains of the three oscillators with equivalent PRCs and the output spike train of the oscillator with feedforward inhibition. j: Cross-correlogram between the output spike train of the oscillator with feedforward inhibition and a periodic spike train with the same frequency. k: Circular variance (see text) for the three oscillators with equivalent PRCs

Computational Simulations
As before, I use computational simulations in order to validate the approximations obtained in this subsection. First, in Fig. 3 I compare the second and third equivalent PRCs from this subsection with the equivalent PRC Δ obtained from Eqs. (15) and (19) in the previous subsection. The performance of the first PRC from this subsection is significantly worse than any of these three and is not shown. We have at this point several methods to obtain an equivalent PRC, and the point of Fig. 3 is to see how those methods measure against each other. To make this comparison, we can simulate an oscillator with feedforward inhibition, and oscillators with the different equivalent PRCs, all of these receiving inputs at the times when the excitatory inputs arrive for the oscillator with feedforward inhibition. We can then see which of the oscillators with an equivalent PRC had the phases and output spike train that best match those of the oscillator with feedforward inhibition.
Panel b of Fig. 3 shows that using Δ leads to a much better agreement in phase PDFs than when using the second or third equivalent PRCs. This is not surprising, since Δ was defined with this purpose in mind. What is remarkable is that despite this, the third PRC provides a greater coherence in the output spike trains, as will be shown below.
In Fig. 3c we can see that Δ and the third equivalent PRC can match the firing rate of the oscillator with two PRCs. In the case of Δ this is the result of adjusting one parameter, but in the case of the third Equiv. PRC this is a consequence of the iterative procedure that leads from the second to the third version of Δ eq . The inability of the second PRC to match firing rates stems largely from the fact that it doesn't account for the correct average phase shifts during the feedforward delay period.
We have stated that an ideal result would be for the oscillator with feedforward delay and the oscillator with the equivalent PRC to produce the same output spike trains. We should therefore compare the output spike trains between the two oscillators. A graphical way of doing this is with cross-correlograms, as in Fig. 2c, d. Panels g-i of Fig. 3 show cross-correlograms between the output spike train of the oscillator with feedforward inhibition and the oscillators with equivalent PRCs. In addition, Fig. 3j shows the cross-correlogram of the spike train from the oscillator with feedforward inhibition, and a periodic spike train of the same frequency.
A perhaps simpler way to visualize the coherence between the output spike trains comes from the phase difference density, shown in Fig. 3d-f. To calculate this, each time the oscillator with feedforward inhibition spiked, the phase of the other oscillator was measured. This resulted in a sample of oscillator phases, whose normalized histogram constitutes the phase difference density plots shown in panels d-f. Notice that if the two oscillators always spiked at the same time, then the phased difference density plot would have a single tall and narrow peak at zero, and if the two oscillators spiked at uncorrelated times the plot would be flat. As with the cross-correlogram, it can be seen that the largest peaks are obtained with the third equivalent PRC. A way to quantify this comes from the circular variance measure in panel k.
Circular variance [60] is a procedure used to obtain the variance in a sample whose data points are periodic. In this case, since the phases 0 and 1 are the same (and the phases 0.99 and 0.01 are very close), using the normal procedures for averaging produces unwanted results. To calculate the circular variance we start with a sample of phases φ i , i = 1, . . . , N, in this case the phases collected to produce the phase difference density plots. Since the phases are periodic, in order to average them they are represented as complex numbers with unit norm. If the number corresponding to φ i is s i , then R{s i } = cos(φ i ), I{s i } = sin(φ i ), and the average of these numbers is The circular variance is defined as Z ≡ 1 − s . From this definition it can be seen that smaller values of Z indicate a larger degree of coherence, and taller, narrow peaks in the phase difference density plot.  Figure 3k shows that the circular variance is significantly smaller when using the third equivalent PRC. This is results shows, at least numerically, that a close match in the phase PDF of two oscillators does not necessarily entail an optimal match in their output spike trains, even when their firing rates are very close.
The next result shown is for the second and third versions of Δ eq . The third version of the equivalent PRC can give accurate results in replicating the response of an oscillator with feedforward inhibition, especially when the feedforward delay is small. If this delay is small enough, even the second version of the equivalent PRC (Eq. (23)) can give reasonable results. This is seen in Fig. 4, where the feedforward delay is taken to be 1 ms. The simulations used to produce Figs. 4, 5, and 6 were similar to those presented in the previous subsection, the only difference is in the equivalent PRCs used, and in the case of Fig. 4, the length of the feedforward delay.
In panel a of Figs. 4, 5, and 6, it is again seen that the matching of stationary phase PDFs is not as good as in the previous case (shown in Fig. 2), when the formulas where derived for this purpose. This is to be expected; for each excitatory input the oscillator with feedforward inhibition advances its phase, whereas the oscillator with  Fig. 4 is the delay value of 5 ms. a: Stationary phase PDFs for the oscillator with feedforward inhibition (red), with the second version of the equivalent PRC (green), and with the third version of the equivalent PRC (blue). b: PRCs used in the simulation. Δ exc (black), Δ inh (red), second version of Δ eq (green), third version of Δ eq (blue). c: Expected value of the inhibition as a function of the phase when the initial excitatory input arrives. The red curve comes from simulating the oscillator with feedforward inhibition, and the black curve comes from the approximation in Eq. (27). d: Cross-correlogram between the output spike trains of the oscillators with feedforward inhibition and with the third equivalent PRC. The horizontal axis indicates time shift, and the vertical axis the normalized spike count. Normalization was done by dividing the bin spike counts between the total number of spikes in the output of the oscillator with two PRCs. e: Firing rate response of the oscillators with feedforward inhibition (red), the second version of Equiv. PRC (green), and the third version of Equiv. PRC (blue) to five different levels of input rates. Input rates range from 240 Hz to 1200 Hz in constant increments of 240 Hz . Magenta dash-dotted lines indicate the phase PDF corresponding to intermediate iterations in the iterative procedure used. b: PRCs used in the simulation. Different colored lines show the equivalent PRCs for the first five synapses. c: Cross-correlogram between the output spike trains of the oscillators with feedforward inhibition and with the equivalent PRCs. The horizontal axis indicates time shift, and the vertical axis the normalized spike count. Normalization was done by dividing the bin spike counts between the total number of spikes in the output of the oscillator with two PRCs. d: Cross-correlograms between the output spike train of the oscillator with feedforward inhibition and a regular spiking oscillator matching its mean frequency. Normalization of the spike counts per bin was done as in c. e: Firing rate response of the oscillators with feedforward inhibition (red), and the equivalent PRCs (blue) to five different levels of input rates. Input rates range from 240 Hz to 1200 Hz in constant increments of 240 Hz the second and third equivalent PRCs experience a shift in phase that includes an estimate of future inhibition. On the other hand, when tracking the timing of output spikes, the oscillator with the third version of the PRC tends to outperform the oscillator with a phase PDF-matching PRC. In panel d of Fig. 4, the normalized crosscorrelogram shows a large similarity between the output spike trains. Simulations with the first equivalent PRC do not produce such a good cross-correlogram, even when the feedforward delay is 1 ms (data not shown). As before, to show that the peaks in the cross-correlogram come from a similar temporal structure in the spikes from both oscillators, and not just from them having similar periods, I also show in panel e the cross-correlogram between the spike train with delayed inhibitory inputs, and a spike train with constant interspike intervals and the same frequency.
Another relevant result is in panel c of Fig. 4. This panel shows a close agreement between the expected inhibitory shift at time t + d obtained from simulations (red curve), and from Eq. (27) (black curve), which approximates the expected inhibition by assuming that during the inhibitory delay all inputs will arrive at regular intervals.
Reported values of delay between an excitatory input and the corresponding feedforward inhibition are usually in the 1-2 ms range [36,59], and have even been described as non-existent [61]. It is nevertheless relevant to test whether the formulas of this paper can be still applicable when the feedforward delay is not as short. The approaches taken here to derive the equivalent PRCs should show their shortcomings as the delay and the input firing rate are increased. With this in mind the simulations in this paper-with the exception of the one in Fig. 4-use a delay of 5 ms, larger than reported values, but still biologically plausible. Figure 5 illustrates simulations done with the second and third versions of Δ eq and a feedforward delay of 5 ms. It is evident in panel e that the second version of the PRC becomes incapable of matching the firing rate of the oscillator with feedforward inhibition as the input rates become larger. On the other hand, the oscillator with the third version of Δ eq still displays similar behavior to the oscillator with feedforward inhibition. It should be noted that the third PRC version only requires to specify three parameters: the maximum number of inputs to be considered during the delay period (limiting the number of terms in Eq. (24)), and for the iterative procedure, the number of iterations M, and the learning rate . In the case shown in Fig. 5, there were up to seven inputs during the delay period, and M = 10, = 0.04. For the case shown in Fig. 3, up to eight inputs were considered during the delay period, with M = 18, = 0.02. The fifth version of the equivalent PRC was also tested, using 60 different inputs, each one with different amplitudes for its excitatory and inhibitory PRCs. The result of the simulations is illustrated in Fig. 6. As can be seen from this figure, and from Fig. 5, the performance for oscillators with homogeneous and heterogeneous inputs is very similar. The fifth PRC version uses the same three parameters as the third version, and optionally, a batch size for the iterative procedure. For the simulation in Fig. 6, M = 10, = 0.04, and up to three inputs were considered during the delay period.

Discussion
Experimental results point to the possibility of Purkinje cells synchronizing due to correlated inputs. Comodulation of firing rates is insufficient to explain this, but we can resort to the emerging mathematical theory of stochastic synchronization. At the outset, however, it is unclear that this theory is applicable here due to the presence of feedforward inhibition. Rather than obtaining a theory of how uncoupled oscillators with feedforward inhibition synchronize, I suggest to define equivalent PRCs meant to show that the systems with feedforward inhibition, which have excitatory and inhibitory PRCs, behave in a similar way to a system with only one PRC. The equivalent PRC is what in practice can be used to study the synchronization properties of oscillators. The first type of equivalent PRC defined here can be used to produce any given phase PDF (assuming Poisson inputs), and in particular, can produce the phase PDF that arises from a system with feedforward inhibition. With the second type of PRC defined here we can obtain an oscillator that, on average, has the same phase as the oscillator with feedforward inhibition after receiving the inhibitory stimulus. Unlike the first type of equivalent PRC, the second one requires knowing the shape of the excitatory and inhibitory PRCs in order to be computed. The phase PDF, excitatory PRC, and inhibitory PRC of Purkinje cells can all be experimentally measured, although the exact methodology is beyond the scope of this paper.
An interesting insight arising from the equivalent PRC is that feedforward inhibition could have the effect of enhancing synchronization. Type II PRCs, which are positive for late phases and negative for early phases [57], are the best at synchronizing due to correlated inputs [41,42,44,46,48]. As can be observed from Figs. 3, 4, 5, 6, feedforward inhibition has the effect of adding a negative region to the equivalent PRC which makes it of type II. Although the shape of inhibitory PRCs is still unknown, it seems reasonable that this result will hold for most reasonable shapes, inhibition delays, and ratios of excitatory to inhibitory amplitudes.
A final observation is that comparing the results from the two types of equivalent PRCs (Fig. 3) leads to the curious conjecture that an optimal equivalent PRC (in the sense of maximizing coherence with the output spike train of the oscillator with feedforward inhibition) need not also produce a matching phase PDF.

Is the PRC Model Adequate?
Purkinje cells are particularly complex, and their computational models tend to be mathematically intractable. The oscillator representation is a mathematically tractable model that allows one to construct complicated hypotheses that may then be addressed by physiological experiments and by detailed models. This extra step in the modeling process is beneficial in the study of synchronization because the effects of changing physiological parameters are not straightforward in this case. Still, modeling Purkinje cells as oscillators and characterizing their inputs through a PRC entails a drastic reduction in complexity, which could obscure important details. Purkinje cells tonically fire simple spikes with frequencies ranging from 30 to 150 Hz; these frequencies are modulated by afferent and efferent inputs (e.g. [62,63]), creating a wide dynamic range. It is likely, however, that Purkinje cells are intrinsically regular, and if they show irregular inter-spikes intervals in vivo, this is because of MLI inputs ( [37]). Moreover, we may not need to consider the complex morphology of PC dendritic trees, as the summation of inputs may happen independently of their location and distribution ( [64]). Indeed, in the dendritic tree there exist voltage-gated calcium channels that amplify distant focal inputs more than proximal ones, canceling cable attenuation and making the somatic response largely independent of input location.
Further considerations regarding the general suitability of the PRC representation are presented in [54]. For the specific case of Purkinje cells, the main concerns when using the theory of stochastic synchronization may be the variability in firing rates, and the fact that the PRCs change depending on this firing rate. These issues are not fully resolved, and a task for the future may be to test possible answers to them.
It is the author's opinion that the results in this paper will become relevant only if synchronization by virtue of correlated inputs is found to exist to a significant degree among Purkinje cells projecting to the same nuclear cell. This may be the next question to be addressed, and it can start by working with more realistic computational models of Purkinje cells. Realistic modeling of cerebellar cells has been taking place for decades, and there are several models that could be used for this endeavor (e.g. [25,38,65,66]). These models can point out possible necessary conditions for the granule layer inputs, Purkinje cells, and molecular layer interneurons in order to express stochastic synchronization. These conditions can then be explored experimentally.

Competing Interests
The author declare that they have no competing interests.