INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS CALCULATION METHODS RESULTS DISCUSSION REFERENCES

Communication between neurons takes place via stochastic spike trains that reflect random synaptic conductance transients. Since synaptic transmission itself is random (Allen and Stevens, 1994; Hardingham and Larkman, 1998), the processing in neuronal systems unavoidably becomes stochastic (Tuckwell, 1988). Indeed, intracellular recordings in vivo have revealed strong stochastic membrane potential fluctuations (Calvin and Stevens, 1968; Stern et al., 1997; Pare et al., 1998). Thus, to understand the principles of neuronal information processing it is necessary to describe the influence stochastic signals have on the statistical properties of membrane potential and firing.

This problem was introduced several decades ago, when simplified integrate-and-fire models of neurons receiving stochastic inputs were analyzed (Gerstein and Mandelbrot, 1964). Later work steadily broadened the scope of analysis by introducing leaky integrate-and-fire models (Stein, 1965), by applying techniques for calculation of firing statistics (Gluss, 1967), and recently by analyzing nonstationary inputs (Burkitt and Clark, 1999). The analytical studies could explain a number of experimentally observed firing statistics (Gerstein and Mandelbrot, 1964; Treves et al., 1999). However, it has been claimed that simple models are not able to explain the high variability of cortical cell firing observed in vivo (Softky and Koch, 1993). It is often assumed that synaptic input consists of independent presynaptic spikes, yet a significant number of neurons in the visual system have spiking patterns that could not be described by Poisson statistics (Reich et al., 1998).

Traditionally, it was believed that neuronal firing rate carries the information for coding and decoding. Recent studies suggest that precise spike timing (Abeles and Prut, 1996; Bair and Koch, 1996) and correlation of firing between different neurons (Gray et al., 1989; Roelfsema et al., 1997) could be implicated in information processing as well. Although questions remain of how precise temporal coding might be (Shadlen and Newsome, 1994), there is growing acceptance that the temporal structure of firing is important (Vaadia et al., 1995; Riehle et al., 1997).

The temporal structure in the synaptic input current is induced by the finite decay time of the synaptic conductance and by the temporal correlation between input event times. The synaptic conductance decay time may vary from <1 ms in auditory neurons (Raman and Trussell, 1992) to several tens of milliseconds for NMDA receptor-mediated currents (Silver et al., 1992). The temporal correlation of the input events may range broadly in time scale and form, e.g., from exponential decay (Weliky and Katz, 1999; Brivanlou et al., 1998) to decaying oscillations (Gray et al., 1989; Roelfsema et al., 1997). The sources of these correlations could be electrical coupling (Brivanlou et al., 1998; Mann-Metzer and Yarom, 1999; Gibson et al., 1999), synaptic interaction, and/or shared input (Brivanlou et al., 1998) of presynaptic neurons.

In this study we apply the theory of random point processes (Stratonovich, 1963; van Kampen, 1992) to study the effects of synaptic input temporal structure on the firing statistics. Our theory, based on small-amplitude inputs, can be applied to describe a broad class of random point series. We assess the approximation's applicability for neurons with spike afterhyperpolarization by comparing our theoretical results with simulations. The membrane potential fluctuations are described accurately as Gaussian when the input event rate exceeds by an order of magnitude or so the reciprocal of the system's slowest time constant. For low-frequency firing the theory adequately estimates the firing rate. Our simulations show that slow temporal correlation in the synaptic input current can bring about high variability of firing by inducing bursts of spikes. Our theory could be used for the analysis of intracellularly recorded membrane potential fluctuations (Lampl et al., 1999; Azouz and Gray, 1999) in order to extract possible temporal structure of the synaptic input.

	THEORETICAL METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS CALCULATION METHODS RESULTS DISCUSSION REFERENCES

We analytically study a simplified one-compartment model neuron where spike generation is associated with a crossing of the threshold level by the fluctuating membrane potential. Our simplified system is governed by equations that describe the dynamics of input synaptic current, I, and membrane potential, V (as deviation from the resting potential):

(1a)

Here, _m is the membrane time constant (in ms), C_N is the neuron's capacitance (µF), and _I(t) is a random process describing the synaptic current that we express as:

(1b)

