Conformational Energy Landscape of the Ritonavir Molecule - The

Article pubs.acs.org/JPCB

Conformational Energy Landscape of the Ritonavir Molecule Debayan Chakraborty,† Neelanjana Sengupta,‡ and David J. Wales*,† †

Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom Physical Chemistry Division, CSIR-National Chemical Laboratory, Dr. Homi Bhaba Road, Pune 411008, India

‡

S Supporting Information *

ABSTRACT: Conformational polymorphism of ritonavir, a well-known pharmaceutical drug, is intricately linked to its efficacy in the treatment of acquired immunodeficiency syndrome (AIDS). Polymorphic transition from the crystalline form I to form II leads to the loss of bioactivity. The constituent ritonavir molecules adopt a trans configuration about the carbamate torsion angle in the form I crystal, and a cis configuration in the form II crystal. Investigating the energetics and mechanistic features of conformational transitions at the single molecule level is a key step toward decoding the complex features of the solid state polymorphism. In this work, we employ the energy landscape framework to investigate the conformational transitions of an isolated ritonavir molecule. The landscape is explored using discrete path sampling (DPS) and visualized in terms of disconnectivity graphs. We identify two distinct funnels corresponding to the two molecular forms that are identified by crystallography. The two regions can be reliably distinguished using the carbamate torsion angle, and the corresponding interconversion rates are predicted to follow Arrhenius behavior. The results provide mechanistic insight into pathways for cis ↔ trans interconversion at the molecular level and may also help in elucidating the polymorphic transitions in the crystal state.

■

INTRODUCTION Polymorphism, the existence of compounds in two or more crystalline states, is a major concern in the chemical and pharmaceutical industries.1,2 Polymorphs usually have distinctive physicochemical properties, and therefore transitions from one form to another can lead to large alterations in their material characteristics, such as solubility or bioactivity. These critical issues may be compounded by the rule of stages (Ostwald’s step rule),3,4 according to which the least stable polymorph is usually the first to emerge during the crystallization process; crystals may therefore be prone to reach more stable polymorphic conformations only on larger time scales. Lately, a number of experimental approaches including crystallization, thermal analysis, and nanoindentation have been employed to investigate small molecule polymorphism.5−7 However, while these methods can isolate and characterize individual polymorphs, they do not generally provide kinetic descriptions of the interconversion. Computational investigations, on the other hand, can potentially elicit molecular details of conformational states and the rates associated with polymorphic interconversions.8 The importance of drug polymorphism was underscored about two decades ago in the case of ritonavir, a protease inhibitor with antiviral properties that is effective as a drug in the treatment of acquired immunodeficiency syndrome (AIDS).9−12 It is noteworthy that only a single crystalline form of ritonavir (form I) was identified during the drug development stage. This polymorph was found to be effective when administered in a semisolid form. However, within about © XXXX American Chemical Society

two years of shelf life, the molecule was found to fail dissolution tests because it transforms to a new crystalline form (form II). This polymorphic transformation severely hindered its bioavailability and the drug had to be eventually withdrawn as a treatment option. Recent crystallization experiments13 have reported additional polymorphic forms of ritonavir, and suggested plausible schemes for achieving “morphology control”. Despite the recent upsurge of interest in crystal energy landscapes,14−18 to the best of our knowledge there are no systematic studies that investigate the conformational polymorphism of ritonavir at the atomistic level. It is worth noting that while most of the existing crystal structure prediction (CSP) methods14−18 aim to identify low energy conformers of a crystal, they do not generally provide comprehensive mechanistic insights into the interconversion between thermodynamically and kinetically relevant conformers. In the case of ritonavir, the two crystal conformers have different values of the carbamate torsion angle with the constituent molecules adopting a trans configuration in the form I crystal and a cis configuration in the form II crystal.10,11 The connection between molecular conformation and crystal structure19,20 suggests that a detailed understanding of the energetics and structural transitions at the level of a single molecule could be very useful. In the present study, we have Received: December 15, 2015 Revised: April 13, 2016

A

DOI: 10.1021/acs.jpcb.5b12272 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

