Letter pubs.acs.org/JPCL
Sublimation Entropy and Dissociation Constants Prediction by Quantitative Evaluation of Molecular Mobility in Crystals Sergiy O. Garbuzynskiy and Alexei V. Finkelstein* Laboratory of Protein Physics, Institute of Protein Research, Russian Academy of Sciences, 4 Institutskaya Street, 142290 Pushchino, Moscow Region, Russia S Supporting Information *
ABSTRACT: Prediction of binding free energies (or dissociation constants) is a crucial challenge for computational biochemistry. One of the main problems here consists in fast and accurate evaluation of binding entropy, which is far more time-consuming than evaluation of binding enthalpy. Here, we offer a fast and rather accurate approach to evaluate the sublimation entropy (i.e., entropy of binding of a vapor molecule to a crystal, taken with the opposite sign) from the average range of molecular movements in the solid state. To estimate this range (and the corresponding amplitude), we considered equilibrium sublimation of small organic molecules from molecular crystals. The calculations were based on experimental data on the sublimation enthalpy, pressure of saturated vapor, and structural characteristics of the molecule in question. The resulting average amplitude (0.17 ± 0.01 Å) of molecular movements in crystals was used to predict sublimation entropies and dissociation constants for sublimation of 28 molecular crystals. The results of these predictions are in close agreement with experimental values.
I
only elastic27 or all28 oscillations of protein and ligand molecules. Here, we present a fast approach to entropy calculation (and, thus, calculation of dissociation constants), which is based on comparison of the mobility of molecules in their bound and free states. For this purpose, we consider equilibrium dissociation of molecules from molecular crystals and analyze the dissociation entropy. In this work, we deal with the simplest kind of dissociation, which is sublimation (when a molecule goes from crystal to vapor). The “biologically interesting” dissociation is, of course, the case when a molecule goes from a solid body (like a crystal) to an aqueous environment. The entropy (or enthalpy) of molecules inside of crystals is the same for vapor or water surrounding the crystal, while the entropy (or enthalpy or free energy) of “free” molecules differs in vapor (which is nearly vacuum) and water. However, the additional free energy of the molecule-to-water interaction is relatively easy to estimate from an implicit accounting of the water surrounding.29−32 Therefore, the results obtained here can be used for other association/dissociation processes, such as dissociation of a protein and a ligand, a monomer protein and an amyloid fibril, and so forth. We calculate mobility (or, rather, the average range of movements) of molecules in crystals from their geometry and two experimentally measured parameters: the sublimation enthalpy and pressure of saturated vapor. Thus, we obtain
ntermolecular interactions are crucial for most processes in nature. Binding affinities of molecules are characterized by dissociation constants, which are exponents of binding free energies. Prediction of these key parameters is of great importance, in particular, for computer-aided drug design. Over the past few decades, many approaches to the problem were undertaken (see, e.g., recent reviews1,2). The most successful ones are based on molecular dynamics using various algorithms for fast sampling of the conformational space, but they are still highly computationally expensive.3−6 To successfully predict dissociation constants, one should accurately evaluate the binding free energy, which consists of the binding enthalpy and binding entropy. Molecular dynamics methods (using force fields that have been being continuously developed5,7−15) are, in principle, able to consider the free energy, which comprises both of these terms.16,17 However, molecular dynamics needs long tracking of molecular complexes and their constituent molecules18 to analyze their equilibrium fluctuating states. The binding enthalpy can be directly evaluated much faster from relatively simple energy minimization; therefore, some researchers focus on enthalpy,5,7,19,20 while others insist that it is entropy that mostly impedes adequate computation of the binding free energy.3,21 As for binding entropy, it depends on fluctuations and therefore can be strictly evaluated by molecular dynamics only after a very long tracking of molecules.18 However, there are a number of other, approximate approaches to entropy evaluation. For instance, the entropy of movements in proteins can be evaluated from dispersion of their observed conformations22−24 or those optimized by different force fields,25 from molecular surfaces hidden from water,26 or by calculating modes of either © 2017 American Chemical Society
Received: April 14, 2017 Accepted: May 30, 2017 Published: May 30, 2017 2758
DOI: 10.1021/acs.jpclett.7b00915 J. Phys. Chem. Lett. 2017, 8, 2758−2763
Letter
The Journal of Physical Chemistry Letters the range (and the amplitude) of molecular movements in a set of crystals at room temperature. We show that these ranges appear quite similar for different organic molecules. This allows us to use the averaged in-crystal amplitude to estimate the entropy of sublimation. (It should be noted that the sublimation entropy is the binding entropy of a molecule coming from vapor to crystal, taken with the opposite sign.) The predicted entropies and dissociation constants are compared with the experiment to evaluate the accuracy of the developed method. To quantitatively evaluate a change in entropy during sublimation, we consider those molecular movements that are different for molecules in a “free” (in vapor) and in a “fixed” state (in crystal). These are the movements having large amplitudes in a free state and thus are strongly hindered by a crystal (Figure 1). In a classic approximation,33 a decrease in entropy of 1 mol of vapor molecules at their fixation by a solid phase (e.g., crystal) can be calculated22,34 as follows: ⎡ δx δx δx ⎤ −ΔSsubl ≡ Scrystal − Svapor = R ln⎢ 1 2 3 ⎥ ⎢⎣ Vvapor·e ⎥⎦ nrotat. ⎡ ⎡ δβ δβ δβ ⎤ ⎤ δφi + R ln⎢ 1 2 2 3 ⎥ + R ∑ ln⎢ rotat ⎥ ⎣ (2π /K i ) ⎦ ⎣ (8π /K 0) ⎦ i=1 n vibr. ⎡ ⎤ δφj ⎥ + R ∑ ln⎢ vibr ⎥ ⎢ j = 1 ⎣ (njΔαj / K j ) ⎦
(1)
Here, R is the gas constant, and the four terms of eq 1 correspond to the four considered types of movements: (A) The first term stands for the loss of translational entropy: δx1, δx2, and δx3 are ranges of movements along three translational degrees of freedom in the solid phase (because all of these movements are restricted by close packing of molecules of a crystal, it is reasonable, to the first approximation, to set all of them equal to δx); here, a molecule is limited to “its” volume of movements Vcrystal = δx1δx2δx3. Vvapor = kBT/Pvapor is a volume per molecule in saturated vapor with pressure Pvapor and temperature T, kB being the Boltzmann constant (the use of Vvapor · e rather than simply Vvapor follows from standard statistical physics, as explained in the Supporting Information). (B) The second term stands for the loss of entropy of rotations of the molecule as a whole: δβ1, δβ2, and δβ3 are ranges (in radians) of angles of all three rotations of the molecule in the solid phase; it is reasonable to set δβk = δx/Ak, where A1, A2, and A3 are three maximal radii of the molecule; A1 · A2 · A3 is calculated from atomic coordinates (see the Supporting Information). K0 is the degeneration due to symmetry of rotations of a free molecule as a whole. (C) The third term stands for the loss of entropy of rotations inside of the molecule: nrotat is the number of free rotations around covalent bonds in the molecule in the free state (these rotations are free due to “low” barriers in the torsional potentials whose heights are 2γ < RT; see the Supporting Information). δφi is the range of angles of rotation around bond i restricted by the solid phase; it is reasonable to set δφi = δx/Bi, where Bi is the maximal radius of the smallest (having the smallest maximal radius) group rotating around the corresponding covalent bond i (see the Supporting Information).
Figure 1. Scheme of equilibrium dissociation of a molecule from a crystal; arrows show those types of movements whose ranges are different in the crystal and vapor: (A) translational and (B) rotational movements of the molecule as a whole; (C) free rotations and (D) “soft” vibrations around covalent bonds. δx is the average range of movements in the crystal (which are restricted by the neighbor molecules); further explanations are given below eq 1.
Krotat is the degeneration due to symmetry of the rotating group i of the corresponding rotation in the free state. (D) The fourth term stands for the loss of vibrational entropy: nvibr is the number of “soft vibrations” around covalent bonds. Vibrations around covalent bonds are rotations hindered, in a free molecule, by high (2γ ≥ RT) barriers of torsional potentials. Some of such vibrations (corresponding to potentials with 2γ > 70 RT; see the Supporting Information) are too small to be additionally hindered by crystals at room temperature; these “rigid” vibrations do not affect the sublimation entropy and are not taken into account here. Other vibrations are soft; they are restricted by crystals and are taken into account. δφj = δx/Bj is the range of angles of vibrations around bond j in the solid phase; here, Bj (like Bi) is 2759
DOI: 10.1021/acs.jpclett.7b00915 J. Phys. Chem. Lett. 2017, 8, 2758−2763
Letter
The Journal of Physical Chemistry Letters
Table S4), and they are excluded from further consideration. Thus, we remain with 28 selected crystals. Thermodynamic properties of the 28 crystals vary significantly; their melting temperatures are between 26.5 and 171 °C, their sublimation enthalpies are between 56 and 101 kJ/ mol, and their equilibrium vapor pressure at 25 °C ranges from 0.0006 to 1730 Pa. Despite this fact, the computed values of ln[δx] for these 28 quite different crystals appeared to be rather close; at room temperature, they vary from −0.94 to 0.03 (Figure 2), with a
the maximal radius of the smallest group rotating around the corresponding covalent bond j. nj is the number of energy minima for torsion in a free molecule around bond j; as a rule, nj = n0j , where n0j is the multiplicity of the torsional potential, but nj < n0j if some energy minima of the torsional potential are greatly elevated by other noncovalent interactions (like the energy of the cis-rotamer of the peptide bond); however, no such cases were met among the molecules considered in this work. Kvibr is the degeneration due to symmetry of the vibrating j group of the corresponding turn in the free state. Δαj is the range of angles of vibrations around bond j in the free state (see the Supporting Information). Thus, a decrease in entropy of 1 mol of molecules upon their fixation in a crystal is ⎡ δx 3 ⎤ ⎡ ⎤ δx 3 ⎥ + R ln⎢ ⎥ −ΔSsubl = R ln⎢ 2 ⎢⎣ Vvapor·e ⎥⎦ ⎣ 8π ·A1A 2 A3 /K 0 ⎦ n vibr ⎡ ⎤ ⎡ ⎤ δx δx ⎥ ⎢ + R ∑ ln⎢ rotat ⎥ + R ∑ ln vibr ⎥ ⎢ ⎣ 2π ·Bi /K i ⎦ ⎦ i=1 j = 1 ⎣ njΔαj · Bj / K j nrotat
(2)
−ΔSsubl = −ΔHsubl/T can be calculated from experimentally measured sublimation enthalpy ΔHsubl at the given temperature T (=298 K everywhere in this work), and Vvapor = kBT/Pvapor can be calculated from the pressure Pvapor of saturated vapor at the same temperature T; Δαj can be obtained from eq S.5 using potentials of rotations around covalent bonds35 (see the Supporting Information Tables S2 and S3); the values A1A2A3, K0, Bi, Bj, Krotat , and Kvibr can be obtained from molecular i j structures (taken from the Cambridge Structural Database, CSD36); the values nrotat, nvibr, and nj (and n0j ) can be obtained from the same structures and torsional potentials; all of these values are collected in Tables S2 and S3. Thus, δx, the range of movements of a molecule in the crystal, can be obtained as follows:
Figure 2. Calculated ranges δx and amplitudes of movements of molecules in the 28 considered crystals, ordered as in Tables S1−S4. The horizontal dashed line with error bars denotes the mean value.
mean of −0.38 and standard deviation of σ = ±0.24. Actually, the mean value ln[δx] is −0.38 ± 0.04; here and below, the standard error of the mean is, as usual, calculated as σ/√n, where σ is the standard deviation of ln[δx] values from the mean ln[δx] and n = 28 is the number of points. The ranges δx of molecular movements in crystals (Figure 2) are from 0.39 to 1.03 Å; their mean value ⟨δx⟩ is 0.70 ± 0.03 Å, and the standard deviation from the mean ⟨δx⟩ is ±0.15 Å. The corresponding amplitudes are (see eq S.6) 2πe ≈ 4.13 times smaller than δx and fall in the range from 0.09 to 0.25 Å, with a mean of 0.17 ± 0.01 Å and standard deviation of ±0.04 Å from the mean. The amplitudes calculated for crystals are close to the average amplitude of thermal movements observed by X-ray scattering in protein crystals.22 Having the estimated ln[δx] and the experimentally measured Vvapor = kBT/Pvapor, one can obtain estimates of sublimation entropy. To this end, eq 2 is rewritten in the form
⎧ ⎡kT ⎤ ⎡ 8π 2eA A A ⎤ ⎪ −ΔS 1 2 3 subl ⎢ B ⎥ + ln⎢ ⎥ ln[δx] = ⎨ + ln ⎪ ⎢ ⎥ R P K ⎣ ⎦ ⎣ ⎦ vapor 0 ⎩ ⎡ 2π ·B ⎤ + ∑ ln⎢ rotati ⎥ + ⎣ Ki ⎦ i=1 nrotat
×
⎡ n Δα ·B ⎤⎫ ⎪ ∑ ln⎢ j vibrj j ⎥⎬⎪ ⎢ ⎥⎦⎭ j = 1 ⎣ Kj n vibr
1 6 + nrotat + n vibr
(3) 3
Here and below, all lengths are in Å and volumes are in Å . Tables S2−S4 contain all data necessary to calculate δx. They are based on our databases13,37 containing information on 61 crystals of small molecules (see Table S1). Thirty-four of them are solid at 25 °C and could be used for our purposes (for others, extrapolations to 25 °C could be a source of errors). The experimental thermodynamic data were collected from numerous literature and web sources (also listed in the Supporting Information). It should be mentioned that the data provided by different reputable sources (CRC Handbook of Chemistry and Physics, PubChem, etc.) differ rather much in some cases. For instance, the sublimation enthalpy at 25 °C for the same compound coming from different resources sometimes differs by 20−28 kJ/mol, and the vapor pressure at 25 °C may differ 100- and even 1000-fold for the same compound! For 6 of our 34 crystals, these discrepancies are too wide (see
⎡kT ⎤ (ΔSsubl /R )theor = −(6 + nrotat + n vibr) ·ln[δx] + ln⎢ B ⎥ ⎢⎣ Pvapor ⎥⎦ ⎡ 8π 2eA A A ⎤ 1 2 3 ⎥+ + ln⎢ K ⎣ ⎦ 0 ⎡ n Δα B ⎤ j j j + ∑ ln⎢ vibr ⎥ ⎢ ⎥⎦ j = 1 ⎣ Kj
nrotat
⎡ 2πBi ⎤ ⎥ ⎣ K irotat ⎦
∑ ln⎢ i=1
n vibr
(4)
The calculated sublimation entropy values are shown in Figure 3 in comparison with the experimentally measured ones. One can see that the sublimation entropy is estimated quite accurately without any molecular dynamics simulations; most points are close to the diagonal that corresponds to ideal 2760
DOI: 10.1021/acs.jpclett.7b00915 J. Phys. Chem. Lett. 2017, 8, 2758−2763
Letter
The Journal of Physical Chemistry Letters
study; ln[A1A2A3/K0] from ≈0 to ≈4, while ΔSsubl/R varies from ≈23 to ≈42 and ln[Pvapor] from ≈−7 to ≈+7; see Tables S2 and S4). Last, the contribution of rotations and vibrations around covalent bonds (for those molecules that have covalent bonds of such type) can sometimes amount to 15% of the total sublimation entropy and sometimes is
