An Electrostatic Funnel in the GABA-Binding Pathway

The γ-aminobutyric acid type A receptor (GABAA-R) is a major inhibitory neuroreceptor that is activated by the binding of GABA. The structure of the GABAA-R is well characterized, and many of the binding site residues have been identified. However, most of these residues are obscured behind the C-loop that acts as a cover to the binding site. Thus, the mechanism by which the GABA molecule recognizes the binding site, and the pathway it takes to enter the binding site are both unclear. Through the completion and detailed analysis of 100 short, unbiased, independent molecular dynamics simulations, we have investigated this phenomenon of GABA entering the binding site. In each system, GABA was placed quasi-randomly near the binding site of a GABAA-R homology model, and atomistic simulations were carried out to observe the behavior of the GABA molecules. GABA fully entered the binding site in 19 of the 100 simulations. The pathway taken by these molecules was consistent and non-random; the GABA molecules approach the binding site from below, before passing up behind the C-loop and into the binding site. This binding pathway is driven by long-range electrostatic interactions, whereby the electrostatic field acts as a ‘funnel’ that sweeps the GABA molecules towards the binding site, at which point more specific atomic interactions take over. These findings define a nuanced mechanism whereby the GABAA-R uses the general zwitterionic features of the GABA molecule to identify a potential ligand some 2 nm away from the binding site.


Introduction
The neurotransmitter γ-aminobutyric acid (GABA) is the brain's major inhibitory neurotransmitter, which binds to the GABA type A receptors (GABA A -Rs). These GABA A -Rs are 'Cysloop receptors' in the pentameric ligand-gated ion channel (pLGIC) superfamily. Cys-loop receptors are so named due to a well-conserved 13-residue loop that is formed between two cysteine (Cys) residues that are connected via a disulfide bond. Upon agonist (GABA) binding, the channel of the GABA A -R opens and increases the intraneuronal chloride ion concentration, hyperpolarizing the cell and inhibiting transmission of the nerve action potential.
GABA A -Rs are heteropentamers that are composed of many different combinations of distinct subunit gene products (α 1-6 , β 1-3 , γ 1-3 , δ, ε, π, and ρ 1-3 ). While the most common GABA A -R subtype in the brain is the α 1 β 2 γ 2 combination (comprising two α 1 -subunits, two β 2 -subunits, and a single γ 2 -subunit) [1,2], this study will focus entirely on a minor subclass of GABA A -Rs that contain a δ-subunit instead of a γ-subunit. The δ-containing GABA A -Rs comprise only 5-10% of the total GABA A -Rs in the brain [3]. They are mostly located away from the synapses [4,5] and are thought to be involved in the constantly active 'tonic' GABAergic current [6,7]. While only comprising a fraction of the total GABA A -Rs in the brain, the α 6 β 3 δ receptor is one of the most highly GABA (and ethanol) sensitive receptors [8], making it the ideal GABA A -R subtype for studying ligand binding.
Existing structural and biochemical data show that the GABA A -R subunits combine to form an ion channel through the membrane via a pore down the center of the pentamer. Currently, no experimental structure of a heteropentameric GABA A -R is available. A structure for a homomeric GABA β 3 pentamer has been recently released [9], but despite being the first (and so far, only) high-resolution structure resolved, it is a non-physiologically occurring construct. This β 3 pentamer possesses the same structural architecture as described in previous extensive comparison studies [10,11]. Each monomer is comprised of three domains; the extracellular ligand-binding domain (LBD) is comprised of a 'β-sandwich' structure; the transmembrane (TM) domain is composed of four helices; and a cytoplasmic domain of relatively unknown structure forms between TM helices 3 and 4. The LBD of each subunit consists of a 'principal' (+) and 'complementary' (-) side. The GABA binding site is formed by a cleft between the αand β+ subunit interfaces [12,13]. The specific GABA binding site within this cleft has been well defined, with explicit key residues identified in experimental literature (α 6 F65, α 6 R132, β 3 L99, β 3 E155, β 3 R196, β 3 Y205, and β 3 R207) [14][15][16][17][18][19][20]. We have previously detected these residues as critical to GABA binding in our own computational studies [21]. While most of the αand β+ subunit interfaces form the 'sides' of the binding pocket, the β-subunit C-loop folds over the top of the pocket to act as a 'cap' or 'roof' to the binding site. The movement of this Cloop has been linked to the activation mechanism of Cys-Loop receptors [22], and a closed Cloop has long been thought necessary for obtaining the active state of Cys-Loop receptors [23][24][25]. Throughout the rest of the paper we will refer to the region of the LBD that is closer to the TM helices than the C-loop as being 'below' the binding site, and the region of the LBD that is further from the TM helices than the C-loop as being 'above' the binding site. A recent gating model in pLGICs suggests that the activation proceeds via a "conformational wave" that starts in the ligand-binding site (notably loops A, B, and C), and then propagates to the LBD/TM interface (the β 1 -β 2 loop and the Cys-Loop) and finally moves to the TM helices (firstly the M2 helix) to cause the ion pore to open [26,27]. This activation model was validated further with coarse-grained normal-mode analyses of the ELIC and GLIC structures to calculate a closed to open state transition pathway [28]. Conversely, one of the latest computational studies [29] observed a pLGIC open to closed transition pathway upon agonist unbinding. The agonist unbinding was mediated by opening of the C-loop and caused a significant reorientation of the β-sandwiches in the LBD that tilted outward. This rearrangement lead to the β 1 -β 2 loop repositioning at the LBD/TM domain interface, resulting in an inward displacement of the M2-M3 loop and an inward tilting of the pore-lining helices to shut the ion pore. This study also suggests that C-loop opening is necessary but not sufficient for agonist unbinding [29].
Consequently, details about the agonist binding site are known, as well as how agonist binding triggers the activation mechanism that leads to channel opening. However, the vast majority of the critical residues that line the binding pocket are obscured by the C-loop in the available structures, leaving very little of the 'binding site' exposed. Thus, questions still remain as to how GABA recognizes this binding site, and what is the pathway by which GABA gets into the binding site.
In order to investigate these open questions, we attempt to map the binding pathway of GABA using unbiased molecular dynamics (MD) simulations. Using extensive MD simulations, GABA molecules are isolated and a consensus pathway to binding the receptor is determined. Fundamental driving forces that control the binding pathway are identified and analyzed, revealing the need to expand our scope for what we consider as important residues and regions for ligand-binding.

