Predicting the Binding Mode of 2-Hydroxypropyl-β-cyclodextrin to

Mar 13, 2018 - ‡Department of Physical Pharmaceutics, Graduate School of Pharmaceutical Sciences and §Department of Clinical Chemistry and Informat...
0 downloads 5 Views 3MB Size
Download PDF
Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Predicting the Binding Mode of 2‑Hydroxypropyl-β-cyclodextrin to Cholesterol by Means of the MD Simulation and the 3D-RISM-KH Theory Yuji Hayashino,† Masatake Sugita,*,† Hidetoshi Arima,‡ Tetsumi Irie,§ Takeshi Kikuchi,† and Fumio Hirata*,∥ †

Ritsumeikan University, College of Life Science, Department of Bioinformatics, 1-1-1, Noji-higashi, Kusatsu, Shiga 525-8577, Japan Department of Physical Pharmaceutics, Graduate School of Pharmaceutical Sciences and §Department of Clinical Chemistry and Informatics, Graduate School of Pharmaceutical Sciences, Kumamoto University, 5-1 Oe-honmachi, Chuo-ku, Kumamoto 862-0973, Japan ∥ Toyota Physical and Chemical Research Institute, 41-1, Yokomichi, Nagakute, Aichi 480-1192, Japan ‡

S Supporting Information *

ABSTRACT: It has been found that a cyclodextrin derivative, 2hydroxypropyl-β-cyclodextrin (HPβCD), has reasonable therapeutic effect on Niemann-Pick disease type C, which is caused by abnormal accumulation of unesterified cholesterol and glycolipids in the lysosomes and shortage of esterified cholesterol in other cellular compartments. We study the binding affinity and mode of HPβCD with cholesterol to elucidate the possible mechanism of HPβCD for removing cholesterol from the lysosomes. The dominant binding mode of HPβCD with cholesterol is found based on the molecular dynamics simulation and a statistical mechanics theory of liquids, or the threedimensional reference interaction site model theory with Kovalenko-Hirata closure relation. We examine the two types of complexes between HPβCD and cholesterol, namely, one-toone (1:1) and two-to-one (2:1). It is predicted that the 1:1 complex makes two or three types of stable binding mode in solution, in which the βCD ring tends to be located at the edge of the steroid skeleton. For the 2:1 complex, there are four different types of the complex conceivable, depending on the orientation between the two HPβCDs: head-to-head (HH), head-to-tail (HT), tail-to-head (TH), and tail-to-tail (TT). The HT and HH cyclodextrin dimers show higher affinity to cholesterol compared to the other dimers and to all the binding modes of 1:1 complexes. The physical reason why the HT and HH dimers have higher affinity compared to the other complexes is discussed based on the consistency with the 1:1 complex. On the one hand, in case of the HT and HH dimers, the position of each CD in the dimer along the cholesterol chain comes right on or close to one of the positions where a single CD makes a stable complex. On the other hand, one of the CD molecules is located on unstable region along the cholesterol chain, for the case of TH and TT dimers.

■

apeutic effect on NPC.8 Although the treatment of HPβCD induces some side effect,9,10 for instance, hearing defect and potential pulmonary toxicity, HPβCD has already passed the first and second phases of the clinical trial.11 A multinational, randomized, controlled trial of intrathecal HPβCD is now ongoing.11,12 However, the functional mechanism of HPβCD for removing the cholesterol is not known very well even though it is suggested that the CDs enhance the cholesterol transport between membranes.13 Efflux process of cholesterol by cyclodextrin is separated at least to two steps, binding

INTRODUCTION Niemann-Pick disease type C (NPC) is a rare disease, medical treatment of which has not been well-developed.1−4 NPC is caused by mutation of either one of two proteins, NPC1 and NPC2, that gives rise to loss or reduction of their function as a transporter of cholesterol: NPC1, a membrane protein; NPC2, a soluble protein.5,6 Mutation of NPC1 or NPC2, which transports unesterified cholesterol from lysosome to endoplasmic reticulum, causes abnormal accumulation of unesterified cholesterol and glycolipids in the late endosome/lysosomes (LE/LY) and shortage of esterified cholesterol in other cellular compartments, and it induces death of the hepatic and/or neuron cells.1,7 Recently, it has been found that a cyclodextrin derivative, 2hydroxypropyl-β-cyclodextrin (HPβCD), has reasonable ther© XXXX American Chemical Society

Special Issue: Ken A. Dill Festschrift Received: March 2, 2018 Published: March 13, 2018 A

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Chemical structure of the (a) HPβCD and (b) cholesterol. In this study, we use a structure where all second hydroxyl groups on the glucopyranose units are substituted by 2-hydroxypropyl groups.

cholesterol, and transporting it to elsewhere. It is crucial to find the binding affinity and binding mode between HPβCD and cholesterol, in order for elucidating the mechanism of the efflux process. However, it is not so easy to determine the binding mode by the experimental method. One of the factors that makes the analysis of the functional mechanism difficult is in the low solubility of cholesterol to aqueous solutions. The low solubility of cholesterol hinders experimental analyses of binding mode in the physiological condition.14,15 Therefore, it is important to develop methodologies that may be able to overcome such difficulty for predicting the binding mode, to apply HPβCD to a treatment of NPC and/or to develop CD homologues, for example, γ-cyclodextrins, that have less side effects.16 One such methodology that has potential to overcome the difficulty may be the in silico evaluation of binding mode and affinity, which has been applied for years to design drug compounds that target a protein.17−19 The method has been employed to analyze the binding mode of some cyclodextrin derivatives and cholesterol.20−22 Recently, a new method for the in silico drug design has been introduced that combines the molecular dynamics simulation (MD) with a statistical mechanics theory of liquids, or the three-dimensional reference interaction site model theory with Kovalenko−Hirata closure relation (3D-RISM-KH). The method, referred to as MM/3DRISM-KH, is attracting increasing attention for predicting the binding mode and binding affinity of the host−guest systems, such as protein and ligands.23−26 The 3D-RISM-KH theory enables us to calculate the solvent distribution around a host molecule, such as protein, and the solvation free energy with reasonable accuracy.23,27−30 The method is free from concerns about so-called “sampling” for the degrees of freedom of the solvents, because it performs the configuration integral of the solvent analytically. Nevertheless, there is some problem that may not be solved just by the method, and that is the conformational fluctuation of a solute molecule. The conformational degrees of freedom of a solute molecule must be sampled by means of another method such as the MD simulation. In the previous report, we developed a protocol for screening cyclodextrin derivatives, based on the MM/3D-RISM-KH method.31 In the protocol, we first find probable binding modes of a ligand in a cyclodextrin derivative from the MD simulation based on the umbrella sampling method, and then we calculate the binding free energy of the ligand to cyclodextrin for each probable binding mode by means of MM/3D-RISM-KH.

