Molecular dynamics simulation of d-Benzedrine transmitting through molecular channels within D3R

Dex-Benzedrine (known as d-Benzedrine or SAT) acts in dopamine receptors of central nerve cell system. In clinic, SAT is used to treat a variety of diseases; meanwhile, it has dependence and addiction. In order to investigate the pharmacology and addiction mechanisms of SAT as a medicine, in this paper, we have studied the structure of D3R complex protein with SAT, and based on which, using potential mean force with umbrella samplings and the simulated phospholipid bilayer membrane (or POPC bilayer membrane), the molecular dynamics simulation was performed to obtain free energy changes upon the trajectories for SAT moving along the molecular channels within D3R. The free energy change for SAT transmitting toward the outside of cell along the functional molecular channel within D3R is 83.5 kJ mol−1. The change of free energy for SAT to permeate into the POPC bilayer membrane along the protective molecular channel within D3R is 87.7 kJ mol−1. Our previous work gave that the free energy for Levo-Benzedrine (RAT) transmitting toward the outside of cell along the functional molecular channel within D3R is 91.4 kJ mol−1, while it is 117.7 kJ mol−1 for RAT to permeate into the POPC bilayer membrane along the protective molecular channel within D3R. The values of free energy suggest that SAT relatively prefers likely to pass through the functional molecular channel within D3R for increasing the release of dopamine molecules resulting in a variety of functional effects for SAT. The obtained results show that the pharmacology and addiction mechanisms of SAT as a drug are closely related to the molecular dynamics and mechanism for SAT transmitting along molecular channels within D3R.


Introduction
Dex-Benzedrine (known as D-Benzedrine or SAT), its molecular structure shown in Figure 1, is clinically used to treat narcolepsy, fatigue, chronic alcoholism, and obesity. SAT is a kind of central stimulants, acting in dopamine receptors of central nerve cell system, which also has dependence and addiction. The pharmacological effects of SAT is considered as the same as Benzedrine, but its excitant effects in central nervous system and its inhibition of appetite both are stronger than those from Benzedrine, equivalent to two times of Benzedrine. The inhibition of appetite by SAT is equivalent to 3-4 times by Levo-Benzedrine (known as L-Benzedrine or RAT).
In central nervous cell system, there are five major DA receptor subtypes (D1-D5 receptors) all belonging to the members of a superfamily of proteins called G-protein-coupled receptors (GPCRs) which are all composed of seven transmembranes (7-GM, or TM). At the beginning of 1990s, with the development of biological cloning technology, five distinct DA receptors (D1-D5 receptors) were cloned and defined as D 1 R, D 2 R, D 3 R, Because there is a cliff gap between nerve cells and nerve cells, like a crack between two cliffs, SAT molecules skip the crack through the nerve cells at a prominent small cliff called synapses for passing through the gap in order to play a functional role. The main ingredients of synapses include all kinds of dopamine subtype receptors. Within the dopamine receptor structures there are molecular channels (Bian, Shi, Chi, & Xu, 2012;Zhang, Bian, Shi, & Xu, 2014). Molecular channels are divided into two categories defined, respectively, as (1) functional molecular channel: to transfer molecules and to play effects of function; (2) protective molecular channel: to prevent excess molecular function from molecular damage and then to play a function of protective body. Since SAT acts in dopamine receptors of the central nervous cell system, it also should pass through molecular functional channel within dopamine receptors to play its functional effects, or pass through protective molecular channel within dopamine receptors to prevent excess molecular function from molecular damage. Therefore, to study SAT transmitting through molecular channels within D 3 R would help to understand the role of SAT's pharmacological mechanisms and kinetics in nerve cell systems.
In this paper, the movement of SAT is characterized by free energy changes for SAT transmitting along molecular channels within dopamine receptors. The free energies for SAT moving along the molecular channels can be experimentally measured or obtained by molecular dynamics simulations. Based on the protein structure of D 3 R and the simulated phospholipid bilayer membrane (or POPC bilayer membrane), we will study and obtain the change of free energy for SAT transmitting along the molecular channels within D 3 R using potential mean force (PMF) with umbrella samplings to probe the mechanism and molecular dynamics of SAT moving along the molecular channels.