Categorization of simulations
The initial criteria used for defining the 'binding state' of GABA in the simulations is the distance between the center of mass (COM) of the GABA molecule and the COM of those residues that line the binding pocket (α 6 F65, α 6 R132, β 3 L99, β 3 E155, β 3 R196, β 3 Y205, and β 3 R20). In order to spatially define when a GABA molecule is within the binding pocket, we ran analysis on a control simulation, where GABA is directly docked into the binding site. This 20 ns simulation shows that the GABA molecule equilibrated to a consistent position within the pocket, with its COM~0.62 ± 0.06 nm from the COM of the binding site. Thus, we define a GABA molecule that has a COM <0.70 nm from the COM of the binding site as being in the pocket. Given that the GABA molecule itself is~0.65 nm long, when the GABA COM is within 1.3 nm of the binding site COM, then at least part of the GABA molecule is within the region of the binding pocket. When the GABA COM is within 1.0 nm of the binding site COM, then the majority of the GABA molecule is within the binding pocket region.
Using this metric for the 'binding state' and cut-off values to indicate binding in the pocket, the simulations fall into four distinct groups (Fig 1). The largest group of simulations is when GABA is 'NON-BINDING' (73 simulations)-defined as those that never get closer than 1.3 nm to the binding site COM. Eight of the simulations are when GABA binds 'NEARBY'defined as GABA getting partially within the binding site (< 1.3 nm) without fully reaching within the binding site (> 0.70 nm). Finally, nineteen simulations that show GABA 'binding'defined as the GABA COM getting closer than 0.70 nm from the binding site COM-are subdivided into two categories: 1) 'BIND' (nine simulations) where GABA remains within the binding site, and 2) 'PARTIAL' (10 simulations) where GABA reaches the 0.70 nm cut-off, but then moves away again. Through the remainder of this article, we refer to the simulations as 'NON-BINDING', 'NEARBY', 'PARTIAL', and 'BIND'. All simulations in 'NEARBY', 'PAR-TIAL', and 'BIND' are used in subsequent analysis (as described in the Methods). A random subset of 20 of the 73 simulations is chosen as a 'NON-BINDING' sample to analyze the non-binders.
It is important that the starting positions of the GABA molecules do not influence their binding, and confirmation of independent starting positions is vital. Thus, three methods were used to verify that the initial starting positions of the GABA molecules in our simulations did not bias the outcome of the results. Firstly, the average starting positions of all the GABA molecules within each category were calculated and are represented relative to the protein and the binding site COM (Fig 2). Not only do the average starting positions all occupy a similar location, but the standard deviations of the GABA position in each category are similar in size and completely overlap one another. This indicates that the quasi-random initial starting positions of the GABA molecules had no bias on the molecule's probability to enter the binding site. Secondly, the average initial positions are all approximately level with the 'height' (z-coordinate position) of the binding site COM, as well as closely aligned with the vector from the binding site COM directly away from the protein; indicating no bias in any particular direction relative to the binding site. Thirdly, the initial movement of the GABA molecules was measured as a movement towards or away from the protein (S1 Fig). This movement shows no bias towards the protein and displays random behavior. All three methods confirm that the initial position of the GABA molecules in each simulation do not bias the results.

