Entropic and enthalpic contributions to stereospecific ligand binding from enhanced sampling methods.

The stereoselective binding of R- and S-propranolol to the metabolic enzyme cytochrome P450 2D6 and its mutant F483A was studied using various computational approaches. Previously reported free-energy differences from Hamiltonian replica exchange simulations, combined with thermodynamic integration, are compared to the one-step perturbation approach, combined with local-elevation enhanced sampling, and an excellent agreement between methods was obtained. Further, the free-energy differences are interpreted in terms of enthalpic and entropic contributions where it is shown that exactly compensating contributions obscure a molecular interpretation of differences in the affinity while various reduced terms allow a more detailed analysis, which agree with heuristic observations on the interactions.

Introduction

The polymorphic cytochrome P450 (CYP) superfamily is a large group of haem-containing mono-oxygenase enzymes which owes its name to its maximum absorption at the wavelength of 450 nm in spectrophotometry. Its members play major roles in the metabolism of numerous compounds of different origins, such as drugs, food, and environmental xenobiotics. Compounds can be a substrate, an inducer, or an inhibitor of a CYP isoform. The isoforms CYP1A2, 2C9, 2C19, 2D6, and 3A4 account for the metabolism of most drugs.[1−3] Each isoform shows specificity to certain molecular features of compounds, which can be individually diversified due to the polymorphic nature of CYPs.[4−6]

A commonly prescribed drug for hypertension, propranolol, is mainly hydroxylated by CYP2D6[7−9] and is a racemic mixture of two enantiomers.[1−3] The S(−)-enantiomer (S-propranol) is more potent than the R(+)-enantiomer (R-propranolol) in blocking beta-adrenergic receptors.[10,11] The metabolism of propranolol by CYPs is also stereospecific and leads to S-propranolol concentrations that are 40–90% higher than the R-propranolol concentrations.[12,13] Experimentally, a single point mutation of phenylalanine at position 483 into an alanine, F483A, was shown to modulate stereoselective binding of propranolol to CYP2D6.[14] This observation made the propranolol–CYP2D6 system interesting.[15,16]

Experimental studies were performed to determine the relative binding free energy (ΔΔGbinding) in vitro, and computational studies have been performed to help elucidate the stereospecific binding of the enantiomers to the CYP2D6 F483A mutant. Using spectroscopic techniques, the difference in binding affinity between R- and S-propranolol was determined to amount to 0.8 kJ·mol–1 in the wild type enzyme, while R-propranolol binds 6.9 kJ·mol–1 less favorably in the F483A mutant.[15] Experimental uncertainties amount to roughly 1 kJ·mol–1. Note that this seems to be in contrast to the observation that R-propranolol is metabolized to a larger extent than S-propranolol.[12,13] The computational methods deployed by de Graaf et al.[15] for free-energy calculations using a homology model could not reproduce the ΔΔGbinding accurately as the thermodynamic cycles did not close due to inadequate sampling of phase space. Furthermore, free energies calculated from forward and backward simulations were inconsistent (hysteresis), which could be interpreted as nonreproducible conformational changes in molecules essential to the calculation, that is, the configurations sampled in forward and backward simulations were on average different. Recently, Nagy et al.[16] made an attempt to calculate the ΔΔGbinding starting from a crystal structure[17,18] using a more expensive enhanced sampling method, Hamiltonian replica exchange,[19] HRE, and longer simulations. These two changes allowed more extensive sampling of phase space than in the previous work.

In this work, the enthalpic and entropic contributions to the stereospecific binding of propranolol to wild-type CYP2D6 (CYP2D6-wt) and the CYP2D6 F483A mutant (CYP2D6-F483A) are considered. The energetic contributions are considered at four different levels, one considering the complete difference in the energy, and three different sets considering the ligands only. This approach allows a better understanding of enthalpy–entropy compensation, a thermodynamic phenomenon in which changes in the enthalpic contribution are canceled by opposing changes in the entropic contribution. Enthalpy–entropy compensation is often observed in drug development[20−24] and has been quantified using isothermal titration calorimetry.[20−23] Understanding the contributions of individual groups of molecules and/or atoms to the protein–ligand interaction and to enthalpy–entropy compensation is highly relevant when applying computational methods in drug design.[23−25] However, contributions of, e.g., the solvent-reorganization exactly cancel out in the enthalpy and entropy but are overwhelmingly large and obscure the contribution of other moieties to the binding[26−32] which are more relevant to the overall binding affinity. This can lead to a different interpretation of the enthalpic and entropic contributions and enthalpy–entropy compensation. Computational methods are complementing methods to experimental measurements and allow a more detailed decomposition of the energetic contributions. Thus, it is possible to limit the contribution of compensating energetic terms to the binding and obtain a more complete view on the binding, which is otherwise inaccessible to experiments.

