Advertisement

Occlusion of the human serotonin transporter is mediated by serotonin-induced conformational changes in the bundle domain

Open AccessPublished:January 20, 2022DOI:https://doi.org/10.1016/j.jbc.2022.101613
      The human serotonin transporter (hSERT) terminates neurotransmission by removing serotonin (5HT) from the synaptic cleft, an essential process for proper functioning of serotonergic neurons. Structures of the hSERT have revealed its molecular architecture in four conformations, including the outward-open and occluded states, and show the transporter’s engagement with co-transported ions and the binding mode of inhibitors. In this study, we investigated the molecular mechanism by which the hSERT occludes and sequesters the substrate 5HT. This first step of substrate uptake into cells is a structural change consisting of the transition from the outward-open to the occluded state. Inhibitors such as the antidepressants citalopram, fluoxetine, and sertraline inhibit this step of the transport cycle. Using molecular dynamics simulations, we reached a fully occluded state, in which the transporter-bound 5HT becomes fully shielded from both sides of the membrane by two closed hydrophobic gates. Analysis of 5HT-triggered occlusion showed that bound 5HT serves as an essential trigger for transporter occlusion. Moreover, simulations revealed a complex sequence of steps and showed that movements of bundle domain helices are only partially correlated. 5HT-triggered occlusion is initially dominated by movements of transmembrane helix 1b, while in the final step, only transmembrane helix 6a moves and relaxes an intermediate change in its secondary structure.

      Keywords

      Abbreviations:

      5HT (serotonin), CG (coarse-grained), hSERT (human serotonin transporter), MD (molecular dynamics), SLC6 (solute carrier 6), TM (transmembrane helix)
      The human serotonin transporter (hSERT) retrieves serotonin (5HT) from the synaptic cleft and transports it into the presynaptic neuron, thereby spatially confining the serotonergic signal (
      • Fuller R.W.
      • Wong D.T.
      Serotonin uptake and serotonin uptake inhibition.
      ,
      • Murphy D.L.
      Serotonin transporter: Gene, genetic disorders, and pharmacogenetics.
      ,
      • Chen N.-H.
      • Reith M.E.A.
      • Quick M.W.
      Synaptic uptake and beyond: The sodium- and chloride-dependent neurotransmitter transporter family SLC6.
      ). The monoamine transporter hSERT is a secondary active transporter, which belongs to the solute carrier 6 (SLC6) family that utilizes the transmembrane concentration gradient of sodium to energize substrate transport. hSERT has been associated with a wide spectrum of neuropsychiatric disorders (
      • Bhat S.
      • El-Kasaby A.
      • Freissmuth M.
      • Sucic S.
      Functional and biochemical consequences of disease variants in neurotransmitter transporters: A special emphasis on folding and trafficking deficits.
      ,
      • Kristensen A.S.
      • Andersen J.
      • Jørgensen T.N.
      • Sørensen L.
      • Eriksen J.
      • Loland C.J.
      • Strømgaard K.
      • Gether U.
      SLC6 neurotransmitter transporters: Structure, function, and regulation.
      ). Among them, depression is listed by the World Health Organization as one of the world’s leading disabilities, affecting over 260 million people (
      • James S.L.
      • Abate D.
      • Abate K.H.
      • Abay S.M.
      • Abbafati C.
      • Abbasi N.
      • Abbastabar H.
      • Abd-Allah F.
      • Abdela J.
      • Abdelalim A.
      • Abdollahpour I.
      • Abdulkader R.S.
      • Abebe Z.
      • Abera S.F.
      • Abil O.Z.
      • et al.
      Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: A systematic analysis for the Global Burden of Disease Study 2017.
      ). Accordingly, hSERT is the target of several drugs that act as inhibitors of 5HT uptake, including clinically relevant drugs such as the selective serotonin reuptake inhibitors and the nonselective tricyclic antidepressants and also drugs of abuse like cocaine. The second group of compounds interfering with normal hSERT function are releasers such as amphetamines and substituted amphetamines, including fenfluramine or methylenedioxyphenethylamine (‘ecstasy'; (
      • Sitte H.H.
      • Freissmuth M.
      Amphetamines, new psychoactive drugs and the monoamine transporter cycle.
      )), which trigger the efflux of 5HT through hSERT. Recently, allosteric inhibitors (
      • Niello M.
      • Gradisch R.
      • Loland C.J.
      • Stockner T.
      • Sitte H.H.
      Allosteric modulation of neurotransmitter transporters as a therapeutic strategy.
      ) and partial releasers (
      • Rothman R.B.
      • Partilla J.S.
      • Baumann M.H.
      • Lightfoot-Siordia C.
      • Blough B.E.
      Studies of the biogenic amine transporters. 14. Identification of low-efficacy “partial” substrates for the biogenic amine transporters.
      ,
      • Hasenhuetl P.S.
      • Bhat S.
      • Freissmuth M.
      • Sandtner W.
      Functional selectivity and partial efficacy at the monoamine transporters: A unified model of allosteric modulation and amphetamine-induced substrate release.
      ) have also been described, highlighting the increasing complexity of hSERT pharmacology.
      The first insights into the structural organization of the SLC6 family were revealed by crystal structures of the small leucine transporter LeuT from Aquifex aeolicus and subsequently confirmed by hSERT structures (
      • Yamashita A.
      • Singh S.K.
      • Kawate T.
      • Jin Y.
      • Gouaux E.
      Crystal structure of a bacterial homologue of Na+/Cl--dependent neurotransmitter transporters.
      ,
      • Krishnamurthy H.
      • Gouaux E.
      X-ray structures of LeuT in substrate-free outward-open and apo inward-open states.
      ,
      • Gotfryd K.
      • Boesen T.
      • Mortensen J.S.
      • Khelashvili G.
      • Quick M.
      • Terry D.S.
      • Missel J.W.
      • LeVine M.V.
      • Gourdon P.
      • Blanchard S.C.
      • Javitch J.A.
      • Weinstein H.
      • Loland C.J.
      • Nissen P.
      • Gether U.
      X-ray structure of LeuT in an inward-facing occluded conformation reveals mechanism of substrate release.
      ,
      • Focke P.J.
      • Wang X.
      • Larsson H.P.
      Neurotransmitter transporters: Structure meets function.
      ,
      • Penmatsa A.
      • Gouaux E.
      How LeuT shapes our understanding of the mechanisms of sodium-coupled neurotransmitter transporters.
      ,
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      ,
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ,
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ) (Fig. 1).
      Figure thumbnail gr1
      Figure 1Sequestering of substrate by hSERT occlusion as part of the transport cycle. A, vertical dissection of the three experimentally solved conformations of hSERT (outward-open state (PDB ID: 5I71; (
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      )), the inhibitor bound outward-occluded conformation (PDB ID: 6DZV, (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      )) and the inward-open conformation (PDB ID: 6DZZ, (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ))). In addition, simulations revealed a new fully occluded structure of SERT, in which 5HT is completely sealed from the extracellular space through a complete closure of the extracellular hydrophobic gate. B, zoom on the substrate binding site S1, highlighting 5HT (green) and the side chains of the gate (salt bridge and hydrophobic gate) of the outward-open state and C, the fully occluded state. 5HT, serotonin; hSERT, human serotonin transporter.
      Conformational changes of substrate transport are associated with the pseudo-2-fold symmetry of the inverted repeat fold (
      • Forrest L.R.
      • Zhang Y.-W.
      • Jacobs M.T.
      • Gesmonde J.
      • Xie L.
      • Honig B.H.
      • Rudnick G.
      Mechanism for alternating access in neurotransmitter transporters.
      ). Substrate transport is enabled by motions of the bundle domain (transmembrane helix [TM]1, TM2, TM6, and TM7) relative to the scaffold domain (TM3, TM4, TM8, and TM9) (
      • Kristensen A.S.
      • Andersen J.
      • Jørgensen T.N.
      • Sørensen L.
      • Eriksen J.
      • Loland C.J.
      • Strømgaard K.
      • Gether U.
      SLC6 neurotransmitter transporters: Structure, function, and regulation.
      ), which allow for alternating the access to the highly conserved substrate binding site (S1), located in the center of the transporter (
      • Yamashita A.
      • Singh S.K.
      • Kawate T.
      • Jin Y.
      • Gouaux E.
      Crystal structure of a bacterial homologue of Na+/Cl--dependent neurotransmitter transporters.
      ). The extracellular and intracellular gates delimit the S1 and determine its exposure to the respective sides of the membrane. These experimental structures represent static snapshots for a highly dynamic transporter (
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      ,
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ,
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ). Complete occlusion of hSERT and sequestering of substrate in the S1 is important for the transport cycle. The occluded state is characterized by two closed gates, which sequester the transporter-bound substrate from both sites of the membrane. This state is very similar to the occluded conformation of 5HT-bound SERT (
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ), but differs from the cyrogenic electron microscopy (cryo-EM) structure of the inhibitor-bound outward-occluded structure of hSERT (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ). The substrate binding site S1 of the inhibitor-bound hSERT remains accessible by water from the extracellular side of the membrane by a narrow path, indicating that the inhibitor-bound outward-occluded cryo-EM structure of hSERT has not reached the fully occluded state.
      Conformational changes of the transport cycle of SLC6 transporters have been previously investigated using molecular dynamics (MD) simulations (
      • Cheng M.H.
      • Bahar I.
      Molecular mechanism of dopamine transport by human dopamine transporter.
      ,
      • Koldsø H.
      • Noer P.
      • Grouleff J.
      • Autzen H.E.
      • Sinning S.
      • Schiøtt B.
      Unbiased simulations reveal the inward-facing conformation of the human serotonin transporter and Na+ ion release.
      ,
      • Koldsø H.
      • Autzen H.E.
      • Grouleff J.
      • Schiøtt B.
      Ligand induced conformational changes of the human serotonin transporter revealed by molecular dynamics simulations.
      ,
      • Hellsberg E.
      • Ecker G.F.
      • Stary-Weinzinger A.
      • Forrest L.R.
      A structural model of the human serotonin transporter in an outward-occluded state.
      ,
      • Sinning S.
      • Musgaard M.
      • Jensen M.
      • Severinsen K.
      • Celik L.
      • Koldsø H.
      • Meyer T.
      • Bols M.
      • Jensen H.H.
      • Schiøtt B.
      • Wiborg O.
      Binding and orientation of tricyclic antidepressants within the central substrate site of the human serotonin transporter.
      ,
      • Thomas J.R.
      • Gedeon P.C.
      • Grant B.J.
      • Madura J.D.
      LeuT conformational sampling utilizing accelerated molecular dynamics and principal component analysis.
      ), typically by applying biasing potential to induce structural changes. These approaches depend on the definition of predefined reaction coordinates or collective variables, which are initially unknown and might therefore deviate from the lowest energy path. We used extensive bias-free equilibrium all atom molecular dynamic simulations (40 μs total simulation time) to investigate the substrate-induced occlusion of hSERT, thereby avoiding the danger of inducing structural changes which could deviate from the lowest energy path. Simulations of sodium- and chloride-bound hSERT show that ion binding stabilizes the outward-open state (
      • Tavoulari S.
      • Margheritis E.
      • Nagarajan A.
      • DeWitt D.C.
      • Zhang Y.-W.
      • Rosado E.
      • Ravera S.
      • Rhoades E.
      • Forrest L.R.
      • Rudnick G.
      Two Na+ sites control conformational change in a neurotransmitter transporter homolog.
      ,
      • Claxton D.P.
      • Quick M.
      • Shi L.
      • de Carvalho F.D.
      • Weinstein H.
      • Javitch J.A.
      • Mchaourab H.S.
      Ion/substrate-dependent conformational dynamics of a bacterial homolog of neurotransmitter:sodium symporters.
      ,
      • Szöllősi D.
      • Stockner T.
      Investigating the mechanism of sodium binding to SERT using direct simulations.
      ), while binding of 5HT to the S1 changes the conformational equilibrium leading to hSERT occlusion. We could identify the sequence of events essential for hSERT occlusion and discovered that TM1b and TM6a move remarkably independently of each other. The initial triggers of hSERT occlusion are electrostatic interactions between the positively charged amino group of 5HT and D98 and the negative helical dipole of TM6a. After initial movements of the bundle domain, interactions between the bundle and scaffold domains increase through close contacts of the side chains of the extracellular hydrophobic gate (I172, Y176, and F335). The subsequent de-wetting of the S1 and the outer vestibule leads to full 5HT shielding from the extracellular side through further movements of TM6a.

      Results

      5HT bound to the substrate binding site S1 induces occlusion

      We studied the 5HT-induced structural changes of hSERT that lead to substrate sequestering using extensive unbiased MD simulations (40 μs total simulation time) of ion-bound hSERT in the absence of 5HT (ten independent simulations) and in complex with 5HT (ten independent simulations) (Fig. 2).
      Figure thumbnail gr2
      Figure 2The outer gate closes in the presence of 5HT. A, graphical legend showing hSERT and the measured distances. The distances between the center of mass of TM1b (residue L99 to Q111) (red), TM6a (residue G324 to L337) (mauve), and the upper part of TM9 (TM9up) (residue F475 to S477) are indicated by yellow lines. Time evolution of distances between (B) TM1b-TM9up and (C) TM6a-TM9up of ten independent simulations of 5HT-free hSERT. Time evolution of distances between (D) TM1b-TM9up and (E) TM6a-TM9up of ten independent simulations of 5HT-bound hSERT. Simulations are colored according to the end state: in green, if a fully occluded state is reached; in purple, if an intermediate state is reached; and in orange, if the hSERT remains in the outward-open state. The horizontal dashed black lines indicate the distances as measured in three known structures of hSERT. 5HT, serotonin; hSERT, human serotonin transporter; TM, transmembrane helix.
      Figure thumbnail gr3
      Figure 3Both gates are completely closed in the fully occluded state. A and B, side view and top view of an overlay of the fully occluded state of replica 10 with the experimental structures of the outward-open (
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      ), outward-occluded (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ), occluded (
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ), and the inward-open (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ) state. For clarity, only TM9, TM1b, and TM6a are shown. All four structures are fitted on the scaffold domain [TM3 (160–185), TM4 (257–269), TM8 (424–449), TM9 (466–477)]. TM, transmembrane helix.
      To quantify vestibule closure, we measured the distances from the bundle helices TM1b and TM6a to the upper part of the scaffold helix TM9 (TM9up). TM9up was selected because it is located across the outer vestibule (Fig. 2A), not directly interacting with elements of the outer vestibule and not changing conformation during the transport cycle. Simulations of ion-bound hSERT in the absence of 5HT (Fig. 2, B and C) showed that distances from TM1b and TM6a to TM9up remain comparable to the distances observed in the outward-open crystal structure. Interestingly, the amplitude of motions and the variance between simulations was larger for the TM1b-TM9up distance as compared to the TM6a-TM9up distance, revealing a larger conformational freedom for TM1b.
      hSERT in complex with 5HT (Fig. 2, D and E) reached the fully occluded state in three out of ten 2-μs-long simulations (green-colored trajectories). The distances from the scaffold domain helix TM9up to the bundle domain helices TM1b and TM6a decreased by 0.32 ± 0.04 nm and 0.34 ± 0.01 nm, respectively, compared to the starting outward-open conformation. The transition led to sequestering of 5HT in the S1 by occlusion of hSERT through a closure of the outer vestibule. In a second group of simulations, hSERT reached an intermediate state (Fig. 2, D and E, purple trajectories), which was partially occluded. These simulations were characterized by an almost complete transition of TM1b, while TM6a moved only halfway toward the fully occluded conformation. This intermediate state was less stable because TM6a continued transitioning between the outward-open and the intermediate conformation in these three simulations. A third group of hSERT simulations (four simulations) in complex with 5HT (Fig. 2, D and E, orange trajectories) remained in the outward-open state. Interestingly, a comparison with the simulations in absence of 5HT showed that TM1b moves by ∼0.1 nm toward the bundle domain in the presence of 5HT, indicating that TM1b senses the presence of 5HT in the S1.
      A comparison between the fully occluded conformation of replica 10 and the experimentally determined outward-open, outward-occluded, occluded, and the inward-open structures (Fig. 3) showed that the outer half of the bundle domain (TM1b and TM6a), starting from the outward-open conformation, reached a conformation that is similar to the conformations observed in the occluded 5HT-bound and in the inward-open cryo-EM structures. In contrast to the transition of the outer vestibule half of the bundle domain, the inner vestibule half of the bundle domain did not move, remaining closed. These results indicate a full occlusion of hSERT that sequesters 5HT in the S1 by closure of both gates.

      Sequential independent movement of TM1b and TM6a

      The comparison of the time evolution of distances between TM1b and TM9up as well as between TM6a and TM9up (Fig. 2, D and E) showed that their motions are not fully coordinated. To unmask these de-synchronized movements, we analyzed the degree of correlation between the movements of TM1b and TM6a by plotting the distances from TM9up to TM1b versus TM9up to TM6a (Fig. 4A).
      Figure thumbnail gr4
      Figure 4Nonsynchronized movements of TM1b and TM6a lead to bundle domain rotations. A, 2D plot of all replicas of hSERT in complex with 5HT, showing all visited conformations (outward-open, intermediate, and fully-occluded) as cumulative counts, visualized by a grayscale map. The paths taken by the individual trajectories (averaged over 150 ns) are highlighted as colored lines, applying the color code as defined in . BD, representative trajectories for an occluding (replica 10) and a non-occluding (replica 6) trajectory, colored in teal and orange, respectively. The thicker line represents running averages of 20 ns, while the thinner lines show all data points with a time resolution of 1 ns. Time evolution of the distance between (B) TM1b and TM9up and (C) TM6a and TM9up. D, time evolution of the angle between TM9up-TM1b-TM6a. E, graphical legend visualizing the angle between TM9up-TM1b-TM6a as measured in panel D. 5HT, serotonin; hSERT, human serotonin transporter; TM, transmembrane helix.
      TM1b moved by 0.1 nm toward TM9up in all simulations of hSERT in complex with 5HT, while TM6a did not show a comparable movement (Fig. 2). This conformation of hSERT stayed outward-open, and the S1 remained fully accessible through the open outer vestibule. hSERT reached the intermediate conformations by synchronized movements of TM1b and TM6a, visible as the diagonal transition from the outward-open to the intermediate conformation (Fig. 4A). Correlated movements of TM1b and TM6a are characteristic for this transition, as all six simulations reaching the intermediate state follow the same path. The final step leading to the fully occluded state consists of movements of TM6a only, which is visible by the vertical transition in Figure 4A.
      These partially nonsynchronized motions of TM1b and TM6a (Fig. 4, B and C), through which the fully occluded state is reached (Fig. 4A), lead to a rotation of the bundle domain relative to the scaffold domain (Fig. 4D), as visualized by the angle between TM9up, TM1b, and TM6a (Fig. 4E). The initial isolated motion of TM1b in response to the presence of 5HT leads to a rotation of TM1b and TM6a by ∼10°, which relaxes to the original value of 80° in the last step, at the point of complete occlusion and closure of the outer vestibule. The bundle domain is therefore contorted in the intermediate state, suggesting that the bundle domain is in a strained, unstable conformation.
      Sequestering of S1-bound 5HT by hSERT occlusion is a global conformational change that represents the first step of the transport cycle. We used principal component analysis to investigate the transitions occurring in the simulations, collectively analyzing the entire dataset of 20 trajectories. Structural analysis of the first eigenvector (principal component 1 [PC1]) showed that PC1 represents (Fig. 5A) the main motion of the bundle domain leading to hSERT occlusion. The second eigenvector (principal component 2 [PC2]) represents a rotation of the bundle domain relative to the scaffold domain with the rotation axis oriented approximately perpendicular to the rotation described by PC1 (Fig. 5B). PC2 describes structural adjustments of the bundle domain relative to the scaffold domain that help the bundle domain to trap the bound 5HT.
      Figure thumbnail gr5
      Figure 5Collective motions of hSERT in response to 5HT binding. Principal component analyses of the 20 trajectories were used to identify conformational changes. A, the largest motion (principal component 1; PC1) is a rotation of TM1b, TM2, and TM6a of the bundle domain. B, PC2 is a rotation of the bundle domain helices, mainly TM1b, TM6a, and TM7, perpendicular to PC1. Projection of all trajectories (color code as in ) onto PC1 and PC2 for (C) the ten 5HT-free simulations, the 5HT-bound hSERT simulations (D) that remain open, (E) that reach the intermediate state, or (F) that fully occlude. 5HT, serotonin; hSERT, human serotonin transporter; TM, transmembrane helix.
      Motions of the individual simulations were projected onto these two largest amplitude eigenvectors (Fig. 5, CF). The four trajectories of 5HT-bound hSERT that remain in the open state (Fig. 5D) are characterized by a large amplitude of motion along PC2, while not moving along PC1. Reaching the intermediate state (Fig. 5E) included motions along PC1 and a restraining of motions in PC2. The trajectories that reached full occlusion (Fig. 5F) are characterized by full amplitude motion along PC1 and transition through the conformation of the intermediate state. Interestingly, movements along PC2 to reach the intermediate state are reversed upon reaching the fully occluded state, which is a property similar to the rotational motion of TM6a described in Figure 4D. In the absence of 5HT, hSERT showed uncorrelated dynamics in PC1 and PC2 (Fig. 5C), which indicates that 5HT binding restrains hSERT motion both along PC1 and PC2.

      Forces between 5HT and residues in the S1 initiate closure of the extracellular vestibule

      Movements of TM1b and TM6a in the presence of 5HT result in the closure of the outer vestibule. Minima on the potential energy surface correspond to stable or metastable states. Forces are derivatives of the potential energy and therefore reflect the steepness of the slope of the energy surface. These forces lead to conformational changes toward minima on the energy surface. To gain an understanding of the driving forces and the interactions that dominate the substrate-triggered global conformational changes, we analyzed the forces acting between 5HT and residues within the S1. Large negative forces represent strong attractions between a residue and 5HT, while repulsive forces originate from atoms that are too close to one another, as when 5HT is pushed toward a residue. An alternation between negative and positive values for forces (e.g., residue 439 in Fig. 6) indicates that 5HT oscillates around the equilibrium distance at the energy minimum.
      Figure thumbnail gr6
      Figure 6Electrostatic and van der Waals forces between 5HT and the residues of the S1. A, visualization of bound 5HT (green) and the surrounding residues which show the largest interaction forces. Residues are colored according to the average forces (electrostatic plus van der Waals) of interaction with 5HT. B, time evolution of interaction forces (electrostatic plus van der Waals) between 5HT and these residues are shown for the occluding replica 10. The horizontal line indicates the time point when the system transitions to the intermediate state. 5HT, serotonin; S1, substrate binding site; TM, transmembrane helix.
      Figure 6A visualizes the average forces from nonbonded interactions (electrostatic and van der Waals forces), mapped in color onto the respective residues, while Figure 6B shows as a heat map the time evolution of nonbonded forces for the fully occluding replica 10.
      Analysis of forces of the other systems are shown in the Fig. S1. Figure 6B shows that the forces which attract some of the crucial residues were weaker in the outward-open state than in the intermediate conformation. The pattern of the time evolution of forces separates the residues into two groups: (i) interactions of the residues that interact with the indole ring of 5HT changed minimally during the transition from the outward-open to the outward-occluded state (Fig. 6B) and (ii) in contrast, residues that mainly interact with the amino group of 5HT showed weaker interactions in the outward-open state, and these interactions became stronger upon reaching the intermediate state. Importantly, the interactions between 5HT and D98 had already increased before the transition.
      The first group of residues are located on the scaffold domain (TM3 and TM8) with Y176 (TM3), N177 (TM3), and T439 (TM8) being of particular importance. Y176, which is part of the outer hydrophobic gate, formed extensive van der Waals interactions with 5HT. The total potential energy of interactions between Y176 and 5HT was negative (−12.74 ± 4.70 kJ/mol), while the total force was positive, indicating that the distance was on average shorter than the equilibrium distance. Residues N177 and T439 formed strong hydrogen bonds with the hydroxyl group of the indole moiety of 5HT in the otherwise largely hydrophobic part of S1, thereby stabilizing 5HT.
      Residues of the second group belong to the bundle domain and the hydrophilic part of the S1. The most important interactions originated from interactions between 5HT and D98 (TM1) as well as F335 (TM6). The negatively charged D98 established a fluctuating salt bridge with the amino group of 5HT, which is indicative of a dynamic interaction and of changes in local geometry. The backbone carbonyl of residue F335 begun interacting with 5HT only during the transition to the intermediate state and was the main interaction partner of 5HT with the helix dipole of TM6a.

      Changes in interaction energy between 5HT and hSERT correlate with the first step of outer vestibule closure

      Next, we investigated the potential energies resulting from key interactions of 5HT with hSERT (Fig. 7) to correlate the time evolution of conformational changes that are associated with the closure of the outer vestibule. A comparison of the position of 5HT in the outward-open and the intermediate conformation showed that 5HT stayed firmly bound to TM3 and TM8, while the bundle domain moved closer (Fig. 7, AC) during occlusion. The analysis of the RMSD of the position of 5HT (Fig. 7E) showed that 5HT slightly changed orientation before TM6a moved toward the scaffold domain and returned to its original position after the transition to the intermediate conformation. The position of 5HT was stabilized by strong hydrogen bonds of the hydroxyl group of 5HT with N177 (TM3) and T439 (TM8) (Fig. 7, F and G).
      Figure thumbnail gr7
      Figure 7Essential interactions between 5HT and the S1 of hSERT. Top view of the central substrate binding site S1 highlights the most important interactions. A, in outward-open conformation; B, in the intermediate conformation. C, side view of intermediate conformation highlighting the key interactions of 5HT with hSERT. Panel DI, shows results from replica 10 and 6 as representative simulations for an occluding and a non-occluding trajectory. D, time evolution of distances between TM6a and TM9up. E, time evolution of the RMSD of 5HT. Time evolution of the potential energy of interactions between functional moieties of 5HT and (F) N177, (G) T439, (H) D98, and (I) F335. Thick lines represent running averages of 20 ns, while the faded lines and circles show all data points with a time resolution of 1 ns. 5HT, serotonin; hSERT, human serotonin transporter; TM, transmembrane helix.
      Interactions of 5HT with D98 (Figs. 6B and 7I) increased early in the simulations and became stronger with the transition of TM1b toward the scaffold domain. The interactions between 5HT and D98 remained dynamic, and the potential energy of interactions fluctuated. The conformation and dynamics of the D98 side chain depended on the presence of 5HT and the degree of occlusion (Fig. S2). The χ1 angle of D98 assumes 60° in the absence of 5HT, 90° in the intermediate and the occluded conformation, and oscillates between the two values in the outward-open conformation.
      The vertical line in Figure 7 indicates the time point when TM6a began to move toward the intermediate state (Fig. 7D) in replica 10. The positively charged amino group of 5HT started to interact with the negative helical dipole of TM6a (Fig. 7H) once the intermediate state was reached, most importantly with the carbonyl group of F335. This interaction was seen in all six trajectories that reached the intermediate state, and it always coincided with the movement of TM6a (Fig. S3). Thus, electrostatic interactions of the amino group of 5HT with hSERT led to both partial occlusion by movements of the bundle domain and to the stabilization of the intermediate state.

      Interactions of the hydrophobic lid led to complete occlusion of hSERT

      Three trajectories (green color) reached full occlusion of 5HT via bundle domain movements. Figure 4 showed that the step leading to the transition to the fully occluded state from the intermediate conformation was dominated by movements of TM6a, while movements of TM1b were very small. The interactions between 5HT and hSERT did not change during this step (Fig. 7). Of particular importance were the interactions between the residues of the outer hydrophobic gate. Figure 8, B and C report on the reduction in distances between the Cα atoms of F335 and I172 as well as between the Cα atoms of F335 and Y176 and show that their motions correlated with the motions of TM6a (Fig. 8A).
      Figure thumbnail gr8
      Figure 8Interactions leading to complete occlusion. Time evolution of distances between (A) TM6a-TM9up (reported from for comparison), (B) F335Cα-I172Cα, (C) F335Cα-Y176Cα, (D) F334H-A330O, and (E) F335H-A331O of the 5HT-bound hSERT and colored according to . The three vertical lines depict the time points at which the systems reached the fully occluded state. All lines represent running averages of 20 ns. Distances between I172, Y176, and F335 (dashed black lines with yellow background) shown for (F) the outward-open state, (G) the intermediate state, and (H) the fully occluded state. Structures showing the α-helical backbone hydrogen bond distances between 330 to 334 and 331 to 335 (green dashed lines with black bordering) for (I) the outward-open state, (J) the intermediate state, and (K) the fully occluded state. 5HT, serotonin; hSERT, human serotonin transporter; TM, transmembrane helix.
      These distances changed gradually during the transition from the outward-open to the intermediate state. In contrast, the transition leading to the fully occluded state was fast, as the distances changed by 0.1 to 0.2 nm in a very short time period (marked by vertical lines). The distances between the Cα atoms of the residues of the outer gate (I172, Y176, F335) are shown in Figure 8, FH and visualize the increasing contacts between these residues. The side chain of F335 interacted mainly with the side chain of I172 in the intermediate state. Upon transition to the fully occluded state, F335 moved closer to Y176, which led to direct interactions between all three residues of the outer hydrophobic gate.
      In the intermediate state, the α-helical structure of TM6a was in a strained conformation that is characterized by a partial loss of secondary structure. Figure 8, D and E depict the distances of two consecutive α-helical backbone hydrogen bonds in TM6a formed between residues A330 and F334 and between residues A331 and F335. These hydrogen bonds broke in the intermediate conformation and reformed upon reaching the fully occluded state. Representative structures of the outward-open, intermediate, and the fully occluded state are shown in Figure 8, I and J. Breaking of these backbone hydrogen bonds led to a partial unwinding of the α-helix, resulting in a bulging of the lowest turn that shortened the distances between F335 and residues I172 and Y176 across the outer vestibule.
      To further elucidate the changes in the outer vestibule during occlusion of hSERT, we examined key interactions of the outer gate and compared the structural changes to the number of water molecules in the outer vestibule. The dynamic changes of TM6a and its transient unwinding reduced the volume of the outer vestibule through increased interaction between residues of the outer hydrophobic gate. The transition to the intermediate state narrowed the outer vestibule, but it remained continuously water filled, thereby connecting the S1 to the extracellular side. The transition from the intermediate to the fully occluded state was accompanied by the closure of the extracellular hydrophobic gate (I172, Y176, and F335). Snapshots of the starting outward-open state and a representative fully occluded state (Fig. 9, A and B, respectively) illustrate that the outer gate salt bridge between R104 and E493 and the hydrophobic gate formed two layers of interactions that stabilized the closed vestibule of the fully occluded state and also separated the S1 from the extracellular environment. The salt bridge provided stability, while the hydrophobic gate formed a permeation barrier for water.
      Figure thumbnail gr9
      Figure 9Formation of the outer gate salt bridge and closure of the hydrophobic gate. Close up and front view of the central substrate binding site S1. A, at the start of the simulation (outward-open conformation); B, at the end of the simulation (fully occluded conformation). The outer gate salt bridge residues (R104 and E493) are highlighted in stick and light green, while the hydrophobic gate residues are colored in brown (I172, Y176, and F335). Water is colored by atom type. TM1, TM6, sodium, and 5HT are colored in red, mauve, dark blue, and by atom type, respectively. C, the radius of the outer vestibule for the outward-open, intermediate, and fully occluded conformation as a function of the Z-coordinate. The vertical lines indicate the size of a water molecule (dash-dot) and the thickness of 5HT (dashed), respectively. The green color indicates the region of the salt bridge, brown the hydrophobic gate. 5HT, serotonin; TM, transmembrane helix.
      Substrate and ions reach the S1 through the outer vestibule, which is first narrowed and then closed by movements of the bundle domain and the hydrophobic lid. We used the program HOLE (
      • Smart O.S.
      • Goodfellow J.M.
      • Wallace B.A.
      The pore dimensions of gramicidin A.
      ) to determine the minimal pore radius. The access path to the S1 was large enough for 5HT in the outward-open conformation (Fig. 9C). In the intermediate state, the pore radius became smaller than 5HT but was still sufficiently large for water to form a connected single file of water molecules reaching the S1. The outer vestibule closed completely in the fully occluded state, as the pore radius decreased further and became smaller than the size of a water molecule. Closing of the hydrophobic gate occurred by movement of the bundle domain and also by changes in side chain orientation, as residues of the hydrophobic lid (I172, Y176 and F335) and salt bridge (R104 and E439) reach over across the vestibule to interact with their respective counterparts.

      Partial dehydration of the S1 in hSERT is associated with occlusion

      Movements of the bundle domain (Fig. 10A) toward the scaffold domain closed the hydrophobic gate and were associated with a reduction in the number of water molecules within 0.64 nm of 5HT (Fig. 10B). The Pearson’s correlation coefficients between the movement of either TM1b-TM9up or TM6a-TM9up and the water count are 0.70 and 0.79, respectively. The other two trajectories reaching the occluded state showed similar behavior and correlation between the number of water molecules in the vestibule, the movements of the bundle domain, and the water densities. Two alternative driving forces may underlie water movements: (i) the bundle helices push the water out of the extracellular gate or (ii) water actively leaves the vestibule due to increasing local hydrophobicity, thereby de-wetting the vestibule and generating space for the subsequent movement of the bundle domain. The former case would be associated with an intermediate higher water density, while the latter would correlate with an intermediate lower water density. We quantified the water density within the outer vestibule, averaged using 1 ns wide time windows. Three representative time points before (533 ns), right at (573 ns), and immediately after the initial large movement of TM6a (585 ns) were selected and are indicated by colored vertical lines (Fig. 10, A and B). The respective structures are shown in Figure 10, CE. Water density is visualized using shades of blue, ranging from light cyan to dark blue representing low to high water density, respectively. Figure 10, CE shows a decrease in water density that preceded the movement of TM6a. The conformation of hSERT was the same at the time points of 533 ns and 573 ns as shown by the structural overlay of Figure 10D, but the water density at 573 ns was lower. Water density remained low during the conformational transition (Fig. 10E) until reaching the fully occluded state. We can infer from these data that water has a tendency to leave the increasingly hydrophobic environment of the outer vestibule. This lower water density then triggers movements of the bundle domain toward the scaffold domain, bringing the residues of the hydrophobic gate close together. These data therefore indicate that partial de-wetting of the outer vestibule contributes to the driving forces for reaching the fully occluded state. Analysis of the change of water density in context of the bundle movement of all three occluding systems are shown in Fig. S4.
      Figure thumbnail gr10
      Figure 10Partial dehydration of the S1 facilitates movements of the bundle domain. Thin lines and circles show all data points with a time resolution of 1 ns, and the thick lines depict running averages. A, distances from TM9up to TM1b and TM6a are colored in red and mauve, respectively. B, number of water molecules within 0.64 nm of 5HT. The three vertical lines colored in magenta, orange, and green indicate the time points (533, 573, and 585 ns, respectively) that are shown in panel CE. C, view into the outer vestibule of hSERT at 533 ns showing the density of water. D, conformation of hSERT and water density at 573 ns (orange). The conformation of hSERT at 533 ns is superimposed (purple) for comparison. E, conformation of hSERT at 585 ns (green) and the corresponding water density is shown. The conformation of hSERT at 573 ns is superimposed (orange). Water density was calculated per volume element using the VMD VolMap tool after fitting the trajectory to hSERT. Water density was averaged using a time window of 1 ns at the depicted time point. Replica 10 is used as a representative example; the data of the other occluding trajectories are similar. 5HT, serotonin; hSERT, human serotonin transporter; TM, transmembrane helix.

      Discussion

      Structures of hSERT have provided the first static snapshots of the transporter in four different conformations, which are described as outward-open, outward-occluded, occluded, and inward-open (
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      ,
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ,
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ). While providing important structural information, these structures lack information on the dynamic conformational changes associated with substrate transport. MD simulations, starting from these structures, can provide high resolution temporal and positional information on the time evolution of molecular processes, if occurring on the nanosecond to microsecond timescale. In this study, we capture the motions leading to the sequestering of substrate 5HT through occlusion of hSERT in a cholesterol-containing phospholipid bilayer that mimics the physiological cell membrane. The unbiased all atom MD simulations confirmed that the outward-open state of ion-bound hSERT is stable in the absence of 5HT, which is consistent with experimental data showing that the outward-facing conformation of SLC6 transporters is stable in the presence of the co-transported ions (
      • Tavoulari S.
      • Margheritis E.
      • Nagarajan A.
      • DeWitt D.C.
      • Zhang Y.-W.
      • Rosado E.
      • Ravera S.
      • Rhoades E.
      • Forrest L.R.
      • Rudnick G.
      Two Na+ sites control conformational change in a neurotransmitter transporter homolog.
      ,
      • Szöllősi D.
      • Stockner T.
      Investigating the mechanism of sodium binding to SERT using direct simulations.
      ,
      • Szöllősi D.
      • Stockner T.
      Sodium binding stabilizes the outward-open state of SERT by limiting bundle domain motions.
      ).
      The unbiased all atom MD simulations showed that occlusion is substrate triggered, as simulations reached full occlusion of three trajectories and partial occlusion in another three replicas. The current study suggests that while the cumulative movement of the bundle is a tilting toward the scaffold domain as described by the rocking bundle model (
      • Forrest L.R.
      • Zhang Y.-W.
      • Jacobs M.T.
      • Gesmonde J.
      • Xie L.
      • Honig B.H.
      • Rudnick G.
      Mechanism for alternating access in neurotransmitter transporters.
      • Forrest L.R.
      • Rudnick G.
      The rocking bundle: A mechanism for ion-coupled solute flux by symmetrical transporters.
      ), the motions of occlusion consist of a sequence of steps that deviate from the simple tilting movement.
      We consistently find that initial motions are limited to TM1b movements which correlate with interactions of the amino moiety of 5HT with D98. This engagement narrows the vestibule and initiates interactions between residues of the outer gate and 5HT. Essential for reaching the intermediate state is partial de-wetting, which is triggered by the increasing hydrophobicity of the outer vestibule and the engagement of the amino moiety of 5HT with the backbone of TM6a. In contrast to the mostly gradual movement of TM1b, the movements of TM6b consist of two large steps, resulting in its partial unwinding, which in turn allow for its increased interaction with the hydrophobic gate. Our data indicate that essential for hSERT occlusion are (i) the positive charge of the amino group of 5HT that attracts the bundle domain by electrostatic interactions and (ii) the hydrophobic part of 5HT that increases the hydrophobicity in the S1, leading to a de-wetting of the S1 and the outer vestibule. Based on these data, it is tempting to speculate whether the intermediate state could participate in interactions that discriminate substrates from nonsubstrates, thereby contributing to the high substrate selectivity of hSERT.
      The alternating access model requires that the extracellular gate closes before the intracellular gate opens (
      • Jardetzky O.
      Simple allosteric model for membrane pumps.
      ). Our simulations show that the occluded cryo-EM structure of hSERT, stabilized by ibogaine and an antibody (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ), represents a partially occluded conformation. It is similar to the intermediate conformation that retains a small access path to the S1. In contrast, the fully occluded hSERT conformation reached by three simulations, similar to the 5HT-bound cryo-EM structure (
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ), represents a state in which both gates are sealed: the extracellular gate reaches a conformation as observed in the inward-open cryo-EM structure of hSERT (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ) and the intracellular gate maintains a conformation as present in the outward-open structure of hSERT (
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      ). We infer from this finding that hSERT has only one entirely occluded state with two fully closed gates, while the outward-occluded (
      • Coleman J.A.
      • Yang D.
      • Zhao Z.
      • Wen P.-C.
      • Yoshioka C.
      • Tajkhorshid E.
      • Gouaux E.
      Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
      ) and the putative inward-occluded states are partially occluded states, stabilized by large inhibitors capable of preventing full occlusion.
      Our data show that the substrate 5HT is essential for inducing hSERT occlusion as it provides several key interactions: (i) strong stabilizing interactions with TM3 and TM8 by hydrogen bonds, (ii) fluctuating electrostatic interactions with D98 to attract TM1b, (iii) strong interactions with the negative helical dipole of TM6a to lock hSERT in the intermediate state, and (iv) increasing hydrophobic surface in the S1 to trigger de-wetting. This enables crossing of the energetic barrier associated with the structural re-arrangement and enables the tilting of the bundle domain. Thus, occlusion can occur as a result of the stochastic de-wetting of the S1 and the essential interaction of 5HT with hSERT. In support of our results, the conservative mutation D98E was shown to reduce 5HT uptake in rat SERT by 70% (
      • Barker E.L.
      • Moore K.R.
      • Rakhshan F.
      • Blakely R.D.
      Transmembrane domain I contributes to the permeation pathway for serotonin and ions in the serotonin transporter.
      ). Inhibition of 5HT uptake by dimetyltryptamine and gramine in WT rat SERT and in the D98E mutant showed that the precise distance between the carboxyl group of residue 98 and the amino group of the substrate is essential. The loss of function caused by the extension of the side chain residue 98 by the insertion of one methylene group (an aspartate to glutamate mutation) could be compensated for by shorting of the acyl chain of the ligand by one methylene group (
      • Barker E.L.
      • Moore K.R.
      • Rakhshan F.
      • Blakely R.D.
      Transmembrane domain I contributes to the permeation pathway for serotonin and ions in the serotonin transporter.
      ). We and others have shown that human and rat SERT can tolerate single methylation of the positively charged amino group of substrates, while double methylation can convert substrates to partial substrates, and triple methylation prevents transport (
      • Gobbi M.
      • Funicello M.
      • Gerstbrein K.
      • Holy M.
      • Moya P.R.
      • Sotomayor R.
      • Forray M.I.
      • Gysling K.
      • Paluzzi S.
      • Bonanno G.
      • Reyes-Parada M.
      • Sitte H.H.
      • Mennini T.
      N,N-dimethyl-thioamphetamine and methyl-thioamphetamine, two non-neurotoxic substrates of 5-HT transporters, have scant in vitro efficacy for the induction of transporter-mediated 5-HT release and currents.
      ,
      • Sandtner W.
      • Stockner T.
      • Hasenhuetl P.S.
      • Partilla J.S.
      • Seddik A.
      • Zhang Y.-W.
      • Cao J.
      • Holy M.
      • Steinkellner T.
      • Rudnick G.
      • Baumann M.H.
      • Ecker G.F.
      • Newman A.H.
      • Sitte H.H.
      Binding mode selection determines the action of ecstasy homologs at monoamine transporters.
      ,
      • Celik L.
      • Sinning S.
      • Severinsen K.
      • Hansen C.G.
      • Møller M.S.
      • Bols M.
      • Wiborg O.
      • Schiøtt B.
      Binding of serotonin to the human serotonin transporter. Molecular modeling and experimental validation.
      ,
      • Saha K.
      • Partilla J.S.
      • Lehner K.R.
      • Seddik A.
      • Stockner T.
      • Holy M.
      • Sandtner W.
      • Ecker G.F.
      • Sitte H.H.
      • Baumann M.H.
      ‘Second-generation’ mephedrone analogs, 4-MEC and 4-MePPP, differentially affect monoamine transporter function.
      ). Mutation of T439A (TM8) reduced the affinity of 5HT 4-fold, while replacing the T439-interacting hydroxyl group of 5HT with a methyl group lowered the affinity 20-fold (
      • Celik L.
      • Sinning S.
      • Severinsen K.
      • Hansen C.G.
      • Møller M.S.
      • Bols M.
      • Wiborg O.
      • Schiøtt B.
      Binding of serotonin to the human serotonin transporter. Molecular modeling and experimental validation.
      ). Structural change of hSERT substrates that alter key interactions by modifications of the substituents of the aromatic ring or methylation of the amine nitrogen frequently displayed a loss in affinity or a change in substrate/blocker profile (
      • Krishnamurthy H.
      • Gouaux E.
      X-ray structures of LeuT in substrate-free outward-open and apo inward-open states.
      ,
      • Sandtner W.
      • Stockner T.
      • Hasenhuetl P.S.
      • Partilla J.S.
      • Seddik A.
      • Zhang Y.-W.
      • Cao J.
      • Holy M.
      • Steinkellner T.
      • Rudnick G.
      • Baumann M.H.
      • Ecker G.F.
      • Newman A.H.
      • Sitte H.H.
      Binding mode selection determines the action of ecstasy homologs at monoamine transporters.
      ,
      • Saha K.
      • Partilla J.S.
      • Lehner K.R.
      • Seddik A.
      • Stockner T.
      • Holy M.
      • Sandtner W.
      • Ecker G.F.
      • Sitte H.H.
      • Baumann M.H.
      ‘Second-generation’ mephedrone analogs, 4-MEC and 4-MePPP, differentially affect monoamine transporter function.
      ,
      • Zeppelin T.
      • Ladefoged L.K.
      • Sinning S.
      • Schiøtt B.
      Substrate and inhibitor binding to the serotonin transporter: Insights from computational, crystallographic, and functional studies.
      ,
      • Niello M.
      • Cintulova D.
      • Hellsberg E.
      • Jäntsch K.
      • Holy M.
      • Ayatollahi L.H.
      • Cozzi N.V.
      • Freissmuth M.
      • Sandtner W.
      • Ecker G.F.
      • Mihovilovic M.D.
      • Sitte H.H.
      para-Trifluoromethyl-methcathinone is an allosteric modulator of the serotonin transporter.
      ,
      • Duart-Castells L.
      • Nadal-Gratacós N.
      • Muralter M.
      • Puster B.
      • Berzosa X.
      • Estrada-Tejedor R.
      • Niello M.
      • Bhat S.
      • Pubill D.
      • Camarasa J.
      • Sitte H.H.
      • Escubedo E.
      • López-Arnau R.
      Role of amino terminal substitutions in the pharmacological, rewarding and psychostimulant profiles of novel synthetic cathinones.
      ).
      In conclusion, our study shows that the bundle domain moves further than observed in the ibogaine-bound cryo-EM structure of hSERT and fully sequesters 5HT in the S1 from the extracellular side. This fully occluded conformation is comparable to the 5HT-bound cryo-EM structure of hSERT. The intracellular gate is closed as observed in the outward-open structure, and the extracellular gate assumes a geometry as present in the inward-open structure of hSERT. The contemporaneous closure of both the extracellular and intracellular gates represents the transition state that is essential for the alternating access model to reach the inward-open state. Water and the hydrophobicity of 5HT play an important role in the process of occlusion. Importantly, we find that the process of occlusion is an orchestrated sequence of events that allows for detection of substrates over nonsubstrates by specific interactions between the ligand and hSERT, thereby limiting the chemical space of potential substrates.

      Experimental procedures

      Model preparation

      For the purpose of this comprehensive work, the outward-open hSERT crystal structure with a resolution of 3.15 Å was used (PDB ID: 5I71, (
      • Coleman J.A.
      • Green E.M.
      • Gouaux E.
      X-ray structures and mechanism of the human serotonin transporter.
      )). Using the MODELLER 9.20 (
      • Webb B.
      • Sali A.
      Protein structure modeling with MODELLER.
      • Shen M.
      • Sali A.
      Statistical potential for assessment and prediction of protein structures.
      ), missing side chains and bound ions were positioned by generating 100 structures and picking out the best 10 based on the Descrete Optimized Protein Energy score and visual inspection. Since SERT is an integral membrane protein, the ten most realistic models with the lowest score previously created were embedded into a 1-palmitoyl-2-oleoylphosphatidylcholine and cholesterol containing membrane with 70:30 mol%.

      Manual docking of 5HT in hSERT

      Since no crystal structure of hSERT bound to 5HT is available, 5HT was docked manually into hSERT according to dopamine bound dDAT (PDB ID: 4XP1, (
      • Wang K.H.
      • Penmatsa A.
      • Gouaux E.
      Neurotransmitter and psychostimulant recognition by the dopamine transporter.
      )). Thus, the protein structures of dDAT and hSERT were aligned, and 5HT was oriented in a way that the ethylamine groups and parts of the ring moieties overlap. Afterward, the coordinates of 5HT were incorporated into the hSERT system. During energy minimization and equilibration, the endogenous substrate will take its favorable position. This pose is similar to a previously proposed model for 5HT binding to hSERT (
      • Celik L.
      • Sinning S.
      • Severinsen K.
      • Hansen C.G.
      • Møller M.S.
      • Bols M.
      • Wiborg O.
      • Schiøtt B.
      Binding of serotonin to the human serotonin transporter. Molecular modeling and experimental validation.
      ) and to the very recent cryo-EM structure of 5HT-bound hSERT (
      • Yang D.
      • Gouaux E.
      Illumination of serotonin transporter mechanism and role of the allosteric site.
      ), which was solved after completion of the simulation.

      Coarse-grained simulations

      The ten protein–membrane systems were converted into coarse-grained (CG) representations by applying the MARTINI force field (
      • Monticelli L.
      • Kandasamy S.K.
      • Periole X.
      • Larson R.G.
      • Tieleman D.P.
      • Marrink S.-J.
      The MARTINI coarse-grained force field: Extension to proteins.
      ,
      • de Jong D.H.
      • Singh G.
      • Bennett W.F.D.
      • Arnarez C.
      • Wassenaar T.A.
      • Schäfer L.V.
      • Periole X.
      • Tieleman D.P.
      • Marrink S.J.
      Improved parameters for the Martini coarse-grained protein force field.
      ,
      • Wassenaar T.A.
      • Ingólfsson H.I.
      • Böckmann R.A.
      • Tieleman D.P.
      • Marrink S.J.
      Computational lipidomics with insane: A versatile tool for generating custom membranes for molecular simulations.
      ). The simulation box was filled with water and 150 mM NaCl, the protein restrained to avoid conformational changes, and the CG systems simulated using Gromacs version 2019.2 (
      • Abraham M.J.
      • Murtola T.
      • Schulz R.
      • Páll S.
      • Smith J.C.
      • Hess B.
      • Lindahl E.
      GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers.
      ) for 1 μs each, yielding equilibrated lipid membranes surrounding the transporters.

      All atom simulations

      Equilibrated membranes, water, and ions were transformed to an atomistic (all atom) representation (
      • Wassenaar T.A.
      • Pluhackova K.
      • Böckmann R.A.
      • Marrink S.J.
      • Tieleman D.P.
      Going backward: A flexible geometric approach to reverse transformation from coarse grained to atomistic models.
      ). The CG-restrained protein structure was removed, and the originally generated transporter structures were incorporated into the membranes using the membed procedure (
      • Wolf M.G.
      • Hoefling M.
      • Aponte-Santamaría C.
      • Grubmüller H.
      • Groenhof G.
      g_membed: Efficient insertion of a membrane protein into an equilibrated lipid bilayer with minimal perturbation.
      ) to avoid overlapping atoms. The amber99sb-ildn force field (
      • Lindorff-Larsen K.
      • Piana S.
      • Palmo K.
      • Maragakis P.
      • Klepeis J.L.
      • Dror R.O.
      • Shaw D.E.
      Improved side-chain torsion potentials for the Amber ff99SB protein force field.
      ) was used to describe the protein, water, and ions properly. 1-Palmitoyl-2-oleoylphosphatidylcholine and cholesterol were described by Slipid (
      • Jämbeck J.P.M.
      • Lyubartsev A.P.
      An extension and further validation of an all-atomistic force field for biological membranes.
      ,
      • Jämbeck J.P.M.
      • Lyubartsev A.P.
      Another piece of the membrane puzzle: Extending slipids further.
      ). All atom MD simulations were performed using Gromacs 2019.3 (
      • Abraham M.J.
      • Murtola T.
      • Schulz R.
      • Páll S.
      • Smith J.C.
      • Hess B.
      • Lindahl E.
      GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers.
      ). After the assembly of the simulations systems, each system was energy minimized and equilibrated using a four-step protocol to slowly release the position restraints (1000, 100, 10, 1 kJ/mol/nm) active on the Cα atoms, the bound ions, and 5HT, if present. The MD production runs were performed for 2000 ns without any restraints for each independent system (10× 5HT-free and 10× 5HT-bound, 40,000 ns in total). To maintain the temperature at 310 K, the v-rescale (τ = 0.5 ps) thermostat (
      • Bussi G.
      • Donadio D.
      • Parrinello M.
      Canonical sampling through velocity rescaling.
      ) was used to independently couple protein, membrane, and solvent. The pressure was maintained at 1 bar by applying the Parrinello-Rahmen barostat (
      • Parrinello M.
      • Rahman A.
      Polymorphic transitions in single crystals: A new molecular dynamics method.
      ) in a semi-isotropic manner and applying a coupling constant of 20.1 ps. The Lennard Jones potentials with a cut-off of 0.9 nm were applied to describe van der Waals interactions. For describing long-range electrostatic interactions, the smooth particle mesh Ewald method (
      • Darden T.
      • York D.
      • Pedersen L.
      Particle mesh Ewald: An N log(N) method for Ewald sums in large systems.
      ) was used with a cut-off of again 0.9 nm. Long-range correction for energy and pressure was applied. Resulting coordinates of all atoms from solving Newton’s equation of motion with the Leap frog integrator were recorded at every 25 ps. The mdp files can be shared upon request.

      Data analysis

      The raw MD trajectories were processed and analyzed using Gromacs version 2019.3, MD Analysis version 0.19.2 (
      • Michaud-Agrawal N.
      • Denning E.J.
      • Woolf T.B.
      • Beckstein O.
      MDAnalysis: A toolkit for the analysis of molecular dynamics simulations.
      ,
      • Gowers R.
      • Linke M.
      • Barnoud J.
      • Reddy T.
      • Melo M.
      • Seyler S.
      • Domański J.
      • Dotson D.
      • Buchoux S.
      • Kenney I.
      • Beckstein O.
      MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations.
      ), the forces were calculated using the Force Distribution Analysis package (
      • Stacklies W.
      • Seifert C.
      • Graeter F.
      Implementation of force distribution analysis for molecular dynamics simulations.
      ), and plots were created using python/mathplotlib and the R package. The principal component analysis was performed with GROMACS 2019.6 on the Cartesian coordinates of all trajectories. The Cα atoms of the transmembrane helices were used as fitting group, and the backbone of transmembrane helices was selected for the covariance analyses.

      Data availability

      The simulations parameter files are contained within the manuscript. Raw data can be shared upon request from the corresponding author Thomas Stockner, email [email protected]

      Supporting information

      This article contains supporting information (Supporting Data, Figures S1–S4).

      Conflict of interest

      The authors declare that they have no conflicts of interest with the contents of this article.

      Author contributions

      R. G., H. H. S., and T. S. conceptualization; R. G. and D. S. methodology; M. N., H. H. S., T. S., R. G., and D. S. writing–original draft; E. L., M. N., H. H. S., T. S., R. G., and D. S. writing-review and editing; R. G. and D. S. software; M. N., H. H. S., T. S., R. G., and D. S. validation; R. G. and D. S. formal analysis, R. G. and D. S. investigation; R. G. and D. S. data curation; E. L. data analysis; H. H. S. and T. S. supervision; H. H. S. and T. S. funding acquisition; H. H. S. and T. S. resources.

      Funding and additional information

      The research underlying the current publication has been supported by the Austrian Science Fund (FWF) stand-alone grant P32017 to T. S. and H. H. S. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 860954. The results presented have been achieved using the Vienna Scientific Cluster.

      Supporting informtaion

      References

        • Fuller R.W.
        • Wong D.T.
        Serotonin uptake and serotonin uptake inhibition.
        Ann. N. Y. Acad. Sci. 1990; 600: 68-80
        • Murphy D.L.
        Serotonin transporter: Gene, genetic disorders, and pharmacogenetics.
        Mol. Interv. 2004; 4: 109-123
        • Chen N.-H.
        • Reith M.E.A.
        • Quick M.W.
        Synaptic uptake and beyond: The sodium- and chloride-dependent neurotransmitter transporter family SLC6.
        Pflugers Arch. 2004; 447: 519-531
        • Bhat S.
        • El-Kasaby A.
        • Freissmuth M.
        • Sucic S.
        Functional and biochemical consequences of disease variants in neurotransmitter transporters: A special emphasis on folding and trafficking deficits.
        Pharmacol. Ther. 2021; 222: 107785
        • Kristensen A.S.
        • Andersen J.
        • Jørgensen T.N.
        • Sørensen L.
        • Eriksen J.
        • Loland C.J.
        • Strømgaard K.
        • Gether U.
        SLC6 neurotransmitter transporters: Structure, function, and regulation.
        Pharmacol. Rev. 2011; 63: 585-640
        • James S.L.
        • Abate D.
        • Abate K.H.
        • Abay S.M.
        • Abbafati C.
        • Abbasi N.
        • Abbastabar H.
        • Abd-Allah F.
        • Abdela J.
        • Abdelalim A.
        • Abdollahpour I.
        • Abdulkader R.S.
        • Abebe Z.
        • Abera S.F.
        • Abil O.Z.
        • et al.
        Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: A systematic analysis for the Global Burden of Disease Study 2017.
        Lancet. 2018; 392: 1789-1858
        • Sitte H.H.
        • Freissmuth M.
        Amphetamines, new psychoactive drugs and the monoamine transporter cycle.
        Trends Pharmacol. Sci. 2015; 36: 41-50
        • Niello M.
        • Gradisch R.
        • Loland C.J.
        • Stockner T.
        • Sitte H.H.
        Allosteric modulation of neurotransmitter transporters as a therapeutic strategy.
        Trends Pharmacol. Sci. 2020; 41: 446-463
        • Rothman R.B.
        • Partilla J.S.
        • Baumann M.H.
        • Lightfoot-Siordia C.
        • Blough B.E.
        Studies of the biogenic amine transporters. 14. Identification of low-efficacy “partial” substrates for the biogenic amine transporters.
        J. Pharmacol. Exp. Ther. 2012; 341: 251-262
        • Hasenhuetl P.S.
        • Bhat S.
        • Freissmuth M.
        • Sandtner W.
        Functional selectivity and partial efficacy at the monoamine transporters: A unified model of allosteric modulation and amphetamine-induced substrate release.
        Mol. Pharmacol. 2019; 95: 303-312
        • Yamashita A.
        • Singh S.K.
        • Kawate T.
        • Jin Y.
        • Gouaux E.
        Crystal structure of a bacterial homologue of Na+/Cl--dependent neurotransmitter transporters.
        Nature. 2005; 437: 215-223
        • Krishnamurthy H.
        • Gouaux E.
        X-ray structures of LeuT in substrate-free outward-open and apo inward-open states.
        Nature. 2012; 481: 469-474
        • Gotfryd K.
        • Boesen T.
        • Mortensen J.S.
        • Khelashvili G.
        • Quick M.
        • Terry D.S.
        • Missel J.W.
        • LeVine M.V.
        • Gourdon P.
        • Blanchard S.C.
        • Javitch J.A.
        • Weinstein H.
        • Loland C.J.
        • Nissen P.
        • Gether U.
        X-ray structure of LeuT in an inward-facing occluded conformation reveals mechanism of substrate release.
        Nat. Commun. 2020; 11: 1005
        • Focke P.J.
        • Wang X.
        • Larsson H.P.
        Neurotransmitter transporters: Structure meets function.
        Structure. 2013; 21: 694-705
        • Penmatsa A.
        • Gouaux E.
        How LeuT shapes our understanding of the mechanisms of sodium-coupled neurotransmitter transporters.
        J. Physiol. 2014; 592: 863-869
        • Coleman J.A.
        • Green E.M.
        • Gouaux E.
        X-ray structures and mechanism of the human serotonin transporter.
        Nature. 2016; 532: 334-339
        • Coleman J.A.
        • Yang D.
        • Zhao Z.
        • Wen P.-C.
        • Yoshioka C.
        • Tajkhorshid E.
        • Gouaux E.
        Serotonin transporter-ibogaine complexes illuminate mechanisms of inhibition and transport.
        Nature. 2019; 569: 141-145
        • Yang D.
        • Gouaux E.
        Illumination of serotonin transporter mechanism and role of the allosteric site.
        Sci. Adv. 2021; 7eabl3857
        • Forrest L.R.
        • Zhang Y.-W.
        • Jacobs M.T.
        • Gesmonde J.
        • Xie L.
        • Honig B.H.
        • Rudnick G.
        Mechanism for alternating access in neurotransmitter transporters.
        Proc. Natl. Acad. Sci. U. S. A. 2008; 105: 10338-10343
        • Cheng M.H.
        • Bahar I.
        Molecular mechanism of dopamine transport by human dopamine transporter.
        Structure. 2015; 23: 2171-2181
        • Koldsø H.
        • Noer P.
        • Grouleff J.
        • Autzen H.E.
        • Sinning S.
        • Schiøtt B.
        Unbiased simulations reveal the inward-facing conformation of the human serotonin transporter and Na+ ion release.
        PLoS Comput. Biol. 2011; 7e1002246
        • Koldsø H.
        • Autzen H.E.
        • Grouleff J.
        • Schiøtt B.
        Ligand induced conformational changes of the human serotonin transporter revealed by molecular dynamics simulations.
        PLoS One. 2013; 8e63635
        • Hellsberg E.
        • Ecker G.F.
        • Stary-Weinzinger A.
        • Forrest L.R.
        A structural model of the human serotonin transporter in an outward-occluded state.
        PLoS One. 2019; 14e0217377
        • Sinning S.
        • Musgaard M.
        • Jensen M.
        • Severinsen K.
        • Celik L.
        • Koldsø H.
        • Meyer T.
        • Bols M.
        • Jensen H.H.
        • Schiøtt B.
        • Wiborg O.
        Binding and orientation of tricyclic antidepressants within the central substrate site of the human serotonin transporter.
        J. Biol. Chem. 2010; 285: 8363-8374
        • Thomas J.R.
        • Gedeon P.C.
        • Grant B.J.
        • Madura J.D.
        LeuT conformational sampling utilizing accelerated molecular dynamics and principal component analysis.
        Biophys. J. 2012; 103: L1-L3
        • Tavoulari S.
        • Margheritis E.
        • Nagarajan A.
        • DeWitt D.C.
        • Zhang Y.-W.
        • Rosado E.
        • Ravera S.
        • Rhoades E.
        • Forrest L.R.
        • Rudnick G.
        Two Na+ sites control conformational change in a neurotransmitter transporter homolog.
        J. Biol. Chem. 2016; 291: 1456-1471
        • Claxton D.P.
        • Quick M.
        • Shi L.
        • de Carvalho F.D.
        • Weinstein H.
        • Javitch J.A.
        • Mchaourab H.S.
        Ion/substrate-dependent conformational dynamics of a bacterial homolog of neurotransmitter:sodium symporters.
        Nat. Struct. Mol. Biol. 2010; 17: 822-829
        • Szöllősi D.
        • Stockner T.
        Investigating the mechanism of sodium binding to SERT using direct simulations.
        Front. Cell. Neurosci. 2021; 15: 673782
        • Smart O.S.
        • Goodfellow J.M.
        • Wallace B.A.
        The pore dimensions of gramicidin A.
        Biophys. J. 1993; 65: 2455-2460
        • Szöllősi D.
        • Stockner T.
        Sodium binding stabilizes the outward-open state of SERT by limiting bundle domain motions.
        Cells. 2022; 11: 255
        • Forrest L.R.
        • Rudnick G.
        The rocking bundle: A mechanism for ion-coupled solute flux by symmetrical transporters.
        Physiology. 2009; 24: 377-386
        • Jardetzky O.
        Simple allosteric model for membrane pumps.
        Nature. 1966; 211: 969-970
        • Barker E.L.
        • Moore K.R.
        • Rakhshan F.
        • Blakely R.D.
        Transmembrane domain I contributes to the permeation pathway for serotonin and ions in the serotonin transporter.
        J. Neurosci. 1999; 19: 4705-4717
        • Gobbi M.
        • Funicello M.
        • Gerstbrein K.
        • Holy M.
        • Moya P.R.
        • Sotomayor R.
        • Forray M.I.
        • Gysling K.
        • Paluzzi S.
        • Bonanno G.
        • Reyes-Parada M.
        • Sitte H.H.
        • Mennini T.
        N,N-dimethyl-thioamphetamine and methyl-thioamphetamine, two non-neurotoxic substrates of 5-HT transporters, have scant in vitro efficacy for the induction of transporter-mediated 5-HT release and currents.
        J. Neurochem. 2008; 105: 1770-1780
        • Sandtner W.
        • Stockner T.
        • Hasenhuetl P.S.
        • Partilla J.S.
        • Seddik A.
        • Zhang Y.-W.
        • Cao J.
        • Holy M.
        • Steinkellner T.
        • Rudnick G.
        • Baumann M.H.
        • Ecker G.F.
        • Newman A.H.
        • Sitte H.H.
        Binding mode selection determines the action of ecstasy homologs at monoamine transporters.
        Mol. Pharmacol. 2016; 89: 165-175
        • Celik L.
        • Sinning S.
        • Severinsen K.
        • Hansen C.G.
        • Møller M.S.
        • Bols M.
        • Wiborg O.
        • Schiøtt B.
        Binding of serotonin to the human serotonin transporter. Molecular modeling and experimental validation.
        J. Am. Chem. Soc. 2008; 130: 3853-3865
        • Saha K.
        • Partilla J.S.
        • Lehner K.R.
        • Seddik A.
        • Stockner T.
        • Holy M.
        • Sandtner W.
        • Ecker G.F.
        • Sitte H.H.
        • Baumann M.H.
        ‘Second-generation’ mephedrone analogs, 4-MEC and 4-MePPP, differentially affect monoamine transporter function.
        Neuropsychopharmacology. 2015; 40: 1321-1331
        • Zeppelin T.
        • Ladefoged L.K.
        • Sinning S.
        • Schiøtt B.
        Substrate and inhibitor binding to the serotonin transporter: Insights from computational, crystallographic, and functional studies.
        Neuropharmacology. 2019; 161: 107548
        • Niello M.
        • Cintulova D.
        • Hellsberg E.
        • Jäntsch K.
        • Holy M.
        • Ayatollahi L.H.
        • Cozzi N.V.
        • Freissmuth M.
        • Sandtner W.
        • Ecker G.F.
        • Mihovilovic M.D.
        • Sitte H.H.
        para-Trifluoromethyl-methcathinone is an allosteric modulator of the serotonin transporter.
        Neuropharmacology. 2019; 161: 107615
        • Duart-Castells L.
        • Nadal-Gratacós N.
        • Muralter M.
        • Puster B.
        • Berzosa X.
        • Estrada-Tejedor R.
        • Niello M.
        • Bhat S.
        • Pubill D.
        • Camarasa J.
        • Sitte H.H.
        • Escubedo E.
        • López-Arnau R.
        Role of amino terminal substitutions in the pharmacological, rewarding and psychostimulant profiles of novel synthetic cathinones.
        Neuropharmacology. 2021; 186: 108475
        • Webb B.
        • Sali A.
        Protein structure modeling with MODELLER.
        Methods Mol. Biol. 2014; 2199: 239-255
        • Shen M.
        • Sali A.
        Statistical potential for assessment and prediction of protein structures.
        Protein Sci. 2006; 15: 2507-2524
        • Wang K.H.
        • Penmatsa A.
        • Gouaux E.
        Neurotransmitter and psychostimulant recognition by the dopamine transporter.
        Nature. 2015; 521: 322-327
        • Monticelli L.
        • Kandasamy S.K.
        • Periole X.
        • Larson R.G.
        • Tieleman D.P.
        • Marrink S.-J.
        The MARTINI coarse-grained force field: Extension to proteins.
        J. Chem. Theory Comput. 2008; 4: 819-834
        • de Jong D.H.
        • Singh G.
        • Bennett W.F.D.
        • Arnarez C.
        • Wassenaar T.A.
        • Schäfer L.V.
        • Periole X.
        • Tieleman D.P.
        • Marrink S.J.
        Improved parameters for the Martini coarse-grained protein force field.
        J. Chem. Theory Comput. 2013; 9: 687-697
        • Wassenaar T.A.
        • Ingólfsson H.I.
        • Böckmann R.A.
        • Tieleman D.P.
        • Marrink S.J.
        Computational lipidomics with insane: A versatile tool for generating custom membranes for molecular simulations.
        J. Chem. Theory Comput. 2015; 11: 2144-2155
        • Abraham M.J.
        • Murtola T.
        • Schulz R.
        • Páll S.
        • Smith J.C.
        • Hess B.
        • Lindahl E.
        GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers.
        SoftwareX. 2015; 1–2: 19-25
        • Wassenaar T.A.
        • Pluhackova K.
        • Böckmann R.A.
        • Marrink S.J.
        • Tieleman D.P.
        Going backward: A flexible geometric approach to reverse transformation from coarse grained to atomistic models.
        J. Chem. Theory Comput. 2014; 10: 676-690
        • Wolf M.G.
        • Hoefling M.
        • Aponte-Santamaría C.
        • Grubmüller H.
        • Groenhof G.
        g_membed: Efficient insertion of a membrane protein into an equilibrated lipid bilayer with minimal perturbation.
        J. Comput. Chem. 2010; 31: 2169-2174
        • Lindorff-Larsen K.
        • Piana S.
        • Palmo K.
        • Maragakis P.
        • Klepeis J.L.
        • Dror R.O.
        • Shaw D.E.
        Improved side-chain torsion potentials for the Amber ff99SB protein force field.
        Proteins. 2010; 78: 1950-1958
        • Jämbeck J.P.M.
        • Lyubartsev A.P.
        An extension and further validation of an all-atomistic force field for biological membranes.
        J. Chem. Theory Comput. 2012; 8: 2938-2948
        • Jämbeck J.P.M.
        • Lyubartsev A.P.
        Another piece of the membrane puzzle: Extending slipids further.
        J. Chem. Theory Comput. 2013; 9: 774-784
        • Bussi G.
        • Donadio D.
        • Parrinello M.
        Canonical sampling through velocity rescaling.
        J. Chem. Phys. 2007; 126014101
        • Parrinello M.
        • Rahman A.
        Polymorphic transitions in single crystals: A new molecular dynamics method.
        J. Appl. Phys. 1981; 52: 7182-7190
        • Darden T.
        • York D.
        • Pedersen L.
        Particle mesh Ewald: An N log(N) method for Ewald sums in large systems.
        J. Chem. Phys. 1993; 98: 10089-10092
        • Michaud-Agrawal N.
        • Denning E.J.
        • Woolf T.B.
        • Beckstein O.
        MDAnalysis: A toolkit for the analysis of molecular dynamics simulations.
        J. Comput. Chem. 2011; 32: 2319-2327
        • Gowers R.
        • Linke M.
        • Barnoud J.
        • Reddy T.
        • Melo M.
        • Seyler S.
        • Domański J.
        • Dotson D.
        • Buchoux S.
        • Kenney I.
        • Beckstein O.
        MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations.
        in: Proceedings of the 15th Python in Science Conference. SciPy.org, 2016: 98-105https://doi.org/10.25080/Majora-629e541a-00e
        • Stacklies W.
        • Seifert C.
        • Graeter F.
        Implementation of force distribution analysis for molecular dynamics simulations.
        BMC Bioinformatics. 2011; 12: 101