A nonequilibrium binary elements-based kinetic model for benzodiazepine regulation of GABA_A receptors

Ion channels are specialized proteins composed of multiple structural domains that regulate the passage of ions across cell membranes. To accomplish this task, they undergo structural rearrangements between conducting and nonconducting states that are initiated or influenced by external signals such as transmembrane voltage or ligand binding. In many cases these signals are known to impinge on structural domains separate from those comprising the channel gate, where they influence channel gating via interactions between domains. Of particular clinical interest is the ability of ion channels, as well as enzymes and many other proteins, to bind ligands that allosterically modulate their activity. This characteristic has made them important drug targets for human disease therapies (Kaczorowski et al., 2008).

Kinetic models have become important tools in the effort to understand how proteins transduce external signals into changes in activity. The standard approach is to posit a limited set of distinct functional states and specific transitions between them (e.g., Goldschen-Ohm et al., 2011). Although the rate constants in any given model encompass the interactions present in the system, without an explicit separation of interactions from the kinetics of distinct functional domains, it can be difficult to interpret whether changes in rate constants reflect changes in interactions between domains or changes in the behavior within a domain. For example, accounting for a localized modulatory stimulus such as drug binding (or any other local perturbation, such as a mutation) with widespread effects on channel function necessitates changing numerous transition rates throughout the model. In contrast, for models where interactions are distinct parameters, one could simply associate the modulator with a local perturbation in a single domain, and then let the interdomain interactions in the system dictate the resulting functional changes in every domain. Furthermore, recent work by several groups has shown that stimulus–response relationships can be used to evaluate the equilibrium energetics of a stimulus pathway (Chowdhury and Chanda, 2012, 2013; Sigg, 2013), observations for which should provide important experimental constraints on model interaction energies.

Numerous clinically important drugs, including barbiturates, benzodiazepines (BZDs), neurosteroids, and general anesthetics, modulate GABAergic inhibition in the brain (Franks and Lieb, 1994). At fast inhibitory synapses these drugs typically lead to a prolongation of inhibitory postsynaptic currents (IPSCs), an effect that has been attributed to drug modulation of postsynaptic GABA_A receptor kinetics because similar changes are observed in excised patches exposed to brief, high-concentration pulses of GABA. The mechanisms by which these effects occur are of considerable interest, and have been investigated using both microscopic single channel and macroscopic ensemble measurements. Different drugs are reported to influence a variety of microscopic kinetic steps, including binding, pore opening, and desensitization (Twyman et al., 1989; Bai et al., 1999). This variety of effects is consistent with the presence of several different drug binding sites inferred through mutational analysis, photolabeling, and molecular modeling.

However, even for a single class of drugs, conclusions from the different approaches have sometimes been contradictory. For example, based upon the observations that the inhaled anesthetics halothane and isoflurane slowed deactivation without altering macroscopic desensitization, we concluded that these drugs slow the microscopic agonist unbinding rate (Li and Pearce, 2000; Benkwitz et al., 2004). However, other general anesthetics, including etomidate and propofol, which bind at the same or nearby transmembrane sites, increase partial agonist efficacy (O’Shea et al., 2000; Topf et al., 2003; Rüsch et al., 2004) and can directly activate receptors, particularly for receptors carrying mutations in the gating apparatus that shift the open/closed equilibrium. These latter observations implicate changes in channel gating. Similar seemingly contradictory observations have been reported for BZDs, which potentiate responses to subsaturating, but not saturating, GABA, and do not affect single channel open dwell times or brief shut times, which is consistent with altered agonist binding (Rogers et al., 1994; Lavoie and Twyman, 1996; Perrais and Ropert, 1999; Goldschen-Ohm et al., 2010). However, BZDs also increase partial agonist efficacy and gate spontaneously open mutants as expected for modulation of channel gating (Downing et al., 2005; Rüsch and Forman, 2005; Campo-Soria et al., 2006; Gielen et al., 2012; Li et al., 2013).

Here, we develop a kinetic model of the GABA_A receptor, a ligand-gated ion channel that is subject to modulation by many different drugs that in clinical practice are often used in combination. Our model explicitly treats the effects of interactions not only on the system’s equilibrium constants, but also on the microscopic transition rates to simulate both equilibrium and nonequilibrium behavior. Such an approach can be crucial for differentiating between possible models based on experimental observations of biologically relevant behavior that for systems such as ion channels often occur far from equilibrium. The model includes two agonist-binding sites that interact directly with two transduction elements, which in turn interact with three gates: one underlying channel activation and two producing desensitization. A BZD-binding site also interacts with one of the transduction elements, but with weaker energetics than the agonist-binding sites.

The model proposed here not only explains the detailed kinetics of activation, desensitization, and deactivation of receptor currents elicited by a variety of stimulus protocols, but also the effect of BZD-positive modulators on all of these kinetic parameters. Furthermore, the model can account for BZD modulation of the affinity and efficacy of peak current dose–response relationships for both full and partial agonists. Thus, this model can explain previous conflicting reports suggesting that BZDs modulate either GABA binding or channel gating, and may represent a more general model of drug modulation at GABA_A receptors.