In the present study, we apply the protocol to predict the most probable binding mode of the cholesterol and HPβCD to elucidate the functional mechanism of the HPβCD to the NPC.

■

MATERIALS AND METHODS

Target Systems. In this study, HPβCD is regarded as a receptor, and a cholesterol is regarded as a ligand. Six-, seven-, and eight-member rings consisting of the glucopyranose units correspond to α-CD, β-CD, and γ-CD, respectively.32 HPβCD is the most popular derivative of the cyclodextrins that have hydroxypropyl groups instead of hydroxyl groups.32 The structures of host and guest molecules are shown in Figure 1. In the structure of a host molecule, in this research, all the second hydroxyl groups on the glucopyranose units are substituted by a 2-hydroxypropyl group. Concerning the structure of HPβCD, there is another complication to be taken into account. It has been found from an experiment that HPβCD seems to bind the cholesterol to make a 2:1 complex when the concentration of HPβCD is high.33,34 In this study, we examine both complexes 1:1 and 2:1, between HPβCD and cholesterol, to find which complex is more favorable. General Procedure for Identifying the Dominant Binding Mode. The protocol to find the dominant binding mode consists of three steps: finding candidates of the binding mode, sampling fluctuation around the possible binding modes, and determining a dominant binding mode. The protocol is based on the previous work that concerns prediction of binding affinity between the HPβCD and seven small molecules.31 The detailed protocol is shown in the following section. Finding Candidates of the Binding Mode. In the present study, the umbrella sampling method is applied for finding possible binding mode.35 The umbrella sampling method requires a reaction coordinate along which the probability distribution of the configuration is calculated. Since we are concerned with the two types of receptor−ligand complex, 1:1 and 1:2, as described in the previous section, we define the reaction coordinate separately for each complex. For the 1:1 complex, we define a reaction coordinate consistent with the previous work, that is, the distance R between the center of mass of the oxygen atoms consisting of glycoside bonds in a cyclodextrin and the 3rd or 25th carbon atom on the cholesterol molecule as shown in Figure 1b. There are two possible mutual orientations for inserting a cholesterol molecule into a cyclodextrin, type I and type II, which are illustrated in Figure 2. In the case of type I, the dimethyl B

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Two types of the 1:1 binding mode of HPβCD and cholesterol. (a) Type I binding mode: the carbon atom labeled as “25” in Figure 1 is located in the substituted side of HPβCD. The reaction coordinate is defined as the distance between the carbon atom and the center of mass of the ring consisting of glycosidic-oxygen atoms. (b) Type II binding mode: the carbon atom labeled as “3” in Figure 1 is located in the substituted side of HPβCD. The reaction coordinate is defined as the distance between the carbon atom and the center of mass of the ring consisting of glycosidicoxygen atoms.

Figure 3. Four variations of the 2:1 binding modes of the HPβCD and cholesterol. These modes are termed as (a) head-to-head (HH), (b) head-totail (HT), (c) tail-to-head (TH), and (d) tail-to-tail (TT).

terminus of cholesterol comes first into the substituted side of cyclodextrin, while the hydroxyl terminus comes first into the substituted side of cyclodextrin in the case of type II. For making the comparison between the location of the cholesterol molecule along the center of geometry of the sugar ring and the termini of substituents of the HPβCD easier, different reaction coordinates or anchor atoms are employed for each of type I and type II complexes. The structure is sampled along the reaction coordinate from R = 1.0 Å to R = 15.5 Å, every 0.5 Å interval. We define another reaction coordinate for the binding mode prediction of the 2:1 complex of the HPβCD and cholesterol. For the 2:1 complex, the reaction coordinate is defined as the distance between the center of mass of the ring consisting of glycosidic-oxygen atoms of one HPβCD and that of the other HPβCD. For convenience, we refer to the substituted side, which has secondary hydroxyl groups, of HPβCD as “Head”,

while the other side, which has primary hydroxyl groups, as “Tail” as shown in Figure 3. Then, four types of the mutual orientation are conceivable between the two HPβCDs and the cholesterol: Head to Head (HH), Head to Tail (HT), Tail to Head (TH), and Tail to Tail (TT). The structure is sampled from R = 5.0 Å to R = 18 Å, every 0.5 Å interval. Five trajectories of 10 ns simulations with a force constant of 2 kcal/mol/Å2 on the NPT ensemble, starting from the different initial coordinates, are produced to explore the phase space in each umbrella to calculate the probability distribution of the distance between the two molecules defined above. Our procedure to obtain initial structures for each of five production runs is explained in detail on the Supporting Information. The number of simulations, simulation length, and force constants of the harmonic restraints are summarized in Table S1. The probability distribution is translated into a Potential of Mean C

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. Standard thermodynamic cycle concerning the host−guest binding process in aqueous solution. Δμrec, Δμlig, and Δμcom are the solvation free energy of the receptor, ligand, and complex structures, respectively.