refinement schemes available within PATHSAMPLE, to locate shorter pathways as well as lower energy barriers, and to remove artificial frustration caused by undersampling.46 DPS runs were deemed to have converged when the rate constants corresponding to the transition between the cis and trans conformers were stable to within an order of magnitude. Dijkstra’s shortest path algorithm47 with suitable edge weighting35 was used to extract the most kinetically relevant pathway connecting the two conformers. The stationary point databases from DPS were employed to estimate the free energies at 298 K using the harmonic superposition approximation.48,49 Here, the molecular partition function is expressed as a sum over the harmonic vibrational densities of states corresponding to the catchment basins of each local minimum. The rate constants for the cis ↔ trans interconversion were estimated at 298 K using the new graph transformation (NGT) method50 in conjunction with a selfconsistent lumping scheme based on free energy barriers.36 The potential and free energy landscapes of ritonavir were visualized using disconnectivity graphs.51−54 This representation is simple but powerful and provides a faithful description of the kinetics by preserving the barriers between different local minima.55 The transition state corresponding to the rate-determining step in the interconversion pathway was reoptimized using density functional theory (DFT), employing the B3LYP exchange-correlation functional56,57 and a 6-31G(d) basis set. The adjoining minima were found by perturbing the transition state geometry in the directions parallel and antiparallel to the eigenvector corresponding to the unique negative eigenvalue. To obtain refined estimates of the barrier heights corresponding to the conformational transition, single point energy calculations in vacuum were carried out on the optimized geometries using the B3LYP56,57 and the MPW91PW9158 exchange-correlation functionals in conjunction with the 6311G++(d,p) basis set. The transition state and the corresponding minima were also reoptimized using dispersion corrected DFT (DFT-D3)59 at the B3LYP/6-31G(d) level to obtain a better description of the aromatic stacking interactions. The dispersion effects were described using version D3 of Grimme et al.’s corrections incorporating Becke−Johnson damping.59,60 The DFT-D3 calculations were performed using the Gaussian 09 code.61 Selected minima from the AMBER database were reoptimized at the B3LYP/6-31G(d) level. All DFT calculations were performed using the Gaussian 0362 package interfaced to OPTIM. These minima were also reoptimized at the DFT-D3 level employing the same combination of functional and basis set as the DFT calculations.

characterized the dominant interconversion pathway between the trans and cis conformers of an isolated ritonavir molecule to seek insight into the complexities of crystal polymorphism. The slow kinetics associated with the cis ↔ trans transition suggest that the corresponding interconversion is unlikely to be captured by unbiased dynamical sampling. However, such “rare events” are ideal targets for our computational energy landscape based approach,21 which is practically independent of barrier heights and time scales. We use discrete path sampling (DPS), a simulation framework based on geometry optimization techniques, to locate stationary points (minima and transition states) on the underlying potential energy landscape for the ritonavir molecule. A kinetic transition network22,23 is constructed from databases of stationary points to provide mechanistic and kinetic insight into the conformational transitions. We find that the potential and free energy landscapes exhibit two principal funnels, corresponding to the cis and trans conformers. The interconversion pathway is consistent with a two-state kinetic picture. In the subsequent sections, our findings are discussed in the context of previous experimental observations of ritonavir polymorphism.

■

COMPUTATIONAL METHODOLOGY The atomic coordinates of trans and cis conformers for the ritonavir molecule are unpublished and were provided to us by Dr. Edward Pyzer-Knapp. The corresponding CSD entries are YIGPIO and YIGPIO01, respectively.10,24 The conformations were modeled using generalized AMBER force-field (GAFF) parameters.25 Solvent effects were accounted for implicitly via the generalized Born model.26 A preliminary exploration of the energy landscape was carried out using the basin-hopping27,28 global optimization algorithm, as implemented in our GMIN code29 interfaced to AMBER9.30 Separate basin-hopping runs were started from both the conformers. In previous work, basin-hopping has been successfully applied in the study of clusters,27 glass formers,31 and biomolecules.32 The algorithm is particularly efficient in locating global minima for landscapes featuring broken ergodicity. Group rotation moves (rotating a group of atoms as a rigid body about a predefined axis) in conjunction with random Cartesian displacements were employed to perturb the current local minimum in the chain. Further exploration of the energy landscape was carried out using discrete path sampling (DPS).33,34 Previously, the DPS procedure has been successfully employed to investigate energy landscapes for a wide range of systems.35−38 In DPS, the objective is to construct a kinetic transition network22,23 in terms of the stationary points on the underlying energy landscape, by mapping out discrete paths between local minima of interest. A discrete path is a sequence of local minima linked by intervening transition states.33 DPS runs were used to connect different pairs of minima located by the basin-hopping. A modified version of the LBFGS algorithm39 was used for all the local minimizations. Initial guesses for transition states between pairs of local minima were found using the doublynudged40 elastic band41,42 method. These transition state candidates were then tightly converged to a root-meansquare-gradient of 10 −6 kcal mol −1 Å −1 using hybrid eigenvector-following.43 The OPTIM code44 interfaced with AMBER9 30 was used to carry out all the geometry optimizations and transition state searches. The PATHSAMPLE code45 was used to distribute parallel OPTIM jobs across compute nodes to connect different pairs of minima. The databases were further expanded using various

■

RESULTS AND DISCUSSION The potential energy and free energy disconnectivity graphs (at 298 K) for ritonavir are shown in Figures 1 and 2, respectively. The graphs share the same topological features. Each energy landscape consists of two major funnels leading to competing structures. The branches leading to the minima are colored according to the values of the carbamate torsion angle. This specific structural order parameter results in a precise segregation of the two major basins and is therefore an excellent diagnostic for the ritonavir landscape. We have confirmed this attribute by constructing disconnectivity graphs where the branches are colored according to committor probabilities63,64 for visiting the global minimum. These graphs B


