In this work we consider simulated macroscopic synaptic currents generated by a varying number of synaptic receptor channels. The channels are assumed to be independent and identical. We assume that the synaptic channel gating is a Markov process and p_ij is the probability of channel transition from state j to state i at time Δt. The rate matrix is an N_s × N_s matrix Q: qij=limΔt→0pijΔt, N_s is the number of states of the synaptic channel model. Each element of the matrix Q gives the rate constant of transition j → i if the transition is allowed by the model and otherwise q_ij = 0. The diagonal elements, q_ii, are set to − ∑jqij, so the sum over each column is zero. Synaptic release of neurotransmitter is modeled as a step pulse of its concentration in the synaptic cleft, which leads to instantaneous change of concentration-dependent transition probabilities, p_ij.

We assume that the kinetic matrix topology (i.e., a set of allowed transitions) and the number of conducting states are known. The required model parameters were arranged into the parameter vector θ = [q, i_ch, N_ch] and they were: rate constants, qij=limΔt→0pijΔt,i≠j, unitary currents, i_ch, and the number of postsynaptic receptors bound with a neurotransmitter right after the concentration transient, N_ch.

The likelihood function, L_θ, that is to be maximized by ML NSFA in order to find the most likely set of parameters is defined as the conditional probability to observe N macroscopic current traces c_i, i = 1 : N, sampled at time points t = [1… N_T] given the model parameters θ (Colquhoun and Hawkes, 1977; Celentano and Hawkes, 2004; Milescu et al., 2005; Stepanyuk et al., 2011):

Here N is the number of synaptic macroscopic current traces c_i (sample size) and N_T is the number of points in each trace; N_chi is a number of channels exposed to neurotransmitter in each current c_i; μ, an N_T × 1 vector and c_m1, an N_T × N_T matrix with elements {c_m1}_t,t′, and denote mean and covariance of single channel current, respectively, and they both are the functions of θ. c_m1 is related to the covariance matrix, c_m, of a macroscopic synaptic current c_i by the following expression: c_m = c_m1 N_chi. Mean and covariance follow equations (Colquhoun and Hawkes, 1977)

Here Q is a rate matrix (Colquhoun and Hawkes, 1977; Celentano and Hawkes, 2004) and p (0) is an initial state vector. The elements of p (0) can be calculated as the equilibrium probabilities determined by the initial experimental conditions, which are assumed to last for sufficiently long time T to allow the channels reach equilibrium

It is generally accepted to maximize the logarithm of the likelihood function logL_θ instead of the likelihood function L_θ itself. Therefore, our objective was to find the most likely model parameter set, θ_ML, i.e., the parameter set that maximized the log-likelihood

The log-likelihood function logL_θ can be efficiently estimated using the fact that c_m1 has a specific structure of semiseparable matrix (DeWilde and van der Veen, 1998; Stepanyuk et al., 2011).

Efficient log-likelihood estimation used in this article is based on our previously described method (Stepanyuk et al., 2011). Briefly, the method was based on the fact that the covariance matrix c_m1, has a specific structure of semiseparable matrix, namely c_m1 can be represented as Stepanyuk et al. (2011).

where

and

where U: e^Qt = Ue^DtU⁻¹ is N_s × N_s matrix of the eigenvectors of Q, and D is a diagonal form of Q provided all the eigenvalues of Q, λ, are different; N_o is a number of open states in the channel model. Efficient linear algebra algorithms for semiseparable matrices (Vandebril et al., 2007; Eidelman and Gohberg, 2008) allowed us to compute the log-likelihood and provided almost linear scaling of its computational cost with the number of states in a kinetic model for the case of sufficiently large number of currents, ensuring fast, and accurate estimation of model parameters. The number of synaptic channels exposed to neurotransmitter was assumed to be the same for all currents. However, in the case of synaptic currents this number could vary between trials due to quantal variability. As a result, logL_θ must be estimated separately for each current, and then summed up, thus increasing the number of operations in N times at least. However, calculations could be substantially simplified if the majority of receptors, which will participate in the current, are found in one particular state immediately after the neurotransmitter concentration transient, as it is expected for the synaptic receptors. To compute logL_θ in this case, let denote noisy macroscopic synaptic current with an N_T × 1 vector c_i and let denote by n_i an N_T × 1 vector of noise imposed on the i-th current. Then −logL_θ of the set of parameters θ given macroscopic synaptic currents without noise imposed on them is (we omit here the constant term NN_T log (2π)/2)