solution. The initial structures for 10 sets of the production run are obtained based on the same procedure for the sampling procedure of the fluctuated ensemble of the complex state but does not include any restraint potential for all the processes. Determination of the Dominant Binding Mode. For determining the dominant binding mode, the MM/3D-RISMKH method is applied for estimating relative binding free energy among the two binding modes for the 1:1 complex, and four binding modes for the 1:2 complex. The binding free energy ΔGbind is estimated by the MM/3D-RISM-KH method according to a standard thermodynamic cycle (Figure 4) described in the previous section, using the equation ΔGbind = ΔGsolute + ΔGdesolve, where ΔGsolute and ΔGdesolve denote the binding free energy in gas phase and the desolvation free energy change in the binding process, respectively. Desolvation free energy is defined by ΔGdesolve = Δμcom − Δμrec − Δμlig, where Δμcom, Δμrec, and Δμlig denote the solvation free energies of the complex, receptor, and ligand, respectively. The binding free energy in gas phase is defined by ΔGgas = ΔEsolute TΔSsolute, where ΔEsolute denotes the change of interaction energies among atoms in the receptor and ligand upon binding, and ΔSsolute represents the corresponding change of the structural entropy. Therefore, the binding free energy is determined by a balance between the binding free energy in gas phase and the desolvation free energy. The detailed explanation of the 3DRISM-KH theory is provided elsewhere.39,40 Standard errors for each thermodynamic variable are evaluated by the bootstrap method.41 Block length is determined as 4000 snap shots corresponding to the trajectory for 80 ns. A thousand sets of bootstrap samples are generated. Standard errors for the difference of thermodynamic values along the binding process

Force (PMF) profile by means of the weighted histogram analysis method (WHAM).36−38 Here it would be worthwhile to clarify the terminologies concerning “binding modes” between host and guest molecules, since the terminology used here is slightly different from that used in our previous study. In the previous study, which focused on the binding affinity between a 1:1 pair of host and guest molecules, we were concerned with a single binding mode that has the deepest well in the PMF along the reaction coordinate. We referred to the binding mode as “most probable binding mode (MPBM).” However, in the present work, we have several different binding modes depending on the number of host molecules, 1:1 and 1:2, and their orientations. Therefore, we use terminologies that are slightly different from the previous work. The structure at each minimum of PMF along the reaction coordinate is defined as a probable binding mode (PBM). After having determined all the PBMs along the reaction coordinate, we screen the MPBM from the PBMs, by calculating the binding free energy based on the MM/3D-RISM-KH method. Sampling of the Fluctuated Ensemble of the Candidates of the Binding Mode. After having determined PBMs for the two types of orientation of the 1:1 complex and the four orientations HH, HT, TH, and TT of the 1:2 complex, we ran further MD simulations for each PBM to explore the structural fluctuation without any restraint potential to avoid the artifact that may be introduced by restraint potential. An ensemble of the bound state was produced by 10 trajectories of 80 ns MD simulations on the NPT ensemble starting from the different initial coordinates, and we got 40 000 snapshots for the free energy calculation. The procedure to obtain initial structures for each of 10 production runs is explained in detail in the Supporting Information. In contrast to the previous work, this procedure does not include the restraint potential, since the PMF of the binding mode of the cyclodextrins and cholesterol has a high energy barrier, which separates the bound state and the unbound state. The structural ensembles of the ligand and the receptor isolated in solution can be readily generated with an ordinary MD simulation, since cyclodextrin does not show large structural change unlike a protein. We regard the 40 000 snapshots generated from 10 sets of 80 ns MD simulations for the ligand and the receptor immersed in a TIP3P octahedral water box as an ensemble of the ligand and receptor isolated in

is estimated as σcom 2 + σrec 2 + σlig 2 . σcom corresponds to the standard error for a thermodynamic variable of the complex state. The free energy difference or the binding free energy ΔGbind consists of three contributions: the potential energy change ΔEsolute of the host and guest molecules, entropy change ΔSsolute of the solute molecule, and the desolvation free energy ΔGdesolve. ΔEsolute is estimated based on the Amber force field. ΔGdesolve is estimated based on the 3D-RISM-KH theory with the partial molar volume correction proposed by Palmer et al.30 ΔSsolute consists of two contributions, the external entropy and the vibrational entropy. The external entropy, or the D

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. PMF and the probability distribution of the binding mode, obtained from the MD simulation, for the type I (a) and II (b) of the 1:1 complex of HPβCD and cholesterol. The blue line corresponds to the PMF calculated by the umbrella sampling. The orange line corresponds to the probability distribution of an ensemble of the candidate of the binding mode explored by the MD simulation. (c, d) Representative structures of the candidate of the type I binding mode around the 5.5 and 12.5 Å on the reaction coordinate. (e) A representative structure of the candidate of the type 2 binding mode around the 9 Å on the reaction coordinate.

representative snapshot around each minimum of the PMF in each system is also shown. The position of the minima in the PMF in Figure 5 corresponds to the peak of the probability distribution of an ensemble of the candidate binding mode. These profiles demonstrate that the protocol based on the umbrella sampling works well to identify the candidate of the most probable binding mode and that the succeeding MD simulation also works well to explore the structural fluctuation of the candidate binding mode. As shown in Figure 5a, two minima are observed in the PMF of the type I complex around R = 5.5 and R = 12.5 Å along the reaction coordinate. There is a barrier between the two minima in the PMF. The region R = 7 to R = 10 Å, along the reaction coordinate, corresponds to where the cyclodextrin ring is located on the steroid skeleton. On the one hand, the results indicate that the radius of the pore of the βCD ring is slightly smaller than that of the steroid skeleton, so that the βCD ring is preferred to stay at the edge of the steroid skeleton. On the other hand, a single minimum is observed in the PMF of the type II complex around R = 9 Å along the reaction coordinate. There is a shallow minimum around R = 1 Å along the reaction coordinate. The minimum corresponds to the structure in which the cyclodextrin ring and the edge of the steroid skeleton make contact, but the substituents of the βCD ring do not make contact with the cholesterol. It is a conceivable reason why only one minimum is observed in the PMF of the type II complex of the HPβCD and cholesterol, in contrast to the type I complex. To check if it is true or not, we initiated the MD simulation from the initial structure, which has a type II binding mode located around R = 1 Å without any constraint. Then, the cyclodextrin ring displaced from the initial structure and moved to the position around R = 9 Å. The results of the binding free energies for all the candidate binding modes of the 1:1 complex between the HPβCD and