Article


Figure 1. Potential energy landscape of ritonavir. The branches are colored according the value of the carbamate torsion angle with blue representing structures corresponding to the trans conformer and red representing structures corresponding to the cis conformer. Some representative structures are also shown.

(see Supporting Information, Figures S1 and S2) are identical to those in Figures 1 and 2. The remarkable agreement is interesting because structural order parameters often fail to fully capture the complexity of the underlying landscape and may merge together regions that are actually well separated in terms of kinetics.65,66 In the disconnectivity graph, minima corresponding to the molecule in the trans conformation constitute the blue funnel, whereas those corresponding to the cis conformation are located at the bottom of the red funnel. The global minimum for the landscape is a trans conformation where the two thiazoyl rings are stacked on top of each other (labeled as f in Figures 1 and 2). In addition, there are subfunnels within the blue section of the disconnectivity graphs, comprising trans conformations in which the thiazoyl rings are partially unstacked (labeled as d in Figure 1, and e in Figure 2), or completely unstacked (labeled as e in Figure 1 and d in Figure 2). In contrast, the red section of the graph has no significant subfunnels. The lowest energy cis structure (labeled as b in Figures 1 and 2) also exhibits stacking of the thiazoyl rings. Interestingly, the conformations located at the bottom of the major funnels have compact structures that are incompatible with the extensive intermo-

lecular hydrogen bonding observed in the crystal state. Conformations that are structurally closest to the crystal conformers (as determined by the all-atom RMSD) have extended structures and are relatively higher in energy. This trend is clearly visible in the map depicting the potential energy of the minima as a function of the RMSD from the cis and trans crystal conformers (Supporting Information, Figure S3). The nodes corresponding to the cis and trans minima that exhibit the lowest RMSD from the crystal conformers are labeled as a and c, respectively, in the disconnectivity graphs. Conformation a is approximately 14 kcal/mol higher in potential energy relative to the lowest energy minimum in the cis basin, and conformation c is approximately 13 kcal/mol higher in potential energy relative to the lowest energy minimum in the trans basin. The cis and trans conformers are separated by high potential and free energy barriers, making the landscape topologically frustrated. To obtain mechanistic insight into conformational switching, the “fastest pathway” determined as the largest contribution to the phenomenological rate constant (Figure 3) between two low-lying potential energy minima in the cis and trans basins was extracted from the stationary point databases. The rate-determining step for the transition from the cis to the C


Article


Figure 2. Free energy landscape of ritonavir at 298 K. The coloring scheme is the same as in Figure 1. The corresponding structures from the potential energy landscape are also shown.

The DFT minima are structurally similar to their AMBER counterparts, except for some loss of stacking interactions between the thiazoyl rings. The transition state has a carbamate torsion angle of 111°, which is slightly higher than the corresponding AMBER value. In addition, the DFT transition state also exhibits unstacking of the thiazoyl rings. The aromatic stacking interactions in the transition state and the corresponding minima are stabilized when the reoptimizations are carried out at the DFT-D3 level (Supporting Information, Figure S4). This observation is consistent with previous studies,67−69 which indicate that inclusion of dispersion corrections is essential to provide a proper description of noncovalent interactions. In fact, the DFT-D3 stationary points correspond quantitatively with their AMBER counterparts, with the transition state exhibiting a carbamate torsion angle of 104° (Supporting Information, Figure S5). The barrier heights computed at different levels of theory are tabulated in Table 1. The inclusion of extra polarization and diffuse functions decreases ΔE2 by approximately 3.2 kcal/mol, whereas there is only a minimal effect on ΔE1. This trend is

trans conformer corresponds to a potential energy barrier of 13.4 kcal/mol. For the reverse transformation, the barrier is 16.6 kcal/mol. This particular transition state exhibits a carbamate torsion angle of 104°, which is approximately midway between the cis and trans configurations. The rest of the pathway is mostly characterized by minor rearrangements of atoms near the thiazoyl and phenyl moieties. Interestingly, the carbamate torsion angle does not change continuously along the pathway but instead switches from around 0° to around 180° via a single jump near the rate-determining step, as shown in Figure 3. To confirm the nature of the rate-determining step, the corresponding transition state was reoptimized at the DFT level, and the adjoining minima were located using the technique described in the Computational Methodology section. We find that approximate steepest-descent paths from the transition state lead to cis and trans minima, respectively. The potential energy profile for the conformational transformation computed at the B3LYP/6-31G(d) level is shown in Figure 4. D


Article


Figure 3. Top panel: The evolution of the total potential energy (V) with integrated path length (s) for the fastest pathway between the cis and trans conformers. Some representative structures encountered at different stages of the conformational transition are also shown. Bottom panel: The variation of the carbamate torsion angle along the pathway. There is a rapid jump from cis (∼0°) to trans (∼180°) via a single transition state.

