Article pubs.acs.org/molecularpharmaceutics
Polymorphism of Paracetamol: A New Understanding of Molecular Flexibility through Local Methyl Dynamics Nikolaos Tsapatsaris,*,† Boris A. Kolesov,‡,§ Jennifer Fischer,∥ Elena V. Boldyreva,§,⊥ Luke Daemen,# Juergen Eckert,#,○ and Heloisa N. Bordallo*,†,◆ †
European Spallation Source ESS AB, P.O. Box 176, Lund 221 00, Sweden Institute of Inorganic Chemistry SB RAS, Novosibirsk 630090, Russia § REC-008 Novosibirsk State University, ul. Pirogova 2, Novosibirsk 630090, Russia ∥ Forschungszentrum Jülich, ICS-3, Jülich 52425, Germany ⊥ Institute of Solid-State Chemistry and Mechanochemistry SB RAS, Kutateladze 18, Novosibirsk 630128, Russia # Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States ○ Department of Chemistry, University of South Florida, Tampa, Florida 33620, United States ◆ Niels Bohr Institute, Universitetsparken 5, Copenhagen 2100, Denmark ‡
S Supporting Information *
ABSTRACT: This study focuses on the interplay of molecular flexibility and hydrogen bonding manifested in the monoclinic (form I) and orthorhombic (form II) polymorphs of paracetamol. By means of incoherent inelastic neutron scattering and density functional theory calculations, the relaxation processes related to the methyl side-group reorientation were analyzed in detail. Our computational study demonstrates the importance of considering quantum effects to explain how methyl reorientations and subtle conformational changes of the molecule are intertwined. Indeed, by analyzing the quasi elastic signal of the neutron data, we were able to show a unique and complex motional flexibility in form II, reflected by a coupling between the methyl and the phenyl reorientation. This is associated with a higher energy barrier of the methyl rotation and a lower Gibbs free energy when compared to form I. We put forward the idea that correlating solubility and molecular flexibility, through the relation between pKa and methyl rotation activation energy, might bring new insights to understanding and predicting drug bioavailability. KEYWORDS: molecular drugs, polymorphism, inelastic neutron scattering, methyl rotation, hydrogen bonding, DFT calculations
■
INTRODUCTION
Moreover, the enzyme activity is modulated by cofactors, coenzymes, and the presence of other molecules. Considering that crystallographic techniques alone are unable to reveal the main enzyme−ligand recognition and that the understanding of the mechanism of the binding processes at the microscopic level is of key importance,3 recent computational studies combined with neutron spectroscopy and diffraction are becoming more common in elucidating biological function through dynamical and static aspects of a biochemical process.4−7 For example, measuring the dynamics of a ligand allows identification of states where significant molecular reorganization takes place. These are associated with entropic changes and alteration of the local force constants that, in turn, can serve as a guide to theoretical models in molecular dynamics simulations or in drug design.8 Indeed, the importance of
Paracetamol, also known under the names acetaminophen, panadol, and p-hydroxyacetanilide, is one of the world’s most widely used medications.1 In his Nobel Prize-winning work on the mechanism of action of aspirin and other nonsteroidal antiinflammatory drugs (NSAIDs), Vane2 demonstrated that this family of molecular drugs inhibits the formation of prostaglandins (PGs), which are associated with pain, fever, and inflammation. Paracetamol, however, did not appear to inhibit PG synthesis despite having effects similar to those of the NSAIDs. In reality, even after extensive computational, structural, and biochemical investigations, the exact molecular mechanism of action of this amino-phenol is not completely understood.1 Currently, there are a number of biochemical studies relating paracetamol activity to the inhibition of the multidomain cyclooxygenase enzymes (COX1, 2, and 3) through interruption of phenoxyl radical formation. Such biochemical processes are complex because they can lower the energy barrier to the creation of products in a number of possible ways. © 2014 American Chemical Society
Received: Revised: Accepted: Published: 1032
November 25, 2013 January 31, 2014 February 7, 2014 February 7, 2014 dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
crystallographic a-axis direction (i.e., along the direction close to the molecular axes and dominated by OH···O). In contrast, contraction is observed along the direction of the molecular layers (parallel to c axis) and in the softest direction of the NH··· O (b axis).16,20,21 The observed anisotropy, which can also be described as a uniaxial negative thermal expansion (NTE), in the crystallographic a axis in form II has a marked temperature dependence, which is accentuated at temperatures below 160 K. Such thermal expansion has been previously reported for molecular crystals such as NMA22 and L-MET23 and has been interpreted as the result of anharmonicity of the interatomic potential. At the same time, experimental as well as computational results report an almost identical C−H bond distance of the methyl rotor in the monoclinic form I but significant bond length differences in the orthorhombic.16 According to ab initio calculations, an individual paracetamol molecule should be flat24 (i.e., the torsion angle equal to 0°), but in the experimental crystal structures, the torsion angle is equal to 21.1° in form I and to 17.7° in form II. This creates an additional degree of conformational freedom, related to the phenyl group, which can be expressed differently in the two polymorphs according to the local intermolecular environment. Such structural differences profoundly affect their solubility, dissolution rate, pressure stability, structure, and lattice dynamics.25,26 The reason for the relative stability of the metastable form II can be sought after in the primary role hydrogen bonding plays in the two polymorphs. These attractive interactions are caused by the dual proton-donor/acceptor nature of the hydrogen in the −OH group, whereas the −NH group is a donor and the −CO group is an acceptor within the molecule (N−H···O−H···OC → N−···H−O···H−O+C). The energy cost of breaking the many hydrogen bonds during a phase change from form II to form I is significant and suggests that the solid-state transformations between the polymorphs during temperature or pressure variations are kinetically hindered. Clearly, the conformation of paracetamol molecules is very sensitive to the local electronic environment and is balanced via van der Waals interactions following a tendency of minimum energy optimization of the NH···O and OH···O intermolecular bonds. Thermodynamic parameters, such as the proton dissociation constant, Ka, are also affected by differences in hydrogen bonding between different structures. The molecular units in forms I and II are arranged in an antiparallel and parallel orientations, respectively, resulting in a different local environment for the terminal methyl group. Therefore, it is expected that their energy barrier to reorientation will be different and may correlate to the free Gibbs energy difference in the two forms. For example, when reducing the thermal energy of the crystal, effectively reducing the inter- and intramolecular space, nonbonding interactions between C(H)3···(H4)C6 hydrogen atoms of the phenyl ring and methyl proton density will affect the local potential. This was corroborated by our earlier study, which showed that a decrease of the energy of the methyl rotation barrier height, V, was observed25 under pressure. Furthermore, methyl groups undergo large amplitude motions, especially at high temperatures,27 creating significant elongations of the apparent methyl group C− H bond distances.28 Interestingly, although this molecular drug has been the focus of considerable research,15−18,20,29−34 information on the side-group dynamics and their relation to the stability of the polymorphic structures has not yet been fully explored. In particular, the low-frequency vibrations of paracetamol polymorphs contain essential information on the coupling between the internal motions of molecular fragments
molecular group dynamics was highlighted in a computational study of galectin9 in which ligand binding accompanied significant loss of translational and rotational freedom of the adjacent hydration layer, facilitating active site−ligand recognition. Moreover, understanding the dynamics of individual molecular fragments to model the dynamics of larger molecules has attracted considerable attention. For example, spectroscopic investigations using neutron scattering as well as quantumchemistry calculations of secondary amide derivatives have been used to understand the force field of proteins and poly amino acids.10 Indeed, compounds such as acetamide,11,12 formamide,13 thioamide,14 and their derivatives served as simple models for investigating the nature of secondary amide functions. In the present study, we focus on the dynamical aspects of the two stable paracetamol polymorphs by the combined use of neutron spectroscopy and ab initio tools. This facilitates the assignment of methyl rotational modes, which are crucial for the elucidation of conformational flexibility in paracetamol. Paracetamol is composed of a benzene ring core, substituted at one terminus by one −OH group and on the other by an amide bonded with a −CH3 group, and can crystallize in three polymorphic states. The monoclinic polymorph, hereafter named form I, crystallizes in the space group P21/a (a = 12.930(40) Å, b = 9.400(10) Å, and c = 7.100(20) Å),15 whereas the orthorhombic polymorph, hereafter named form II, crystallizes in the space group Pcab (a = 7.406(1) Å, b = 11.837(7) Å, and c = 17.162(20) Å),16 and the highly unstable form III, in space group Pca2117 (Figure 1). Considering the
Figure 1. Molecular structure of paracetamol forms I and II after performing periodic zero-point structural refinement using VASP. The marked intramolecular C−H distances are calculated computationally and show that in the relaxed structures there are significant conformational differences and a significant rotor asymmetry in the orthorhombic form II (distances are given in angströms; blue, nitrogen and red, oxygen).
difference in the free Gibbs energy between forms I and II, form II is metastable at all temperatures at ambient pressure and was shown to transform to form I upon heating or cooling and ambient pressure over time and in the presence of humidity.18,19 In addition, the structures of forms I and II are different: form I is composed of corrugated layers containing four asymmetric units, whereas form II is composed of parallel planes containing eight asymmetric units. There is marked anisotropy of the structural strain in forms I and II, both upon cooling and with increasing pressure. Upon cooling, the interlayer compression in form II is larger than in form I, and the structure expands slightly in the 1033
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
Figure 2. (a) Inelastic neutron scattering intensity as a function of energy measured at 2 K at FDS and (b) compared to density functional theory calculations using VASP in the Brillouin center and 0 K for both forms. Gray dotted line indicates additional dispersion-analyzed vibrations for the monoclinic form using a 4 × 4 × 4 k-mesh. The computational resources were insufficient for dispersion analysis of the orthorhombic form II. (c, d) Experimental incoherent inelastic scattering collected from the NEAT instrument and converted to χ″(ω) as a function of temperature for forms I (c) and II (d). Insets: Temperature evolution of the intensity of the barycenter of the strongest spectral bands at ∼10 and ∼14 meV for forms I and II, respectively (denoted by asterisk). The intensities of both forms have been fitted using an exponential dependence of the form ∼1 − eW/kBT. Note that 1 meV ≈ 8 cm−1.
range −1.7 to 500 meV, covering 40 ps to 8 fs in the time domain. Details on the instrumental setup are given in the Supporting Information. Because of the simplicity of the neutron−nucleus interaction, the observed spectra can be easily compared with computed results, as the intensities need to be scaled only by the neutron scattering cross sections of the nuclei. Computational Details. The published experimental crystal structures15,16 were taken as the reference for the atomic coordinates used in our ab initio calculations. Calculations on gas-phase paracetamol polymorphs at 0 K were carried out at the B3LYP/6-31G** level of theory with the GAMESS-UK package.36 Periodic calculations on crystalline paracetamol forms I and II were performed with the VASP37 package using the PBE functional along with Vanderbilt ultrasoft pseudopotentials38 with a plane wave kinetic energy cutoff of 450 meV. Additional periodic calculations were carried with a 4 × 4 × 4 Monkhorst−Pack mesh of k-points on the monoclinic form to improve agreement at lower frequencies. The rotational profiles around the C−CH3 group were determined by varying the corresponding dihedral angle in steps of 10° from 0 to 120° for the monoclinic form and 0 to 360° for the orthorhombic form, respectively. The rotational energy barrier for the methyl group was then derived from single-point energy calculations in which the methyl group orientation about the C−CH 3 axis was incremented in 10° steps in a frozen atom approach.
and the lattice modes. Here, by combining neutron spectroscopy with ab initio tools, we investigated the behavior of weak intraand intermolecular interactions on different time scales as a function of temperature to assess this missing information (i.e., how the anisotropy of the hydrogen bonds affects the conformational freedom of paracetamol forms I and II). The combination of density functional theory (DFT) and incoherent inelastic neutron scattering (IINS) allowed for the assignment of the librational modes, which are crucial in the identification of anomalous behavior because they are related to the degree of rigidity of the amide chain. This, in turn, can be associated with the anisotropy of the structure and consequently to its degree of conformational flexibility.
■
MATERIALS AND METHODS Preparation of Paracetamol Polymorphs. Form I was obtained after recrystallization of commercial monoclinic paracetamol from an ethanol solution and subsequent drying for 1 day in a vacuum oven. Paracetamol form II was prepared by recrystallization of a melt. The phase purity of the purchased and synthesized paracetamol polymorphs was assessed using a Bruker D8 Advance diffractometer with Cu Kα radiation (1.5418 Å) in θ/θ geometry. The diffraction patterns, Figure S1 in the Supporting Information, are in full agreement with previous results corresponding to from I (monoclinic, P21n)35 and form II (orthorhombic, Pcab).16 The samples were kept dry during the experiments. Incoherent Inelastic Neutron Scattering Experiments (IINS). The dynamic behavior of polycrystalline paracetamol forms I and II was studied by means of incoherent inelastic neutron scattering (IINS) using spectrometers NEAT and FDS at the HZB (Berlin, Germany) and LANSCE (Los Alamos, USA), respectively. The combination of these two spectrometers allowed for the analysis of motions observed over the energy
■
RESULTS AND DISCUSSION Lattice and Internal Vibrations Assignment. To identify the main vibrational modes for both polymorphs, we started our analysis with the IINS data measured using FDS together with the temperature-dependent susceptibility data from NEAT. Neutrons scattered incoherently by hydrogen atoms dominate
1034
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
Figure 3. (a,b) Elastically scattered intensity integrated over the elastic line obtained from the NEAT spectrometer as a function of temperature and momentum transfer for forms I and II. (c, d) Temperature evolution of the quasielastic broadening measured from NEAT. (e) Mean-square displacements, ⟨u2⟩, derived from the elastic-scattered intensity and modeled according to a harmonic Debye model. (f) Quasielastically scattered intensity, obtained at 240 K, fitted using the phenomenological model of eq 4 at an average momentum transfer of 1.35 Å−1 for form I (blue) and form II (red).
energy band between 18 and 25 meV, we have identified CH3 methyl librations that appear to be mixed predominantly with −C6H4 librations/rotations around its −N−C−C− center symmetry axis (solid red arrow). The vibrational band at approximately 40 meV (Figure 2a) has been correctly reproduced for both forms by our DFT calculations (Figure 2b). These can be attributed to mixed modes involving most atoms in the lattice. However, the vibrational density of this band also includes significant bending and translational (C−CH3) modes in the form of out-of-plane bending and translations as well as significant −OH bending modes, which are directly related to the strength of hydrogen bonding across the crystallographic axis a. In the inset of Figure 2a,b, the modes located at approximately 95 and 100 meV for the monoclinic form and at 105 and 120 meV for the orthorhombic form can be assigned to the out-ofplane bending of the N−H amide group.40 These bands also contain O−H and N−H bending modes, highlighting differences in the strength and symmetry of the hydrogen bonding of the donor −NH and −OH groups. The most striking finding in this spectral region is that the −OH bend modes occur at a much higher energy in form II (+8 meV). The result indicates that form II is less flexible than form I along axis a and b, which are dominated by hydrogen bonds formed by O···HO. In the spectral region of 125−410 meV, not shown in Figure 2, characteristic of internal vibrations, we find agreement between
the total scattered intensity; thus, one can assume that only the self-dynamics of these atoms are being probed. The spectral assignment of the main vibrational bands was achieved by visualizing the molecular vibrational frequencies created by VASP using Jmol39 and by comparison with literature studies. Bands in the low-energy part of the vibrational data in molecular crystals correspond largely to intermolecular collective vibrations and torsional modes ( TD, ⟨u2⟩ increases linearly with temperature, and conversely, at low temperatures, the meansquare fluctuations become temperature-independent, approaching a plateau.44 Assuming such a model and within the statistical accuracy of our data, on the measured time scale, form II can be described as an harmonic system up to 140 K, whereas for form I, anharmonic contributions are evident above 60 K. For monoclinic form I, a deviation from the Debye model happens at a lower temperature than form II. This suggests that the proton motions in form I exhibit higher disorder and are less hindered than in form II. This is in full agreement with the anisotropic thermal parameters deduced from previous crystallographic study.45 To quantify the energetics of the different thermal-activated dynamical processes, the quasielastic part of the spectra was analyzed. The accessible dynamical phenomenon in the time window corresponding to the resolution of the instrument is approximately 1.3−40 ps (ΔE ∼ 100 μeV). The simplest model to analyze such data is to use a single Lorentzian function that corresponds to the population statistics of one relaxation process. This approach is represented by the convolution of a delta function with the instrumental resolution and a Lorentzian function (eq 4).
Γ(T) = Γ0 e−Eact / kBT
Γ0 is the rotational constant and describes the attempt frequency of the stochastic process at infinite temperature. In the case of paracetamol, we observed that the activation energy of the methyl group reorientation corresponds to Eact form II = 35 meV (4.2 ± 0.08 kJ/mol) and Eact form I = 29.5 meV (2.8 ± 0.16 kJ/mol) (rotational constant Γo = 0.5 kJ/mol). The methyl group rotation seems to be more hindered in form II, suggesting stronger intermolecular forces than form I. Naively assuming that the dynamics of the methyl reorientation are contained in a harmonic potential with three-fold symmetry, it is possible to estimate the V3 potential barrier and the first librational substate from the activation energy of the methyl reorientation using the simplified numerical relations E0−1 (meV) = 0.47[V3(K)]0.548 and Eact(K) = 0.598[V3(K)]1.05, which have been described previously.47 The obtained values are given in Figure 4c and d. Our single-point energy calculations, fitted by a V3 potential, agree well with the QENS-derived V3 barriers for form I (V3 = 30 meV QENS vs V3 = 33 meV SP DFT) (Figure 4c,e) but disagree markedly for form II (V3 = 43 meV QENS vs V3 = 62 meV SP DFT) (Figure 4c). Considering the DFT-derived data, the difference in the height of the potential barrier between forms II and I is ΔV = V3 formII − V3 formI = 62 − 33 = 29 meV or ∼2.8 kJ/mol. With the aim of relating the findings from the IINS experiments and DFT calculations to well-known thermodynamic parameters, we now turn to comparing the potential barrier height of the methyl rotation with the free Gibbs energy, G, in forms I and II. The difference of G in the two polymorphs can be calculated using ΔG = ΔH − TΔS using thermodynamic parameters derived from calorimetry,18 as ΔG = GformII − GformI = 32.2 − 18.8 = 13.4 kJ/mol at 300 K. This result indicates that a positive difference in the free Gibbs energy is accompanied by a positive difference in the height of the potential barrier between forms II and I. This demonstrates that the higher thermodynamic stability of monoclinic form I compared to orthorhombic form II is juxtaposed by a more unhindered, higher-entropy, low-barrier methyl rotation. Conversely, the lower thermodynamic stability of orthorhombic form
S(Q , ω) = DWF[(A 0(Q )δ(ω) + (1 − A 0(Q )L(Γ, ω)) ⊗ R(Q , ω)] + B(Q )
(5)
(4)
DWF is the Debye−Waller factor, A0(Q) or EISF is the timeaverage spatial distribution of all scatterers, L(Γ,ω) is the Lorentzian representing the broadened energy distribution because of neutron−nucleus collision, R(Q,ω) is the instrumental resolution function, B(Q) is a background term, and δ(ω) is a delta function. Clearly, the quasielastic contribution observed in both polymorphs, Figure 3f, can be well-described using this simple 1037
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
Figure 5. Elastic incoherent structure factors extracted after the quasielastic fitting of the data and subsequently fitted using Bessel functions corresponding to different symmetry operations. (a) Data fitted using a three-fold rotational symmetry for form I. Inset: Radius of rotation, R, of the protons about their symmetry axis. (b) Data fitted using two-jump + three-jump model. Inset: a pure three-jump model.
A′ 0(Q ) (EISF) = p + (1 − p)[1 + 2j0 ( 3 QR CH3)]/3
II is reflected by a higher potential barrier, a more hindered, lower-entropy methyl rotation. This can be related to the tendency of form II to transform over time to form I at all temperatures and constant volume. Furthermore, it is of interest to relate our findings with the well-known proton dissociation constant, Ka, widely used to indicate solubility and bioavailability of pharmaceuticals. The magnitude of K a is directly related to the underlying thermodynamics of the dissociation reaction in solution and is proportional to ΔGsolution. In Figure 4c, a comparison of the potential barrier of the methyl group with proton dissociation values (in the form of the dissociation logarithm constant pKa = −log10 Ka) for different molecules is shown. We observed that paracetamol form I, a phenol-based acid, has a higher pKa = 9.948 than aspirin (pKa = 3.5),49 a carboxylic-based acid. Although the high pKa of form I is associated with low aqueous solubility and low bioavailability, its lower potential energy barrier in the solid state is associated with a higher thermodynamic stability, implying a lower Glattice. A crystal can be dissolved if solvation is accompanied by a negative change of the Gibbs free energy or if the Gsolution < Gsolvent + Glattice. Moreover, if we consider that the height of the potential barrier is a measure of entropy in the system and that Ssolution > Ssolvent + Slattice (not considering enthalpy), then the high entropy contribution of the methyl reorientation in form I suggests that the interaction of the solvent with the lattice is not favorable and the dissociation is hindered. In contrast, the asymmetrically deformed methyl group in aspirin50 has a higher potential barrier and its rotation is severely hindered with a lower contribution to the Slattice. As a result, the dissolving process is preferred as a means to increase the entropy of the combined (Ssolvent + Slattice) system. Nevertheless, albeit large, the entropic contributions are only part of the system’s thermodynamical picture, and to generalize this idea, a significant number of pKa values and QENS from other methyl-containing molecular crystals combined with considerations of noncovalent interactions, enthalpic contributions, and computational molecular calculations would be needed for a quantitative comparison. Geometry of the Methyl Reorientation and Mixed Motions. To clarify the geometry of the methyl rotation and the apparent discrepancy between the single-energy point DFT calculations and our QENS-estimated three-fold potential V3 barrier of form II, we have modeled the elastic incoherent structure factor (EISF). As seen in Figure 5, the elastic intensity does not disappear at high momentum transfer. This means that there is a significant immobile population on the time scale of the experiment. To account for these populations, the EISF is written as
(6)
where Á 0(Q) is the apparent EISF and p and (1−p) denote the immobile and mobile hydrogen-population fractions, respectively. Details of the analysis are given in the Supporting Information. The reorientation of the methyl rotor in monoclinic form I seems to follow well the three-fold symmetry as shown in Figure 5a. Above 140 K, the EISF signal changes only slightly as the motions become too fast to be observed in the energy window of the measurement. Moreover, the symmetry of the methyl reorientation does not seem to be affected as a function of temperature. In addition, the radius of rotation, RCH3, is approximately 1 ± 0.1 Å, which corresponds well with the reported31 C−H bond distances of 1.09 Å. From eq 6, the jump distance, 31/2RCH3 = 1.73 ± 0.09 Å, and the immobile fraction, pCH3 = 0.79 ± 0.04, were calculated. This analysis confirms that for form I, the observed motions are primarily due to the methyl protons rotating about a C3v symmetry axis. The analysis of the EISF obtained for the orthorhombic form, however, required an additional component in the form of a three-jump combined with a two-jump model to allow a reasonable fit (Figure 5b and inset). Our EISF was modeled with the assumption that these two motions happen on the same time scale. This is not an unreasonable assumption because, as we have evidenced from IINS and DFT, the methyl and phenyl motions are coupled. ⎛ A′ 0(Q ) (EISF) = p + (1 − p)⎜⎜(1 − pCH ) 3 ⎝ [1 + j0 (2QR C6H4)] 2 [1 + 2j0 ( 3 QR CH3)] ⎞ ⎟⎟ + pCH 3 3 ⎠
(7)
The p fraction represents the sum of the immobile populations of both two-jump and three-jump reorientation, whereas the fractions pCH3 and (1−pCH3) represent the portion of the elastic intensity corresponding to the three-jump methyl reorientation and two-jump phenyl reorientation. RC6H4 is the distance of the phenyl protons to a two-fold symmetry axis, C2v. From this analysis, it is clear that ∼90% of the observed QE component is due to methyl protons, whereas the remaining ∼10% originates from phenyl protons. Both methyl and phenyl proton motions 1038
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
lower, mainly methyl, proton disorder than form I (Figure 3e). The interplay between disorder and compressibility can be sourced in the anisotropy of the intermolecular environment of the methyl group. Form II is mostly compressible along crystallographic axis c (where no hydrogen bonding is present). Therefore, in form II, the methyl rotations are severely hindered because of strong intermolecular forces along both axes a and b. Conversely, although form I is on average slightly less compressible (flexible) than form II, both axes b and c are largely and equally compressible and the methyl group rotation is significantly unhindered. Moreover, in form II, the substantial compressibility of the c axis, where weak electrostatic interactions are dominant, provides additional evidence for phenyl reorientational jumps. Paracetamol Activity: Conformational Flexibility and Induced Fitting. Paracetamol boasts significant conformational flexibility, which is demonstrated through its ability to (a) exist in completely different crystallographic environments, (b) perform low-energy relatively methyl reorientations (form I 2.8 kJ/mol and form II 4.2 kJ/mol), and (c) exhibit phenyl−methyl mixed motions. This leads us to evoke the model of induced fitting58 and conformational proofreading59,60 whereby the enzymatic active site is rearranged and a structural mismatch is introduced in the ligand so that energetically favorable binding can take place with high specificity. Paracetamol (and other similar phenols) has a degree of conformational flexibility that allows intermediate conformational transitions that may consequently enhance induced fitting. This is similar to the changes that take place in molecular fragments during subtle polymorphic transformations in the solid state. Lastly, we would like to put forward the idea that because the weakly dissociating monoclinic paracetamol with pKa = 9.9 is accompanied by an extremely low V3 potential barrier (Figure 5c) establishing a correlation between solubility and potential barrier may be a promising avenue in improving predictions of bioavailability.
will contribute to the Arrhenius-estimated activation energy, and as a result, our calculated Eact represents a mixture of these motions. The fitting values are found in the Supporting Information Table 1. The apparent asymmetry and mixed motions of the methyl group in orthorhombic form II, which is corroborated by our selfconsistent structural optimizations and X-ray studies,45 has several important implications. Quantum Mechanical Tunnelling. Tunnelling in molecular crystals containing methyl groups is observed because of the identical nature (degeneracy) of the wave functions of the hydrogen protons in −CH3. The local potential of −CH3 protons is affected by the intermolecular structure, or lattice modes, thus modifying the proton wave functions. As the proton wave functions become different, nondegenerate states are formed. This results in a significantly lower probability of proton tunnelling and a substantially reduced attempt tunnel frequency.31 Our computational calculations indicate that for form I a potential barrier of V3 = 33 ± 0.5 meV results in a tunnel splitting of 8 μeV. This apparent disagreement with the measured value of 3 μeV42 can be related to quantum effects.51 However, for form II, the calculated barrier (V3 = 62 ± 0.5 meV) results in a tunnel splitting of 0.4 μeV, which lies below the instrumental resolution of our previous measurements.25 Prediction of the First Librational Substate E01. Excellent agreement was found between the QENS-predicted V3 barrier and E01 and our ab initio single-point (SP) energy calculations and DFT for monoclinic form I (E(01) EXPERIM fi = 12.1 meV and E(01) COMPUT fi = 13 meV). In contrast, significant differences between experimental QENS (E(01) EXPERIM fii = 14 meV) and computational results (E(01) COMPUT = 18 meV) were observed in the orthorhombic form. To rationalize these differences, we must consider that (a) the coupled methyl−phenyl motions affect the methyl rotational symmetry and consequently it cannot be described by our simple numerical models52 and (b) methyl reorientational dynamics are partly of quantum mechanical nature53,54 even at room temperature. The partial loss of quantum coherence, which is expected to occur as the temperature is increased, makes the rotation increasingly incoherent,55 and these processes can be governed by different rate constants. As a result, the potential barriers, calculated through the activation energy using QENS, will always be lower than the actual barriers because of quantum effects. It is important to consider this effect, especially in coupled-motion systems. Local Frustration, Flexibility, and Proton Disorder. The significant asymmetry of the potential surface of the methyl rotor leads us to discuss concepts of frustration and anisotropy and to relate this to the relatively common effect in molecular crystals of negative thermal expansion (NTE). Upon cooling to 100 K, a small expansion of the crystallographic axis a (Δd/d = 0.2%, approximately in the direction of the −N−C−C−H3− chain) in orthorhombic form II has been observed in a temperature dependent XRD study.45 The low compressibility in this crystallographic direction can be explained by electrostatic repulsions of the methyl protons and the phenyl C−H groups. Most importantly, the slight expansion along axis a can be elucidated by considering that the terminal O···HO hydrogen bonds can bend out of plane, effectively increasing the intermolecular space.23,56,57 It is interesting to consider that, although form II is overall slightly more compressible (flexible) than form I, it demonstrates
■
CONCLUSIONS By combining IINS with ab initio and DFT computational tools, this work brings new insight to how the flexibility of the CH3 group and the mixing of motions are manifested in the molecular structure of paracetamol. Furthermore, the observed methyl− phenyl coupling illustrates the importance of considering quantum effects in the analysis of the potential barrier of the methyl group rotation. In addition, we examined the idea of whether the high dissociation constant logarithm pKa can be connected to flexibility and bioavailability or even toxicity. More results from similar molecules need to be gathered to create a coherent global picture of molecular flexibility. Lastly, we conclude that studying molecular drugs in the solid-state form provides for a controlled environment in which the behavior of molecular fragments can be explored without the effects of solutes or the noisy intra- and extracellular environment.
■
ASSOCIATED CONTENT
■
AUTHOR INFORMATION
S Supporting Information *
Detailed description of INS data analysis and fitting results for the different EISF models. This material is available free of charge via the Internet at http://pubs.acs.org. Corresponding Authors
*E-mail:
[email protected] (N.T.). 1039
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
*E-mail:
[email protected] (H.N.B.).
(14) Velkov, Z.; Velkov, Y.; Balabanova, E.; Tadjer, A. First principle study of the structure of conjugated amides and thioamides. Int. J. Quantum Chem. 2007, 107, 1765. (15) Haisa, M.; Kashino, S.; Kawai, R.; Maeda, H. Monoclinic form of para-hydroxyacetanilide. Acta Crystallogr., Sect. B: Struct. Sci. 1976, 32, 1283. (16) Drebushchak, T. N.; Boldyreva, E. V. Variable temperature (100− 360 K) single-crystal X-ray diffraction study of the orthorhombic polymorph of paracetamol (p-hydroxyacetanilide). Z. Kristallogr. 2004, 219, 506. (17) Perrin, M. A.; Neumann, M. A.; Elmaleh, H.; Zaske, L. Crystal structure determination of the elusive paracetamol Form III. Chem. Commun. 2009, 3181. (18) Boldyreva, E. V.; Drebushchak, V. A.; Paukov, I. E.; Kovalevskaya, Y. A.; Drebushchak, T. N. DSC and adiabatic calorimetry study of the polymorphs of paracetamol. An old problem revisited. J. Therm. Anal. Calorim. 2004, 77, 607. (19) Espeau, P.; Ceolin, R.; Tamarit, J. L.; Perrin, M. A.; Gauchi, J. P.; Leveiller, F. Polymorphism of paracetamol: Relative stabilities of the monoclinic and orthorhombic phases inferred from topological pressure-temperature and temperature-volume phase diagrams. J. Pharm. Sci. 2005, 94, 524. (20) Boldyreva, E. V.; Shakhtshneider, T. P.; Vasilchenko, M. A.; Ahsbahs, H.; Uchtmann, H. Anisotropic crystal structure distortion of the monoclinic polymorph of acetaminophen at high hydrostatic pressures. Acta Crystallogr., Sect. B: Struct. Sci. 2000, 56, 299. (21) Boldyreva, E. V.; Shakhtshneider, T. P.; Ahsbahs, H.; Sowa, H.; Uchtmann, H. Effect of high pressure on the polymorphs of paracetamol. J. Therm. Anal. Calorim. 2002, 68, 437. (22) Bordallo, H. N.; Argyriou, D. N.; Barthès, M.; Kalceff, W.; Rols, S.; Herwig, K. W.; Fehr, C.; Juranyi, F.; Seydel, T. Hydrogen in Nmethylacetamide: Positions and dynamics of the hydrogen atoms using neutron scattering. J. Phys. Chem. B 2007, 111, 7725. (23) Fischer, J.; Lima, J. A.; Freire, P. T. C.; Melo, F. E. A.; Havenith, R. W. A.; Filho, J. M.; Broer, R.; Eckert, J.; Bordallo, H. N. Molecular flexibility and structural instabilities in crystalline l-methionine. Biophys. Chem. 2013, 180−181, 76. (24) Binev, I. G.; Vassileva-Boyadjieva, P.; Binev, Y. I. Experimental and ab initio MO studies on the IR spectra and structure of 4hydroxyacetanilide (paracetamol), its oxyanion and dianion. J. Mol. Struct. 1998, 447, 235. (25) Tsapatsaris, N.; Landsgesell, S.; Koza, M. M.; Frick, B.; Boldyreva, E. V; Bordallo, H. N. Polymorphic drugs examined with neutron spectroscopy: Is making more stable forms really that simple? Chem. Phys. 2013, 427, 124. (26) Paukov, I. E.; Kovalevskaya, Y. A.; Arzamastcev, A. E.; Pankrushina, N. A.; Boldyreva, E. V. Heat capacity of methacetin in a temperature range of 6 to 300 K. J. Therm. Anal. Calorim. 2011, 108, 243. (27) Wilson, C. C. Variable temperature study of the crystal structure of paracetamol (p-hydroxyacetanilide), by single crystal neutron diffraction. Z. Kristallogr. 2000, 215, 693. (28) Wilson, C. C. Zero point motion of the librating methyl group in p-hydroxyacetanilide. Chem. Phys. Lett. 1997, 280, 531. (29) Haisa, M.; Kashino, S.; Maeda, H. Orthorhombic form of parahydroxyacetanilide. Acta Crystallogr., Sect. B: Struct. Sci. 1974, 30, 2510. (30) Nichols, G.; Frampton, C. S. Physicochemical characterization of the orthorhombic polymorph of paracetamol crystallized from solution. J. Pharm. Sci. 1998, 87, 684. (31) Wilson, C. C. Neutron diffraction of p-hydroxyacetanilide (paracetamol): Libration or disorder of the methyl group at 100 K. J. Mol. Struct. 1997, 405, 207. (32) Naumov, D. Y.; Vasilchenko, M. A.; Howard, J. A. K. The monoclinic form of acetaminophen at 150K. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1998, 54, 653. (33) Wilson, C. C. Variable temperature study of the crystal structure of paracetamol (p-hydroxyacetanilide), by single crystal neutron diffraction. Z. Kristallogr. 2000, 215, 693. (34) Kolesov, B. A.; Mikhailenko, M. A.; Boldyreva, E. V. Dynamics of the intermolecular hydrogen bonds in the polymorphs of paracetamol in
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS
■
REFERENCES
H.N.B. and N.T. acknowledge the support of the HelmholtzCenter Berlin (Helmholtz Zentrum Berlin, HZB) for providing some of the neutron research facilities used in this work. This work has also benefited from the use of the Manuel Lujan, Jr. Neutron Scattering Center at Los Alamos National Laboratory and funding from the Department of Energy’s Office of Basic Energy Sciences. Los Alamos National Laboratory (LANL) is operated by Los Alamos National Security LLC under DOE contract DE-AC52-06NA25396. E.V.B. acknowledges the Russian Academy of Sciences and the financial support provided by the Russian Ministry of Science and Education (project 14.B37.21.1093). J.E. thanks the Physics and Chemistry of Materials Group (T-1) at LANL for making computing resources available. Last, N.T. thanks Andrew Jackson, Paul Henry, Esko Oksanen, and Hanna Wacklin at the European Spallation Source for many fruitful discussions.
(1) Graham, G. G.; Davies, M. J.; Day, R. O.; Mohamudally, A.; Scott, K. F. The modern pharmacology of paracetamol: therapeutic actions, mechanism of action, metabolism, toxicity and recent pharmacological findings. Inflammopharmacology 2013, 21, 201. (2) Vane, J. R. Inhibition of prostaglandin synthesis as a mechanism of action for aspirin-like drugs. Nature-New Biol. 1971, 231, 232. (3) Bronowska, A. Thermodynamics of ligand-protein interactions: Implications for molecular design. Thermodynamics−Interaction Studies−Solids, Liquids and Gases; Moreno-Pirajan, J. C., Ed.; InTech: Rijeka, Croatia, 2011. (4) Fitter, J. A measure of conformational entropy change during thermal protein unfolding using neutron spectroscopy. Biophys. J. 2003, 84, 3924. (5) Setny, P.; Baron, R.; Kekenes-Huskey, P. M.; McCammon, J. A.; Dzubiella, J. Solvent fluctuations in hydrophobic cavity−ligand binding kinetics. Proc. Natl. Acad. Sci. U.S.A. 2013, 1197. (6) Sakai, V. G.; Arbe, A. Quasielastic neutron scattering in soft matter. Curr. Opin. Colloid Interface Sci. 2009, 14, 381. (7) Chu, X.; Mamontov, E.; O’Neill, H.; Zhang, Q. Temperature dependence of logarithmic-like relaxational dynamics of hydrated tRNA. J. Phys. Chem. Lett. 2013, 4, 936. (8) Ladbury, J. E. Just add water! The effect of water on the specificity of protein-ligand binding sites and its potential application to drug design. Chem. Biol. 1996, 3, 973. (9) Saraboji, K.; Hakansson, M.; Genheden, S.; Diehl, C.; Qvist, J.; Weininger, U.; Nilsson, U. J.; Leffler, H.; Ryde, U.; Akke, M.; Logan, D. T. The carbohydrate-binding site in galectin-3 is preorganized to recognize a sugarlike framework of oxygens: Ultra-high-resolution structures and water dynamics. Biochemistry 2012, 51, 296. (10) Krishnan, M.; Smith, J. C. Reconstruction of protein side-chain conformational free energy surfaces from NMR-derived methyl axis order parameters. J. Phys. Chem. B 2012, 116, 4124. (11) Hase, Y. Theoretical study of the force field and vibrational assignments of acetamide and deuterated analogues. Spectrochim. Acta, Part A 1995, 51, 2561. (12) Uno, T.; Machida, K.; Saito, Y. Out-of-plane vibrations of acetamide and partially N-deuterated acetamide. Spectrochim. Acta, Part A 1971, 27, 833. (13) Tam, C. N.; Bour, P.; Eckert, J.; Trouw, F. R. Inelastic neutron scattering study of hydrogen-bonded solid formamide at 15 K. J. Phys. Chem. A 1997, 101, 5877. 1040
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
Molecular Pharmaceutics
Article
relation to crystal packing and conformational transitions: a variabletemperature polarized Raman spectroscopy study. Phys. Chem. Chem. Phys. 2011, 13, 14243. (35) Haisa, M.; Kashino, S. The monoclinic form of p-hydroxyacetanilide. Acta Crystallogr., Sect. B: Struct. Sci. 1976, 58. (36) Guest, M. F.; Bush, I. J.; Van Dam, H. J. J.; Sherwood, P.; Thomas, J. M. H.; Van Lenthe, J. H.; Havenith, R. W. A.; Kendrick, J. The GAMESS-UK electronic structure package: Algorithms, developments and applications. Mol. Phys. 2005, 103, 719. (37) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169. (38) Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B 1990, 41, 7892. (39) Jmol: An Open-Source Java Viewer for Chemical Structures in 3D; http://www.jmol.org. (40) Barthes, M.; Bordallo, H. N.; Eckert, J.; Maurus, O.; Nunzio, G. De. Dynamics of crystalline N-methylacetamide: Temperature dependence of infrared and inelastic neutron scattering spectra. J. Phys. Chem. B 1998, 5647, 6177. (41) Schober, H. Neutron spectroscopy: An ideal tool for the materials scientist. J. Phys. 2003, 103, 173. (42) Johnson, M. R.; Prager, M.; Grimm, H.; Neumann, M. A.; Kearley, G. J.; Wilson, C. C. Methyl group dynamics in paracetamol and acetanilide: Probing the static properties of intermolecular hydrogen bonds formed by peptide groups. Chem. Phys. 1999, 244, 49. (43) Bordallo, H. N.; Zakharov, B. A.; Boldyreva, E. V.; Johnson, M. R.; Koza, M. M.; Seydel, T.; Fischer, J. Application of incoherent inelastic neutron scattering in pharmaceutical analysis: relaxation dynamics in phenacetin. Mol. Pharm. 2012, 9, 2434. (44) Bordallo, H. N.; Boldyreva, E. V.; Buchsteiner, A.; Koza, M. M.; Landsgesell, S. Structure property relationships in the crystals of the smallest amino acid: An incoherent inelastic neutron scattering study of the glycine polymorphs. J. Phys. Chem. B 2008, 112, 8748. (45) Drebushchak, T. N.; Boldyreva, E. V. Variable temperature (100− 360 K) single-crystal X-ray diffraction study of the orthorhombic polymorph of paracetamol (p-hydroxyacetanilide). Z. Kristallogr. 2004, 219, 506. (46) Laidler, K. J. Chemical Kinetics, 3rd ed.; Harper & Row: New York, 1987. (47) Smuda, C.; Busch, S.; Wagner, B.; Unruh, T. Methyl group dynamics in glassy, polycrystalline, and liquid coenzyme Q10 studied by quasielastic neutron scattering. J. Chem. Phys. 2008, 129, 074507. (48) Rhodes, H. J.; Denardo, J. J.; Bode, D. W.; Blake, M. I. Differentiating non-aqueous titration of aspirin acetaminophen and salicylamide mixtures. J. Pharm. Sci. 1975, 64, 1386. (49) Schanker, L. S.; Shore, P. A.; Brodie, B. B.; Hogben, C. A. M. Absorption of drugs from the stomach 0.1. The rat. J. Pharmacol. Exp. Ther. 1957, 120, 528. (50) Johnson, M. R.; Frick, B.; Trommsdorff, H. P. A high resolution, inelastic neutron scattering investigation of tunnelling methyl groups in aspirin. Chem. Phys. Lett. 1996, 258, 187. (51) Prager, M.; Heidemann, A. Rotational tunneling and neutron spectroscopy: A compilation. Chem. Rev. 1997, 97, 2933. (52) Johnson, M. R.; Kearley, G. J. Quantitative atom-atom potentials from rotational tunneling: Their extraction and their use. Annu. Rev. Phys. Chem. 2000, 51, 297. (53) Ratajczyk, T.; Szymanski, S. Theory of damped quantum rotation in NMR spectra. I. Fundamental aspects. J. Chem. Phys. 2005, 123, 204509. (54) Szymanski, S. Nuclear magnetic resonance line shapes of methyllike quantum rotors in low-temperature solids. J. Chem. Phys. 1999, 111, 288. (55) Limbach, H. H.; Ulrich, S.; Grundemann, S.; Buntkowsky, G.; Sabo-Etienne, S.; Chaudret, B.; Kubas, G. J.; Eckert, J. NMR and INS line shapes of transition metal hydrides in the presence of coherent and incoherent dihydrogen exchange. J. Am. Chem. Soc. 1998, 120, 7929. (56) Bordallo, H. N.; Boldyreva, E. V.; Buchsteiner, A.; Koza, M. M.; Landsgesell, S. Structure-property relationships in the crystals of the
smallest amino acid: an incoherent inelastic neutron scattering study of the glycine polymorphs. J. Phys. Chem. B 2008, 112, 8748. (57) Weller, M.; Henry, P.; Wilson, C. An analysis of the thermal motion in the negative thermal expansion material Sc2 (WO4) 3 using isotopes in neutron diffraction. J. Phys. Chem. B 2000, 2, 12224. (58) Koshland, D. E. Application of a theory of enzyme specificity to protein synthesis. Proc. Natl. Acad. Sci. U.S.A. 1958, 44, 98. (59) Park, I. H.; Li, C. L. Dynamic ligand-induced-fit simulation via enhanced conformational samplings and ensemble dockings: A survivin example. J. Phys. Chem. B 2010, 114, 5144. (60) Savir, Y.; Tlusty, T. Conformational proofreading: The impact of conformational changes on the specificity of molecular recognition. PLoS One 2007, 2, e468.
1041
dx.doi.org/10.1021/mp400707m | Mol. Pharmaceutics 2014, 11, 1032−1041