translational and rotational entropies, is estimated by the classical statistical mechanical formulas, based on a rigid rotor assumption.42 The vibrational part is estimated by means of the method proposed by the Chong and Ham.43 Parameters. All simulations are performed using the AMBER 14 software package. Production runs are conducted at 300 K and 1 bar using periodic boundary conditions and NPT ensemble. In this study, the generalized Amber force field (GAFF)44 with am1-bcc charge is employed for the HPβCD and cholesterol.45,46 The SHAKE algorithm47 was applied to keep the bond length involving the hydrogen atoms constant. The particle mesh Ewald algorithm48 was applied to handle the long-range electrostatic interactions accurately. All the 3DRISM-KH calculations were performed in a 64 × 64 × 64 Å3 cubic box with1283 grids. Convergence of the 3D-RISM-KH calculation was judged by ∑s (hs i + 1 − hs i)2 /N < 1 × 10−7, where hsi is the total correlation functioin at grid s in the ith iteration, and N is the number of grids. The TIP3P model is chosen for water with the modified Lennard-Jones (LJ) parameters for hydrogen atoms: ε = 0.0152 kcal mol−1 and σ = 1.236 Å. The calculations are conducted in pure water. It is our future concern to extend the study to aqueous solutions including cosolvents such as eletctrolytes.

■

RESULTS Candidates of the Probable Binding Mode for the 1:1 Complex. The PMF and the probability distribution of the binding mode, obtained from the MD simulation, for the 1:1 complex of HPβCD and cholesterol are depicted in Figures 5. The blue line corresponds to the PMF calculated by the umbrella sampling. The orange line corresponds to the probability distribution of an ensemble of the candidate binding mode produced by the MD simulation. The E

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 1. Binding Free Energy and its Componentsa Obtained from MM/3D-RISM-KH Method for Each Candidate of the Binding Mode of the 1:1 Complex HPβCD type1 (12.5 Å) HPβCD type1 (5.5 Å) HPβCD type2

ΔEsolute (kcal/mol)

−TΔSsolute (kcal/mol)

ΔGdesolve (kcal/mol)

ΔGbind_cal (kcal/mol)

−46.72 ± 0.34 −48.03 ± 0.35 −50.45 ± 0.34

25.71 ± 3.66 18.29 ± 3.83 20.60 ± 3.71

27.84 ± 0.19 28.79 ± 0.19 30.45 ± 0.20

6.83 ± 3.69 −0.95 ± 3.87 0.60 ± 3.76

ΔGbind_cal, ΔEsolute, TΔSsolute, and ΔGdesolve denote the binding free energy, the potential energy change of the solute molecule, the entropy change of the solute molecule, and the desolvation free energy, respectively. Each component of the translational, rotational, and vibrational part of the entropy change is shown in Table S3.

a

Figure 6. PMF and the probability distribution of the binding mode, obtained from the MD simulation, for the (a) HH, (b) HT, (c) TH, and (d) TT types of the 2:1 complex of HPβCD and cholesterol. The blue line corresponds to the PMF calculated by the umbrella sampling. The orange line corresponds to the probability distribution of an ensemble of the candidate of the binding mode explored by the MD simulation. Representative structures of the candidate of the binding mode around a peak of the probability distribution are also shown in each figure.

Table 2. Binding Free Energy and Its Componentsa Obtained from MM/3D-RISM-KH Method for the Four Types of the Candidate of the Binding Mode, HH, HT, TH, and TT target complex

ΔEsolute (kcal/mol)

HH HT TH TT

−104.58 ± 0.48 −114.54 ± 0.44 −99.37 ± 0.44 −99.96 ± 0.43

−TΔSsolute (kcal/mol) 29.26 32.66 35.20 39.71

± ± ± ±

5.95 6.03 5.84 6.12

ΔGdesolve (kcal/mol)

ΔGbind_cal (kcal/mol)

± ± ± ±

−3.18 ± 5.98 −6.41 ± 6.07 1.13 ± 5.87 1.47 ± 6.16

72.13 75.47 65.30 61.72

0.29 0.24 0.24 0.25

ΔGbind_cal, ΔEsolute, TΔSsolute, and ΔGdesolve denote the binding free energy, the potential energy change of the solute molecule, the entropy change of the solute molecule, and the desolvation free energy, respectively. Each component of the translational, rotational, and vibrational part of the entropy change is shown in Table S4.

a

cholesterol are shown in Table 1. The type I complex that corresponds to the minimum around R = 5.5 Å in the PMF shows the best binding affinity. The type II complex around R = 9 Å along the reaction coordinate shows the second-best binding affinity. Since these values of the binding free energy corresponding to the two types of the complex are comparable, those complexes may coexist in the solution. The type I complex corresponding to R = 12.5 Å along the reaction coordinate is less stable compared with the other binding modes. As shown in Figure S1a, the auto-correlation function of the potential energy of the HPβCD−cholesterol complex decays almost instantaneously to 0. This profile demonstrates that our simulation time is enough to estimate the reliable potential energy change. However, the change in the structural entropy of the solute molecules still has a certain amount of error. Longer simulation might have given better results. However, there is an unsettled question concerning the convergence of the structural entropy due to the molecular simulation: how long does it take to converge the structural entropy.26,49

Therefore, the results concerning the entropy should be considered as preliminary. Candidates of the Probable Binding Mode for the 2:1 Complex. The PMF and the probability distribution for the four types of the 2:1 complex of HPβCD and cholesterol, HH, HT, TH, and TT, obtained from the MD simulation, are depicted in Figure 6. The orange line corresponds to the probability distribution of an ensemble of the binding mode sampled by the MD simulation, while the blue line is the corresponding PMF. The representative snapshots around the minimum of the PMF in each system are also shown. Note that the position of the minima in the PMF corresponds to the peak of the probability distribution of the candidate binding mode. These profiles demonstrate that the protocol based on the umbrella sampling works well to identify the candidate of the most probable binding mode, and the successive MD simulation also works well to explore the structural fluctuation of the candidate binding mode. All the representative structures show that the cyclodextrin ring tends to be located at the edge of the steroid skeleton. The F

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 7. Projected structural distribution of each binding mode of the 2:1 complex onto the reaction coordinate defined for the 1:1 complex. A, B, and D are the position that showed high peaks in the distribution of the 1:1 complex along the reaction coordinate. The position C corresponds to the unstable shallow minimum of the PMF of the type II of 1:1 complex. The green line describes the distribution of the HH complex along the reaction coordinate. The yellow line describes the distribution of the HT complex on the reaction coordinate. The red line describes the distribution of the TH complex along the reaction coordinate. The purple line describes the distribution of the TT complex on the reaction coordinate.