Human embryonic kidney (HEK) 293 cells were cultured in minimum essential medium with Earle’s salts (Corning) containing 10% bovine calf serum (Sigma-Aldrich) in a 37°C incubator under a 5% CO₂ atmosphere. Cells were transfected using Lipofectamine 2000 (Invitrogen) using the prescribed protocol, with 1–4 µg total of human α₁β₂γ₂ (1:1:1 ratio) GABA_A receptor subunit cDNAs in vector pcDNA3.1 (Invitrogen). Recordings were performed 24–80 h after transfection.

Recordings from outside-out patches excised from HEK 293 cells were made at room temperature using borosilicate glass pipettes filled with (in mM): 140 KCl, 10 EGTA, 2 MgATP, 20 phosphocreatine, and 10 HEPES, pH 7.3, with KOH; osmolarity was 315 mosM, adjusted with sucrose. Patches were voltage clamped at −60 mV and placed in the stream of a multibarreled flowpipe array (Vitrodynamics) mounted on a piezoelectric bimorph (Morgan Electro Ceramics). GABA and diazepam were dissolved in the perfusion solution, which contained (in mM): 145 NaCl, 2.5 KCl, 2 CaCl₂, 1 MgCl₂, 10 HEPES, and 10 dextrose, pH 7.3, with NaOH; osmolarity was 320 mosM, adjusted with sucrose. All reagents were from Sigma-Aldrich. A computer-controlled constant current source (WPI) drove the bimorph to move solution interfaces over the patch with 10–90% exchange times of <200 µs, as measured by the liquid junction current at the open pipette tip after each experiment. Junction currents were generated by altering the ionic strength with an additional 5 mM NaCl or 1% H₂O in solutions containing GABA or diazepam, respectively.

For a given stimulus protocol, current recordings in the presence and absence of diazepam were interleaved for the same patch to allow a direct comparison of current kinetics. The same approach was used for recordings at 30 µM and 10 mM GABA to estimate the channel’s concentration–response relationship. Currents were low-pass-filtered at 5 kHz with a four-pole Bessel filter and acquired at 20 kHz using an amplifier (Axopatch 200B) and digitizer (Digidata 1320A; both from Molecular Devices), controlled by AxoGraph software (Axograph Scientific) running on a Macintosh G4 (Apple Inc.).

A system is modeled as a collection of individual elements, each of which represents a physical or functional feature. For example, in a system describing an ion channel, one element might represent an agonist binding site, another a gate, another a voltage sensor, etc. Each element can exist in either of two metastable configurations, hence their designation as “binary elements.” The configuration of an element denotes its physical or functional status, such as an agonist- or drug-binding site that is unoccupied or occupied, a gate that is closed or open, or a voltage-sensing element that sits in one position or another. For a system of N elements, there will be 2^N states representing all of the possible unique combinations of element configurations.

In the absence of any interactions, each element can flip between its configurations according to intrinsic rate constants γ and δ, which are defined by the height of a single transition state energy barrier between configurations (Fig. 1 A), according to Eyring’s transition rate theory (Eyring, 1935), as

where G₀ and G₁ are the free energies of the element in configuration 0 or 1, respectively; G_‡ is the free energy of the element’s transition state between configurations; k_B is Boltzmann’s constant; T is temperature; h is Planck’s constant; and κ is a proportionality constant.

Binary elements as described above are analogous to subunits in Koshland, Nemethy, and Filmer’s (KNF) description (Koshland et al., 1966) of interacting subunits. Furthermore, each ligand binding step in Monod, Wyman, and Changeux’s (MWC) description (Monod et al., 1965) of cooperative ligand binding can be thought of as an additional binding site element, all of which are coupled to a single element that can exist in either a relaxed or tense configuration.

Interactions between binary elements are state dependent. That is, a particular interaction will only occur when both interacting elements adopt the appropriate configuration. Thus, for each pair of elements there are four possible interactions, one for each unique pair of configurations. Each interaction is associated with an energetic perturbation of not only those metastable states where the interaction is occurring, but also of each of the two interacting element’s transition state intermediates. For example, given an interaction between element i in configuration 1 and element j in configuration 0, we define $ΔGi1j0$ as the energetic perturbation caused by the interaction of those states where elements i and j adopt configurations 1 and 0, respectively. Furthermore, we define $ΔGij0‡$ as the change in free energy of element i’s transition state intermediate due to its interaction with configuration 0 of element j, and $ΔGji1‡$ as the analogous change in element j’s transition state when element i adopts configuration 1 (Fig. 1 B). In the absence of a change in configuration of element j, the effect of its interaction with element i’s transition state is independent of whether or not element i is transitioning from configuration 0 to 1, or 1 to 0 (i.e., interactions with an element’s transition state depend only on the nature of the transition state itself, not on how it got there). Thus, for example, $ΔGij0‡$ will be the same for interactions in both the i₀j₀ and i₁j₀ states (Fig. 1 B). Finally, interactions are assumed to be independent of each other, such that the energies of multiple simultaneous interactions are additive.

