{'fraction_non_alphanumeric': 0.030975409317908844, 'fraction_numerical': 0.027403754978190784, 'mean_word_length': 4.760924981791697, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Free energies of binding of R- and S-propranolol to wild-type and F483A mutant cytochrome P450 2D6 from molecular dynamics simulations

Detailed molecular dynamics (MD) simulations have been performed to reproduce and rationalize the experimental finding that the F483A mutant of CYP2D6 has lower affinity for R-propranolol than for S-propranolol. Wild-type (WT) CYP2D6 does not show this stereospecificity. Four different approaches to calculate the free energy differences have been investigated and were compared to the experimental binding data. From the differences between calculations based on forward and backward processes and the closure of thermodynamic cycles, it was clear that not all simulations converged sufficiently. The approach that calculates the free energies of exchanging R-propranolol with S-propranolol in the F483A mutant relative to the exchange free energy in WT CYP2D6 accurately reproduced the experimental binding data. Careful inspection of the end-points of the MD simulations involved in this approach, allowed for a molecular interpretation of the observed differences. Electronic supplementary material The online version of this article (doi:10.1007/s00249-006-0126-y) contains supplementary material, which is available to authorized users.


Introduction
The superfamily of cytochrome P450 (CYP) enzymes constitutes the most important phase I metabolic enzymes. CYPs generally detoxify pharmacologically active and potentially hazardous compounds, but in a number of cases non-toxic parent compounds are bioactivated into pharmacologically or toxicologically active metabolites, or procarcinogens into their ultimate carcinogens (Sheweita 2000;Vermeulen 1996). Although the expression levels of CYP2D6 represent only 2% of all hepatic CYPs, it is the second most important drug metabolizing enzyme, after CYP3A4, and involved in the metabolism of about 15-30% of the currently marketed drugs (Bertilsson et al. 2002;Zanger et al. 2004). Large interindividual differences exist in CYP2D6 activity, due to gene multiplicity and polymorphisms further emphasizing its clinical importance (Ingelman-Sundberg 2004;Oscarson and Ingelman-Sundberg 2002;Wormhoudt et al. 1999). The early identification of potential CYP2D6 substrates and prediction of their metabolism is therefore advantageous in the discovery and development of new drugs.
Experimentally, the active site of CYP2D6 can be probed by site directed mutagenesis studies as described in e.g. Flanagan et al. (2004), Keizers et al. (2004Keizers et al. ( , 2005b, Lussenburg et al. (2005) and Paine et al. (2003). Such experiments can be designed from structural models of enzyme-substrate complexes and the observed effect of mutations on metabolism and inhibition can subsequently be rationalized at an atomic resolution by modeling approaches. One mutant of CYP2D6 that was proposed based on computational models, is the F483A mutant (Lussenburg et al. 2005).
The effect of this mutation on substrate binding affinities and orientations was shown to be strongly dependent on the substrate. For example, the metabolism of MAMC was no longer observed in the F483A mutant and metabolism of bufuralol was reduced. On the other hand no effects on the MDMA metabolism were found, while for dextrometorphan several alternative metabolites were observed. The binding affinity of these compounds did not seem to be influenced significantly by the mutation. For the substrate R-propranolol, however, the affinity is reduced by a factor 20 upon mutation. Interestingly, a large stereospecific affinity is observed in the mutant CYP2D6, but not in wild-type (WT). The metabolites formed from both R-and S-propranolol remain the same upon mutation of the enzyme (unpublished data).
In this study, we explore the stereospecific binding of propranolol to wild-type CYP2D6 and the F483A mutant. The spectral dissociation constants obtained were computationally rationalized. Computational tools may offer an atomic resolution often lacking in experiments. In the present case, one expects that subtle differences in the protein-ligand interactions in WT CYP2D6 or in the F483A mutant have a determining influence on the observed affinities. In order to catch such subtle differences a static approach based on protein-ligand structures only will not likely suffice, but rather the structure and dynamics of protein-ligand interactions have to be taken into account, for instance by molecular dynamics (MD) simulation. In CYPs it was shown that docking and scoring strategies can be applied to accurately predict binding modes and to select high affinity compounds from chemical databases (e.g. CYP101, CYP2D6 and CYP3A4) (de Graaf et al. 2006;DeVoss and Demontellano 1995;DeVoss et al. 1997;Kemp et al. 2004;Keseru 2001;Zhang et al. 1998). However, the scoring functions available in docking programs show little or no correlation with (experimental) binding affinity for these CYPs Keseru 2001). In the case of MDMA binding to CYP2D6 (WT and F120A mutant) it was shown that a dynamic view on the interaction was required to rationalize experimental findings (Keizers et al. 2005a). Detailed MD simulations also offer the possibility to calculate the free energy difference between complexes of different proteins and ligands.
A remaining hurdle in the accurate calculation of CYP ligand binding affinity, however, is the diversity of binding modes of substrates, stemming from a flexible binding cavity that is large in relation to the size of the substrates. Subsequently, for many substrates the binding mode in CYPs is not as strictly defined as is the case for many other more substrate specific enzymes. This introduces sampling problems, prohibiting the application of MD based methods that can provide the most accurate free-energy calculations.
In the case of R-and S-propranolol binding to WT CYP2D6 and the F483A mutant, several approaches can be made, making use of different thermodynamic cycles (Tembe and McCammon 1984). We distinguish four of such approaches, which are exemplified in approach 1: Fig. 1 Four different approaches to calculate the differences in binding affinities of R-propranol (R) and S-propranolol (S) to wild-type (WT) and F483A mutant (F483A) CYP2D6, making use of different thermodynamic cycles (described in the text). Simulations at the initial and final stages along the solid arrows are labeled with a small letter (also referred to in the text and in Table 2) binding affinities of R-and S-propranolol in protein (WT or F483A) can be obtained ðDDG WT bind ðR; SÞ and DDG F483A bind ðR; SÞÞ: 2. Calculate the free energy of inversion of R-propranolol into S-propranolol when bound to WT and when bound to the F483A mutant: DG WT inv and DG F483A inv : From this the relative free energy of transfer from WT to F483A for R-and S-propranolol can be obtained ðDDG transfer ðR; SÞÞ: 3. Calculate the free energy of mutating the phenylalanine sidechain into an alanine at position 483 in a substrate free protein model and a substrate (R-or S-propranolol) bound protein model DDG unbound mut ; DG R mut and DG S mut : From this the relative binding affinities of the ligand (R-or S-propranolol) to WT and F483A can be obtained (DDG R bind ðWT; F483AÞ and DDG S bind ðWT; F483AÞÞ: 4. Calculate the free energy of mutating the phenylalanine side chain into an alanine at position 483 in a R-propranolol bound protein model and in a S-propranolol bound protein model: DG R mut and DG S mut : From this the relative free energy of exchanging R-propranolol by S-propranolol in WT and the F483A mutant can be obtained (DDG exchange (WT,F483A)).
Each of these approaches yields information that can be compared directly to the experimentally determined differences in binding affinities of R-and S-propranolol to WT and F483A CYP2D6. DG free inv ; representing the inversion of the stereocenter of propranolol free in solution, should equal zero. Because free energy is a state function, combinations of individual free energy terms defined above can be found that should sum up to exactly zero ðDG cycle Þ: Comparisons of calculated to experimental values and to theoretically defined values will be used to determine which of the approaches is most accurate and can be used to rationalize the stereospecificity of the F483A mutant.

Optical titrations
Dissociation constants of the substrates to the enzymes were determined by spectral titration, according to the method of Jefcoate (1978). Spectra were taken at room temperature on a Pharmacia Ultrospec 2000 spectrometer. In short, 1 ml of 0.5 lM purified enzyme in a 50 mM potassium phosphate buffer, pH 7.4 containing 10% glycerol, was divided over two cuvettes, to the first 5 ll of a 1 mM solution of R-or S-propranolol in the same buffer was added, to the second the same volume of just buffer. Difference spectra from 350 to 450 nm were taken before, and after each addition of compound. The substrates themselves did not show absorbance in this spectral range at the concentrations used. The difference between the absorbance at 390 and 425 nm was plotted against the substrate concentration to estimate the dissociation constant (K s ) using the following equation: where B is the absorbance difference (390-425) and B max is the absorbance difference when [S] fi ¥.

CYP2D6-propranolol model
Coordinates for a ligand bound model of cytochrome P450 were obtained from homology modeling as described in (de Graaf et al. 2007;Keizers et al. 2005a). Very recently a model based on the crystal structure of this isoenzyme has been published, but this involves a ligand-free conformation of the protein (Rowland et al. 2005). Our homology model was found to agree with most of the details of the recently solved substrate-free CYP2D6 crystal structure. Structural differences between the homology model and crystal structure were the same differences observed between substrate-free and substrate-bound structures of other CYPs (e.g. 2C5; Wester et al. 2003aWester et al. , 2003bWilliams et al. 2000), suggesting that these conformational changes are required upon substrate binding (de Graaf et al. 2007). For this reason we chose to rather use the ligand-bound homology model, which was validated in earlier MD simulations (Keizers et al. 2005a). Initial coordinates for R-propranolol in complex with the protein model were obtained by an automated docking approach as described in (de Graaf et al. , 2006.

Molecular dynamics simulations
Molecular dynamics simulations and free energy calculations were performed along all solid angles in the thermodynamic cycles in Fig. 1. All molecular dynamics simulations were performed using a preliminary version of the GROMOS05 programme (Christen et al. 2005). Interaction parameters were taken from the GROMOS 43A1 force field (Van Gunsteren et al. 1996). For propranolol, parameters were determined according to analogy to protein functional groups (see supplementary material). The GROMOS force field is a united atom force field, treating aliphatic hydrogen atoms together with the carbon atoms to which they are bound as a single interaction site. The wild-type protein was Eur Biophys J (2007) 36: 589-599 591 solvated in a rectangular periodic box containing 5,314 (ligand free simulation) or 5,306 (with R-propranolol bound) SPC water molecules (Berendsen et al. 1981).
For the simulations of propranolol free in solution, the simulation box contained 2,367 water molecules. After an initial energy minimization in which the protein was positionally restrained, simulations were started at a temperature of 300 K, maintaining harmonic restraints with force constant 2,500 kJ mol -1 nm -2 on all solute atoms.
The system was carefully heated up and equilibrated as described in (Keizers et al. 2005a). Bond lengths were constrained at their minimum energy values using the SHAKE algorithm (Ryckaert et al. 1977) allowing for a time step of 2 fs. Nonbonded interactions were calculated using a triple range cutoff scheme. Interactions within a short range cutoff of 0.8 nm were calculated at every time step from a pairlist that was generated every five steps. Longer range interactions within a cutoff of 1.4 nm were calculated at pairlist generation and kept constant between updates. To account for electrostatic interactions beyond the long range cutoff, a reaction field contribution to the energies and the forces was calculated (Tironi et al. 1995). The temperature was kept constant by coupling solute and solvent degrees of freedom separately to a temperature bath of 300 K and using a relaxation time of 0.1 ps (Berendsen et al. 1984). Constant pressure was maintained by isotropically coupling to a weak pressure bath of 1 atm using a relaxation time of 0.5 ps and an estimated isothermal compressibility of 4.575 · 10 -4 (kJ mol -1 nm -3 ) -1 .

Free energy calculations
Free energies were calculated along the solid arrows in Fig. 1, using the thermodynamic integration method (Beveridge and DiCapua 1989;Kirkwood 1935). The arrows in this figure connect an initial (A) and a final (B) state, which can be defined as: (1) bound to wild-type and F483A mutant CYP2D6, respectively; (3) unbound wild-type (A) versus unbound F483A mutant (B) CYP2D6; (4) wild-type (A) versus F483A mutant (B) binding R-and S-propranolol, respectively. The Hamiltonians describing the initial and final states of each process were related using a coupling parameter k, such that at k = 0 the Hamiltonian described the initial state and at k = 1, it described the final state of the process. Simulations were performed at a number of discreet k-values and the free energy was obtained by numerical integration of the derivative of the Hamiltonian with respect to k: where AE ¶H(r, p, k)/ ¶kae k represents an ensemble average obtained at the specified value of k, r and p are the positions and momenta of all particles in the system. Free energies of mutating residue 483 and of R/S inversion of propranolol were calculated as shown Fig. 2. Processes involving the inversion of the chiral center of propranolol were performed in a two-step process (Zhou et al. 2005(Zhou et al. , 2006 Fig. 2 Schemes describing the (reversible) inversion of R-into S-propranolol and the (reversible) mutation of the residue on position 483 in CYP2D6 (F (WT) into A (mutant)). The planar propranolol intermediate is indicated by ''I'', noninteracting dummy atoms in the A483 residue are indicated by '' du''. An asterisk indicates atom C5 of propranolol was changed from a CH2 group into a CH3 group. The nonbonded interactions at intermediate k-values were described using the soft-core approach to avoid singularities in the derivative (Beutler et al. 1994). No changes were made in the covalent interactions. Twenty ps of equilibration were followed by 200 ps of data collection at every k-value. Statistical error estimates on the ensemble averages were obtained using block averages (Allen and Tildesley 1987). Simulations with relatively large statistical errors were prolonged for another 200 ps, and additional points were added wherever needed to obtain a smooth profile of AE ¶H/ ¶kae as a function of k. A more empirical error estimate was obtained by comparing calculated free energy values associated with a forward process and the corresponding reverse process, also called hysteresis. An average free energy estimate from the forward and backward processes can be calculated using Boltzmann weighted averages using (Jarzinsky 1997;Lu et al. 2003): A total of 76 simulations (6 · 2 · 6 + 4 added kpoints) were performed involving R/S inversions and 72 simulations (6 · 11 + 6 added k-points) were performed involving F/A mutations. This leads to an overall simulation time of 33 ns.

Experimental results
Both the propranolol enantiomers displayed a type I binding spectrum when added to wild-type CYP2D6, as derived from a decrease of the Soret maximum at 417 nm, and the appearance of a shoulder at 380 nm, which is typical for type I binding of a substrate (Jefcoate 1978). The affinity of the wild-type enzyme for the R-and S-enantiomers was similar (Table 1). The addition of S-propranolol to the F483A mutant gave the same results as obtained for wild-type CYP2D6. The addition of R-propranolol to the F483A mutant of CYP2D6 did also lead to a type I spectrum, however, the affinity was 20-fold lower compared to the wild-type enzyme.

Computational approaches
Simulations at the initial and final stages along the solid arrows in Fig. 1 have been analyzed in detail. For every combination of R-or S-propranolol and WT or F483A mutant, four of such simulations of 200 ps were performed, labeled with small letters in Fig. 1. The four simulations included the end points of at least two different processes (R/S-inversion or F/A mutation) allowing for a comparison of convergence at every state.
The atom-positional root-mean-square deviation of C a atoms with respect to the initial homology model, remained stable throughout the simulations, with maximum values of 0.25 nm. For a 464 residue protein (the N-terminus of the homology model is truncated by 33 residues) such values are reasonable and indicate simulations of a stable protein structure.

Binding mode orientation
The active site of CYP2D6 is formed by substrate recognition sites (SRS) 1 (B-C loop), 2 (F-helix), 4 (I-helix), 5 (b 1-4 ) and 6 (b 4-1 ). Figure 3 displays a top view of the active site with interacting residues displayed in stick. The protonated nitrogen atom and the hydroxyl group of propranolol can form a hydrogen bond network with the negatively charged carboxylate of E216 (not visible in Fig. 3, but shown in Fig. 5). Two binding modes for the naphthalene moiety can be distinguished (see Fig. 3). In WT CYP2D6, propranolol occupies a hydrophobic cavity formed by F120, V308, V370, and F483 (binding mode 1). In the F483A mutant, a larger cavity becomes available and SRS5 slightly shifts. The naphthalene moiety of propranolol now occupies a cavity formed by F120, V374, V370 and A483 (binding mode 2). In order to distinguish these binding modes, an atomic distance and an angle are defined in Fig. 3. The distance between the heme NB atom and the C5 atom of propranolol will be larger in binding mode 2. The angle defined by the C n atom in F120, the C5 atom of propranolol and the C b in V370 will approach 180 in binding mode 2. All simulations originated from the WT protein and R-propranolol occupying binding mode 1. During the free energy calculations, the ligand was seen to move to binding mode 2. Table 2 presents the average values of the NB-C5 distance and the C n -C5-C b angle over the end-point All values are the means of at least two independent experiments; ±SD as described in the experimental procedures Eur Biophys J (2007) 36: 589-599 593 simulations, together with the average number of hydrogen bonds between propranolol and E216. The presence of a hydrogen bond was determined using geometric criteria: a donor-acceptor distance of at most 0.25 nm and a H-bond angle of at least 135°. Alternative hydrogen bonds partners of propranolol were not observed for significant periods of time.

Free energy calculations
The values of AE ¶H/ ¶kae as a function of k are displayed in Fig. 4 for all thermodynamic integration simulations. Values for the backward processes are multiplied by -1, for ease of comparison to the forward processes. Free energy estimates for all processes were obtained by a simple trapezoid integration of these curves and the resulting free energies are presented in Table 3. Using the definitions in Fig. 1, the different relative free energies are calculated and compared to experimental values.

Discussion
The primary aim of this study was to reproduce and rationalize the experimental finding that the F483A mutant of CYP2D6 has lower affinity for R-propranolol than for S-propranolol, using detailed molecular dynamics (MD) simulations. From these simulations, free energy estimates for all processes in Fig. 1 have Propranolol is depicted in yellow sticks and hydrophobic substrate binding amino acid residues are depicted in grey ball-and-sticks. Two geometric features used to discriminate between the two binding modes are displayed: the distance between the nitrogen atom of heme pyrrole ring B (NB) and the C5 atom of R-propranolol and the angle defined by the C n atom in F120, the C5 atom in R-propranolol and the C b in V370 been derived, which are presented in Table 3. The error estimates in Table 3 should be considered as statistical error estimates, indicating the statistical precision of the calculation. The absolute accuracy, however, is more difficult to determine. Three criteria that may give an indication of the accuracy can be distinguished. First, the hysteresis, defined as the difference between performing a process in the forward and backward directions, yields information on the convergence of the directions. At infinite sampling times, the hysteresis will be equal to zero. In Table 3, two cases can be observed with relatively large hystereses, the R/S inversion of propranolol in WT CYP2D6 and the F/A mutation in S-propranolol bound protein. In the other cases the hysteresis falls within the statistical uncertainty. Second, using the studied processes, one can obtain cyclic pathways along  Fig. 3 together with the average number of hydrogen bonds between propranolol and E216 for the different MD-simulations at the initial and final stages of the inversion and mutation processes defined in Fig. 1 The results belonging to approach 2 (defined in Fig. 1  which the sum of free energies should be equal to zero. In Table 3, such a pathway is presented as DG cycle : It is clear that the present calculations do not meet this criterion, with values of DG cycle larger than the statistical uncertainty. This indicates that not all free energy calculations have converged properly. The third criterion is the ability to reproduce experimental data. Knowing that not all simulations have reached convergence, this last criterion allows us to determine if any of the four approaches mentioned in the introduction can still be used to rationalize the experimental data at an atomic level.
Free energies of binding: approaches 1 and 3 In the first approach, the relative free energies of binding R-and S-propranolol are determined for both proteins. This approach uses a simulation of the R/Sinversion free in solution, of which the free energy value should fundamentally equal zero. As can be seen in Table 3, indeed DG free inv is calculated to be 0 within the error estimate. However, the relative free energies DDG WT bind ðR; SÞ and DDG F483A bind ðR; SÞ do not correspond to the experimentally determined value. The inversion process of propranolol when bound to the protein seems to be too unfavourable. Still, this may be due to a systematic contribution in the protein simulations, which would cancel in approach 2, where only free energies resulting from simulations of the propranolol inversion when bound to the protein are compared.
In the third approach, the relative free energies of binding R-and S-propranolol to the WT and mutant CYP2D6 protein are calculated and compared. The mutation in the ligand-free protein shows the smallest hysteresis in Table 3. The obtained value of DDG R bind ðWT; F483AÞ is within 4 kJ mol -1 of the experimental value, which is within the statistical uncertainty, but not satisfactory. The value of DDG S bind ðWT; F483AÞ is even further from the experimental value. Note that the value of DG S mut shows the largest hysteresis, indicating that particularly these simulations may not have converged properly. An indication of the convergence can be obtained by comparing simulations k and l to simulations w and x in Fig. 1 in terms of the structural properties. From the average distance and angle in Table 2, we can see that in simulations k and l, the naphthalene moiety of S-propranolol binds in binding mode 2, while in simulations w and x, Spropranolol still orients according to binding mode 1 (see Fig. 2). This difference indicates that the processes (R/S inversion in the mutant protein, or F/A mutation with S-propranolol bound) do not converge to the same end point. This is also reflected by significant differences in H-bond counts (see Table 3): S-propranolol forms more H-bonds with E216 in the F483A mutant in simulations k and l than it does in simulations w and x. A similar difference in the average number of hydrogen bonds between S-propranolol and E216 in WT is observed when comparing simulations g and h with simulations s and t. This indicates that S-propranolol occupies a different area in phase space during the F/A mutation than during the other simulations leading to the same end-points. With these observations, we also find an explanation for the non-zero value of DG cycle : Since not all end-points of the simulations are connected, the calculated cycle is not complete.
Free energies of transfer and exchange: approaches 2 and 4 Approach 4 uses the same mutational processes with propranolol bound to the protein as were used in approach 3. From these values the free energies of exchanging R-propranolol by S-propranolol in the mutant relative to WT (DDG exchange (WT,F483A)) can be obtained. As was discussed above, the F/A mutation with S-propranolol bound does not converge to the same end-points. For this reason it should not be expected that the value of DDG exchange (WT,F483A) will correspond to the experimental data (Table 3). The only approach for which the calculated free energies match the experimental data appears to be approach 2. In this approach, the free energy of transferring R-propranolol from WT CYP2D6 to the F483 mutant is calculated relative to the free energy of transferring S-propranolol ðDDG transfer ðR; SÞÞ: The obtained value lies only 0.2 kJ mol -1 from the experimental value, which is well within the statistical error estimate of the calculations. It is also within k B T (~2.5 kJ mol -1 ), corresponding to thermal fluctuations, which is the accuracy that can be expected from free energy calculations (Brandsdal et al. 2003;van Gunsteren et al. 2002).
In summary, the approaches using protein mutations (i.e., 3 and 4) do not converge sufficiently to reproduce the experimental data. A systematic error in the R/S inversion for protein-bound propranolol disqualifies approach 1, but still allows for a comparison with approach 2 (the inversion in the different proteins), which does yield free energy values in agreement with the experimental data. From Fig. 4 it can be seen that AE ¶H/ ¶kae takes much larger values in the processes involving a mutation than in the processes involving an inversion. It has been shown earlier that free energy calculations involving smaller structural changes and smaller energy barriers yield more accurate free energy estimates than integrations involving large changes in AE ¶H/ ¶kae (Mark et al. 1991;Zhou et al. 2005).
Possibly, better results with lower hysteresis and cycle closure, could be obtained by further extending the length of the simulations and by adding more k-points (Mark et al. 1991;Oostenbrink and Van Gunsteren 2006). Another possible source of error could be the fact that a (substrate-bound) homology model was used, which has slightly different properties than the recently solved substrate-free crystal structure. However, we conclude that the changes required in the active site of CYP2D6 upon ligand binding are so large (de Graaf et al. 2007;Rowland et al. 2005) that the use of the homology model is justified. In addition, CYP2D6 will remain a challenging target for computational approaches, because of the broad range of inhibitors and substrates and the multiple binding modes observed for these . Now that we have determined that approach 2 reproduces the experimental data best, it is of interest to investigate the molecular basis for the calculated free energies in this approach.

Structural rationalization
In order to rationalize the difference in affinity for Rand S-propranolol upon mutation of residue 483, we can compare the observed binding modes in simulations e -l in Fig. 1. Representative orientations of propranolol in the binding cavities of WT and mutant protein in these simulations are presented in Fig. 5. In all simulations propranolol is involved in a salt bridge between the protonated nitrogen atom and the negatively charged E216 in the protein. Together with the propranolol hydroxyl group, these three moieties form a dynamic hydrogen bond network. The average number of hydrogen bonds observed throughout the simulations, however varies for different combinations of R/S-propranolol and WT/F483A proteins. Already in WT (simulations e -h), we observed fewer hydrogen bonds for R-propranolol than for S-propranolol (see Table 2 and Fig. 5).
For both stereoisomers, the naphthalene moiety orients itself according to binding mode 1, it fits snugly into a hydrophobic pocket formed by V308, V370, and F483 (pocket 1, see Figs. 3,5). Upon the F483A mutation, the size of the pocket increases, which opposes tight hydrophobic interactions. The naphthalene moiety now occupies the pocket formed by F120, V370, V374 and A483 (pocket 2), according to binding mode 2, as indicated by a larger value of the heme-NB to propranolol-C5 distance (see Table 2). A large value of Eur Biophys J (2007) 36:589-599 597 the C n -C5-C b angle indicates that a position between F120 and V370 is occupied. In the F483A mutant, R-propranolol binds in pocket 2 without significant decrease in the average number of hydrogen bonds (simulations i and j), compared to the simulations performed in WT (e and f). When S-propranolol binds in pocket 2 in F483A, hydrogen bond interactions are lost, in particular the propranolol-OH to E216-O e hydrogen bond is observed only for 50% of the time in simulation l. This hydrogen bond can be regained by moving away from F120 and V370 (angle C n -C5-C b decreases) in simulation k. In fact the average number of hydrogen bonds in this simulation is even larger than for S-propranolol in WT, compensating the lack of tight hydrophobic interactions at the naphthalene moiety. R-propranolol cannot compensate for this loss. In summary, the following picture arises. R-propranolol forms fewer hydrogen bonds than S-propranolol, both in WT and in the F483A mutant. The mutation causes a loss of favourable hydrophobic interactions, which can be compensated by increased hydrogen bond formation by S-propranolol, but not by R-propranolol.

Conclusions
Detailed molecular dynamics (MD) simulations have been performed to reproduce and rationalize the experimental finding that the F483A mutant of CYP2D6 has lower affinity for R-propranolol than for S-propranolol. Wild-type (WT) CYP2D6 does not show this stereospecificity. Four different approaches to calculate the free energy differences have been investigated and were compared to the experimental binding data. From the differences between calculations based on forward and backward processes and the closure of the thermodynamic cycles, it was clear that not all simulations have converged sufficiently. The approach that calculates the free energies of exchanging R-propranolol with S-propranolol in the F483A mutant relative to the exchange free energy in WT accurately reproduced the experimental binding data. Careful inspection of the end-points of the simulations involved in this approach, allowed for a molecular interpretation of the observed differences.