Average binding pathway
To determine the distinct pathway for GABA to bind to the GABA binding site, average positions of GABA relative to the protein within each simulation category are calculated (as described in Fig 3A). These average GABA positions indicate that all of the 'binding' (BIND, The 100 independent GABA simulations are broken down into four categories based upon their distance from the GABA binding site. Those that enter the binding site (< 0.70 nm from the binding site COM) and remain there fall into the 'BIND' category (red). Those that enter the binding site and leave again fall into the 'PARTIAL' category (orange). Those that partially enter the binding site (< 1.3 nm from the binding site COM) fall into the 'NEARBY' category (green). Those that do not even partly enter the binding site (never < 1.3 nm from the binding site COM) fall into the 'NON-BINDING' category (blue). The radial representation on the right shows the populations of each of the categories (with the corresponding color) and the relative proximity they reach to the binding site COM (shown as a red circle). The number next to the category is the number of simulations that are in that category. PARTIAL, NEARBY) simulations follow a similar pathway ( Fig 3B) with comparable characteristics: 1) GABA approaches the binding site from the membrane side of the C-loop (below) before reaching the protein and then moves ('flips up') up into the binding site, and 2) the distribution of the GABA positions within the pathway is very narrow. By contrast, the NON-BINDING simulations appear to have GABA positioned further from the membrane, 'above' the binding site and have a much wider distribution of GABA positions with larger standard deviations.
Although a consistent pathway for GABA molecules to approach the binding site from the membrane side of the C-loop was identified, further analysis shows this pathway is determined by influences of the protein. A theoretical non-biased "random" distribution of GABA locations was artificially generated and used to calculate the hypothetical positional standard deviations expected if GABA molecules were to approach the binding site in a completely random manner. This standard deviation data calculated from an artificially generated random distribution was used as a metric for random/non-biased GABA dispersal. The NON-BINDING simulations do indeed have GABA distributions that are comparable to these hypothetical 'random' distributions ( Fig 4A). The GABA molecules in these NON-BINDING simulations begin to adopt this quasi-random and mostly-unbiased distribution behavior once they reach a distance of 2.7 nm from the binding site COM (red line, Fig 4A). At this distance, the GABA molecules are beyond the forcefield cut-off distance to be influenced by van der Waals interactions with the protein. As such, these GABA molecules are essentially experiencing bulk solution-like conditions with only weak, long-range electrostatic effects from the protein and thus are more randomly distributed. When the hypothetical 'random' distribution is compared to the average distributions seen in all the binding (BIND, PARTIAL, and NEARBY) simulations, the pathway by which GABA molecules approach the binding site is indeed far from random, as shown in Fig 4B. The presence of the protein greatly influences the GABA molecules to become more concentrated at specific positions as they approach the binding site. In particular, the GABA molecule positions are definitively 'focused' or 'funneled' once they reach a distance of~1.9-2.0 nm from the binding site COM (red line, Fig 4B). Furthermore, deconstructing the average GABA positions into the angles they make with both the XY-plane and the YZplane shows that for all three of the binding categories (BIND, PARTIAL, and NEARBY), both the XY-plane angle and the YZ-plane angle converge to almost the same value when the GABA molecules are at a distance of~1.9 nm from the binding site COM (XY-plane angle =~95°, YZ-plane angle =~110°) ( Fig 5). Thus, at a distance of 1.9-2.0 nm from the binding site COM, the GABA molecules of all the binding simulations converge to this same point and remain narrowly distributed as they approach the protein past this location. We have termed this The positions of all the GABA molecules within a particular bin (ii) are used to calculate an average position and standard deviation of GABA at that particular distance from the binding site COM (iii). These data are represented graphically as a disc with a radius that is proportional to the standard deviation of the GABA positions within that bin, and with its axis aligned to the vector between the average GABA position and the binding site COM (iv). (B) The average binding pathways calculated for each of the simulation categories (L-R: BIND, PARTIAL, NEARBY, NON-BINDING) are shown. The radius of the disc is proportional to the standard deviation of GABA molecule positions at that distance from the binding site. The ligand-binding domain of the β 3 -subunit is also illustrated in a grey cartoon format. A red circle represents the binding site COM. The discs used for these figures were constructed by using the draw feature of VMD to create cylinders with defined centers along the vector between the ligand binding site COM and the ligand position COM. The radius of the cylinder is defined as the measured standard deviation at that position.

