Inclusion of noise in iterated firing time maps based on the phase response curve

Fred H. Sieling,* Carmen C. Canavier, and Astrid A. Prinz Wallace H. Coulter Department of Biomedical Engineering, Georgia Institute of Technology and Emory University, Atlanta, Georgia 30332, USA Neuroscience Center for Excellence and Department of Ophthalmology, Louisiana State University Health Sciences Center, New Orleans, Louisiana 70112, USA Department of Biology, Emory University, Atlanta, Georgia 30322, USA Received 16 December 2009; revised manuscript received 10 February 2010; published 25 June 2010

The problem of synchronization and patterning of synchronization of oscillators is a problem that is of general interest in physics ͓1-6͔.More specifically, the study of pulse-coupled oscillators has commanded substantial interest ͓7-14͔.Networks of pulse-coupled oscillators are used to study many things ͓15͔-e.g., plate tectonics ͓16͔, heart rhythms ͓17͔, and neural networks ͓18͔.Synchrony is a dynamical feature of such networks ͓19-21͔.These systems consist of oscillators coupled by discrete pulses.Phase response theory is often used to study network synchrony.The phase response curve ͑PRC͒ is tabulated as the shift in phase of an oscillator in response to a stereotyped input pulse ͓22͔.When the pulse is sufficiently weak, interactions can be added linearly so the infinitesimal PRC ͑iPRC͒, or response to weak perturbations, which is the foundation for a powerful technique-the theory of weakly coupled oscillators-can be used to predict network activity ͓23͔, however, many real systems are coupled too strongly for weak coupling approaches ͓24,25͔.
Phase response theory has a rich history in biological oscillators-especially in cardiac electrophysiology and neurophysiology.In the latter field, phenomena emerging from networks of two neurons are often used to help understand behaviors of networks containing large populations of neurons ͓18͔, however, this paradigm is limiting-the complexity seen in biology is far beyond that explained by existing models, which are largely limited to networks of identical neurons coupled identically and weakly with no noise included.Conversely, real neural systems often contain heterogeneous neural types coupled heterogeneously, strongly, and with noise.Several recent studies have pushed one or more these boundaries, e.g., showing noise-stabilized antiphase oscillations ͓26͔, synchronization in a network despite noise ͓27͔, optimal PRC shape for noise-driven synchrony ͑stochastic synchrony͒ ͓28͔, and bistability induced by shared noise ͓29͔.Although noise has been included in maps of pulse-coupled neurons ͓30͔, no study has yet formalized the effects of noise in nonweakly coupled pulse-coupled networks.In ͓31͔ we briefly proposed a method to address these issues simultaneously, and we applied it in one example as a proof of concept.
Here, we present our method formally and apply it generally to describe 78 experimental networks from ͓31͔.We explicitly incorporate noise measured during experimental determination of the PRC into the prediction method.In a hybrid network of two bursting neural oscillators coupled via synaptic excitation ͑ syn Ϸ 10 ms with a network period of about a second͒, we found that this method could accurately predict phase locking in the presence of 10% variability in the intrinsic period.We have shown that networks including a biological oscillator may undergo noise-induced transitions, and that including measured noise in the firing time map may capture the effect of noise, giving results that are qualitatively and quantitatively similar to experiment ͓31͔.
Following the biological tradition ͓32,33͔, in a two neuron network our phase response curves ͑PRCs, F qi below͒ describe the effect of a stereotyped input from the partner neuron j on neuron i as a function of the timing of input, or phase i .The qth order PRC is tabulated experimentally by measuring the qth cycle lengths P qi and P 0i , where first order resetting q = 1 measures response in the perturbed cycle, and second order resetting q = 2 measures response in the next subsequent cycle.As such, F 1i measures the change of the cycle containing the start of input from neuron j If the input is an infinitesimal voltage deflection, then our F 1i is equivalent to PRCs used in weak coupling theory ͓34,35͔.Second order resetting F 2i has been previously shown to be important in bursting neurons ͓24,25,31͔ and was recently modeled as leftover charge in dendrites resulting in partial reset ͓36͔.
We use these PRCs to construct a general model of pulsecoupled phase oscillators with no delays * fred.sieling@bme.gatech.edu˙i = where i ͓0,1͒ is the phase of neuron i modulo 1.A threshold event such as a spike or burst occurs when the phase reaches one, and the phase is immediately reset to zero afterwards.The coupling function is given by F 1i , ␦ is the Dirac delta function, and t j ‫ء‬ is the time of a spike or burst initiation in neuron j.Because we use ideas from phase resetting theory, we constrain the behavior of our uncoupled oscillators to autonomous periodic activity, i.e., spiking or bursting, where for each neuron i there is an intrinsic period P 0i such that i ͑t + P 0i ͒ = i ͑t͒.Note that the oscillators may be heterogeneous, each having its own intrinsic period P 0i P 0j .
The major drawback to using PRCs as coupling functions in our model is that we must assume the system returns to its limit cycle between pulses.In the phase space this means that when a relatively strong input moves the coupled oscillator far from its uncoupled trajectory, it must return before the next input arrives.In general, this is not a problem for weakly coupled oscillators.Returning to the limit cycle between inputs also implies that in the case of a neuron receiving multiple inputs in between two firing events, the second order effects F 2i must be zero for all inputs except the last one received before the second event.For two neurons, this is rarely a problem.
To expand Eq. ͑1͒ for nonweak inputs, we include F 2i by adding a second term, where i last is the value of i at the last t j ‫ء‬ if an input was received by neuron i in the previous cycle.
The switch A i is one if an input was received in the previous cycle and zero otherwise.Additionally, each i must now be defined for small negative values because if F 2i is a delay, i will be reset to a negative value.A geometrical interpretation of the negative branch in i is given in ͓37͔-it "corresponds to an isochron that curls around the limit cycle, intersecting it at a position that is retrograde to the peak of the action potential."For i Ͻ 0, F qi are undefined, however, this is not a problem in our analysis, where we assume that for i Ͻ 0, F qi ͑ i ͒ = F qi ͑0͒ because few inputs ͑if any͒ were received on this branch.
To account for noise in our experimentally measured PRC data, we simply include it directly in the PRCs, where F qi fit is a piecewise polynomial fit to the experimental data, qi represents unitary Gaussian noise, and qi is the standard deviation of the measured noise as a function of i , found by binning the data and measuring in each bin.All PRCs used in this study are the same as in ͓31͔.
Consider a network of two heterogeneous coupled neural oscillators as in Eq. ͑2͒.To find modes of rhythmic network activity, we modify the iterative pulse coupled map with no predetermined firing order from ͓38͔ to include noise in the PRCs as in Eq. ͑3͒.We iterate the map as follows, where i ͓m͔ is the phase of neuron i after all effects of the previous pulse event m are added, S͓m͔ is the set of neurons that emit a pulse m, and R͓m͔ is the set that receives a pulse.
We initialize each neuron with a phase i ͓0͔ and a memory of second order resetting F 2i ͑ i last ͓0͔͒.Then, we poll for which neuron͑s͒ will fire next: if P 0i ͑1− i ͓m͔͒ = min͕P 0i ͑1− i ͓m͔͖͒ over all i, then i S͓m͔, else i S͓m͔.If i sends a pulse then j receives one: for i S͓m͔, j R͓m͔ else j R͓m͔, where i j.Now update each neuron according to whether it fires a pulse and whether it receives a pulse: If i S͓m͔, i R͓m͔, then reset i to zero plus F 2i from previous cycle: i ͓m +1͔ =0−A i F 2i ͑ i last ͓m͔͒, and reset the switch A i to zero.If i S͓m͔, i R͓m͔, then set and reset the switch A i to zero.If i S͓m͔, i R͓m͔, then advance phase by elapsed time plus 1st order resetting from current cycle: Update the memory of second order resetting: for i R͓m͔, i last ͑m +1͒ = i ͑m͒ + ͓t͑m +1͒ − t͑m͔͒ / P 0i , and set the switch A i to one.Now start again-poll for which neuron͑s͒ will fire next, etc.
To summarize: This map is difficult to analyze because of the discontinuities introduced by the min function.First, there is no assumption of firing order.It is even possible that the neurons will fire synchronously or that one neuron will not fire at all.Next, a neuron's phase is updated differently depending on whether it fires a pulse, receives a pulse, or does both simultaneously.Note that the phase is updated each time that an input is received, in contrast to the weak coupling assumption where phase is not updated when each pulse is received based on resetting that occurs earlier in the cycle.
We have ignored conduction delays here, although they may be added to the map ͓39͔.In our experimental system consisting of bursting neurons, we defined an event as the start of a burst in the presynaptic neuron.Using the map, we can iterate a number of cycles N map roughly matching the number of observed cycles N obs and compare results using circular statistics as in ͓31͔.The inclusion of noise in the firing time map is similar to ͓30͔; the major difference is the use of F 2 in addition to F 1 in our method.Further, we explicitly restrict the tails of i according to their causal limits.
In ͓31͔, these limits translate to recovery interval t r Ն 0 and stimulus interval t s Ն 0 as defined below.