In addition to decomposition of the enthalpic and entropic contributions to binding at different levels, the suitability of two enhanced sampling methods for this system, namely one-step-perturbation (OSP)[33] and one-step-perturbation with local elevation umbrella sampling (OSP+LEUS)[34] are validated against the HRE method that Nagy et al.[16] used to calculate the relative free energy (ΔΔGbinding,) associated with the perturbation of R-propranolol into S-propranolol, in CYP2D6-wt and CYP2D6-F483A. By comparing OSP, LEUS, and HRE, we compare three distinctly different methods: where OSP tries to modify the potential energy landscape such that barriers disappear and different minima can be readily sampled, LEUS adds a biasing potential to cross existing barriers in the physical potential energy landscape. HRE uses multiple replicas to mix in sampling of the conformational space from one replica to all other replicas and uses an unmodified potential (other than the modifications used for the free-energy calculation).[35] Furthermore, these three methods are commonly applied and readily available or easily implemented.

Methods

One-Step-Perturbation and Local Elevation Umbrella Sampling

The one-step-perturbation method[33] uses the Zwanzig perturbation formula (FEP)[36] to estimate the free-energy difference from a single reference simulation, that should sufficiently sample all relevant configurations for the real states. Here, the real states are R- and S-propranolol which differ by the zero-energy value of the improper dihedral E in Figure 1 (+35° and −35°, respectively). The reference state was created by setting the force constant of the dihedral angle D and improper dihedral E to zero, see Figure 1, to remove the preference for one of the stereoisomers.[37,38] The free-energy difference between the reference state, represented by the artificial Hamiltonian A, and one of the real states, represented by R or S, is calculated using Zwanzig’s perturbation formulawhere kB is the Boltzmann constant, T is the absolute temperature and ⟨⟩A indicates an ensemble average obtained using A. The local elevation umbrella sampling (LEUS)[39,40] method consists of an LE simulation followed by an US simulation. The LE simulation uses a time-dependent biasing potential for a generalized coordinate Q in one dimension, defined aswhere k is the force constant, d is the distance between grid points, F is defined as a truncated polynomial,[41] and N(Q) is the number of grid points for the coordinate Q, while n(t) is the number of times Q takes on a value corresponding to grid point b. Adding the LE bias of eq 2 to the physical Hamiltonian (q; p) yields the Hamiltonian LE, which becomes time-dependent,and leads to nonequilibrium simulations. Therefore, an equilibrium biased US simulation usingis performed, during which the values of n(t) are no longer updated.

Figure 1

Angles, dihedral angle, and improper dihedral angle of the propranolol molecule that were modified in the Hamiltonian replica exchange (HRE) and in the one-step perturbation (OSP) simulations. In HRE, the simulations were performed via a planar intermediate state, in which the lowest energy values of angle A, B, and C was changed from 111° to 120°, the force constant of dihedral angle D was set to zero, and the lowest energy value of improper dihedral angle E was changed from 35° to 0°. In OSP the force constants for D and E were set to zero. A local-elevation bias was added to E in the OSP+LUES simulations.

As demonstrated by Garate et al.[34] a suitable local elevation biasing potential can be used to ensure adequate sampling of the reference state in OSP. Considering the biasing potential as part of the artificial reference state in OSP, the corresponding Zwanzig perturbation formula now readswhere ALE is the Hamiltonian of the biased artificial reference state.

Enthalpic and Entropic Contribution Decomposition