where μ is an expectation of synaptic current without noise and logL_θ (c − n) denotes the required log-likelihood given the set of macroscopic synaptic currents without noise. Equation (8) can be rewritten as

Since the background noise and the macroscopic current are uncorrelated the last term in Equation (9) can be neglected without loss in accuracy given the number of currents, N, is large enough. Therefore, we have

To quickly evaluate the last term in Equation (10), let us approximate ∑i=1N(niT)k(ni)jNchiby{cnoise}kjN〈1Nch〉. Hence,

Finally, from Equation (11) we obtain

where ° denote Hadamard multiplication.

Keeping in mind that ∑∑cm1−1°cnoise=tr(cm1−1cnoise), we rewrite Equation (10) for logL_θ (c − n) as

where

is the log-likelihood function of macroscopic synaptic currents with noise.

To quickly evaluate tr(cm1−1cnoise) we note that matrices c⁻¹_m1 and c_noise are quasiseparable (as an inverse of semiseparable matrix, Vandebril et al., 2007) and semiseparable matrice, respectively. Semiseparability of noise covariance matrix, c_noise, follows from the fact that experimental background noise can be well approximated by a stationary Gaussian process, and the covariance matrix of such process is semiseparable (DeWilde and van der Veen, 1998). Then, the computation of trace of the product of such matrices can be accelerated by representing it as tr(F · C)=2tr(H · B)+∑k=1NTFkkdk, where H is (N_T − 1) × N_S matrix, F is symmetric N_T × N_T semiseparable or quasiseparable matrix and B is defined by Equation (7) (see also Equations A1.26–A1.35 from Text S1 in Appendix in Stepanyuk et al., 2011).

To estimate the number of channels N_chi (see Results for further definition), we re-write Equation (8) for a single macroscopic synaptic current:

In the last expression we neglect the 3-d term, as it was done in Equation (9), and the 5-th and the last terms does not depend on N_chi at all. Leaving terms that depend on N_chi only we obtain log-likelihood as a function of the number of channels:

The number of channels, N_chi, can be approximated for each macroscopic synaptic current c_i as a number that gives maximum of the likelihood function when being substituted into Equation (16). Therefore, after differentiation of Equation (16) with respect to N_ch

we find an approximation for the number of channels, N_chi, for each macroscopic synaptic current

Here c_i is not the whole decaying phase of each current but only those part where signal-to-noise ratio is high and therefore noise term can be neglected (usually from peak of the current to 0.1–0.5 of the peak). Therefore, before calculating the log-likelihood, we first estimate N_chi for each macroscopic synaptic current, c_i, then substitute N_chi into Equations (13) and (14) and calculate the log-likelihood of the set of parameters θ given the set of simulated macroscopic synaptic currents. Accordingly, N_ch is estimated automatically when the maximization is finished.

The peak open probability, P(o, peak), was defined as a probability that a channel is opened at the peak of the macroscopic current given that this channel was bound with a neurotransmitter immediately after the end of concentration transient, which was assumed to be sufficiently short (0.1–0.2 ms) with respect to the time interval (1–4 ms) from the moment of stimulation to the starting point of the analyzed fragment of current. P(o, peak) can be expressed as a function of rate constants: P(o, peak) = max(e^Qtp(0)), where p (0) is an initial state probability vector assumed to be zero for all states except for the RG2 state in the case when currents were simulated with 7-state GABA_AR scheme or RL state in the case when currents were simulated with simple 3-state kinetic scheme (see descriptions of both schemes below).

Summing up, ML NSFA can be used for the fast estimation of rate constants, unitary current of synaptic ion channel, the number of synaptic channels bound with a neurotransmitter right after the concentration transient for each synaptic current and peak open probability from the set of macroscopic synaptic currents under Gaussian colored background noise.

We search for the log-likelihood global maximum to obtain the required model parameters from a set of macroscopic synaptic currents. In order to do this, we minimize the negative log-likelihood with a variant of graduated optimization technique using SQP algorithm embedded in fmincon function in MATLAB Optimization toolbox. Initial estimates of each parameter were chosen randomly and uniformly from the logarithmic scale interval, [θ₀/10, θ₀ · 10], where θ₀ is a vector composed of the true values of each parameter (rate constants and unitary current), i.e., of values utilized by the macroscopic current generator (see below). During the search of a minimum, all parameters were bounded within the interval [θ₀/50, θ₀ · 50].