I. COMPARISON OF MAP TO EXPERIMENTAL DATA AND NOISE-FREE METHOD
To evaluate our methods on experimental networks, we applied them to an existing data set of 78 coupled networks ͓31͔.These networks are constructed from one biological neuron and one model neuron by using the dynamic clamp ͓40,41͔ to implement artificial synapses.The biological neuron used ͑N =5͒ was the AB/PD complex from the pyloric network in the Stomatogastric ganglion of Homarus Americanus, which is a regularly bursting group of 3 electrically coupled neurons that act as a single oscillator.This group was isolated pharmacologically using standard methods ͓31͔.The model neurons used were from a heterogeneous set of four regularly bursting pyloric model neurons from a database ͓42͔.
We compared experimental observations to predictions made using a noise-free method of pulsatile coupling ͓25,31͔ and to the predictions of our iterated map using steady state stimulus intervals t s for each neuron in the network, where t s is the interval between the time of an event in one neuron and the time of the next input it receives from its partner neuron ͑Fig.1͒.For the case of 1:1 locking this is equivalent to t͓M͔ above, where M is the set of odd or even whole numbers.Using both methods, we predicted qualitative behavior-whether a mode of 1:1 phase-locking would be observed-and quantitative behavior-the stable phase difference and strength of phase-locking.We quantify the strength of phase-locking using the circular statistic R 2 .When each n of N phase differences ⌽͓n͔ =2 ‫ء‬ t s1 ͓n͔ / P net ͓0,1͒, P net = 1 N ͚ n=1 N ͑t s1 ͓n͔ + t s2 ͓n͔͒ is plotted on the unit circle, then R ͓0,1͔ is the length of the vector from the origin to the average position of the plotted points ͓43,44͔.We used a threshold R 2 Ͼ 0.7 to define phase-locking as in ͓31͔.
In the 78 networks tested, the noise-free method was qualitatively incorrect 16 times and the iterative map was able to make the correct prediction for 8 of those networks ͑see ͓31͔ for full analysis of the noise-free method in these networks͒ as well as 61 that were correctly predicted by the noise-free methods, for a total success rate of 69 out of 78, an improvement of 10%.We excluded eight networks that exhibited continuous spiking due to positive excitatory feed-back as a result of the coupling, breaking the assumption of a limit cycle burster under which the PRCs were generated ͓31͔.Five failures of the map were due to overestimation of noise.Over the course of hours, noise levels in the biological neuron varied.In all cases, removing noise from the firing time map by setting = 0 in Eq. ͑3͒ produced a match to the predictions of the noise-free method.
To best compare predictions made using our firing time map to those made using the noise-free methods, we focus on networks that were experimentally observed to phaselock.In Fig. 2, we compare methods for 28 such networks by showing steady state t s intervals for each neuron in the network.Each network is represented by an ordered group of symbols separated by vertical dotted lines, showing t s values obtained for each of the methods used.Where a phase-locked mode was not observed or predicted by a method, the appropriate symbol is plotted in red below.
When F 2 was removed from the firing time map ͑not shown͒, F 2i = 0 in Eq. ͑2͒, prediction accuracy dropped from 85% to 45%, as we would expect from ͓24,25,31͔.In B, we compare prediction methods in 8 networks previously shown to exhibit noise-induced transitions ͓31͔.The analytical method, which is noise free, was qualitatively incorrect in each of these cases, and the map was correct in each case.Out of the 8 networks with such behavior, only 1 was observed to phase-lock and the map method was qualitatively accurate ͑network 8, shaded͒.Note that in this network, the analytical methods did not predict a phase-locked mode, so adding noise to the PRCs caused a new stable fixed point to appear, showing that in the firing time map, noise can both generate and destroy modes of phase locking.
We have verified that variability in a PRC does modulate the strength of phase locking in a network of 2 neurons.When noise was added to the map, the differential effects of noise on structurally stable and unstable systems was captured qualitatively and quantitatively using circular statistics.This simple idea is a powerful tool for understanding the predictive utility of PRCs because it simulates the measured variability of the biological neuron and gives outputs that are readily compared to experimental data.
Bendels and Leibold ͓27͔ note a threshold in noise level above which "low-synchrony" phase-locking occurs but is not predicted by their mean-field approximation method due to occasional slipping.This threshold can be explained by the noise level at which the envelope of phase-dependent noise around the PRC disrupts the basin of attraction of the noise-free fixed point.We do not analyze this threshold here, but our data give examples where a stable fixed point in the noise-free system transitions to a quasistable fixed point when noise is added, seen as stable phase-locking with occasional slipping.
An input can never advance the event time to a time before the input was applied: this is the causal limit.If no input is given at phase , then the next event is expected to follow at a time interval of P 0 ͑1−͒, the remaining interval.Therefore the observed phase resetting can never drop below the line ͓F 1i ͑͒ = −1͔, which biases the noise observed in first order PRCs.For example, in Tsubo et al. ͓45͔ the authors note that they could fit the distributions of the fluctuations in the phase responses at phases less than 0.5 with Gaussian distributions with a variance that was nearly constant.However, at phases greater than 0.5, the variance decreased.Figure 2͑E͒ of that paper makes it clear that the causal limit is responsible, but the role of the causal limit was not addressed.Phoka et al. show that noise in Purkinje cells introduces a bias in the mean measured PRC due to the causal limit ͓46͔.To correct for the bias, they use all spikes in a train as reference, one at a time, rather than using only the spike before and after.Polhamus et al. proposed a different solution, which was to average only ISIs longer than the stimulus interval to determine the average period to use to calculate the phase.Our method accounts for the missing variance in the first order resetting by treating second order resetting separately.
Netoff et al. ͓30͔ also observed PRCs in which the noise was limited by causality.In their studies of excitatory inputs, the scatter decreased inversely proportional to the square root of the phase, but for inhibitory inputs causality is not an issue and variance was independent of phase.Netoff et al. ͓30͔ also incorporated noise into a PRC-based map.They chose firing times for each cell by adding the spike time response curve ͑STRC͒ value ͑mean+ normally distributed random component scaled by the time-dependent standard deviation, each measured from the STRC͒ to the unperturbed firing period of the cell in question ͓30͔.
Since our map is different from that of Netoff et al. ͓30͔ only in that we used second order resetting, it is interesting that our maps perform more poorly if F 2 is ignored.This is likely because the neurons in ͓30͔ are spiking, rather than the bursting neurons that we used, and spiking neurons tend to have small F 2 amplitudes.The new method incorporates the true variance rather than underestimating it at late phases.To our knowledge, this point has not been made previously.For example, after commenting on the phase dependence of the variance observed by Netoff et al., Ermentrout and Sanders state that "An adequate theory of this dependence remains to be resolved" ͓26͔.In our method, the variance in F 2 in the subsequent cycle compensates for reduction in variance at late phases in F 1 .
We acknowledge Robert Butera, Andrey Olifer, Srisairam Achuthan, and Lakshmi Chandrasekaran for their helpful comments.This work was supported by NIH ͑Grant No. 2R01NS054281-5͒.Comparison of observed behavior from ͓31͔ to analytical and firing time map predictions including F 2 in the firing time map.Empirical observation ͑o and error bars͒, analytical prediction ͑x͒, and prediction using the firing time map ͑square and error bars͒ are shown for mdl ͑model neuron, blue͒ and bio ͑biological neuron, green͒ neurons in each hybrid network, separated by dotted lines and sorted by qualitative match for clarity.Modes that were not observed or predicted to exist are marked by a red symbol ͑o, x, or square͒ in the row labeled DNE.For clarity, error bars are shown only in one direction.Analytical predictions are noise-free so they do not have error bars.The x axis is the arbitrary experimental network number.A. For networks observed to phase lock, the noisy map produced results that were generally in agreement with observations.In one case ͑red x͒ the noisy map succeeded where the noise free method failed, but in five cases ͑red squares͒ the noisy map failed presumably due to overestimation of noise.͑b͒ In cases where noise was previously shown to cause the analytical prediction to fail, the noisy map correctly predicted the failure to lock in seven cases, and the noise-induced ability to lock, albeit with some error, in one case.

