Efficient Calculation of Microscopic Dissolution Rate Constants: The

Letter pubs.acs.org/JPCL

Efficient Calculation of Microscopic Dissolution Rate Constants: The Aspirin−Water Interface Julian Schneider and Karsten Reuter* Chair of Theoretical Chemistry and Catalysis Research Center, Technische Universität München, Lichtenbergstrasse 4, D-85747 Garching, Germany S Supporting Information *

ABSTRACT: We present a molecular simulation approach to efficiently determine quantitative rate constants for rare dissolution events from organic crystals. Metadynamics is employed to generate a tailored bias potential that accelerates the escape from the bound state in subsequent hyperdynamics simulations. The robustness and acceleration obtained in the application to kink site dissolution at the aspirin(001)/water interface suggests this technique as a suitable tool to unravel the atomic-scale mechanisms of crystal dissolution.

SECTION: Kinetics and Dynamics

T

trajectory approaches like transition path sampling.11 (ii) The rotational degrees of freedom of the noncentrosymmetric molecular dissolution unit adds additional complexity to the reaction coordinate.12,13 This prevents a low-dimensional description of the reaction pathway and concomitantly regular transition-state theory (TST)-type approaches14 that are based on free energy profiles along corresponding reaction coordinates.9 (iii) The individual transition between bound and dissolved states involves substantially changing entropic contributions,15 excluding simple reaction barrier estimates based on bond-breaking energies that otherwise sufficiently determine, for instance, relative ratios of different rate constants.5,13 Aiming for a method that is both quantitatively reliable and numerically efficient, we here propose to overcome these challenges by combining hyperdynamics16 with metadynamics 17 simulations. We thereby exploit the strong MD acceleration achieved through the use of an optimized bias potential Vbias that, in turn, is custom-tailored to the dissolution process by preceding metadynamics simulations. Complementing the on-the-fly combination of both techniques recently suggested by Tiwary and Parrinello to simulate back-and-forth transitions between bound states18 and over well-defined free energy barriers, we demonstrate, by means of its application to kink site dissolution at the aspirin/water interface, that our two-

he dissolution of crystal surfaces at a given solid/liquid interface plays a central role in processes as diverse as the weathering of materials1,2 or the delivery of active pharmaceutical ingredients (API) from a solid drug compound.3 With the vision of a controlled deceleration or acceleration of such processes, there is a growing demand for tools that provide insight down to the underlying molecular-level mechanisms. Especially with respect to the increasing importance of computer-aided drug design, this includes in silico schemes for the quantitative prediction of macroscopic dissolution rates based on accurate microkinetic models. However, in contrast to the broad attention that, e.g., corresponding theoretical modeling of the “reverse” crystallization from solution has received over the past decades,4,5 only few molecular simulation studies have aimed to go beyond prevalent effective models6,7 by evaluating the detachment of individual atoms or molecules at the solid/liquid interface.2,8,9 As such, there is to date no established framework for the efficient calculation of quantitative molecular rate constants that in turn would serve as the basis for predictive-quality microkinetic models of the dissolution process. In particular, for the dissolution of prototypical hydrophobic API molecules like aspirin, whose dissolution rates play an essential role in their pharmacokinetics,10 this lack is largely due to the concurrence of several characteristics that each alone already challenge existing simulation techniques: (i) Similar to the dissolution of ionic crystals,9 the dissolution process is rare, and simultaneously likely to proceed along rough pathways involving long-lived intermediate states, e.g., due to consecutive breaking of bonds. It is thus neither efficiently accessible by regular molecular dynamics (MD) simulations, nor by reactive © 2014 American Chemical Society

Received: September 12, 2014 Accepted: October 21, 2014 Published: October 21, 2014 3859

dx.doi.org/10.1021/jz501939c | J. Phys. Chem. Lett. 2014, 5, 3859−3862

The Journal of Physical Chemistry Letters

Letter

is known to take place primarily via detachment and attachment of molecules at kink sites.4 As a showcase, we correspondingly focus on the dissolution of a molecule from a selected kink defect at a (100) step shown in Figure 1.

step protocol robustly tackles escape processes via rough and diffusive transition regions at the solid/liquid interface. According to the hyperdynamics formalism, the escape rate constant kHD from an initial state A can be obtained as16 k HD =

k bias k bias = facc ⟨exp(βVbias)⟩A ,bias

(1)