For any pair of binary elements i and j, there are four possible states (i₀j₀, i₁j₀, i₀j₁, i₁j₁) and eight possible transitions between states, each of which reflects the change in configuration of one of the elements (Fig. 1 C). In the absence of any interactions, the transition rates will be given by each element’s intrinsic rates γ_i, δ_i, γ_j, and δ_j. When elements i and j interact, each transition rate associated with a change in configuration of element i or j will be modified according to the interaction’s energetic perturbation of both the state from which the transition originates and the transition state intermediate of the element undergoing the change in configuration. For example, the height of the energy barrier governing the transition rate from state i₁j₀ to i₀j₀ will by perturbed by the amount $ΔGij0‡−ΔGi1j0.$ For each element i, adding the appropriate state-dependent interaction energies with every other element j to the free energies of activation in Eqs. 1 and 2 gives the effect of interactions on element i’s intrinsic transition rates for any given state of the system as

where the product over j reflects the sum over any number of simultaneous interactions with each other element j; p ∈{0,1} and q∈{0,1} denote the configuration of element i and j, respectively; $ΔGipjq$ is the interaction energy for those states where elements i and j adopt configurations p and q, respectively; and $ΔGijq‡$ is the energetic perturbation of element i’s transition state when element j is in configuration q. The $θipjq$ terms reflect the effect of the interaction on the system’s metastable states analogous to the equilibrium parameterization used by Chowdhury and Chanda (2010), and the $λijq$ terms reflect the interaction’s effect on the transition state intermediates of the interacting elements. Fig. 1 (B and C) illustrates all of the possible state-dependent interaction factors for a two-element system. However, Eqs. 3–6 are valid for an element involved in any number of simultaneous interactions.

As can be seen in Fig. 1 C, traversing around the two-element cyclic state loop for any pair of elements returns the system to the same energy at which it started, and preserves the ratio of the products of transition rates in either direction around the loop. Thus, interactions as defined above maintain microscopic reversibility by construction. This is true for an arbitrary number of binary elements with any number of interactions.

In the previous section, the system of binary elements was characterized by specifying the intrinsic transition rates γ_i and δ_i for each element i, and all of the pairwise state-dependent interaction factors $θipjq,$$λijq,$ and $λjip$ for each pair of elements i and j. However, this is more information than is required to define the resulting kinetic scheme. That is, the system is “over-parameterized.” This can be understood intuitively by considering an element that interacts in both of its configurations, such that it will never behave according to its intrinsic transition rates γ and δ. Indeed, it is often the case experimentally that the behavior of an element in the absence of any interactions is unknown, and it is sufficient to describe only how each element behaves in some reference state of the system, and how this behavior is altered when interaction energies change upon changes in element configuration.

Without loss of generality, we can characterize a system of binary elements using the minimum number of independent parameters by defining two transition rates α_i and β_i for each element i, and three interaction energies $ΔGij,$$ΔGij‡,$ and $ΔGji‡$ for each pair of interacting elements i and j. We select the state where all of the elements are in configuration 0 as the zero energy reference state. The rate constants α_i and β_i are the rates at which element i flips configurations when all other elements are in configuration 0 (i.e., the transition rates of each element with respect to the reference state; Fig. 2 A). Thus, α_i and β_i are not necessarily the intrinsic rates γ_i and δ_i; rather, they include all of the interactions between element i and configuration 0 of each other element j as

The interaction energy $ΔGij‡$⁠, which imparts a “catalytic” effect, is the difference in the energy barrier governing the transition of element i when element j adopts configuration 1 as compared with configuration 0. $ΔGji‡$ is the analogous change for transition of element j when element i adopts configuration 1. Finally, the interaction energy $ΔGij,$ which imparts the classical “allosteric” effect embodied in the MWC and KNF formulations, is the energy difference between elements i and j both being in configuration 1 as compared with the sum of the energies for either element i or j to be in configuration 1, but not both. These three interaction energies define how the rates α_i and β_i for each element i are perturbed when element j adopts configuration 1 as compared with configuration 0 according to

where the product over j is only over those elements in configuration 1. Fig. 2 A illustrates all of the possible generalized interaction factors for a two-element system. However, Eqs. 9–12 are valid for an element involved in any number of simultaneous interactions. The relationship between the generalized and state-dependent interaction factors for each pair of elements (compare Fig. 1 C and Fig. 2 A) is given by