Materials and methods
Due to difficulty of crystallization and problem of stability for dopamine receptor proteins, among them only a crystal protein structure of D 3 R currently was reported with a protein code of 3PBL (Chien et al., 2010;De Paulis, Hall, & Ogren, 1985;Griffon, Pilon, Sautel, Schwartz, & Sokoloff, 1996;Rosenbaum et al., 2007;Roth, Hanson, & Stevens, 2008). In fact, the 3PBL crystal protein can be considered as a mutated crystal protein, because the existence of multiple point mutations in this protein is designed in order to obtain a stable D 3 R protein crystal. In addition, a series of hydrophilic residue groups designed with T4-lysozyme for feasible crystallization is outreached to supersede most of the third cytoplasmic loop (ICL3) (Arg222-Arg318). Further its stabilization is carried out with the antagonist eticlopride, a potent D 2 R/ D 3 R antagonist introduced into the purified protein.
In our previous work, based on the 3PBL protein structure, we had studied and built an original human D 3 R complex structure with dopamine (DA) by docking technology and molecular dynamics simulation (Jin et al., 2011). Because five subtype dopamine receptors have high homology, D 3 R protein structure can be used as a representative of dopamine receptors to investigate free energy for SAT moving along the molecular channels within dopamine receptor structures. Therefore, the built D 3 R-DA complex structure reported in our previous work is exploited as a material of D 3 R protein to work in this paper (Jin et al., 2011).