The free-energy difference between R-propranolol and S-propranolol (ΔG), in either bound or unbound state, is calculated using different free-energy calculation methods, namely OSP, OSP+LEUS, and HRE. For the HRE, the enthalpy, ΔH, is calculated by subtracting the ensemble average of the Hamiltonian of corresponding end-states of the free-energy simulations. Hence, the entropy, ΔS, can be calculated from ΔH, and ΔG, using the Gibbs equation,The enthalpic and entropic contributions to the free energy contain large compensating terms, which may obscure a detailed analysis of the binding free energy. In the context of a free-energy perturbation, the Hamiltonian can be split into a λ-dependent perturbed part, d, and a λ-independent unperturbed part, i,It was previously shown[26,27,31,42,43] that the ensemble average of i cancels out exactly in the enthalpy and entropy, such that the Gibbs equation can also be written aswhere ΔHd = ⟨d(1)⟩ – ⟨d(0)⟩. Here, the thermodynamic quantities are studied by using different definitions of d that include different energy terms. For every choice, the ΔH is reduced by more contributions that exactly cancel out in the ΔS, and eq 8 is used to calculate the corresponding reduced entropy term. The first level (1) of reduction includes the full enthalpy, ΔH, and entropy, ΔS, of the system, and the Gibbs equation remains unchanged (eq 6). The second level (2) includes the ligand-surrounding enthalpy, ΔHls, and entropy, ΔSls, which were introduced in our previous work as generalized solvent–solute terms derived from solvation studies and which are known to converge readily.[31] Ligand-surrounding interactions consist of the nonbonded and bonded interactions within the ligand and between the ligand and the protein, and the ligand and the solvent. This approach improves convergence that is otherwise limited by contributions and errors from significantly larger surrounding–surrounding energy terms, such as solvent–solvent, protein–solvent, and protein–protein interactions. The Gibbs equation becomesThe third level (3) includes only the bonded energy terms of the ligands in the enthalpy, ΔHcov, and entropy, ΔScov. This level of reduction is valid, because all perturbed properties are reflected by bonded energy terms. The Gibbs equation becomesIn the fourth (4), and probably the most descriptive level with regards to the binding, only the energetic contributions from the sampling of all angles, dihedral angles, and improper dihedral angles that were perturbed in the HRE calculations are considered. These energetic contributions are λ-dependent and were recalculated from the coordinate trajectories of the HRE simulations at the end-states. Now, the Gibbs equation becomeswhere ΔHpt, the perturbed enthalpy term, is calculated from the sum of the energetic contributions of all perturbed angles, dihedral angles, and improper dihedral angles, and ΔSpt is the corresponding perturbed entropy term.

Analogously, enthalpic and entropic terms were calculated from the trajectories of the OSP and OSP+LEUS simulations and the relevant term of the observed Hamiltonian, term, was reweighted from the artificial reference state to the real state, which is either R-propranolol or S-propranolol, using

Simulation Setup

HRE data were taken from previous simulations by Nagy et al.,[16] and raw trajectories of the energy, free energy, and coordinates were reanalyzed.

All OSP and OSP+LEUS simulations were performed using the GROMOS11 simulation package.[44] The GROMOS++ software for analysis of biomolecular trajectories was used to assist in setting up and analyzing the simulations.[45] Force-field parameters for the protein were taken from the GROMOS 45A4 united-atom force field.[46,47] Propranolol was described using the parameters of Hritz et al.[18] Initial configurations and corresponding velocities of a propranolol molecule solvated in 1812, simple point charge (SPC) water molecules,[48] or in complex with CYP2D6 solvated in 12019 SPC and 5 Na+ molecules were taken from the end-states of the HRE simulations. Figure 2 shows propranolol in the active site of the enzyme. Simulations were performed with a constant number of particles, at a constant pressure of 1 atm and at a constant temperature of 298 K. Solvent and solute degrees of freedom were coupled separately to two temperature baths with a relaxation time of 0.1 ps using the weak-coupling method.[49] Pressure was also kept constant using the weak-coupling scheme, with a relaxation time of 0.5 ps and an estimated isothermal compressibility of 4.575 × 10–4 (kJ·mol–1·nm–3)−1. The leapfrog algorithm[50] with a time step of 2 fs was used. All bonds were constrained to their minimum energy values using the SHAKE algorithm.[51] Center of mass translation was removed every 1000 steps. Nonbonded interactions were calculated using a triple range cutoff scheme. Interactions up to a short-range distance of 0.8 nm were calculated at every time step from a pairlist that was updated every 5 steps. At pairlist construction,[52] interactions up to an intermediate range of 1.4 nm were also calculated and kept constant between updates. A reaction field contribution[53] was added to the forces and energies to account for a dielectric continuum with relative permittivity of 61 beyond the cutoff sphere of 1.4 nm.[54]