result indicates again that the radius of the pore of the βCD ring is slightly smaller than the radius of the steroid skeleton, and the βCD ring is preferred to stay on the edge of the steroid skeleton. These four types of the 2:1 complex seem to have the cyclodextrin ring at the edge of the steroid skeleton and to form a stable complex. However, the predicted binding affinity shows an interesting feature described in the following. The predicted binding free energies, or the affinity, of all the possible binding modes for the 2:1 complex between the HPβCD and cholesterol are shown in Table 2. The results indicate that the HT dimer shows the highest affinity to cholesterol, −6.42 kcal/mol, and the HH dimer shows the second-best affinity, −3.13 kcal/mol. On the one hand, since the difference of these values is not so large, it is very likely that these complexes may coexist in the solution. On the other hand, the TH and TT dimers show rather large values for the binding free energy compared with the HT complex, indicating that the complex may not exist in solution so much. As shown in Figure S1b, the relaxation time of the autocorrelation function of the potential energy for the HH 2:1 complex rapidly converges to small value; −0.2 to 0.2. This profile demonstrates that our simulation time is enough to estimate the reliable potential energy change. However, there are small oscillations in the correlations, the periodicity of which is 10−40 ns. This result indicates that there is a largescale fluctuation in the complex of the three solute molecules, the periodicity of which is 10−40 ns. This large scale of fluctuation is caused by the change in the intermolecular distance among the three molecules. Such large-scale fluctuation will make contribution to the structural entropy. We believe that is the most important reason why the convergence of the structural entropy requires a much longer simulation time compared to the other quantities such as energy.26,49 Shown in Tables S5 and S6 is the binding free energy in which ΔGdesolve is estimated from GBSA or PBSA.50−55 The results for the binding free energy due to both GBSA and PBSA

substantially shift toward positive side from the MM/3DRISM-KH methods especially for the 2:1 complex, and the results are worse than those from the MM/3D-RISM-KH methods. The results demonstrate that only the 1:1 complex is observed in the solution, and it is inconsistent with the experimental observation. Intrinsic problem of the continuous models is the lack of the microscopic structure of the solution.56 The results from the continuum models may be improved if one adjusts the empirical parameters, the dielectric constant, the size of atoms, etc., included in the model so as to produce a better agreement with the experiments. However, there will be no unequivocal results from the continuum models, since the adjustable parameters used in the models are equivocal.56 To unveil the reason why the HT and HH dimers make distinctly stable complex with cholesterol, compared with the other complexes, we project the structural distribution of each binding mode of the 2:1 complex onto the reaction coordinate defined for the 1:1 complex. The reaction coordinate can be defined by hypothetically removing one of the HPβCDs from the 2:1 complex and by using the same definition for the distance between HPβCD and cholesterol with that used for the case of the 1:1 complex. Since there are two types of the mutual orientation for the 1:1 complex of HPβCD and cholesterol, Type I and Type II, we have eight different combinations for the projection: HH, HT, HT, and TT for Type I and HH, TH, TH, and TT for Type II. The results are shown in Figure 7, along with the PMF, which is calculated from the umbrella sampling method. In the figure, the structures of the 1:1 complex, which show high peaks in the distribution along the reaction coordinate, are marked by A, B, C, and D (see Figure 5). As is obvious from the figure, the HT complex (yellow) has peaks in the distribution at the two positions along the reaction coordinate, which correspond to the stable structures of the 1:1 complex, A and B. The two peaks in the PMF also overlap largely with the two minima in the profile of the PMF. The HH dimer (green) G

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

stable complex with the whole position of the steroid skeleton, in contrast to our finding. We think this difference originates from the difference of the structure of CDs, between βCD and HPβCD. In the work of Lopez et al., only the HH dimer of βCDs can make a stable binding mode with a cholesterol. On the one hand, the CD dimer is stabilized by the 14 hydrogen bonds between the two CD molecules. On the other hand, HPβCDs in our study cannot make such hydrogen bonds, since the second hydroxyl group is substituted by a 2-hydroxypropyl group. The substitution would induce change of the conformation to the cyclodextrin ring, since the hydrophobic moiety of the 2-hydroxypropyl groups tends to interact with each other. The difference in the structures between βCD and HPβCD would have made the difference in the binding modes of the cholesterol and βCDs. Some experimental papers suggest that the HH, HT, and TH dimers of the methylated CDs might make 2:1 complexes with cholesterol.57−59 Our study, too, demonstrates that the HT and HH dimers make stable complexes with cholesterol but not the TH dimer. Two reasons are conceivable for the contradiction: the difference in the length of the substituent, and the CD molecule in our calculation in which all the second hydroxyl groups on the glucopyranose units are substituted by a 2hydroxypropyl group. When HPβCD forms a TH or TT type (2:1) complex, the cyclodextrin substituents located near the hydroxyl group of the cholesterol cannot make strong interactions with the molecule. Our results indicate that these substituents contribute to destabilize the TH and TT types of the complex. If the cyclodextrin has small or no substituent, the TH and TT dimers may exist with a reasonable propensity. In contrast with our work, HPβCD used in experiments is a mixture of the cyclodextrins, the substitution degree and position of which are different. If the second and sixth hydroxyl groups are substituted, the difference of the binding affinity among HH, HT, TH, and TT should be smaller than the case in which the second hydroxyl groups are substituted. Effect of the substitution degree and position on the binding mode and affinity will be studied in our future work.

has a rather broad peak in the distribution at the position corresponding to the structure B in the type I orientation. It has another broad peak in the distribution at the position R = 13− 14 Å in the type II orientation, but the position is slightly off from the stable position D for the 1:1 complex due to the collision of the hydroxypropyl groups of the two cyclodextrin rings. The TH dimer (red) has two distinct peaks in the distribution, one of which corresponds to the stable structure of the 1:1 complex, D. However, the position C is in the region where the HPβCD is difficult to stay. The TT dimer (purple) also has two distinct peaks in the distribution, one in the curve corresponding to the Type I orientation and the other to the Type II orientation. One of the two peaks comes right on the positions corresponding to the stable structures of the 1:1 complex, A. However, the structure corresponding to C is in the region where the HPβCD is difficult to stay. To summarize, the HT and HH dimers satisfy the criteria for the stability; that is, the two peaks in the distribution correspond to the two stable structures of the 1:1 complex, although one of the distributions for the HH dimer is slightly off from the stable position. We believe it is the reason why the HT and HH pairs show distinctly high affinity with cholesterol.