Fig. 2 B illustrates the energy landscape correlate of the kinetic scheme for a pair of interacting elements i and j as depicted in Fig. 2 A. It is clear that the eight transition rates are completely defined by the relative energies between the four possible states and four transition state barriers between states, where we take the energy of the i₀j₀ state as our zero energy reference. The parameters α_i, β_i, α_j, and β_j define the relative energies of the i₁j₀ and i₀j₁ states as well as the transition state barriers between them and the i₀j₀ state. The transition state barriers between the i₁j₀ and i₀j₁ states and the i₁j₁ state are defined relative to the first two transition state barriers by the interaction energies $ΔGij‡$ and $ΔGji‡,$ and $ΔGij$ sets the energy of the i₁j₁ state relative to the sum of the energies of the i₁j₀ and i₀j₁ states. This is analogous to defining seven of the eight transition rates, where the eighth rate is constrained by microscopic reversibility (Colquhoun et al., 2004).

For a system of N binary elements, there are

possible pairs of interacting elements, and the maximum number of system parameters n_dof-max accounts for two transition rates per element and three interaction energies per pair of elements according to

All of the transition rates between system states are defined by Eqs. 7–12. For a system of N elements, there will be 2^N states, but only N allowed transitions into or out of each state corresponding to the change in configuration of one of the elements (i.e., transitions between states that differ in the configuration of more than one element are disallowed; see Fig. 3). Each system state can be represented by a unique 0-based integer index ∈ [0,2^N − 1] corresponding to the binary representation of each of the elements (Fig. 3 A). For example, “101” is denoted by the index 5 and represents the state of a three element system where the first and third elements are in configuration 1 and the second element is in configuration 0. For a transition from state index u to state index v, the index i of the single element undergoing a change in configuration between states u and v satisfies (1 << i) = (u ^ v), where << and ^ are the bitwise shift and exclusive OR operators, respectively. Furthermore, if (1 << i) & u = 0, then element i transitioned from configuration 0 to configuration 1, otherwise vice versa, where & is the bitwise AND operator. For example, consider the transition from “101” (u = 5) to “111” (v = 7). In this case, the element having the 0-based index i = 1 flipped from configuration 0 to 1, and it holds that (1 << 1) = (5 ^ 7) and (1 << 1) & 5 = 0. Therefore, the transition rate k_uv from state u to state v is given by

where i is the index of the element changing configurations, the product over element indices j is over only those elements that are both in configuration 1 (i.e., (1 << j) & u = 1) and are interacting with the element having index i, and the rate constants and interaction factors refer to Eqs. 7–12.

Channels were modeled as collections of interacting “binary elements” as described above. All model parameters are referenced with respect to the symbol of their associated elements such that α_i and β_i are the rate constants for element i, and Θ_ij, Λ_ij, and Λ_ji are the interaction factors between elements i and j (see Eqs. 7–12). The sets of differential equations describing the time-dependent probability of state occupancies are defined by the transition rates between states as given in Eq. 18, and were solved numerically, either using an eigensystem decomposition of the matrix of transition rates (Colquhoun and Hawkes, 1995), a fourth order Runge-Kutta differential equation solver with error checking and adaptive step sizes (Press et al., 2007), or by Monte Carlo techniques. Programs were written in C++ (see “Availability of software” below).

For each of schemes I–III (Fig. 4 A), the agonist binding sites were coupled directly to the main gate, such that bound agonist speeds main gate opening by the factor Λ_mA. For simplicity, we set Λ_mA = Θ_Am and Λ_Am = 1, which results in no net change to the main gate closing rate upon agonist binding. The main gate closing rate β_m was set to 380 s⁻¹ based on observed single channel open times (Keramidas and Harrison, 2008), and the opening rate α_m was constrained according to α_m = β_mP_o-min/(1 − P_o-min), where P_o-min = 1e⁻⁵ is the channel’s estimated spontaneous opening probability (Campo-Soria et al., 2006). Finally, the interaction energy between each agonist binding site and the main gate was obtained from the constraint Λ_mA²α_m = β_mP_o-max/(1 − P_o-max), where P_o-max = 0.8 is the single channel intraburst open probability in the presence of saturating GABA (Keramidas and Harrison, 2010). Thus, Λ_mA = Θ_Am = 632, which at a temperature of 298.15 K corresponds to ΔG_Am = −3.8 kcal/mol.

We coupled the agonist binding sites directly to either one or both desensitization gates in an analogous fashion to the main gate (i.e., Λ_dA = Θ_Ad and Λ_Ad = 1, where d represents either d_f or d_s). Fits to our experimental data were most sensitive to the combination of desensitization gate transition rates and their energetic coupling, rather than either one of these parameters alone. Thus, we arbitrarily fixed the energy of each interaction with the desensitization gates, and only allowed the rates of desensitization/resensitization to vary. For scheme I, we set $ΘAdf=ΘAds=10,$ which corresponds to $ΔGAdf=ΔGAds=−1.4 kcal/mol.$ For schemes II and III, we set $ΘAds=7.26$ and $ΘAdf=7.26,$ respectively, for a net agonist binding energy of ΔG_Am + ΔG_Ad = −5 kcal/mol (d represents either d_f or d_s), similar to observations of the thermodynamics of ligand binding (Jackson, 1993; Maksay, 1994).

