Biophysical Journal 85:2170-2185 (2003)
© 2003 The Biophysical Society
Wavelet Analysis of Nonstationary Fluctuations of Monte Carlo-Simulated Excitatory Postsynaptic Currents
F. Aristizabal * and
M. I. Glavinovic 
Departments of * Chemical Engineering and Physiology, McGill University, Montreal, Canada
Correspondence: Address reprint requests to Dr. M. I. Glavinovic, Dept. of Physiology, McGill University, 3655 Sir William Osler Promenade, Montreal, Quebec H3G 1Y6, Canada. Tel.: 514-398-6002; Fax: 514-398-7452; E-mail: mladen.glavinovic{at}mcgill.ca.
 |
ABSTRACT
|
|---|
Tracking spectral changes of rapidly varying signals is a demanding task. In this study, we explore on Monte Carlo-simulated glutamate-activated AMPA patch and synaptic currents whether a wavelet analysis offers such a possibility. Unlike Fourier methods that determine only the frequency content of a signal, the wavelet analysis determines both the frequency and the time. This is owing to the nature of the basis functions, which are infinite for Fourier transforms (sines and cosines are infinite), but are finite for wavelet analysis (wavelets are localized waves). In agreement with previous reports, the frequency of the stationary patch current fluctuations is higher for larger currents, whereas the mean-variance plots are parabolic. The spectra of the current fluctuations and mean-variance plots are close to the theoretically predicted values. The median frequency of the synaptic and nonstationary patch currents is, however, time dependent, though at the peak of synaptic currents, the median frequency is insensitive to the number of glutamate molecules released. Such time dependence demonstrates that the "composite spectra" of the current fluctuations gathered over the whole duration of synaptic currents cannot be used to assess the mean open time or effective mean open time of AMPA channels. The current (patch or synaptic) versus median frequency plots show hysteresis. The median frequency is thus not a simple reflection of the overall receptor saturation levels and is greater during the rise phase for the same saturation level. The hysteresis is due to the higher occupancy of the doubly bound state during the rise phase and not due to the spatial spread of the saturation disk, which remains remarkably constant. Albeit time dependent, the variance of the synaptic and nonstationary patch currents can be accurately determined. Nevertheless the evaluation of the number of AMPA channels and their single current from the mean-variance plots of patch or synaptic currents is not highly accurate owing to the varying number of the activatable AMPA channels caused by desensitization. The spatial nonuniformity of open, bound, and desensitized AMPA channels, and the time dependence and spatial nonuniformity of the glutamate concentration in the synaptic cleft, further reduce the accuracy of estimates of the number of AMPA channels from synaptic currents. In conclusion, wavelet analysis of nonstationary fluctuations of patch and synaptic currents expands our ability to determine accurately the variance and frequency of current fluctuations, demonstrates the limits of applicability of techniques currently used to evaluate the single channel current and number of AMPA channels, and offers new insights into the mechanisms involved in the generation of unitary quantal events at excitatory central synapses.
|
INTRODUCTION
|
|---|
The analysis of current fluctuations or noise analysis has long been an important tool for the study of ion channels (Derksen and Verveen, 1966
; Verveen and Derksen, 1968; Katz and Miledi, 1970, 1972; Anderson and Stevens, 1973). Although single-channel recording has replaced noise analysis as the prevalent method of evaluating single channel properties, noise analysis is still needed when: a), the single channel conductance is low, b), channels are inaccessible, or c), rapid kinetics of activation and/or desensitization complicate single-channel analysis. Such is the case with glutamate activated fast synaptic currents in the central nervous system. Current fluctuation analysis was adapted for nonstationary conditions by Sigworth (1981). Nonstationary mean-variance current analysis, now used in several laboratories (Robinson et al., 1991; Traynelis et al., 1993; De Koninck and Mody, 1994; Barbara and Takeda, 1996; Nusser et al., 1997, 2001; Traynelis and Jaramillo, 1998; Wierenga and Wadman, 1999), provides estimates of the number of receptors and their saturation, but the interpretation is not simple. The spatial inhomogeneity of transmitter concentration in the synaptic cleft alters mean-variance plots (Kruk et al., 1997). Even when spatially uniform, the instantaneous concentration is highly variable, especially at low concentrations (Glavinovic, 1999). Finally the number of receptors present postsynaptically may also vary, especially if the synaptic morphology is highly heterogeneous (Trommershauser et al., 2001). The degree of postsynaptic receptor saturation during fast miniature excitatory postsynaptic currents (mEPSCs) could be evaluated more accurately if the spectral changes of current fluctuations could be determined. The spectra are expected to shift to higher frequencies at high saturation levels (Silberberg and Magleby, 1993). The power spectra of the current fluctuations have been used in the past to study the ion channel kinetics, but until now only as composite spectra-spectra of the synaptic current fluctuations that ignore the time-dependent changes (Conti et al., 1980; Sigworth, 1981; Robinson et al., 1991; Otis and Mody, 1992; Traynelis et al., 1993; De Koninck and Mody, 1994; Traynelis and Jaramillo, 1998).
It is not an easy task to evaluate the spectral changes of synaptic currents. The spectral estimates provided by Fourier analysis are accurate only for current segments much longer than the mean channel open time or effective mean open time (Silberberg and Magleby, 1993). Traditional Fourier analysis is thus poorly suited for nonstationary spectral analysis. It has only one set of bases for signal decomposition (sines and cosines) that cannot be tailored to the nature of the signal analyzed, and it localizes only in frequency, but not in time. With small sections of the signal analyzed at a time ("Short-Time Fourier Transform" (Gabor, 1946)), the spectral resolution is the same for all frequencies (frequency-independent windowing), resulting in poor accuracy of spectral estimates at low frequencies.
The recently developed wavelet technique offers several advantages: a), an infinite number of bases for signal decomposition, b), localization of spectral changes in frequency and in time, and c), frequency-dependent windowing of the signals (Daubechies, 1991, 1998; Coifman and Wickerhauser, 1992; Mallat, 1989, 1992). The technique has been used recently to study oscillatory brain activity (Senkowski and Herrmann, 2002). Wavelet analysis, and especially its generalized version—wavelet packet analysis—can successfully track spectral changes of Monte Carlo-simulated fast mEPSCs (Wahl et al., 1996; Kruk et al., 1997; Glavinovic and Rabie, 1998; Glavinovic, 1999), and thus promises to be an exciting new tool for the analysis of receptor saturation and nonstationary synaptic current fluctuations. A preliminary account has appeared as an abstract (Glavinovic and Aristizabal, 2000).
|
METHODS
|
|---|
Wavelet fundamentals
Since the wavelet decomposition of signals is a comparatively recent technique (Teolis, 1998), we will present a brief theoretical introduction of its fundamentals, compare it with the traditional Fourier transform technique, and discuss its advantages and disadvantages.
According to Fourier analysis, the signal f(t) can be decomposed into a sum (infinite) of sine and cosine functions given by the Fourier transform. In the complex domain, the equations for the decomposition are:
 | (1) |
which is the sum over all time of the signal f(t) multiplied by a complex exponential, which breaks down into real and imaginary sinusoidal components.
The wavelet transform of the signal f(t) is defined as:
 | (2) |
where * denotes complex conjugation, a represents the scale parameter, and b the translation parameter (time shifting), whereas 
a,b(t) is calculated by scaling the basic wavelet at time b and scale a
 | (3) |
