Sodium-assisted formation of binding and traverse conformations of the substrate in a neurotransmitter sodium symporter model.

Therapeutics designed to increase synaptic neurotransmitter levels by inhibiting neurotransmitter sodium symporters (NSSs) classify a strategic approach to treat brain disorders such as depression or epilepsy, however, the critical elementary steps that couple downhill flux of sodium to uphill transport of neurotransmitter are not distinguished as yet. Here we present modelling of NSS member neuronal GAT1 with the substrate γ-aminobutyric acid (GABA), the major inhibitory neurotransmitter. GABA binding is simulated with the occluded conformation of GAT1 homodimer in an explicit lipid/water environment. Simulations performed in the 1-10 ns range of time elucidated persistent formation of halfextended minor and H-bridged major GABA conformations, referred to as binding and traverse conformations, respectively. The traverse GABA conformation was further stabilized by GAT1-bound Na(+)(1). We also observed Na(+)(1) translocation to GAT1-bound Cl(-) as well as the appearance of water molecules at GABA and GAT1-bound Na(+)(2), conjecturing causality. Scaling dynamics suggest that the traverse GABA conformation may be valid for developing substrate inhibitors with high efficacy. The potential for this finding is significant with impact not only in pharmacology but wherever understanding of the mechanism of neurotransmitter uptake is valuable.

INTRODUCTION

Cellular uptake of a large variety of essential amino acids, biogenic amines, related nutrients and osmolytes against their concentration gradients in cells is facilitated by Na+ ion symporters coupling entropy increase to uphill solute translocation on the expense of transmembrane electrochemical gradient [1]. Uptake of inhibitory neurotransmitter amino acids gamma-aminobutyric acid (GABA) and glycine as well as monoamine neurotransmitters serotonin, dopamine and norepinephrine is regulated by the SLC6 family subclass neurotransmitter sodium symporters (NSSs) that are also important targets for central nervous system therapeutics and illicit drugs.

A basic problem in understanding the mechanism of NSSs-mediated cellular uptake is how uphill neurotransmitter translocation is spatially and temporally linked to the downhill passage of Na+ ions. In this paper we shall therefore investigate events evoked by substrate binding and the role for Na+ ions, taking the major human GABA transporter subtype (hGAT1) as an example. Molecular modelling of hGAT1 subtype and simulation of the GABA/Na+ symport (co-transport) mechanisms has become possible by the high-resolution structure elucidation of the outward open, occluded [2] and inward open [3] conformations of the bacterial orthologue leucine symporter (LeuT) from Aquifex aeolicus. Although GABA transporter is functionally active as a monomer [4], proper trafficking requires its dimerisation providing GABA transporters as dimers or clusters of dimers at the cell surface [5]. In order to disclose functional consequences of dimerisation, a homodimeric model of the occluded hGAT1 was constructed using the LeuT template [2].

MATERIALS AND METHODS

Modeling hGAT1 Homodimer from LeuT Dimeric Crystal Structure

The model was built by the Modeller 9.10 program [6], using the alignment of Beuming et al., [7]. In detail, the dimeric structure of 2A65 was downloaded from the PDB database, and the residues Leu, beta octylglucoside and water were cut out from the structure, and the protein chain and the four sodium ions (two on each side) were retained. The sequence of hGAT1 (SC6A1_HUMAN) was downloaded from SwissProt (UniProtKB/Swiss-Prot code: P30531). This sequence in duplicate was continuously aligned to chain A and chain B of 2A65. The alignment was set manually according to Beuming et al., [7]. The position of the two Na+ ions Na+(1) and Na+(2) was set to be directly taken from the pdb template, to occupy the same position (coordinates) in the hGAT1 model as in 2A65. Putative disordered regions of hGATs were searched using the IUPRED server [8, 9], which revealed none or negligible disorder even in the flexible EL2 loop, therefore we found it reliable to include it in the automated modelling procedure. However, to restrict the number of possible conformations, we took advantage of the fact, that the EL2 loop for NSS transporters contains a disulfide bond formed between C164 and C173 [10], therefore we integrated Modeller’s “DISU patch” (disulfide modelling option) into our modelling routine and set these residues to form disulfide bonds in both chains. The resulting structure from Modeller was obtained in pdb format. Hydrogen atoms and MMFF94 charges were put on this model in SYBYL, where also acetyl and N-methyl groups were placed on the N and C termini. In addition, Cl- ions were added in SYBYL to the suggested binding site [11, 12], namely among residues Tyr86, Ser295, Asn327 and Ser331. To evaluate molecular motions that can be considered specific for cellular uptake of GABA, we also built an inward open hGAT1 model based on the inward open LeuT variant (pdb code: 3TT3, [3]). In contrast to the occluded state LeuT, however, the inward-facing variant was crystallised as a monomer. Therefore the dimeric inward open model was constructed in two steps. First, a dimeric template was built by fitting two 3TT3 monomers on the 2A65 dimer followed by modelling the open-to-in conformation of hGAT1 based on this dimeric template. Neither Na+ and Cl- ions, nor GABA were included during the model building.

