maximum entropy deconvolution, and allows us to proceed to refine our estimate of f. Estimation of the Impulse Response The quality of the estimate ^s computed in the previous step depends crucially on the estimate of f we supplied. In this section, we show how to recursively improve the estimate of f over the initial guess. We have formulated the problem such that the impulse response–the fundamental EPSC shape–is invariant, although the amplitude and location in time of each event differs. Conditional on ^s, we estimate the impulse response non- parametrically by ^f~ arg max PðfDx,^sÞ, in close analogy to ^s. f Because f has finite support, its degrees of freedom are far fewer than provided by x and ^s. In this well-posed situation, if we consider the likelihood on its own, ^f is the solution to the Wiener- Hopf equations [7]: 1 {s22+f x2~Cxs{f Css~0, 5 ð Þ in which Cxs and Css are cross-correlation functions defined by Cxy(t)~ P x(tzt’)y(t’) and * is the convolution operator. t’ We may simply use the sparse ^s we inferred in the previous step and optimize f, taking care to enforce the negativity of f as we do so. When the events are all well-separated, this process is similar to identifying the isolated events and finding the impulse response that best fits them in a self-consistent manner. Convergence of Iteration Successive applications of the deconvolution and impulse- response estimation steps quickly converge to the optimal signal and impulse response. Our computational approach is described in Appendix S1, with the nested nature of the resulting iterative loops shown in Figure S2. Despite the availability of a sparse representation following the second step above, the overlap among events implies that directly optimizing the amplitudes and times of events could converge to the wrong result. In general, it may also be prudent to run the algorithm from a variety of starting conditions in order to verify that the phenomenon under study is not strongly dependent on artifacts that have been introduced during deconvolution. When we explored different initial condi- tions with our data sets, we found that the algorithm always converged to the same solutions. Unlike typical software packages, our method is defined not by an algorithm but by the generative model in the likelihood and the priors over the variables of interest. We can take advantage of additional capabilities by analyzing the posterior; we can, for instance, generate the error of an estimate. As better approxima- tion techniques improve the quality of the solutions, the limit set by the signal-to-noise ratio may be approached to arbitrary proximity. Although many aspects of EPSC deconvolution are analogous to image deconvolution, there are several features that are unique to electrophysiology. The additive noise in experimental record- ings is often non-white and records can exhibit baseline drift over time. Furthermore, the impulse response is always negative. Because the AMPA receptors that generate EPSCs are fluctuat- ing–unlike most point-spread functions–the impulse response is nondeterministic. The impulse response may vary owing to the presence of more than one synapse. Parameters of the model whose values we assume to be constant for the duration of a Blind Deconvolution of Postsynaptic Currents record, such as the power spectrum of the additive noise and the membrane time constant involved in the recording, may in fact exhibit non-stationarity. Incorporating the effects of AMPA receptors, multiple synapses, and non-stationarity would be a useful goal of future research. Performance on Physiological Data To test the performance of our algorithm, we analyzed data recorded from the bullfrog’s amphibian papilla, an auditory end- organ whose spontaneous exocytosis reveals complex, highly overlapping EPSCs [1]. The experiments were conducted under protocols approved by the Institutional Animal Care and Use Committee of The Rockefeller University. In a first example, we sought to demonstrate the ability of our algorithm to identify individual EPSCs from complex records. From long, stable records (Fig. 2A) the algorithm readily recovered the characteristic sum-of- exponentials shapes of the constituent EPSCs (Fig. 2B). The distributions of EPSC amplitudes (Fig. 2C) were compatible with distributions found elsewhere [11]. The distributions of inter- EPSC time intervals were well approximated by exponential functions (Fig. 2D). We also characterized the data by examining correlations between event amplitudes and timing. Little or no correlation was apparent between amplitudes and interevent intervals in our recordings (Fig. 2E–G). In a second example, we sought to use our algorithm to answer a physiologically important question: are stimulus-dependent increases in the postsynaptic current the consequence of an increase in the rate of EPSCs or an increase in the mean amplitude of EPSCs? Without deconvolving the signal, it is not possible to answer this question. We explored the complex waveforms resulting from an experiment in which the presynaptic cell was given a series of depolarizing steps between 265 mV and 220 mV. This procedure was a rigorous test of our algorithm, for large depolarizations produce postsynaptic currents that are highly complex and overlapping. Using our algorithm, we conclude that the increase in total postsynaptic current stems primarily from an increase in the rate of EPSCs (Fig. 3A) rather than an increase in their amplitudes (Fig. 3B). Performance on Synthetic Data We also estimated the performance of our algorithm by analyzing synthetic data whose EPSC amplitudes, timings, shapes, and noise spectrum were fully controlled. We tested our algorithm’s ability to extract from a complex record the characteristic EPSC impulse response as well as the timing and amplitude of individual events. We also informally compared the performance of our algorithm to that of MiniAnalysis [12], a popular software program for the analysis of physiological recordings. MiniAnalysis represents a heuristic approach that provides flexibility at the cost of requiring operator assistance. By contrast, our global deconvolution approach is self-consistent and allows the program to run with less intervention. Synthetic data were constructed by simulating AMPA receptors to generate EPSCs: at an instant drawn from a Poisson process, a fraction 0.1–0.5 of 150 receptors was selected to be open instantaneously by drawing from a binomial distribution. The receptors were allowed to reclose stochastically with a timescale of about 1 ms. The result was low-pass filtered with a time constant of approximately 0.25 ms and sampled at 50 ms intervals, parameters selected to match the membrane time constants and sampling frequencies of typical recordings, such as those in Figure 2. Additive Gaussian noise was filtered with a low-pass filter fitted to the shape of a typical noisy experiment and adjusted to a root-mean-square amplitude of PLoS ONE | www.plosone.org 4 June 2012 | Volume 7 | Issue 6 | e38198