When a becomes large, a,b becomes a stretched version of the basic wavelet (Fig. 3), which is a useful property when analyzing low-frequency components of the signal. In contrast, when the scale parameter is small, the basis function a,b becomes contracted, as needed for the analysis of the high-frequency components of the signal. Wavelet analysis can thus be viewed as an analysis consisting of band-pass filters having a constant relative bandwidth. Given a signal of n data points and with j levels of decomposition, a is 2j; b, which is a multiple of an integer and a, determines the number of translations and thus the time resolution.
View larger version (25K):
[in this window]
[in a new window]
|
FIGURE 3 A comparison of a Fourier and a wavelet transform. Original signals, top traces. (A) Eight terms of the Fourier decomposition series having the highest magnitude spectral coefficients are shown below the signal trace (middle traces). Numbers on the left are Fourier coefficients in polar coordinates (magnitude and phase). The sum of the eight terms (bottom traces) gives an approximation (thick line) of the original signal (thin line). (B and C) The same signal (top trace) but decomposed using a wavelet method. Eight terms of the series were chosen on the basis of the amplitude of their real value coefficients (given on the left of each individual decomposition term) but are presented in an inverse frequency order (i.e., with the lowest frequencies at the top). Coefficients, which are localized in time, are given on the right of each decomposition term. Filled symbol shows the coefficient used for the calculation of the corresponding decomposition term. The level, position, and the indication whether the decomposition term belongs to a detail or approximation are indicated on the left of each set of coefficients (see text). Note that using the Haar instead of Daubechies-3 (db3) wavelet leads to a better reconstruction because of the greater similarity between the signal being analyzed and the wavelet (see text).
|
|
Signal decomposition
The analysis begins with a decomposition of the signal into an approximation and a detail (Fig. 3). The approximations are computed by convolving the signal with the wavelet scaling function ("father wavelet"). Convolving the signal with the wavelet filter ("mother wavelet") yields the details. The decomposition process is iterative, but in each successive step it is only the approximations that are decomposed. At the end of such a procedure the signal is broken down into many lower-resolution components. In wavelet packets (Coifman and Meyer, 1990; Coifman et al., 1992), which are a generalization of the wavelet transform, both the approximations and details are decomposed at each step. Wavelet packet (WP) analysis is a richer method for signal processing, and is used in this study to evaluate current fluctuations. WP analysis is especially useful when better frequency localization is needed, because of its ability to partition not only the low frequency of the signal but also the high frequency.
Fig. 1, A and B (top), show the time course of a scaling function and a wavelet function together with their Fourier transforms at several resolution levels (as indicated) illustrating their filtering properties. As can be clearly seen from the power spectra of the scaling and wavelet functions, the scaling function acts as a low-pass filter, and the wavelet function as a high-pass filter (for the first level of decomposition) and a band-pass filter for all other levels of decomposition (Fig. 1, C and D).
View larger version (31K):
[in this window]
[in a new window]
|
FIGURE 1 Wavelet function is both temporally and frequency localized. (A) Scaling function ("father wavelet"), which produces approximations of wavelet decomposition, acts as a low-pass filter. (B) Wavelet function ("mother wavelet"), which produces the details of wavelet decomposition, acts as a high-pass or (band-pass) filter. Daubechies-3 (db3) wavelet is chosen as an example. Four levels of decomposition are shown (1–4). (C and D) Fourier transform of the scaling function and wavelet function. In the first level, first approximation and first detail are calculated. As the Fourier transform of the wavelet function demonstrates at the first level of decomposition, the wavelet function acts as a high-pass filter, filtering out the highest frequencies, whereas its counterpart scaling function filters out the lowest frequencies. For higher decomposition levels (2–4), the scaling function acts again as a low-pass filter, filtering out progressively lower frequencies. In contrast, for higher decomposition levels the wavelet function acts as a band-pass (and not as a high-pass) filter, but it also filters out progressively lower frequencies. Note that all amplitudes of the Fourier transforms of both wavelet and scaling functions are normalized to 1.0.
|
|
Scale pseudofrequency relationship
Because wavelet or WP analyses do not decompose signals into sines and cosines, a purely periodic signal of frequency Fp (pseudofrequency) is defined to allow wavelet and WP decomposition of the signals to be related to the more traditional Fourier decomposition. This is a convenient way of characterizing the dominant frequency of the wavelet (Yu and Karlsson, 1997; Karlsson et al., 2000). Pseudofrequency corresponding to the scale of wavelet (or wavelet packet) decomposition is calculated as:
 | (4) |
