Revisiting Partition in Hydrated Bilayer Systems - Journal of Chemical

Article pubs.acs.org/JCTC

Revisiting Partition in Hydrated Bilayer Systems Joaõ T. S. Coimbra, Pedro A. Fernandes, and Maria J. Ramos* UCIBIO, REQUIMTE, Departamento de Química e Bioquímica, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, s/n, 4169-007 Porto, Portugal S Supporting Information *

ABSTRACT: We discuss potential-of-mean force convergence issues involving the methods for generating starting structures of umbrella sampling simulations for the study of the water−bilayer permeation of drug-like compounds. We provide a proof of concept of the benefits of potential-energy modification methods, in particular the flooding technique, to enhance the sampling and improve the convergence of the potential-of-mean force. In addition, we also discuss how these potential modifiers could introduce the necessary bias to induce the correct description of other relevant internal degrees of freedom for the compounds, thus providing a more accurate description of the permeation process. Effective methodological approximations are taken into consideration, more specifically the inclusion of atomic polarization, and experimental data are used in the validation of computed values whenever possible. The parametrization of the atomic point charges is always a relevant issue in biomolecular simulations. We discuss the influence of a molecule’s conformationdependent atomic point charges and their impact on the predicted potential-of-mean force.

■

INTRODUCTION To know how compounds permeate biomembranes is very important in many fields of research. From basic pharmacology to drug design, knowledge of how drugs permeate biomembranes (be it passive diffusion or active transport) and how fast this permeation is, is a crucial aspect when one wishes to develop new drug candidates. The importance mainly resides on the relation with ADMET properties, since many drug candidates fail if they do not present good bioavailability and pharmacokinetics properties.1 Numerous approximations and empirical methods employed in “predictive” chemistry (here, we refer mainly to quantitative structure-permeability relationship (QSPR) models) have provided a way of attaining computational speed with relatively good accuracy.2 These are often employed in industrial applications. However, they entail a loss of atomistic detail, and often they consider a nonpathway dependent approach as they frequently contemplate only discrete phases or single conformations for the compounds. For biological membranes, and due to their heterogeneous nature, more dynamic methods are crucial to provide answers on the chemical-physics of compounds permeating the bilayer. In addition, empirical methods could eventually neglect the serendipity of drug development, as numerous candidates are filtered with established correlations between physical properties and successful drug development analysis (as for instance Lipinski’s rule-of-five,3 one of the most commonly employed drug-like character predictors). Physics-based methods could grant a more directed approach for drug design, despite the computing time drawbacks inherent to them. Obviously, physics-based © XXXX American Chemical Society

methods, such as molecular dynamics (MD), cannot address the level of complexity that a full biological system entails, or that of empirical models, when they are constructed from more complex systems (caco2-cells, PAMPA assays, etc.). In addition, the determination of relevant experimental quantities (such as permeability coefficients), without enough sampling, is sometimes difficult to achieve.4,5 However, the combination of both methods could eventually optimize current protocols of screening and drug design. With the exponential growth in computational power, ranging from better hardware to more efficient algorithms, we believe computational chemistry is coming of age. Is it possible that MD simulation could take its place in predictive chemistry and biochemistry? Is it capable of substituting faster (but less accurate) approaches in virtual screening campaigns? This is why, in this paper, we were interested in studying the permeation process of compounds through the bilayer. This has already been widely described in the literature, and thus we could sense the impact of these predictions.6−15 However, the literature normally suggests that microecond time-scales are often necessary for attaining convergence in the potential-ofmean force (PMF) describing permeation. More special cases also need enhanced sampling methods such as metadynamics16 for describing multiple collective variables,6 which often require a level of computational exigence that takes time. These are also not often trivial to determine, since the larger the compounds, the larger the number of degrees of freedom (DOF). Taking Received: December 29, 2016 Published: April 7, 2017 A

DOI: 10.1021/acs.jctc.6b01258 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 1. List of Compounds Evaluated in This Studya compound

three-letter code

log Kmem

T/K

bilayer system

method

ref

amlodipine dibutyl succinate 5-phenylvaleric acid

AML DBS 5PA

ibuprofen

IBF

3.75 2.40 2.94 3.17 3.80

310.15 308.15 310.15 298.15 298.15

DMPC DMPC DMPC DOPC DOPC

ultracentrifugation nondepletion PA-SPME ultracentrifugation pH metric titration pH metric titration

30 31 30 32 32

a

A three-letter code was attributed to each compound. Experimental log Kmem data are provided, together with the experimental conditions in which they were obtained.