In our version of graduated optimization technique, the whole minimization procedure was divided into sequential minimization steps. On the first step the negative log-likelihood was minimized given the first 2 or 3 currents regularly sampled at 50 points each. On each consequent minimization step the number of points and the number of currents was increased. The parameter estimates, θ_ML, obtained on each previous step were taken as initial parameters θ₀ for each next minimization step. For all calculations in this work each minimization was rerun 5 (3-state scheme) or 10 (7-state scheme) times, each time starting from the different initial parameter set.

First series of simulations of macroscopic synaptic currents were based on experimentally derived 7-state kinetic scheme for GABA_A receptor that had one unliganded state, R, two liganded closed states (RG, singly-liganded and RG2, doubly-liganded) and the respective open (O1 and O2) and desensitized (D1 and D2) states (Mozrzymas et al., 2003). The following rate constants were adapted from Mozrzymas et al. (2003): k_off = 0.13 ms⁻¹, d₁ = 0.14 ms⁻¹, d₂ = 1.5 ms⁻¹, r₁ = 0.02 ms⁻¹, r₂ = 0.12 ms⁻¹, a₁ = 1.5 ms⁻¹, a₂ = 1 ms⁻¹, b₁ = 0.15 ms⁻¹, b₂ = 8 ms⁻¹; k_on1 = 4 ms⁻¹ mM⁻¹, k_on2 = 8 ms⁻¹ mM⁻¹; Unitary current was the same for singly- and doubly-liganded states and was set to 1 pA. Variability in the amplitude of macroscopic postsynaptic responses was achieved by trail-to-trial Gaussian variation of the number of available synaptic channels (mean = 250; SD = 50). Simulation time step was Δt = 0.2 ms. Synaptic vesicle release was modeled as a square pulse of saturating agonist concentration with a duration equal to the single simulation time-step (Δt = 0.2 ms), which caused transition of all available channels from R to RG2 state. A total of 1000 macroscopic synaptic currents were simulated and colored noise that resembled baseline noise of experimentally recorded IPSCs was added to each current. Colored noise (SD = 3 pA) was modeled as a sum of 4 AR(1) processes (Qin et al., 2000; Venkataramanan and Sigworth, 2002):

with parameters ϕ = [0.0067, 0.61, 0.96, 0.999] and σ = [0.32, 1.0, 1.42, 0.72], pA that were obtained from the approximation of autocorrelation function of experimentally recorded (whole-cell patch clamp) background noise by the sum of 4 exponentials (see Equations 23, 24 in Stepanyuk et al., 2011). The decaying phases of the responses (starting in 1 ms after the end of stimulation pulse) were taken for the consequent log-likelihood maximization.

In a second series of simulations we have used simple 3-state kinetic scheme of an abstract synaptic receptor. The scheme consisted of unliganded state, R, singly-liganded state, RL, and open state, O and had the following rate constants: binding rate, k_on = 6 mM⁻¹ms⁻¹, unbinding rate, k_off = 0.025 ms⁻¹, opening rate, b = 0.25 ms⁻¹. The closing rate constant, a, was 2.5 ms⁻¹ for Model R and Model N and 1.25 ms⁻¹ for Model A (see Section ML NSFA Distinguishes Between Changes in the Channel Gating and Changes in the Number of Receptors Bound with a Neurotransmitter in Results). Unitary current was set to 1 pA. Variability in the amplitude of macroscopic postsynaptic responses was achieved by trail-to-trial Gaussian variation of the number of available synaptic channels (mean = 400; SD = 50 for Models R and A; mean = 800; SD = 71 for Model N). Simulation time step was Δt = 0.1 ms. Synaptic vesicle release was modeled as a square pulse of saturating agonist concentration with a duration equal to two simulation time-steps (0.2 ms), which caused transition of all channels from R to RL state.