Electrostatics drive the binding pathway
Given that the GABA molecule is zwitterionic, with charged termini connected by a short hydrocarbon linker, the driving forces behind the binding pathway may be electrostatic in nature. We calculated the electrostatic potential surface of the protein to test this hypothesis. These data are also used to construct a representation of the electrostatic field that surrounds the GABA receptor ( Fig 6). Visualization of the field lines provides an intuitive approach to identify the regional intensity and gradient of electric fields in relation to the GABA A -R structure. In other protein systems, such as acetylcholinesterase, the field lines around the protein are often used to interpret binding mechanisms for the positively charged acetylcholine [30,31].
Analysis of the electrostatic field lines around the GABA receptor reveals that the strongest, most persistent areas of the field converge at two regions of the GABA receptor (Fig 6A-6C) that correspond to the two GABA binding sites at the α+βsubunit interfaces. Specifically, the electrostatic potential surface shows a highly electronegative region in the α+βcleft just below the GABA binding site (Fig 6D), where the electrostatic field lines converge. As a comparative assessment, this same electronegative region has been observed in previously published GABA models for both the α 6 β 3 and α 1 β 2 clefts [21,32,33], as well as test models constructed using the recently published GluCl [34] and GABA β 3 [9] homopentamer crystal structures (S2 Fig). The influence of the electrostatic field on GABA is measured by calculating the net strength and direction of the dipole on the GABA molecules as they approach the binding site (Fig 7). All of the GABA binding pathways follow the electrostatic field lines with the dipoles of the GABA molecules strongly aligned within the field. It appears that the GABA molecules are being 'swept along' within the electrostatic field. By stark contrast the GABA positions in the NON-BINDING simulations do not appear to correlate with the electrostatic field lines, and  Fig 7B). At these positions the GABA molecules are randomly oriented, and as such any dipole on the molecules will cancel out.
By examining the overall properties of the GABA molecules, we observe that a large change in behavior occurs at the pathway 'midpoint'~1.9-2.0 nm from the binding site COM. It is at Finally, the electrostatic surface of a β 3 -subunit and a α 6 -subunit is presented to illustrate the dense electronegative region just below the binding site. The binding site is highlighted by a white dashed circle, and the electronegative region is highlighted by a red circle.

Statistical analysis of the pathway
One of the critical junctures in our pathway is the midpoint~1.9-2.0 nm from the binding site COM, where the GABA molecules enter the electrostatic field. To emphasize the essential importance of this point in the pathway, the distance between GABA and the midpoint was measured for each of the 100 independent simulations (S4 Fig). GABA molecules are considered as having reached the midpoint when the distance between the GABA COM and the midpoint coordinates is <1.0 nm (as the midpoint is actually a 'region' with a variability of~0.5-0.6 nm). Using this metric, 28 of the 100 simulations reach the midpoint, whereas 72 do not. Of the 28 simulations that reach the midpoint, 24 of them (86%) go on to enter the binding region. Of the 72 simulations that do not reach the midpoint, only 3 enter the binding region. Thus, of the 27 total simulations (BIND, PARTIAL, and NEARBY) where GABA enters the binding region, 24 (89%) proceed via the midpoint, indicating that the midpoint is indeed a crucial checkpoint that must be reached before progressing to the binding site.
Further analysis was carried out to look at the probability of GABA reaching the midpoint. The midpoint is~1.9-2.0 nm from the binding site COM, but~2.5 nm from the average GABA starting position. Thus, to test the initial dispersal of the GABA molecules, each of the simulations was monitored until the GABA molecule reached a distance of 2.5 nm from the average starting position (Fig 8A-8C). The coordinates at which the GABA molecules first reach this '2.5 nm spherical shell' were measured and cross-referenced with those that are within 1.0 nm of the midpoint (Fig 8D). Of the 100 simulations, six GABA molecules first reach the 2.5 nm shell within the vicinity of the midpoint. Thirty-two equally distributed . The red hashed area is devoid of GABA molecules. This is due to the region being occluded by part of the GABA A -R itself (C), specifically, the α 6 -subunit that comprises part of the binding site. This α 6 -subunit is shown in dark grey. The relative position of the GABA binding site COM is illustrated as a red circle, and the rest of the protein is shown in light grey, from a viewpoint looking down on the protein. Thus, excluding this protein-occluded region, the calculated accessible area on this 2.5 nm shell is~73 nm 2 . The number of GABA molecules that first reach the 2.5 nm shell within 1.0 nm of the midpoint was calculated (D). The binding site COM is shown as a red circle, and the β 3 -subunit of the binding site is shown in grey. We are essentially measuring the number of GABA molecules within a specific circle of radius 1 nm on the surface of a sphere of radius 2.5 nm. 32 additional, 'random' overlapping sampling points were measured (E), and represented as colored patches on the surface of the 2.5 nm sphere (F and G). The chance of a molecule randomly reaching any sample point is~4.3%. This is the area of sample point (approximately a circle of radius 1 nm-3.14 nm 2 ) as a fraction of the available area (calculated as~73 nm 2 ). The numbers measured for the sample points (4.55 ± 1.80) portrayed in (F) and (G) indicate that there is no preference for the GABA molecules to initially travel to the midpoint, and that the positions of the GABA molecules are consistent with a random distribution. doi:10.1371/journal.pcbi.1004831.g008 An Electrostatic Funnel in the GABA-Binding Pathway 'random' sample points on the 2.5 nm shell were also analyzed (Fig 8E-8G). These points are positioned on the opposite side of the 2.5 nm shell to the midpoint (and are away from the protein), and are~1.0 nm from each other. The average number of GABA molecules within 1.0 nm of each 'random' sample point is 4.55 ± 1.8. Thus, there does not appear to be a significantly increased population of GABA molecules initially moving towards the midpoint region (6 vs 4.55 ± 1.8). The probability of a GABA molecule randomly reaching the 2.5 nm shell within a specific region of radius 1 nm is~4.3% (the percentage of available sampling space that is within 1 nm of a specific point, see Fig 8). Given that the average distribution we observe is~4.55%, the initial movement of the GABA molecules is indeed almost completely random/ non-biased in nature.
However, once the GABA molecules 'randomly' reach the vicinity of the midpoint, the influence of the electrostatic field takes over, and the molecules are funneled and focused into the binding pocket. This funnel acts as a sink on the population of randomly moving GABA molecules; even though the GABA molecules are unbiased and can reach any point on the 2.5 nm shell. If the GABA molecules reach a point on the 2.5 nm shell that is near the vicinity of the midpoint, the molecule becomes 'trapped' and is concentrated into the binding region. The time-dependent density distributions of the GABA molecules further illustrate this hypothesis (Fig 9 and S5 Fig). The molecules that arrive at the midpoint/binding site regions persist for substantial periods of time, and thus these areas become more densely occupied by GABA molecules (Fig 9B-9F), while the other GABA molecules randomly disperse from their starting locations. The GABA-midpoint distance calculations also support the idea of the focusing of the GABA molecules toward the binding site with 86% (24 of 28) of the simulations successfully proceeding to the binding region once they have reached the midpoint.