The drug unbinding rate in schemes IV–VI (Fig. 5, A and C; and Fig. 6 A) was set to β_D = 0.3 s⁻¹ based on measurements of the kinetics of diazepam unbinding (Goldschen-Ohm et al., 2010), and the binding rate was estimated as α_D = β_D/K_D = 1 × 10⁸ M⁻¹s⁻¹, where K_D is the affinity constant found in radioligand binding and current potentiation experiments (Sieghart, 1995; Boileau et al., 1998). For scheme VI, the interactions between the intermediate elements and the main and desensitization gates is analogous to those between the agonist binding sites and channel gates in scheme I (i.e., Θ_mF = Θ_mA and Θ_dF = Θ_dA, where d represents either d_f or d_s). The direct interaction with each agonist binding site was constrained such that Θ_AF = 4,590, which corresponds to ΔG_AF = −5 kcal/mol at 298.15 K similar to observations of the free energy of GABA binding (Maksay, 1994). For simplicity, one of the two catalytic interaction factors was further constrained by requiring that Θ_AF = Λ_AFΛ_FA with Λ_AF ∈ [1,Θ_AF], which is equivalent to splitting the interaction energy ΔG_AF between each of the interacting element’s transition state intermediates (i.e., $ΔGAF=ΔGAF‡+ΔGFA‡$⁠). The drug-binding site in scheme VI was coupled to one of the intermediate elements in an analogous fashion to agonist binding (i.e., Θ_DF = Λ_DFΛ_FD), although the interaction energy Θ_DF was allowed to vary.

Models were optimized by minimizing the weighted sum of squared errors between simulated and observed GABA_A receptor currents in response to rapid application of GABA onto outside-out patches from HEK 293 cells expressing human α₁β₂γ₂ receptors (solution exchange <200 µs for accurate measurement of fast receptor kinetics; see “Patch clamp electrophysiology” above). Time points were weighted according to a decaying exponential plus a constant for each period of constant stimulus, such that time points immediately after a change in stimuli (e.g., application of GABA) where currents are changing most rapidly are given the most weight. Current responses to saturating (10 mM) GABA were scaled to a peak open probability of 0.7 based on noise analysis (Goldschen-Ohm et al., 2010), and current responses to 30 µM GABA were scaled based on the relative peak current between interleaved 10 mM and 30 µM GABA pulses for individual patches (see Fig. 4 C). Current responses to 30 µM GABA were scaled based on preliminary model simulations after fitting to 10 mM and 30 µM GABA responses. Optimization involved simultaneously fitting exemplary GABA_A receptor current responses to multiple stimulus protocols, a stringent constraint designed to ensure that the models can describe numerous aspects of channel behavior. Current recordings from a minimum of three separate patches were used for each stimulus protocol.

The software used to simulate and optimize the kinetic models are available as Supplemental material. To get the latest version of the software, we recommend checking for future updates that may be available at http://www.sourceforge.net/projects/modbuilder. The software includes a graphical interface for constructing new models and simulating their activity in response to user-defined stimuli.

The code for the software used to simulate and optimize the kinetic models is provided, as well as a user’s guide and a Mac OS X disk image for simple drag-and-drop installation.

We constructed a simple GABA_A receptor model incorporating two agonist binding site elements (A), a main gate (m), and fast and slow desensitization gates (d_f, d_s), such that all three gates (m, d_f, d_s) must be open for the channel to conduct current. We explored three different sets of interactions between these elements in which at rest the main and desensitization gates are closed and open, respectively, and agonist binding promotes main gate opening and desensitization gate closing for a desensitizing current response (Fig. 4 A, schemes I–III). Schemes I–III each have 2⁵ = 32 states, but assuming equivalent agonist binding sites and applying several constraints based on observations of single channels and agonist binding thermodynamics (see Materials and methods) yields a total of only six free parameters describing the transition rates of agonist binding/unbinding and opening/closing of each desensitization gate. The models were optimized by fitting to observed GABA_A receptor currents in response to rapid application of GABA onto outside-out patches from HEK 293 cells expressing human α₁β₂γ₂ receptors (solution exchange <200 µs for accurate measurement of fast receptor kinetics; see Materials and methods). The final model parameters across patches (mean ± SEM) are given in the legend for Fig. 4.

Each of schemes I–III was able to quantitatively account for all of the stimulus protocols shown in Fig. 4, with nearly indistinguishable fits to the observed current time courses. Fig. 4 B illustrates both simulated and observed GABA_A receptor currents in response to a brief 4-ms pulse of 10 mM GABA, similar to the stimulus conditions at a synapse. During prolonged 500-ms and 3-s pulses of 10 mM GABA at least two kinetic components of desensitization become apparent (Fig. 4, C and D, respectively), which are accounted for in schemes I–III by the fast and slow desensitization gates. The time course of recovery from desensitization was assayed with pairs of 5-ms pulses of 10 mM GABA separated by a variable recovery period (Fig. 4 E). Finally, the agonist concentration dependence of the kinetics of both activation and desensitization was estimated by measuring the peak response to 500-ms pulses of both saturating 10 mM and subsaturating 30 µM GABA in the same patch (Fig. 4 C; 10 mM and 30 µM pulses were interleaved).