Accuracy of kinetic rates, unitary current, number of liganded channels, and peak open probability estimates was estimated using bootstrap. To this end, N = 5, 10, 20, 30, 40, or 100 currents were randomly sampled with replacement from the initially generated set of 1000 macroscopic current traces. This procedure was repeated until that 30–40 bootstrap samples were generated. For each bootstrap sample we rerun likelihood maximization m = 5 or 10 times (for currents generated with 3- or 7-state scheme, respectively), starting m-1 times from different randomly generated initial parameter sets (see Section The Log-Likelihood Maximization Procedure above) and mth start was done from θ₀. The estimated model parameters, θ_ML, were obtained from a maximization trial that resulted in the best log-likelihood, which was considered to be a global maximum. The accuracy of estimated model parameters was assessed as a deviation of these parameters (θ_ML) from those (θ₀) used for the generation of the currents, (θML−θ0)2θ0 (hereinafter relative error). The algorithm was implemented in MATLAB.

Accuracy of single-channel current estimates obtained with ML NSFA method presented here was compared to those obtained by PS NSFA. In PS NSFA, variance in currents arising from the stochastic gating of the ion channel is isolated from variance arising from sources such as quantal variability by scaling the mean simulated current waveform to the peak amplitude of each individual simulated current and then subtracting the two waveforms (Traynelis et al., 1993).

To estimate the accuracy of unitary current estimates with PS NSFA, it was applied to n = 1000 bootstrap samples each of which contained either N = 5, 10, 20, 30, 40, or 100 currents simulated with a 7-state GABA_A receptor scheme (see Section Simulation of Macroscopic Synaptic Currents above). For each bootstrap sample the ensemble variance, σ², was plotted against the ensemble mean, 〈I〉, and then fitted with parabola:

where σ²₀ is the variance of the background noise. Accuracy of unitary current estimates was calculated as described above, and was then compared with the accuracy of estimates obtained with ML NSFA. To ensure the best accuracy possible with PS NSFA, the ensemble mean current 〈I〉 and variance, σ², were calculated for each data point and the rising phase of variance vs. mean curve was fitted with parabola using weighted (with weights ω_i = 1/var(σ²_i)) least squares method.

