The Structural Basis of the Activity Cliff in Modafinil-Based Dopamine Transporter Inhibitors

Abstract

Modafinil analogs with either a sulfoxide or sulfide moiety have improved binding affinities at the human dopamine transporter (hDAT) compared to modafinil, with lead sulfoxide-substituted analogs showing characteristics of atypical inhibition (e.g., JJC8-091). Interestingly, the only distinction between sulfoxide and sulfide substitution is the presence of one additional oxygen atom. To elucidate why such a subtle difference in ligand structure can result in different typical or atypical profiles, we investigated two pairs of analogs. Our quantum mechanical calculations revealed a more negatively charged distribution of the electrostatic potential surface of the sulfoxide substitution. Using molecular dynamics simulations, we demonstrated that sulfoxide-substituted modafinil analogs have a propensity to attract more water into the binding pocket. They also exhibited a tendency to dissociate from Asp79 and form a new interaction with Asp421, consequently promoting an inward-facing conformation of hDAT. In contrast, sulfide-substituted analogs did not display these effects. These findings elucidate the structural basis of the activity cliff observed with modafinil analogs and also enhance our understanding of the functionally relevant conformational spectrum of hDAT.

1. Introduction

The primary function of the dopamine transporter (DAT) is to regulate dopamine levels in the synaptic cleft by transporting the released dopamine back into the presynaptic dopaminergic neurons, thereby terminating dopamine neurotransmission [1,2]. Disruption of dopamine uptake and signaling has been associated with various psychiatric disorders, including schizophrenia, Parkinson’s disease, bipolar disorder, depression, attention deficit hyperactivity disorder, and substance use disorders [3,4,5]. Cocaine binds to all three monoamine transporters, but its psychostimulant actions have been primarily attributed to inhibition of DA reuptake via DAT [6,7,8]. Consequently, this inhibition results in a rapid accumulation of dopamine in the synaptic cleft, which is widely regarded as the primary cause of the reinforcing effects of cocaine, leading to addiction in humans [9]. It has been hypothesized that all compounds that inhibit DAT would induce behaviors similar to those caused by cocaine [10,11,12]. However, many studies have demonstrated that certain DAT inhibitors, including benztropine, modafinil, and their analogs, exhibit diminished rewarding effects, in humans and experimental animals, and thus reduced liability for misuse and addiction [13,14,15,16]. The term “atypical DAT inhibitors” has been designated to refer to such compounds, whereas substances such as cocaine and cocaine-like compounds are categorized as “typical DAT inhibitors” [13]. Atypical DAT inhibitors can mitigate the behavioral effects of cocaine and other typical DAT inhibitors in experimental animals, and may have therapeutic potential for treating cocaine use disorder [14].

DAT belongs to the neurotransmitter:sodium symporter family, the members of which translocate substrates across the cellular membrane by utilizing the energy stored in the transmembrane (TM) Na⁺ gradient [8]. In DAT, there are two sodium binding sites, each playing a distinct role. The first sodium ion (Na1) is closely associated with substrate binding, while the bound second sodium ion (Na2) is crucial in stabilizing the outward-facing state [17,18,19,20]. In the transport process, DAT undergoes conformational transitions by alternating between outward-facing (including outward-open and outward-occluded) and inward-facing (including inward-open and inward-occluded) states [21,22]. It has been found that the typical DAT inhibitors trap the transporter in the outward-facing states, while the atypical inhibitors prefer the inward-facing states [23].

X-ray crystallography has captured the outward-open, outward-occluded, and inward-open conformations of the Leucine transporter (LeuT), a bacterial homolog of DAT [24,25,26,27]. Various substrates and inhibitors have been co-crystallized with Drosophila melanogaster DAT (dDAT), which has a sequence identity of more than 50% with human DAT (hDAT) [28,29]. However, the available dDAT structures are all in outward-facing conformations. Based on the crystal structures of LeuT and dDAT, computational modeling and simulations have been instrumental in gaining insights into the molecular mechanisms of ligand binding and conformational dynamics in DAT [22,30,31,32,33,34,35,36,37]. Specifically, the mechanisms through which atypical inhibitors affect the DAT and the resulting conformational rearrangements in the protein have been explored through molecular dynamics (MD) simulations. Our previous study indicated that JHW007, a benztropine analog, preferred an inward-occluded DAT conformation [23,38]. The atypical effects of modafinil and its derivatives have been the subject of studies exploring the correlation between these effects and conformational changes in DAT [39,40]. JJC8-091, a modafinil analog that does not display any cocaine-like behaviors in rodents, favored a more occluded DAT conformation in comparison to its derivative, JJC8-088, a modafinil analog that is cocaine-like in rats and nonhuman primates [41,42]. Thus, MD simulations have proven valuable in unveiling these conformational changes induced by different inhibitors, providing insights into the molecular intricacies.