The function s_I(t  t_j) describes the unitary synaptic current time course, e.g., as an exponential, exp((t  t_j)/_s) or as an alpha function, (t  t_j)exp((t  t_j)/_s), where _s is the synaptic time constant; s_I is equal to 0 if t < t_j. The synaptic current's random amplitude is a_j, and t_j is the arrival time of the synaptic event. If we replace s_I(t  t_j) with a function, v(t  t_j), for the unitary synaptic potential which is the solution of Eq. 1a, we obtain the equation for the membrane potential as a stochastic process:

(1c)

A more realistic description of the synaptic input is to account for changes in the synaptic conductance G. In this case the equation for the model neuron changes to:

(2)

where E_s is the synaptic current's reversal potential (deviation from the resting membrane potential) and _G(t) is a process for synaptic conductance amplitude and is described as in Eq. 1b.

Our analysis is based on the characteristic functional for a continuous stochastic process (van Kampen, 1992)
where braces indicate ensemble averaging. From this we may obtain the correlation functions, k_n(t₁, ... , t_n), of the random membrane potential, _V, or synaptic current, _I. A set of these functions completely describes the stochastic process. The characteristic functional is a generalization of the multivariable characteristic function (u₁, ... , u_n) = exp(iu₁₁ + ··· + iu_nn) for a continuous stochastic process. By using functional derivatives, the characteristic functional allows for calculating correlation functions
and, thus, could be expressed as follows (Stratonovich, 1963; van Kampen, 1992):

(3)

Stochastic continuous processes _I(t) and _V(t) are generated by series of synaptic events viewed as random point processes. Assuming independence between arrival times t_j and of synaptic amplitudes a_j, the characteristic functional for the net synaptic current is

(4)

where w(a_j) is the probability density for synaptic amplitudes and braces indicate averaging with respect to the arrival times t_j. Let us denote the terms inside the braces by W[u(t_i)]:

(5)

Then the characteristic functional of the random point process acquires the simple form

(6)

In order to average over arrival times we must first specify the statistical properties of the time series of synaptic events. For a complete description of these input trains, we use the correlation functions, g_n(t₁, ... , t_n) (Stratonovich, 1963; van Kampen, 1992), for the trains as random point processes. These correlation functions are related to the probabilities, f_n(t₁, ... , t_n)dt₁ ... dt_n, of having one event in each of the time intervals: dt₁ around t₁, dt₂ around t₂, etc. The relation is similar to that between cumulants (semi-invariants) and moments of random variables (Risken, 1989). For example, the first function g₁(t₁) = f₁(t₁) is the expected rate of events at time t₁. The second function g₂(t₁, t₂) = f₂(t₁, t₂)  f₁(t₁)f₁(t₂) is the cross-correlation between presynaptic spikes similar to the joint peristimulus histogram (Aertsen et al., 1989) used to study neuronal interaction by correlating spike trains from neuron pairs (Aertsen et al., 1989; Vaadia et al., 1995). Thus, if g_n(t₁, ... , t_n) = 0 for n > 1, the process would be a non-homogeneous Poisson point process. Using the correlation functions, g_n(t₁, ... , t_n), we can perform the averaging over the arrival times of synaptic events (Stratonovich, 1963; van Kampen, 1992), and obtain the characteristic functional in the explicit form:

(7)

To calculate the integral W[u(t)] in Eq. 5 we notice that it is equal to the characteristic function, _a(z(t_j)), of the amplitudes of synaptic events, where the variable z(t_j) = ₀^T u(t)s_I(t  t_j)dt. Let us assume that the probability density for the synaptic amplitudes is either a Gaussian w_a(a) = exp((a_j  A)²/2²)/ (Wahl et al., 1997; Hardingham and Larkman, 1998) or an exponential w_a(a) = 1/ · exp(a_j/) (Matsui et al., 1998; Hardingham and Larkman, 1998), where A and  represent mean and  is standard deviation of the synaptic amplitude. Then, assuming that inputs are small, we expand with respect to z(t_j), obtaining

and
for the Gaussian and exponential probability densities, respectively. After substituting this expansion into Eq. 7 and changing the order of integration, the correlation functions k_n(t₁, ... , t_n), can be obtained by equating coefficients of like powers of u(t) in Eq. 3.