A family of simulated responses to 500-ms pulses of 10 nM–100 mM GABA are shown in Fig. 4 F, which describes a peak dose–response relationship with a median ligand concentration of 50.3 µM GABA (Fig. 4 G; see Chowdhury and Chanda [2013] for a description of median ligand concentration). Although the constraint on the dose–response relationship is largely caused by the relative kinetics and peak responses elicited by 10 mM versus 30 µM GABA, the resulting model prediction is consistent with previous observations of the concentration dependence of peak currents in response to rapid application of GABA (<200 µs solution exchange time) over a larger concentration range (e.g., Mozrzymas et al., 2003). Model simulations also predict that increased desensitization during prolonged GABA application results in the prolongation of deactivation after GABA removal, as observed in GABA_A receptors studied using rapid solution exchange methods (Fig. 4 H; Jones and Westbrook, 1995).

Despite having 32 states, open and closed dwell time distributions from Monte Carlo simulations had small numbers of components: one open and three closed (Fig. 4 I). In the presence of a saturating concentration of GABA (10 mM), thereby simulating a “cluster analysis” type of experiment (Keramidas and Harrison, 2008), the single apparent open state was largely governed by the main gate closing rate (β_m⁻¹ = 2.6 ms), whereas the three apparent closed states reflected the opening rates of the main and desensitization gates in the presence of two bound agonists ((α_mΛ_ma²)⁻¹ = 0.66 ms; $βdf−1=11.6 ms;$$βds−1=391 ms$⁠). Burst lengths separated by closures longer than 15 ms were concentration dependent (Fig. 4 J) similar to observations in single channel records (Keramidas and Harrison, 2010).

Although numerous clinically important drugs exert their effects by modulating ligand-gated ion channels (Kaczorowski et al., 2008), it is often unclear whether these drugs function by altering channel gating or agonist binding, both of which contribute to shaping postsynaptic current kinetics. In general, relating drug effects on observed physiological responses to modulation of receptor microscopic properties depends on the structure of the kinetic model used. Thus, we explored the predicted effects of drug binding in scheme I (Fig. 4 A) for either modulation of the main gate (Fig. 5 A, scheme IV) or an agonist binding site (Fig. 5 C, scheme V). In scheme IV, the drug binding site (D) was coupled to the main gate analogous to an agonist binding site, albeit with different energetics. In scheme V, drug binding directly affected the agonist binding rate and the stability of the bound agonist. See the legend for Fig. 5 for parameter values.

Drug binding in both schemes IV and V had nearly identical effects on receptor kinetics, slowing deactivation without altering desensitization (Fig. 5, B and D). However, drug binding in schemes IV and V had different effects on the channel’s peak dose–response relationship (Fig. 5 E), knowledge of which can potentially discriminate between the two schemes (Wagner et al., 2004; Goldschen-Ohm et al., 2011). Although peak dose–response curves were left-shifted in the presence of drug for both schemes IV and V, only for scheme IV did drug binding increase the peak open probability in the presence of saturating GABA. Also, only scheme IV allows for drug binding to promote channel activation in the absence of agonist (i.e., increase the resting state open probability).

BZDs are an important class of GABA_A receptor modulators that bind at an interface between subunits (Pritchett et al., 1989; Wieland et al., 1992; Boileau and Czajkowski, 1999). They are used in clinical practice for procedural sedation and to treat a wide variety of conditions including disorders of mood and anxiety, epilepsy, and many others (Rudolph and Knoflach, 2011). However, their kinetic effects on the channel remain a matter of contention, with multiple studies reporting nearly exclusive changes in either binding (Vicini et al., 1987; Rogers et al., 1994; Lavoie and Twyman, 1996; Mellor and Randall, 1997; Perrais and Ropert, 1999; Goldschen-Ohm et al., 2010) or gating (Downing et al., 2005; Rüsch and Forman, 2005; Campo-Soria et al., 2006; Gielen et al., 2012; Li et al., 2013), the totality of which cannot be explained by either scheme IV (Fig. 5 A) or V (Fig. 5 C). For example, scheme IV predicts that drug binding will increase the peak open probability in response to a saturating agonist, whereas BZDs do not potentiate responses to saturating GABA. In contrast, drug binding in scheme V will not enhance the peak open probability at a saturating concentration of either a full or partial agonist, whereas BZDs can potentiate the peak open probability in response to saturating partial agonists, as opposed to full agonists.