In this study, we carried out quantum mechanical calculations and extensive MD simulations to comparatively characterize the interactions between hDAT and two pairs of modafinil analogs, JJC8-091 and RDS04-010, with the sulfoxide substitution and their sulfide analogs JJC8-089 and RDS03-094, respectively. We performed extensive analysis of our MD simulation results to identify the different impacts of sulfoxide and sulfide substitutions on the conformational equilibrium of hDAT.

2. Methods

2.1. hDAT Binding Assay

hDAT-HEK293 cells were grown as previously described [43]. Briefly, HEK293 cells stably expressing human DAT were grown in Dulbecco’s modified Eagle medium (DMEM), supplemented with 5% fetal bovine serum, 5% calf bovine serum, 100 units/mL penicillin/100 μg/mL streptomycin, and 2 μg/mL puromycin and kept in an incubator at 37 °C and 10% CO₂. Upon reaching 80−90% confluence, cells were harvested using premixed Earle’s balanced salt solution with 5 mM ethylenediaminetetraacetic acid (EDTA) (Life Technologies, Carlsbad, CA, USA) and centrifuged at 3000 rpm for 10 min at 21 °C. The supernatant was removed, and the pellet was resuspended in 10 mL of hypotonic lysis buffer (5 mM MgCl₂, 5 mM Tris, pH 7.4 at 4 °C) and centrifuged at 14,500 rpm (∼25,000× g) for 30 min at 4 °C. The pellet was then resuspended in binding buffer (50 mM Tris, 120 mM NaCl, pH 7.4). Bradford protein assay (Bio-Rad, Hercules, CA, USA) was used to determine the protein concentration and the membranes were diluted to 1000 μg/mL and stored in a −80 °C freezer for later use. Radioligand binding assays were conducted similar to those previously described [44,45]. Experiments were conducted in 96-well polypropylene plates containing 50 μL of various concentrations of the test compound, diluted using 30% DMSO vehicle, 300 μL of binding buffer (50 mM Tris, 120 mM NaCl, pH 7.4), 50 μL of [³H]WIN35,428 [46] (final concentration 1.5 nM; K_d = 38.1 nM; NOVANDI Chemistry AB, 78 Ci/mmol SA), and 100 μL of membranes 30 μg/well). All compound dilutions were tested in triplicate and the competition reactions started with the addition of tissue; the plates were incubated for 120 min at 4 °C. Nonspecific binding was determined using a 10 μM final concentration of indatraline. Incubations were terminated by rapid filtration through PerkinElmer Uni-Filter-96 GF/C presoaked in 0.2% polyethylenimine, using a Brandel 96-Well Plates Harvester manifold (Brandel Instruments, Gaithersburg, MD, USA). The filters were washed a total of three times with 3 mL (3 × 1 mL/well or 3 × 1 mL/tube) of ice-cold binding buffer. After drying, 65 μL PerkinElmer MicroScint20 Scintillation Cocktail was added to each filter well. Plates were counted using a PerkinElmer MicroBeta Microplate Counter. For each experiment, aliquots of the prepared radioligand solutions were measured to calculate the exact amount of radioactivity added, taking in account the experimentally determined top-counter efficiency for the radioligand. K_i values have been extrapolated by constraining the bottom of the dose−response curves (=0% residual specific binding) in the nonlinear regression analysis. K_i values were calculated using GraphPad Prism 8 version 8.4.0 for Macintosh (GraphPad Software, San Diego, CA, USA) utilizing the One site-Fit K_i model. K_d values for the radioligands were determined via separate homologous competitive binding or radioligand binding saturation experiments. K_i values were determined from at least three independent experiments performed in triplicate and are reported as mean ± SEM.

2.2. Conformational Search

The conformational search for each compound was carried out with the Conformational Search module implemented in the Schrödinger suite (version 2022-3) using the OPLS4 force field [47] and customized parameters for the missing torsions of RDS03-094 and RDS04-010.

The conformation search was carried out with water as the solvent model, and the search was limited to 1000 steps. During the energy minimization stage, a Polak–Ribière conjugate gradient (PRCG) was performed, with a cutoff distance of 8.0 Å for van der Waals interactions, 20 Å for electrostatic interactions, and 4 Å for hydrogen bonds. The energy minimization step was performed with 100,000 iterations. The optimal pose for each compound was chosen as the conformation that exhibited the lowest energy after minimization.

2.3. Quantum Mechanical Calculations