where kbias is the biased rate constant, facc is the so-called acceleration factor, β = 1/kBT, kB is the Boltzmann constant, and T is the temperature. The accuracy of this approach relies crucially on an exhaustive sampling of the exponential average ⟨···⟩A,bias, in fact over the entire biased initial state A. A major challenge is thus to construct a suitable bias potential, which fulfills the requirement of being zero at the barrier region and which facilitates the escape, while not pushing the system out of state A too quickly as to impede a sufficient sampling. Different methods have been proposed,16,19,20 which all essentially aim to modify a smooth potential energy landscape, for instance, governed by atomic bonds or dihedral torsions. Applying such bias to nonbonded potential energies, is, particularly in an explicit solvent environment, prevented by the large fluctuations affecting these energies. In this work we therefore follow a different approach for the construction of the bias potential, consisting of two subsequent steps: 1. An initial metadynamics simulation is run to adaptively build a suitable bias potential acting on previously selected collective variables (CVs). The time-dependent evolution of the bias potential is truncated before the first escape event out of state A, to fulfill the requirement of zero bias at the barrier region.16 2. Using the resulting static bias potential, the hyperdynamics simulation step then generates an ensemble of biased trajectories starting from different conditions in state A, and determines kHD via eq 1, the recorded mean first-passage-times (FPTs), and through the evaluation of the acceleration factor facc. In contrast to the on-the-fly combination of both steps, as suggested in ref 18, the proposed approach does not rely on the same collective variables accelerating both the forward and reverse transition in order to observe a large number of transitions. Instead we are able to focus only on the forward reaction. This aspect becomes particularly beneficial for the present bound-to-free transition and the associated largely changing entropic contributions.15 Moreover, the separation of the two steps allows us to control that the bias potential always vanishes in the barrier region. Therefore, the trajectory segments that follow the first barrier crossing can be exploited to explicitly calculate a transmission coefficient κ to additionally account for recrossing effects.14 Thus, by augmenting the TST approximation toward a more general reactive flux rate constant, k = κkHD, the requirement of having a sharp, welldefined transition barrier (which the original hyperdynamics scheme inherits from TST) can be substantially weakened. In our aspirin/water application, this led to an impressive robustness of k with respect to the exact choice of the transition boundary (cf. Supporting Information). We correspondingly expect the suggested algorithm to be more generally applicable to systems with complex transition regions. We demonstrate and validate our approach with the application to the aspirin (001)/water interface, which is one of the dominant facets of aspirin crystals grown from solution.21 At near-equilibrium conditions, crystal dissolution and growth

Figure 1. Side (a) and top (b) view of the considered kink defect at the aspirin(001)/water interface. The top view features only the upper terrace layer and labels the selected kink site molecule by bold atomic spheres.

Pursuing the simulation scheme suggested above, we first prepare a suitable bias potential using metadynamics. The interactions of aspirin and water molecules are described by the generalized AMBER force field (GAFF)22 in combination with the TIP3P water model.23 All simulations are carried out in the NVT-ensemble and are performed with the GROMACS package,24 using the PLUMED plugin25 to include the bias potentials. Computational details are provided in the Supporting Information (SI), along with an extensive description of the employed supercell surface model containing 2707 water and 334 aspirin molecules. The center-of-mass (COM) vector of the kink site molecule is employed as the CV that triggers the escape of the molecule from its original location. To comply with the translational invariance of the chosen CV, the atoms in the bottom molecular layer of the crystal slab are tethered to their crystal positions by harmonic position restraints. In contrast to, for instance, the situation reported recently for the detachment of an ion from an inorganic surface,9 no unique, reproducible escape pathway is obtained from repeated initial metadynamics runs. Instead, the aspirin molecule is observed to leave its crystal position only temporarily, typically followed by an immediate reattachment in either the original orientation or in a rotated, less stable arrangement. The resulting distribution of molecular orientations, spanning the ensemble of adsorbed kink site conformations, is displayed in the SI. The final detachment trajectory can start from any of these possible orientations, thus offering various different escape pathways, instead of exhibiting one unique molecular dissolution mechanism. The metadynamics trajectory of the COM components, used to generate the bias potential for the subsequent hyperdynamics simulations, is shown in Figure 2. During the simulation, the COM vector starts to exhibit increased fluctuations after the addition of about 600 hills, and eventually, after approximately 1300 hills, the molecule escapes from the kink site location. To ensure a spatially confined potential, that is zero at the barrier region, only the first 600 hills are included to generate the hyperdynamics bias potential, as the subsequent fluctuations might possibly extend into the transition region. The boundary of the thus constructed bias potential is depicted 3860

dx.doi.org/10.1021/jz501939c | J. Phys. Chem. Lett. 2014, 5, 3859−3862

The Journal of Physical Chemistry Letters

Letter