where a is a scale, 
is the sampling period, and Fw is the center frequency of the wavelet in Hz (i.e., the frequency corresponding to the spectral peak of the wavelet). The sampling period and the scale a take into account dilations and contractions of the wavelet function. Fig. 2 shows an example how such a conversion is made with db3 (Daubechies-3) as a wavelet function.
View larger version (28K):
[in this window]
[in a new window]
|
FIGURE 2 Scale to pseudofrequency conversion. (A) Schematic diagram of the wavelet packet decomposition tree. A, approximations; D, details. Three levels of decomposition yield 23 or eight nodes and eight sets of coefficients at the third (and final) decomposition level. (B) Wavelet (db3, thick line) was reconstructed from a single coefficient at the node 2, whose magnitude was set to 1; it was positioned in the middle to avoid the problems occurring at the edges, and corresponds to the set of coefficients with the second lowest frequencies. Single sinusoid (thin line) was reconstructed from the largest Fourier coefficient of the decomposed wavelet. It characterizes well the dominant frequency of the wavelet. (C) Calibration line of the frequency of the dominant sinusoid ("pseudofrequency") versus "node number". The vertical bars indicate the frequency range within which 50% of the energy of the wavelet lies (see text). (D) Wavelet packet spectral densities using single pseudofrequencies (crosses) are very similar to those using all frequencies of the wavelet (open circles; see text).
|
|
We take a zero signal and decompose it using a WP analysis. Out of all coefficients, we arbitrarily take those from the second node third level of decomposition (or DAA3 at the decomposition tree; Fig. 2 A) and set one coefficient to 1 (all others are equal to zero). We reconstructed that coefficient and obtained a wavelet shown in Fig. 2 B. The chosen coefficient is typically located in the middle of the window so that we can have the complete time course of its reconstruct. This reconstruct is then decomposed using Fourier analysis. All coefficients of the Fourier decomposition are then set to zero except the one with the largest magnitude. Its reconstruction is a sinusoid shown in Fig. 2 B. Note that on visual inspection the sinusoid characterizes well the dominant frequency of the wavelet. The same procedure is then repeated at other nodes (not shown). In each case a sequence of decompositions and recompositions yields a dominant or a pseudofrequency. Fig. 2 C shows the relationship between the node number of wavelet packet decompositions and the corresponding dominant frequency. This relationship links the timescale decomposition provided from the wavelet or WP analysis to the time-frequency decomposition available from the traditional Fourier decomposition. The line calculated using the least-squares fit was (±SD): f = -0.057 (±0.007) + 0.068 (±0.001) x N, where f is pseudofrequency in Hz. Note that N corresponds to the frequency ordered nodes. Node ordering is needed because the dominant frequency of the wavelets at individual nodes does not increase monotonically when moving from left to right. The scatter of the individual pseudofrequency estimates about the best-fitted line is very small. The coefficient of correlation R was 0.999, and P < 0.0001. Vertical bars indicate the limits within which 50% of the energy of the wavelet is contained for wavelets at each node.
The spectra were estimated in several steps. A discrete WP transform of the signal was first calculated. The WP coefficients were then taken to approximate the coefficients of the Fourier waves, whose frequencies are the same as the "pseudofrequencies" of the wavelets used. All possible combinations of low- and high-pass filters, at the last decomposition level, were used. However, the pseudofrequency is only an approximate measure of the spectral content of the wavelet. To assess how choosing such an approximation affects the accuracy of the spectra, we compared two spectral estimates of the same signal. First, and as described above, it was assumed that a wavelet is well represented by a single frequency—pseudofrequency. Second, all spectral contributions of individual wavelets were taken into account. At each node the individual coefficients were reconstructed, their individual spectra were calculated and added up. The spectra are essentially identical (Fig. 2 D), demonstrating that the pseudofrequency approximates well the spectral content of the wavelet.
Comparison of Fourier and wavelet analysis
The advantages of the wavelet transform over Fourier transform can be best understood by presenting an example of signal decomposition using both transforms and by comparing their results. In Fig. 3, a 64-sample signal is shown (top trace). The Fourier decomposition series of eight of the most prominent spectral coefficients are shown below (A, middle traces). With each trace, the corresponding amplitude and phase of Fourier coefficients are shown on the left. The approximation of the original signal provided by the sum of the decomposition series is seen (A, bottom trace). We decomposed the same signal using the wavelet transform with two short wavelets—Haar (B) or Daubechies (db3, C). Haar wavelet can be expressed as a function 
(t) = 1 for 0 
t < 1/2; -1 for 1/2 t < 1, and 0 otherwise. The corresponding scaling function is 
= 1 for 0 t 1 and 0 otherwise. However, db3 wavelet has no explicit expression. Eight shifted and scaled decomposing wavelets with the most prominent real value coefficients (middle traces) and their sum (bottom trace) are also shown. The prominent coefficients were obtained by comparing the real value of all coefficients of the decomposition pyramid (i.e., all coefficients of all details and all coefficients of the last approximation) and taking the eight that were the largest. For each trace, we give the amplitude of its coefficient (on the left of the decomposition trace). On the right of the trace, we labeled whether it came from the approximation (indicated as small A) or detail (indicated as small D), and gave the decomposition level as well as its index at that level. For example D5,2 indicates that the coefficient came from the detail, the decomposition level was 5, and its index at that level was 2. Finally, further to the right, we give all the coefficients from that decomposition level. Note that: a), a filled symbol indicates the coefficient used for a reconstruction of the decomposition term shown on the left, and b), the number of coefficients decreases due to down-sampling as the decomposition level rises.
The sums of the decompositions are generally good for both Fourier and wavelet analysis. This is not surprising because both Fourier and wavelet transforms map a time-domain function into the frequency domain. However, there are important differences. The decomposition terms of wavelet transform contain information localized in time. In contrast, in the Fourier series terms the information about the time localization is spread over the whole duration of the signal, and although not lost, it is not explicitly available and is hidden in the phase(s) of the decomposition series. Therefore, Fourier analysis isolates well in frequency, but not in time. This is owing to the nature of support of the basis functions. The support is infinite for Fourier transforms because sines and cosines are of infinite duration. The support is finite (i.e., compact) for wavelet analysis because wavelets are localized waves, although some wavelets (Morlet wavelet, Meyer wavelet, Gauss wavelets, Mexican hat wavelet) do not have compact support. In this case, the performance of Haar decomposition is exact due to a great similarity between the analyzed signal and the wavelet itself. This points to the importance of choice of wavelet for the decomposition of the signal. Generally, choosing longer wavelets leads to a better frequency resolution but a worse time resolution.
Wavelet analysis provides information that is localized in frequency and in time, which makes it highly suitable for analysis of nonstationary signals. Fourier analysis can also be used if the signal is segmented, and if it is assumed stationary per segment. The limitation of this technique ("Short-Time Fourier Transform" (Gabor, 1946)) is that once a choice of the size of the window is made, it is the same for all frequencies. Using a shorter window results in better time, but lower frequency resolution. In contrast, wavelet transform provides a multiresolution analysis of nonstationary signals with frequency-dependent windowing (long windows for low frequencies and short windows for high frequencies).
Monte Carlo simulations and graphical presentation
Synaptic transmission in the excitatory synapse in the hippocampus was simulated using a Monte Carlo method (Bartol et al., 1991; Wahl et al., 1996), because the method enables one to study the fluctuations due to the stochastic nature of the diffusion of the individual transmitter molecules, their interactions with the receptor molecules, and the stochastic nature of the gating of the associated ion channels. At each discrete time step: a), every transmitter molecule has a position (x,y,z), and is flagged either as free or bound—free molecules move randomly in all three dimensions or interact with the receptors; and b), every receptor has a fixed position and is considered to be in one of the states of a given kinetic scheme—every receptor has a finite probability of changing to another state according to the same kinetic scheme. The changes of the receptor states were assumed to be Markovian (i.e., they were assumed to depend only on their present position and not on their previous history). Since the rise time of mEPSCs or patch currents is of the order of 
100 µs, we chose the time step in our simulation to be 0.1 µs (see below (Wahl et al., 1996)).
The distance traveled by a transmitter molecule (in each of three dimensions) was chosen randomly from a Gaussian distribution with mean of 0 and a standard deviation 
given by
 | (5) |