Sampling from a likelihood distribution of model parameters that were estimated from a single macroscopic synaptic current was done by the slice sampling Markov chain Monte Carlo algorithm (Neal, 2003), implemented in “MCMC Methods for MLP and GP and Stuff” toolbox by Toni Auranen and Aki Vehtari (available at http://www.lce.hut.fi/research/compinf/mcmcstuff/).

Postsynaptic architecture restricts direct electrophysiological access to individual receptors in native synaptic environments, with only occasional exceptions when channel openings and closings can be resolved on the very tail of postsynaptic currents (Silver et al., 1992). Both these exceptions and application of PS NSFA (Traynelis et al., 1993) do not allow estimating any parameters of synaptic receptors except their unitary conductance and the number of receptors open at the peak of synaptic current (Hartveit and Veruki, 2006).

In this part of the work we tested how ML NSFA estimates the unitary current and kinetic constants of GABA_A receptors from stochastically simulated macroscopic currents having a trial-to-trial Gaussian variation in the number of available receptors (N_ch = 250 ± 50). Currents were simulated with a 7-state model of this receptor (Mozrzymas et al., 2003, see Methods) having one unbound, two liganded closed, two open and two desensitized states (Figure 1A). Synaptic release of GABA was modeled as a brief (0.2 ms) step of saturating GABA concentration resembling GABA release in real synaptic connections (Perrais and Ropert, 1999, 2000; Hájos et al., 2000; Nusser et al., 2001; Biró et al., 2006; Scimemi and Beato, 2009). 1000 macroscopic currents generated in response to this stimulation had the mean amplitude of 184 ± 35 pA and decay kinetics of 43.4 ± 3.6 ms (Figure 1B) and resembled postsynaptic currents routinely recorded in cortical GABAergic synapses (Nadkarni et al., 2010). Background colored noise (SD = 3 pA, see Section Simulation of Macroscopic Synaptic Currents in Methods) was added to the simulated currents.

Samples consisting of N = 5, 10, 20, 30, or 40 macroscopic currents were randomly chosen from initially generated set of 1000 currents and analyzed using ML NSFA. In order to assess the accuracy of estimates for the unitary current and kinetic rates, parameter search was performed for 60 samples obtained in such a way and log-likelihood maximization was run 10 times for each sample in order to achieve the global maximum (see Section Accuracy of the Estimates in Methods). For each run, the initial parameter values were chosen randomly and uniformly in the logarithmic scale from the range [θ₀/10, θ ₀ · 10], where θ₀ denotes true parameter values, i.e., those used for simulation of currents.

The unitary current was estimated with good accuracy even from samples consisting of only 10 simulated postsynaptic currents (Figure 1C, 8.1% relative error) whereas it was estimated with almost 2-fold better accuracy when the number of currents in the sample was increased from 10 to 40 (4.3% relative error). Three rate constants: unbinding rate (k_off), desensitization (d₂) and resensitization (r₂) rate from double-liganded state could also be estimated (Figures 1D–F). For samples consisting of 10 and 40 currents the relative errors of estimates were: k_off —49.0% and 19.1%; d₂—28.3 and 14.6%; r₂—8.9 and 4.7%, respectively. Some of rate constants associated with single-liganded states were estimated in order of magnitude (a₁) or bounded from below (b₁, d₁).

Thus, we demonstrate that ML NSFA could reliably estimate the unitary current of synaptic receptor channel and several kinetic constants of synaptic receptor model from the very limited number of postsynaptic currents (5–40). These results indicate that ML NSFA may allow analysis of kinetic models of synaptic receptors in their native biochemical environment using routinely recorded macroscopic postsynaptic currents.

The number of currents necessary for a particular algorithm to secure a given accuracy of unitary current estimate is an important practical issue. With many hundreds or even thousands of simulated macroscopic currents accuracy of PS NSFA in estimating the unitary current is fairly good (Markova et al., 2005; Hartveit and Veruki, 2006). At the same time it is hard to collect more than about 100 of evoked postsynaptic currents in steady-state conditions in routine electrophysiological experiments.

Thus, to see whether ML NSFA gives any advantage with respect to the number of required traces we calculated a relative error of unitary current estimates obtained with ML NSFA from the above described samples of different sizes (5, 10, 20, 30, 40, and 100 currents; 60 samples were analyzed in each case to estimate the error) and compared this error with one estimated with PS NSFA applied to 1000 samples of similar sizes.

Figure 2A demonstrates that the error of unitary current estimates obtained with both methods decreases with the number of currents taken for the analysis. However, the unitary current can be estimated with as low as 10.8, 8.1, and 4.9% relative error from only 5, 10, and 20 simulated synaptic currents, respectively, whereas PS NSFA resulted in about 2-fold lower accuracy (23.0, 14.7, and 10.4 relative error, respectively). The estimates obtained with ML NSFA from the analysis of samples of 30 and 40 currents had relative error of 4.6 and 4.3%, whereas PS NSFA gave 8.6 and 7.2% error for these cases. At last, accuracies of unitary current estimates obtained from 100 simulated currents were high for both methods and were comparable (Errors: 2.9% for ML NSFA vs. 4.5% for PS-NSFA; Figure 2A).

Thus, for some complex models ML NSFA allows evaluation of the unitary current with good accuracy using experimentally realistic number of macroscopic currents and substantially outperforms PS NSFA in terms of accuracy when only a few (5–30) postsynaptic currents are available for estimating the unitary current.

PS NSFA was specifically designed to be independent of variations in the number of postsynaptic receptors exposed to neurotransmitter and peak open probability for the sake of more accurate estimation of a unitary current (Silver et al., 1996) from postsynaptic current fluctuations. Unfortunately, this method could not be used for the estimation of the total number of receptors in the synapse. To the contrary, ML NSFA presented here allows estimation of the number of receptors bound with neurotransmitter by the end of neurotransmitter concentration transient in each macroscopic current (liganded channels, N_ch). It is assumed that this transient time course is known or sufficiently brief, meaning that it could be approximated by delta function in the latter case. Indeed, GABA concentration in the synaptic cleft decreases by a factor of 10 during less than 0.1 ms after synaptic vesicle release (Scimemi and Beato, 2009) resulting in almost instantaneous concentration transient. For such a brief concentration transient and for a given GABA receptor model (Figure 1A) ML NSFA would estimate the number of receptors bound with two GABA molecules by the end of concentration transient in all synapses of particular synaptic connection independently upon receptor saturation in the case when most of the current is mediated by the receptors in double-liganded states.

The open probability P(o) at any given time is determined as a mean current divided by a product i_chN_ch, and it is a function of rate constants: P(o) = e^Qtp(0) (see Section Estimation of the Number of Channels and Peak Open Probability in Methods). Thus, P(o) as a function of time is automatically estimated at the end of log-likelihood maximization procedure. The peak open probability is simply a maximum of this function, P(o, peak) = max (e^Qtp(0)). Peak open probability estimated by ML NSFA is, thus far, a ratio of the number of receptors open at the peak of postsynaptic current to the number of double-liganded receptors by the end of neurotransmitter concentration transient.

Figures 2C,D demonstrate that the error of N_ch and P(o, peak) estimates obtained with ML NSFA decreases with the number of currents taken for the analysis. The number of liganded receptor channels, N_ch, was calculated as an average over all currents in the sample and was estimated with 24.5 and 12.4% relative error from samples consisting of only 5 and 10 simulated macroscopic synaptic currents, respectively. The respective estimates of accuracy for the peak open probability, P(o, peak), had 14.4 and 9.8% relative error, respectively. Both N_ch and P(o, peak) were estimated with even better accuracy from samples consisting of 100 simulated currents (10.0 and 4.3% relative error, respectively).

At the same time PS NSFA applied to the same samples resulted in a variance vs. mean curve that was profoundly skewed (Figure 2B, gray dots) and, therefore, could not give an estimate of the number of liganded channels, N_ch.

Most ligand-gated channels are described by kinetic schemes with multiple, non-identical open states often having different conductance levels (Jin et al., 2003; Mozrzymas et al., 2003; Robert and Howe, 2003; Wyllie et al., 2006; Keramidas and Harrison, 2010; Mortensen et al., 2010). In practice some open states should be considered rare and excluded from the fitting of experimental results in order to estimate at least some parameters of receptor kinetic schemes (Mortensen et al., 2010). Unfortunately, PS NSFA is also not applicable to examination of receptors having multiple conductance levels giving values of unitary current and channel number having no obvious physical interpretation (Hartveit and Veruki, 2006). Thus, at the present moment single-channel recordings are virtually the only approach that allows identifying multiple conductance levels of ligand-gated receptors and this approach is also not applicable for studying of synaptic receptors.

We next wanted to investigate if ML NSFA suggested in this work is applicable to analysis of ion channels and ligand-gated receptors with multiple conductance levels, described by kinetic schemes with non-identical open states. 7-state kinetic model of GABA_A receptor (Mozrzymas et al., 2003) having two open states O1 and O2 with identical unitary current (i₁ = i₂, see Figure 1A) was modified to have the unitary current i₁ = 2 pA and i₂ = 1 pA for the states O1 and O2, respectively (Figure 3). Rate constants of the model were modified in such a way that the contribution of single- and double-liganded open states to the total macroscopic current became comparable. Modified constants were (in ms⁻¹): b₂ = 4, b₁ = 1.2, d₁ = 1, r₁ = 1, d₂ = 0.15, r₂ = 1. Colored background noise with SD = 3 pA was added to the simulated currents (Figure 3A, upper panel).

Representative examples of the simulated current components associated with either state O1 or state O2 are shown in Figure 3A, lower panel, by blue and green lines, respectively. When 250 samples consisting of N = 50 simulated currents (Figure 3A, upper panel) were analyzed by PS NSFA the unitary current estimates were close to the unitary current of single-liganded open state O1 (1.86 ± 0.03 pA vs i₁ = 2 pA for the state O1). At the same time, ML NSFA gave reasonable estimates for both conductance levels (Mean ± SE i₁ = 2.00 ± 0.11 pA and i₂ = 0.89 ± 0.08 pA; n = 15 samples of N = 50 simulated currents; Figure 3B). ML NSFA also reliably estimated kinetic rates for single-liganded state transitions (k_off = 0.13 ± 0.01, d₁ = 0.89 ± 0.34, r₁ = 1.02 ± 0.08, a₁ = 1.55 ± 0.05, b₁ = 1.17 ± 0.24 ms⁻¹, Figure 3C) and the mean number of liganded channels (N_ch = 557 ± 53 vs. 500 ± 50 used in simulation).

Thus, contrary to PS NSFA, ML NSFA can reliably estimate kinetic schemes with several open states having different conductance levels and gives precise values of unitary currents, some kinetic rates, and the mean number of liganded receptors in a given synaptic connection.

We next attempted to explore ML NSFA capability to identify which postsynaptic parameters were changed in the case when mean amplitude of simulated currents was increased without changes in macroscopic current waveform and unitary current.

To this end three distinct groups of 1000 macroscopic currents were generated using a simple 3-state scheme of synaptic channel (Figure 4A, see Section Simulation of Macroscopic Synaptic Currents in Methods). A similar increase in mean current amplitude was achieved by changes in either receptor gating or receptor number. A reference kinetic scheme (Model R; Figure 4A, red) had the closing rate, a = 2.5 ms⁻¹ and the total number of channels N_ch = 400 ± 50 and was used to generate a group of macroscopic currents before putative remodeling of synaptic connection (Figure 4B). In the second kinetic scheme (Model A; Figure 4A, blue) mimicking remodeling of receptor gating the closing rate, a, was changed from 2.5 ms⁻¹ to 1.25 ms⁻¹ resulting in almost 2-fold increase of average current amplitude (Figure 4C, blue) without substantial changes in current waveform (Figure 4D, blue vs. red). Conversely, in the third model (Model N; Figure 4A, black) the number of available channels, N_ch, was increased from 400 ± 50 to 800 ± 71 without any changes in the kinetic constants, which led to similar changes in current amplitude (Figure 4C, black) as for Model A without any changes in current waveform (Figure 4D, black vs. red). Therefore, currents generated with Models A and N had similar amplitudes and when normalized, appeared to have the same waveforms as reference currents generated by Model R (Figures 4B–D).

ML NSFA was run with n = 20 bootstrap samples consisting of N = 100 currents (see Section Accuracy of the Estimates in Methods) for each of the 3 groups of simulated currents in order to evaluate the receptor model parameters and the respective errors. Log-likelihood maximization was run 5 times for each bootstrap sample in order to achieve the global maximum. When the parameter estimates obtained from currents generated with Model R were compared to those obtained from currents generated with Model A (Figure 4E, red vs. blue boxes) the difference, Δ_RA, between mean values of each parameter estimates except the closing rate, a, and peak open probability, P(o, peak), was small and was within the standard error (SE) range of the respective estimates: k_off: Δ_RA = 1.7% (SE = 2.2%), b: 8.5% (13.0%), i_ch: 0.6%(2.4%), N_ch: 8.8% (14.0%). At the same time, Δ_RA was 49.7% for the closing rate, a and 70.1% for the peak open probability, P(o, peak) and did not fall within the narrow ranges of the respective SE's (2.4% and 14.1%, respectively). The mean values of the respective estimates were a = 2.51 ± 0.04 for Model R and 1.26 ± 0.03 for Model A, P(o, peak) = 0.08 ± 0.01 for Model R and 0.14 ± 0.02 for Model A. Therefore, we could infer that these were the parameters that altered. These results directly indicate that ML NSFA may reliably determine changes in receptor gating, which leads to an increase in peak open probability.

When estimates obtained from currents generated with Model R and Model N were compared, we observed insufficient differences, Δ_RN, between mean values of all parameter estimates except the number of receptors, N_ch, which was changed from 419 ± 62 for Model R to 942 ± 228 for Model N (Figure 4E, compare red vs black boxes). Δ_RN for N_ch was 124.9% and did not fall within the range of its SE (24.2%). At the same time, Δ_RN for other parameters fell within the respective standard error (SE) range: k_off: Δ_RN = 0.3% (SE = 1.7%), a: 1.3% (1.7%), b: 7.1% (18.0%), i_ch: 0.8%(1.6%), P(o, peak): 8.3% (18.9%) and it was possible to conclude that the number of receptors was the only altered parameter in this case.

Thus, with ML NSFA it becomes possible to distinguish between alteration in receptor channel gating and receptor number, which nonetheless resulted in visually indistinguishable postsynaptic currents.

The key assumption of the PS NSFA is that all receptors in a particular synaptic connection under study have identical kinetic properties (Silver et al., 1996). As a result, all variance in the currents could be attributed to the stochastic nature of the channel gating rather than to the variability in their kinetics. In fact, this assumption could be violated since receptors in the synaptic connection could have different subunit composition or could be differentially modulated (Popescu and Auerbach, 2003) and a set of receptors contributing to each postsynaptic current could vary from trial to trial. In this case PS NSFA overestimates the unitary current and this overestimation could be quite significant even if the difference between receptor kinetic rates is so small that it could be hardly noticed from the observation of synaptic currents (see Figure 5A and below).

Using likelihood approximation it is possible in principle to estimate unitary current and other parameters independently for each individual synaptic current. To test this possibility we have conducted a series of computational experiments. A group of 1000 synaptic currents was simulated using 7-state kinetic scheme of GABA_A receptor channel (Mozrzymas et al., 2003, see scheme in Figure 1A) and the other 1000 currents were simulated using similar scheme in which several parameters (closing rate, k_off, desensitization rate, d₂ and resensitization rate, r₂) varied between trials randomly and uniformly in the range of ±20% of parameter values that were used to generate the first group of currents. In both cases the unitary current was set at 1 pA and colored background noise (SD = 3 pA) was added to the simulated currents (see Section Simulation of Macroscopic Synaptic Currents in Methods for details). Figure 5A demonstrates that both groups of currents had similar waveforms and their decay times were almost identical although variability of decay times in the second group was slightly higher (Mean ± SD: 43.6 ± 3.7 ms vs. 43.9 ± 6.1 ms, N = 1000 currents). Nevertheless, variance vs. mean curves for these two groups of currents differed significantly (Figure 5B) and for the second group unitary current appeared to be 1.9-fold overestimated by PS NSFA (Mean ± SE was 1.01 ± 0.03 pA for the group of currents without variation of parameters vs 1.92 ± 0.05 pA for the group of currents with variation of k_off, d₂, and r₂; N = 250 currents; true value was 1 pA).

To the contrary, when ML NSFA was applied to the group of currents with varying rate constants and log-likelihood of each current in the group was optimized independently, a reasonably accurate estimate of unitary current was obtained (Mean ± SD = 0.89 ± 0.23 pA, N = 50 currents). Standard error of mean unitary current estimate was very low (SE = 0.033 pA), but bias from the true value (1 pA) was significant. We have noticed that the cause of this bias is the skewed shape of the likelihood distribution of a single simulated synaptic current, which means that for the case of single current the maximum likelihood value of unitary current is not the most common value. An example of the typical distribution of unitary current obtained by sampling from the likelihood distribution for a single simulated macroscopic current using the slice sampling Markov chain Monte Carlo method (MCMC, 2000 samples) is shown in Figure 5C (upper panel). It can be seen that the distribution maximum significantly differs from the distribution mean (red vertical line). Therefore, in order to obtain “typical” values of unitary current, mean values of unitary current were also estimated by slice sampling from the likelihood distributions obtained for individual synaptic currents (1000 MCMC samples for each current) from the same group of 50 currents. The resulting distribution of unitary current estimates obtained by accumulation of all 50 distributions for individual currents is represented in Figure 5C (lower panel). The final estimate of unitary current was obtained by averaging over N = 50 mean unitary currents and was in perfect agreement with its true value (Mean ± SD = 0.97 ± 0.39 pA, red vertical line in Figure 5C, lower panel; SE = 0.056 pA). Figure 5D shows statistics of the mean unitary current estimates obtained with MCMC applied to likelihood distributions of individual currents (left box, N = 50) in comparison with the same statistics obtained with PS NSFA applied to individual currents as described above (right box, N = 250, n = 50 bootstraps). It is clearly seen that, contrary to MCMC, PS NSFA significantly overestimates unitary current (green line indicates true value, 1 pA). Among the other model parameters only the number of liganded channels, N_ch, and the resensitization rate, r₂, were estimated with MCMC with relatively high accuracy. The desensitization rate, d₂, and GABA unbinding rate, k_off, were estimated in order of magnitude. The median of the absolute difference between estimates of model parameters and their true values for k_off, d₂, r₂, i_ch, and N_ch was 191, 188, 22, 31, and 35% of their true values, respectively.

We conclude that the mean values for several parameters of the synaptic receptor model, such as the unitary current, the number of channels and the peak open probability, can be estimated with a reasonable accuracy using ML NSFA or MCMC sampling from the likelihood distribution of each individual current in the group of currents even if these currents were mediated by receptors having different kinetic models.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.