Entropy and Polarity Control the Partition and Transportation of Drug-like Molecules in Biological Membrane

Partition and transportation of drug in the plasma membrane of a mammalian cell are the prerequisite for its function on target protein. Therefore, comprehensive understanding of the physicochemical properties and mechanism behind these complex phenomena is crucial in pharmaceutical research. By using the state-of-art molecular simulations with polarization effect implicitly or explicitly included, we studied the permeation behavior of 2-aminoethoxydiphenyl borate (2-APB), a broad-spectrum modulator for a number of membrane proteins. We showed that the protonation state and therefore the polarity of the drug is critical for its partition, and that the drug is likely to switch between different protonation states along its permeation pathway. By changing the degrees of freedom, protonation further affects the thermodynamic of the permeation pathway of 2-APB, leading to different entropic contributions. A survey on 54 analog structures with similar backbone to 2-APB showed that delicate balance between entropy and polarity plays an important role in drugs’ potency.

. The initially assembled system. The yellow represents the neutral 2-APB, and the green denotes the charged one.  Systems are assumed to be symmetric, one half of the profile is calculated by averaging the two halves, the other half is the mirror of it. By comparing with the logP determined from experiments, the 2-APB was predicted to be less lipophilic by using NPA charge. Therefore, the performance of ESP charge is much better.  Computational Methods: (1) Parameterization for 2-APB with conventional CHARMM force field. To get the initial force field parameters for 2APB, CHARMM General Force Field (CGenFF) parameters were developed in a selfconsistent fashion following previously described protocols 1-2 and initial parameters were generated from similar fragments. All empirical force field (MM) calculations were performed using the CHARMM program 3 and Quantum mechanical (QM) calculations were performed using the Gaussian09 program. 4 To prepare for QM target data, geometry optimizations and vibrational calculations were done at the MP2/6-31G(d) level. [5][6] The target QM frequencies were scaled by 0.943 to account for limitations in the level of theory. 7 QM level water interactions were used as target data for charge optimization, and water molecules in TIP3P geometry 8 were individually placed at optimized geometries interacting with hydrogen bond donors and acceptors in the model compound at its monomer MP2/6-31G(d) conformation. The interaction distance was optimized at the HF/6-31G(d) level with all other degrees of freedom fixed for each interacting pair. In order to yield parameters appropriate for the condensed phase simulation, the HF/6-31G* interaction energies were scaled by a factor of 1.16 and a difference of -0.2 Å in interaction distances between QM and MM was targeted for moderate to strong model compound-water interactions. 9-10 Relaxed potential energy scans (PES) were used to optimize dihedral parameters and were performed at the MP2/6-31G(d) level with the target dihedral angle(s) constrained at an increment step of 10°, while other degrees of freedom were fully optimized.
Gas phase MM energy minimizations were conducted using the conjugate gradient minimization method followed with the Newton-Raphson minimization algorithm to a gradient tolerance of 10 -5 kcal/mol/Å using an infinite nonbonded cutoff distance. Vibrational calculations were done using the VIBRAN module and corresponding potential energy distribution (PED) analyses were carried out using MOLVIB module with internal valence coordinates as proposed by Pulay et al. 11 Such PED analyses were used to optimize force constants in bonded parameters. MM PES were performed by reading the QM geometries of all the scan points into CHARMM and harmonic restrains with force constant of 10,000 kcal/mol/radian were placed on the target dihedral angle(s) while all other degrees of freedom were minimized.
The parametrization started from optimizing atomic charges using QM water interactions and then move forward to reproduce minimized geometry and vibration analysis, and finally trying to reproduce QM PES profile in a self-consistent way. Due to lack of related experimental data, the van der Waals (VDW) parameters (epsilon and Rmin) of the B atom were not specifically adjusted against the bulk properties of Bcontaining compounds. Fortunately in this special case, the B atom is buried at the center of bulky 2APB molecule. Even a wild range of change of the VDW parameters of B showed very limited affection on other properties like molecule-water interactions according to our test work.
(2) Explicit polarizable model. In this model, the general idea is that, the partial charges of a molecule were not fixed but were changeable in response to the local condition represented with background point charges. [12][13][14] As shown in Scheme 1, the procedure of the explicit-polarizable model is as follows: (i) Run MD simulations with the conventional force field, and relax the system; (ii) Based on the relaxed structure of the system, calculate the partial charge of a 2-APB molecule with QM methods, and the rest of the systems is represented as background point charge to mimic the polarization effect of the environment. For systems with multiple 2-APB molecules, each of them is treated with QM with all of the others being the environment. Considering the computational cost and accuracy, two cycles were performed to get the charge information of the each 2-APB molecules; (iii) Update the partial charge of all 2-APB molecules, calculate the total energy with the new parameters, and compare it with the value obtained from the former parameters. If the energy difference is within 10%, the systems is believed to be stable. Otherwise, the former parameters will be used.
(iv) With a given interval, MD simulations are carried out with the parameters determined from the above step. Such cycle is repeated until the desired total time length of simulations is achieved.
The above procedure was implemented using a homemade Fortran code and TCL script embedded in the VMD. 15 ab initio quantum mechanical (QM) was used to get new partial charges. All the QM calculations were performed using the Gaussian 09 program package 4 at the M06-2X 16 /6-311++G(d,p) level of theory. Seemingly, two factors are critical in this explicit polarizable model. One is the partial charge, and the other is the interval used for updating charges. Extensive exams on these two factors were carried out.
Two kinds of population analysis, the natural population analysis (NPA) and the electro-static-potential (ESP) were tested. And seven independent MD simulations using different interval time to renew the atomic charge, i.e. 200ps, 400ps, 600ps, 800ps, 1000ps, 2ns, and 10ns, were carried out. The result was shown in Figures S2 and S3. Consequently, the combination of ESP charge with the interval of 10 ns was used.
(3) Implicit polarizable model. Different from the explicit polarizable model, in the present implicit one, the partial charge of the drug is fixed after parameterization. During the parameterization, the polarizable continuum model (PCM) was used to take the surrounding environment into consideration, and partial charge of the drug was derived from restrained fit of electrostatic potential 17 calculated with QM calculations (RESP charge), and the AMBER 16 package suite 18 was used. The protocol used to derive ESP partial charge is same as that we mentioned above. In addition, different values of ε were tested in this S8 work, as shown in Figure S4. We found that charges are nearly independent of the dielectric constant. Therefore, the ε value of water (78.4) was used to in this model. As shown in Scheme 2, the procedure contains the following steps: (i) Run MD simulations with the conventional force field, and relax the system; (ii) Based on the relaxed structure of the system, calculate the partial charge of a 2-APB molecule with QM methods for the subsequent RESP calculations, and the PCM was used to take the environment into account. Importantly, ensemble average of the trajectory was used to account for the redistribution of molecular dipole.
(iii) Use the partial charges obtained from step (ii) to run MD simulations.
(4) System setup. The lipid bilayer containing 91 POPC molecules was set up using the Package psfgen supplied by NAMD. 19 The membrane normal was oriented parallel to the z-axis of the simulation cell. 8 charged and 8 neutral forms of drugs were distributed evenly in the system ( Figure S1), where half were in the membrane and half in the bulk water. To mimic the real physiological environment, 150 mM NaCl was added to the lipid bilayer. The total number of water, sodium, and chloride were 4281, 15, and 23 respectively, where the extra 8 anions were used to neutralize the system.