The quantum mechanical (QM) calculations were performed using Jaguar in the Schrödinger suite. The optimized ligand geometries of each compound from the conformational search were used as inputs in QM calculations. The geometries were first optimized using the B3LYP-D3 theory and the 6-31G** basis set in Jaguar. The Poisson–Boltzmann finite element method (PBF) with water as solvent was also applied in the calculation. The same theory and basis set were used for the electrostatic potential calculation and the results were mapped onto the surfaces of constant electron density (Figure 1).

2.4. Optimization of Small-Molecule Parameters

To accurately simulate protein–ligand interaction without the missing torsions for compounds used in the MD simulation, the ligand parameters were systematically scanned using the force field builder in the Schrödinger suite. The identified missing torsion parameters of RDS03-094, RDS04-010, and JJC8-089 were then calculated using the force field builder with QM calculations.

2.5. Modeling and Molecular Docking

As RDS04-010, RDS03-094, and JJC8-089 share the same modafinil scaffold as JJC8-091, to establish the hDAT/RDS04-010, hDAT/RDS03-094, and hDAT/JJC8-089 models, we docked the corresponding compound to a well-equilibrated hDAT/JJC8-091 model from previous simulations [42]. The induced-fit docking (IFD) protocol from the Schrödinger suite (version 2022-3) was utilized to dock RDS04-010, RDS03-094, and JJC8-089 to the central binding pocket in a previously equilibrated hDAT/JJC8-091 model [42]. The center of mass of the bound JJC8-091 was selected as the center of the docking box. The standard IFD protocol, which can generate up to 20 poses for each docked compound, was performed, and the final poses were selected by considering the following criteria: (i) in the same orientation as the equilibrated JJC8-091, (ii) forming the ionic interaction between the protonated N atom of the molecule and Asp79, and (iii) the IFD score. The resulting complex models were then immersed back into the original lipid bilayer water box of hDAT/JJC8-091 to construct corresponding new simulation systems.

2.6. Molecular Dynamics Simulations

The MD simulations were performed using a Desmond MD engine (D. E. Shaw Research, New York, NY, USA) with our previously established protocols for the DAT simulation systems [23]. The OPLS3e force field [48] and customized ligands force field parameters were used in the MD simulation. Langevin dynamics was performed with an NPγT ensemble at a constant temperature (310 K) and 1 atm of constant pressure with the hybrid Nosé–Hoover Langevin piston method [49] on an anisotropic flexible periodic cell. The initial complexes systems after docking were first minimized and equilibrated with restraints on the ligand heavy atoms and protein backbone atoms. The restraints were removed during the stage of the production simulations (Table S1).

2.7. Structural Element Definitions

The structural elements of hDAT were adapted from our previous simulation [23] with the following definitions: NT (N-terminus, residues 58–64), TM1i (the intracellular segment (i) of TM1, residues 65–74), TM1m (the middle segment (m) of TM1, residues 75–82), TM1e (the extracellular segment (e) of TM1, residues 83–92), EL1 (the extracellular loop 1, residues 93–95), TM2e (residues 96–101), TM2m (residues 102–111), TM2i (residues 112–124), IL1 (residues 125–135), TM3i (residues 136–151), TM3m (residues 152–156), TM3e (residues 157–172), EL2 (residues 173–237), TM4e (residues 238–246), TM4i (residues 247–255), IL2 (residues 256–261), TM5i (residues 262–266), TM5m (residues 267–274), TM5e (residues 275–284), EL3 (residues 285–307), TM6e (residues 308–316), TM6m (residues 317–328), TM6i (residues 329–335), IL3 (residues 336–342), TM7i (residues 343–352), TM7m (residues 353–359), TM7e (residues 360–374), EL4a (residues 375–386), EL4b (residues 387–404), TM8e (residues 405–417), TM8m (residues 418–426), TM8i (residues 427–436), IL4 (residues 437–441), TM9i (residues 445–454), TM9e (residues 455–465), EL5 (residues 466–469), TM10e (residues 470–478), TM10m (residues 479–483), TM10i (484–496), IL5 (residues 497–518), TM11i (residues 519–529), TM11e (residues 530–540), EL6 (residues 541–557), TM12e (residues 558–572), TM12i (residues 573–584), and CT (C-terminus, residues 585–595). Note that we redefine the segment range of CT from residues 585–600 to residues 585–595 to avoid the impact of the flexible tail of the CT in the PIA analyses.

2.8. Conformational Analysis

The processing of trajectories and calculation of geometric measurements were performed using a combination of in-house Python scripts, MDAnalysis [50], and VMD [51].

The water molecules near the sulfur atoms of the ligands were analyzed using two steps. For each MD frame, the first step identified the water molecules within 3.5 Å of the sulfur atom of the ligand and counted their numbers. The second step calculated the minimum distance between any oxygen atom of the water molecules and the sulfur atom. Then the water counts and the distribution of the minimum distances from the representative MD frame ensemble for each simulated condition were plotted (Figure S4B).