FIG. 1 .
FIG.1.Definition of terms, A: stimulus time ͑t s ͒ and response time ͑t r ͒ during 1:1 phase locking are shown for bursting neurons 1 and 2. The shaded regions correspond to the burst duration, which can be a substantial fraction of the cycle period.
FIG.2.͑Color͒ Map with noise.Comparison of observed behavior from ͓31͔ to analytical and firing time map predictions including F 2 in the firing time map.Empirical observation ͑o and error bars͒, analytical prediction ͑x͒, and prediction using the firing time map ͑square and error bars͒ are shown for mdl ͑model neuron, blue͒ and bio ͑biological neuron, green͒ neurons in each hybrid network, separated by dotted lines and sorted by qualitative match for clarity.Modes that were not observed or predicted to exist are marked by a red symbol ͑o, x, or square͒ in the row labeled DNE.For clarity, error bars are shown only in one direction.Analytical predictions are noise-free so they do not have error bars.The x axis is the arbitrary experimental network number.A. For networks observed to phase lock, the noisy map produced results that were generally in agreement with observations.In one case ͑red x͒ the noisy map succeeded where the noise free method failed, but in five cases ͑red squares͒ the noisy map failed presumably due to overestimation of noise.͑b͒ In cases where noise was previously shown to cause the analytical prediction to fail, the noisy map correctly predicted the failure to lock in seven cases, and the noise-induced ability to lock, albeit with some error, in one case.