between both curves starts only well beyond the transition region, as the bias potential has not been constructed for an enhanced sampling of the delocalized detached state. Plotting the accessed Vbias values in Figure 3b equally as a function of d confirms that the bias potential is indeed nonzero only within the free energy well immediately around the adsorbed state. In the major part of the transition region and, more importantly, at the barrier, Vbias vanishes, justifying its construction and allowing one to extract unbiased transition segments from the hyperdynamics trajectories for the calculation of the transmission coefficient κ. For the calculation of the molecular dissolution rate constant, we use a value of d = 19 for the boundary between kink site ensemble and free states, noting that the exact choice does not effect the final results to a significant extent (cf. SI). The hyperdynamics simulations then yield a total of 1071 escape trajectories, including detachment processes that followed a reattachment process after recrossing back into the inner region of the kink site ensemble. The histogram of the individual FPTs confirms that the biased system still obeys the correct first-order kinetics (cf. SI). From the decay constant or the mean FPT we therefore obtain a biased kbias = 8.7 × 108 s−1. The time-dependent transmission coefficient likewise exhibits the expected quick initial decay before it plateaus at a constant value of κ = 0.08 after approximately 200 ps. Together with an acceleration factor of facc = 89, this finally leads to a true rate constant and mean FPT of k = 7.7 × 105 s−1 and 1.3 μs, respectively. This is consistent with the estimate obtained from the unbiased PRD simulations, which exhibit one true escape attempt throughout the 128 trajectories with an accumulated total simulation time of 2 μs. If we employ, for comparison, the computed barrier height of 23 kJ/mol within a harmonic TST estimate (using a generic vibrational prefactor kBT/h), we arrive at a rate constant of 6 × 108 s−1, which overestimates the hyperdynamics rate constant kHD = 9.8 × 106, calculated without the transmission coefficient, by almost 2 orders of magnitude, and thus reveals considerable deviations of the dissolution kinetics from TST. In conclusion, we have presented an efficient algorithm to quantitatively determine dissolution rate constants for molecular crystals. An initial metadynamics step generates a bias potential that is specifically adapted to match the free energy minimum around the initial state, thus efficiently boosting subsequent hyperdynamics simulations. The explicit and computationally unintensive inclusion of the transmission coefficient renders the approach robust and applicable to complex, diffusive barrier regions, as well as allows to drop the strict prerequisite of a well-defined free energy barrier. Already for the comparably fast kink site dissolution at the aspirin (001)/water interface, the approach shows a remarkable acceleration factor and thus speed-up in CPU time of almost 2 orders of magnitude. For an optimized bioavailability of an API like aspirin, its crystal dissolution has to be fast against fluid diffusion. Addressing this, computational drug design requires reliable dissolution rate constants as the central parameter within a microkinetic model, e.g., a modified version of the well-known spiral growth model.4,12 We correspondingly plan to employ the calculated rate constants, along with recently obtained accurate defect formation free energies,15 in order to predict macroscopic interfacial dissolution rates. The benefit of the proposed rate constant calculation within such a scheme becomes obvious, considering that in the absence of established free energy molecular simulation protocols, engineering approaches often rely on bond-breaking energies for simple

Figure 2. (a) Trajectory of the COM components (x, y, z) of the kink molecule during the metadynamics simulation that generates the bias potential. The truncation point is marked by the dotted line (see text). (b) Isosurface corresponding to Vbias = 0.1 kJ/mol, depicting the outer boundary of the bias potential distribution.

in Figure 2b, confirming that all nonzero values concentrate exclusively within a small region around the kink site location. To initiate the hyperdynamics simulations based on this bias potential, we employ an ensemble of 128 independent starting conformations. Each hyperdynamics simulation is then run for 15 ns to also capture comparably long biased escape times. Aiming to acquire an unbiased reference, the same set of simulations is additionally run without invoking any bias potential, then in the spirit of parallel-replica dynamics (PRD).26 Both ensembles of trajectories are first analyzed to establish the free energy profile, applying a reweighting scheme to the biased ensemble.19 As a suitable dimensionless order parameter describing the dissolution process, we thereby consider the covariance-weighted distance from the average equilibrium kink site position d 2 = (rCOM − rav)T · C −1·(rCOM − rav)

where C is the covariance matrix of the fluctuations of the molecule COM from its average position, as calculated from an equilibrium simulation. Both profiles shown in Figure 3a agree very well and exhibit a deep free energy well at small distances d that characterizes the perfect kink site location. This is followed by a transition region with a reduced slope that at around d = 19 leads to a free energy barrier, separating the adsorbed ensemble from free states. Both curves give almost exactly the same well-to-barrier height of 23 kJ/mol. The disagreement

Figure 3. (a) Free energy profile associated with the dimensionless distance d, obtained by reweighting the hyperdynamics trajectories (black) and from the unbiased PRD trajectories (gray/red). The dotted line indicates the location of the barrier separating adsorbed and free states (see text). (b) Scatter plot of the encountered Vbias values as a function of d. 3861

dx.doi.org/10.1021/jz501939c | J. Phys. Chem. Lett. 2014, 5, 3859−3862

The Journal of Physical Chemistry Letters