Below, we will write the mean, k₁(t₁), and the two-time correlation function, k₂(t₁, t₂), for the synaptic current; the corresponding expressions for the membrane potential can be obtained if v(t  t') is used instead of s_I(t  t'). Thus, for amplitudes with the Gaussian distribution the mean of the process and the two-time correlation function are given by:


For the synaptic current process with amplitudes obeying the exponential distribution:

(8)

(9)

In these equations, the integration upper limit, T, should be defined according to the principal of causality, i.e., the upper limit cannot exceed observation time t₁ or t₂. If the higher-order correlation functions are small compared to the first two, the synaptic current or membrane potential will have Gaussian probability distributions.

Notice that, according to Eq. 9, the auto-correlation has two contributions: the first term is defined mostly by the shape of the unitary synaptic current, while the second term includes the temporal correlation structure of the input trains. Suppose the input train is stationary, g₁ = const. Then, if the synaptic current and input train's temporal correlation decay exponentially, the influence of these two contributions are similar. Indeed, if synaptic current decays very fast and the correlation between inputs decays exponentially, then only the second term contributes and k₂(t₁, t₂) ~ exp(|t₁  t₂|/), where  is the characteristic time for the decay. However, if synaptic currents decay exponentially and the correlation between inputs is very small or decays very fast, the auto-correlation function has similar time dependency because only the first term is important. Thus, it is enough to study only the influence of synaptic current decay on the statistics of membrane potential or firing to gain an understanding of the qualitative effects.

	CALCULATION METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS CALCULATION METHODS RESULTS DISCUSSION REFERENCES

Simulations were done in order to check applicability of our analytical results, as approximations based on assuming small amplitude inputs. For the statistics of subthreshold membrane potential fluctuations, the firing of the model neuron was disallowed. Equations 1a and 2 (see Results) were solved numerically using the implicit trapezoidal scheme (Kloeden and Platen, 1992) with a time step of 0.005 ms. If the smaller time step of 0.001 ms was used, the statistics were the same. The occurrence times and synaptic amplitudes were generated using algorithms from Press et al., 1992. Although our theory handles both Gaussian and exponential densities for the synaptic current amplitudes (see above), we used only exponential distributions to evaluate the theoretical results.

For estimating neuronal firing rate we did not directly simulate synaptic input trains, but instead generated the stochastic synaptic current as a continuous stochastic process with Gaussian probability density and exponentially decaying correlation function (Risken, 1989). The parameters of the distribution and correlation functions were calculated according to the theoretical description (Eqs. 8 and 9) of the process for the net synaptic input. In these simulations, a spike event is registered when the membrane potential crosses a threshold level. This also evokes a hyperpolarizing current by activating a conductance that decays exponentially with time constant _h. The current's reversal potential E_h = 90 mV and the conductance amplitude is equal to the leak conductance. Runs of 20 or 4000 s were simulated to collect statistics for the membrane potential and firing, accordingly.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS CALCULATION METHODS RESULTS DISCUSSION REFERENCES

We studied the influence of the synaptic input's temporal structure on the statistical properties of membrane fluctuations and firing in the neuron model (see Theoretical Methods). Although the theory is not restricted to stationary inputs or to particular temporal structure, we studied the effects of stationary uncorrelated series of events (homogeneous Poisson point process, where g₁ = const, g₂(t₁, t₂) = 0) with exponentially decaying synaptic current. As mentioned in the Methods section, our results are also applicable for the case when synaptic currents decay very fast but temporal correlation of the input trains decays exponentially.