parameters, as is customary in studies involving lipid bilayers. These have shown a very satisfactory agreement to the values found in the available experimental literature (see SI, section 0). Molecular Dynamics Simulations of Bilayer PMFs. The simulations for the PMFs for drug permeation across the bilayer were performed with the GROMACS 5.1 software.33−36 A native GPU acceleration was employed, in which the most computationally demanding calculations were performed in the GPUs (nonbonded force calculation), whereas the rest was conducted in the CPUs (bonded forces and PME electrostatics). The compounds were parametrized with the General Amber Force Field (GAFF),37 and lipid molecules were described using the SLIPIDS FF,38,39 since it has been consistent in providing very accurate results for bilayer systems, both structurally and energetically.40 The GAFF parameters were obtained with the ANTECHAMBER program41 and were converted to the Gromacs topology through the acpype.py script.42 The TIP3P water model43 was chosen, as it is most commonly employed with the Amber family of force fields, and was also employed in the SLIPIDS’s original validation. We have reduced the cutoff value of long-range interactions (to 1.2 nm), since it has been shown that no significant effects were noted with lower cutoffs than the recommended 1.4 nm (more specifically, 1.0 for the short-range electrostatic cutoff and 0.9 for the Lennard-Jones cutoff).40 The temperature was set to 310.15 K (physiological temperature) for AML and 5PA, or 308.15 K for DBS, the same temperature at which the experimental log Kmem values were determined (see Table 1). IBF was simulated at 310.15 K. The umbrella sampling (US) simulations44 were conducted for the four compounds, presented in Table 1. After the molecule’s parametrization (protocol is described in detail below), the molecules were placed in the water slab of the previously equilibrated DMPC bilayer system. The system was again equilibrated (for an additional 20 ns) in the presence of the compound with its atomic positions restrained (at ca. 3.5 nm from the bilayer’s center of mass). A constant-force pulling simulation toward the bilayer center was subsequently performed, for a maximum of 10 ns and with a force constant ranging from 50 to 75 kJ·mol−1·nm−1. At −0.2 nm from the bilayer’s center of mass, each compound was constrained with an umbrella potential, and a flooding simulation45 (also not exceeding 10 ns) was conducted, until a rotation of the compound was detected within the bilayer. This rotation implied that the compounds (more specifically AML, 5PA, and IBF) repositioned themselves within the nonpolar region of the bilayer at −0.2 nm, orienting their polar groups closer to the other bilayer leaflet (see SI, Scheme S1, for schematics of this reorientation). The details of the flooding simulations will be discussed in the following section. After we reached this transition, the compounds were pulled in the opposite direction (toward the water slab of the system), under similar conditions

into account that time is of the essence in computational screening campaigns, we were interested in the convergence of the permeation process. Another crucial aspect in MD is the atomic charge parameters for molecular mechanics force fields. This issue is widely discussed in the literature.11,17,18 From the introduction of polarizable force fields to different charge schemes, the choice is not straightforward. One has to balance the time requirements of improved charge models (such as polarizable charge schemes, or even improved solvent charge models) and the best fitted charge model per force field family. In the case of membrane permeation, the atomic point charges are another crucial aspect, since we are dealing with two (or more) distinct phases with different dielectric constants.19 In this paper, we describe a conformational flooding-based approach to achieve the convergence of the standard binding free energy of the solutes to the lipid bilayer (ΔGbinding°) of three compounds, for which experimental data exist. We discuss the importance of the protocol for generating umbrella sampling starting structures. Finally, embarking on an already known difficult case, we apply the same flooding potential to ibuprofen, to study the carboxylic acid conformation (cis/trans). An accurate portrayal of this distribution is essential for the correct description of its permeation process. With this work, we hope to contribute even more to the field of water/ membrane drug permeation, as it is intimately connected to drug design and development, and to the comprehension of the mechanism of action of an array of compounds, such as antimicrobial peptides.

■