Figure 2

S-Propranolol at its binding site in CYP2D6. S-Propranolol is shown in ball and stick representation, with the heme directly located below it. Hydrogen bonds are shown as blue dashed lines. All residues of CYP2D6 within 0.45 nm of the propranolol molecule are shown and labeled. Secondary structure elements are shown in tube representation.

The free-energy calculations using HRE by Nagy et al.[16] were performed via a planar intermediate propranolol structure in which three angles, one dihedral angle and one improper dihedral were modified; see Figure 1. Transitions of dihedral angle D were seen to influence the convergence of the calculations significantly. In the OSP and OSP+LEUS simulations, the force constant of the dihedral angle D and the improper dihedral angle E were set to 0, to allow adequate sampling of the R- and S-configurations.

A 5 ns OSP simulation of propranolol in solvent was performed. Four 4 ns OSP simulations were started from configurations of R- and S-propranolol in complex with CYP2D6-wt or CYP2D6-F483A in solvent extracted from the end-states of different HRE simulations. Two 10 ns OSP+LEUS simulations of propranolol in complex with either CYP2D6-wt or CYP2D6-F483A in solvent were performed. A bias along the improper dihedral angle E with a force constant of 0.001 kJ·mol–1 (CYP2D6-wt) or 0.00034 kJ·mol–1 (CYP2D6-F483A) and a bin width of 4° (90 bins) was build up during the LE stage of the OSP+LEUS simulations (300 ps). The US stage was subsequently performed for 10 ns.

Error estimates for the averages obtained from simulations were determined from block averaging and extrapolation to infinite block length.[55] Error estimates in the thermodynamic terms are subsequently obtained from standard propagation of the error estimates on the simulation averages.[56]

Results and Discussion

OSP and OSP+LEUS Sampling

The sampling of different configurations of propranolol by OSP or OSP+LEUS is shown in Figure 3, where time series and corresponding histograms show the sampling of improper dihedral angle E during simulations, which distinguishes R- and S-propranolol. Figure 3A shows the OSP simulation of the artificially modified propranolol in solvent. The possible configurations of improper dihedral angle E were sampled with equal probability. This indicates that intramolecular interactions or the solvent are not limiting factors to the configurational sampling of the reference state. However, OSP simulations started from different end-states of the HRE simulations which contain either R- or S-propranolol in complex with the CYP2D6-wt and CYP2D6-F483A, Figure 3B–E, do not always show sufficient configurational sampling. Figure 3D shows an ideal case where both R- and S-propranolol configurations are sampled with equal and high probability. The other three cases are simulations where either R- or S-propranolol is sampled predominantly. This bias toward a single configuration hampers the calculation of free energies from these simulations. However, Figure 3F and G show that OSP+LEUS can be used to improve the sampling of R-propranolol and S-propranolol configurations during the simulations, allowing both configurations to be sampled more frequently. Remarkably, the flat propranolol was sampled significantly more by the OSP+LEUS than OSP alone, which, like the increased sampling of the R- and S-propranolol configuration, is the direct result of adding LEUS to OSP and may indicate a slight overbuilding of the biasing potential. Note however, that this state is also significantly sampled in panel A where no bias was added.

Figure 3

Time series and histograms of improper dihedral angle E (ζ) indicated in Figure 1. (A) OSP simulation of propranolol in solvent. (B and D) OSP simulations of propranolol with CYP2D6-wt using different starting configurations. (C and E) OSP simulations of propranolol with CYP2D6-F483A using different starting configurations. (F) OSP+LEUS simulations of propranolol with CYP2D6-wt protein. (G) OSP+LEUS simulations of propranolol with CYP2D6-F483A.

Free Energy between R-Propranolol and S-Propranolol