(Hille, 1992; Wahl et al., 1996; Glavinovic and Rabie, 1998), where the 
t is the length of the time step and D is the diffusion coefficient of the transmitter molecule. The random numbers were obtained from a random number generator. The diffusion in a restricted space was simulated by assuming that the transmitter molecules collide elastically with the "walls" of the space (presynaptic or postsynaptic membranes (Wahl et al., 1996)). Finally, all glutamate molecules that reach the edge of the model synapse and diffuse away into the "infinite" space (or return from it into the synaptic cleft) were also followed throughout simulation.
Synaptic transmission was simulated with glutamate released from an instantaneous point source facing the synaptic cleft, in the center of the synapse, on the membrane opposite to where the postsynaptic AMPA receptors lie (Fig. 4 A). The synaptic dimensions were 200 x 200 x 15 nm; 14 x 14 (or 196) AMPA receptors were distributed equidistantly on the postsynaptic surface (Wahl et al., 1996; Glavinovic and Rabie, 1998). The kinetic scheme used was as given by Jonas et al. (1993) (Fig. 4 B). In all simulations, the temperature was assumed to be 37°C. All rates were adjusted by assuming Q10 = 3.0. The diffusion coefficient of glutamine (D = 760 µm2/s measured in water at 25°C (Longsworth, 1953)), was taken to be the diffusion coefficient of glutamate. The adjustment for the temperature difference was done by assuming Q10 = 1.3 (Hille, 1992). The single channel current amplitude was 1 pA. In the simulations of voltage-clamp experiments on an outside-out membrane patch, a 30 nm box was constructed above the 200 x 200 nm postsynaptic membrane containing also 14 x 14 or 196 AMPA receptors. We fixed the concentration of glutamate at a given level while keeping the size of the box constant. To maintain the concentration of glutamate, any binding or unbinding of glutamate molecules was not associated with removal or addition of glutamate molecules from or to the "box." The glutamate molecules were distributed randomly within the box at the beginning of the simulation. To simulate abrupt changes in the glutamate concentration (pulses), we distributed additional molecules at random positions in the box or removed a certain number of free molecules.
View larger version (27K):
[in this window]
[in a new window]
|
FIGURE 4 Models of glutamate diffusion in the synaptic cleft and the kinetics of gating of AMPA receptors. (A) Diagram (approximately to scale) depicts how glutamate molecules, after release from an instantaneous point source, diffuse in synaptic cleft. (B) Kinetic scheme of gating of AMPA receptor-channels by glutamate, used for the Monte Carlo simulations. U, SB, DB, and O indicate the unbound, singly bound, doubly bound, and open state, respectively. D1, D2, and D3 are three desensitized states. The rate constants, taken from Jonas et al. (1993), were adjusted for the temperature of simulations (37°C), assuming a Q10 of 3.0. They are: K+1 = 2.07 x 107 M-1s-1, k-1 = 1.92 x 104 s-1, K+2 = 12.9 x 107 M-1s-1, k-2 = 1.47 x 104 s-1, K+3 = 1.14 x 107 M-1s-1, and k-3 = 207 s-1 for glutamate binding;  = 1.91 x 104 s-1 and ß = 4050 s-1 for channel opening; and 1 = 1.30 x 104 s-1, ß1 = 176 s-1, 2 = 774 s-1, ß2 = 3.27 s-1, 3 = 79.7 s-1, ß3 = 32 s-1, 4 = 75.6 s-1, and ß4 = 855 s-1 for the desensitization pathway.
|
|
According to the kinetic scheme of channel gating (Fig. 4 B), a receptor can be unbound (U), in a single or in a double bound state (SB and DB, respectively; these are also called activatable states), open (O), or in one of the three desensitized states (D1, D2, and D3). Each state of the receptor is associated with a surface area and a probability of binding, given that a transmitter molecule "hits" this receptor surface (Bartol et al., 1991; Wahl et al., 1996). The inverse of the receptor surface area is taken to be the density of the receptor molecules (r) at the postsynaptic membrane. The probability (Pb) that a transmitter molecule, after hitting the receptor surface, will bind in a given time step (t) is related to the macroscopic rate constant by
 | (6) |
where Na is Avogadro number and 
is the appropriate rate constant of binding (in M-1 s-1). For steps in the kinetic scheme that did not involve binding of the transmitter, the probability (p) that a receptor will move to a new state in a given time step (t) is related to the macroscopic rate constant by
 | (7) |