For an exponential unitary synaptic current s_I(t  t') = exp((t  t')/_s) Eq. 1a can be solved to obtain the unitary postsynaptic membrane potential: v(t  t') = _sm(exp[(t  t')/_s]  exp[(t  t')/_m)]/[C_N(_s  m)] (if _m  s). Since the mean injected current I₁ = _sg₁, the mean membrane potential V₁ = _sg₁R_N, where  is mean synaptic current amplitude (nA), and R_N is neuron's input resistance (M). By using Eq. 9, the auto-correlation function evaluates to:

(10)

where  = t₁  t₂, and the standard deviation (SD) of the membrane potential

(11)

In order to check the accuracy of our theory, simulations were performed for different rates of synaptic inputs. The synaptic amplitudes had exponential distribution and their mean value, , was adjusted in all simulations to cause membrane potential fluctuations with SD of 5 mV. The mean membrane potential was kept equal to the resting potential by injection of a steady current opposite to the derived mean synaptic current (Eq. 8). When the arrival rate for synaptic events was higher that 5 events/ms, the probability density for the membrane potential approached a Gaussian function, which is satisfactorily described by the theory (Fig. 1).

View larger version (18K): [in this window] [in a new window]   FIGURE 1   Dependence of membrane potential probability density on the rate of (stationary Poisson) input events in the case of synaptic input modeled as current injection. For low input rate, the density was skewed. As the rate was increased, the probability density approached a Gaussian function calculated analytically (see text). The average exponentially distributed amplitude of the synaptic current, , was calculated a priori in each case to ensure a prescribed SD for the membrane potential according to Eq. 11. A steady current equal and opposite in sign to this average was added as a counterbalancing input to keep the mean membrane potential near rest. The membrane time constant was 5 ms; the synaptic current decay time was constant at 2.5 ms;  was equal to 15.5 pA and 5 pA for the rate of 0.5 and 5 events/ms, respectively.

In real neurons the injected synaptic current is due to transient changes in synaptic conductance, which is accounted for in Eq. 2. Although analytic expressions for Eq. 2 could not be obtained, it is possible to estimate the membrane potential statistics. For a stationary process, the mean of the membrane potential, V₁, satisfies the equation V₁ = G_s1E_s/(G_s1 + 1/R_N), where G_s1 = _sg_1G is the mean of the net synaptic input conductance and _G is mean unitary synaptic conductance (µS). Thus, V₁ = E_s/[1 + 1/(_sg_1GR_N)]. In the case of membrane potential fluctuations with SD small relative to E_s, we can approximate the synaptic current as I(t)  _G(t)(V₁  E_s), and use the same equation (Eq. 1a) for the voltage dynamics. By rescaling the membrane time constant to '_m = _m/(1 + _sg_1GR_N) and scaling potentials by (V₁  E_s) in Eqs. 10 and 11, we obtain an approximate correlation function and SD of the membrane potential:

(12)

(13)

To check these results, we performed simulations for different mean values of the synaptic conductance amplitude. The rate of incoming synaptic events was set to 10 events/ms (see Discussion) and the synaptic current's reversal potential, E_s, was equal to 50 mV above the resting potential. The synaptic amplitudes had exponential distribution and their mean value was calculated from Eq. 13 to ensure a prescribed membrane potential standard deviation. As shown in Fig. 2 A, the probability density of the membrane potential is almost Gaussian and is well described for a broad range of standard deviations despite up to fivefold decrease of effective membrane time constant, '_m. The Gaussian shape is well preserved for the membrane fluctuations that had standard deviations <E_s/10. The correlation function also follows the predicted shape (Fig. 2 B).

View larger version (20K): [in this window] [in a new window]   FIGURE 2   Statistics of the membrane potential for the case of synaptic conductance inputs. (A) Theoretical and simulated probability density for the membrane potential. Analogous to Fig. 1, the synaptic conductance amplitude, G, was adjusted to achieve the appropriate value of V; a compensatory steady input was included. The input rate was fixed at 10 ms1. Note that for high SD the density in the tails deviates from the Gaussian form. (B) The autocorrelation function of the membrane potential closely matches the theoretically predicted curves. The membrane time constant without synaptic input was 5 ms; the synaptic current decay time constant was 2.5 ms; G had values of 0.05, 0.07, 0.1, and 0.5 nS for SD of 2.5, 3.3, 5, and 6.6 mV, respectively.

Having described the statistical properties of membrane potential fluctuations, we next estimate the frequency with which the membrane potential crosses some threshold level. Although the first-passage time problem is usually solved by using the Fokker-Planck equation (Tuckwell, 1988; Plesser and Tanaka, 1997), we applied another method that uses the correlation functions of membrane potential (Stratonovich, 1967). The probability to cross a given threshold level, V_th, with some rate of potential change,  dV/dt, is equal to w_j(, V_th)V, where w_j is the joint probability density. After noticing that V = t, we can write for the rate of threshold crossings from below

(14)

For stationary Gaussian-distributed membrane potential fluctuations, dV/dt and V are statistically independent, thus: n₀ = w(V_th)|dV/dt|/2, where w is the probability density of the membrane potential. Since the rate, being a linear transformation of membrane potential, also is Gaussian-distributed, the mean value for the absolute rate, |dV/dt|, can be calculated knowing only the standard deviation of the rate, which equals d²k_V2()/d²|₌₀:

(15)

After substitution of Eqs. 10 and 11 into 15, the firing rate in the model neuron can be estimated as

(16)

The latter estimate was checked by simulations (see Methods) for different decay times of the synaptic current and different threshold levels (Fig. 3 A). In order to directly compare with analytical results from Eq. 16, we changed the mean amplitude  of the unitary synaptic current according to Eq. 11 to have the same SD for membrane fluctuations in cases with different decay time, _s. The synaptic current decay time constant varied from 2.5 to 40 ms, and threshold level for the spike generation was set to 1_V, 2_V, and 3_V from the mean level of membrane fluctuations.

View larger version (18K): [in this window] [in a new window]   FIGURE 3   Statistics of firing for the integrate-and-fire-like model. (A) Theoretical values (open symbols) satisfactorily describe the firing rate for different synaptic current decay times if the spike threshold is >2 SD of the membrane potential fluctuations. For lower thresholds the rate is distorted by temporal correlation (see text). (B) ISI histograms are compared for fast, s = 2.5 ms, and slow, s = 40 ms, synaptic decay times when the threshold level Vth = V. Theoretical curves (solid lines) fit the distribution for longer ISIs. The number of short ISIs is reduced due to spike afterhyperpolarization. Note, for slowly decaying synaptic current the histogram's sharp peak for short-to-medium ISIs reflects bursts of spikes evoked by broad fluctuations of the synaptic current. The membrane time constant was 10 ms and spike afterhyperpolarization decayed with a time constant of 5 ms.

For very small fluctuations when _V < (V_th  V₁)/2, our estimate coincides with the firing rate of the model neuron (Fig. 3 A). For a given _V, the rate of firing decreased proportionally to 1/ as the decay time of the synaptic current was increased. However, for high firing rates there was a discrepancy between the predicted firing rate and the simulated one. For fast-decaying synaptic currents the simulated rate was less than estimated, and for slowly decaying synaptic currents the simulated rate was higher than predicted.

The reasons for these differences could be understood from the interspike interval (ISI) histogram (Fig. 3 B). For long interspike intervals the distribution is well predicted by the exponential density with 1/n₀ as the decay constant (Fig. 3 B), as one expects for independent events when temporal correlation does not play any role. However, the number of short interspike intervals is much less than predicted. For slow- and fast-decaying synaptic currents the histogram peaks at almost the same ISI value, indicating that spike afterhyperpolarization is responsible for this effect. Thus, for fast-decaying currents, afterhyperpolarization reduces the rate of firing (Figs. 3 B and 4 A).

View larger version (17K): [in this window] [in a new window]   FIGURE 4   Effects of spike afterhyperpolarization on the number of spikes and coefficient of variation. (A) Simulated traces illustrate a difference between the number of threshold crossings and spikes predicted, when spike afterhyperpolarization is or is not included. The dotted trace is from a simulation without evoked afterhyperpolarization. The continuous lineshows that spike afterhyperpolarization reduces the number of spikes (arrows) in this case of fast synaptic current decay (s = 2.5 ms). (B) When synaptic current decays slowly (s = 40 ms) very long peaks cause bursts of spikes (arrows) and increase the rate above that predicted theoretically. The same bursts are responsible for the excess number of short-to-medium ISIs (Fig. 3 B) and for the increase of the coefficient of variation (panel C).

For slowly decaying synaptic currents there is an excess of interspike intervals with medium duration (Fig. 3 B) and the rate of firing is increased relative to the prediction (Fig. 3 A) despite the reduced number of short ISIs. As the correlation decay is slowed (_s increased) the fluctuations become smoother and the duration of the input peaks increases. Some suprathreshold excursions are longer than the spike afterhyperpolarization, as can be seen in the trace of the membrane process (Fig. 4 B). Therefore, bursts of spikes can occur during a single excursion caused by the slower decaying synaptic current. Also, the ISI's coefficient of variation, CV, increased with synaptic current decay and could be higher than 1 (Fig. 4 C). The histogram (Fig. 3 B) has a long, slow tail but with a sharp peak at short-to-medium ISIs (because of the bursts), leading to large CV.

Assuming that burst-like spiking is responsible for the increased CVs, it is possible to estimate the ISI CV from the calculated, n₀, and simulated, n_s, firing probabilities. Since ISIs are relatively small for intraburst spiking, the ISI variance, which depends on the squared ISIs, reduces only slightly if we omit this contribution. The interburst intervals have approximately exponential distribution (Fig. 3 B). Furthermore, we took into account the fact that the mean-squared interburst interval is two times the mean interval squared, 2/n₀², for an exponential distribution. Since interburst intervals constitute a fraction, n₀/n_s, of all ISIs, the variance can be approximated as (2/n₀²) · (n₀/n_s)  1/n_s². Dividing ISI's variance by the squared mean interval, 1/n_s², we get:

(17)

If the difference between estimated and simulated firing probabilities is small, then CV  1 + (n_s  n₀)/n₀. Despite the underestimation due to the reduced variance, Eq. 17 predicts quite well the observed coefficient of variation (Fig. 4 C). For fast-decaying synaptic currents n_s is less than n₀ because spike afterhyperpolarization prevents some spikes from occurring (Fig. 4 A). Again, this simple equation provides a good estimate for the coefficient of variation (Fig. 4 C) showing that the effect of spike afterhyperpolarization can account for the changes in ISI CV.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS CALCULATION METHODS RESULTS DISCUSSION REFERENCES

In this study we used the theory of random point processes to describe the statistics of membrane potential and neuronal firing rate. For a neuron that receives random transient synaptic conductance inputs the theory allows one to calculate the distribution of membrane potential and the correlation function. The firing rate can be estimated analytically for low firing rates when spike afterhyperpolarizations do not change the rate significantly. For lower spike thresholds and higher firing rates spike afterhyperpolarization reduces the rate for fast-decaying synaptic inputs. However, the rate is increased if synaptic input correlations decay slower than spike afterhyperpolarizations. In this case spike bursts may occur during broadly peaked (single or composite) inputs, leading to a high coefficient of variation for ISIs, possibly exceeding one.

The derived relation (Eq. 9) between the correlation functions of membrane potential and synaptic input indicates the importance of spike synchrony in the efficiency of synaptic input to evoke a spike. For example, if synaptic input is completely decorrelated (Poisson train), only the rate of the incoming spikes defines the amplitude of membrane fluctuations. However, if a correlation between presynaptic spikes is induced in the train with the same rate, then fluctuations of membrane potential become stronger due to the additive term in Eq. 9. It could be that this effect of spike correlation is used to gate signals in the nervous system, for example, as suggested in a recent study (Steinmetz et al., 2000), which found that attention increases the spike synchrony in neurons of somatosensory cortex in behaving monkeys.

Although we used only exponentially decaying synaptic currents and conductance, our theory could be used for any dynamically structured synaptic conductance. Thus, the alpha function that is widely used to simulate synaptic conductance (Koch and Segev, 1998) can also be incorporated. We did not explore this particular case since our goal was to capture the essential influence of the correlation decay time on the membrane potential and firing statistics.

As shown in Fig. 1 the membrane potential distribution approaches a Gaussian for high input rates. For these simulations the input rate was 5 events per ms and the membrane time constant was equal to 5 ms. Thus, it can be concluded that for the distribution to be similar to a Gaussian the rate should be an order of magnitude higher than the reciprocal of the system's slowest time constant. Could these conditions be fulfilled in vivo? For cortical neurons, which have ~10⁴ synaptic contacts, this rate could be achieved because of the spontaneous activity of neurons (a few Hz) and the constant spontaneous generation of synaptic events. More direct evidence for the theory's applicability is seen in intracellular recordings in vivo, where the membrane potential fluctuations are reported as Gaussian-distributed (Calvin and Stevens, 1968; Stern et al., 1997; Pare et al., 1998; Azouz and Gray, 1999). Even if the input rate is not very high, the theory can describe the fluctuations, but then higher-order correlation functions from Eqs. 3 and 7 are needed for a more accurate description. However, for high-rate net synaptic inputs with large unitary conductance the theory satisfactorily describes membrane fluctuations (Fig. 2 A) even when average synaptic conductance causes severalfold decrease in the effective membrane time constant, which was observed experimentally in vivo (Pare et al., 1998; Destexhe and Pare, 1999).

In all simulations we used a single excitatory stochastic input process in order to study membrane potential fluctuations, and we compensated by subtracting a steady current equal to the average net random input current. We did so in order to avoid introducing additional parameters that could make our understanding of the phenomena less transparent. In principle, Gaussian distributed amplitudes of synaptic conductance could be used to represent excitation and inhibition. For synaptic inputs with different current decay rates, for example mediated by NMDA and AMPA receptors, separate processes could be used since the equations describing the system are linear. However, inhibitory inputs are more difficult to describe due to the small difference between resting membrane potential and synaptic reversal potential. In this case our theory could be applied only for small mean synaptic amplitudes causing small membrane fluctuations near the spike threshold.

Equation 16 gives the firing rate only in the case when the threshold for spike generation is higher than the average membrane potential. This assumption could be applicable at least for cortical neurons. High variability of firing in cortical neurons and intracellular recordings in vivo suggest that cortical neurons could operate under a balance of excitation and inhibition (Ferster, 1986; Shadlen and Newsome, 1998). In this case of constant integration of excitatory and inhibitory inputs neurons could remain near, but just below, the threshold, and firing could be due to random crossings of the threshold for spike generation (Shadlen and Newsome, 1998). Indirect evidence suggests that this could also be true for some other neural systems, too. In the spinal cord, stimulation of the pyramidal tract evokes both excitation and inhibition of almost equal strength in some types of motoneurons (Binder et al., 1998).

Although Eq. 16 provides a good estimate only for low firing rates, it could also be used to analyze cases with high firing frequency. In this case Eq. 16 provides only an estimate for the rate of long interspike intervals (Fig. 3 B), but not for intraburst firing rate. High firing rates, which are observed in experiments, could be accounted for by slow temporal correlation among spikes in presynaptic neurons (Brivanlou et al., 1998; Weliky and Katz, 1999; Gray et al., 1989; Roelfsema et al., 1997) and/or slow synaptic currents (Fig. 4 B). Usually, the NMDA receptor-mediated current decays over several tens or hundreds of milliseconds (Silver et al., 1992; Lester et al., 1990) and, according to Eq. 9, could be partially responsible for slow membrane fluctuations observed experimentally (Azouz and Gray, 1999; Lampl et al., 1999). As can be seen from Figs. 3 B and 4 B, firing frequency inside a burst is 40 Hz. If afterhyperpolarization decays faster, it is possible to observe even higher intraburst firing rates, which causes high variability of ISIs. Previously, high ISI coefficient of variation was explained by the interaction of synaptic input with potential dependent currents in a postsynaptic neuron (Softky and Koch, 1993; Wilbur and Rinzel, 1983). Since ISI CVs much higher than 1 are often observed in in vivo recordings (Victor and Purpura, 1998), possibly very general mechanisms, like the proposed slow membrane fluctuations, are responsible for this variability.

The theory as used here to estimate firing rate is limited to neurons with only simple fast spike generation mechanisms. Many different voltage-dependent currents that may have slow time scales start to activate below the firing threshold. Our theory provides only qualitative understanding in such cases; it would have to be extended to be considered as quantitative. However, our description of the synaptic input as a stochastic process (Eqs. 8 and 9) from the theory of random point processes is quite general. The mean, k₁(t), and correlation function, k₂(t₁, t₂), of the membrane potential recorded in vivo (Lampl et al., 1999; Azouz and Gray, 1999) could be used to extract mean synaptic input rate, g₁(t), and correlation between presynaptic spikes, g₂(t₁, t₂), by using Eqs. 8 and 9. Also, the theory will help in generating temporally structured synaptic input in studies of complex nonlinear models of neurons.