Letter

estimates.5,13 For the present system and using a standard generalized Born implicit solvation model27 the potential energy difference between free and kink site aspirin amounts to 108 kJ/mol. Without any further treatment of entropy, evaluating this thermodynamic barrier at the TST level would already underestimate the true dissolution rate constant by more than 10 orders of magnitude.

■

(14) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. Dynamics of Ion Pair Interconversion in a Polar Solvent. J. Chem. Phys. 1990, 93, 7137−7147. (15) Schneider, J.; Zheng, C.; Reuter, K. Thermodynamics of Defects at the Aspirin/Water Interface. J. Chem. Phys. 2014, 141, 124702− 124713. (16) Voter, A. F. A Method for Accelerating the Molecular Dynamics Simulation of Infrequent Events. J. Chem. Phys. 1997, 106, 4665−4677. (17) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−12566. (18) Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys. Rev. Lett. 2013, 111, 230602−230605. (19) Hamelberg, D.; Mongan, J.; McCammon, J. A. Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120, 11919−11929. (20) Miron, R. A.; Fichthorn, K. A. Accelerated Molecular Dynamics with the Bond-Boost Method. J. Chem. Phys. 2003, 119, 6210−6216. (21) Kim, Y.; Matsumoto, M.; K, M. Specific Surface Energies and Dissolution Behavior of Aspirin Single Crystal. Chem. Pharm. Bull. 1985, 33, 4125−4131. (22) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (23) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (24) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (25) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (26) Voter, A. F. Parallel Replica Method for Dynamics of Infrequent Events. Phys. Rev. B 1998, 57, R13985−R13988. (27) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129.

ASSOCIATED CONTENT

S Supporting Information *

The supporting material contains computational details, as well as additional investigations regarding convergence and sampling. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We gratefully acknowledge funding by the German Research Council, DFG, under Grant No. RE 1509/18-1.

■

REFERENCES

(1) Lüttge, A.; Arvidson, R. S.; Fischer, C. A Stochastic Treatment of Crystal Dissolution Kinetics. Elements 2013, 9, 183−188. (2) Lasaga, A. C.; Lüttge, A. Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model. Science 2001, 291, 2400− 2404. (3) Dokoumetzidis, A.; Macheras, P. A Century of Dissolution Research: From Noyes and Whitney to the Biopharmaceutics Classification System. Int. J. Pharm. 2006, 321, 1−11. (4) Burton, W. K.; Cabrera, N.; Frank, F. C. The Growth of Crystals and the Equilibrium Structure of their Surfaces. Philos. Trans. R. Soc. 1951, 243, 299−358. (5) Lovette, M. A.; Robben Browning, A.; Griffin, D. W.; Sizemore, J. P.; Snyder, R. C.; Doherty, M. F. Crystal Shape Engineering. Ind. Eng. Chem. Res. 2008, 47, 9812−9833. (6) Nernst, W. Theorie der Reaktionsgeschwindigkeit in Heterogenen Systemen. Z. Phys. Chem. 1904, 47, 52. (7) Edwards, L. The Dissolution and Diffusion of Aspirin in Aqueous Media. Trans. Faraday Soc. 1951, 47, 1191−1210. (8) Piana, S.; Gale, J. D. Understanding the Barriers to Crystal Growth: Dynamical Simulation of the Dissolution and Growth of Urea from Aqueous Solution. J. Am. Chem. Soc. 2005, 127, 1975−1982. (9) Stack, A. G.; Raiteri, P.; Gale, J. D. Accurate Rates of the Complex Mechanisms for Growth and Dissolution of Minerals Using a Combination of Rare-Event Theories. J. Am. Chem. Soc. 2012, 134, 11−14. (10) Amidon, G.; Lennernäs, H.; Shah, V.; Crison, J. A Theoretical Basis for a Biopharmaceutic Drug Classification: The Correlation of in Vitro Drug Product Dissolution and in Vivo Bioavailability. Pharm. Res. 1995, 12, 413−420. (11) Dellago, C.; Bolhuis, P. G.; Geissler, P. L. Transition Path Sampling. Adv. Chem. Phys. 2002, 123, 1−78. (12) Liu, X. Y.; Boek, E. S.; Briels, W. J.; Bennema, P. Understanding the Barriers to Crystal Growth: Dynamical Simulation of the Dissolution and Growth of Urea from Aqueous Solution. Nature 1995, 374, 342−345. (13) Kuvadia, Z. B.; Doherty, M. F. Spiral Growth Model for Faceted Crystals of Non-Centrosymmetric Organic Molecules Grown from Solution. Cryst. Growth Des. 2011, 11, 2780−2802. 3862

dx.doi.org/10.1021/jz501939c | J. Phys. Chem. Lett. 2014, 5, 3859−3862