where k is the appropriate rate constant (in s-1). The time step (chosen to be 0.1 µs in all our calculations) was such that the probability of a receptor changing state twice within one time step was <1 in 100. It has been shown previously that the results do not become more accurate by reducing the time step further (Wahl et al., 1996). Equation 7 also applies to the situation when a transmitter molecule unbinds from a receptor. We have chosen to move a molecule the mean length of the random jump (0.67 ) perpendicularly away from the receptor surface at unbinding: 1), to ensure that the probability of a given receptor making a transition to a bound state does not depend on the receptor's previous history; and 2), because such physical separation between the neurotransmitter and the receptor accurately reproduces the macroscopic unbinding rate constant (Wahl et al., 1996).
All graphics were done using Origin (Microcal Software, Northampton, MA). Skew was calculated as the ratio of the third sample moment about the mean (m3) and the SD3 (SD, standard deviation of the sample; m3 = (
(xi - xmean))3/n; xi, individual values of the spectral densities; n, number of spectral density points. All simulations were done using MATLAB (MathWorks, Natick MA) on a Pentium IBM computer. The programs for wavelet analysis were locally made using the operators available in MATLAB Wavelet Toolbox.
|
RESULTS
|
|---|
Wavelet analysis of stationary fluctuations
Fig. 5 A shows glutamate-activated AMPA patch currents simulated using the Monte Carlo method. Glutamate pulses were 10 ms long and were spatially uniform. The glutamate concentration was constant but differed from run to run, ranging from 63 to 500 µM. After an initial rise, each patch current reaches a steady state. On visual inspection, the frequency of the current fluctuations changes with glutamate concentration and current level. To evaluate these changes, we determined the power spectral density of the current fluctuations using Fourier and WP analysis from the zero-mean steady-state segments of current fluctuations—residuals (8 ms long, starting 2.0 ms from the beginning of the simulation and lasting until the end; Fig. 5 B). Note that WP spectral densities are shown versus frequency, calculated by creating a link between the scale description of wavelet (and WP) analysis and the frequency description of Fourier analysis at each decomposition level (see Methods). The analysis performed throughout this study requires very good time resolution, whereas the frequency resolution is less critical, because the spectral densities change rapidly with time but do not have sharp peaks. We thus used two short wavelets, Daubechies-2 (db2, Fig. 5 B, filled circles) and Haar (open upward triangles), as the bases for decomposition. The number of decomposition levels influences also the time and frequency resolutions. We chose nine levels of decomposition that yield 29 or 512 spectral points. Raising the number of decomposition levels provides a greater number of spectral points and expands the examined spectral range at low frequencies, but it also lowers the time resolution. Finally note that if the length of the signal is not an exact multiple of a, the signal is extended on both ends by adding zeros (zero padding), which may lead to edge artifacts (see below). Both WP spectral densities are similar to each other and to the Fourier spectral density, but the WP spectral densities are less scattered. Finally all spectral densities agree well with the theoretical steady-state spectral density (Fig. 5 B, line). The theoretical spectral density was calculated using the eigenvalues and eigenvectors that can be obtained from the Q-matrix of the kinetic scheme of gating of AMPA channels given in Fig. 4 B (Colquhoun and Hawkes, 1977). A good agreement between the theoretical and WP (but also Fourier) spectral estimates raises our confidence in the validity of the present analysis.
View larger version (33K):
[in this window]
[in a new window]
|
FIGURE 5 (A) Time course of Monte Carlo simulated patch currents in the absence of desensitization (all rates leading to and from the desensitized states were set to zero). The glutamate pulses (concentration as indicated) started at zero time, were 10 ms long, constant, and spatially uniform. All currents reached a steady state in <2 ms. (B) The wavelet packet-based spectral density distributions of the stationary current fluctuations (db2, filled circles; Haar, open upward triangles) resemble the Fourier-based distribution (crosses) but are less scattered. Furthermore, all three spectra are close to the theoretical spectra calculated from rates of the kinetic scheme (see text). Glutamate concentration was 250 µM. (C) The mean-variance plot of the current fluctuations calculated using the wavelet packet analysis (WP-db2-l9, filled circles) is well fitted by the parabola and is close to the true parabola (see text). Fourier mean-variance plot of the nonsegmented signal (Fou, open circles) is also close to the true parabola, but not if the signal is segmented into 18 segments (Fou-Seg, filled upper triangles). If the individual segments are further made zero-mean (Fou-Seg-Zero Mean, filled downward triangles), the variance is also reduced, making the mean-variance plot even more inaccurate. (D) Coefficient of variation and skew of the spectral densities of the current fluctuations are not constant, which indicates a changing shape and lower dispersion of spectra. (E) Median frequency (f50, the frequency at which half of the energy is below and half above) versus mean current (and AMPA receptor saturation) relationship is nonlinear, and similar irrespective of whether Fourier (nonsegmented signal) or WP analysis is used. If the signal is divided into 18 segments, as for WP analysis, f50 (Fourier) is similar to f90 (WP). Mean current of 196 pA indicates 100% saturation (all channels open). The relationship between 1/f50 and the mean current, however, appears to be linear irrespective of whether Fourier (nonsegmented signal) or WP analysis is used but the fits are different (see text).
|
|
Because it is of uppermost concern, the accuracy of Fourier and WP estimates of the spectral densities was also evaluated by comparing the mean-variance plots calculated from such spectral densities with true mean-variance plots. Fig. 5 C shows the true mean-variance plot (thick line) calculated using the equation 2 = iI - I2/N (Robinson et al., 1991) together with the Fourier and WP estimates (as indicated). 2 and I are the variance and the mean of the current, i is the single channel current, and N is the number of channels. The variance (or total energy of the signal) of the current fluctuations was calculated as an integral of the power spectral density. The mean-variance plot from WP spectra approximates well the true plot. Those estimated by Fourier analysis, but using a nonsegmented signal, are similarly close to the true plot. The estimates of the single channel current were 0.92 pA and 0.79 pA, whereas those of N are 225 and 263 (WP and Fourier, respectively). However, Fourier estimates provide only one value of variance for the whole 8 ms signal, whereas WP analysis provides several (in this case 18) values (only their mean value is shown).
Having only a single value is a major drawback for the analysis of nonstationary signals such as deactivation or desensitization patch currents and synaptic currents (see below). The Fourier estimates of the variance, with the same degree of segmentation as that of WP analysis, were poor and mean-variance plots inadequate. This is not surprising because the individual segments were not zero-mean, and even small errors of the mean can lead to large errors of the estimated variance (Neher and Stevens, 1977; Conti et al., 1980; Silberberg and Magleby, 1993). When the individual segments were made zero-mean, the variance was reduced to values below true values, and the plot was even less accurate. The estimates of the single channel current were 1.50 pA and 0.27 pA, whereas those of N were 123 and 1625 (Fourier segmented nonzero mean and zero-mean data, respectively). Nevertheless, and in agreement with earlier studies (Sigworth, 1980; Traynelis and Jaramillo, 1998; Wierenga and Wadman, 1999), the mean-variance plot could in all cases be well fitted by a parabolic function (Fig. 5 C).
To evaluate how the frequency of the current fluctuations changes with glutamate concentration, we calculated the median frequency f50 (the frequency at which 50% of energy of current fluctuations is below and 50% above it) or f90 (the frequency at which 90% of energy of current fluctuations is below and 10% above it). f50 was used instead of mean frequency, as it describes more reliably the location of the main frequency components. f90 was also used as it provided less scattered estimates than f50, especially when a very high fraction of the total variance was attributable to a few spectral values (or even a single one) at the low frequency end. The conventional way of characterizing the frequency of the current fluctuations is to calculate the critical frequency fc (the frequency at which the power declines to a half of its asymptotic zero-frequency value). However, we decided not to use it (except for a comparison) because such a characterization is model bound. In the present studies, the kinetic scheme of AMPA channel gating was in most cases complex; moreover, given that three kinetic rates were concentration dependent, it changes with time. The glutamate concentration was usually time-dependent and spatially nonuniform (see below).
f50 rises as the mean current (and saturation level) increases irrespective of what method is used to calculate the spectral densities (Fig. 5 E). WP estimates (which are segmented) are similar to those of Fourier analysis (nonsegmented signal), but Fourier f50 estimates of the segmented signal are much higher and resemble f90 WP estimates. The relationship between the inverse of the median frequency (1/f50) and the mean current for the WP and Fourier nonsegmented estimates could be fitted by a straight line (Fig. 5 F). The best-fitted lines are: a), 1/f50 = 2.200 - 0.012 x I (Fourier), and b), 1/f50 = 1.542 - 0.007 x I (WP), where I is the mean current. This agrees well with previous studies which demonstrated that: a), fc of the current fluctuations rises nonlinearly as the probability of opening of ion channels and saturation levels rise, and b), the relationship between the observed mean burst length (
d; calculated from fc as d = 1/(2
fc) and the saturation level is negative but linear (Silberberg and Magleby, 1993). We also calculated the spectral coefficient of variation (the standard deviation/mean of the spectral densities) and spectral skew (see Methods). Both diminish with greater saturation, which precludes a straightforward comparison between fc and f50 and thus a direct calculation of d from f50 (Fig. 5 D). Coefficient of variation describes the dispersion, and the skew the shape of the spectral densities.
Spectral density and mean-variance plots estimated using a modified Fourier analysis
The inaccuracy of the spectral density estimates of short signals can be partially remedied by a modification of traditional Fourier analysis, which better estimates a contribution of the baseline to the power spectra. If two signals are identical but their baselines are different, their spectral densities will be identical (except at the lowest frequency end), because all spectral values result from zero-mean operations—band-pass and high-pass filtering. Despite the extremely localized effect of the baseline on the power spectra, the estimates of the parameters characterizing the power spectra, such as the variance, f50, f90, and fc, will nevertheless differ, often greatly, owing to the fact that the spectral values at the low frequency end can be comparable or greater than all other points combined. Although the true value of the spectral density at the low frequency end is not known especially for short signals, it can be inferred from the general trend of the spectral estimates by fitting the Lorentzians that exclude the spectral value at the lowest frequency end.
Fig. 6 A compares the true mean-variance plots (thick line) with the plots obtained from hybrid spectra calculated using the modified Fourier analysis at different levels of segmentation (as indicated) and with various types of fitting (filled symbols, on linear scales; open symbols, on logarithmic scales). The mean-variance plot from the nonsegmented data approximated well to the true plot when the data were fitted on linear scales, but less well when fitted on the logarithmic scales. The variances calculated from the spectral densities (and fits) on the logarithmic scale were below the true values and below those estimated on the linear scales by up to 20%, and the errors were greater for currents induced by glutamate pulses of low concentration, presumably because of the greater low-frequency spectral content. This is as expected. Fitting the Lorentzians to data on the logarithmic scale attributes a disproportionally high importance to the small values of the spectral densities. Increasing the segmentation of data leads to progressively greater discrepancies between the true and observed mean-variance plots, but the plots remained good with three segments used (i.e., for currents 2.67 ms long). These discrepancies are reflected in the estimates of the single channel currents and the number of AMPA receptors postsynaptically. The single channel and N estimates were: a), 1.11 pA and 172 (nonsegmented), b), 1.10 pA and 179 (three segments), c), 0.65 pA and 370 (nine segments), d), 0.39 pA and 1086 (18 segments), and e), 0.90 pA and 227 (nonsegmented). All variance estimates were calculated from the hybrid power spectra and used the Lorentzian fits, which were done on the linear (a–d) or logarithmic scales (e).
View larger version (25K):
[in this window]
[in a new window]
|
FIGURE 6 The length of current segments affects the spectral density and variance estimates. (A) The mean-variance plots calculated using the improved Fourier analysis (see text) and Lorentzian fitting on linear (filled symbols) or logarithmic scales (empty symbols). In both cases, the parabolic fits are close to the true parabola for nonsegmented data (labeled as 1), for the segmented data with a few (3) but long segments, but not for segmented data with many (9 or 18) but short segments. (B) Mean-f50 plots. (C) Mean-fc plots. Vertical bars are standard errors.
|
|
Both f50 and fc rose as mean current increased (Fig. 6, B and C). Similarity of the mean-f50 plots, and to a lesser extent of the mean-fc plots, for signals >2.67 ms (duration of the signal segmented into three sections) suggests that they represent the true plots. If so, the spectral densities of the signals whose fc is >0.5 kHz can be determined accurately. The length of the signal that is needed for an accurate estimation of the variance, but also f50, f90, and fc, is thus <8 x d (d – mean burst length). Finally note that, as observed for mean-variance plots, mean-f50 and especially mean-fc plots were also less accurate when obtained from spectra fitted on logarithmic scale.
Wavelet analysis of patch currents produced by long and short glutamate pulses
A decrease in responsiveness (desensitization) of the evoked currents can be observed when prolonged, though constant and spatially uniform transmitter pulses are applied to the transmitter-activated channels with a complex gating mechanism, such as ACh or AMPA channels (Katz and Thesleff, 1957; Trussell and Fischbach, 1989; Raman and Trussell, 1995). As single-channel studies of the acetylcholine receptor have shown, in such cases the stochastic nature of gating changes and displays a progressively changing pattern of bursts and clusters (Sakmann et al., 1980). Because the opening and closing of ion channels is nonstationary, the spectra of the current fluctuations should also be nonstationary. Fig. 7. A shows a simulated patch current resulting from a spatially uniform 1 mM glutamate pulse, 10 ms long, starting at zero time. The trace was fitted by the product of two exponentials: I = 11.3 + [126.1 exp(-t/3.04)] [1 - exp(-t/0.138)], where I is the current in pA and t time in ms. Fig. 7 B shows the residuals. Both the frequency of current fluctuations and their variability appear to change in a time-dependent manner. Indeed, the variance (Fig. 7, C and D), f50 (Fig. 7, E and F), and f90 (Fig. 7 G) were time-dependent and their time courses similar to that of desensitization current. As expected, f90 values were higher and less variable than those of f50. The corresponding db2 and Haar estimates (i.e., variance, f50, and f90) were similar except for the discrepancies at the start of the signal. The edge effects degraded the db2 estimates of both the variance and f90 at the start due to zero-padding. The Haar wavelet is shorter and less prone to such problems. The estimates of variance with seven decomposition levels were similar to those with nine levels but were more variable. They are based on the spectra with fewer values. In contrast, f50 estimates with seven decomposition levels were not only more variable but were also generally somewhat higher. It is difficult to evaluate f50 accurately with few spectral values and with f50 near the low frequency end.
View larger version (35K):
[in this window]
[in a new window]
|
FIGURE 7 (A) Simulated patch current during desensitization. Constant and spatially uniform pulse of 1 mM glutamate was applied throughout. After a rapid onset, the current decreased, but more slowly, to a low maintained level. (B) Residuals. (C and D) Variance of the current fluctuations changes with time and the estimates using Haar and db2 wavelets are very similar except at the beginning. Throughout the time course, the values of variance were very similar irrespective of whether nine-level and seven-level decompositions were used, but nine-level decompositions provided less variable values but with poorer time resolution. Note the break on the ordinate axis. (E–G) f50 and f90 are time dependent and similar for Haar and db2 wavelets except at the beginning, where db2 inaccurately yielded very low values. Variance estimates were more variable and consistently, though moderately, higher for seven-level decompositions. Open and filled symbols, seven- and nine-level decompositions, respectively. Averages of 10 runs. Vertical bars are standard errors.
|
|
Fig. 8, A and B, depicts the mean-variance plots with the corresponding best parabolic fits for the decay phase of desensitization currents (A, nine decomposition levels; B, seven decomposition levels). The fits are similar irrespective of the decomposition level or whether db2 or Haar wavelet was used. The estimates of the single-channel current were 0.45 pA (individual run; db2 wavelet), 0.43 pA (average of 10 runs; db2 wavelet), and 0.41 pA (average of 10 runs; Haar wavelet) for nine-level decomposition, whereas for the seven-level decomposition they were 0.44 pA (average of 10 runs; db2 wavelet) and 0.41 pA (average of 10 runs; Haar wavelet). The number of AMPA receptors with nine-level decomposition was estimated to be 909 (individual run; db2 wavelet), 455 (average of 10 runs; db2 wavelet), and 526 (average of 10 runs; Haar wavelet), whereas for seven-level decomposition they were 435 (average of 10 runs; db2 wavelet) and 526 (average of 10 runs; Haar wavelet). The mean-f50 relationship is positive irrespective of whether db2 or Haar wavelet or whether seven- or nine-level decompositions were used, but the hysteresis could be assessed accurately only when using Haar wavelet (see arrows; Fig. 8, C–E).
View larger version (35K):
[in this window]
[in a new window]
|
FIGURE 8 (A) Mean-variance plots of the desensitizing current fluctuations during the decay phase (same data as in Fig. 7) with parabolic fits. Using Haar and db2 wavelets gives similar plots and fits; nine levels of decomposition. (B) Same mean-variance plots and parabolic fits but with seven levels of decomposition. Note greater number (and better time resolution) but more variable mean-variance data pairs, and parabolic fits that are similar to those shown in A. (C) Mean-f50 relationship is positive, but the hysteresis is revealed with Haar but not with db2 wavelet (see text). (D and E) Similarly positive relationship between current and f50 is observed with seven-level decomposition irrespective of whether Haar or db2 wavelets are used, but the hysteresis is clearly visible only with Haar and not with db2 wavelet. Open circles, individual estimates; filled symbols, averages of 10 runs; upward triangles, Haar; downward triangles, db2. Vertical bars are standard errors.
|
|
Short glutamate pulses are often applied to excised macropatches serving as surrogate synapses (Colquhoun et al., 1992; Hestrin, 1992; Jonas et al., 1994; Raman and Trussell, 1995). Simulated deactivation current induced by a spatially uniform 2 mM, 0.5 ms glutamate pulse starting at zero time is shown together with the best fit (Fig. 9 A). At the peak, a saturation level of 70% was reached (135 AMPA receptors were open out of 196). The deactivation current was fitted by the following function: I = 0.041 + [226.9 exp(-t/0.671)] [1 - exp(-t/0.134)], where I is the current in pA and t is the time in ms. Fig. 9 B depicts the residuals. Note that though the kinetic scheme of gating of AMPA channels is the same as that used for the simulations of desensitization currents (Fig. 7), the deactivation current clearly decays faster than the desensitization current, as observed experimentally (Colquhoun et al., 1992; Jonas and Spruston, 1994; Edmonds et al., 1995; Raman and Trussell, 1995; Hausser and Roth, 1997). Faster time course of the current, but especially the abrupt change of three concentration-dependent kinetics rates of gating, are expected to lead to rapid spectral changes of the current fluctuations. Indeed, f50 and f90 both change markedly with time, reaching levels clearly above those observed for the desensitization currents, although the saturation level was approximately the same (Fig. 9, C–E). The variance, however, changed less. Note that the db2 and Haar wavelet estimates of variance, f50, and f90 were similar. The parabolic fit to the mean-variance data pairs (points are averages of 10 runs) yields the estimate of the single channel current of 0.41 pA and 0.29 pA (db2 and Haar wavelets, respectively) and of the number of AMPA receptors postsynaptically of 500 and 1136 (db2 and Haar wavelets, respectively; Fig. 9 F). Irrespective of whether db2 or Haar wavelet were used, the mean-f50 and mean-f90 relationships are both positive, and both plots reveal clear hysteresis with f50 and f90 values that are higher during the rise phase (see arrows, Fig. 9, G and H).
View larger version (32K):
[in this window]
[in a new window]
|
FIGURE 9 (A) Simulated deactivation current. Constant and spatially uniform pulse of 2 mM glutamate was applied for 0.5 ms. A rapid rise of current is followed by a rapid decay. (B) Residuals. (C–E) The variance, f50, and f90 of the current fluctuations. (F) Mean-variance plots with the parabolic fits. (G and H) Mean-f50 and mean-f90 plots both reveal a positive relationship and presence of hysteresis. All values of variance, f50, and f90 are averages of 10 runs. Vertical bars are standard errors. Using db2 wavelet (filled circles) and Haar wavelet (crosses) provided very similar results.
|
|
Wavelet analysis of mEPSCs
As shown above, the estimates of the parameters characterizing the fluctuations of deactivation and desensitization currents such as variance, f50, and f90, though good, are not excellent presumably because: a), a high fraction of AMPA receptors enters into desensitized states from which they do not return or only very infrequently, and b), the glutamate concentration changes rapidly. Estimating the same parameters from the synaptic currents is more difficult because the glutamate pulse in the synaptic cleft is spatially nonuniform and changes rapidly with time. Fig. 10 A shows a simulated mEPSC after a release of 1500 glutamate molecules. The best fit to the fast-mEPSC calculated using the least-squares method (smooth curve) was: I = 0.046 + [33.3 x exp(-t/0.518)] x [1 - exp(-t/0.046)], where I is the current in pA, and t is the time in ms. Overall saturation of AMPA receptors at the peak of the mEPSC was 15%. The residuals suggest that both the amplitude and the frequency of current fluctuations are nonstationary (Fig. 10 B). Indeed, the variance was time dependent and increased more when more glutamate molecules were released (Fig. 10 C). However, the variance changed less whereas f90 changed more for the same current than observed for patch currents (stationary or nonstationary). Finally, f90 changed largely independently of the number of glutamate molecules released (Fig. 10 D). All variance, f50, and f90 estimates were very similar irrespective of whether db2 or Haar wavelet was used.
View larger version (33K):
[in this window]
[in a new window]
|
FIGURE 10 (A) Simulated miniature fast excitatory postsynaptic current (mEPSC) with the best fit by the sum of two exponentials (see text). (B) Residuals. (C and D) Time course of variance and f90 of the current fluctuations. (E) Mean-variance plots with the best parabolic fits. (F) Similar mean-f90 plots are obtained with db2 wavelet (filled symbols) or Haar wavelet (open symbols). Upward and downward triangles, 1500 and 6000 glutamate molecules released, respectively. Note the hysteresis of mean-f90 plots and the similarity of f90 values at the peak. All mean-variance and mean-f90 data pairs are averages of 10 runs. Vertical bars are standard errors.
|
|
Fig. 10 E shows the mean-variance data pairs with parabolic fits. The single-channel currents estimated from the fits were 0.20 pA (N – 1500; db2), 0.27 pA (N – 1500; Haar), 0.33 pA (N – 6000; db2), and 0.36 pA (N – 1500; Haar), whereas those of the number of AMPA receptors postsynaptically were 667 (N – 1500; db2), 131 (N – 1500; Haar), 435 (N – 6000; db2), and 357 (N – 6000; Haar; N – number of glutamate molecules released). Note that the estimates of the single-channel currents and of N are usually better when more glutamate molecules are released. The saturation levels reached are higher when more molecules are released. The mean-variance data pairs thus have values over a greater part of parabola and yield more accurate estimates of the single-channel current, and of the number of AMPA receptors postsynaptically. Furthermore, using Haar wavelet provided overall better estimates. Nevertheless the estimates of the single-channel current are clearly below the true values, whereas those of the number of AMPA receptors are typically above. Fig. 10 F shows mean-f90 plots. Because f90 at the current peak is essentially independent of the number of glutamate molecules released, the slope of the mean-f90 plot is steeper when fewer molecules are released. Hysteresis of the mean-f90 plots was, however, clear irrespective of the number of glutamate molecules released.
The current fluctuations of experimentally recorded mEPSCs are not solely due to the opening of AMPA channels by the glutamate. A variety of electronic or biological sources also contribute as is evidenced by the baseline that is never noiseless. We now evaluate how reliable wavelet analysis is under such clearly more realistic conditions. mEPSCs were made "noisy" by adding a white Gaussian random noise with a zero mean and whose standard deviation was 0.5 pA. The baseline fluctuations are often similar in amplitude to the amplitude of single-channel currents. Fig. 11 A shows a "noisy" mEPSC, with the best fit. The amplitude histogram of the added noise is given in Fig. 11 B with the corresponding Gaussian probability density function (line). On visual inspection, the residuals (Fig. 11 C) reveal a change of the variance but not of the frequency of current fluctuations, and the unitary events are not clearly discernible at the tail of the mEPSC. The variance changes of the noisy mEPSCs (filled triangles; Fig. 11 D) generally resemble those of the noiseless mEPSCs shown in Fig. 10 C, and are not very much affected by the contribution of the baseline noise. As can be seen very clearly from Fig. 11 E, the relationship between the variance of the noiseless mEPSCs (abscissa) and the variance difference between noisy mEPSCs and the baseline (ordinate) is very close to the line with the unity slope. The variance contribution of the baseline noise can be successfully accounted for.
View larger version (32K):
[in this window]
[in a new window]
|
FIGURE 11 (A) Simulated noisy mEPSC with the best fit by the sum of two exponentials (see text). Gaussian random noise with a zero mean and a standard deviation of 0.5 pA was added to the simulated mEPSC. (B) Amplitude histogram of the added noise with the line showing the corresponding Gaussian probability density function. (C) Residuals. (D) The variance of the current fluctuations when 1500 and 6000 glutamate molecules were released (downward and upward triangles, respectively). Filled symbols, noisy mEPSCs; open symbols, variance difference between noisy mEPSCs and baseline; crosses, added Gaussian noise (baseline). (E) Variance of the noiseless mEPSCs (abscissa) versus variance of noisy mEPSCs but with the contribution of the baseline current fluctuations accounted for. Filled and open symbols, 1500 and 6000 glutamate molecules released, respectively. The line has unity slope. (F) f90 of current fluctuations of noise-free (open symbols) and noisy mEPSCs (filled symbols) are different, especially when 1500 glutamate molecules (downward triangles) were released. f90 of baseline fluctuations (crosses) simulated as white Gaussian noise is very high. f90 of noisy mEPSCs decreases initially before increasing as the baseline contribution becomes more important. (G) f90 of the noise-free mEPSCs (open symbols) and noisy mEPSCs (filled symbols) but corrected for the baseline contribution agree well, except when mEPSC fluctuations become comparable or smaller than baseline fluctuations. The agreement was further improved when f90 was calculated from spectral averages (crosses; see text). All variance and f90 values are averages of 10 runs. Vertical bars are standard errors. Only Haar wavelet analysis results are shown to avoid crowding.
|
|
However, f90 of noisy and noiseless mEPSCs were in this case very different, especially with 1500 glutamate molecules released and when mEPSCs decayed to low levels (Fig. 11 F). This is not surprising. The baseline was simulated as a white Gaussian noise and its f90 is much higher than that of the noiseless mEPSCs. When the contribution of the baseline current fluctuations to noisy mEPSCs is taken into account, the agreement between f90 of the noiseless mEPSCs and f90 of the corrected noisy mEPSCs improved. The agreement was further improved when f90 was calculated from the averaged spectra rather than as an average of f90s, each evaluated from individual mEPSCs (Fig. 11 G).
Origin of the hysteresis of current-median frequency relationship
The existence of the hysteresis in the mean-f50 and mean-f90 plots demonstrates that f50 and f90 do not simply reflect the level of saturation of AMPA receptors. However, the origin of the hysteresis remains unclear. Two possible sources of the hysteresis are considered: a), the time course of the occupancy of the doubly bound state, and b), the temporal changes of the spatial distribution of AMPA receptor activation during the time course of an mEPSC.
Currents fluctuate because the number of channels open at any time varies stochastically. Given a very simple kinetic scheme (one unbound and one open state) and stationary conditions, the observed time constant of the return of the current to the mean current level is given by a simple relation (Hill, 1909; Hodgkin and Huxley, 1952):
Under such restricted conditions, the time constant (and the frequency) of current fluctuations depends only on two kinetic rates, and ß. This is owing to the fact that the occupancy of two states is determined exclusively by the two rate constants. The relationship between saturation and observed is negative and linear (Silberberg and Magleby, 1993). If the probability of opening (and thus the level of saturation) is low, the observed is determined by the closing rate (ß) only.
When the kinetic scheme of gating is complex and the conditions are not stationary, the frequency of the current fluctuations depends not only on the kinetic rates but also on the occupancy of various states. The occupancy of the doubly bound state preceding the open state is clearly the most important, whereas the occupancy of the desensitized state directly linked to the open state is as a rule very low and can be disregarded (Wahl et al., 1996; Glavinovic and Rabie, 1998). However, during deactivation and synaptic currents, the occupancy of the doubly bound state rises very rapidly, and is clearly higher during the rise than during the decay time for the current level (Glavinovic and Rabie, 1998; see also Fig. 10 F). This leads to higher frequencies during the rise phase for the same current producing hysteresis of the mean-f50 plot.
We also explored whether the changes of the spatial distribution of open AMPA receptors could be a contributing factor. Fig. 12 A shows the spatial distribution of the occupancy of the open state at different times (as indicated; 1500 glutamate molecules were released) at the start of glutamate release. Corresponding spatial distributions of the double bound, unbound, and desensitized (combined) states are shown in Fig. 12, B–D, respectively. Fig. 12 E shows how the radii of the areas, where 30%, 50%, and 90% of activation occurs, change during the time course of an mEPSC. There is surprisingly little if any spatial spread of the activation disk. Furthermore, the radius where 30% of the activation occurs is similar throughout the time course of an mEPSC, irrespective of whether 1500, 3000, or 4500 glutamate molecules are released. Fig. 12 F shows the time course of the probability of opening of the disk where 30% and 50% of activation occurs together with the probability of opening for the whole synapse. As expected, the probability of opening is higher for the activation disk comprising 50% and especially 30% of all open channels. Finally note that though the occupancy of the doubly bound state is lower during the decay phase of an mEPSC, it is not constant but gradually decreases.
View larger version (39K):
[in this window]
[in a new window]
|
FIGURE 12 (A) Spatial profiles of the occupancy of the open state of AMPA channels estimated at different times (as indicated) after the instantaneous release of 1500 glutamate molecules (from a point source). (B–D) Spatial profiles of the occupancy of the doubly bound, unbound, and fully desensitized states, respectively. Note that the activation is initially highly localized but is associated with a high occupancy of doubly bound state and subsequently open state. In contrast, close to the edges of the synapse, binding of glutamate is initially associated with an entry into the desensitized states. Unless they are unbound, AMPA receptors tend to be in desensitized states, and remain so throughout the duration of an mEPSC. (E) Absence of spread of the radius of the activation disk containing 30%, 50%, and 90% of the open AMPA channels. The radii of the activation disks associated with a release of 3000 and 4500 glutamate molecules are similar to those resulting from a release of 1500 glutamate molecules. (F) Time course of the occupancy for the activation disk containing 30% and 50% of AMPA receptor openings, as well as for the whole synapse, and of the doubly bound state (DB). Note its rapid time course.
|
|
|
DISCUSSION
|
|---|
Evaluation of stationary current fluctuations by wavelet packet and by modified Fourier analyses
In this study we explore, on Monte Carlo-simulated excitatory patch and synaptic currents (Wahl et al., 1996; Glavinovic and Rabie, 1998), the use of wavelet analysis for the evaluation of the spectral changes of current fluctuations. We assessed first the spectral properties of excitatory patch currents under steady-state conditions. In this case the theoretical spectra can be calculated and the binomial theory predicts a parabolic mean current-variance relationship. A comparison between the theoretical and estimated spectra and mean-variance plots provides a test of the accuracy of WP and Fourier methods.
Glutamate pulses were long, constant, and spatially uniform; the kinetics of gating of AMPA channels were simple and had only an open state, an unbound state, and a bound but a closed state. There were no desensitized states. The probability of opening of AMPA channels in different runs varied over a wide range. Both WP and Fourier estimates of the spectral densities were similar and close to the theoretical predictions, but WP estimates were less scattered. Furthermore, WP analysis yields many sets of spectral densities, an important asset when determining the time course of nonstationary signals (see below).
Wavelet analysis of nonstationary fluctuations of AMPA patch currents
The real importance of the wavelet analysis lies in its ability to assess the fluctuations of the nonstationary signals (see Introduction and Methods). We tested its capabilities first on the desensitizing AMPA patch currents. The kinetic scheme of the gating kinetics of AMPA channels was complex and consisted of an open state, three closed states (an unbound, a single bound, and a double bound state), and three desensitized states as described previously (Wahl et al., 1996; Glavinovic and Rabie, 1998). Given a complex pattern of opening and closing of single channels during desensitization, characterized by bursts and clusters (Sakmann et al., 1980), the spectra of the current fluctuations ought to be time dependent, and indeed they are. Although f50 and f90 changed rapidly, their time course could be determined. Haar and db2 estimates of f50 and f90 were similar except at the beginning, where db2 estimates were less accurate owing to zero padding (db2 wavelet is longer than Haar wavelet). f50 and f90 estimates were more variable, and their values were modestly higher when the number of decomposition levels decreased to seven from nine but the time resolution improved. This is as expected because the values are based on fewer spectral points. Mean-f50 and mean-f90 relationships were positive except for the hysteresis shown by a higher f50 (or f90) values during the rise phase (for the same current). The hysteresis could be well tracked with Haar, but not db2 wavelets. The Haar but not db2 wavelet was similarly able to track variance changes at the beginning. Reducing the level of decomposition from nine to seven made variance estimates more variable, but the mean<
TOP
ABSTRACT
INTRODUCTION
