Abstract
Schmich, Robert M. and Michael I. Miller. Stochastic threshold characterization of the intensity of active channel dynamical action potential generation. J. Neurophysiol. 78: 2616–2630, 1997. This paper develops a stochastic intensity description for action potential generation formulated in terms of stochastic processes, which are direct analogues of the physiological processes of the pre and postsynaptic complex of the cochlear nerve: 1) neurotransmitter release is modeled as an inhomogeneous Poisson counting process with release intensity μ_{t}, 2) the excitatory postsynaptic conductance (EPSC) process is modeled as a marked, linearly filtered Poisson process resulting from the linear superposition of standard shaped postsynaptic conductances of size G, and 3) action potential generation is modeled as resulting from the EPSC exceeding a random threshold determined by active channel dynamics of the HodgkinHuxley type. The random threshold is defined to be the least upper bound in the size of a standardshaped neurotransmitter release injected at time t given the previous action potential time and the number of releases occurring in a short preconditioning time increment. The action potential process is modeled as a selfexciting point process with stochastic intensity resulting from the probability that the random threshold process crosses the threshold in some small time increment that is a function of time since previous action potential, release intensity, and the probability that a single synaptic event exceeds the stochastic threshold. The stochastic intensity model is consistent with a direct simulation of the nonlinear HodgkinHuxley differential equations over a variety of parameters for the vesicle release intensity, vesicle size, vesicle duration, and temperatures. Results are presented showing that the regularity properties seen in the vestibular primary afferent in the lizard, Calotes versicolor, associated with a slowtoactivate potassium channel resulting in a long afterhyperpolarization can be accommodated directly by the stochastic intensity description. The stimulus dependence of the model is attributed to synaptic transmission and the probabilistic nature to the threshold conductance process, which is dependent upon the EPSC process. The stochastic intensity is seen to have a form consistent with the phenomenologically based SiebertGaumond model, a stimulusrelated function of time multiplied by a refractoryrelated function of time since previous action potential.
INTRODUCTION
This paper develops a stochastic intensity description for action potential generation in the cochlear nerve associated with active channel dynamics of the HodgkinHuxleytype. The model is formulated in terms of coupled stochastic processes that are direct analogues of the following physiological processes: 1) vesicle release is modeled as an inhomogeneous Poisson counting process {M_{t} , t ≥ 0} with random event times {u_{i} , i = 1, . . . , M_{T} } and release intensity μ_{t} = lim_{Δ→0} Pr{M _{t,t+Δ} = 1}/Δ, 2) excitatory postsynaptic conductance (EPSC) {G _{t} , t ≥ 0}, which is modeled as a marked, linearly filtered Poisson process consisting of a superposition of standard shaped neurotransmitter releases of size G, and 3) action potential generation is a counting process {N _{t}, t ≥ 0} with event times {w _{i}, i = 1, . . . , N _{T}}. These are depicted in the Fig. 1, left.
Our approach for modeling spike generation is to characterize the threshold properties of nonlinear active channel dynamics. We model action potential generation as resulting from the random EPSC process passing threshold for the nonlinear active channel dynamical systems of the HodgkinHuxley (HH) type. Action potential generation is modeled as a selfexciting point process with intensity determined by the neurotransmitter release intensity μ_{t}, stochastic threshold, and synaptic conductance conditioning history. Fundamental to our model is the notion that action potential generation results from the random EPSC process driving the nonlinear channel dynamics to cross threshold on random times {w _{i}, i = 1, . . . , N _{T}}, which form the action potential process. To calculate the intensity of discharge, we define the random synaptic threshold conductance function θ(t − w _{Nt } = τ, m) as the justnecessarysized neurotransmitter quanta injected at time t required to cause an action potential at time t as a function of the time since previous action potential t − w _{Nt } = τ and the number of preconditioning vesicles M _{t−H,t} = m occurring in a small history of the synaptic conductance process [t − H, t). By calculating the probability laws on the sequence of thresholds θ(τ, 0), θ(τ, 1), . . . , we can derive the stochastic intensity of discharge conditioned on the history of the action potential times. The parameters of the model are the temperature, vesicle duration, vesicle intensity and size, correlation, and channel characteristics of the postsynaptic membrane.
For characterizing a particular active channel system, we empirically calculate the stochastic threshold distribution Pr{θ(τ, m) ≤ G} as a function of time since previous action potential τ and preconditioning events m by adapting the channels and temperatures of a single compartmental HH model to the system of interest and empirically calculating the stochastic threshold for that system. For simulating the threshold characteristics of the afferent terminal, excitatory synaptic currents, passive membrane components, and the voltagegated currents generating the action potential are included. These are depicted in Fig. 2, left, which shows the idealized model with membrane capacitance in parallel with the sodium, potassium, leak, and synaptic conductances G ^{Na} _{t}, G ^{K} _{t}, G ^{L}, and G _{t}, having reversal potentials V ^{Na}, V ^{K}, V ^{L}, V ^{S}, and postsynaptic membrane voltage V _{t}. Throughout the paper, a variety of stochastic threshold functions have been generated from the nonlinear HH systems, each one reflecting the range of systems that are of current interest to us, including fish vestibular neurons at 17–27°C (Boyle and Highstein 1990; Rabbitt et al. 1994, 1995) as well as bird (Boyle et al. 1994) and, of course, mammalian auditory nerve at 37°C. As well, we present results that model the vestibular primary afferent in the lizard Calotes versicolor at 32°C as described by Highstein and coworkers (Schessel et al. 1991). For examining the regularity properties of neural discharge in the vestibular system, we include a slowtoactivate potassium channel causing long afterhyperpolarizations. HH computation results are presented that illustrate excellent agreement between the stochastic intensity of discharge model and the HH spike generation model.
A strong departure of this work from previous work in Miller and Wang (1993) is the predominant emphasis on the stochastic intensity associated with time since previous action potential strictly larger than dead time. Herein we do not explore the nonzero probability behavior associated with the discontinuity in threshold as active channel systems release from dead time. As well, we parameterize threshold in terms of synaptic conductance directly, not requiring the small signal approximations assumed previously. This requires a characterization of the active channel systems in terms of a stochastic threshold measured in synaptic conductance units. As demonstrated, the stochastic synaptic conductance threshold determines the stochastic intensity of action potential generation. Interestingly, we show that the model is consistent with the product model of Siebert and Gray (Gray 1967; Siebert and Gray, 1963) and Gaumond, Kim, and Molnar (Gaumond et al. 1982, 1983). The stochastic intensity derived here has a natural factorization in terms of synaptic vesicle release intensity μ_{t} and the threshold characteristics of the postsynaptic active channel dynamical system having a clear functional dependence on the time since previous action potential. However, as discussed below, the stochastic intensity can be a nonlinear function of vesicle release intensity and can diverge from the linear model at varying stimulus levels.
METHODS
Stochastic HodgkinHuxley threshold characteristics
Neurotransmitter release is modeled as a Poisson counting process {M_{t}
, t ≥ 0} with intensity μ_{t} = lim_{Δ→0} Pr{M
_{t,t+Δ}
= 1}/Δ. At this time, we take this as an ordinary inhomogeneous Poisson process; this has been extended as described in Wang et al. (1994a,b) and Schmich (Schmich 1996; Schmich and Miller 1997), including shortterm adaptation reflected by the synaptic state. For our purposes here for accounting for steady state responses such as spontaneous activity, μ_{t} = μ is a fixed deterministic, unknown parameter. The resulting EPSC process, {G_{t}, t ≥ 0}, is modeled as a marked, linearly filtered Poisson process (Snyder and Miller 1991, chapts. 4 and 5) resulting from each neurotransmitter release adding an independent, short time increment of synaptic conductance of some random height in the postsynaptic cleft independent of others. Linear superposition of quantal conductances is assumed. The single transmitter response is idealized to a fixed waveform of standard shape Π(t), t ≥ 0. A response shifted to random time u
_{i} with height G
_{i} is denoted as G
_{i}Π_{ui
} (t). The synaptic conductance process {G
_{t}, t ≥ 0}, due to a number of such conductance enhancements occurring on Poisson random times {u
_{1}, . . . , u
_{MT
}}, becomes G
_{t} =
For calculating the random threshold properties associated with particular active channel dynamics, we use the standard HHtype differential equation, supplemented with the standard kinetic equations for both G
^{Na}
_{t} and G
^{K}
_{t}, for describing the input and output relationship between postsynaptic conductance G
_{t} and membrane voltage V
_{t}
The action potential process {N_{t} , t ≥ 0} and its random times {w_{i} , i = 1, . . . , N_{T} } are generated by solving the stochastic Eq. 1 and defining the action potential times to be the time when the voltage waveform increases rapidly toward V _{Na} and obtains its peak value. These are depicted in the Fig. 1, bottom right.
To calculate the intensity of discharge of action potential generation as a function of time since previous spike t − w_{Nt} , we must be able to calculate to o(Δ), Pr{N_{t,t+Δ} = 1‖t −w _{Nt } = τ}, where the random measure on the spike generation process is induced by the measure on the synaptic conductance processes {G _{t}, t ≥ 0}. For this, the HH generation of action potentials is characterized as a nonlinear thresholder. The action potential times {w _{i}, i = 1, . . . , N _{T}} result from the random threshold crossing of the synaptic conductance process {G _{t}, t ≥ 0}.
Begin by defining threshold, a random process totally determined by the synaptic conductance sample paths. Define the random element ω_{t} ≡ {G _{σ}, −∞ < σ < t} a sample path of the EPSC process in the sample space 𝒢_{t} = {ω_{t} = {G _{σ}, −∞ < σ < t}} of conductance processes. Also define the subset of elements 𝒢_{t}(τ) ⊂ 𝒢_{t} of sample paths having no action potentials in [t − τ, t), 𝒢_{t}(τ) = {ω_{t} ∈ 𝒢_{t} : t − w _{Nt } = τ}.
The random threshold process {θ (t − w _{Nt } = τ, ω_{t}) t > −∞} is taken as a function of ω_{t} ∈ 𝒢_{t}(τ), the random sample conductances for which the time since previous action potential is t − w _{Nt } = τ. Definition 1: The random threshold function θ(τ, ω_{t}) at time t is defined as the justnecessary standardshaped vesicle of size ϑ injected at time t required to cause an action potential with preconditioning ω_{t} ∈ 𝒢_{t}(τ) where t − w _{Nt } = τ. See appendix a and Eq. EA3 for a formal definition.
Figure 2, right, shows an example of the threshold function computed throughout the conductance sample path ω _{t} = {G _{σ}, 0 ≤ σ < t} with an action potential assumed to reset the system at time t = 0. Figure 2, right, shows θ(t − w _{Nt }, ω_{t}) computed from the HH dynamical system according to Eq. EA3 generated by solving the HH differential equation at 17°C. Notice how the threshold resets dramatically after an action potential; also notice the absolute dead time T _{D} = 3.0 ms.
We have assumed implicitly that the threshold function θ(t − w _{Nt }, ω_{t}) is not an explicit function of time t only in its dependence on the time since previous action potential and on the conductance sample path. As well, we assume θ(t − w _{Nt }, ω_{t}) is continuous in its argument t, i.e., θ(t + ζ − w _{Nt+ζ}, ω_{t+ζ}) → θ(t − w _{Nt }, ω_{t}) as ζ ↓ 0. This assumes t − w _{Nt } > T _{D}, the absolute dead time! Not surprisingly, the threshold function is discontinuous at absolute dead time τ = T _{D}. As shown in Miller and Wang (1993), this must be handled separately when describing stochastic intensities. Our analysis throughout assumes t − w _{Nt } > T _{D} so that continuity holds.
Notice that in defining the threshold function, we have conditioned explicitly on the time since previous action potential precisely because channel kinetics change dramatically as a function of time since previous action potential. This allows for the efficient encoding of the conditional intensity of discharge. This is obvious as shown in Fig. 2, right. The selfexciting point process conditioned on conductance sample paths becomes
To examine the history parameter somewhat more quantitatively, we have computed the conditional covariance between the HH action potential process {N_{t}
, t ≥ 0} and the input synaptic conductance process {G_{t}
, t ≥ 0}
With the notion of equivalence defined through finite history, we compute the probability law on stochastic thresholds by dividing the sample paths into subsets 𝒢_{t}(τ), τ ≥ 0, of sample paths for which there are no action potentials in [t − τ, t) with t − w
_{Nt
} = τ, and subsets 𝒢_{t}(τ, m) of elements having M
_{t−H,t} = m vesicles released in [t − H,t)
Empirical generation of stochastic threshold
Our strategy is to characterize each active channel system by empirically computing the stochastic threshold from a nonlinear dynamical system with parameters chosen to reflect that of the system being characterized. The probabilistic and average properties of threshold can be calculated straightforwardly for subsets 𝒢_{t}(τ, m), m = 0, 1, 2, . . . . Shown in Fig. 5 are the mean thresholds, denoted as θ̄(τ, m), characterizing the average threshold dynamics for the m = 0, 1, 2 subsets. The mean thresholds are calculated by initially causing the system to spike at time −τ, driving the system forward to time 0 with M _{−H,0} = m synaptic conductances as the preconditioning stimulus and then calculating the conductance required at t = 0 to elicit an action potential by altering the height of a standardshaped tansmitter release Π_{0}(⋅). The preconditioning is varied over all Poisson release times for M _{−H,0} = m, m = 0, 1, 2. Figure 5, top, shows a pictorial representation for generating these characterizations. Notice, there is no randomness for them = 0 subset because there are no preconditioning stimuli in the preconditioning interval [−H, 0). Figure 5, bottom, shows the associated HH mean threshold functions θ̄(τ, m) for the m = 0, 1, 2 subsets as a function of the time since previous spike τ, with varying parameters for temperature (17, 27, and 37°C). The mean threshold curves θ̄(τ, m) illustrate the discharge history dependence in the model: during the absolute dead time interval, the mean threshold is effectively infinite, after which there is a decrease to the steady state threshold conductance with various time courses depending upon temperature and membrane currents. Also notice that both θ̄(τ, 1) at 27 and 37°C and θ̄(τ, 2) at 17, 27, and 37°C decrease to 0 because the preconditioning stimuli themselves are sufficient to generate action potentials.
The HH mean threshold function θ̄(τ, 0) shown in Fig. 5, bottom left, also explains the difference in the covariance curves seen in Fig. 4 for 17°C (left) and 27°C (right). Notice that the covariance curves for τ = 4, 12 ms at 17°C are not equivalent, whereas the covariance curves for τ = 2, 12 ms at 27°C are equivalent. This difference is explained by the shapes of the mean threshold curves shown in Fig. 5 with the threshold curve for 27°C reaching steady state level near τ = 2 ms; the threshold curve for 17°C does not reach steady state until τ = 13 ms. Thus the mean thresholds for τ = 2, 12 ms at 27°C in Fig. 4 are equal, causing the respective covariance curves to be equivalent. For 17°C, the mean thresholds at τ = 4, 12 ms are not equivalent, implying the respective covariance curves will be different.
Stochastic intensity of discharge
To study the global probabilistic behavior of the action potential process, it should be described via its intensity, which is the local probability of discharging in a small increment conditioned on the history of the process (chapter 6, Snyder and Miller 1991) and is given by
To calculate conditional probabilities over waveforms containing m vesicles with no action potentials N_{t−τ,t}
= 0 for time duration τ = t − w
_{Nt
}, define γ(τ, m) ∈ [0, 1] to be the probability of the subset 𝒢_{t}(τ, m) conditioned on 𝒢_{t}(m)
Equation 6 is a stochastic intensity description of action potential generation associated with active channel dynamics of HH type models. Essentially, this is the superposition of all events associated with m = 0, 1, 2, . . . preconditioning transmitter quanta, multiplied by the probability of one synaptic event in [t, t + Δ] causing an action potential, which goes as μ_{t} Pr{θ(τ, m) ≤ G}.
Figure 6 shows the calculation of the probabilities Pr{θ(τ, m) ≤ G} as a function of τ for m = 1, 2, 3, . . . . Note that because there is no randomness in the preconditioning interval for θ(τ, 0), Pr{θ(τ, 0) ≤ G} 1 if θ̄(τ, 0) ≤ G, and 0 otherwise. For the mth subset, the probability law is calculated for the justnecessary conductance required to elicit an action potential at time t. The exact location of the m preconditioned vesicles in [t − H,t) is varied according to the Poisson occurrence times, with the probability curves Pr{θ(τ, m) ≤ G} computed for m = 1, 2, 3, . . . as a function of time since previous spike τ by varying the placement of the preconditioning transmitter quanta in [t − H,t). These curves were generated by uniformly releasing m transmitter quanta in the interval [t − H,t) over n trials and empirically calculating the associated probability function Pr{θ(τ, m) ≤ G} for each τ,
RESULTS
Stochastic intensity model of HH dynamics
We use hazard rate histograms to compare the accuracy of the stochastic intensity model, λ_{t}(τ) of Theorem 1, to the output of the actual HH dynamical system driven by the EPSC process. Hazard rate functions (chapter 6, Snyder and Miller 1991) characterize the discharge intensity for such selfexciting point processes and have the following functional form: h[i] = (n[i])/(N _{T} − n[i]), where n[i] is the number of interarrivals in the ith bin and N _{T} is the total number of action potentials. This form for the hazard rate function assumes that the discharge intensity is stationary, i.e., λ_{t}(τ) = λ(τ), or equivalently that the stimulus function is stationary, μ_{t} = μ. The general nonstationary form for the hazard rate function for λ_{t}(τ) is derived in Mark and Miller (1992). The hazard function is essentially the maximum likelihood estimate of the intensity λ(τ) with stimulus μ_{t} = μ (Mark and Miller 1992). The simulation results are generated assuming a stationary form for the stimulus function and, subsequently, the discharge intensity.
Begin by comparing the stochastic intensity description and the HH active channel model using parameters from the classic paper of Frankenhaueser and Huxley (1964) at 20°C. Figure 7, left, shows the ionic currents resulting from the current stimulation at 20°C corresponding to Frankenhaueser and Huxley (1964), their Fig. 6. Figure 7 (right, dotted line) shows the hazard discharge rate generated from solving the differential Eq. 1 for action potential generation and computing h[i] = (n[i])/(N_{T} − n[i]). Superimposed via the solid line is the stochastic intensity λ(τ) as a function of τ, the time since previous action potential. Notice the accurate fit of the stochastic intensity model of Theorem 1 with the HH generated data. The dead time parameter T_{D} is replicated exactly; the rapid rise in hazard intensity as well as the neural discharge passing into its relative refractory period are all duplicated as well.
Figure 8 shows a comparison of the stochastic intensity model λ(τ) ( ) and the HH dynamical equation's hazard rates (⋅⋅⋅) at 17°C. The vesicle size was increased fromG = 6.8 mS/cm^{2} (Fig. 8, top left) to G = 9 mS/cm^{2} (top right), and the release intensity was increased from μ = 200 s (bottom left) to μ = 300 s (bottom right). The single transmitter width was W = 0.5 ms (top) and W = 1.0 ms (bottom). Conditioning history was selected to be the vesicle width H = W in all four plots. Notice that the firing rate increases with transmitter quantal size G (top).
Shown in Fig. 9 are the discharge hazard rates generated from solving the HH model at 27°C (⋅⋅⋅) as a function of time τ since previous action potential. Shown superimposed ( ) is the stochastic intensity λ(τ). The stochastic models are shown with transmitter quantal size G increased from 6 to 8 in Fig. 9, top, and transmitter release intensity μ decreased from 500 to 250 s from the left to right panels in Fig. 9, bottom. The width, W = 0.5 ms in Fig. 9, top, was decreased to 0.3 ms in Fig. 9, bottom. The preconditioning interval H was selected to be the transmitter width H = W in all four plots.
Figure 10 examines the rapid dynamics of hightemperature systems such as is appropriate for mammalian auditory nerve. The stochastic intensity, λ(τ) ( ) and hazard rate intensities (⋅⋅⋅) are shown at 37°C with μ = 250 s (top left), μ = 375 s (top right), μ = 500 s (bottom left), and μ = 625 s (bottom right), with H = W = 0.5 ms. Notice, that the discharge rates for 37°C almost immediately reach their steady state level for τ > T_{D} ≈ 0.6 ms. This is consistent with the mean threshold curve dynamics θ̄(τ, m) shown in Fig. 5 for the 37°C curves (solid line).
Figures 810 show that the HH hazard rates and the stochastic intensity model decrease linearly with the release intensity and quantal size. As well, notice the agreement between the HH hazard rates and stochastic intensity model λ(τ). Figure 11 examines the model as the preconditioning increment H is varied. Figure 11, left, shows the stochasticintensity for width W = 0.5 and H = 0.5 ms (bold line),H = 0.7 ms (thin line), and the stochastic intensity withH = 0.3 ms (dashed line) all at 17°C. The HH hazard rates are shown via the dotted lines. Figure 11, right, shows results for transmitter width W = 0.3 ms and for the stochastic intensity for H = 0.3 ms (bold line), the stochastic intensity with H = 0.5 ms (thin line), and the stochastic intensity with H = 0.1 ms (dashed line) at 27°C. These panels illustrate a divergence of the stochastic intensity model from the solution to the HH differential equation when H < W, with the intensity model significantly below the HH hazard rates. The postsynaptic conductance history is not being sufficiently included in the calculation. Notice that the model's performance is slightly better for H > W (thin line) compared with when H = W (bold line), especially for the 17°C case (left). Figure 11, bottom, shows the covariance curves between the EPSC process G_{t} and the action potential process N_{t} with the same parameters for 17 and 27°C as in Fig. 11, top. The covariance curves demonstrate that the stochastic intensity model fits the hazard rate better for the 17°C case when H is increased because there is a covariance between G_{t} and N_{t} for s > W. The covariance plots for 27°C illustrate that G_{t} and N_{t} are uncorrelated for s > W and thus, increasing the history H has little effect on the stochastic intensity model.
It is interesting to point out that the stochastic intensity model for auditory nerve discharge predicts an intensity form that is similar to the SiebertGaumond product model type (Gaumond et al. 1982; Siebert and Gray 1963). The stochastic intensity model in Eq. 6
can be rewritten according to
In our simulations, we have found a mismatch between the stochastic intensity model λ(τ) and the solution of the HH dynamical system for extremely high discharge rates. Figure 13 demonstrates that the discharge rates obtained from the solution of Eq. 1 (⋅⋅⋅) diverge from the stochastic intensity λ(τ) ( ) at the highest discharge rates. In Fig. 13, μ is increased steadily from 500 s (left) to 2,000 s (right), with the result that the stochastic intensity model and the solution to the differential equation diverging at 500–1,000 spikes/s steady state.
Application to long afterhyperpolarization data
Figure 14 shows results that model the vestibular primary afferent in the lizard, C. versicolor, as described by Highstein and coworkers (Schessel et al. 1991). These neurons are examples of fibers that contain channels that produce large afterhyperpolarizations (AHP) due to a slowly activating potassium channel. The AHP results in regular discharge properties with preferred intervals of fixed duration, not exhibited in systems we have examined previously. To accommodate the regular discharge properties of these neurons, the HH differential Eq. 1 was modified to include a channel of functional form taken from Goldberg and coworkers (Smith and Goldberg 1986) being I_{k0} = g_{k0}e ^{− τ/τk}(V _{t} − V _{K}), where I _{k0} is the channel current, τ is the time since previous spike, g _{k0} and τ_{k} are adjustable parameters. The C. versicolor lizard data of Schessel et al. (Schessel et al. 1991) are exquisite because they provide excellent values for the parameters W, G, τ_{k}, g _{k0}. Results were generated at a temperature of 32°C, with width W = 0.5 ms, size G = 0.46 mS/cm^{2} (which corresponds to an excitatory postsynaptic potential (EPSP) with a size of ∼1 mV), vesicle release rate μ = 25,000 s, g _{k0} = 0.03, and τ_{k} = 8 ms. Figure 14, top left, shows a sample waveform for the postsynaptic membrane voltage V _{t} for Eq. 1, with the parameters in the FrankenhaueserHuxley model chosen to duplicate the voltage waveforms depicted in Fig. 17 of Schessel et al. (1991). Figure 14, top right, shows the interval histogram demonstrating the regularity in firing induced by the AHP channel. Notice the preferred interval of duration ∼20 ms. Figure 14, bottom left, shows the covariance curves between the synaptic conductance process G _{t} and the action potential process N _{t} plotted as a function of correlation distance s for τ = 20 ms ( ) andτ = 25 ms (⋅⋅⋅). Figure 14, bottom right, shows the discharge rates with the added channel in the HH differential equation (⋅⋅⋅) and the stochastic intensity model λ(τ)(– – –, ) as a function of time since previous spike τ. The stochastic intensity model was generated for two values of preconditioning interval: H = 0.5 ms (– – –) and H = 1.5 ms, with width W = 0.5 ms ( ). Notice that for the H = W = 0.5 ms case, there is a significant mismatch between the stochastic intensity model (– – –) and the hazard rate (⋅⋅⋅); this can be explained by inspecting the covariance curves in Fig. 14, bottom left. The covariance curves illustrate that there is a correlation between the conductance process G _{t} and the action potential process N _{t} for s > w, no doubt due to the additional slow potassium channel.
DISCUSSION
The stochastic intensity model was developed by direct modeling of various channel dynamics, as well as additional channels for various physiological systems such as the slowtoactivate potassium channel producing large afterhyperpolarizations in the vestibular system of lizards. Simulation results presented show that the stochastic intensity model λ_{t}(τ) fits the discharge rates obtained from solving the HH differential equation for a wide variety of values for the transmitter quantal conductance, release intensity, and channel dynamics. For simulations, parameters of the EPSC and HH system were chosen to parallel the various components of several wellstudied systems. Siegel and Dallos have recorded EPSPs from the guinea pig organ of Corti, where most afferent fibers make a single synaptic contact with one inner hair cell (Siegel 1992; Siegel and Dallos 1986). The single event EPSPs that did not cause a spike had a duration of ≤1 ms, with both symmetric rising and falling phases. This should be contrasted to the unsymmetric onsets and decays seen in other systems (Siegel 1992; Siegel and Dallos 1986). The simulations idealize these symmetric EPSPs and study the HH system for transmitter quantal durations of 0.5–1.0 ms, consistent with these time scales. Synaptic transmitter quantal widths were chosen to be comparable with the measurements in the guinea pig organ of Corti (Siegel 1992; Siegel and Dallos 1996) and recordings by Highstein and coworkers (Schessel et al. 1991) measuring miniature EPSPs (mEPSPs) in the primary afferents of the horizontal semicircular canal in the lizard at 32°C.
By comparison, in Yamada et al. (1989), fast synaptic inputs for the bullfrog sympathetic ganglion cells are modeled with postsynaptic waveforms of the form const ×te ^{−t/tpeak } with t _{peak} = 2.5 ms to fit experimental data measured by Kuba and Nishi (1979). To examine the dependence of threshold conductance results on waveform shape, we have examined various alternative forms for the synaptic conductance including const × te ^{−t/tpeak } as well as a triangular shape. By proper normalization of the calculation of the conductance threshold in units of total area for each conductance waveform, the threshold results were found to be not significantly affected by the exact shape of waveform used for the conductance. We have chosen to use the idealized rectangular waveform for the conductance in our modeling.
The amplitudes of the idealized synaptic conductance waveform also were chosen based on such experimental recordings and by the synaptic conductance threshold curves. As discussed in Siegel and Dallos (1986), primary auditory neurons are exquisitely sensitive, requiring on the order of one to several synaptic releases to cause a postsynaptic action potential. Several of the simulations studied the active channel dynamics at 37°C for synaptic conductance values on this order. As well, for systems with the slowtoactivate potassium channel (see Fig. 14), parameters were selected to be consistent with experimentally measured primary afferents in the lizard, which had mEPSPs with a mean amplitude of ≈0.9 mV (Schessel et al. 1991). In such systems, upwards of 10 mEPSPs are required to pass steady state threshold.
Throughout this paper, we have assumed that the time since previous action potential t − w _{Nt } = τ is always greater than the absolute refractory period or dead time of the system T _{D}, i.e., τ > T _{D}. In previous work (Miller and Wang 1993), it was shown that the selfexciting pointprocess model with bounded intensity does not hold for the time since previous action potential τ = T _{D}, but it instead has a nonzero probability of discharge that is not modeled here. The probability of discharge at τ = T _{D} predicts the nonmonotonic behavior of discharge rates as evidenced by an initial peak in the hazard rate histogram at τ = T _{D} seen in high spontaneous neurons and high driven rate neurons (Gaumond 1980; Miller and Wang 1993). As seen in Fig. 10, bottom, such nonmonotonic peaks of discharge rate are found in direct simulation of the HH dynamical system at high temperatures such as 27 and 37°C and in particular for large transmitter release rates μ or large quantal sizes G. We in fact have shown this in HH simulation (Schmich 1996). As this paper does not focus on the absolute dead time discharge intensity, we have chosen not to focus on these issues.
Another important aspect of the stochastic intensity model is the preconditioning increment parameter. The preconditioning interval H of the system is on the order of membrane integration for passive dynamics and expresses the fact that synaptic conductance history before time t − H has little effect on the system at time t conditioned on there being no action potential in [t − H, t). This allows a characterization via equivalence classes according to the number of synaptic events m released in the time increment [t − H, t). The preconditioning increment is determined by the correlation distance between the action potential process and the synaptic conductance. The increment was shown to be related inversely to the system's temperature (see Figs. 4 and 11) and, as shown in Fig. 14, a function of the channel dynamics. Figure 11 also showed the important role correlation history has in the prediction of the stochastic intensity.
The stochastic intensity model λ_{t}(τ) also predicts a product form for the intensity that is similar in form to the SiebertGaumond model: μ_{t} r _{t}(τ; μ). The refractory function in the SiebertGaumond model is not a function of the stimulus. However, as first argued in Miller and Wang (1993), the stochastic intensity is a nonlinear function of the vesicle release rate μ_{t} because of the nonlinear dependence of r _{t}(τ; μ) on vesicle release intensity. This is illustrated in Fig. 15, which shows the discharge rates obtained from solving the HH differential equation (⋅⋅⋅). Figure 15, left, shows the monotonic recovered probability MLE for the SiebertGaumond model (– – –); Fig. 15, right, shows the stochastic intensity λ(τ) model ( ) with the HH hazard rates superimposed in both panels. The three models are shown as functions of τ and release intensity μ at 37°C (top) and 27°C (bottom). Notice that in all plots, doubling release intensity more than doubles discharge rates, demonstrating that the intensities are nonlinear functions of the release rate. The curves were plotted on a log scale in the spirit of Gaumond et al. (1982), who obtained the product form for the intensity of discharge as s _{t} r(τ) by essentially plotting the measured discharge rates on a log scale for the y axis and noticing that the discharge rates differed by a multiplicative constant—the stimulus term s _{t}.
As depicted in Fig. 13, it may be significant that at the highest discharge rates the simulations diverge from the stochastic intensity model. A crucial assumption in the model is the conditioning on time since previous action potential. As first demonstrated by Gaumond et al. (1982) there is a divergence of the conditional independence for the shortest interval increments when conditioning on only the most previous action potential, suggesting dependence on multiple previous interspike times τ_{1}, τ_{2}. This would suggest a conditional intensity of the form λ_{t}(τ_{1}, τ_{2}), which is a function of the threshold θ(τ_{1}, τ_{2}, m), and its probability distribution P[θ(τ_{1}, τ_{2}, m) ≤ G]. It is possible that this is related to the divergence seen between the HH simulations and the stochastic intensity model.
An interesting aspect of the discharge rates at high temperatures, such as 37°C, is the prediction by the stochastic intensity model of the rapid onset to the steady state level (see Fig. 10). This rapid onset is consistent with results predicted by others (Lane et al. 1995). There also exists data (Bosch 1990; Gaumond 1980; Li and Young 1993) that show a slow component of recovery at 37°C that is not exhibited by a purely postsynaptic phenomenon in our model. This has caused various investigators (Li and Young 1993; Miller and Wang 1993) to conjecture that a possible cause for the slow component is a slowtoactivate potassium channel seen by SantosSacchi and others in spiral ganglion cells (SantosSacchi 1993; Siegel 1992). Our recent work demonstrates that by adding a slowtoactivate potassium channel to the HH model, the slow components of recovery can be duplicated for the specialized case of small conductance sizes such that four or more vesicles need to overlap to reach threshold; this is somewhat inconsistent with findings by Siegel (Siegel 1992; Siegel and Dallos 1986). Instead, we have found that the slow component of recovery as well as the fast component of recovery as referred to by Li and Young (1993) can be generated by altering time constants and parameters in the presynaptic models developed by Miller et al. (Schmich and Miller 1997; Wang et al. 1994a,b) for synaptic vesicle regeneration and transfer to the docking sites (Schmich 1996). As demonstrated in Schmich (1996), the slow and fast components of recovery may be associated with such presynaptic models.
Acknowledgments
We acknowledge F. Yang for help in the covariance simulations as well as E. Young, who originally encouraged us to generate the model in terms of synaptic conductance directly so that small signal approximations such as in Miller and Wang (1993) are not necessary. We also acknowledge S. Highstein for encouraging us to model the lizard system, Calotes versicolor, where physiological parameters have been characterized. We thank thereviewers for excellent suggestions for improving the paper as well as the suggestions for the possible explanation of the divergence of the stochastic intensity model from the HH solution at high spike rates. Finally, we thank J. A. O'Sullivan for helping to simplify the notational definition of the random threshold.
This research was supported by National Institute of Deafness and Other Communications Disorders Grants R01DC00333011 and 1P01DC0183701A2.
Footnotes

Address for reprint requests: M. I. Miller, Dept. of Electrical Engineering, Washington University, Campus Box 1127, St. Louis, MO 63130.
 Copyright © 1997 the American Physiological Society
APPENDIX A: A FORMAL DEFINITION OF THRESHOLD
Formally, construct a sample conductance ω′ ≡ {G_{σ}′, −∞ < σ < ∞} from ω_{t} ≡ {G
_{σ}, −∞ < σ < t} by injecting a standardsized vesicle of height ϑ at time t + Δ
APPENDIX B: PROOF OF THEOREM 1
Define P(⋅), the probability measure on the sample space 𝒢_{t} = {ω_{t} = {G
_{σ}, −∞ ≤ σ < t}} of conductance processes,
Proof
Noting
The last equation follows because θ(τ,0) is assumed to be continuous for all τ > T_{D} and θ(τ,0) > G, resulting in the first term Pr{N_{t,t+Δ} = 1‖M_{t,t+Δ} = 0, ω_{t}, 𝒢_{t}(τ)} = 0 for all τ > T _{D}.
Then removing the conditioning on ω_{t} gives the intensity
This gives