Why do we get three 'binding' outcomes?
All of the trajectories that fall into the three different simulation categories (BIND, PARTIAL, and NEARBY) have been thoroughly analyzed. We find that the GABA molecules in these different categories all follow similar pathways, all converge at the same point~1.9-2.0 nm from the GABA binding site COM, and importantly, all at least mostly enter the GABA binding region. However, the question remains, if the pathways are so similar, and they all reach this converged midpoint, then why do we observe these three distinct outcomes?
The orientation of the net dipole on the GABA molecules was decomposed into two components; the angle it makes with the XY plane, and the angle it makes with the YZ plane (Fig 10). When the molecules reach the midpoint at 1.9 nm from the binding site COM, the orientations of the dipole in both the BIND and PARTIAL categories are virtually identical (angles of~65°a nd~60°, respectively). However, the orientation of the dipole in the NEARBY category is nearly orthogonal to the BIND and PARTIAL dipoles (Fig 10C). Thus, the GABA molecules in the NEARBY category enter the electrostatic field with a different orientation.
One reason for this alternate orientation is the presence of a charged amino acid sidechain near the GABA molecule. The Arg207 residue from the β 3 subunit is close to the midpoint and may influence the GABA alignment. Visual assessment agrees with this hypothesis (Fig 11A,  inset). Furthermore, analysis of the number of contacts that Arg207 makes to GABA indicates that prior to the midpoint, GABA-Arg207 contacts are only present in the NEARBY simulations ( Fig 11A). A possible cause for these 'premature' GABA-Arg207 contacts is revealed in the distribution of arginine sidechain rotamers (Fig 11B). Of the nine possible rotamers formed by the N-CA-CB-CG and CB-CG-CD-NE dihedrals, only one allows consistent positioning to a location where GABA contacts are possible (Rotamer #3, Fig 11B, where the N-CA-CB-CG and CB-CG-CD-NE dihedrals are both between 240°and 360°). This rotamer accounts for a significant (~20%) proportion of Arg207 conformations as GABA enters the electrostatic field during the NEARBY simulations (Fig 11C), but is not present at the equivalent position in any of the BIND and PARTIAL simulations.
The increased population of this specific arginine rotamer in the NEARBY simulations may form additional contacts to GABA molecules, which subsequently influences their orientation as they enter the electrostatic field. Thus, the GABA molecules in the NEARBY simulations are in a different orientation and appear to be 'swept along' to a slightly different position, leaving them in a sub-optimal orientation to then 'flip up' into the binding pocket. Notably, most of the GABA molecules in this category actually persist near, or just inside the binding region for the duration of the simulations, and it may be the case that given extended simulation time they would indeed eventually reorient into the GABA binding site. Our statistical analysis does not reveal the cause-and-effect relationship between the Arg207 rotamer and the GABA orientation. We have hypothesized that the possible increased presence of Arg207 rotamer #3 induces the altered orientation of the GABA molecule. However, it is also plausible that the GABA molecule may have already adopted that orientation, and it is the occurrence of this orientation that causes the increase in the Arg207 rotamer #3 population.
Therefore, a differing orientation may account for the varying behavior seen in the NEARBY category of simulations. In contrast, the BIND and PARTIAL simulations both occupy the 'correct' orientation and as such, both proceed into the binding site. Once in the binding site, the GABA molecules from the BIND trajectories remain there for the remainder of the simulation (an average of almost 7 ns-70% of the simulation time), whereas molecules from the PARTIAL trajectories leave after an average of only~2.4 ns.
To suggest a potential basis for this differing behavior, components of the binding site were investigated. As the system was fully solvated and had an effective ionic concentration of 0.15 M, there are 135 Clions present in the simulation. Analysis of the Cldensity shows that during the PARTIAL simulations (where GABA leaves the binding pocket), there is a higher density of Clions in the core of the binding pocket compared to the BIND simulations. Early GABA A -R studies showed that higher concentrations of Clions modulate GABA binding [35,36], as well as inhibiting muscimol binding into the GABA site [37]. Moreover, anions in general have been reported to perturb GABA binding [38,39], suggesting that they may prevent/ shield key interactions between the carboxyl end of the GABA molecule and the basic binding residues.

Discussion
Through unbiased molecular dynamics simulations we have constructed a possible pathway by which GABA molecules enter the binding site of the GABA receptor. We also hypothesize an electrostatically driven mechanism for these GABA molecules to be 'detected'. The GABA A -R system represents an ideal test case to investigate ligand binding, as the binding site is well defined and the ligand interactions are primarily electrostatic, which have long-range effects.  ) and CB-CG-CD-NE (χ 1 dihedral) dihedrals of Arg207 were analyzed in terms of their distance to the average GABA position (for the bin when GABA is 2 nm from the binding site). These distances were calculated and standardized for each rotamer during all of the NEARBY simulations. This shows that the conformation of rotamer #3 allows Arg207 to get significantly closer to GABA than any other rotamer. This is illustrated by a visual comparison between rotamer #3 and rotamer #4 (B-inset, where the position of GABA is represented by an olive green sphere). (C) The percentage population of Arg207 rotamers for each simulation category was calculated as function of the GABA distance from the binding site. This shows that in the NEARBY simulations, Arg207 occupies rotamer #3 for a significant percentage of time as GABA is entering the electrostatic field (indicated by the highlighted red region), whereas in the BIND and PARTIAL simulations Arg207 does not occupy this rotameric conformation at all during the time GABA is entering the electrostatic field. This methodology may be applicable to other protein systems in terms of characterization of the ligand-binding pathway.
Our findings indicate that all the GABA molecules that at least partly enter the binding pocket follow a very similar pathway, whereby they approach the protein from below the binding site, before progressing up behind the C-loop and into the binding pocket. The pathway passes through a 'midpoint'. This midpoint is a region~1 nm below the binding site, and~1 nm laterally out away from the binding site, where the combined electrostatic potential surface of the protein creates a very strong electrostatic field.
The midpoint is a critical decision point, where the GABA molecule is 'captured' by the receptor and its far-reaching electrostatic field. The behavior and distribution of the GABA molecules up to this point is almost random/unbiased, as illustrated by their indiscriminate initial dispersion (Fig 8) and distribution density (S5 Fig). However, once the GABA molecules reach the midpoint they get 'swept along' by the electrostatic field, which funnels them to a position much nearer the binding site, whereupon more specific (but shorter distance) interactions can take over. The 'capturing' nature of the midpoint is highlighted by the fact that 86% (24 of 28) of the simulations where GABA reaches the midpoint then proceed on to the binding region. This is further illustrated by accumulating density of GABA molecules in this area (Fig  9 and S5 Fig).
Thus, the binding pathway may operate in a two-step, two-resolution manner. The longerrange electrostatic field interactions are used to 'pull in' molecules that match the general properties of the GABA A -R agonist-small and zwitterionic. Once these molecules have been brought closer to the protein, the precise sidechain-ligand contacts could determine if the molecule is suitable for binding. Previous studies suggest that the electrostatic interactions between glutamate and the glutamate receptor LBD become significant at~5 Å [40]. Thus, diffusing glutamates within about 5 Å of the protein are readily drawn in to the binding site through electrostatic interactions. We suggest that the effects on GABA may be significant at even longer distance.
This hypothesis of a non-specific electrostatic interaction is further highlighted by the fact that the electronegativity is a general property of the region, rather than a specific residue, and that the midpoint coordinates are actually~0.6 nm from the surface of the protein. Experimental investigation into the residues of this region may require more than single-point mutagenesis in order to alter the overall electronegative nature of the area. A fundamental difficulty in the assessment of mutation effects is that if ion flow through the channel is the endpoint measurement, then it can be problematic to distinguish a change in ligand binding versus a change in channel gating.