GABA Docking and MD Simulations

GABA docking to primary site S1 [2, 13-17] in the occluded conformation of the homodimeric hGAT1 model was performed using the GOLD 4.0 docking program (Cambridge Crystallographic Data Centre Software Ltd., Cambridge, UK). Monomers A and B refer to the chain names in 2A65 dimer. The search space for docking was a 10 Å sphere around Na+(1) on each side. GABA docking was repeated 10 times and the average of the GOLD scores in 10 runs was used for further analysis. Docked GABA residues were merged to the protein in SYBYL and the structure obtained so far was subjected to molecular dynamics (MD) simulations using the CHARMM modelling package [18] in a 1-palmitoyl-2-oleoylphosphatidylcholine(POPC) membrane bilayer. In the MD simulations hydrogen atoms were automatically generated and the structures were energy minimized in a protein only environment in 100 steps steepest descent cycles. During the minimization procedure, positions of heavy atoms were kept in place with a harmonic restraining force (first 50 then 10 kcal/mol/Å). Scripts from the CHARMM-GUI [19] membrane builder module were used for the generation, minimization and preproduction equilibration calculations of the hGAT1 in POPC membrane bilayer system. The simulation system comprised 154963 (occluded) and 157649 (inward-facing) atoms consisting of hGAT1, GABA, POPC, TIP3 water molecules and sodium, potassium and chloride ions. During constant temperature, pressure and fixed lipid area MD simulations Langevin piston algorithm was used to maintain the temperature at 303 K and a pressure of 1 atm. The pressure was controlled by the Langevin barostat with a piston collision frequency of 20 ps-1 and the target temperature was kept with the thermal piston mass of 2000 kcal*mol-1*ps2. Electrostatic interactions were evaluated using the PME method with a grid spacing below 1 Å. The bonds between hydrogen and heavy atoms were constrained with the SHAKE algorithm. In all simulations the all-atom force field CHARMM22 for the protein and CHARMM36 for the lipid atoms [20] were used. The integration time-step was 2.0 fs. Figures were generated using the enhanced perspective option of Pymol.

MD simulation time was 13.641 ns. The hGAT1 backbone atoms RMSD plots during the last 5 ns of the simulation did not show any significant change. For further analysis, structures from the last 1 ns simulation were used. Two additional similar MD computations were also performed. In the second simulation (run 2) only the initial random number used at the initial step of the MD simulation was changed. In the third simulation (run 3) the values of the constraining force used during the minimization procedure were doubled compared to original.

RESULTS

Construction of the hGAT1 Homodimer Model