■

DISCUSSION For the case of the 1:1 complex, our results indicate that both the type I and II binding modes can be formed, since the difference in the binding free energies between the two types is marginal. On the contrary, in the case of the 2:1 complex, our study predicts that the HT and HH dimers have greater affinity to cholesterol compared with the other two dimers. Our results also indicate that the βCD ring prefers to stay at the edge of the steroid skeleton, irrespective of the type of binding mode. Unfortunately, there have not been many experimental studies performed concerning the binding mode between cholesterol and HPβCD that would be directly compared with the results in the present work. Nonetheless, there are experimental and computational studies for some other βCD derivatives, βCD itself, dimethyl(DM)βCD, and trimethyl(TM)βCD, that are closely related to the present work. HPβCD in our study has all the second hydroxyl groups on the glucopyranose units substituted by a 2-hydroxypropyl group. DMβCD is a CD derivative, the second and sixth hydroxyl groups of which are substituted by a methyl group, while TMβCD is a CD derivative in which the second, third, and sixth hydroxyl groups are substituted by a methyl group. In the following, we make contact our results for HPβCD with those for the different βCD derivatives. The 2:1 complex of βCD, dimethyl(DM)βCD, and trimethyl(TM)βCD with a cholesterol have been studied experimentally to elucidate the binding mode.57−59 The studies commonly conclude that the cyclodextrins interact preferentially with the 18th, 19th, 26th, and 27th methyl groups on a cholesterol molecule. These results suggest that the alkyl-tail side of the cholesterol is enclosed by one of the DM or TM cyclodextrins, while the hydroxyl side of the steroid skeleton is enclosed by the other cyclodextrin. Those findings indicate that the βCD rings prefer to stay on the edge of the steroid skeleton, in harmony with our results. The computational study by Lopez et al. deduced a conclusion apparently different from that described above.21 They found that tightly dimerized βCDs can extract a cholesterol molecule with a barrierless process from a model lipid membrane. Their results indicate that the CD rings make a

■

CONCLUSION In the present study, we have identified the dominant binding mode of the cholesterol and HPβCD by applying the protocol that has been developed in previous work. Two types of complexes of HPβCD with cholesterol are considered: one-toone (1:1) complex and two-to-one (2:1) complexes. It is predicted that the 1:1 complex makes three types of stable binding mode in solution, depending on the mutual orientation, in which the βCD ring tends to be located at the edge of the steroid skeleton. For the 2:1 complex, there are four different types of the complex conceivable depending on the orientation between the two HPβCDs: head-to-head (HH), head-to-tail (HT), tail-to-head (TH), and tail-to-tail (TT). The HT and HH cyclodextrin dimers show higher affinity to cholesterol compared to the other dimers and to all the binding modes of 1:1 complexes. The physical reason why the HT and HH dimers have higher affinity compared to the other complexes is considered as follows. There are three positions along the cholesterol chain at which a single cyclodextrin is bound to make a stable complex, depending on their orientation. In case of the HT and HH dimers, the position of each CD in the dimer along the cholesterol chain comes right on or close to one of the positions where a single CD makes a stable complex. Namely, the HT and HH dimers are stabilized cooperatively at H

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(3) Poorthuis, B. J. H. M.; Wevers, R. A.; Kleijer, W. J.; Groener, J. E. M.; de Jong, J. G. N.; van Weely, S.; Niezen-Koning, K. E.; van Diggelen, O. P. The Frequency of Lysosomal Storage Disease in The Netherlands. Hum. Genet. 1999, 105, 151−156. (4) Pinto, R.; Caseiro, C.; Lemos, M.; Lopes, L.; Fontes, A.; Ribeiro, H.; Pinto, E.; Silva, E.; Rocha, S.; Marcão, A.; et al. Prevalence of Lysosomal Storage Diseases in Portugal. Eur. J. Hum. Genet. 2004, 12, 87−92. (5) Carstea, E. D.; Morris, J. A.; Coleman, K. G.; Loftus, S. K.; Zhang, D.; Cummings, C.; Gu, J.; Rosenfeld, M. A.; Pavan, W. J.; Krizman, D. B.; et al. Niemann-Pick C1 Disease Gene: Homology to Mediators of Cholesterol Homeostasis. Science 1997, 277, 228−231. (6) Naureckiene, S.; Sleat, D. E.; Lackland, H.; Fensom, A.; Vanier, M. T.; Wattiaux, R.; Jadot, M.; Lobel, P. Identification of HE1 as the Second Gene of Niemann-Pick C Disease. Science 2000, 290, 2298− 2301. (7) Ramirez, C. M.; Liu, B.; Aqul, A.; Taylor, A. M.; Repa, J. J.; Turley, S. D.; Dietschy, J. M. Quantitative Role of LAL, NPC2, and NPC1 in Lysosomal Cholesterol Processing Defined by Genetic and Pharmacological Manipulations. J. Lipid Res. 2011, 52, 688−698. (8) Camargo, F.; Erickson, R. P.; Garver, W. S.; Hossain, G. S.; Carbone, P. N.; Heidenreich, R. A.; Blanchard, J. Cyclodextrins in the Treatment of a Mouse Model of Niemann-Pick C Disease. Life Sci. 2001, 70, 131−142. (9) Irie, T.; Uekama, K. Pharmaceutical Applications of Cyclodextrins. III. Toxicological Issues and Safety Evaluation. J. Pharm. Sci. 1997, 86, 147−162. (10) Davidson, C. D.; Fishman, Y. I.; Puskás, I.; Szemán, J.; Sohajda, T.; McCauliff, L. A.; Sikora, J.; Storch, J.; Vanier, M. T.; Szente, L.; et al. Efficacy and Ototoxicity of Different Cyclodextrins in NiemannPick C Disease. Ann. Clin. Transl. Neurol. 2016, 3, 366−380. (11) Ory, D. S.; Ottinger, E. A.; Farhat, N. Y.; King, K. A.; Jiang, X.; Weissfeld, L.; Berry-Kravis, E.; Davidson, C. D.; Bianconi, S.; Keener, L. A.; et al. Intrathecal 2-Hydroxypropyl-β-Cyclodextrin Decreases Neurological Disease Progression in Niemann-Pick Disease, Type C1: A Non-Randomised, Open-Label, Phase 1−2 Trial. Lancet 2017, 390, 1758−1768. (12) Mallinckrodt Pharmaceuticals, http://www.mallinckrodt.com/ research/technology/pipeline, accessed Mar 19, 2018. (13) McCauliff, L. A.; Xu, Z.; Storch, J. Sterol Transfer between Cyclodextrin and Membranes: Similar but Not Identical Mechanism to NPC2-Mediated Cholesterol Transfer. Biochemistry 2011, 50, 7341− 7349. (14) Saad, H. Y.; Higuchi, W. I. Water Solubility of Cholesterol. J. Pharm. Sci. 1965, 54, 1205−1206. (15) Haberland, M. E.; Reynolds, J. a. Self-Association of Cholesterol in Aqueous Solution. Proc. Natl. Acad. Sci. U. S. A. 1973, 70, 2313− 2316. (16) Soga, M.; Ishitsuka, Y.; Hamasaki, M.; Yoneda, K.; Furuya, H.; Matsuo, M.; Ihn, H.; Fusaki, N.; Nakamura, K.; Nakagata, N.; et al. HPGCD Outperforms HPBCD as a Potential Treatment for Niemann-Pick Disease Type C during Disease Modeling with iPS Cells. Stem Cells 2015, 33, 1075−1088. (17) Gallicchio, E.; Levy, R. M. Advances in All Atom Sampling Methods for Modeling Protein-Ligand Binding Affinities. Curr. Opin. Struct. Biol. 2011, 21, 161−166. (18) Chodera, J. D.; Mobley, D. L.; Shirts, M. R.; Dixon, R. W.; Branson, K.; Pande, V. S. Alchemical Free Energy Methods for Drug Discovery: Progress and Challenges. Curr. Opin. Struct. Biol. 2011, 21, 150−160. (19) Abel, R.; Mondal, S.; Masse, C.; Greenwood, J.; Harriman, G.; Ashwell, M. A.; Bhat, S.; Wester, R.; Frye, L.; Kapeller, R.; et al. Accelerating Drug Discovery through Tight Integration of Expert Molecular Design and Predictive Scoring. Curr. Opin. Struct. Biol. 2017, 43, 38−44. (20) Yu, Y.; Chipot, C.; Cai, W.; Shao, X. Molecular Dynamics Study of the Inclusion of Cholesterol into Cyclodextrins. J. Phys. Chem. B 2006, 110, 6372−6378.