To reconcile these seemingly disparate observations, we extended scheme I (Fig. 4 A) to allow both GABA and BZD-positive modulator binding to similarly influence the configuration of an intermediate element involved in transducing ligand binding to channel gating, although with different energetics for GABA versus BZDs (Fig. 6 A, scheme VI). Here, the ligand-binding sites and channel gates do not directly interact, but instead ligand binding promotes activation of an intermediate transduction element (F), which is coupled to the main and desensitization gates analogously to how the agonist-binding sites were so coupled in scheme I (Fig. 4 A). We treat BZD-positive modulator binding at the drug-binding site element (D) in an analogous fashion to agonist binding, albeit with a different energy. Scheme VI has 2⁸ = 256 states, but applying the same constraints used for scheme I and constraining the agonist binding energy at 298.15 K to ΔG_AF = −5 kcal/mol, similar to observations for the thermodynamics of GABA binding (Maksay, 1994), yields only 11 free parameters (see Materials and methods). The final optimized model parameters based on fits to observed GABA_A receptor currents in response to rapid application of GABA (solution exchange <200 µs; see Materials and methods) across patches (mean ± SEM) are given in the legend for Fig. 6.

Drug binding in scheme VI stabilizes one of the intermediate elements in its activated configuration, thus slowing deactivation kinetics as observed in the presence of 1 µM of the BZD-positive modulator diazepam after both brief (4 ms) and prolonged (500 ms) pulses of saturating (10 mM) GABA (Fig. 6 B), similar to observations for inhibitory postsynaptic currents (IPSCs; Mellor and Randall, 1997; Perrais and Ropert, 1999). As observed during 500-ms pulses of 10 mM GABA, desensitization is not affected by bound drug (Fig. 6 B, right). Consistent with previous reports, diazepam potentiated peak current responses to rapid application of subsaturating GABA (10 µM; Fig. 6 C), but not saturating GABA (10 mM; Fig. 6 B), observations that were previously interpreted to suggest predominant changes in binding only (Lavoie and Twyman, 1996; Mellor and Randall, 1997; Perrais and Ropert, 1999; Goldschen-Ohm et al., 2010). In scheme VI (Fig. 6 A), this reflects the relatively much greater ability of bound GABA to activate its associated intermediate element than bound drug (ΔG_DF = −1.7 kcal/mol at 298.15 K), such that drug binding while GABA is bound provides very little additional enhancement of current activation. In contrast, for subsaturating GABA concentrations drug binding allosterically promotes GABA binding at unoccupied GABA binding sites (and stabilizes GABA once bound) via its stabilization of the activated intermediate element. The simulated peak open probability dose–response relationship was left shifted from a median ligand concentration of 30.8 µM to 14.7 µM GABA in the presence of 1 µM drug (Fig. 6, D and E), similar to previous observations of the concentration dependence of BZD potentiation (Downing et al., 2005; Campo-Soria et al., 2006; Gielen et al., 2012; Li et al., 2013).

Drug binding in scheme VI can also account for the observation that BZD-positive modulators increase the efficacy of partial, but not full agonists (e.g., GABA; Fig. 6 E; Downing et al., 2005; Rüsch and Forman, 2005; Campo-Soria et al., 2006; Gielen et al., 2012; Li et al., 2013). Here, partial agonism was modeled by decreasing the interaction energy between the agonist binding sites and intermediate elements to ΔG_AF = −3 kcal/mol at 298.15 K, although similar effects could be obtained by a combination of reductions in binding energy ΔG_AF and binding rate α_A. In scheme VI (Fig. 6 A), the drug-enhanced efficacy of partial agonists reflects the decreased disparity between drug and agonist binding on activation of the intermediate element. Scheme VI (Fig. 6 A) also allows for drug binding to influence the channel opening probability in the absence of bound agonist by directly promoting activation of the associated intermediate element (Fig. 6 F), observations previously ascribed to exclusive effects on channel gating (Downing et al., 2005; Rüsch and Forman, 2005; Campo-Soria et al., 2006; Gielen et al., 2012). Thus, the major observations for BZD-positive modulation of GABA_A receptors including detailed kinetics, direct gating, and enhanced efficacy of partial agonists can be explained by positing that BZD-positive modulators themselves act as weak partial agonists whose effects on downstream gating elements are additive to those of agonist binding at one of the two agonist-binding sites.

The models shown in Figs. 4 A and 6 A reproduce many aspects of both equilibrium and nonequilibrium GABA_A receptor behavior. In addition, scheme VI (Fig. 6 A) accounts for the effects of drug modulation by BZD-positive modulators on both the detailed kinetics of receptor behavior in response to multiple stimulus protocols, as well as dose–response relationships for both full and partial agonists, and also allows drug binding to itself modulate the open/closed equilibrium as observed in spontaneously open mutants. Despite having 256 states, by postulating a model “structure” based upon our intuition of how the elements interact with each other, and applying experimental constraints on both the rate constants and interaction energies as described, we are left with only 11 free parameters. This compares favorably with other kinetic models of GABA_A receptor gating, most of which incorporate between 10 and 24 parameters and do not include drug binding or modulation (Jones and Westbrook, 1995; Haas and Macdonald, 1999). However, our models remain oversimplified. For example, single channel recordings exhibit more open/closed states than are apparent in our simulated dwell time distributions, and also show a redistribution of channel openings among observed open states across agonist concentrations (Haas and Macdonald, 1999; Keramidas and Harrison, 2008). Nonetheless, schemes I–III provide a starting point for more complicated models such as scheme VI, which provides an intuitive mechanistic description of drug modulation that can explain the kinetic changes underlying drug action at synapses.