Model construction and validation
The model used in this work was derived from previously published studies, [21,32] which contain more detailed method descriptions. In brief, the main template of the GABA A -R model was the Torpedo marmorata nAChR structure (PDB [41] ID: 2BG9). [42] This template was used as a scaffold upon which other better-resolution, higher-homology protein sections were incorporated. This protocol was used for regions that were missing in the structure, lacked alignment to the template, or were poorly conserved. The DIG (Deletions Insertions Gaps) tool from the LGA (Local-Global Alignment) [43] program was used for this additional alignment and modeling. The DIG tool analyzes protein structures, or fragments of protein structures, and can complete the missing sections by searching for areas of similar sequence from a database of structural folds or a manually selected library of appropriate structural regions. In this study, as we were particularly interested in GABA binding, we concentrated on the LBD. Thus, we focused on gaps and regions that were poorly modeled/aligned in this domain. In order to fill these gaps and increase the accuracy of these regions, we searched homologous sections from all available pLGIC structures and the LBD-analogous AChBPs (such as PDB IDs: 2BYN [23] and 1UX2 [44]). Thus, the models of the LBD domains were completed and refined using these additional structural data, resulting in an overall model that has a better resolution and is more closely aligned to the GABA A -R sequence. In order to increase the likelihood of observing a GABA-binding event, the starting LBD structure was modeled using the available apo structures, rather than structures with a ligand-bound, 'closed' C-loop conformation, such as the glutamate-bound GluCl [34]. We also chose to model the α 6 β 3 δ receptor as it has the highest affinity [8] for GABA and thus the greatest chance of binding success. Consistent with published modeling studies of GABA A -R, [33,45] the cytoplasmic domain was not included due to high structural uncertainty. The integrity of these models was assessed using PROCHECK [46,47]. Most parameters were typical of a structure of 1.5-2.5 Å resolution, an improvement on the main template resolution, and an enhancement of the overall quality.
The protein model was inserted into a preformed and equilibrated POPC (palmitoyl-oleoyl phosphatidylcholine) bilayer that was used for previous GABA A -R simulations. [21,33] The system was solvated and had counter-ions added to neutralize the charge. Further Na+/Cl − ions were added to create an effective concentration of 0.15 M. The final system (consisting of~210,000 atoms) was energy minimized. Following minimization, the system was simulated for 2 ns with harmonic positional restraints on the protein, allowing the relaxation of lipid and water molecules. After allowing the packing of the lipids around the protein, the area per lipid was calculated and found to be representative of POPC. Subsequently, the area of the XY-plane was fixed (to maintain the area per lipid), positional restraints were removed, and the entire apo system was run for a further 20 ns to fully equilibrate the system. The root-mean-square deviation (RMSD) of the Cα atoms during this equilibration compared favorably to previous GABA A -R simulations [33]. DSSP (Define Secondary Structure of Proteins) [48] was used to measure the secondary-structure conformation of the GABA A -R. This confirmed that the structural composition of the GABA A -R model remained consistent throughout the 20 ns equilibration. Finally, there are critical GABA A -R salt bridges that are known to be important for gating dynamics [49] or maintaining the binding pocket structure [19,25]. Analysis showed that these salt bridges were present in the homology model and remained intact during equilibration.