the two positions along the cholesterol chain. However, HPβCD is stabilized just at one of the three positions along the cholesterol chain in the case of the (1:1) complex. We believe the insights obtained here concerning the binding mode and binding affinity between the cholesterol and cyclodextrin derivative will be quite useful to analyze a therapeutic effect of cyclodextrin derivatives on the NiemannPick disease type C, and to find a better drug, which has less toxicity.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b02098. Procedure to obtain initial structures for production runs of the umbrella sampling. Summary of the number of simulations, simulation length, and force constants of the harmonic restraints for the umbrella sampling method. Procedure to obtain initial structures for production runs for sampling fluctuation around the possible binding modes. Summary of the number of simulations, simulation length, and force constants of the harmonic restraints for the sampling procedure of the fluctuated ensemble. The contributions of each entropic component of TΔSsolute for each candidate of the binding mode of the 1:1 complex. The contributions of each entropic component of TΔSsolute for each candidate of the binding mode of the 2:1 complex. Binding free energy and its components obtained from MM/GBSA and MM/PBSA method for each candidate of the binding mode of the 1:1 complex. Binding free energy and its components obtained from MM/GBSA and MM/PBSA method for each candidate of the binding mode of the 2:1 complex. Auto correlation function of the 10 sets of 80 ns simulation for (a) the structures around 12.5 Å on the reaction coordinate of the type 1of the 1:1 complex, and (b) HH type of the 2:1 complex (PDF) Representative structures for the complex, ligand and host molecules. (ZIP)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. (M.S.) *E-mail: [email protected]. (F.H.) ORCID

Masatake Sugita: 0000-0002-6298-4055 Hidetoshi Arima: 0000-0001-5498-4234 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS F.H. and M.S. are grateful to Molecular Design Frontier Co. Ltd. for a financial support. F.H. is a fellow of Toyota Physical and Chemical Research Institute.

■

REFERENCES

(1) Vanier, M. T. Niemann-Pick Disease Type C. Orphanet. J. Rare. Dis. 2010, 5, 16. (2) Meikle, P. J.; Hopwood, J. J.; Clague, A. E.; Carey, W. F. Prevalence of Lysosomal Storage Disorders. J. Am. Med. Assoc. 1999, 281, 249−254. I

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(42) McQuarrie, D. A. Statistical Mechanics, 2nd ed.; University Science Books: Sausalito, CA, 2000. (43) Chong, S.-H.; Ham, S. Protein Folding Thermodynamics: A New Computational Approach. J. Phys. Chem. B 2014, 118, 5017− 5025. (44) Wang, J.; 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. (45) Jakalian, A.; Bush, B. L.; Jack, D.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21, 132−146. (46) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623−1641. (47) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (48) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An W log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (49) Sawle, L.; Ghosh, K. Convergence of Molecular Dynamics Simulation of Protein Native States: Feasibility vs Self-Consistency Dilemma. J. Chem. Theory Comput. 2016, 12, 861−869. (50) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129. (51) Schaefer, M.; Karplus, M. A Comprehensive Analytical Treatment of Continuum Electrostatics. J. Phys. Chem. 1996, 100, 1578−1599. (52) Onufriev, A.; Bashford, D.; Case, D. A. Modification of the Generalized Born Model Suitable for Macromolecules. J. Phys. Chem. B 2000, 104, 3712−3720. (53) Onufriev, A.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-Scale Conformational Changes with a Modified Generalized Born Model. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (54) Tan, C.; Tan, Y.-H.; Luo, R. Implicit Nonpolar Solvent Models. J. Phys. Chem. B 2007, 111, 12263−12274. (55) Rastelli, G.; Rio, A. d.; Degliesposti, G.; Sgobba, M. Fast and Accurate Predictions of Binding Free Energies Using MM-PBSA and MM-GBSA. J. Comput. Chem. 2009, 31, 797−810. (56) Hirata, F.; Redfern, P.; Levy, R. M. Viewing the Born Model for Ion Hydration through a Microscope. Int. J. Quantum Chem. 1988, 34, 179−190. (57) Ravichandran, R.; Divakar, S. Inclusion of Ring A of Cholesterol Inside the Cyclodextrin Cavity: Evidence from Oxidation Reactions and Structural Studies. J. Inclusion Phenom. Mol. Recognit. Chem. 1998, 30, 253−270. (58) Nishijo, J.; Moriyama, S.; Shiota, S.; Kamigauchi, M.; Sugiura, M. Interaction of Heptakis (2,3,6-Tri-O-Methyl)-β-Cyclodextrin with Cholesterol in Aqueous Solution. Chem. Pharm. Bull. 2004, 52, 1405− 1410. (59) Castagne, D.; Dive, G.; Evrard, B.; Frédérich, M.; Piel, G. Spectroscopic Studies and Molecular Modeling for Understanding the Interactions between Cholesterol and Cyclodextrins. J. Pharm. Pharm. Sci. 2010, 13, 362−377.