To compute the water counts in the intracellular and extracellular vestibules, the MD frames were aligned to the same reference frame, using the stable transmembrane (TM) segments (TM3i, TM3m, TM3e, TM5m, TM5e, TM6e, TM6m, TM6i, TM7m, TM8e, TM8m, and TM8i) [23]. The volumes of the extracellular and intracellular vestibules were measured by computing the number of water molecules in each vestibule as described in our previous studies [23].

The quantitative analyses of the MD simulations are based on the bootstrapping from MD frame ensembles for each of the four simulated conditions. Specifically, for the two outward-facing conditions, hDAT/JJC8-089 and hDAT/RDS03-094, we pooled the frames from all trajectories for each condition to assemble the corresponding ensemble. For the two inward-facing conditions, hDAT/RDS04-010 and hDAT/JJC8-091, we identified the frames after the hDAT conformation had transitioned to the inward-facing state to assemble the corresponding ensemble. The inward-facing state was defined by large intracellular vestibule water counts. A total of ten bootstrapping samples, each containing 1000 frames, were generated from each ensemble for the PIA and other analyses shown in Figure 2, Figure 3 and Figures S2–S6, Table S2.

3. Results and Discussion

3.1. The Sulfoxide and Sulfide Modafinil Analogs Have Distinct Physical–Chemical Properties

The only difference between RDS03-094 and RDS04-010 (referred to as RDS pairs in this study) is the substitution of sulfoxide (-S(=O)-) versus sulfide (-S-). Similarly, JJC8-091 and JJC8-089 (referred to as JJC pairs) also have the same difference. The results of our binding assay showed that the sulfoxide analogs in both pairs have ~30-fold lower affinities at hDAT compared to the corresponding sulfide analogs (

Table 1. Binding affinities of the modafinil analogs.

Compound	Structure	hDAT Ki ± SEM (nM) a (n)
JJC8-089		69.1 ± 10.5 (4)
JJC8-091		2140 ± 110 (4)
RDS3-094		50.1 ± 6.51 (3)
RDS4-010		1570 ± 186 (3)
Modafinil		2280 ± 107 (3)

^a Each K_i value represents data from at least three independent experiments, each performed in triplicate. K_i values were obtained utilizing GraphPad Prism. Binding assay procedures are described in detail in Section 2. (n) = number of experiments.

). Compared to the JJC pairs, the RDS pairs have a 2,6-dimethyl substitution on the piperazine ring, which had no significant effect on DAT binding affinities (Table 1). To compare the impact of the sulfoxide and sulfide substitution on the compound conformation and electrostatic properties, we first carried out a conformational search of each compound to identify their conformers with the lowest energy. Our results showed that although the lowest-energy conformers of these two pairs of analogs share significant similarities in their extended conformations, they demonstrate noticeable differences near their aligned sulfoxide and sulfide moieties.

Using these lowest-energy conformers as the inputs, we then further optimized the geometry of these conformers and characterized their electrostatic potential surfaces (EPSs) with quantum mechanical (QM) calculations. Using the RDS pair as an example, the lowest-energy conformers show a relatively extended conformation similar to those resulting from the conformational search (Figure 1C,D). As expected, the most prominent difference on their EPS is the negative partial charge distributed around the oxygen atom of the sulfoxide (RDS04-010), which is absent in the sulfide (RDS03-094) (Figure 1A,B). These different charge distributions lead to varied orientations between the bisphenyl and 2,6-dimethylpiperazine moieties of these two compounds. Based on the optimized geometries of RDS04-010 and RDS03-094 from the QM optimization, the dihedral angle between those two moieties can vary from −35.8° to −46.0° due to the presence of sulfoxide substitution (Figure 1C,D). The differences in both charge distribution and ligand conformation may result in a varied propensity to form either ligand–water or ligand–protein interactions (see below).

Figure 1. The addition of the sulfoxide moiety significantly alters the physical–chemical properties of the ligand. The calculated electrostatic potential surfaces for RDS04-010 (A) and RDS03-094 (B) show that the sulfoxide introduces a partial negative charge (red) near the bisphenyl moiety. As a result, the dihedral angle between the bisphenyl and dimethylpiperazine moieties has >10° difference between RDS04-10 (C) and RDS03-094 (D).

Figure 1. The addition of the sulfoxide moiety significantly alters the physical–chemical properties of the ligand. The calculated electrostatic potential surfaces for RDS04-010 (A) and RDS03-094 (B) show that the sulfoxide introduces a partial negative charge (red) near the bisphenyl moiety. As a result, the dihedral angle between the bisphenyl and dimethylpiperazine moieties has >10° difference between RDS04-10 (C) and RDS03-094 (D).