COMPUTATIONAL DETAILS Setting Up the Bilayer System. We have chosen the hydrated DMPC bilayer as our membrane model system. For this study, a small scale model, and representative of common experimental setups, was employed to determine the thermodynamics of the permeation of four compounds through the bilayer system. The compounds are presented in Table 1, together with experimental water/membrane partition coefficients (log Kmem)20 and the conditions in which they were determined. The bilayer was prepared with the CHARMMGUI interface.21,22 A total of 64 lipids (32 per leaflet) were packed as a solvated bilayer phase, together with a 0.15 M concentration of NaCl. A hydration of ca. 50 water molecules per lipid was considered. After standard minimization, NVT and NPT pre-equilibrations protocols, an NPT conventional MD of 100 ns was run (see details in the Supporting Information (SI), and a careful analysis of important structural data parameters was performed taking into account experimental reference values (at the simulated temperature, when possible).23−29 These results can be consulted in the SI. We have analyzed the following properties over the last 20 ns: area per lipid, bilayer’s thickness, electron density profiles, and order B


Article

Journal of Chemical Theory and Computation

PMF near the point where the atomic charges are modified (more specifically between 1.2 and 0.7 nm):

to those above. Conformations spaced at 0.1 nm were extracted from this pulling simulation (comprising the center-of-mass bilayer distance slab of [−0.2; 3.5] nma total of 38 umbrella sampling windows), and these were later used as the starting structures for the umbrella sampling simulations. For each umbrella window starting conformation, an equilibration stage of 5 ns was performed, in the NPT ensemble. After the equilibration stage, two protocols were tested for AML, 5PA, and DBS. The first protocol (no f lood) comprised an umbrella simulation of 60 ns for each window, whereas the second protocol (f lood) considered a flooding simulation of 10 ns (for each window) prior to a 50 ns umbrella sampling production simulation per window, in which the flooding potential was removed. The total simulation time per compound reached 2.47 μs. IBF was grouped differently as the protocol comprising an additional flooding simulation stage was built from the last structure of the first protocol. In addition, it only comprised a flooding simulation of 200 ps for each window, preceding the umbrella sampling production simulation. The umbrella force constant applied ranged from 1000 to 1500 kJ·mol−1·nm−2, and it was adapted according to the superposition of the umbrella windows’ histograms. Detailed technical aspects for the umbrella sampling simulations are presented in the SI (see SI, section 1).46−52 The convergence of the profiles was evaluated from comparison of the free-energy cumulative profiles, by considering increasing sets of data from the umbrella sampling simulation, from the last to the first structure. The profile from the data sets of 20−50 ns or 30−60 ns (depending on the employed protocol) was used as the reference profile for the calculation of the free energy of binding (ΔGbinding°). We have decided to discard the initial 20 ns of each sampling distribution window, since it has been explained that equilibration takes place over this period, considering a large amphiphilic compound.12 Polarization Corrections to the PMF. Recently, the community has become more aware of the relevance of polarizability phenomena and their importance to PMF prediction.11 However, the applicability of polarizable force fields still represents a considerable increment to total computation times. In this study, we have accounted for polarizability effects in the determined PMFs, in a close resemblance to the model discussed by Jämbeck and Lyubartsev.11 Our approach is resumed subsequently. The reader is directed to the detailed theory in more relevant literature as well.17 All compounds were geometry optimized with B3LYP/6311+G(2d,2p) in gas-phase. The optimized structures were then used to produce polarized atomic point charges for the compounds in two polarizable continuum solvents, mimicking water (ε = 78.4) and n-hexadecane (ε = 2.04). Partial atomic charges were computed from the fitting of HF/6-311+G(2d,2p) electrostatic potentials (ESP) using the RESP scheme.53 The ESPs were obtained at 10 concentric layers around each atom and with 2500 points per atom, as described elsewhere.11 The water polarized point-charges were employed from 3.5 to 0.8 nm of the bilayer’s COM, whereas n-hexadecane polarized charges were applied from 0.7 to −0.2 nm. The polarization correction to the PMF is added subsequently by subtracting the difference in the free energy cost of polarization for condensed phases wat and bil (ΔΔGpol (wat)→(bil)) from the

pol pol pol ΔΔG(wat) → (bil) = Wd(bil) − d(gas) − Wd(wat) − d(gas)

(1)

where Wpol d(i)−d(gas) is the dipole−dipole polarization cost. This is approximated according to the following expression: pol Wd( i) − d(gas) =

1 (μ ⃗ − μgas ⃗ )T (α −1)T (μi⃗ − μgas ⃗ ) 2 i

(2)

in which μ⃗i is the dipole moment for condensed phases wat or bil (here, water and n-hexadecane, mimicking the bilayer tails), μ⃗gas is the dipole moment for the solute in the gas-phase, and α is the dipole−dipole polarizability tensor. Gas-phase dipole moments and polarizabilities were determined at the LCBLYP/311+G(2d,2p) level of theory (at the gas-phase optimized geometry), and the polarized RESP charges were employed for the determination of the dipole moments of the solutes in condensed phases. This differs slightly from the protocol described by Jämbeck and Lyubartsev, mainly due to computational costs associated with their choice of methodology (MP2/aug-cc-pVTZ for polarizabilities and dipole moments). Nonetheless, our choice was defined in view of the work by Maekawa and Moorthi,54 in which they present a benchmark for polarizability estimates from long-range corrected DFT functionals. LC-BLYP was among the best behaving functionals for 105 medium-sized organic compounds. Polarization costs per compound are presented in the SI (see SI, Table S2). All these calculations were performed with the Gaussian 09 software55 in its default settings, unless otherwise specified. Flooding Simulations. In a flooding simulation,45,56,57 a flooding potential is added to the force field to accelerate the escape from initial energy minima without affecting the reaction pathway. Slow events are then approachable with this technique that otherwise, with conventional MD, would only be accessible with much larger time scales. One of the main advantages relies on the nonrequirement of a priori knowledge of the final state, which sometimes could be a difficult, or even impossible, challenge. The flooding potential only acts on the destabilization of an educt conformation, without affecting surrounding energy barriers or product conformations derived from the previous state. To generate the collective variables that comprise the most significant global motions of each compound, principal component analysis (PCA) was applied to the configurational space ensemble derived from the first 5 ns of an umbrella sampling simulation (for each umbrella window). The PCA technique then identifies the eigenvectors of collective motions associated with the largest variance and lower frequencies (showing the largest eigenvalues). In this work, the flooding potential and PCA analysis were solely performed taking into account the tested compounds, and using the entire system for a least-squares fit. For AML, 5PA, and DBS, only the first six eigenvalues were considered for the flooding potential simulation. This is obviously a debatable issue, and in principle more essential coordinates could be chosen, depending on the system. Usually, a small number of modes suffice to describe the majority of the total fluctuation.58 In proteins, for instance, up to 20 modes have been reported to be sufficient to describe the “essential space” that incorporates the motions governing biological function (even for large proteins).59 With this protocol, we aimed at spatially limiting the bias to regions of C


Article

Journal of Chemical Theory and Computation higher interest. For IBF, since we were interested in observing the cis/trans conformational distribution of the carboxylic group, the entire set of essential coordinates (3N − 6) was chosen. By choosing only the first six PCA modes, fewer umbrella sampling window transitions occurred for the cis/trans distribution (data not shown). Another essential parameter in flooding simulations is the respective flooding strength (Efl), which is proportional to the destabilization free energy introduced by the flooding potential. For AML, 5PA, and DBS, a 300 kJ·mol−1 constant flooding strength was employed, and for IBF a 350 kJ·mol−1 maximum constant flooding strength was defined. Flooding simulations were conducted both in the pull-simulations and umbrella sampling windows as explained above. On a final note, this protocol was defined with a priori knowledge of the distributions we wished to observe. More specifically, these included a reorientation of the compound within the bilayer and the trans/cis distribution of the carboxylic group of IBF. We advise the reader not to apply the defined protocol to every simulation, and a more specific definition of the flooding parameters, Efl and choice of essential modes of PCA, should be considered in a per system fashion. This work only provides, in our belief, a proof of concept of the potentialities of flooding simulations in respect to bilayer permeation processes and convergence of umbrella sampling simulations. Additional information is provided in the SI (see SI, section 1). Another possible issue arising from the flooding simulation is the choice of reference structure for the PCA fitting. The full system comprising the initial structure for the NPT equilibration was chosen. Since the bilayer is a highly flexible system, this could lead to structural transition artifacts.58 By simulating the system without the flooding potential, during the production umbrella sampling simulations, we should account for the relaxation of eventual artifacts that could have surfaced during the flooding simulations. By also having an experimental checkpoint, we should be aware of potential limitations arising from simulations. Dihedral PMFs. We performed four additional simulations, in this case with the AMBER 12 simulation package.60 These comprised the determination of the PMFs for the trans to cis transition of the carboxylic group in both IBF and 5PA in an explicit water environment, and considering two sets of point charges: (i) charges derived from the compound with the carboxylic group conformation in cis and (ii) from the compound with the carboxylic group conformation in trans. The scanned dihedral in both molecules is presented in the SI (see SI, section 2). The two ESPs used to calculate the charges were polarized by a continuum with ε = 78.4 to mimic water (the same procedure as described above was performed). Umbrella sampling simulations were then conducted in a consecutive protocol, in which the final structure of the previous umbrella window was used as the starting point for the next window. The full details are presented in the SI (see SI, section 1).50,51,61−63

Figure 1. Molecular representation of the compounds evaluated in this study. The compounds’ names and three-letter codes are also depicted.

carboxylic acid group, which is known to be slow and to influence the PMF of the bilayer permeation process.6 We intend to study not only the convergence of the PMFs for this set of compounds but also the applicability of an enhanced sampling method that modifies the potential energy in the generation of the same PMFs (in this case, the addition of a flooding potential). The results section is thereby organized in three sections. The first one describes the importance of a careful production of umbrella sampling starting configurations considering solute distributing along hydrated bilayers. The second analyzes the accuracy of the PMFs in terms of validation with experimental partition results. It also discusses the convergence of the PMFs. The final section examines the importance of exploration of the internal degrees of freedom and issues arising from the choice of specific, nonexplicitly polarizable atomic point charges in the force field. Compounds’ Orientation within the Bilayer. One important issue that affects the convergence of the bilayer permeation PMFs is the starting conformation for the drugs at each umbrella sampling window. Numerous approaches are used to set up these configurations. Frequently, the compounds are inserted at the desired positions (relative to the bilayer’s COM) by gradually switching the intermolecular interactions between the solute and the rest of the system, or they are pulled through the bilayer (with steered MD for instance) to generate a pathway from which structures are extracted along the pulling trajectory. The latter approach is susceptible to the pulling strength and to the starting position of the compounds. It is also prone to bilayer defects that can affect the PMFs (such as formation of water pores).4,12 Other strategies that are also more demanding would be either the simulation of a hydrated bilayer/solute system, until the desired configurations are sampled (simulations could reach microsecond time scales or higher),15 or the consecutive umbrella sampling simulations, in which the previous umbrella sampled structure is used to generate the next umbrella sampling window structure. The last protocol, although conceptually desirable, is not often employed, since the previous protocols allow all umbrella sampling window simulations to be performed simultaneously,

■

RESULTS The tested compounds are represented in Figure 1. These were chosen based on the availability of experimental data, size, and molecular properties. Both 5PA and AML are polar molecules, and DBS contains a high number of rotatable bonds, which are known to be closely connected to the bioavailability of compounds.64 IBF is a difficult case to study with MD simulations, due to the cis to trans transition occurring at the D


Article

Journal of Chemical Theory and Computation and not sequentially (as in the latter case), which is much more efficient computationally. In this paper, we employ a constant-force pulling of the compound, for generation of initial umbrella sampling conformations. Despite being prone to convergence difficulties, it is one of the least computer demanding methods. We have tried to apply the smallest possible force in the solute in order to preserve the hydrated bilayer structure, and to minimize protrusions arising from the pulling of the compounds. Subsequently, we will further explore the convergence of the PMF considering our data set. For AML, we have tested two types of pulling. The first one comprised a pulling from the water slab to the bilayer center (defined as W_to_B); the second considered the pulling from the proximity of the bilayer center (more specifically from −0.2 nm) to the water slab of the system (B_to_W). The structure at −0.2 nm was generated from the W_to_B simulation, and pre-equilibrated as described in the Computational Details section: with a flooding potential until the desired reorientation was produced. For the orientation analysis of AML, we have produced three additional umbrella sampling windows (in the range [−0.5; −0.3] nm). This was performed in order to enhance the observed differences between the employed protocols. As we can see through Figure 2, both protocols give distinct orientations of AML within the bilayer. The orientation was determined from the angle between a vector defined with two atoms of the compound (see SI, Figure S2) and the bilayer normal (the z axis). In an equilibrium state, the PMF profile should be symmetric at 0.0 nm, as we are dealing with a homogeneous bilayer system, composed solely of DMPC phospholipids. If PMF symmetry is also dependent on the drugs’ orientation within the bilayer (near 0.0 nm there should be an inversion of the compound’s orientation relative to the bilayers’ normal), we see that for the W_to_B umbrella sampling, the AML orientation (and thus the PMF) is not symmetric in 0.0 nm (see Figure 2). This is supported also by the PMF shown in the SI (see SI, Figure S3), where we can see that AML has an inversion of the PMF at ca. 0.0 nm considering the B_to_W protocol and continues to increase in the W_to_B umbrella sampling protocol (until ca. −0.3 nm). For the other compounds, we have not tested the W_to_B pulling (as it has already been proven to be a slower converging approach).12 The orientation for the other compounds is shown in the SI (see SI, Figures S4 and S5). The orientation for 5PA and IBF is also inverted near 0.0 nm. For DBS, since this is a highly symmetric compound, no distinct orientation is detected, as the compound is shown to indistinctively rotate within the bilayer’s tail slab. Still, some discrepancies could be detected in the case of AML. In Figure 2, we see in the middle graph (plot B) that in the region of 0.2 nm (to the bilayer COM), there is a window in which the AML orientation is inverted. As discussed in the methodology section, a 10 ns umbrella sampling simulation with an additional conformational flooding potential was conducted for each window. In this case, shown in the bottom graph in Figure 2 (plot C), we see a more coherent distribution of the orientation of AML within the bilayer. The hysteresis of the permeation process has already been discussed for amphiphilic compounds, comparing both types of pulling discussed here: starting from the bilayer and water slabs, respectively.12 However, our studies in the AML compound support that this issue is also relevant for smaller compounds,

Figure 2. Amlodipine orientation relative to the bilayer’s normal axis (z axis) for each window. This analysis comprised the last 30 ns of each window. (A) Umbrella sampling simulation for amlodipine’s W_to_B protocol. (B) Umbrella sampling simulation considering the B_to_W protocol. (C) Umbrella sampling simulation considering the B_to_W protocol and comprising an additional 10 ns constant flooding simulation. On the left of each plot, we present a snapshot65 of the bilayer’s z axis projection (from −0.5* to 3.5 nm of the bilayers’ center-of-mass distance). See text for details on the different protocols. *The analysis of the orientation for AML was extended until −0.5 nm, so as to clearly observe the reorientation of the vector.

which should still exhibit high autocorrelation times for rotation throughout the bilayer region,10 and that problem can be solved or ameliorated through the use of enhanced sampling methods. Approaching Experiment and Convergence. In an attempt to evaluate the effects of a flooding potential on hidden E


Article


In terms of binding free energy, we observe no relevant differences, in the case of 5PA, or an improvement toward experimental determined values, for AML and DBS, when comparing the f lood and no f lood protocols. Specifically in terms of convergence, simulations where the conformational flooding potential is employed often exhibit faster convergence of the resulting PMFs. This was observed from the analysis of PMFs, considering increasing time sets from the last structure (see SI, section 3). The differences in ΔGbind° between protocols and the experimental data are not substantial for these compounds, despite the better description of the experimental reference data by the simulations subjected to a flooding stage. This approach might work better for more complex cases (showing increased degrees of freedom and slower relaxation rates), and as stated by others, the applicability to peptides should be interesting in future studies.6 Sampling Difficult Degrees of Freedom. Finally, we have looked at IBF, which represents a challenge in terms of PMF prediction. In Figure 4, we depict the carboxylic acid dihedral angle distribution for the umbrella sampling simulations conducted. We see that, without resorting to the conformational flooding potential, the dihedral angle remains mainly in the cis conformer in the majority of umbrella sampling windows and throughout the last 30 ns of each window (see Figure 4, plot A). When we apply a flooding potential considering 3N − 6 PCA eigenvectors, we observe a transition of the carboxylic acid conformer dependent on the position within the bilayer. In the majority of the umbrella sampling windows comprising the polar region, we have detected a transition to the trans conformer, whether in the nonpolar region the conformer remained in the cis conformer. This behavior has already been observed by others.6 This is related to the contrast interactions through hydrogen bond networking with water or 1,4 intramolecular interactions with the carbonyl. For some umbrella sampling windows, which comprised the polar region of the bilayer, this transition did not occur. This is a phenomenon that we attribute to metastable states. However, when we extend the simulations for umbrella sampling analysis (and without the flooding bias) starting from the conformations generated with the flooding potential, we have observed a reversion to the cis conformer for many of these umbrella windows (see Figure 4, plot B). In addition, the convergence (see SI, section 3) and binding free energy (see Table 2) were poorly predicted for this PMF. This was an unexpected result, taking into account previous studies.6 As a result, we have decided to look for the root of the problem. First, we performed a geometry optimization of the IBF structure, with the carboxylic acid in its trans conformer. We carried on by producing polarized charges on this conformer (in the same fashion as for cis). We then performed PMFs for the complete rotation of its carboxylic group dihedral angle (see SI, Figure S8). These PMFs were produced for both charge distributions and comprising a system of the compounds inserted in a TIP3P water box. For comparison purposes, we have performed the same calculations for the 5PA compound, since it also has a carboxylic group, but its dihedral angle distribution during the umbrella sampling simulation showed an expected behavior (see SI, Figure S7). Very interestingly, the PMFs were very similar for both charge distributions in 5PA, whether for IBF the cis and trans PMFs were distinct. In particular, the trans charges stabilized in ca. 0.6 kcal·mol−1

convergence problems of the PMFs, for three of the compoundsAML, 5PA, and DBSwe have tested the ability to correctly describe the experimental partition on the same experimental model and under the same pressure/temperature conditions. We have considered a protocol in which the polarization is taken into account in a very simplistic (but effective) way.11 This protocol considers the inclusion of polarization effects without resorting to polarizable force fields, which are still computationally demanding for large biological systems. Two approaches were considered: (i) during the first 10 ns umbrella sampling simulation, a conformational flooding potential was applied to the compounds (comprising the first six PCA eigenvalues generated from a 5 ns simulation) in each windowthe f lood protocoland (ii) under the same conditions of the previous simulations but without resorting to the conformational flooding potentialthe no f lood protocol. To evaluate the PMFs’ accuracy, we have determined the binding free energy (ΔGbinding°), which can be determined from experimental partition coefficients. Computationally, it can be calculated according to eq 3: ⎛ 1 ° = −kBT ln⎜ ΔGbind ⎝ 2zb

zb

∫−z

b

⎞ e−βw(z) dz⎟ ⎠

(3)

where zb = 3.5 nm and represents the absolute distance at which w(z) = 0.0 kcal·mol−1 and β = 1/kBT, where T is adjusted to the temperature in which the PMFs were produced and kB represents the Boltzmann constant. Partition coefficients can also be determined according to eq 4: ° P = e−β ΔG bind

(4)

In Table 2, we present experimentally determined binding free energies (ΔGbinding°), together with computed estimates. In Table 2. Experimental (exptl.) and Computed (sim.) ΔGbinding° and P Values for the Four Compounds Depicted in This Studya ΔGbinding°/kcal·mol−1 compound

protocol

AML

flood no flood WtoB flood no flood flood no flood flood (cis) flood (trans) no flood (cis)

DBS 5PA IBF

sim. −5.50 −4.73 −4.43 −3.65 −4.42 −3.14 −3.16 −6.93 −5.47 −6.19

± ± ± ± ± ± ± ± ± ±

0.59 0.43 0.50 0.42 0.40 0.39 0.47 0.73 0.42 0.44

P

exptl.

sim.

exptl.

−5.32

7556 2158 1321 387 1353 164 167 76376 7125 23186

5623

−3.38 −4.71

251 871

a

The different protocols in which these values were computed are also portrayed, and details are provided in the main text. Simulation contemplating derived charges taking into account different isomers of IBF (cis or trans) are also highlighted.

Figure 3, we also depict the PMFs from which the data were generated. We obtain a very good agreement with experimental findings, with a deviation of circa 0.3 kcal·mol−1 for AML and DBS, and around 1.5 kcal·mol−1 for 5PA. This is not far from 1kBT, which should be expected from molecular mechanics based methods. The larger differences observed for 5PA are probably related to force field parameters and not on the methodology (as both protocols produce the same results). F


Article


Figure 3. PMFs comprising the four tested compounds. Different simulation protocols are depicted and are discussed in the text. The cis or trans are related to the charge distributions generated from the cis or trans conformation of the carboxylic acid group of IBF. The compounds’ name is presented on the top right corner above each plot. In the AML plot, we have introduced a molecular representation of the bilayer system in the background,65 in order to visually support the distance to the bilayer COM.

(∼1kBT) the trans conformer, while for the cis charges the difference was close to null. This result could eventually justify the improved distribution of the dihedral angles observed for the umbrella sampling simulation conducted with the same trans charges (see Figure 4, plot C). We see that the population of the cis conformer is reduced in the polar regions of the hydrated bilayer system, a clear improvement over the cis charge’s f lood and no f lood simulations. These results support that a careful inspection of the atomic point charges is always necessary, to account for the correct conformational sampling within distinct phases of a hydrated bilayer system. In principle, it should be easily corrected through the inclusion of polarization within the force field, describing the system. However, this should be computationally demanding, albeit probably providing the most accurate

estimation of the PMF for the permeation of solutes through hydrated bilayers.

■

CONCLUSIONS We have presented a study of the influence of a biasing potential and of the atomic point charges over the convergence of the PMFs for drug permeation across water−bilayer permeation. We have applied a conformational flooding potential prior to an umbrella sampling simulation for the permeation for three drugs (amlodipine, dibutyl succinate, and 5-phenylvaleric acid), to accelerate the convergence of the PMFs. The results with a flooding potential reproduced slightly better the experimental partition coefficients. The same flooding potential was applied to ibuprofen, a compound known for difficulties inherent to the inefficient G


Article


umbrella sampling simulation, starting from the same structure as the above case). This highlights the necessity of a careful inspection of the atomic charges. It also indicates that polarizable force fields may represent the best option for this case, as different conformations have a significant effect in atomic charge, and vice versa. Despite all technical issues, the computational methods used here reproduce well the basic molecular processes occurring during permeation. The flooding potential was able to improve the accuracy of the umbrella sampling simulated PMFs, by accelerating the conformational transitions in each umbrella window, in particular in terms of the compounds’ orientation within the bilayer. In addition, it was also able to enhance the conformational description of DOF orthogonal to the bilayer insertion reaction coordinate (in this case, a cis/trans conformation transition of a carboxylic group). Full exploration of these DOFs is crucial for the correct description of the minimum energy pathway.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b01258. Membrane’s structural validation (Table S1 and Figure S1), additional protocol specifications (Table S2), and supporting results (Scheme S1 and Figures S2−S11) (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Joaõ T. S. Coimbra: 0000-0001-9138-7498 Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. Rui Neves for his meaningful discussions of the work. We also acknowledge Ó scar Passos for providing the necessary technical assistance in terms of computing hardware. This work received financial support from the following institutions: European Union (FEDER funds POCI/01/ 0145/FEDER/007728) and National Funds (FCT/MEC, Fundaçaõ para a Ciência e Tecnologia and Ministério da Educaçaõ e Ciência) under the Partnership Agreement PT2020 UID/MULTI/04378/2013; NORTE-01-0145-FEDER-000024, supported by Norte Portugal Regional Operational Programme (NORTE 2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF). J.T.S.C. further thanks the FCT for grant SFRH/BD/87434/2012.

Figure 4. IBF carboxylic acid dihedral distribution considering three umbrella sampling simulations. (A) no f lood protocol, (B) f lood protocol with the cis charge distribution, and (C) f lood protocol with the trans charge distribution.

exploration of the conformational distribution of its carboxylic group. With the flooding potential, we were able to observe (on the picosecond scale of the “flooded” simulation time) the desired distributionbetween cis and trans conformersover the full hydrated bilayer system. Still, when the flooding potential was removed, and an umbrella sampling simulation was conducted for 50 ns/window, the desired trans conformation reverted back to cis in some windows. This led us to study the effects of conformation-adaptable atomic point charges of the carboxylic acid in the permeation of ibuprofen through the lipid bilayer. The trans to cis transition occurred less when the charges were calculated from the trans conformation (considering an additional 50 ns/window

■

REFERENCES

(1) Gleeson, M. P.; Hersey, A.; Montanari, D.; Overington, J. Probing the links between in vitro potency, ADMET and physicochemical parameters. Nat. Rev. Drug Discovery 2011, 10, 197−208.

H


Article

Journal of Chemical Theory and Computation (2) Leung, S. S. F.; Mijalkovic, J.; Borrelli, K.; Jacobson, M. P. Testing Physical Models of Passive Membrane Permeation. J. Chem. Inf. Model. 2012, 52, 1621−1636. (3) Lipinski, C. A. Lead- and drug-like compounds: the rule-of-five revolution. Drug Discovery Today: Technol. 2004, 1, 337−41. (4) Lee, C. T.; Comer, J.; Herndon, C.; Leung, N.; Pavlova, A.; Swift, R. V.; Tung, C.; Rowley, C. N.; Amaro, R. E.; Chipot, C.; Wang, Y.; Gumbart, J. C. Simulation-Based Approaches for Determining Membrane Permeability of Small Compounds. J. Chem. Inf. Model. 2016, 56, 721−733. (5) Swift, R. V.; Amaro, R. E. Back to the Future: Can Physical Models of Passive Membrane Permeability Help Reduce Drug Candidate Attrition and Move Us Beyond QSPR? Chem. Biol. Drug Des. 2013, 81, 61−71. (6) Jambeck, J. P. M.; Lyubartsev, A. P. Exploring the Free Energy Landscape of Solutes Embedded in Lipid Bilayers. J. Phys. Chem. Lett. 2013, 4, 1781−1787. (7) Chew, C. F.; Guy, A.; Biggin, P. C. Distribution and Dynamics of Adamantanes in a Lipid Bilayer. Biophys. J. 2008, 95, 5627−5636. (8) Bemporad, D.; Essex, J. W.; Luttmann, C. Permeation of small molecules through a lipid bilayer: A computer simulation study. J. Phys. Chem. B 2004, 108, 4875−4884. (9) MacCallum, J. L.; Tieleman, D. P. Computer simulation of the distribution of hexane in a lipid bilayer: Spatially resolved free energy, entropy, and enthalpy profiles. J. Am. Chem. Soc. 2006, 128, 125−130. (10) Neale, C.; Bennett, W. F. D.; Tieleman, D. P.; Pomes, R. Statistical Convergence of Equilibrium Properties in Simulations of Molecular Solutes Embedded in Lipid Bilayers. J. Chem. Theory Comput. 2011, 7, 4175−4188. (11) Jambeck, J. P.; Lyubartsev, A. P. Implicit inclusion of atomic polarization in modeling of partitioning between water and lipid bilayers. Phys. Chem. Chem. Phys. 2013, 15, 4677−86. (12) Filipe, H. A. L.; Moreno, M. J.; Rog, T.; Vattulainen, I.; Loura, L. M. S. How To Tackle the Issues in Free Energy Simulations of Long Amphiphiles Interacting with Lipid Membranes: Convergence and Local Membrane Deformations. J. Phys. Chem. B 2014, 118, 3572− 3581. (13) Ma, J.; Domicevica, L.; Schnell, J. R.; Biggin, P. C. Position and orientational preferences of drug-like compounds in lipid membranes: a computational and NMR approach. Phys. Chem. Chem. Phys. 2015, 17, 19766−19776. (14) Orsi, M.; Essex, J. W. Permeability of drugs and hormones through a lipid bilayer: insights from dual-resolution molecular dynamics. Soft Matter 2010, 6, 3797−3808. (15) Paloncyova, M.; Berka, K.; Otyepka, M. Convergence of Free Energy Profile of Coumarin in Lipid Bilayer. J. Chem. Theory Comput. 2012, 8, 1200−1211. (16) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−6. (17) Swope, W. C.; Horn, H. W.; Rice, J. E. Accounting for Polarization Cost When Using Fixed Charge Force Fields. II. Method and Application for Computing Effect of Polarization Cost on Free Energy of Hydration. J. Phys. Chem. B 2010, 114, 8631−8645. (18) Jambeck, J. P. M.; Mocci, F.; Lyubartsev, A. P.; Laaksonen, A. Partial Atomic Charges and Their Impact on the Free Energy of Solvation. J. Comput. Chem. 2013, 34, 187−197. (19) Tanizaki, S.; Feig, M. A generalized Born formalism for heterogeneous dielectric environments: Application to the implicit modeling of biological membranes. J. Chem. Phys. 2005, 122, 124706. (20) Endo, S.; Escher, B. I.; Goss, K. U. Capacities of Membrane Lipids to Accumulate Neutral Organic Chemicals. Environ. Sci. Technol. 2011, 45, 5912−5921. (21) Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; DavilaContreras, E. M.; Qi, Y. F.; Lee, J. M.; Monje-Galvan, V.; Venable, R. M.; Klauda, J. B.; Im, W. CHARMM-GUI Membrane Builder Toward Realistic Biological Membrane Simulations. J. Comput. Chem. 2014, 35, 1997−2004.

(22) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: a webbased graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859−65. (23) Douliez, J. P.; Leonard, A.; Dufourc, E. J. Restatement of Order Parameters in Biomembranes - Calculation of C-C Bond Order Parameters from C-D Quadrupolar Splittings. Biophys. J. 1995, 68, 1727−1739. (24) Kucerka, N.; Nieh, M. P.; Katsaras, J. Fluid phase lipid areas and bilayer thicknesses of commonly used phosphatidylcholines as a function of temperature. Biochim. Biophys. Acta, Biomembr. 2011, 1808, 2761−2771. (25) Kucerka, N.; Liu, Y. F.; Chu, N. J.; Petrache, H. I.; TristramNagle, S.; Nagle, J. F. Structure of fully hydrated fluid phase DMPC and DLPC lipid bilayers using X-ray scattering from oriented multilamellar arrays and from large unilamellar vesicles. Biophys. J. 2005, 88, 2626. (26) Petrache, H. I.; Tristram-Nagle, S.; Nagle, J. F. Fluid phase structure of EPC and DMPC bilayers. Chem. Phys. Lipids 1998, 95, 83−94. (27) Almeida, P. F. F.; Vaz, W. L. C.; Thompson, T. E. Lateral Diffusion in the Liquid-Phases of Dimyristoylphosphatidylcholine Cholesterol Lipid Bilayers - a Free-Volume Analysis. Biochemistry 1992, 31, 6739−6747. (28) Oradd, G.; Lindblom, G.; Westerman, P. W. Lateral diffusion of cholesterol and dimyristoylphosphatidylcholine in a lipid bilayer measured by pulsed field gradient NMR spectroscopy. Biophys. J. 2002, 83, 2702−2704. (29) Filippov, A.; Oradd, G.; Lindblom, G. Influence of cholesterol and water content on phospholipid lateral diffusion in bilayers. Langmuir 2003, 19, 6397−6400. (30) Austin, R. P.; Davis, A. M.; Manners, C. N. Partitioning of Ionizing Molecules between Aqueous Buffers and PhospholipidVesicles. J. Pharm. Sci. 1995, 84, 1180−1183. (31) Vaes, W. H. J.; Urrestarazu Ramos, E.; Verhaar, H. J. M.; Cramer, C. J.; Hermens, J. L. M. Understanding and estimating membrane/water partition coefficients: Approaches to derive quantitative structure property relationships. Chem. Res. Toxicol. 1998, 11, 847−854. (32) Avdeef, A.; Box, K. J.; Comer, J. E. A.; Hibbert, C.; Tam, K. Y. pH-metric logP 10. Determination of liposomal membrane-water partition coefficients of ionizable drugs. Pharm. Res. 1998, 15, 209− 215. (33) Berendsen, H. J. C.; Vanderspoel, D.; Vandrunen, R. Gromacs a Message-Passing Parallel Molecular-Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (34) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (35) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (36) Pall, S.; Abraham, M. J.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. Lect. Notes Comput. Sc. 2015, 8759, 3−27. (37) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157−1174. (38) Jambeck, J. P. M.; Lyubartsev, A. P. Derivation and Systematic Validation of a Refined All-Atom Force Field for Phosphatidylcholine Lipids. J. Phys. Chem. B 2012, 116, 3164−3179. (39) Jambeck, 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. (40) Paloncyova, M.; Fabre, G.; DeVane, R. H.; Trouillas, P.; Berka, K.; Otyepka, M. Benchmarking of Force Fields for MoleculeI


Article

Journal of Chemical Theory and Computation Membrane Interactions. J. Chem. Theory Comput. 2014, 10, 4143− 4151. (41) Wang, J. M.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (42) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 367. (43) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (44) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187−199. (45) Lange, O. E.; Schafer, L. V.; Grubmuller, H. Flooding in GROMACS: Accelerated barrier crossings in molecular dynamics. J. Comput. Chem. 2006, 27, 1693−1702. (46) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (47) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (48) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (49) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 126. (50) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (51) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. Multidimensional Free-Energy Calculations Using the Weighted Histogram Analysis Method. J. Comput. Chem. 1995, 16, 1339−1350. (52) Hub, J. S.; de Groot, B. L.; van der Spoel, D. g_wham-A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates. J. Chem. Theory Comput. 2010, 6, 3713−3720. (53) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges - the Resp Model. J. Phys. Chem. 1993, 97, 10269−10280. (54) Maekawa, S.; Moorthi, K. Polarizabilities from Long-Range Corrected DFT Calculations. J. Chem. Eng. Data 2014, 59, 3160− 3166. (55) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (56) Bouvier, B.; Grubmuller, H. A molecular dynamics study of slow base flipping in DNA using conformational flooding. Biophys. J. 2007, 93, 770−86. (57) Grubmuller, H. Predicting Slow Structural Transitions in Macromolecular Systems - Conformational Flooding. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 52, 2893− 2906.

(58) Hayward, S.; de Groot, B. L. Normal modes and essential dynamics. Molecular Modeling of Proteins, Methods in Molecular Biology; Humana Press: Totowa, NJ, 2008; Vol 443, pp 89−106. (59) David, C. C.; Jacobs, D. J. Principal component analysis: a method for determining the essential dynamics of proteins. Methods Mol. Biol. 2014, 1084, 193−226. (60) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Götz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco, 2012. (61) Pastor, R. W.; Brooks, B. R.; Szabo, A. An Analysis of the Accuracy of Langevin and Molecular-Dynamics Algorithms. Mol. Phys. 1988, 65, 1409−1419. (62) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (63) Grossfield, A. WHAM: the weighted histogram analysis method. http://membrane.urmc.rochester.edu/content/wham. (64) Veber, D. F.; Johnson, S. R.; Cheng, H. Y.; Smith, B. R.; Ward, K. W.; Kopple, K. D. Molecular properties that influence the oral bioavailability of drug candidates. J. Med. Chem. 2002, 45, 2615−2623. (65) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38.

J


Revisiting Partition in Hydrated Bilayer Systems - Journal of Chemical

Recommend Documents