Simulation system setup
Four separate, random frames were taken from the last 5 ns of the apo GABA A -R equilibration simulation. In each of these four frames, a GABA molecule was placed in a random orientation at a different location~2.50 nm from the center of the binding site. These four systems were used as 'seed points' and each underwent molecular dynamics simulations for one ns. From each of these four seed simulations, 25 frames (output at every 2 ps) were randomly chosen, producing 100 different starting positions for GABA in various conformations and orientations that are~2.46 ± 0.50 nm from the center of the binding site (Figs 2 and 12). Thus, unbiased, randomly oriented starting positions for GABA molecules were generated. Okada et al. [50] recently demonstrated that they could achieve binding of a ligand within~2 ns when it was placed~0.5 nm from the binding site. Thus, for our systems with a starting distance of~2.5 nm to the center of the binding site, a timescale of 10 ns was an appropriate choice to observe binding. Thus, each of the 100 systems was simulated for 10 ns, to produce 1 μs of total data for subsequent analysis.
All molecular dynamics (MD) simulations were run using CHARMM [51] (v27 for lipid and protein) in NAMD [52]. As GABA is a simple amino acid, we used the PARAtool [53] VMD [54] plugin to carry out force field parameterization by analogy to existing residues in the CHARMM force field to determine the bonded interactions. The partial charges are obtained using the RESP procedure with GAUSSIAN03 [55], at the Hartree-Fock 6-31G Ã [56][57][58] level of theory. GABA was modeled as a zwitterionic molecule. All systems use NPT and a nonbonded vdw cutoff of 1.2 nm. Constant pressure and temperature were maintained with Langevin pistons set at 1 atm and 310 K, using an oscillation time constant of 200 fs, a damping time constant of 50 fs, and a temperature damping coefficient of 1/ps. As the CHARMM force field parameters used do not reproduce the observed area per lipid over long MD trajectories, the area of the system in the XY-plane was kept constant. All simulations were run with particle mesh Ewald electrostatics, SHAKE [59], and TIP3P water [60], and with a 2 fs time step. All the simulations were run on the Sierra Linux cluster at the Livermore Computing Center (LC), consisting of 23,328 Intel Xeon EP X5660 cores.