3.2. The Sulfoxide- and Sulfide-Containing Ligands Have Different Binding Poses at hDAT

To evaluate the differential impact of the sulfoxide versus sulfide substitution on the binding of the ligands with the same scaffold at hDAT, we docked the ligands into an hDAT model in an outward-facing conformation, immersed the resulting hDAT/ligand complexes individually in an explicit lipid bilayer/water environment, and performed prolonged molecular dynamics (MD) simulations (see Section 2, Table S1, and Figure 2A). One limitation of MD simulations is that they can be trapped in local energy minima along the evolution of the trajectory. To address this challenge, we collected multiple trajectories for each complex (see Table S1).

Figure 2. The bound of RDS04-010 cannot stably interact with Asp79, resulting in the dissociation of Na2. (A) hDAT in complex with RDS04-010 is embedded in a lipid bilayer (pale cyan) with the central ligand binding pocket indicated by a red box. The initial binding poses of RDS04-010 and RDS03-094 from docking are shown in (B) and (E), respectively, and the representative snapshots of their equilibrated binding poses are shown in (D) and (G), respectively. During the MD simulations, as shown in panels (B–D), from a representative trajectory, RDS04-010 gradually loses its interaction with Asp79, which is present in the starting model (B), and gets close to Asp421 at the end of the trajectory (D); the evolution of the distances between the pyramidal nitrogen of the ligand and these two residues are shown in panel (C). The dissociation of Na2 from its binding site, which is demonstrated by its elongated distance to Asp421, is a prerequisite for the transition to happen. In comparison, such a transition does not happen in the simulations of hDAT bound with RDS03-094 (E–G). The sodium ions (Na1 and Na2) are shown in purple. For simplicity, the stably bound chloride ion is not illustrated. Note that the distance between the ligand pyramid nitrogen and Asp421 or Asp79 was the minimum distance between the pyramidal nitrogen and the sidechain oxygen atoms. Similarly, the distance between Asp421 and Na2 was the minimum distance between Na2 and the sidechain oxygen atoms of Asp. The ionic interactions are indicated by the dotted lines; the red arrow indicates the conformational change of Phe320 from (B–D).

Figure 2. The bound of RDS04-010 cannot stably interact with Asp79, resulting in the dissociation of Na2. (A) hDAT in complex with RDS04-010 is embedded in a lipid bilayer (pale cyan) with the central ligand binding pocket indicated by a red box. The initial binding poses of RDS04-010 and RDS03-094 from docking are shown in (B) and (E), respectively, and the representative snapshots of their equilibrated binding poses are shown in (D) and (G), respectively. During the MD simulations, as shown in panels (B–D), from a representative trajectory, RDS04-010 gradually loses its interaction with Asp79, which is present in the starting model (B), and gets close to Asp421 at the end of the trajectory (D); the evolution of the distances between the pyramidal nitrogen of the ligand and these two residues are shown in panel (C). The dissociation of Na2 from its binding site, which is demonstrated by its elongated distance to Asp421, is a prerequisite for the transition to happen. In comparison, such a transition does not happen in the simulations of hDAT bound with RDS03-094 (E–G). The sodium ions (Na1 and Na2) are shown in purple. For simplicity, the stably bound chloride ion is not illustrated. Note that the distance between the ligand pyramid nitrogen and Asp421 or Asp79 was the minimum distance between the pyramidal nitrogen and the sidechain oxygen atoms. Similarly, the distance between Asp421 and Na2 was the minimum distance between Na2 and the sidechain oxygen atoms of Asp. The ionic interactions are indicated by the dotted lines; the red arrow indicates the conformational change of Phe320 from (B–D).

The inspection of the equilibrated hDAT/ligand complex models from the MD simulations indicates that these two substitutions in both the RDS and JJC pairs result in similarly different interactions within the binding pocket. While these four ligands commonly interact with 26 residues in TMs1, 3, 6, 8, and 10 (Table S2), overall, the ligand binding pose tends to be more extended for the ligands with the sulfoxide compared to those with the sulfide. Consequently, Gly425 of TM8 interact more frequently with the sulfoxide-containing ligands, while Ala81 of TM1, Asn157 of TM3, and Ala479 of TM10 interact more frequently with the sulfide-containing ligands. In addition, a direct ionic interaction was formed persistently between the pyramidal nitrogen of the ligands (L_N) and the negatively charged Asp79, in hDAT/RDS03-094 and hDAT/JJC8-089, throughout the simulations (Figure 2E–G and Figure S1D–F), while this interaction was gradually lost during the hDAT/RDS04-010 and hDAT/JJC8-091 simulations (Figure 2B–D and Figure S1A–C). The divergence at this critical interaction is highly likely associated with the higher affinities of RDS03-094 and JJC8-089 compared to RDS04-010 and JJC8-091, respectively (Table 1).