energy trans conformer (LMtrans AMBER) in the AMBER database. These two minima are separated by a relatively low barrier (around 1.5 kcal/mol), and after recursive grouping36 of the stationary points in the database based on free energy barriers at 298 K they become members of the same macrostate. The lowest energy cis conformer (LMcis DFT) found from the DFT calculations can be mapped back to a local minimum on the AMBER landscape, which is approximately 4.1 kcal/mol higher in energy than the lowest energy minimum in the cis funnel cis (LMcis AMBER). We find that LMDFT is identified as the secondlowest minimum after reoptimization and is approximately 0.3 kcal/mol higher in energy relative to LMtrans DFT, which is in contrast to the energy difference of about 3.6 kcal/mol between trans LMcis AMBER and LMAMBER. The structures corresponding to trans cis cis LMDFT and LMDFT, superimposed on LMtrans AMBER and LMAMBER respectively, are shown in Figure 5c,d. As illustrated in the trans has a more open structure compared to figure, LMDFT cis LMtrans . Both LMtrans AMBER DFT and LMDFT exhibit loss of stacking interaction between the thiazoyl rings. The lowest energy DFT conformers are also shown superimposed on their corresponding minima from the AMBER database in Figure 5e,f, respectively.

consistent among the different functionals and basis sets. The inclusion of dispersion interactions also has a more pronounced effect on ΔE2 compared to ΔE1. In addition, the difference between the ΔE2 and ΔE1 values predicted by DFT-D3 is closer to the AMBER estimate than DFT. The estimated barrier heights from DFT and DFT-D3 indicate that the energy landscape consists of two major basins corresponding to cis and trans conformers, in line with the predictions using the AMBER force-field. Minima from the kinetic transition network that are connected to at least two transition states and lie at the bottom of monotonic sequences70,71 were also reoptimized at the B3LYP/6-31G(d) and B3LYPD3/6-31G(d) levels without inclusion of solvent effects. This choice is particularly effective as the resultant subnetwork preserves the organization of the landscape (Supporting Information, Figure S6), thereby ensuring that the selected minima are thermodynamically as well as kinetically relevant. The lowest energy trans conformer (LMtrans DFT) identified by DFT corresponds to a low-lying local minimum on the AMBER potential energy surface. This local minimum is destabilized by only 0.8 kcal/mol with respect to the lowest E


Article


Figure 4. Relative potential energy as a function of the path length for the rate-determining step, computed at the B3LYP/6-31G(d) level. Positive and negative values of the path length indicate different directions from the transition state, toward the reactant and product minimum, respectively.

Table 1. Barrier Heights Corresponding to the Cis ↔ Trans Transition, Computed at Different Levels of Theorya

a

level of theory

ΔE1(cis→trans)

ΔE2(trans→cis)

AMBER B3LYP/6-31G(d) B3LYPD3/6-31G(d) B3LYP/6-311++G(d,p) MPW91PW91/6-311++G(d,p)

13.40 17.11 18.95 16.65 17.08

16.60 25.33 21.50 22.15 22.15

All energies are in kcal/mol.

Reoptimizations with DFT-D3 indicate that a poor description of dispersion interactions in conventional DFT leads to destabilization of the stacking interactions between the thiazoyl rings. As shown in Figure 6a, in the lowest energy trans conformer (LMtrans DFT−D3) identified by DFT-D3 the thiazoyl rings exhibit perfect stacking. This structure corresponds to a low-lying local minimum on the AMBER potential energy surface, which is destabilized by approximately 1 kcal/mol with respect to LMtrans AMBER. The aromatic stacking interactions are also maintained in the lowest energy cis conformer (LMcis DFT−D3), as shown in Figure 6b. This particular minimum can be mapped back to an AMBER local minimum, which is approximately 0.3 kcal/mol higher in energy than LMcis AMBER. The energy cis difference between LMtrans and LM DFT−D3 DFT−D3 is 3 kcal/mol, which is much closer to the corresponding AMBER estimate (3.6 kcal/mol) than DFT. Interestingly, both LMtrans DFT−D3 and LMcis DFT−D3 map onto minima on the AMBER landscape that are lumped together into the same macrostate with LMtrans AMBER and cis LMAMBER respectively, after recursive regrouping.36 The cis structures corresponding to LMtrans DFT−D3 and LMDFT−D3, supercis imposed on LMtrans and LM , respectively, are shown in AMBER AMBER Figure 6c,d respectively. Figures 6e,f illustrate the lowest energy DFT-D3 conformers superimposed on their corresponding minima from the AMBER database. The relative energies of local minima with respect to the global minimum are within the same range of values for AMBER, DFT, and DFT-D3. As shown in Figure 7, the distribution progressively shifts toward the low energy side on