The free-energy difference between R- and S-propranolol (ΔG) is given in Table 1. ΔG calculated for propranolol in water, in complex with CYP2D6-wt, or CYP2D6-F483A using HRE, and OSP, or OSP+LEUS simulations match with differences of less than 1 kJ·mol–1, well within statistical error estimates. The excellent agreement between the different methods to calculate ΔG can be expected from robust free-energy calculation methods and suggests that sampling is adequate.

Table 1

Free Energy, Enthalpy, and Entropy Differences between R- and S-Propranolol (kJ·mol–1) Calculated Using HRE and OSP+LEUS Simulations at Different Levels of Reduction

HRE	water	wt	mutant
ΔGR→S	–0.3 ± 0.4	7.1 ± 1.1	1.5 ± 1.0
level 4: perturbed terms
ΔHpt	–0.3 ± 0.2	0.6 ± 0.1	–0.3 ± 0.2
TΔSpt	0.0 ± 0.5	–6.5 ± 1.1	–1.8 ± 1.1
level 3: covalent terms
ΔHcov	0.4 ± 0.4	3.7 ± 1.0	6.9 ± 0.6
TΔScov	0.7 ± 0.5	–3.4 ± 1.5	5.4 ± 1.1
level 2: ligand-surrounding terms
ΔHls	2.5 ± 1.4	3.2 ± 7.0	2.6 ± 2.0
TΔSls	2.8 ± 1.4	–3.9 ± 7.1	1.1 ± 2.3
level 1: full terms
ΔH	14.5 ± 6.9	–16.6 ± 34.4	–69.6 ± 24.1
TΔS	14.8 ± 6.9	–23.7 ± 34.4	–71.1 ± 24.1

Changing the configuration from R-propranolol to S-propranolol in water using HRE or OSP resulted in ΔG = −0.3 ± 0.4 and 0.2 ± 0.4 kJ·mol–1, respectively. Both estimates are zero within the statistical uncertainly as is appropriate for two enantiomers in an achiral environment. Changing of the configuration of the propranolol molecule from R to S in CYP2D6-wt is unfavorable by 7.1 ± 1.1 or 6.9 ± 1.2 kJ·mol–1, an observation that has been discussed by Nagy et al. and confirmed again by the OSP or OSP+LEUS simulations in this work. These observations seem to be in agreement with the experimental finding that R-propranolol is metabolized more efficiently than S-propranolol.[12,13] The value of ΔG is more favorable and closer to zero in the mutant, which is in agreement with the experimental observation that the mutation strongly influences the stereoselectivity.[15] Note, however that the spectroscopically determined binding affinities rather suggest values of 0.8 kJ·mol–1 for CYP2D6-wt and −6.9 kJ·mol–1 for CYP2D6-F483A.[15] The OSP+LEUS calculations reported in this work represent the third independent computational estimate of the binding affinities, strongly suggesting that the discrepancy between the computed and experimental data is not due to limited sampling. Rather, the shift by about 7 kJ·mol–1with respect to the experimental data could be due to inappropriate force-field parameters or systematic errors in the experiments. We repeat that the experimentally determined binding affinities do not agree with the observed rates of metabolism either.[12,13]

Enthalpic and Entropic Contributions to the Difference between R- and S-Propranolol

The enthalpic (ΔH) and entropic (TΔS) terms at each level of reduction are shown in Figure 4 with exact values, and estimates of the statistical error given in Table 1. The statistical error estimates are calculated from block averaging and extrapolation to infinite block length on the individual ensemble averages, followed by a proper error propagation.[55] In particular, the various enthalpy terms calculated from eq 12 suffer from large statistical uncertainties. As the size of the term increases, the error estimate increases as well. Equation 12 reweights the data, with significant weights remaining for a reduced set of the data, further increasing the uncertainty, which subsequently remains in the enthalpy difference. As such, the use of OSP and related methods for the calculation of the full enthalpic and entropic contributions to the binding process is in this particular case, and possibly in other cases, of limited applicability.

Figure 4