In DAT, there are two sodium binding sites, each playing a distinct role. The first sodium ion (Na1) is closely associated with substrate binding, while the bound second sodium ion (Na2) is crucial in stabilizing the outward-facing state [17,18,19,20]. The hDAT/RDS04-010 and hDAT/JJC8-091 simulations showed a common tendency for the more intracellular-bound Na⁺ ion near the ligand binding pocket (Na2) to escape to the intracellular milieu (Figure 2B,D,E,G), which exposed a negatively charged Na2 coordinating residue, Asp421, to potentially attract and interact with the L_N of ligands. Note that such a Na2-escaping tendency has been previously associated with the transition to the inward-facing conformations in hDAT and homologous transporters [17,19,31,32].

To further understand the coordination among these structural elements, we compared the evolution of the distances between the L_N of the ligands and Asp79 and Asp421, as well as that between Na2 and Asp421 (Figure 2C,F). In the hDAT/RDS04-010 simulations, in eight out of nine trajectories, Na2 was not stable and the ionic interaction between Na2 and Asp421 gradually dissociated (i.e., the distance > 5 Å). The escape of Na2 from the binding pocket resulted in an opening at the bottom of the ligand binding site (Figure 2), while the L_N of RDS04-010 gradually dissociated with Asp79 before forming a water-mediated interaction with Asp421. In contrast, in the hDAT/RDS03-094 simulations, Na2 remained stably bound in the Na2 binding site, while the ionic interactions of Na2-Asp421 and L_N-Asp79 persistently formed. The same coordination among these elements was also observed in the hDAT/JJC8-091 but not in the hDAT/JJC8-089 simulations (Figure S1B,E).

3.3. The More Extended Pose of RDS04-010 Is Associated with the More Polar Nature of the Sulfoxide

The conformational transition between the inward-facing and outward-facing conformations of transporter proteins could be characterized by the changes in the volumes of the intracellular and extracellular vestibules (IV and EV, respectively). Previously, we developed a protocol to quantitatively estimate these volumes by the numbers of water molecules detected in these vestibules (referred to as “water count”, see Section 2) [17,23,52]. Using this protocol, our analysis showed that the representative frame ensemble from the hDAT/RDS03-094 condition has an average IV water count of 31.3, while that of hDAT/RDS04-010 has an average IV water count of 54.2. Therefore, the hDAT/RDS03-094 condition exhibited a notably smaller IV volume compared to hDAT/RDS04-010. On the extracellular side, hDAT/RDS03-094 showed a significantly higher average EV water count of 73.7, in contrast to the 62.0 observed in hDAT/RDS04-010 (Figure 3). A similar trend was also evident in the JJC pair (Figure S2). These findings further confirmed that RDS03-094 stabilized the outward-facing state in hDAT, while RDS04-010 promoted a transition towards the inward-facing state.

Figure 3. The outward-facing and inward-facing hDAT conformations stabilized by RDS04-010 and RDS03-094, respectively, show drastically different extracellular and intracellular vestibule volumes. We measured the volumes by counting the water molecules in the vestibules (see Section 2) and found that hDAT/RDS03-094 (orange) has ~11.7 more waters in the extracellular vestibule (A) and ~22.8 fewer waters in the intracellular vestibule (B) than hDAT/RDS04-010 (blue). These differences are demonstrated in the hDAT/RDS04-010 (C) and hDAT/RDS03-094 (D) models, with the water molecules shown in red spheres.

Figure 3. The outward-facing and inward-facing hDAT conformations stabilized by RDS04-010 and RDS03-094, respectively, show drastically different extracellular and intracellular vestibule volumes. We measured the volumes by counting the water molecules in the vestibules (see Section 2) and found that hDAT/RDS03-094 (orange) has ~11.7 more waters in the extracellular vestibule (A) and ~22.8 fewer waters in the intracellular vestibule (B) than hDAT/RDS04-010 (blue). These differences are demonstrated in the hDAT/RDS04-010 (C) and hDAT/RDS03-094 (D) models, with the water molecules shown in red spheres.

By extracting and comparing the frame ensembles corresponding to the inward-facing and outward-facing states, we further compared the key structural features of the RDS04-010 and RDS03-094 molecules bound in their preferred inward-facing and outward-facing hDAT conformations. We found that the distribution of the ligand end-to-end distances has the peak value of 5.6 Å for RDS03-094 in the outward-facing hDAT conformational ensemble and 7.3 Å for RDS04-010 in the inward-facing one (Figure S3). This difference is consistent with a more extended binding pose of RDS04-010, which has the sulfoxide substitution, than that of RDS03-094.