Figure 5. (a) The lowest energy trans conformer, LMtrans DFT, identified by DFT. (b) The lowest energy cis conformer, LMcis DFT, identified by DFT. (c,e) Panels show the trans conformer, LMtrans DFT, superimposed on LMtrans AMBER, and its mapped AMBER minimum, respectively. (d, f) cis Panels show the cis conformer, LMcis DFT, superimposed on LMAMBER, and its mapped AMBER minimum, respectively.

going from DFT-D3 to AMBER. This trend is in agreement with previous work on cyclic peptides32 where the energy gaps predicted by AMBER were found to be lower than those from DFT. Significant reordering of the lowest energy structures occurs after reoptimization with DFT and DFT-D3. This effect is illustrated in Figure 8. Each AMBER minimum was assigned a rank based on its energy (with lower energy having a higher rank) and the change in its ranking, ΔR AMBER→DFT F


Article


Figure 6. (a) The lowest energy trans conformer, LMtrans DFT−D3, identified by DFT-D3. (b) The lowest energy cis conformer, LMcis DFT−D3, identified by DFT-D3. (c,e) Panels show the trans conformer, trans LMtrans DFT−D3, superimposed on LMAMBER and its mapped AMBER minimum, respectively. (d,f) Panels show the cis conformer, cis LMcis DFT−D3, superimposed on LMAMBER and its mapped AMBER minimum, respectively.

(ΔRAMBER→DFT‑D3), after reoptimization with DFT (DFT-D3). Negative values of ΔRAMBER→DFT (ΔRAMBER→DFT‑D3) in Figure 8 indicate a more favorable ranking after reoptimization, whereas positive values indicate destabilization with respect to other minima. In Figure 8, the blue and red bars correspond to trans and cis minima, respectively. Most of the cis conformers are stabilized with respect to other minima after reoptimization with DFT as well as DFT-D3. In fact, in the AMBER database, the lowest energy cis conformer (LMcis AMBER) is assigned a rank of 50, whereas after reoptimization with DFT 9 cis conformers ranked within the top 50. Eight cis conformers are assigned ranks within the top 50 after reoptimization with DFT-D3. However, compared to DFT, DFT-D3 supports a larger number of cis conformers within the top hundred. Twenty-one

Figure 8. Change in ranking with respect to other minima after reoptimization with DFT-D3 (top) and DFT (bottom). A positive value indicates that the structure is destabilized with respect to other minima, causing a shift toward lower ranks. A negative value indicates a more favorable ranking after reoptimization. The red and blue bars correspond to cis and trans minima, respectively.

Figure 7. Distribution of energy gaps relative to the global minimum predicted by DFT-D3 (top panel), DFT (middle panel), and AMBER (bottom panel). G



■

cis conformers are ranked within the top hundred after reoptimization with DFT-D3, compared to 15 with DFT. A complete mapping of the energy landscape at the DFT or DFT-D3 level will require more extensive sampling of the underlying potential energy surface, and is beyond the scope of the current work. Nonetheless, we believe that the doublefunnel nature of the landscape will be preserved, and the key mechanistic and kinetic details will not change. The slow transition between the two conformers was quantified by estimating the corresponding rate constants at 298 K, employing the self-consistent recursive lumping scheme.36 This procedure alleviates bias arising due to the original choice of end points. The rate constant estimates are 64 s−1 for cis → trans and 4.8 × 10−3 s−1 for the reverse transformation. A lumping threshold of 3 kcal/mol proved to be appropriate, as quantitatively similar results were obtained in terms of rate constants and disconnectivity graphs for a range of values around this threshold. The activation energy for the cis to trans transformation was calculated using the Arrhenius equation, employing rate constants in the temperature range of 280 to 400 K obtained using the NGT procedure.50 For each temperature, an appropriate lumping threshold was chosen to converge the rate constant. The Arrhenius expression is quantitatively accurate (Supporting Information, Figure S7) with an activation energy of 14 kcal/mol and a preexponential factor of 1. 35 × 1012 s−1. The estimated activation energy is in excellent agreement with the barrier height corresponding to the rate-determining step for the cis to trans conformational change. The appearance of the underlying potential and free energy landscapes, as well as the mechanistic and kinetic analysis, indicate that the cis ↔ trans interconversion is a two-state kinetic process. Furthermore, relaxation from a high energy minimum is much more likely to lead to the trans form due to its greater entropy. However, once kinetic trapping occurs in the cis funnel, subsequent equilibration between the conformations is slow, leading to broken ergodicity. This separation of time scales is a characteristic property of multifunnel potential energy landscapes.72 Our observation agrees with previous studies, which reported that the form I polymorph (with a trans configuration about the carbamate torsion), suddenly “disappeared” from the solution, and all attempts to synthesize it were unsuccessful.9,11,73

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b12272. Additional data related to this publication is available at the Cambridge Data Repository. The appropriate link is https:// www.repository.cam.ac.uk/handle/1810/252982. Figures showing potential and free energy disconnectivity graphs colored according to committor probabilities; potential energy map depicting the energy of the minima as a function of their RMSD from the cis and trans crystal conformers; stationary points reoptimized at the DFT-D3 level; DFT-D3 stationary points overlaid on the corresponding AMBER counterparts; potential energy disconnectivity graph for the minima connected to at least two transition states and defining monotonic sequences; Arrhenius plot showing the variation of rate constants with temperature. (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge Dr. Edward Pyzer-Knapp for providing us with the atomic coordinates of the ritonavir molecule. We are also grateful to Mr. Kyle Sutherland-Cash, Dr. Chris Whittleston, Dr. Cheng Shang, and Dr. Joanne Carr for fruitful discussions. We thank the ERC for financial support. D.C. acknowledges the Cambridge Commonwealth, European and International Trust for financial support. N.S. thanks Council for Scientific and Industrial Research, Government of India, for her 2014-15 Raman Research Fellowship award.