We modelled the major human GABA transporter subtype hGAT1. The occluded conformation of the homologous bacterial LeuT homodimeric structure (pdb code: 2A65) was used as a template for model building. Since LeuT formed dimers in the crystal in the occluded state, this structure can be directly used to model hGAT1 in the dimeric state (Fig. top). To evaluate molecular motions that can be considered specific for cellular uptake of GABA, we also built an inward open hGAT1 model based on the inward open LeuT variant (pdb code: 3TT3, (3) (Fig. bottom). The main difference between the inward-facing and the occluded states of the hGAT1 homodimers was the tilted position of the TM1 helices, similarly to the template LeuT structure [3].

A slightly asymmetrical structure, most pronounced at the EL2 loop regions (Fig. top) was obtained for the homodimer hGAT1, which is possibly the consequence of the lack of a complete loop in the LeuT template that allows a high degree of flexibility in this region. The transmembrane parts of the protein, exhibited only subtle differences in side chain organization. In contrast to our hGAT1 monomer model [13, 14], the OH side chain of Tyr60 points towards the S1 binding crevice in the homodimer. It is to note, however, that the use of different modelling program (Modeller [15, 16) instead of Swiss PDB Viewer [13, 14]), can be accounted for the difference.

Early Formation of H-Bridged GABA Assisted by hGAT1-Bound Na+(1)

After the construction of the hGAT1 homodimer model, GABA was docked to the S1 substrate binding crevice. The docking was made separately for the two monomer units (monomer-A and monomer-B) using the GOLD 4.0 docking program. Notably, GABA docked to the S1 substrate binding crevices with high docking scores [48-50]. GABA oxygen atoms turned towards Na+(1) forming a GABA(O)-Na+(1) complex (d = 3 ± 0.2 Å) in both monomers without the need to apply any other constraints. The highest scoring GABA molecules were then merged to the protein chains and these complexes were later used in molecular dynamics (MD) runs. Since Tyr60 points to the binding crevice, the docked GABA is able to contact its side chain, in line with the models of Wein and Wanner [15] and Skovstrup et al. [16] (see also the previous section). Following GABA docking procedures, the structures obtained so far were subjected to short MD simulations using the CHARMM modelling package in a POPC membrane bilayer [18-20].

We have observed that in MD runs 1 and 3, the GABA(O)-Na+(1) complex formed a ring-like structure with a characteristic 3 Å CD-N distance of GABA in both monomer-A and monomer-B. To test whether the formation of a ring-like GABA is energetically favourable, "in vacuo" Hartree-Fock and B3LYP density functional theory calculations (6-31G** and 6-311++G** basis sets) were performed by using commercial software Gaussian E01 (MO calculations were performed with the DFT method at the B3LYP/6-31G* level using the Gaussian 03 program package). It was found that GABA formed into a ring-like structure even when the calculations started from the linear GABA conformation. The ring-like GABA structure was further stabilized by the addition of Na+(1) ion. In MD run 2 that is basically identical to MD run 1 and 3, however, GABA(O)-Na+(1) complex showed different conformations in the two monomers, namely the ring-like structure in monomer-A and a half-extended GABA(O)-Na+(1) complex with d(CD-N) = 4.5 Å in monomer-B (Fig. ). Formation of the ring-like Na+(1)-GABA complex can be characterised by an apparent rate constant around 3 x 108 s-1 (Fig. ). These findings indicate that the occurrence of the ring-like GABA(O)-Na+(1) complex as well as the co-existence of half-extended binding [13] and ring-like [14, 16] traversing [21] conformations of GABA in the hGAT homodimer are statistically relevant events in the hGAT1 homodimer.

When the structures of hGAT1 monomer-A and monomer-B were compared, one of the MD simulations (run 1) showed a striking difference between the overall structures of the two monomers. Monomer-A remained in the initial occluded state conformation, while monomer-B showed large-scale molecular rearrangement including the opening of a water-permeable pore that finally arranged a structure resembling the inward-facing conformation (c.f. Fig. bottom). These MD data suggest that despite the sequence identity of monomer-A and monomer-B in the hGAT1 homodimer, asymmetries have evolved between the two monomers of the GAT1 homodimer, i.e. the two monomers can potentially behave differently.

To get insights into the molecular motions underlying functional transitions, we analyzed the movements of individual residues. We calculated the RMS difference of each residue between monomer-A and monomer-B in MD run 1. Besides the extra- and intracellular loops, the highest ΔRMS values were observed for residues in the transmembrane helices near the GABA binding crevice (TM1, TM6, TM7 and TM8), but interestingly, high residue mobility was found for TM5 and TM11 that are relatively far from the substrate binding site (Fig. ). According to their flexible nature, the extra- and intracellular loops, especially the highly flexible EL2 loop showed prominent movements. The rearrangements, however, also involved the transmembrane helices close to the S1 GABA binding site (TM1, TM5, TM6, TM7, TM8 and TM11, Fig. ).

Detailed analysis of the structural changes in monomer-B revealed functionally essential molecular transitions. GABA in monomer-B lost its coordination to Na+(1), instead, it participated in H-bonding with two water molecules that were penetrated into the binding cavity during MD (Fig. ). Concurrently, Na+(1) formed salt bridge with bound Cl- ion. Na+(1) ion was coordinated by residues from TM1 (Ala61, Asn66), TM6 (Ser295) at start, and gained another coordinating residue from TM7 (Asn327) after 13 ns. Most importantly, the salt bridge formation with bound Cl- ion pulled Na+(1) away from GABA concurrently. Na+(2) was coordinated by less residues after 13 ns MD collecting a water molecule into its binding crevice (Fig. 5). These energetically favourable changes shape monomer-B in a conformation clearly different from monomer-A, also characterised by the unbinding tendency of Na+(2) in monomer-B (Fig. ). In monomer-A, however, GABA remained in complex with Na+(1) ion, i.e. Cl- ion did not interact with Na+(1) ion after MD simulation. In addition, neither GABA nor Na+(2) ion became hydrated (data not shown). These binding interactions are similar to those observed in the occluded conformation [2], supporting that monomer-A mimics the occluded state.

DISCUSSION

Temporal analysis of global residue movements throughout the MD simulation has shown that the molecular rearrangement of the transporter was initiated at the extracellular surface mostly involving helices implicated in GABA, Na+(1) and Cl- binding TM1, TM6, TM7 and TM11 [21] (Fig. ). Eventually, these movements were conducted through the transmembrane helices to the intracellular side. Importantly, after approximately 4-6 ns, the TM8 became involved in the molecular restructuring also initiating its movement at the extracellular surface. The involvement of TM8 might be of paramount importance, because its residues are associated to GABA and Na+(2) binding and also to the formation of intracellular gate [21] and the water penetration pathway, through which water accesses GABA and Na+ ions to get them hydrated (Fig. ). Local events occurring in binding crevices of the homodimer disclosed that the two monomers may undergo different fluctuations, resulting in distinguishable conformations allowing traverse of GABA in one of the monomers (Fig. ). GABA prefers the H-bridged ring-like (traversing) conformation in complex with Na+(1), although either half-extended or Na+(1) released conformations possibly co-exist in the two monomers. H-bridged ring-like GABA coordinates Na+(1) in monomer-A, while in monomer-B, it is hydrated by water molecules (Fig. ).

In the same run, Na+ ions are found in a more polar environment bridging Cl- ion (Na+(1)) and water (Na+(2)) (Fig. ). The asymmetric fluctuations observed in one of the MD simulations of the occluded state hGAT1 homodimer model, therefore, may provide useful information to better understand the mechanistic clues underlying substrate traverse through NSSs also featuring the LeuT-fold symmetry and a break in TM1 and TM6 helices forming Na+/solute binding crevice, like hGAT1 does.

We found structural fluctuations distinguishing monomers by short-term MD simulations of the occluded hGAT1 homodimer. Notably, water was able to penetrate into the binding crevice only in monomer-B, but not in monomer-A. In monomer-B water appeared next to GABA and Na+(2), like in the inward-facing homodimer model (Fig. bottom). Being persistent in the nanosecond range of time, statistical time events in the membrane symporter such as helix-coil transitions collecting solute binding, intramolecular hydrogen bonding, hydration and charge translocation [22] may spread throughout the monomer-B towards the intracellular membrane face similarly to the opening of the acetylcholine receptor channel [23]. In this way, monomer-B is ready to transit to the inward-facing state, therefore the structural fluctuations may be considered as functional, representing an intermediate state during the solute translocation process. Thus, the initial LeuT-fold asymmetry of the occluded state monomer, furthering in the homodimer after short MD may be indicative of forthcoming inward substrate release by hGAT1. We conjecture substrate/ligand-induced asymmetry of membrane protein oligomers including NSSs (this work) or ligand-gated ion channel subtypes [23-27] as a structural prerequisite for facilitating symport, channel opening or desensitisation.

In support of the notion that the lifetime of the traversing conformation that was modelled is in the order of nanoseconds we refer to data collected in Table . Plotting of log(kapp) vs. log(sampling time) reveals linear correlation (R2=0.95) indicating that the apparent rate constant obtained from macroscopic GABA uptake data scales right down to the microscopic elementary steps, such as the formation of traversing conformation of GABA described in the 1-10 ns range of time. It is to note, that opening of ionic channels also acts upon scaling dynamics [33].

CONCLUSION

It appears likely that time events appearing within the simulation time chosen in the homodimer of neuronal hGAT1 subtype reflect mechanistically relevant molecular events underlying complex kinetics of GABA-sodium symport [29 and references cited, 30-32, 34]. We extend our findings to astroglial hGAT2 and hGAT3 subtypes in particular and NSSs homodimers in general. Thus, deeper insights into binding/transport phenomena of NSS oligomers through coarse-grain simulation of the whole alternating access transport cycle in the future may open up new avenues in NSS-based CNS drug design and discovery [1, 21, 35, 36 and reference cited].

CONFLICT OF INTEREST

The authors confirm that this article content has no conflict of interest.

Table 1

Apparent rate constant for GABA uptake correlates the sampling time of techniques applied.

Sampling Time (s)	kapp (s-1)	References
10-11	3 x 1010	this work
2x 10-10	107	28
2x 10-10	7x106	17
4x10-5	3x102	29
8x10-5	70	30
10-4	10	31
10	5.4x10-2	32
100	4.44 10-3	32