Free energy (red), enthalpy (green), and entropy (blue) differences between R- and S-propranolol at different levels of reduction of the enthalpic and entropic contributions. The panels left from the middle show data from the HRE simulations, and the panels the right from the middle show data from the OSP+LEUS simulations. Whiskers on each bar indicate the statistical uncertainties of the calculations; see also Table 1. (1) First level, the full enthalpy and entropy. (2) Second level, only ligand–ligand and ligand-surrounding energy terms are included. (3) Third level, only bonded interactions within the ligand are included. (4) Fourth level, only the energetic contribution of the perturbed angles, dihedral angles. and improper dihedral angles are included.

From Figure 4, one can see that the calculations using HRE or OSP+LEUS mostly follow similar trends, with notable exceptions. At the highest level of reduction, enthalpic contributions (ΔHpt) calculated using HRE and OSP+LEUS are in good agreement, and consequently, the entropic contributions (TΔSpt) also are. Altogether, this observation suggests that HRE and OSP+LEUS are both suitable for studying reduced enthalpic and entropic contributions, while explicit simulations at the end-states are required for the full enthalpy (ΔH) and entropy (TΔS).

In CYP2D6-wt, the enthalpic and entropic components of ΔG of propranolol show that across all levels of reduction, the change of R-propranolol into S-propranolol in CYP2D6-wt is entropically unfavorable. The enthalpic contribution becomes unfavorable only if the surrounding–surrounding energetic terms are excluded. This observation emphasizes once more the fact that changes in interactions between the solvent–solvent, solvent–protein, or protein–protein interactions, which are exactly compensated in the entropy, obscure the nature of the protein–ligand interactions.

In fact, some enthalpy–entropy compensation remains at the second and third level of reduction and it is only at the fourth level of reduction that no compensating contributions remain. At this level, it is clear that the free-energy difference is mostly determined by the entropy (TΔSpt).

The enthalpic components of ΔG for propranolol in CYP2D6-F483A follow a similar trend as seen in the CYP2D6-wt calculations. However, changes between levels are significantly larger, which indicates that the contributions from the exactly canceling surround–surrounding terms are more significant. Again, the entropic contribution is unfavorable at the first and fourth levels, while it changes sign between the levels. The change in the enthalpy seems to suggest that the interactions of the ligand with its environment rather than the covalent ligand–ligand interactions, make the binding enthalpically favorable. The enthalpy follows the same trend as in the wild-type over all levels of reduction, while the entropic contributions change more and become favorable at the second and third level. The increase in entropy at these levels suggest increased mobility of the propranolol molecule in the mutant and seems to originate from parts of the ligand which were not perturbed.

Relative Enthalpic and Entropic Contributions

From Table 1 it is possible to calculate ΔΔH and TΔΔS by subtracting the values of the change from R- to S-propranolol in water from the values in the protein, or by considering the differences between CYP2D6-wt and CYP2D6-F483A directly. All values of ΔΔH at the first three levels of reduction contain contributions that are exactly compensated by an opposing term in ΔΔS, which suggests that it is possible to obtain a clearer description of the protein–ligand interactions by focusing on the essential changes in the system rather than the entire system.

Focusing on the HRE data once again, the full enthalpy and entropy (ΔH and TΔS) show that perturbation of R- into S-propranolol is enthalpically more favorable for propranolol in complex with CYP2D6-F483A than for propranolol in complex with CYP2D6-wt, which in turn is enthalpically more favorable than in solvent. From the entropic point of view the situation is reversed. The preference of S-propranolol over R-propranolol binding to CYP2D6-wt is favored by enthalpy and disfavored by entropy. Both effects were more pronounced in CYP2D6-F483A simulations. However, Nagy et al. suggested, based on hydrogen bond analysis that the enthalpic contribution did not change significantly, but rather that the number of unique hydrogen bonds is higher in CYP2D6-F483A than in CYP2D6-wt, suggesting that the difference between binding to CYP2D6-wt or to CYP2D6-F483A is mostly of entropic nature. Additionally, it was observed that S-propranolol has more direct protein–ligand interactions than R-propranolol which is mostly interacting through the water network. A recent experimental observation[32] suggested that interactions with such a water network can obscure the proper interpretation of the protein–ligand interaction through compensating terms. These observation together suggest that the surrounding–surrounding terms are the major contributors to the full enthalpy and entropy terms (ΔH and TΔS).