Trajectory preparation
After the completion of the simulations, the trajectories were categorized depending on the distance of GABA relative to the binding site. As we were investigating the binding pathway, further preparation was needed for some of the simulations. In all of the PARTIAL and some of the NEARBY trajectories, the GABA molecules reach a position within or near the binding site COM, but subsequently leave and diffuse away. These 'leaving' sections of the simulation may obscure the analysis results and were not of interest. In those instances, the trajectories are split An Electrostatic Funnel in the GABA-Binding Pathway at the point where the GABA molecule reaches its minimum distance to the binding site COM. The second 'leaving' section of the trajectory was discarded, and only the first 'binding' section was used for analysis.
After refining the simulations by removing the 'leaving' sections, all the trajectories within a particular category were combined. In order to standardize the molecule positions and orientations, the frames (output at every 20 ps) within each group were aligned to the Cα atoms of the protein backbone by a least-squares fitting method. These frames were sorted based upon the distance between the GABA COM and the binding site COM. These sorted frames were binned into windows defined by 0.1 nm increments of this GABA COM to binding site COM distance. The overall average number of descriptors (such as position, standard deviation of position, dipole direction, and dipole strength) of the GABA molecules within these distance-defined bins were calculated for each bin and compared between the four categories; BIND, PARTIAL, NEARBY, and RANDOM.
While we have defined a GABA molecule as 'binding' when within a distance cutoff to the binding site COM, these compounds still may not occupy the correct binding orientation. Indeed, the objective of this study is to determine the route and mechanism by which the GABA molecules get into the binding pocket, but not the specific details of sidechain contacts or subtle molecule rearrangements once within the pocket. In fact, formation of these precise interactions may be beyond the timescale of our study. Despite this, however, we are confident that our 'binding' simulations do indeed represent the initial stages of true GABA binding. The BIND simulations reach a final position that almost overlaps the binding site COM (Fig 13A), and once they have reached that position, they persist there for the remainder of the simulation (only undergoing minor fluctuations in position with the binding site; as possible indication of GABA reorientation into an optimal binding conformation). Furthermore, in three of the nine BIND simulations, after GABA has reached the binding site, closure of the C-loop around the GABA molecule was observed ( Fig  13B). This process is indicative of the early stages of LGIC activation by a ligand and may represent GABA reaching the appropriate binding site position.

Analysis and visualization
The PDB2PQR [61,62] and Adaptive Poisson-Boltzmann Solver [63][64][65] software packages were used for the modeling of the biomolecular solvation through the solution of the Poisson-Boltzmann equation. VMD [52,54] was used to map the electrostatic potential to the biomolecular surface to produce an electrostatic potential surface and visualization of the electrostatic field lines. All further analysis was carried out using GROMACS [66] functions, VMD, and locally written scripts. Figure preparation was done using VMD. The electrostatic surfaces were calculated for various different GABA A -R homology models that were constructed using either different GABA A -R subunit sequences or different structural templates. They are (L-R): model of the alternative α 1 β 2 γ 2 isoform, α 6 β 3 δ model using a complete AChBP for the LBD template, α 6 β 3 δ model using GluCl as a template, α 6 β 3 δ model using the new GABA A -R β 3 homopentamer as a template.