■

REFERENCES

(1) Hilfiker, R. Polymorphism: in the pharmaceutical industry; John Wiley & Sons: New York, 2006. (2) Desiraju, G. R. Polymorphism: The Same and Not Quite the Same. Cryst. Growth Des. 2008, 8, 3−5. (3) Ostwald, W. Studien uber die Bildung und Umwandlung fester Korper. Z. Phys. Chem. 1897, 22, 289−330. (4) Ng, J. D.; Lorber, B.; Witz, J.; Théobald-Dietrich, A.; Kern, D.; Giegé, R. The crystallization of biological macromolecules from precipitates: evidence for Ostwald ripening. J. Cryst. Growth 1996, 168, 50−62. (5) Hamilton, B. D.; Ha, J.-M.; Hillmyer, M. A.; Ward, M. D. Manipulating Crystal Growth and Polymorphism by Confinement in Nanoscale Crystallization Chambers. Acc. Chem. Res. 2012, 45, 414− 423. (6) Giron, D. Applications of Thermal Analysis and Coupled Techniques in Pharmaceutical Industry. J. Therm. Anal. Calorim. 2002, 68, 335−357. (7) Varughese, S.; Kiran, M. S. R. N.; Solanko, K. A.; Bond, A. D.; Ramamurty, U.; Desiraju, G. R. Interaction anisotropy and shear instability of aspirin polymorphs established by nanoindentation. Chem. Sci. 2011, 2, 2236−2242. (8) Sanderson, K. Model predicts structure of crystals. Nature 2007, 450, 771−771. (9) Bauer, J.; Spanton, S.; Henry, R.; Quick, J.; Dziki, W.; Porter, W.; Morris, J. Ritonavir: an extraordinary example of conformational polymorphism. Pharm. Res. 2001, 18, 859−866.

■

CONCLUSIONS We have investigated the mechanistic and kinetic details of the conformational transitions for the ritonavir molecule in terms of its underlying energy landscape. Discrete path sampling simulations were used to explore the landscapes, which were visualized in terms of disconnectivity graphs. The energy landscape has double-funnel character and can be accurately segregated in terms of the carbamate torsion angle. We find that the cis and trans conformers are separated by a high energy barrier and interconvert slowly following Arrhenius kinetics. In the present study, we have not considered intermolecular interactions between the different monomers that are present in the crystal. It would be interesting to evaluate how these additional noncovalent as well as hydrogen-bonding interactions reshape the landscape and influence specific energetic barriers. These effects will be systematically investigated in future work. H


Article