Furthermore, since the sulfoxide substitution carries more negative charge compared to sulfide (see Figure 1), more water molecules were attracted to the sulfoxide of the ligands (<3.5 Å) in the simulations than to sulfide, as the sulfoxide but not the sulfide can form a direct H-bond with these water molecules (Figure S4). Thus, the presence of oxygen in the sulfoxide not only impacts its charge distribution and overall ligand conformations, but also enhances the tendency to foster a polar environment in the binding pocket, which ultimately contributes to the more extended ligand binding poses in hDAT (Figure 2 and Figure S1).

3.4. The RDS04-010 Stabilized hDAT Conformation Shows Key Features of Inward-Facing State

We then further characterized the overall conformations of the TM domain of hDAT using the protein interaction analyzer (PIA) [53,54], which provides a quantitative assessment of conformational differences between two conditions without relying on protein superposition. To achieve this goal, we have previously defined extracellular, middle, and intracellular subsegments for each TM segment of DAT (see Section 2 for subsegment definitions) (Figure 4A) [23,53].

Figure 4. The binding of RDS04-010 induces hDAT to transition from the outward-facing to the inward-facing conformation. Panel (A) shows the division of TM segments of hDAT into the extracellular (“e”), middle (“m”), and intracellular (“i”) subsegments. The differences in the distances among the subsegments shows that compared to the hDAT model bound with RDS03-094 (hDAT/RDS03-094), hDAT/RDS04-010 has many subsegments moving closer to each other on the extracellular side (reddish color in panel (B)) and moving away from each other on the intracellular side (blueish color in panel (C)), thus transitioning to an inward-facing conformation. On panels (D,E), the equilibrated hDAT/RDS04-010 (light blue) and hDAT/RDS03-094 (light pink) complexes are superimposed. On the extracellular side (D), TM1e and TM6e and their neighboring subsegments move closer to TM10e, while on the intracellular side (E), TM1i swings out towards the membrane, with a coordinated rearrangement of TMs5i, 6i, and 7i. Panel (F) shows the Cα–Cα distance between Ile67 of TM1i and Leu583 of TM12i (red in panel (G)). Leu583 is relatively stable and selected as a reference for distance measurement. This distance was increased in the hDAT/RDS04-010 simulation, but not in the hDAT/RDS03-094 one.

Figure 4. The binding of RDS04-010 induces hDAT to transition from the outward-facing to the inward-facing conformation. Panel (A) shows the division of TM segments of hDAT into the extracellular (“e”), middle (“m”), and intracellular (“i”) subsegments. The differences in the distances among the subsegments shows that compared to the hDAT model bound with RDS03-094 (hDAT/RDS03-094), hDAT/RDS04-010 has many subsegments moving closer to each other on the extracellular side (reddish color in panel (B)) and moving away from each other on the intracellular side (blueish color in panel (C)), thus transitioning to an inward-facing conformation. On panels (D,E), the equilibrated hDAT/RDS04-010 (light blue) and hDAT/RDS03-094 (light pink) complexes are superimposed. On the extracellular side (D), TM1e and TM6e and their neighboring subsegments move closer to TM10e, while on the intracellular side (E), TM1i swings out towards the membrane, with a coordinated rearrangement of TMs5i, 6i, and 7i. Panel (F) shows the Cα–Cα distance between Ile67 of TM1i and Leu583 of TM12i (red in panel (G)). Leu583 is relatively stable and selected as a reference for distance measurement. This distance was increased in the hDAT/RDS04-010 simulation, but not in the hDAT/RDS03-094 one.

In comparing the hDAT conformation stabilized by RDS04-010 and that stabilized by RDS03-094, we found that the hDAT/RDS04-010 condition has many smaller pairwise distances among the extracellular subsegments (reddish pixels in Figure 4B), while having many larger distances on the intracellular side compared to hDAT/RDS03-094 (Figure 4C). Specifically, in hDAT/RDS04-010, when compared to hDAT/RDS03-094, the TM subsegments TM1e, TM6e, TM7e, TM11e, EL3, and EL4 moved closer to TM10e on the extracellular side, while TM1i swung out towards the lipid membrane with coordinated rearrangements of TM5i, TM6i, and TM7i on the intracellular side (Figure 4D,E). Thus, hDAT/RDS04-010 has clearly transitioned into an inward-facing conformation, while hDAT/RDS03-094 remained in an outward-facing conformation. A similar trend of these rearrangements was also observed between the hDAT/JJC8-091 and hDAT/JJC8-089 conditions (Figure S5).

