Blind Deconvolution of Postsynaptic Currents 1=p kskp~ P DsiDp . There are efficient optimization techniques i for deconvolution under certain conditions, for instance when the regularizer is an L1 norm [6] or L2 norm [7]. The L1 norm ensures sparseness; neither enforces positivity of the signal. Although our algorithm works with a variety of priors, for instance a Gaussian or a Burg entropy, we primarily use the so- called quantified maximum entropy (QME) prior [8] defined by the negative of the entropy, the negentropy H: H~ X ½mi{sizsi logðsi=miÞ , i in which si are the elements of s, and mi are the elements of the default signal m, usually set to the same constant value for all times. H has useful regularization properties by enforcing both positivity and sparsity in accordance with our stated prior assumptions about s. QME approximates the true maximum entropy prior while remaining numerically favorable due to the differentiability and convexity of H. In the following paragraph we outline the treatment of the QME prior implemented in the commercial software package MemSys5 [9]. When the filter is kept constant, the posterior is given by Pðs,fDhÞ. The signal is sparse, meaning that s is zero over most of its length. At any time when s is non-zero an event may have occurred; these abrupt spike-like events are spread along the length of s (Fig. 1C). In contrast, f is dense, much shorter than s, and mostly smooth (Fig. 1B). Lastly, we know that our signal must be everywhere positive and the impulse response must be everywhere negative. This problem is ill-posed because the number of degrees of freedom in the variables we wish to infer, s and f, exceeds that of the data x. There is an infinity of ways of maximizing the likelihood by fitting the noise, so a regularizer is needed to guide the algorithm toward the type of sparse solution that we require. A specific aim of our algorithm is to analyze electrophysiological records with complex EPSCs. We therefore developed a method that takes EPSC kinetics into account to accurately disentangle most complex EPSCs. Fitting complex EPSCs by multiple events cannot be achieved by simply extending the process used to fit a single event: the nature of the task changes dramatically when events are allowed to overlap. For instance, although fitting a single Gaussian to data is trivial, fitting a sum of two or more Gaussians is challenging and is an active area of research in the field of clustering [2]. A special case of deconvolution occurs when the likelihood is separable in the data, all events being separated by a time greater than the length of the filter. The equations then simplify substantially and each event can be fitted individually (Appendix S1). But unless the data are separable in the sense described here–that is, containing no complex EPSCs–an algo- rithm must overcome the challenge of complex EPSCs by optimizing the signal globally rather than on an event-by-event basis. A consequence of our sparseness prior is a non-linear filtering process, described below as the deconvolution step, which we use to separate s from f. Linear deconvolution filters, such as the Wiener filter, do not break the symmetry between s and f, and hence algorithms that use linear filters cannot be blind but require f to be fully specified a priori [3–5]. With the prior and likelihood defined, the posterior must be characterized. A point estimate, in the form of the maximum or mean of the posterior, is a common and convenient way to convey information about the posterior, and is also practical from the point of view of an experimentalist seeking an optimal way to visualize the signal. However, it is difficult to solve for signal and filter simultaneously. Although Monte Carlo techniques are applicable to the problem, we opted to make use of approximation techniques that are nearly as accurate but considerably faster. With this approach, the run-time of the algorithm scales approximately linearly with the length of the data. We decompose the problem into two parts. During the first step, deconvolution, we find the optimal signal given the data and a preliminary estimate of the filter. During the second step, filter estimation, we improve our estimate of the filter given the data and our current estimate of the signal. We alternate between these two steps in a way analogous to expectation maximization. Deconvolution Starting with the data x and a filter f, we find a sparse representation by optimizing the strength of a maximum- entropy regularizer. A regularizer is a cost function introduced into the optimization procedure to ensure the smoothness, magnitude, or some other property of the representation being sought, and often corresponds from the Bayesian perspective to the logarithm of a prior. A common family of regularizers seeks to minimize the Lp norm of the signal vector s: PLoS ONE | www.plosone.org 3 June 2012 | Volume 7 | Issue 6 | e38198 3 ð Þ 4 ð Þ PðsDx,f,a,hÞ~ exp {1x2{aH 2 Z , 1 in which x2~s2kx{s fk2 is the sum of squared errors. We have introduced the regularization strength a, which determines the influence of the regularizer. The normalization constant or evidence, Z~ Ð ds exp {1x2{aH , which depends on f, a 2 and additional parameters h, ensures that the posterior is a normalized probability density. Finding the best value of the regularization parameter is frequently a challenge in deconvolu- tion problems, as sparsity comes about only at the correct value of a. When a is too large, the prior dominates and the inferred signal lacks the full resolution available from the data. When a is too small, x corresponds to fitting the noise. Because of our Bayesian formulation, we can use the evidence framework [10] to optimize a. Instead of fully marginalizing over a, we maximize the evidence, Z, with respect to a, because fluctuations in a around its optimal LZ value are very small. At the maximum, ~0 defines the optimal La regularization a^. Due to the form of x2, both the noise strength s2 and the regularization constant a may be optimized at the same time, drawing attention to the fundamental importance of the signal-to-noise ratio in determining the resolution. The point estimate of the signal is found by maximizing the posterior probability PðsDx,f,a^Þ with respect to s, for ^s~ arg max PðsDx,f,a^Þ. The presence of noise and the approx- s imate nature of the prior means that ^s is nearly, but not exactly, zero in regions where there are no events but there is noise. Any point estimate suffers from the presence of spurious tiny events in ^s; if we were to inspect the full posterior, we would find that their location and presence is highly uncertain. The process of obtaining a sparse representation by thresholding ^s reduces the degrees of freedom significantly, typically by three orders of magnitude on our experimental data. This sparse representation is the key contribution of this first step, the