(5) Unbiased MD simulations.
The molecular dynamics (MD) simulation was performed with the NAMD package suite. 19 CHARMM 36 force field 20 was used for lipids and ions and the TIP3P model for water. 8 The system was performed in the NPT ensemble at 310 K with 1 atm pressure. The electrostatic calculations were handled by the particle-mesh Ewald (PME) 21 with a 12 Å cutoff. Integration time step was set to be 2 fs with all bonds containing hydrogen held rigid using SHAKE algorithm 22 . Langevin dynamics was used to maintain temperature, while the pressure of the system was kept by a Nose-Hoover-Langevin piston, 23 except for the hydrogen atoms. During the simulations, the temperature was rescaled every 1000 integration steps, and the barostat oscillation time scale was set to be 100 fs.
For the system by using explicit polarizable model, 4 ns/day could be achieved by using NVIDIA GTX980Ti + Interl(R) Core(TM) i7-4790K(@4.0GHz). A toal of 200 ns trajectory was accumulated. For the implicit one, the performance is 12 ns/day by using the same machine, and we ran the MD simulations for 1 μs.
(6) Free energy calculations. Due to the relative long trajectory needed to get converged free energy profile, a system containing 59 POPC molecules, 2810 TIP3P water molecules was constructed for the biased MD simulation. In the system with neutral 2-APB, 7 sodium and 7 chlorides were added to mimic the physiological environment, while in the charged one, 7 sodium and 8 chlorides were added. The same force field parameters used for the unbiased MD simulations were applied for the calculation of potential of mean force. The simulation was performed with the NVT ensemble by using the adaptive biasing force method. [24][25] The distance between the center of mass of drug and all the phosphate atoms in the lipid (therefore the membrane center) was defined as a reaction coordinate. As the two leaflets are symmetric in membrane, the permeation free energy on one side was carried out. 6 non-overlapping windows ranging from 0 to 30 Å with an increase of 5 Å were set up. A bin width was set to 0.1 Å, and 200 steps for force samples was collected.

S9
For the free energy decomposition, we run additional ABF calculations at 305 K and 315 K, respectively, with the same protocol mentioned above. To get converged free energy profiles, a total simulation time of 4.8 μs were collected for each temperature. For two protonation states at three different temperatures, a total of 28.8 μs trajectory was collected.

(7) Partition coefficient.
The definition of partition coefficient is shown in the following equation: (9) Properties of 54 stuided drug-like molecules. Their partition coefficient (logP) and pKa were all collected from the website of the sci-finder 27 , the values of dipole moment were calculated using the G09 4 at the theory level of M06-2X/6-311++G(d,p), and values of Solvent Accessible Surface Area (SASA) were calculated by software VMD 15 . To classify the inherent relationship between these values, density-based clustering method 28 was used. Clustering by using K-means implemented in the scikit-learn 29 gave consistent results.