SAT molecule docked with D 3 R
For docking work, we used similar methods and steps reported in literatures (Chi, Xie, Zhang, & Xu, 2015;Shi et al., 2012;Wang et al., 2011;. Using MP2/6-31G (d,p) method to optimize the molecular structure of SAT with no imaginary frequency as the optimization standard (Frisch et al., 2003), SAT was then docked into D 3 R to replace DA molecule by the Dock6 program (Lang et al., 2006).

Molecular dynamics simulation for D 3 R-SAT complex protein
For the D 3 R-SAT protein obtained by docking, we used the same steps reported in the literature NH 2 Figure 1. Molecular structure of D-Benzedrine (SAT). (Jin et al., 2011) to inlay it into the POPC membrane model. Using the Gromacs 4.5.3 software package, the MD simulation was performed with visual molecular dynamics (VMD) program for analyzing simulated results (Berendsen, Van der Spoel, & Van Drunen, 1995;Hub, de Groot, Grubmüller, & Groenhof, 2014;Humphrey, Dalke, & Schulten, 1996;Van der Spoel et al., 2005). The steps of MD simulation are similar to those in our previous work (Chi et al., 2015;Xu et al., 2009). The force field parameters of POPC originated from the reported work were rewritten to be in consistence with the format of gmx force fields (Hoff et al., 2005). The force field parameters of SAT were individually defined using gmx force fields, while the parameters of residues for D 3 R were taken from the Gromacs package. The initial structure of SAT was optimized at the level of MP2/6-31G(d,p) with no imaginary frequency. The gmx-united atom types with Mulliken atomic charges were applied to define the parameters of SAT (Daura, Mark, & Van Gunsteren, 1998;Van Gunsteren et al., 1996), where the CH, CH 2 , and CH 3 groups were used as single atoms. The parameters of SAT are available in the Supplementary Material with a proof obtained in a box of 3 × 3 × 3 nm for SAT with H 2 O molecules performed by 50 ns of MD simulation. The complex protein was put into a box (+y axis corresponding to extra-cellular orientation), where a buffer of 2.0 nm was given between protein atoms and the box edge. Using Na + ion to neutralize the system, therefore, in the system there are two Na + ions, 8207 water, 84 POPC, 279 residues, and one SAT, total 31,743 atoms. The complex protein geometry was minimized to meet a standard of 1000 kJ mol −1 nm −1 using the steepest descent method to remove bad van der Waals contacts between atoms, and then the position-restrained MD simulation was run for 200 ps. Water molecules move freely, while the LINCS algorithm was applied to restrain the equilibrium positions of solute bonds, in which at the condition of one standard atmosphere pressure, the temperature rose slowly from 0 to 310 K. Next, using 2 fs of time step, 50 ns of MD simulation was carried out and output into a file of coordinates written in every 2000 steps (4 ps) for calculating various energies and analyzing system structures. Non-bonding interactions in every 10 steps for a pair list of neighbors were calculated. Along the simulation box, periodic boundary conditions were applied with a cut-off distance (1.0 nm) both for shortrange electrostatic interactions and for Lennard-Jones non-bonding interactions. A PME procedure was used to treat the long-range electrostatic effect with default parameters. The Berendsen thermostat for 310 K temperature was coupled with the time constant .1 ps. The Berendsen exponential relaxation pressures were coupled with the time constant of .5 ps.

Free energy surfaces for SAT moving along the trajectories within D 3 R
Based on the simulated D 3 R-SAT-POPC-H 2 O system, using PMF with umbrella samplings (Hess, Kutzner, Van der Spoel, & Lindahl, 2008), the molecular dynamics simulation was performed to obtain free energy changes upon the trajectories for SAT moving along the molecular channels within D 3 R. The principle of PMF is based on the formula: ΔG = −RTlnK, in which ΔG is free energy and K is equilibrium constant obtained by molecular simulation (Bian, Zhang, Wang, & Xu, 2014;Van der Spoel et al., 2013;Zhang et al., 2014). The initial structure of trajectory comes from the file of the simulated structure, and the topology necessary for the trajectory simulation is similar to its structural simulation. Along x, z, and y axes with positive and negative directions, six trajectories were obtained by simulating 500,000 steps in 2 fs every step at the simulated temperature of 310 K. SAT molecule was set up as a moving group, on which is exerted an external force of smaller than 2000 kJ mol −1 , with less than 10 nm ns −1 moving speed. A residue in the very opposite direction of movement was selected as a reference group. After the simulation of trajectory, the appropriate umbrella samplings were selected from the data file of trajectory by a general rule of .05-nm distance. To generate the trajectories for SAT moving along the molecular channels, it is necessary that external forces are provided to make SAT move, which results in the molecular configuration and the D 3 R system to be deviated from their equilibrium. Therefore, it is necessary to perform the simulations of umbrella samplings, i.e. the distances between the moving and reference groups are limited to their original mass central distance, and then perform MD simulations to let the molecular configuration and the D 3 R system rebalanced. For the simulation of umbrella samplings, there are three files including the selected umbrella sampling files, the topology file and a new pointing parameter file. The parameters for simulations of umbrella samplings are almost the same as those of trajectory simulations, just only the velocity set to 0 nm ns −1 indicating that SAT molecule do not move along the trajectory, but is limited within the given original mass central distance of two groups. The simulations were set to 400,000 steps (800 ps) in 2 fs every step for umbrella samplings rebalanced, which is determined by root mean square deviation (RMSD). After the simulations of umbrella samplings were completed, using Weighted Histogram Analysis Method (WHAM) (Hub, de Groot, & Van der Spoel, 2010) of the Gromacs 4.5 programs, the biased sampling results were converted to unbiased sampling of statistical results, and the PMF is drawn out from a series of umbrella samplings of statistical results. PMF surface is the free energy surface, on the basis of which the changes of binding free energy (ΔG bind ) can be calculated. The other files output from the g_wham calculations are the data files of histograms for umbrella samplings. The results of umbrella histograms show the overlapping degree of umbrella samplings, and the better overlapping implies the reliability of ΔG obtained by the calculation through WHAM.

SAT molecule docked with D 3 R
In Figure 2 is the docking result. The main and top views both display that SAT molecule is docked into the cavity surrounded by seven transmembrane helices near middle but close to the extra-cellular position. This position for SAT is similar as dopamine is docked within D 3 R (Jin et al., 2011). The docking energy is −29.1 kJ mol −1 for SAT molecule docked with D 3 R. If the positive-ion state of SAT is used to dock with D 3 R, its docking energy is 1.04 × 10 4 kJ mol −1 , a positive energy to mean a repulsive energy, which indicates that molecular channels are not suitable for ion-state molecule, they are molecular channels, but not ion molecule channels. Therefore, in this paper, the neutral molecular state of SAT is used, while an ionic state of SAT is used in the system without protein to permeate through the bilayer membrane (Wang, 2015). Figure 3 is a plot of RMSD vs. 50 ns MD simulation for the D 3 R-SAT complex protein within the bilayer POPC-H 2 O membrane. The values of RMSD demonstrate that after 25 ns of MD simulation, the D 3 R-SAT structure reaches a balance of stability with the background values of ±.01 nm located at .46 nm for D 3 R and 4.66 nm for the whole system. On the view of molecular dynamics simulation, the values of RMSD show that the system achieves a balance that can be used as the initial structure to study the trajectories of SAT moving within D 3 R.

Molecular dynamics simulation of D 3 R-SAT complex protein
3.3. Free energy potential surfaces for SAT moving along the trajectories within D 3 R 3.3.1. Free energy potential surfaces for SAT moving along the trajectories of functional molecular channels within D 3 R For exogenous molecules, it is well known that they usually depend on their permeability through cell membrane arriving inside cell, and under normal conditions, any molecules with a strong ability to penetrate through the cell membrane would have the potentials as medicines.
SAT has a strong ability to penetrate through cell membranes (Wang, 2015) corresponding to its functional role in clinic. On the basis of the refined D 3 R complex protein with SAT, we have studied six trajectories for SAT moving along molecular channels within D 3 R shown in Figure 4. The +y axis is supposed for SAT moving from the interior of D 3 R toward the extra-cellular orientation, while the −y axis is designed for SAT moving from the interior of D 3 R toward the intra-cellular direction. The +x, −x, +z, and −z axes are considered for SAT moving from the interior of D 3 R toward the middle space of bilayer membranes of cell. Considering four directions, it may cover the most likely path for SAT entering the space of bilayer membrane. Along the +y axis, Leu129 is selected as a reference group for SAT moving. The selection of the reference group is upon the reverse direction as near as possible close to SAT. Because water and phospholipid molecules have relatively large movability in the system, they are not suitable to be used as a reference. Along the −y axis, Ser242 is selected as the reference group for SAT. Their trajectories are shown in Figure 4, which are generated by a suitable external force of less than 2000 kJ mol −1 allotted from the simulation program. Along the +y axis, 66 samplings were selected from the trajectory data file by a rule of .05-nm mass central distance and used to simulate the free energy of SAT. In the trajectory file, there are 1001 trajectory points preserved along each direction, without necessity for all used to perform simulation. The 1001 trajectory data points correspond to 1000 ps time of MD simulation, i.e. every 1 ps keeping a trajectory point of system conformation data. In Supplementary Tables S1 and S2, besides the sampling numbers indicating the scope to calculate PMF, they also mean the simulation time length (in ps unit). Using the trajectory point systems, we calculate the mass central distances between the reference and moving groups for umbrella sampling. Due to the change of external force in different regions of system, SAT has not only rotational motion but also translational motion, leading to the mass central distance changes to be larger, so that even if two neighboring samplings, the changing value of their mass central distance is larger than .05 nm. It is a general rule using .5 nm to select umbrella samplings listed in Supplementary Tables S1 and S2.
Along the −y axis, 48 samplings were selected. These samplings, limited to their original mass central distance, are simulated to achieve both SAT molecule and D 3 R structural system rebalanced, which are determined by RMSD. For each axis, we chose four samplings as the representative to analyze the values of RMSD. The four representative sampling systems are shown in Figure 4 and labeled in o, a, b, c four trajectory points, also listed in Supplementary Tables S1 and S2 in bold type of letter. In Supplementary Figure S1, along the +y and −y axes, from 700 to 800 ps, the RMSD values of four representative samplings are maintained at the background values of ±.01 nm. The values of RMSD demonstrate that after 400,000 steps of MD simulation, the sampling systems have been rebalanced. For other samplings, the values of RMSD are similar, but their figures and data are not given in this paper. On the basis of rebalanced sampling systems, using WHAM, the biased sampling results are converted to unbiased sampling of statistical results, and PMF is drawn out from the statistical results. Although many umbrella samplings are picked out, they still belong to the biased sampling results. In Figure 5 are the results of PMF and the corresponding weighted histograms of umbrella samplings showing that the sampling windows with good overlapping to account for the reliability of PMF. The VMD program is used to display and analyze the SAT trajectories, as well as to identify and distinguish the various positions of SAT situating within D 3 R. For the sake of analysis, the images of four key positions for SAT moving along each channel are shown in Figure 5. Their labels in bold line and the distances of center mass are listed in Supplementary Tables S1 and S2. The distances of center mass shown in Figure 5 gained by WHAM have some deviations from the values calculated by the GROMACS tool. For example, along the +y axis, the minimum distance of .86 nm is obtained by WHAM in Figure 5, corresponding to the minimum distance of 1.02 nm from sampling 2 calculated by the GROMACS tool. In the trajectory of Figure 4, along the +y axis, SAT starts from the O point (.86 nm shown in Figure 5), through the a point (2.43 nm) and then passing through the b point (3.29 nm) arriving at the c point (3.56 nm), where it leaves out of the D 3 R interior structure. The result shown in Figure 5 gives that the free energy change is 83.5 kJ mol −1 for SAT moving along the +y axis. Along the −y axis in the trajectory of 3.3.2. Free energy potential surfaces for SAT moving along the trajectories of protective molecular channel within D 3 R For the free energy surfaces for SAT moving along protective molecular channels within D 3 R, we have studied four trajectories along the +x, −x, +z, and −z axes all possible for SAT transferring from the interior structure of D 3 R into the space of bilayer membranes. Figure 6 is the trajectory of SAT moving along the +x axis into the space of bilayer membrane, where Trp24 acts as a reference group. Figure 7 is the trajectory of SAT moving along the −x axis into the space of bilayer membrane, where Ala137 acts as a reference group. Figure 8 is the trajectory of SAT moving along the +z axis into the space of bilayer membrane, where Cys172 acts as a reference group. Figure 9 is the trajectory of SAT moving along the −z axis into the space of bilayer membrane, where Val56 acts as a reference group.
From the trajectory along the +x axis, 64 samplings are picked up and listed in Supplementary Table S3 up and listed in Supplementary Table S4. From the trajectory along the +z axis, 67 samplings are chosen and listed in Supplementary Table S5. From the trajectory along the −z axis, 55 samplings are picked up and listed in Supplementary Table S6. For all the selected samplings, the distances between the SAT moving and reference groups are limited to their original mass central distance, and using MD simulation they are rebalanced, which are characterized by the system's RMSDs from four representative samplings. The respective samplings are shown in Figures  6-9 and labeled in o, a, Figures  S2-S5, from 700 to 800 ps, all the RMSD values from respective samplings are kept at the background values within the vibrating value of ±.01 nm to demonstrate after through 400,000 steps of MD simulation the sampling systems all are rebalanced. For other samplings, the similar results are not given in the paper. Also shown in Figures 6-9, respectively, are PMFs drawn out from the statistical results using WHAM, and the corresponding weighted histograms of umbrella samplings to respective SAT reaction coordinates with good overlapping of umbrella samplings to show the reliability of free energy.
In Figure 6, SAT moves along the +x axis, from the o point (.60 nm), through the a point (1.12 nm), and then the b point (1.72 nm) arriving at the c point (2.57 nm), where it seems to be left out of the D 3 R interior structure. On the views of free energy curve with the main and top views for SAT positioning on the system, from the gap between TM2 and the TM3, SAT leaves out of D 3 R. It can be judged that at the b point (1.72 nm) SAT actually and basically breaks away from the control of D 3 R with ΔG of 87.7 kJ mol −1 . From the b point to the c point, actually SAT is moving within the POPC bilayer membrane, during which ΔG is 33.1 kJ mol −1 in agreement with 30.3 kJ mol −1 obtained by MD simulation for SAT in ion molecule state moving within pure POPC bilayer membrane (Wang, 2015).
In Figure  and top views for SAT positioning in the system, from the gap between TM6 and the TM7 SAT leaves out of D 3 R, and furthermore it can be judged that at the a point (2.34 nm) SAT basically breaks away from the control of D 3 R, during which ΔG is 309.4 kJ mol −1 . From the b point to the c point, SAT is moving within the POPC bilayer membrane, in which its ΔG is 42.0 kJ mol −1 to agree with the result obtained for SAT in ion state moving within POPC bilayer membrane. In Figure 8, SAT moves along the +z axis, from the o point (.42 nm), through the a point (1.62 nm), and then the b point (2.57 nm) arriving at the c point (3.48 nm), where it completely leaves out of the D 3 R interior structure. On the views of free energy curve with the main and top views for SAT positioning in the system, from the gap between TM1 and the TM2 SAT leaves out of D 3 R. From the view of free energy curve, the free energy is rising in uniform, until the b point, the free energy curve reaches a platform with the value of 224.3 kJ mol −1 . On the main view for SAT in the system, the b point just is ready for SAT to leave out of D 3 R from the gap between TM1 and TM2. After that, from the point of b to reach the c point it is actually to move within the POPC bilayer space with the free energy of 29 kJ mol −1 , also in agreement with 30.3 kJ mol −1 obtained by MD simulation for SAT within POPC bilayer membrane.
In Figure 9, SAT moves along the −z axis, from the o point (.63 nm), through the a point (1.93 nm), and then the b point (2.36 nm) arriving at the c point (3.49 nm) where it completely leaves out of D 3 R. On the main and top views for SAT in the system, from the gap between TM3 and the TM5, SAT leaves out of D 3 R. On the views of free energy curve for SAT in the system, the free energy is basically rising in uniform, until the c point, the free energy curve reaches a platform with ΔG of 221.7 kJ mol −1 .
From the trajectory pictures, it is seen that SAT moves along the molecular channels within D 3 R with a certain degree of flexibility. SAT does not keep moving along a single axis direction, but starts to move along an axis as the initial direction, and transfers always along the direction of the minimum energy, which are verified by four trajectory figures, if necessary, sometimes it performs slanted movements. When SAT is blocked by front of a transmembrane helix, it will slightly change the direction of movement, then from the transmembrane helix column gap move toward the lowest energy direction. If SAT moves along the Y axis, the most likely direction for SAT moving is not to transfer through the column gap of transmembrane helices. Now we have studied four directions for SAT moving along +x, −x, +z, −z axes, it should cover the most probable trajectory for SAT to transfer into the space of POPC bilayer membrane. For SAT moving along the protective molecular channels within D 3 R, the most probable trajectory of channel is to leave through the gap between TM2 and TM3 from D 3 R into the space of POPC bilayer membrane with ΔG of 87.7 kJ mol −1 . By MD simulation, we can even see SAT moving along different trajectories through the different TM gaps within D 3 R. If using the experimental method, we cannot peer the real-time trajectories of SAT moving along the molecular channels within D 3 R. Furthermore, the changes of free energy obtained by MD simulation for SAT moving along the molecular channels within D 3 R should be similar to those determined with the experimental method, although there is no direct evidence to support it, we can compare the related experimental data to acquire part of support for our results.

Conclusions
The change of free energy is 83.5 kJ mol −1 for SAT moving along the functional molecular channel within D 3 R toward the extra-cellular direction of trajectory, and it is 261.8 kJ mol −1 for SAT toward the intra-cellular direction of trajectory. The changes of free energy are 87.7, 309.4, 224.3, and 221.7 kJ mol −1 , respectively, for SAT moving along the different protective molecular channels within D 3 R into the space of cell bilayer membrane. It is likely that through the column gap between TM2 and TM3, SAT moves with ΔG of 87.7 kJ mol −1 along the protective molecular channel into the space of POPC bilayer membrane. Our previous work presented that the free energy for Levo-Benzedrine (RAT) to transmit toward the outside of cell along the functional molecular channel within D 3 R is 91.4 kJ mol −1 , while it is 117.7 kJ mol −1 for RAT to transmit into the space of POPC bilayer membrane along the protective molecular channel within D 3 R (Xie, Xu, Wang, & Xu, 2016). The values of free energy imply that SAT relatively easily passes through the functional molecular channel within D 3 R for exerting its molecular functions and to increase the release of functional dopamine molecules resulting in a variety of functional effects from SAT. Therefore, the obtained results show that the pharmacology and addiction mechanisms of SAT as a medicine are closely related to the molecular dynamics and mechanism for SAT to transmit along molecular channels within dopamine receptors.

Disclosure statement
No potential conflict of interest was reported by the authors.