(21) López, C. A.; de Vries, A. H.; Marrink, S. J. Molecular Mechanism of Cyclodextrin Mediated Cholesterol Extraction. PLoS Comput. Biol. 2011, 7, e1002020. (22) López, C. A.; de Vries, A. H.; Marrink, S. J. Computational Microscopy of Cyclodextrin Mediated Cholesterol Extraction from Lipid Model Membranes. Sci. Rep. 2013, 3, 2071. (23) Yamazaki, T.; Kovalenko, A. Spatial Decomposition Analysis of the Thermodynamics of Cyclodextrin Complexation. J. Chem. Theory Comput. 2009, 5, 1723−1730. (24) Genheden, S.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Ryde, U. An MM/3D-RISM Approach for Ligand Binding Affinities. J. Phys. Chem. B 2010, 114, 8505−8516. (25) Yesudas, J. P.; Blinov, N.; Dew, S. K.; Kovalenko, A. Calculation of Binding Free Energy of Short Double Stranded Oligonucleotides Using MM/3D-RISM-KH Approach. J. Mol. Liq. 2015, 201, 68−76. (26) Hasegawa, T.; Sugita, M.; Kikuchi, T.; Hirata, F. A Systematic Analysis of the Binding Affinity Between the Pim-1 Kinase and Its Inhibitors Based on the MM/3D-RISM/KH Method. J. Chem. Inf. Model. 2017, 57, 2789−2798. (27) Imai, T.; Hiraoka, R.; Kovalenko, A.; Hirata, F. Water Molecules in a Protein Cavity Detected by a Statistical-Mechanical Theory. J. Am. Chem. Soc. 2005, 127, 15334−15335. (28) Yoshida, N.; Phongphanphanee, S.; Maruyama, Y.; Imai, T.; Hirata, F. Selective Ion-Binding by Protein Probed with the 3D-RISM Theory. J. Am. Chem. Soc. 2006, 128, 12042−12043. (29) Imai, T.; Hiraoka, R.; Kovalenko, A.; Hirata, F. Locating Missing Water Molecules in Protein Cavities Bythe Three-Dimensional Reference Interaction Site ModelTheory of Molecular Solvation. Proteins: Struct., Funct., Genet. 2007, 66, 804−813. (30) Palmer, D. S.; Frolov, A. I.; Ratkova, E. L.; Fedorov, M. V. Towards a Universal Method for Calculating Hydration Free Energies: A 3D Reference Interaction Site Model with Partial Molar Volume Correction. J. Phys.: Condens. Matter 2010, 22, 492101. (31) Sugita, M.; Hirata, F. Predicting the Binding Free Energy of the Inclusion Process of 2-Hydroxypropyl-β-Cyclodextrin and Small Molecules by Means of the MM/3D-RISM Method. J. Phys.: Condens. Matter 2016, 28, 384002. (32) Kurkov, S. V.; Loftsson, T. Cyclodextrins. Int. J. Pharm. 2013, 453, 167−180. (33) Williams, R. O., III; Mahaguna, V.; Sriwongjanya, M. Characterization of an Inclusion Complex of Cholesterol and Hydroxypropyl-β-Cyclodextrin. Eur. J. Pharm. Biopharm. 1998, 46, 355−360. (34) Frijlink, H. W.; Eissens, A. C.; Hefting, N. R.; Poelstra, K.; Lerk, C. F.; Meijer, D. K. F. The Effect of Parenterally Administered Cyclodextrins on Cholesterol Levels in the Rat. Pharm. Res. 1991, 8, 9−16. (35) Torrie, G. M.; Valleau, J. P. Monte Carlo Free Energy Estimates Using Non-Boltzmann Sampling: Application to the Sub-Critical Lennard-Jones Fluid. Chem. Phys. Lett. 1974, 28, 578−581. (36) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for Freeenergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (37) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo Data Analysis. Phys. Rev. Lett. 1989, 63, 1195−1198. (38) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo Technique for Studying Phase Transitions. Phys. Rev. Lett. 1988, 61, 2635−2638. (39) Kovalenko, A.; Hirata, F. Three-Dimensional Density Profiles of Water in Contact with a Solute of Arbitrary Shape: A RISM Approach. Chem. Phys. Lett. 1998, 290, 237−244. (40) Kovalenko, A.; Hirata, F. Self-Consistent Description of a Metal−water Interface by the Kohn−Sham Density Functional Theory and the Three-Dimensional Reference Interaction Site Model. J. Chem. Phys. 1999, 110, 10095−10112. (41) Efron, B.; Tibshirani, R. Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy. Stat. Sci. 1986, 1, 54−77. J

DOI: 10.1021/acs.jpcb.8b02098 J. Phys. Chem. B XXXX, XXX, XXX−XXX