Fluctuation Analysis of Nonideal Shot Noise - BioMedSearch
Fluctuation Analysis of Nonideal Shot Noise Application to the Neuromuscular Junction R . FESCE, J . R . SEGAL, and W . P . HURLBUT From The Rockefeller University, New York 10021, and the Biophysics Laboratory, Veterans Administration Medical Center, New York 10010 Procedures are described for analyzing shot noise and determining the waveform, w(t), mean amplitude, (h), and mean rate of occurrence, (r), of the shots under a variety of nonideal conditions that include : (a) slow, spurious changes in the mean, (b) nonstationary shot rates, (c) nonuniform distribution of shot amplitudes, and (d) nonlinear summation of the shots . The procedures are based upon Rice's (1944 . Bell Telephone System journal. 23 : 282332) extension of Campbell's theorem to the second (variance), X2, third (skew), X3, and fourth, X4, semi-invariants (cumulants) of the noise . It is shown that the spectra of X2 and X3 of nonstationary shot noise contain a set of components that are proportional to (r) and arise from w(t), and a set of components that are independent of (r) and arise from the temporal variations in r(t) . Since the latter components are additive and are limited by the bandwidth of r(t), they can be removed by appropriate filters ; then (r) and (h) can be determined from the X2 and X3 of the filtered noise . We also show that a factor related to the ratio (X3)2/(\2)(X4) monitors the spread in the distribution of shot amplitudes and can be used to correct the estimates of (r) and (h) for the effects of that spread, if the shape of the distribution is known and if r(t) is stationary . The accuracy of the measurements of X4 is assessed and corrections for the effects of nonlinear summation of A2, X3, and X4 are derived . The procedures give valid results when they are used to analyze shot noise produced by the (linear) summation of simulated miniature endplate potentials, which are generated either at nonstationary rates or with a distribution of amplitudes . ABSTRACT
INTRODUCTION
Quantal secretion at neuromuscular junctions is traditionally studied by measuring the amplitudes of indirectly evoked endplate-potentials (EPPs) or currents Dr . Fesce's present address is Dept. of Medical Pharmacology and Center for the Study of Peripheral Neuropathies and Neuromuscular Diseases, University of Milan, CNR Center of Cytopharmacology, 20129 Milan, Italy . J.
GEN . PHYSIOL.
© The Rockefeller University Press - 0022-1295/86/07/0025/33 $1 .00
Volume 88 July 1986
25-57
25
26
THE JOURNAL OF GENERAL PHYSIOLOGY " VOLUME 88 - 1986
and by measuring the spontaneous rate of occurrence of miniature endplate potentials (MEPPs) or currents . Continuous measurements of the MEPP rate and amplitude during prolonged periods of intense secretion provide data that can be used to test the vesicle hypothesis and to characterize the kinetic properties of processes, such as vesicle recycling, that sustain high rates of quantal release. Very high MEPP rates cannot be measured by simple counting techniques because the individual MEPPs cannot be recognized, and other means must therefore be used . Since the MEPPs recorded at a single junction are similar in amplitude and time course, their mean rate and amplitude can, in principle, be estimated by applying the classical procedure of shot noise analysis to the fluctuations in the membrane potential of the endplate of a vigorously secreting junction . This procedure has been used to analyze shot noise in a number of physiological preparations (Katz and Miledi, 1972 ; Anderson and Stevens, 1973 ; Wong and Knight, 1980 ; Finger and Stettmeier, 1981). Ideal "shot noise" is a fluctuating signal produced by the linear summation of uniform elementary events (shots) that occur randomly at a constant mean rate . The standard theory of the statistics of shot noise states that its power spectrum is shaped like the spectrum of the shot waveform and that the mean, X,, and variance, X2, of the noise are related to the mean shot rate, (r), amplitude, h, and waveform, w(t), by the equations X, = (r)h fw(t)dt and a2 = (r)h 2 fw 2 (t)dt (Campbell, 1909 ; Rice, 1944). However, serious errors occur when these equations are applied to the neuromuscular junction, where the shots are the MEPPs, because (a) the mean membrane potential of the endplate is affected by many factors in addition to the summation of MEPPs, (b) the MEPPs may not occur at a stationary random rate, (c) MEPPs are not uniform in amplitude, and (d) MEPPs do not sum linearly (Martin, 1955). We show that these errors can be avoided or corrected for if the analysis is based upon the higher semi-invariants of appropriately filtered records of the noise. Courtney (1978) has previously pointed out the usefulness of the higher moments of the fluctuations in the amplitude of the EPP to the statistical analysis of quantal secretion . Campbell's theorem can be extended to higher semi-invariants (cumulants) of the noise according to the general relation (Rice, 1944) X" = (r)h"I", where I = f [w(t)]"dt, so that, in principle, any pair of semi-invariants can be used to compute (r) and h . If X3 (skew) and X2 are used, then errors arising from slow, spurious changes in membrane potential are avoided (because the mean is not used) and errors arising from nonlinear summation are greatly reduced. Segal et al . (1985) used X3 and X2 to estimate (r) and h at La 3+ -stimulated frog neuromuscular junctions, where this rate is quasi-stationary, changing so slowly that the power spectrum of the fluctuations has the shape of the spectrum of the MEPP waveform . However, under many experimental conditions (e .g ., in hypertonic solutions or in the presence of black widow spider venom [BWSV]), the power spectra of the endplate noise differ from the spectrum of the MEPP waveform in that the former are not white (flat) at low frequencies . These extra lowfrequency components in the noise spectrum could arise from a variety of phenomena that include nonstationary MEPP rates, correlations among MEPP waveform parameters, or extraneous sources of noise; they must be eliminated
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
27
or corrected for before Xs and a 2 can be used to calculate (r) and h . This article derives the frequency composition of a 2 (power spectrum) and a y (skew "bispectrum") of nonstationary shot noise and demonstrates that the effects of this phenomenon on these spectra can be virtually eliminated by appropriately filtering the data before analysis, as long as the deviations of the noise spectrum from the expected stationary shape are confined to a limited region of the spectrum . The semi-invariants of the filtered records can then be used to compute (r) and h . The power spectrum of nonstationary shot noise has been derived before by others (Rice, 1944 ; Schick, 1974 ; Sigworth, 1981), but the extension to the skew bispectrum is new. We also show that the ratio R = (X3/I8)2/(X2/I2)(X4/I4), which monitors the spread in the distribution of the shot amplitudes, can be used to correct the estimates of (r) and (h) for the effects of that spread . The errors in the measurement of X4 (fourth semi-invariant) are treated theoretically, and the corrections for the second-order effects of nonlinear summation on the Xn ' s are derived. The procedures are tested by analyzing nonideal shot noise generated by computer simulation, and they give valid results. In the companion article (Fesce et al ., 1986), the procedures are used to measure the MEPP rate and amplitude at neuromuscular junctions treated with BWSV, and to correct previous results obtained with La s+ . A synopsis of the derivation of the spectral components of a 2 and a s has been published previously (Fesce, 1986). THEORY
Nonstationary Shot Rate
Frequency composition ofvariance and skew of a f uctuating signal. Let a signal, V(t), be recorded over the period of time [t = 0, T] . The finite Fourier transform of V(t) is defined as : T
-2arint/T v(n) = f V(t) e dt
{n integer) .
(1)
From the theory of Fourier series, it follows that : V(t)
=
00 1 v(n) e2zint/T = ~ ) + 1 T n-_T T
= (V) + v(t),
v(m) e2,rimt/T
{m 0 01
(2)
where v(0)/T = (V) is the mean of V(t) over T, and v(t) is the departure of V(t) from (V) . The variance of V(t), (v 2(t)), can also be expressed as a series in v(n) : (V2(t)) ~ T 1
fT
fT
v2(t)dt
dt v(n)v(m) e2r=(n+m)t/T n=-oo - Ts m=-m
In 00,m00}.
28
THE JOURNAL OF GENERAL PHYSIOLOGY " VOLUME
88 - 1986
This equation can be evaluated by reversing the order of summation and integration ; the integrals with m 0 -n equal zero, whereas those with m = -n equal T. Therefore, (v2(t)) ~'
v(n)v( -n) = 12 v2(n) = 22 Iv(n)I 2, T12 n =_~ T n=_ T n =,
(3)
where the last step can be taken because v(-n) is the complex conjugate of ,6(n) . This result indicates that each value of n (n > 0) contributes a factor 21v(n) 12 IT2 to the variance . Since 1 /T is the interval between successive Fourier frequencies, the factor 2Iv(n)I 2/T is interpreted as the power density per unit bandwidth, G,(n), at the frequency f = n/T. When expressed in these terms, Eq . 3 becomes: (v 2(t)) =
T
Gv(n) . E Iv(n)I2/T = T E~ i
( 3 a)
This is the classical expression equating the variance of a fluctuating signal to the area under its power spectrum (Bendat and Piersol, 1971). A similar procedure applied to the skew of V(t) gives: (v3(t))=
T3
v(n)v(m)v(-n - m)
n=-oo
= 13
v3(n, m)
{n, m 54 0} .
The right-hand side of Eq . 4 gives the frequency composition of the skew, and it can be rewritten using only positive values for n and m. Notice that vs(n, m), unlike v2(n), is a complex quantity and that v3(-n, -m) is the complex conjugate of v3(n, m) . Therefore, iv 3(n, m) + v3(-n, -m) = 2Real[v3(n, m)], and v3(n, -m) + v3(-n, m) = 2Re[v3(n, m - n)] = 2Re[v3(m, n - m)]. If now we consider only positive values [ 1, oo] for n, m, n - m, and m - n, we have : (u3 (t)) = 62 T
tt=~
m=1
Re[v3(n, m)]/T .
(4a)
By analogy with the power density, the factor 6Re[v3(n, m)]/T is defined as the skew "bispectral" density, Bv(n, m), and we have : (v 3(t)) = 12 T
Bv(n, m) .
(4b)
The factor 1 /T 2 in Eq . 4b is the product of two frequency intervals between successive Fourier components ; the surface generated when B (n, m) is plotted against n and m is called the skew bispectrum and the volume under it equals the skew (Subba Rao and Gabr, 1984). Examples of bispectra are given in Fig. 7. Frequency composition of shot noise . Let V(t) be produced by the linear summation of randomly occurring elementary events, each of amplitude h and waveform w(t), which is limited in time such that w(t < 0) = w(t > -r) = 0 . If a
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
29
record extends over the interval [t = 0, T], then only those events occurring between the times [-r] and [T] will contribute to it. If K is the number of such events, then V(t) is given by: x V(t) _ E hw(t - B% where j is an arbitrary event index (not implying time sequence), and B is a event. random variable, representing the time of occurrence of the If the process is Poissonian, the individual events occur independently and the probability density functions of the B 's are all equal . Let p(t)dt be the probability for each of the K events to occur in the infinitesimal time interval between t and t + dt . If an ensemble of records is available, each with the same time course for p(t), then the expected value of V(t) for the ensemble is: j
jt'
j
E[V(t)]=E
[
K
E hw(t-B) =hEKEwt-B
hEK
T
t'wt-t'dt' .
Since E[K]p(t') is the expected rate of occurrence of events, r(t'), at time t', we have : EVt
h
T rt wt-t dt
h
f
rt-uwudu,
where u = t - t', and the change in the limits of integration is justified because all the values of t' for which w(t - t') 0 0 are included in the first integral . The ensemble expected value of the time average of V(t), (V), is computed by averaging E[V(t)] over time : 1-7[(V)] =
fT T
dt
f
0o
= h
f
r(t - u)w(u)du I
w(u)du
7,
T v
f
u
r(t')dt'
{t' = t - u} .
The integrands in Eq . 7 are zero when u > T, since w(u > T) = 0. Therefore, if T >> T, and if r(t') does not change rapidly and progressively (see Appendix), then changing the limits of integration of t' to [0, T] causes negligible errors, and we can write : w fT h w(u)du T r(t')dt' = (r)hl,, (7a) E[(V)] f
where (r) is the average value of r(t) over the interval [0, T], and 1, is the integral of w(t) . Note that the expected value of the mean signal is proportional to the mean rate and is independent of temporal variations in r(t) . The expected value of v(n) is, from Eqs. I and 6:
30
THE JOURNAL OF GENERAL PHYSIOLOGY " VOLUME
F.[v(n)] = f =h
T
.fo
Ch f
r(t - u)w(u)du
- w(u)du
J
T-u
e 2ainu/T fu
88 " 1986
dt e -2 rint/T r(t')dt e-taint'/T
(8)
it' =t- u,n00} . Neglecting errors at the edges and changing the limits of integration of t' to [0, T] gives: E[v(n)] ^' h f w(u)du e-2ainu/T
T
fo
r(t')dt' e -taint'/T = hzvo(n)r"(n),
(8a)
where z'a(n) and ~(n) are the finite Fourier transforms of w(t) and r(t), respectively
(Eq. 1).
Frequency composition of the variance of shot noise . The average power density in the bandwidth 1/T centered on the frequency n/T is, from Eq . 3a, G (n) _ 2v 2(n)/T = 2/T dt ds V(t)V(s) e2asn(t-s)/T. 1 n order to compute its expected value, we need the expected value of V(t)V(s) :
f f
K lI E[V(t)V(s)] = E [i i h2w(t - Bj)w(s - Bk) .
k
(9)
J
The K 2 terms of this product of two series can be grouped into two sets : one containing K terms with j = k and the other containing K(K - 1) terms with j54k : E[V(t)V(s)] = h 2E[K]E[w(t - Oj)w(s - Bj)] + h 2E[K(K - 1)] lk 0j) .
.E[w(t - Oj)w(s - Bk)]
(9a)
Now: T
E[w(t - Oj)w(s - B;] = f p(t')w(t - t')w(s - t')dt' ; T
(96)
and, since j and k are independent, E[w(t - Oj)w(s - Bk)] = E[w(t - Bj)]E[w(s - Bk)] = Therefore :
E[V(t)V(s)]
= E[K]h2
T
fT
f
T T
(9c)
p(t')w(t - t')dt' fT p(s')w(s - s')ds' . T
p(t')w(t - t')w(s - t')dt'
+ F,[K(K - 1)]h2
p(t')w(t - t')dt'
T
f,
p(s')w(s - s')ds' .
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
31
If the events occur independently, then the distribution of the K's is Poissonian so that E[K(K - 1)] = (E[K]) 2. Since E[K]p(t') = r(t'), then E[K(K - 1)]p(t') . p(s') = r(t')r(s'), and we have : EVtVs
h`
rt wt-t ws-t dt
Tf,
T
+ h2
f,
r(t')w(t - t')dt'
T
r(s')w(s - s')ds'
ft
(10)
= h2 f r(t')w(t - t')w(s - t')dt' + E[V(t)] E[V(s)], where the limits of the first integral have been extended to integrands are zero for t' < -T or t' > T. If Eq. 10 is substituted into Eq. 9, we get: T
±00,
since the
T
E[v2(n)] = f dt f ds E[V(t)V(s)] e2*in(t-s)/T T =
f
dt T
"
f 0
T
fo
ds E[V(t)]E[V(s)] e2sin(t-s)/T + h 2 f.' r(t')dt' T
dt
f 0
ds h2w(t - t' )w(S
- t')
(11)
e2xin(t-s)/T .
When edge errors are neglected, the first term of Eq. 11 gives: E[v(n)] E[v(-n)] = h2zo(n)zvu(-n)r"(n)i(-n)
=
h2za 2 (n)r2 (n).
(11a)
The second term can be evaluated over three different regions of t' : (a) t' < -T or t'>T,(b)020 Hz .
FESCE ET AL.
Fluctuation Analysis of Nonideal Shot Noise
47
computed from the mean and variance, but all errors decreased as the filter RC decreased ; when RC = 1 ms, the random errors were reduced to -10% (compare with stationary conditions in Fig . 3A) . Similar results were obtained when r(t) was varied in a continuous random manner . The fractional errors (mean ± SD) in the estimates obtained from mean A
___________________
+0.1 -------------i
-0 .1-____________
B
-0.1 ____
o
f 4-------------
.________
___-----
_
______________
-------------___1l
. ;_____1Y-----------------
C -0 .1 -------------- _ _________________
s ____________________________ --------------------------------------------
N
Z
__ 6______
+0 .1-------------_ _o,,_-------
_______________
T --------
______I_j.f______________ -
-
__________________________
E
+0 .1 -------------- _ .----------------0.1 ______________ _
RC
10
.________________
5
3 .5
------------------________
2
___________________
1
0.5
ms
Effect of filter time constant on the estimates of (r) and h from nonstationary records : superimposed volleys with the parameters reported in the text (A-E) . Fractional errors (mean ± SD, 10 independently simulated records) in the estimates of MEPP rate (triangles) and amplitude (squares), from mean and variance (solid symbols), or from variance and skew (open symbols) . Variance and skew are those of data filtered with the time constant indicated on the abscissa, while the mean was computed before filtering. FIGURE 6 .
(unfiltered) and variance were : -43 ± 2% for (r) and 77 ± 7% for h (no filter) ; -12 ± 2% for (r) and 14 ± 3% for h (RC - 5 ms); -2.5 ± 2% for (r) and 2.7 ± 2% for h (RC = 1 ms). The errors in the estimates obtained from variance and skew were: -13 ± 23% for (r) and 47 ± 21% for h (no filter) ; -5 ± 10.5% for (r) and 10 ± 6% for h (RC ^- 5 ms); -0.8 ± 8.5% for (r) and 2 ± 5% for h
FIGURE 7 . (Left) Bispectral density of the skew . To obtain this plot, a record was divided into 12 sections of 2,048 points, Fourier transforms were performed on each section, and the bispectral densities were averaged over groups of neighboring frequencies (15 x 15 for low frequencies, 21 x 21 for intermediate, and 71 x 71 for high frequencies), and also over the 12 sections . Bispectral densities were further smoothed once over neighboring points (each point is half its original value plus 1/8 of each of the four adjacent points) . The two horizontal axes (X-Y) are Fourier frequencies on logarithmic scales, and the vertical axis is the bispectral density on a logarithmic scale . Positive values that would fall below the X-Y plane (i .e ., densities between 0 and 1) are not plotted ; negative values are plotted below the X-Y plane (negative logarithm of the absolute value) . Top panel : simulated stationary record (same as Fig . 2) ; second panel : analytical bispectrum of the simulated MEPP waveform ; third panel : simulated nonstationary record (same as Figs . 4B and 5) ; note the deviations from the analytical shape ; fourth panel : same record as in the third panel, filtered through a 1-ms RC filter; fifth panel : analytical bispectrum of the simulated MEPP waveform, filtered through a 1-ms RC filter . Differences are still apparent beween the densities in the two bottom panels, but their contributions to the skew are negligible (right) . (Right) Integrals of the skew bispectra in the corresponding panel on the left . Horizontal axes as on the left . Vertical axes are linear (arbitrary units) . Notice the good agreement between the two top panels (unfiltered simulated stationary record and analytical plot for stationary data), and between the two bottom panels (filtered simulated nonstationary record and analytical plot of filtered stationary data) .
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
49
(RC = 1 ms). The random errors in the estimates of (r) and h were reduced further when the moments of three independent filtered records were averaged . This averaging procedure is equivalent to increasing by threefold the duration of a data-collection interval, and the resultant reduction in the random errors agrees with the predictions of the theory (Appendix) . The skew bispectra of the simulated records shown in Fig. 7 support the theoretical prediction that the contributions of nonstationarities to the skew are limited to the bandwidth of r(t) and are effectively removed by appropriate highpass filtering.
FIGURE 8 . Histograms of the distribution of h as determined from 'y functions with various values of ,B and y. Ordinates : number of events. Abscissae: value of h (arbitrary units) . For one set of histogram (5,000 events/histogram), ,B = 1 and (h) (indicated by vertical arrows) decreased as y was decreased from 10 to 1 . For the other set (10,000 events/histogram), both 3 and ti were varied such that (h) was constant .
The Dependence of X4 upon (r) and the Effects ofDistributed MEPP Amplitudes When r(t) is stationary, the X 's are proportional to (r) (Rice, 1944). The random error in their estimates should increase with (r) (for n > 2), and should decrease
with high-pass filtering (Segal et al ., 1985). We examined the behavior of X4 of simulated stationary records . For unfiltered data, the slope of the log-log plot of X4 vs. (r) was 0 .998 (correlation coefficient, 0.993) over the range from 10/s to 1,000/s; the coefficient of variation (CV = mean/SD) of the estimates (15 at each (r)) rose from 0 .07 to 0.41 over this range, and it exceeded 1 at 3,000/s. For filtered data (r = 1 ms), the slope was 0.968 (correlation coefficient, 0.983) over the range from 10/s to 10,000/s, and CV increased from 0 .11 to 0.49 over this more extensive range. Thus, a4 behaves as expected when (r) is stationary, and it can be estimated with reasonable accuracy at rates approaching 10,000/s, if the data are appropriately filtered .
50
THE JOURNAL OF GENERAL PHYSIOLOGY - VOLUME
88 - 1986
We also simulated stationary ((r) = 932/s) records of filtered MEPPs whose h factors were distributed in accordance with y distributions . In one set of simulations, the ,# parameter of the distributions was fixed at 1, while y varied from oo (all h's equal) to 1 . In a second set, both (3 and y were varied such that (h) remained constant [(h) = K(y + 1)/a, where K is a scale factor). Fig. 8 shows the resulting distributions of h. The distribution is approximately Gaussian (CV -0 .3) when y = 10 and is highly skewed when y = 1 . Table I shows the values of R and the apparent (r) and (h) calculated from X2 and X3 of the simulated records . As y falls, the apparent (r) falls and the apparent (h) rises; they approach values of -1/3 and 2, respectively, of the applied values when y nears 1 . The table also shows the results obtained when the measured values of R are used to correct the apparent (r) and (h) for the distribution of h . The corrected values lie within 9% of the applied ones when y
INTRODUCTION
Quantal secretion at neuromuscular junctions is traditionally studied by measuring the amplitudes of indirectly evoked endplate-potentials (EPPs) or currents Dr . Fesce's present address is Dept. of Medical Pharmacology and Center for the Study of Peripheral Neuropathies and Neuromuscular Diseases, University of Milan, CNR Center of Cytopharmacology, 20129 Milan, Italy . J.
GEN . PHYSIOL.
© The Rockefeller University Press - 0022-1295/86/07/0025/33 $1 .00
Volume 88 July 1986
25-57
25
26
THE JOURNAL OF GENERAL PHYSIOLOGY " VOLUME 88 - 1986
and by measuring the spontaneous rate of occurrence of miniature endplate potentials (MEPPs) or currents . Continuous measurements of the MEPP rate and amplitude during prolonged periods of intense secretion provide data that can be used to test the vesicle hypothesis and to characterize the kinetic properties of processes, such as vesicle recycling, that sustain high rates of quantal release. Very high MEPP rates cannot be measured by simple counting techniques because the individual MEPPs cannot be recognized, and other means must therefore be used . Since the MEPPs recorded at a single junction are similar in amplitude and time course, their mean rate and amplitude can, in principle, be estimated by applying the classical procedure of shot noise analysis to the fluctuations in the membrane potential of the endplate of a vigorously secreting junction . This procedure has been used to analyze shot noise in a number of physiological preparations (Katz and Miledi, 1972 ; Anderson and Stevens, 1973 ; Wong and Knight, 1980 ; Finger and Stettmeier, 1981). Ideal "shot noise" is a fluctuating signal produced by the linear summation of uniform elementary events (shots) that occur randomly at a constant mean rate . The standard theory of the statistics of shot noise states that its power spectrum is shaped like the spectrum of the shot waveform and that the mean, X,, and variance, X2, of the noise are related to the mean shot rate, (r), amplitude, h, and waveform, w(t), by the equations X, = (r)h fw(t)dt and a2 = (r)h 2 fw 2 (t)dt (Campbell, 1909 ; Rice, 1944). However, serious errors occur when these equations are applied to the neuromuscular junction, where the shots are the MEPPs, because (a) the mean membrane potential of the endplate is affected by many factors in addition to the summation of MEPPs, (b) the MEPPs may not occur at a stationary random rate, (c) MEPPs are not uniform in amplitude, and (d) MEPPs do not sum linearly (Martin, 1955). We show that these errors can be avoided or corrected for if the analysis is based upon the higher semi-invariants of appropriately filtered records of the noise. Courtney (1978) has previously pointed out the usefulness of the higher moments of the fluctuations in the amplitude of the EPP to the statistical analysis of quantal secretion . Campbell's theorem can be extended to higher semi-invariants (cumulants) of the noise according to the general relation (Rice, 1944) X" = (r)h"I", where I = f [w(t)]"dt, so that, in principle, any pair of semi-invariants can be used to compute (r) and h . If X3 (skew) and X2 are used, then errors arising from slow, spurious changes in membrane potential are avoided (because the mean is not used) and errors arising from nonlinear summation are greatly reduced. Segal et al . (1985) used X3 and X2 to estimate (r) and h at La 3+ -stimulated frog neuromuscular junctions, where this rate is quasi-stationary, changing so slowly that the power spectrum of the fluctuations has the shape of the spectrum of the MEPP waveform . However, under many experimental conditions (e .g ., in hypertonic solutions or in the presence of black widow spider venom [BWSV]), the power spectra of the endplate noise differ from the spectrum of the MEPP waveform in that the former are not white (flat) at low frequencies . These extra lowfrequency components in the noise spectrum could arise from a variety of phenomena that include nonstationary MEPP rates, correlations among MEPP waveform parameters, or extraneous sources of noise; they must be eliminated
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
27
or corrected for before Xs and a 2 can be used to calculate (r) and h . This article derives the frequency composition of a 2 (power spectrum) and a y (skew "bispectrum") of nonstationary shot noise and demonstrates that the effects of this phenomenon on these spectra can be virtually eliminated by appropriately filtering the data before analysis, as long as the deviations of the noise spectrum from the expected stationary shape are confined to a limited region of the spectrum . The semi-invariants of the filtered records can then be used to compute (r) and h . The power spectrum of nonstationary shot noise has been derived before by others (Rice, 1944 ; Schick, 1974 ; Sigworth, 1981), but the extension to the skew bispectrum is new. We also show that the ratio R = (X3/I8)2/(X2/I2)(X4/I4), which monitors the spread in the distribution of the shot amplitudes, can be used to correct the estimates of (r) and (h) for the effects of that spread . The errors in the measurement of X4 (fourth semi-invariant) are treated theoretically, and the corrections for the second-order effects of nonlinear summation on the Xn ' s are derived. The procedures are tested by analyzing nonideal shot noise generated by computer simulation, and they give valid results. In the companion article (Fesce et al ., 1986), the procedures are used to measure the MEPP rate and amplitude at neuromuscular junctions treated with BWSV, and to correct previous results obtained with La s+ . A synopsis of the derivation of the spectral components of a 2 and a s has been published previously (Fesce, 1986). THEORY
Nonstationary Shot Rate
Frequency composition ofvariance and skew of a f uctuating signal. Let a signal, V(t), be recorded over the period of time [t = 0, T] . The finite Fourier transform of V(t) is defined as : T
-2arint/T v(n) = f V(t) e dt
{n integer) .
(1)
From the theory of Fourier series, it follows that : V(t)
=
00 1 v(n) e2zint/T = ~ ) + 1 T n-_T T
= (V) + v(t),
v(m) e2,rimt/T
{m 0 01
(2)
where v(0)/T = (V) is the mean of V(t) over T, and v(t) is the departure of V(t) from (V) . The variance of V(t), (v 2(t)), can also be expressed as a series in v(n) : (V2(t)) ~ T 1
fT
fT
v2(t)dt
dt v(n)v(m) e2r=(n+m)t/T n=-oo - Ts m=-m
In 00,m00}.
28
THE JOURNAL OF GENERAL PHYSIOLOGY " VOLUME
88 - 1986
This equation can be evaluated by reversing the order of summation and integration ; the integrals with m 0 -n equal zero, whereas those with m = -n equal T. Therefore, (v2(t)) ~'
v(n)v( -n) = 12 v2(n) = 22 Iv(n)I 2, T12 n =_~ T n=_ T n =,
(3)
where the last step can be taken because v(-n) is the complex conjugate of ,6(n) . This result indicates that each value of n (n > 0) contributes a factor 21v(n) 12 IT2 to the variance . Since 1 /T is the interval between successive Fourier frequencies, the factor 2Iv(n)I 2/T is interpreted as the power density per unit bandwidth, G,(n), at the frequency f = n/T. When expressed in these terms, Eq . 3 becomes: (v 2(t)) =
T
Gv(n) . E Iv(n)I2/T = T E~ i
( 3 a)
This is the classical expression equating the variance of a fluctuating signal to the area under its power spectrum (Bendat and Piersol, 1971). A similar procedure applied to the skew of V(t) gives: (v3(t))=
T3
v(n)v(m)v(-n - m)
n=-oo
= 13
v3(n, m)
{n, m 54 0} .
The right-hand side of Eq . 4 gives the frequency composition of the skew, and it can be rewritten using only positive values for n and m. Notice that vs(n, m), unlike v2(n), is a complex quantity and that v3(-n, -m) is the complex conjugate of v3(n, m) . Therefore, iv 3(n, m) + v3(-n, -m) = 2Real[v3(n, m)], and v3(n, -m) + v3(-n, m) = 2Re[v3(n, m - n)] = 2Re[v3(m, n - m)]. If now we consider only positive values [ 1, oo] for n, m, n - m, and m - n, we have : (u3 (t)) = 62 T
tt=~
m=1
Re[v3(n, m)]/T .
(4a)
By analogy with the power density, the factor 6Re[v3(n, m)]/T is defined as the skew "bispectral" density, Bv(n, m), and we have : (v 3(t)) = 12 T
Bv(n, m) .
(4b)
The factor 1 /T 2 in Eq . 4b is the product of two frequency intervals between successive Fourier components ; the surface generated when B (n, m) is plotted against n and m is called the skew bispectrum and the volume under it equals the skew (Subba Rao and Gabr, 1984). Examples of bispectra are given in Fig. 7. Frequency composition of shot noise . Let V(t) be produced by the linear summation of randomly occurring elementary events, each of amplitude h and waveform w(t), which is limited in time such that w(t < 0) = w(t > -r) = 0 . If a
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
29
record extends over the interval [t = 0, T], then only those events occurring between the times [-r] and [T] will contribute to it. If K is the number of such events, then V(t) is given by: x V(t) _ E hw(t - B% where j is an arbitrary event index (not implying time sequence), and B is a event. random variable, representing the time of occurrence of the If the process is Poissonian, the individual events occur independently and the probability density functions of the B 's are all equal . Let p(t)dt be the probability for each of the K events to occur in the infinitesimal time interval between t and t + dt . If an ensemble of records is available, each with the same time course for p(t), then the expected value of V(t) for the ensemble is: j
jt'
j
E[V(t)]=E
[
K
E hw(t-B) =hEKEwt-B
hEK
T
t'wt-t'dt' .
Since E[K]p(t') is the expected rate of occurrence of events, r(t'), at time t', we have : EVt
h
T rt wt-t dt
h
f
rt-uwudu,
where u = t - t', and the change in the limits of integration is justified because all the values of t' for which w(t - t') 0 0 are included in the first integral . The ensemble expected value of the time average of V(t), (V), is computed by averaging E[V(t)] over time : 1-7[(V)] =
fT T
dt
f
0o
= h
f
r(t - u)w(u)du I
w(u)du
7,
T v
f
u
r(t')dt'
{t' = t - u} .
The integrands in Eq . 7 are zero when u > T, since w(u > T) = 0. Therefore, if T >> T, and if r(t') does not change rapidly and progressively (see Appendix), then changing the limits of integration of t' to [0, T] causes negligible errors, and we can write : w fT h w(u)du T r(t')dt' = (r)hl,, (7a) E[(V)] f
where (r) is the average value of r(t) over the interval [0, T], and 1, is the integral of w(t) . Note that the expected value of the mean signal is proportional to the mean rate and is independent of temporal variations in r(t) . The expected value of v(n) is, from Eqs. I and 6:
30
THE JOURNAL OF GENERAL PHYSIOLOGY " VOLUME
F.[v(n)] = f =h
T
.fo
Ch f
r(t - u)w(u)du
- w(u)du
J
T-u
e 2ainu/T fu
88 " 1986
dt e -2 rint/T r(t')dt e-taint'/T
(8)
it' =t- u,n00} . Neglecting errors at the edges and changing the limits of integration of t' to [0, T] gives: E[v(n)] ^' h f w(u)du e-2ainu/T
T
fo
r(t')dt' e -taint'/T = hzvo(n)r"(n),
(8a)
where z'a(n) and ~(n) are the finite Fourier transforms of w(t) and r(t), respectively
(Eq. 1).
Frequency composition of the variance of shot noise . The average power density in the bandwidth 1/T centered on the frequency n/T is, from Eq . 3a, G (n) _ 2v 2(n)/T = 2/T dt ds V(t)V(s) e2asn(t-s)/T. 1 n order to compute its expected value, we need the expected value of V(t)V(s) :
f f
K lI E[V(t)V(s)] = E [i i h2w(t - Bj)w(s - Bk) .
k
(9)
J
The K 2 terms of this product of two series can be grouped into two sets : one containing K terms with j = k and the other containing K(K - 1) terms with j54k : E[V(t)V(s)] = h 2E[K]E[w(t - Oj)w(s - Bj)] + h 2E[K(K - 1)] lk 0j) .
.E[w(t - Oj)w(s - Bk)]
(9a)
Now: T
E[w(t - Oj)w(s - B;] = f p(t')w(t - t')w(s - t')dt' ; T
(96)
and, since j and k are independent, E[w(t - Oj)w(s - Bk)] = E[w(t - Bj)]E[w(s - Bk)] = Therefore :
E[V(t)V(s)]
= E[K]h2
T
fT
f
T T
(9c)
p(t')w(t - t')dt' fT p(s')w(s - s')ds' . T
p(t')w(t - t')w(s - t')dt'
+ F,[K(K - 1)]h2
p(t')w(t - t')dt'
T
f,
p(s')w(s - s')ds' .
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
31
If the events occur independently, then the distribution of the K's is Poissonian so that E[K(K - 1)] = (E[K]) 2. Since E[K]p(t') = r(t'), then E[K(K - 1)]p(t') . p(s') = r(t')r(s'), and we have : EVtVs
h`
rt wt-t ws-t dt
Tf,
T
+ h2
f,
r(t')w(t - t')dt'
T
r(s')w(s - s')ds'
ft
(10)
= h2 f r(t')w(t - t')w(s - t')dt' + E[V(t)] E[V(s)], where the limits of the first integral have been extended to integrands are zero for t' < -T or t' > T. If Eq. 10 is substituted into Eq. 9, we get: T
±00,
since the
T
E[v2(n)] = f dt f ds E[V(t)V(s)] e2*in(t-s)/T T =
f
dt T
"
f 0
T
fo
ds E[V(t)]E[V(s)] e2sin(t-s)/T + h 2 f.' r(t')dt' T
dt
f 0
ds h2w(t - t' )w(S
- t')
(11)
e2xin(t-s)/T .
When edge errors are neglected, the first term of Eq. 11 gives: E[v(n)] E[v(-n)] = h2zo(n)zvu(-n)r"(n)i(-n)
=
h2za 2 (n)r2 (n).
(11a)
The second term can be evaluated over three different regions of t' : (a) t' < -T or t'>T,(b)020 Hz .
FESCE ET AL.
Fluctuation Analysis of Nonideal Shot Noise
47
computed from the mean and variance, but all errors decreased as the filter RC decreased ; when RC = 1 ms, the random errors were reduced to -10% (compare with stationary conditions in Fig . 3A) . Similar results were obtained when r(t) was varied in a continuous random manner . The fractional errors (mean ± SD) in the estimates obtained from mean A
___________________
+0.1 -------------i
-0 .1-____________
B
-0.1 ____
o
f 4-------------
.________
___-----
_
______________
-------------___1l
. ;_____1Y-----------------
C -0 .1 -------------- _ _________________
s ____________________________ --------------------------------------------
N
Z
__ 6______
+0 .1-------------_ _o,,_-------
_______________
T --------
______I_j.f______________ -
-
__________________________
E
+0 .1 -------------- _ .----------------0.1 ______________ _
RC
10
.________________
5
3 .5
------------------________
2
___________________
1
0.5
ms
Effect of filter time constant on the estimates of (r) and h from nonstationary records : superimposed volleys with the parameters reported in the text (A-E) . Fractional errors (mean ± SD, 10 independently simulated records) in the estimates of MEPP rate (triangles) and amplitude (squares), from mean and variance (solid symbols), or from variance and skew (open symbols) . Variance and skew are those of data filtered with the time constant indicated on the abscissa, while the mean was computed before filtering. FIGURE 6 .
(unfiltered) and variance were : -43 ± 2% for (r) and 77 ± 7% for h (no filter) ; -12 ± 2% for (r) and 14 ± 3% for h (RC - 5 ms); -2.5 ± 2% for (r) and 2.7 ± 2% for h (RC = 1 ms). The errors in the estimates obtained from variance and skew were: -13 ± 23% for (r) and 47 ± 21% for h (no filter) ; -5 ± 10.5% for (r) and 10 ± 6% for h (RC ^- 5 ms); -0.8 ± 8.5% for (r) and 2 ± 5% for h
FIGURE 7 . (Left) Bispectral density of the skew . To obtain this plot, a record was divided into 12 sections of 2,048 points, Fourier transforms were performed on each section, and the bispectral densities were averaged over groups of neighboring frequencies (15 x 15 for low frequencies, 21 x 21 for intermediate, and 71 x 71 for high frequencies), and also over the 12 sections . Bispectral densities were further smoothed once over neighboring points (each point is half its original value plus 1/8 of each of the four adjacent points) . The two horizontal axes (X-Y) are Fourier frequencies on logarithmic scales, and the vertical axis is the bispectral density on a logarithmic scale . Positive values that would fall below the X-Y plane (i .e ., densities between 0 and 1) are not plotted ; negative values are plotted below the X-Y plane (negative logarithm of the absolute value) . Top panel : simulated stationary record (same as Fig . 2) ; second panel : analytical bispectrum of the simulated MEPP waveform ; third panel : simulated nonstationary record (same as Figs . 4B and 5) ; note the deviations from the analytical shape ; fourth panel : same record as in the third panel, filtered through a 1-ms RC filter; fifth panel : analytical bispectrum of the simulated MEPP waveform, filtered through a 1-ms RC filter . Differences are still apparent beween the densities in the two bottom panels, but their contributions to the skew are negligible (right) . (Right) Integrals of the skew bispectra in the corresponding panel on the left . Horizontal axes as on the left . Vertical axes are linear (arbitrary units) . Notice the good agreement between the two top panels (unfiltered simulated stationary record and analytical plot for stationary data), and between the two bottom panels (filtered simulated nonstationary record and analytical plot of filtered stationary data) .
FESCE ET AL .
Fluctuation Analysis of Nonideal Shot Noise
49
(RC = 1 ms). The random errors in the estimates of (r) and h were reduced further when the moments of three independent filtered records were averaged . This averaging procedure is equivalent to increasing by threefold the duration of a data-collection interval, and the resultant reduction in the random errors agrees with the predictions of the theory (Appendix) . The skew bispectra of the simulated records shown in Fig. 7 support the theoretical prediction that the contributions of nonstationarities to the skew are limited to the bandwidth of r(t) and are effectively removed by appropriate highpass filtering.
FIGURE 8 . Histograms of the distribution of h as determined from 'y functions with various values of ,B and y. Ordinates : number of events. Abscissae: value of h (arbitrary units) . For one set of histogram (5,000 events/histogram), ,B = 1 and (h) (indicated by vertical arrows) decreased as y was decreased from 10 to 1 . For the other set (10,000 events/histogram), both 3 and ti were varied such that (h) was constant .
The Dependence of X4 upon (r) and the Effects ofDistributed MEPP Amplitudes When r(t) is stationary, the X 's are proportional to (r) (Rice, 1944). The random error in their estimates should increase with (r) (for n > 2), and should decrease
with high-pass filtering (Segal et al ., 1985). We examined the behavior of X4 of simulated stationary records . For unfiltered data, the slope of the log-log plot of X4 vs. (r) was 0 .998 (correlation coefficient, 0.993) over the range from 10/s to 1,000/s; the coefficient of variation (CV = mean/SD) of the estimates (15 at each (r)) rose from 0 .07 to 0.41 over this range, and it exceeded 1 at 3,000/s. For filtered data (r = 1 ms), the slope was 0.968 (correlation coefficient, 0.983) over the range from 10/s to 10,000/s, and CV increased from 0 .11 to 0.49 over this more extensive range. Thus, a4 behaves as expected when (r) is stationary, and it can be estimated with reasonable accuracy at rates approaching 10,000/s, if the data are appropriately filtered .
50
THE JOURNAL OF GENERAL PHYSIOLOGY - VOLUME
88 - 1986
We also simulated stationary ((r) = 932/s) records of filtered MEPPs whose h factors were distributed in accordance with y distributions . In one set of simulations, the ,# parameter of the distributions was fixed at 1, while y varied from oo (all h's equal) to 1 . In a second set, both (3 and y were varied such that (h) remained constant [(h) = K(y + 1)/a, where K is a scale factor). Fig. 8 shows the resulting distributions of h. The distribution is approximately Gaussian (CV -0 .3) when y = 10 and is highly skewed when y = 1 . Table I shows the values of R and the apparent (r) and (h) calculated from X2 and X3 of the simulated records . As y falls, the apparent (r) falls and the apparent (h) rises; they approach values of -1/3 and 2, respectively, of the applied values when y nears 1 . The table also shows the results obtained when the measured values of R are used to correct the apparent (r) and (h) for the distribution of h . The corrected values lie within 9% of the applied ones when y