Removing the surrounding–surrounding terms, that is, moving to the second level of reduction, reveals enthalpic and entropic contributions (ΔHls and TΔSls) that are more in line with previous qualitative observations. The enthalpic contribution in CYP2D6-wt and CYP2D6-F483A are practically identical (ΔΔHls = −0.6 kJ·mol–1), which means that surrounding–surrounding interactions, account for a large portion of the full enthalpy (ΔH). The relative entropic contribution (TΔΔSls) is 5.2 kJ·mol–1, which suggest that the favorable interaction of S-propranolol with CYP2D6-F483A is entropy driven. The differences in enthalpy (ΔHls), between the complexes and propranolol in solvent is comparable at this level, further emphasizing that removal of the exactly compensating surrounding–surrounding terms assists in clearer interpretation of the ligand–surrounding interactions.

Moving onto the third level reveals an even stronger entropy driven process (TΔΔScov = 8.8 kJ·mol–1), and the enthalpic contribution becomes more unfavorable (ΔΔHcov = 3.1 kJ·mol–1), which is not entirely unexpected. By taking only covalent ligand–ligand interactions into consideration, the rather promiscuous hydrogen-bond forming behavior of S-propranolol is excluded from the analysis, and this is reflected in the enthalpic contribution (ΔHcov). Hydrogen bonds with water molecules reflect indirect interactions of propranolol with CYP2D6 and may be related to the increase in enthalpy, assuming that hydrogen bonds between propranolol and water molecules are more transient than hydrogen bonds between propranolol and CYP2D6.

At the last level of reduction, the difference between enthalpic contributions (ΔHpt) is nearly indistinguishable (ΔΔHpt = 0.9 kJ·mol–1), but the difference in entropy remains (TΔΔSpt = 4.8 kJ·mol–1). The change of configuration from R- to S-propranolol is entropically (TΔSpt) most unfavorable in CYP2D6-wt and less unfavorable in CYP2D6-F483A, which is in line with previous computational and experimental observations, which stated that the mutation in CYP2D6-F483A affects stereoselectivity for S-propranolol.

Comparison to Configurational Entropy

In an attempt to understand the entropic contribution to the binding of S- and R-propranolol, Nagy et al. estimated configurational entropies using Schlitter’s equation[57] and found for the free energy of binding of S-propranolol in favor of R-propranolol by CYP2D6-F483A, a configurational entropy contribution of 4.8 kJ·mol–1 due to the conformations of F483, −0.9 kJ·mol–1 due to the conformations of propranolol, and 2.6 kJ·mol–1 due to the conformations of propranolol including translation and rotation within the active site. Configurational entropies are heuristic estimates based on observed (co)variances in the simulations and ignore contributions from desolvation or solvent-reorganization. Typically, only contributions for a reduced set of degrees of freedom are considered.[58] As enthalpy and entropy at the third and fourth level are not directly considering interactions with the surrounding, they may be most comparable with the configurational entropy. A TΔΔScov of 8.8 kJ·mol–1 was found at the third level, which is within kBT of 7.4 kJ·mol–1 that is the sum of the configurational entropy of F483 and propranolol with translation and rotation. At the highest level, the entropic contribution (TΔSpt) calculated from the HRE amounts to 4.9 kJ·mol–1. Although both approaches quantify entropy using different assumptions, they agree quite closely in terms of a molecular interpretation of the thermodynamic contributions to the binding.

Conclusion

OSP and OSP+LEUS were used to calculate free energies and enthalpic and entropic contributions to the binding of R- or S-propranolol to CYP2D6-wt or CYP2D6-F483A. The results showed that OSP and OSP+LEUS can be used for free-energy calculations, while it remains difficult to calculate the full enthalpy and entropy from a single simulation of the reference state, mainly due to error propagation in the reweighting process.

The enthalpic and entropic contributions were considered at four levels of reduction by excluding different energetic terms at each level. This approach further emphasizes that enthalpy–entropy compensations severely hamper the thermodynamic interpretation of protein–ligand interactions at a molecular level and shows that the reduced terms for protein–ligand interactions are more suitable tools for improving the understanding of the underlying processes in protein–ligand binding and enthalpy–entropy compensation.