The Journal of Physical Chemistry B (10) Galek, P. T. A.; Allen, F. H.; Fabian, L.; Feeder, N. Knowledgebased H-bond prediction to aid experimental polymorph screening. CrystEngComm 2009, 11, 2634−2639. (11) Nangia, A. Conformational Polymorphs, Multiple Z’Crystal Structures and Phase Transformations. J. Indian Inst. Sci. 2012, 87, 133. (12) Chemburkar, S. R.; Bauer, J.; Deming, K.; Spiwek, H.; Patel, K.; Morris, J.; Henry, R.; Spanton, S.; Dziki, W.; Porter, W.; Quick, J.; Bauer, P.; Donaubauer, J.; Narayanan, B.; Soldani, M.; Riley, D.; McFarland, K. Dealing with the Impact of Ritonavir Polymorphs on the Late Stages of Bulk Drug Process Development. Org. Process Res. Dev. 2000, 4, 413−417. (13) Morissette, S. L.; Soukasene, S.; Levinson, D.; Cima, M. J.; Almarsson, Ö Elucidation of crystal form diversity of the HIV protease inhibitor ritonavir by high-throughput crystallization. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 2180−2184. (14) Price, S. L. Predicting crystal structures of organic compounds. Chem. Soc. Rev. 2014, 43, 2098−2011. (15) Price, S. L. Computed Crystal Energy Landscapes for Understanding and Predicting Organic Crystal Structures and Polymorphism. Acc. Chem. Res. 2009, 42, 117−126. (16) Ismail, S. Z.; Anderton, C. L.; Copley, R. C. B.; Price, L. S.; Price, S. L. Evaluating a Crystal Energy Landscape in the Context of Industrial Polymorph Screening. Cryst. Growth Des. 2013, 13, 2396− 2406. (17) Kendrick, J.; Stephenson, G. A.; Neumann, M. A.; Leusen, F. J. J. Crystal Structure Prediction of a Flexible Molecule of Pharmaceutical Interest with Unusual Polymorphic Behavior. Cryst. Growth Des. 2013, 13, 581−589. (18) Price, S. L. From crystal structure prediction to polymorph prediction: interpreting the crystal energy landscape. Phys. Chem. Chem. Phys. 2008, 10, 1996−2009. (19) Cruz-Cabeza, A. J.; Bernstein, J. Conformational Polymorphism. Chem. Rev. 2014, 114, 2170−2191. (20) Bernstein, J.; Hagler, A. Conformational Polymorphism. The Influence of Crystal Structure on Molecular Conformation. J. Am. Chem. Soc. 1978, 100, 673−681. (21) Wales, D. J. Energy Landscapes; Cambridge University Press: London, 2003. (22) Wales, D. J. Energy landscapes: some new horizons. Curr. Opin. Struct. Biol. 2010, 20, 3−10. (23) Noe, F.; Fischer, S. Transition networks for modelling the kinetics of conformational change in macromolecules. Curr. Opin. Struct. Biol. 2008, 18, 154−162. (24) Bauer, J.; Spanton, S.; Henry, R.; Quick, J.; Dziki, W.; Porter, W.; Morris, J. Ritonavir: An Extraordinary Example of Conformational Polymorphism. Pharm. Res. 2001, 18, 859−866. (25) 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. (26) Onufriev, A.; Bashford, D.; Case, D. A. Modification of the Generalized Born Model Suitable for Macromolecules. J. Phys. Chem. B 2000, 104, 3712−3720. (27) Doye, J. P. K.; Wales, D. J. Global optimization by basinhopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. J. Phys. Chem. A 1997, 101, 5111−5116. (28) Wales, D. J.; Scheraga, H. Global Optimization of Clusters, Crystals, and Biomolecules. Science 1999, 285, 1368−1372. (29) Wales, D. J. GMIN: A program for finding global minima and calculating thermodynamic properties; http://www-wales.ch.cam.ac.uk/ software.html. (30) Case, D. A.; Darden, T. A.; Cheatham, T.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Goetz, A. W.; Kolossváry, I. AMBER 9, 2006. http://ambermd.org/. (31) Middleton, T. F.; Hernández-Rojas, J.; Mortenson, P. N.; Wales, D. J. Crystals of binary Lennard-Jones solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 64, 1842011−1842017.

(32) Oakley, M. T.; Johnston, R. L. Exploring the Energy Landscapes of Cyclic Tetrapeptides with Discrete Path Sampling. J. Chem. Theory Comput. 2013, 9, 650−657. (33) Wales, D. J. Discrete Path Sampling. Mol. Phys. 2002, 100, 3285−3305. (34) Wales, D. J. Some Further Applications of Discrete Path Sampling to Cluster Isomerization. Mol. Phys. 2004, 102, 891−908. (35) Evans, D. E.; Wales, D. J. Folding of the GB1 Hairpin Peptide from Discrete Path Sampling. J. Chem. Phys. 2004, 121, 1080−1090. (36) Carr, J. M.; Wales, D. J. Folding Pathways and Rates for the Three-Stranded beta-sheet Peptide Beta3s Using Discrete Path Sampling. J. Phys. Chem. B 2008, 112, 8760−8769. (37) Farrell, J. D.; Lines, C.; Shepherd, J. J.; Chakrabarti, D.; Miller, M. A.; Wales, D. J. Energy landscapes, structural topologies and rearrangement mechanisms in clusters of dipolar particles. Soft Matter 2013, 9, 5407−5416. (38) Chakraborty, D.; Collepardo-Guevara, R.; Wales, D. J. Energy Landscapes, Folding Mechanisms, and Kinetics of RNA Tetraloop Hairpins. J. Am. Chem. Soc. 2014, 136, 18052−18061. (39) Liu, D.; Nocedal, J. On the Limited Memory Method for Large Scale Optimization. Math. Program. 1989, 45, 503−528. (40) Trygubenko, S. A.; Wales, D. J. A Doubly Nudged Elastic Band Method for Finding Transition States. J. Chem. Phys. 2004, 120, 2082− 2094. (41) Henkelman, G.; Jònsson, H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys. 1999, 111, 7010−7022. (42) Henkelman, G.; Uberuaga, B. P.; Jònsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 2000, 113, 9901−9904. (43) Munro, L. J.; Wales, D. J. Defect Migration in Crystalline Silicon. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 3969− 3980. (44) Wales, D. J. OPTIM: A program for optimizing geometries and calculating pathways; http://www-wales.ch.cam.ac.uk/software.html. (45) Wales, D. J. PATHSAMPLE: A program for generating connected stationary point databases and extracting global kinetics; http://wwwwales.ch.cam.ac.uk/software.html. (46) Strodel, B.; Whittleston, C. W.; Wales, D. J. Thermodynamics and Kinetics of Aggregation for the GNNQQNY Peptide. J. Am. Chem. Soc. 2007, 129, 16005−16014. (47) Dijkstra, E. W. A Note on Two Problems in Connexion with Graphs. Numer. Math. 1959, 1, 269−271. (48) Hoare, M. R.; McInnes, J. J. Statistical mechanics and morphology of very small atomic clusters. Faraday Discuss. Chem. Soc. 1976, 61, 12−24. (49) Strodel, B.; Wales, D. J. Free energy surfaces from an extended harmonic superposition approach and kinetics for alanine dipeptide. Chem. Phys. Lett. 2008, 466, 105−115. (50) Wales, D. J. Calculating Rate Constants and Committor Probabilities for Transition Networks by Graph Transformation. J. Chem. Phys. 2009, 130, 2041111−2041117. (51) Becker, O. M.; Karplus, M. The topology of multidimensional potential energy surfaces: theory and application to peptide structure and kinetics. J. Chem. Phys. 1997, 106, 1495−1517. (52) Wales, D. J.; Miller, M. A.; Walsh, T. R. Archetypal Energy Landscapes. Nature 1998, 394, 758−760. (53) Krivov, S. V.; Karplus, M. Free energy disconnectivity graphs: Application to peptide models. J. Chem. Phys. 2002, 117, 10894− 10903. (54) Evans, D. A.; Wales, D. J. Free energy landscapes of model peptides and proteins. J. Chem. Phys. 2003, 118, 3891−3897. (55) Krivov, S. V.; Karplus, M. Hidden complexity of free energy surfaces for peptide (protein) folding. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 14766−14770. (56) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (57) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism I