The swinging-out of TM1i in the transition towards the inward-facing state, which has been previously observed in the LeuT simulations [55], prepares the transporter to release the substrate to the intracellular side. To quantitively evaluate the potential impact of this TM1i rearrangement, we measured the Cα-Cα distance between Ile67 of TM1i and Leu583 of TM12i, and found it fluctuated near ~45 Å in the hDAT/RDS03-094 condition but increased to ~55 Å in hDAT/RDS04-010 (Figure 4F,G). This ~10 Å outward movement of TM1i in hDAT/RDS04-010 was coordinated with the detachment of the N-terminal loop from the TM domain, and created an opening below the ligand binding pocket to the intracellular water milieu (Movie S1). Such a process was associated with the tendency of RDS04-010 to be in a more extended conformation (Figure 2). In line with this study, our prior long-timescale MD simulations demonstrated that a cocaine analog, CFT, stabilized hDAT in the outward-facing conformation, while an atypical DAT inhibitor, JHW007, steered hDAT towards an inward-facing conformation [23]. However, in that study, the transition to the inward-facing state did not progress far enough to allow for the outward swinging of the TM1i.

On the extracellular side of the ligand binding pocket, the gating residue, Phe320, demonstrated disparate behaviors in the inward-facing and outward-facing conformations. Its χ1 rotamer in the hDAT simulations bound with the sulfide-substituted ligands (RDS03-094 and JJC8-089) was populated more in gauche+, while in hDAT/RDS04-010 and hDAT/JJC8-091, this rotamer was in trans and closed the binding pocket from the extracellular side (Figure S6).

In addition, in the simulations of both the hDAT/RDS04-010 and hDAT/JJC8-091 conditions, similarly to what has been observed in LeuT [55], we found that the release of Na2 was associated with an outward movement of the TM1i during the transition to the inward-facing conformation.

Taken together, our simulations suggest that ligands with the sulfoxide substitution (RDS04-010 and JJC8-091) induce and stabilize hDAT in an inward-facing conformation, predicting that these compounds would be atypical DAT inhibitors. In contrast, those analogs with the sulfide (RDS03-094 and JJC8-089) prefer an outward-facing conformation and may have a more typical cocaine-like behavioral profile. In the case of JJC8-091, rodent studies support these predictions [42]. Extensive behavioral investigation of the other analogs is underway.

4. Conclusions

This study focused on two pairs of modafinil analogs, JJC8-091 and JJC8-089, and RDS04-010 and RDS03-094, aiming to investigate the impacts of the sulfoxide versus sulfide groups in affecting conformational changes to hDAT that translate into divergent behavioral profiles. Our ligand conformational search and QM calculations revealed differences in charge distribution and ligand conformation due to these sulfoxide or sulfide moieties. The results of our long MD simulations of hDAT/ligand complexes illustrated that these substitutions prompted the ligand to engage in distinct interactions with hDAT. We found that sulfoxide analogs RDS04-010 and JJC8-091 tend to promote an inward-facing conformation compared to the sulfide analogs of RDS03-094 and JJC8-089. In the transition to the inward-facing state, the initial events involved the Na2 escaping from its binding site, followed by the dissociation of the interaction between Asp79 and the pyramidal nitrogen of the ligand. Subsequently, the pyramidal nitrogen was found to establish an interaction with Asp421. This process is associated with an outward movement of the TM1i that opens the intracellular vestibule to the water milieu.

Thus, the seemingly subtle substituent difference between the sulfoxide versus sulfide in these modafinil analogs strongly affects their preference for either inward-facing or outward-facing DAT conformation, leading to what is known as an activity cliff. Indeed, it has been found in many instances that similarly minor modifications of small-compound ligands can lead to drastic changes in function at the target proteins, namely, activity cliffs. For example, we have reported that introducing an additional methyl group to a synthetic cannabinoid receptor agonist leads to a substantial improvement in its efficacy at cannabinoid receptor 1 [56]. A comparable subtle alteration in characteristics has also been observed among some 4-phenylpiperazine ligands at the dopamine D3 receptor, wherein the (R)-enantiomers are antagonists or weak partial agonists, while the (S)-enantiomers exhibit considerably greater efficacy [54].

Taken together, our findings of the disparate characteristics of sulfoxide versus sulfide analogs, along with their atomistic interaction details with hDAT, provide insights into the molecular mechanisms underlying the distinct pharmacological properties of these inhibitors. While our computational approach allows exploration of dynamic behavior and interactions over time, providing insights into the sequence of events in their mechanisms of action, the demanding nature of the extensive simulations and analysis required for robust conclusions highlights the non-trivial task of predicting whether a compound functions as a typical or atypical DAT inhibitor. Nonetheless, the gained insights will deepen our mechanistic understanding of the DAT function at the atomistic level by fully revealing its functionally relevant conformational spectrum, which will provide the groundwork for the rational design of drugs targeting DAT.