Neither direct modulation of either the main gate or one of the agonist binding sites was alone sufficient to explain the totality of observations for BZD-positive modulators (Fig. 5, schemes IV and V). In contrast, allowing drug binding to modulate an intermediate element involved in transducing ligand binding to channel gating accounts for both equilibrium and nonequilibrium observations for BZD-positive modulators, including all of the effects previously ascribed to exclusive changes in either binding or gating (Fig. 6 A, scheme VI). Scheme VI posits a direct interaction between the drug-binding site and only one of the two intermediate transduction elements based on the proximity of BZD and GABA binding sites in a structural model (Bergmann et al., 2013). Interestingly, a similar conclusion that BZDs directly modulate only one of the two GABA binding events was reached based on single channel observations (Macdonald and Olsen, 1994).

The inclusion of an intermediate element between agonist binding and channel gating builds upon a previous model proposed by Gielen et al. (2012), hereafter referred to as the GLS model, where BZDs influence a preactivation step intermediate to ligand binding and channel gating. There are several significant enhancements and important differences between scheme VI and the GLS model. First, the GLS model was parameterized for equilibrium conditions only, whereas scheme VI describes both the equilibrium and nonequilibrium kinetic effects of BZD-positive modulators. Second, the GLS model, like the allosteric schemes on which it was based (Downing et al., 2005; Rüsch and Forman, 2005; Campo-Soria et al., 2006), does not consider desensitization, a potentially confounding complication when comparing model simulations to peak current observations, which are shaped by both activation and desensitization processes. Third, for scheme VI, BZD-positive modulator binding affects both agonist binding and gating via the propagation of interactions from the BZD binding site to both the agonist-binding site and channel gates. In contrast, BZD binding in the GLS model does not affect the agonist-binding step itself, but only a subsequent gating step. Because peak current is a nonequilibrium measure depending on both binding and gating, additional evidence such as a measure of binding alone is necessary to distinguish between these two model predictions. Finally, the GLS model postulates different sets of gating equilibrium constants for agonists versus partial agonists, whereas scheme VI provides an explanation based exclusively on differences in interactions between ligands and their binding sites.

Scheme VI (Fig. 6 A) predicts that GABA binding should reciprocally stabilize BZD-positive modulator binding via the stabilization of the activated intermediate transduction element. Such reciprocal interactions between the BZD and GABA binding sites have been observed, as changes in the rates of accessibility of an introduced cysteine residue to methanethiosulfonate reagents at either the GABA or BZD binding site upon ligand binding at the other site (Kloda and Czajkowski, 2007; Sharkey and Czajkowski, 2008).

In summary, we provide a GABA_A receptor model that describes both the detailed kinetics of receptor behavior and the mechanism of BZD-positive modulation in the presence of intact desensitization. Given that the intravenous anesthetic etomidate has been similarly shown to increase partial agonist efficacy and directly activate the channel at high concentrations (Rüsch et al., 2004), this model may also be useful for describing the mechanism of various anesthetics. Finally, the binary elements allosteric-kinetic model representation used here provides a comprehensive approach for building and constraining nonequilibrium kinetic models incorporating multiple allosteric sites. Such models are becoming increasingly important for understanding the behavior of systems such as ion channels that are subject to modulation by numerous allosteric stimuli, including both physiological modulators and therapeutic agents.

This work was supported by the American Epilepsy Society and the Lennox Trust Fund (to M.P. Goldschen-Ohm), the Betty J. Bamforth Research Professorship, the Ralph M. Waters Distinguished Professorship (to R.A. Pearce), National Institutes of Health grant GM55719 (to R.A. Pearce), NIH grant NS046378 (to M.V. Jones), and the Department of Anesthesiology at the University of Wisconsin-Madison. A. Haroldson is presently at Dart Neuroscience LLC.

The authors declare no competing financial interests.

Author contributions: R.A. Pearce conceived the study, M.P. Goldschen-Ohm and R.A. Pearce developed the binary elements model parameterization, M.P. Goldschen-Ohm and R.A. Pearce developed the simulation software, and M.P. Goldschen-Ohm provided the graphical interfaces to the software. R.A. Pearce performed initial simulations for a precursor to scheme I, and M.P. Goldschen-Ohm performed all subsequent simulations for schemes I–VI. M.P. Goldschen-Ohm and R.A. Pearce wrote the manuscript. M.P. Goldschen-Ohm and A. Haroldson performed the experiments in the laboratory of M.V. Jones.

Kenton J. Swartz served as editor.