Article

The Journal of Physical Chemistry B Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (58) Adamo, C.; Barone, V. Exchange functionals with improved long-range behavior and adiabatic connection methods without adjustable parameters: The mPW and mPW1PW models. J. Chem. Phys. 1998, 108, 664−675. (59) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104−154104−20. (60) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456−1465. (61) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B. et al. Gaussian 09; Gaussian Inc.: Wallingford, CT, 2009. (62) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J. et al. Gaussian 03; Gaussian, Inc.: Wallingford, CT, 2004. (63) Onsager, L. Initial Recombination of Ions. Phys. Rev. 1938, 54, 554−557. (64) Du, R.; Pande, V. S.; Grosberg, A. Y.; Tanaka, T.; Shakhnovich, E. S. On the transition coordinate for protein folding. J. Chem. Phys. 1998, 108, 334−350. (65) Wales, D. J. Perspective: Insight into reaction coordinates and dynamics from the potential energy landscape. J. Chem. Phys. 2015, 142, 130901. (66) Krivov, S. V.; Karplus, M. One-dimensional Free-Energy Profiles of Complex Systems: Progress Variables that Preserve the Barriers. J. Phys. Chem. B 2006, 110, 12689−12698. (67) Moellmann, J.; Grimme, S. DFT-D3 Study of Some Molecular Crystals. J. Phys. Chem. C 2014, 118, 7615−7621. (68) Mukherjee, S.; Kailasam, S.; Bansal, M.; Bhattacharyya, D. Energy hyperspace for stacking interaction in AU/AU dinucleotide step: Dispersion-corrected density functional theory study. Biopolymers 2014, 101, 107−120. (69) Huber, R. G.; Margreiter, M. A.; Fuchs, J. E.; von Grafenstein, S.; Tautermann, C. S.; Liedl, K. R.; Fox, T. Heteroaromatic π-Stacking Energy Landscapes. J. Chem. Inf. Model. 2014, 54, 1371−1379. (70) Kunz, R. E.; Berry, R. S. Statistical interpretation of topographies and dynamics of multidimensional potentials. J. Chem. Phys. 1995, 103, 1904−1912. (71) Berry, R.; Breitengraser-Kunz, R. Topography and Dynamics of Multidimensional Interatomic Potential Surfaces. Phys. Rev. Lett. 1995, 74, 3951−3954. (72) Doye, J. P. K.; Miller, M. A.; Wales, D. J. The double-funnel energy landscape of the 38-atom Lennard-Jones cluster. J. Chem. Phys. 1999, 110, 6896−6906. (73) Bucar, D.-K.; Lancaster, R. W.; Bernstein, J. Disappearing Polymorphs Revisited. Angew. Chem., Int. Ed. 2015, 54, 6972−6993.

J


Conformational Energy Landscape of the Ritonavir Molecule - The

Recommend Documents