Article pubs.acs.org/molecularpharmaceutics
Molecular Dynamics Simulation of Amorphous Indomethacin Tian-Xiang Xiang and Bradley D. Anderson* Department of Pharmaceutical Sciences, College of Pharmacy, University of Kentucky, Lexington, Kentucky 40536, United States ABSTRACT: Molecular dynamics (MD) simulations have been conducted using an assembly consisting of 105 indomethacin (IMC) molecules and 12 water molecules to investigate the underlying dynamic (e.g., rotational and translational diffusivities and conformation relaxation rates) and structural properties (e.g., conformation, hydrogenbonding distributions, and interactions of water with IMC) of amorphous IMC. These properties may be important in predicting physical stability of this metastable material. The IMC model was constructed using X-ray diffraction data with the force-field parameters mostly assigned by analogy with similar groups in Amber-ff03 and atomic charges calculated with the B3LYP/ccpVTZ30, IEFPCM, and RESP models. The assemblies were initially equilibrated in their molten state and cooled through the glass transition temperature to form amorphous solids. Constant temperature dynamic runs were then carried out above and below the Tg (i.e., at 600 K (10 ns), 400 K (350 ns), and 298 K (240 ns)). The density (1.312 ± 0.003 g/cm3) of the simulated amorphous solid at 298 K was close to the experimental value (1.32 g/cm3) while the estimated Tg (384 K) was ∼64 degrees higher than the experimental value (320 K) due to the faster cooling rate. Due to the hindered rotation of its amide bond, IMC can exist in different diastereomeric states. Different IMC conformations were sufficiently sampled in the IMC melt or vapor, but transitions occurred rarely in the glass. The hydrogen-bonding patterns in amorphous IMC are more complex in the amorphous state than in the crystalline polymorphs. Carboxylic dimers that are dominant in α- and γ-crystals were found to occur at a much lower probability in the simulated IMC glasses while hydrogen-bonded IMC chains were more easily identified patterns in the simulated amorphous solids. To determine molecular diffusivity, a novel analytical method is proposed to deal with the non-Einsteinian behavior, in which the temporal evolution of the apparent diffusivity Dt is described by a relaxation model such as the KWW function and extrapolated to infinite time. The diffusion coefficient found for water diffusing in amorphous indomethacin at 298 K (2.7 × 10−9 cm2/s) compares favorably to results obtained in experimental IMC glasses (0.9−2.0 × 10−9 cm2/s) and is mechanistically associated with β-relaxation processes that are dominant in sub-Tg glasses. KEYWORDS: indomethacin, molecular dynamics simulation, amorphous solids, hydrogen bonding, diffusion, mobility, relaxation
■
of continuing structural relaxation after glass formation.11,12 Phase separation of amorphous mixtures13 or drug crystallization to form stable or metastable polymorphs10,14 may also occur after glass formation. These tendencies may pose significant challenges during drug product development, delay regulatory approval, and increase risks of subsequent product failure. Rates of drug crystallization from amorphous formulations are known to depend on molecular mobility.10,15−17 In amorphous drug−excipient mixtures, translational diffusion of drug molecules is a prerequisite for nucleation and crystal growth while rotational motions are necessary to optimize intermolecular interactions at the nuclear−amorphous interface.18,19 Detecting low-temperature mobility such as βrelaxation (or Johari−Goldstein20) processes and the corresponding kinetics in pharmaceutical glasses21 is difficult.
INTRODUCTION Indomethacin (1-(p-chlorobenzoyl)-5-methoxy-2-methylindole-3-acetic acid, IMC) is a powerful antipyretic, analgesic, and anti-inflammatory drug. However, in its most stable crystalline γ-form, IMC has a very low water solubility of ∼1−2 μg/mL.1 Amorphous solid dispersions of poorly soluble drugs such as IMC or its mixtures with certain polymer carriers may be useful to enhance solubility.2−5 This solubility advantage may translate into improved dissolution rates and corresponding enhancements in oral bioavailability.6 However, despite progress in the development of more sensitive analytical techniques to characterize amorphous dispersions, many of their thermodynamic and kinetic properties are not predictable and remain poorly understood. Moreover, these properties often depend on the method of preparation (e.g., melt-cooling, milling, spray drying, and vapor deposition).7−10 A particular concern associated with amorphous solid dispersions is that they are high-energy states relative to crystalline forms and therefore metastable. Thermodynamic properties such as density, enthalpy, free energy, etc. gradually drift toward their equilibrium values in such systems as a result © 2012 American Chemical Society
Received: Revised: Accepted: Published: 102
February 7, 2012 October 25, 2012 November 1, 2012 November 1, 2012 dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
history. These differences, in turn, may contribute to different stability outcomes. In this context, it is noteworthy that the γcrystalline polymorph forms from amorphous indomethacin at temperatures near or below the Tg while higher temperatures favor the metastable α-crystalline form.22 Water sorption facilitates crystallization, favoring the γ-crystalline form at low water contents and the α-polymorph at higher water content.26 Karmwar et al. found that the indomethacin polymorphs formed during crystallization were dependent on the method of preparation of the amorphous material (i.e., melt quenching, spray drying, ball milling, and cryo-milling) and that structural features identified from principal component analyses of the Raman spectra correlated with the specific polymorph that recrystallized.10 Unfortunately, unraveling molecular structure in amorphous drug solids remains a challenge as little detailed information can be gained from X-ray diffraction, the technique that is commonly used to determine crystal structures. While other techniques, including vibrational spectroscopy, NMR, shear viscosity measurements, differential scanning calorimetry, isothermal microcalorimetry, dielectric relaxation, and thermally stimulated depolarized current (TSDC),5,18,27,28 have provided useful information, a mechanistic understanding of the underlying structural properties and dynamic processes that may depend on the method of preparation and aging history of these metastable materials29 is still lacking. Molecular dynamics (MD) simulation is ideally suited to explore various atomic-level structures and certain dynamic processes that may be important in predicting properties of amorphous drug dispersions that cannot be easily investigated experimentally. Using this computational method, significant progress has been made in recent years toward understanding glass transition phase behaviors,30,31 miscibility,32−34 and molecular motions accompanying relaxation and diffusion processes in amorphous solids.35,36 In the present MD simulation study, IMC assemblies containing a small amount of water (0.6% w/w) were equilibrated in their molten state and cooled through the glass transition temperature to form amorphous solids. The newly formed IMC glasses were then subjected to long aging dynamic runs (240 ns) to investigate the underlying dynamic (e.g., rotational and translational diffusivities and conformation relaxation rates) and structural properties (e.g., conformation, hydrogen-bonding distributions, and interactions of water with IMC).
Nevertheless, residual molecular mobility persists well below the glass transition temperature (Tg) and is considered to be the reason for the slow crystallization of IMC at temperatures as low as 20 °C, well below its Tg of 42 °C.15,22 On a molecular level, the driving force for a given phase transformation is at least partly determined by differences in the molecular packing motifs and associated interaction patterns. IMC possesses a largely hydrophobic indole and chlorobenzyl group and several hydrophilic groups (i.e., the amide, methoxyl, carboxylic acid) that can act as hydrogen bond (HB) donors and/or acceptors. The formation of HB dimers between the carboxylic acid groups of adjacent IMC molecules in γ-crystals is considered to be a major factor contributing to their lower free energy while the existence of multiple HB patterns in metastable α-crystals, in which two IMC molecules form a carboxylic acid dimer while the carboxylic acid of a third molecule is hydrogen bonded to one of the amide carbonyls of the dimer, may contribute to their unusually high density.23 Due to the hindered rotation of the partial double bond between its N1 and C9 atoms (Figure 1) characterized by the
■
COMPUTATIONAL METHODS IMC molecules with different conformations were built using xLeap/Amber with the atomic numbering and associated atomic types listed in Table 1. Table 2 presents the forcefield parameters for bond length and angle associated with different atomic types that were determined based on available X-ray diffraction crystal-structure data (Cambridge Structural Database) and utilized in place of the corresponding default Amber force-field parameters (ff03). Other force-field parameters for various torsional motions and nonbonded L-J interactions were assigned by analogy with existing parameters in the Amber database (ff03) except that the torsional barrier for rotation of the chlorobenzyl group (i.e., C9−C10 bond in Figure 1) was lowered (from 29.0 to 8.0 kcal/mol) to facilitate transition between and therefore sampling of various IMC conformations. Partial charges for IMC (cf. Table 1) were obtained by first calculating the electrostatic potentials (ESPs) at the level of B3LYP/ccpVTZ30 using the IEFPCM
Figure 1. Chemical structure of IMC with atomic numbering.
angle ϕ1, IMC can exist as different geometric isomers. Different conformations predominate in different crystalline forms (e.g., γ- and α-polymorphs) suggesting that the polymorphism may stem from this conformational diversity.23−25 In γ-crystals, IMC exists exclusively as the Z isomer while three different conformations (Z, E, and α#3; labeled as M1, M4, and M2 by Aubrey-Medendorp24) coexist in metastable α-crystals.23 The disordered molecular organization and conformational diversity in amorphous IMC may lead to complex interaction patterns and conformational distributions that vary with the method of glass preparation and thermal 103
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
initially equilibrated temperature to a final temperature of 200 K at a cooling rate of 0.03 K/ps. After the completion of the cooling process, a microstructure for the newly formed glassy solid was used as a restarting configuration for a prolonged aging dynamic run (240 ns) at 298 K and 1 bar. To compare the diffusion and relaxation behaviors of the molecules in the simulated IMC glass at 298 K with those at temperatures near and above the Tg, dynamic runs were also carried out at 600 K (10 ns) and 400 K (350 ns) from equilibrated assemblies at the corresponding temperatures. System trajectories were acquired during the above dynamic runs at 0.01−0.1 ns intervals for subsequent analyses. The minimization and dynamic runs were performed using Sander 11 in which Newton’s equations of motion were evolved using the Verlet leapfrog algorithm42 with a time step of 1 fs. The dielectric constant was set at 1.0. Constant temperature and pressure were maintained by coupling to an external thermal bath. 43 Electrostatic interactions were calculated via the particle mesh Ewald method44 with a cutoff of 10 Å. The SHAKE algorithm was used to constrain all covalent bonds involving hydrogen. The ptraj programs in Amber 11, VMD 1.8.7,45 and numerous computer programs developed in-house were utilized to calculate numerically and analyze visually various structural and dynamic properties of the simulated IMC system. Calculations were performed on a Lipscomb HPC cluster (Dell Inc.) rated at >40 teraflops at the High Performance Computing complex, University of Kentucky, and an SGI workstation and PCs in the authors’ laboratory.
Table 1. Atomic Numbers, Types, and Partial Charges in IMC number
type
charge/e
number
type
charge/e
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14 C15 C16 C17 C18 C19 H1 H2
CW C* CB CA CA CA CA CN C CA CA CA CA CA CA CT CT C CT HO HA
0.1333 −0.1811 0.0755 −0.3321 0.3645 −0.3412 −0.1605 0.0232 0.7459 −0.2449 −0.0774 −0.0319 −0.0081 −0.0319 −0.0774 −0.3241 −0.0272 0.7025 0.0973 0.4666 0.1821
H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 H13 H14 H15 H16 Cl1 N1 O1 O2 O3 O4
HA HA HA HA HA HA HC HC HC HC HC HC HC HC Cl NA O OS OH O
0.1837 0.1816 0.1223 0.1184 0.1184 0.1223 0.1242 0.1242 0.1242 0.0604 0.0604 0.0481 0.0481 0.0481 −0.0987 −0.2105 −0.5422 −0.3810 −0.6002 −0.6049
Table 2. Bond Lengths and Angles for Different Atomic Types bond
length (Å)
angle (deg)
OS−CA CT−CW NA−CW C*−CW CA−HA Cl−CA CA−CA C−CA C−O C−NA CA−CB CB−CN NA−CN CT−HC CT−OS HO−OH C−OH C−CT C*−CT C*−CB
1.38 1.50 1.42 1.35 0.93 1.73 1.38 1.48 1.21 1.42 1.40 1.40 1.42 0.97 1.42 1.06 1.30 1.50 1.49 1.43
CA−CA−CA CA−CA−HA CA−C−NA O−C−OH C−OH−HO O−C−NA CA−C−O CA−CA−OS OS−CT−HC CW−CT−HC NA−CN−CB CN−NA−CW NA−CW−C* CW−C*−CB C*−CB−CN C*−CB−CA NA−CN−CA C*−CT−C C*−CW−CT CT−CW−NA CN−NA−C CW−NA−C CA−OS−CT
■
120.00 120.00 116.70 123.05 110.30 120.82 122.42 120.00 109.50 109.50 106.99 108.08 108.52 108.73 107.64 131.37 132.14 112.80 128.96 122.37 126.67 124.76 117.92
RESULTS AND DISCUSSION Formation and Aging of Amorphous IMC Glasses. The glass transition can be experimentally brought about by cooling of the melt. For example, an amorphous IMC glass has been prepared by heating the γ-form up to the temperature of 210 °C and quenching the sample into liquid nitrogen.46 A similar technique was employed in the present MD simulation, in which an IMC glass was generated by fast cooling of the IMC melt from 600 to 200 K at a cooling rate of 0.03 K/ps. As an approximation that is only empirically justified, the glass transition temperature, Tg, can be estimated from the point of intersection of straight lines drawn through the curve at temperatures above and below Tg as illustrated in Figure 2.20
continuum-solvent model37−39 in Gaussian 03, followed by fitting the ESPs with the RESP method.40,41 An assembly consisting of 105 IMC molecules combined with 12 TIP3P water molecules (0.6% w/w) was then constructed in an enlarged cubic box followed by energy minimization (100/200 iterations of the steepest descent and conjugate gradient) to eliminate bad contacts and equilibration for 10 ns at 600 K and 1 bar under periodic boundary conditions. The equilibrated assemblies were then subjected to a dynamic run during which the systems were cooled from the
Figure 2. Density−temperature phase diagram for the simulated IMC assembly. The system was cooled from 600 to 200 K at a rate of 0.03 K/ps. 104
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
Linear sections of the density−T curve were identified by conducting lack-of-fit analyses using the F-test at p < 0.05 for linear fits of the data starting at either 200 or 600 K. To obtain an estimate of pure error for the lack-of-fit tests, 2 successive data points were binned to produce 60 determinations of average density for the 60 temperature points plotted in Figure 2. This procedure produced an estimated Tg of 384 K, ∼64 deg higher than the experimental value of ∼320 K.27 Apart from the much faster (∼1011-fold) cooling rate employed here (0.03 K/ ps) compared with the experimental cooling rate of 20 K/ min,27 which has been known to increase Tg,35,36 the slight bending of the phase diagram may also contribute to the Tg disparity. Similar bending phenomena are also found in MD simulations of other materials47−49 and are perhaps attributable to (1) mismatch between the relaxation frequency of IMC and the cooling rate, (2) the presence of additional phases (e.g., rubbery phase) between the melting and glass transition temperatures, and/or (3) statistical uncertainty in the calculations. A representative 3-D image of a simulated IMC glass formed at 298 K is shown in Figure 3 (upper panel). The density of the glass at 298 K, 1.312 ± 0.003 g/cm3, is very close to the reported experimental value of 1.32 g/cm3,50−52 suggesting that by this criterion the present IMC assembly is a realistic representation of amorphous IMC glass. In the simulation, bond lengths and angles determined from X-ray diffraction crystal-structure data for IMC were used rather than the existing Amber force-field parameters (ff03). The latter values resulted in a significantly lower density (1.25 g/cm3) for a simulated IMC glass and were therefore deemed inadequate. Also shown in Figure 3 (lower panel) is the distribution of the 12 water molecules that were included in the simulated IMC assembly. At this concentration (0.6% w/w) most, but not all, of the water molecules exist as monomers. Apart from the thermodynamic changes, molecular mobility was also monitored during the cooling process as highlighted by two representative trajectories of the distance (r) between the carboxylic acid proton (−OH) in one IMC molecule and the carbonyl oxygen atom (OC) in the carboxylic acid of another IMC molecule as illustrated in Figure 4. Over the 12 ns in which the temperature was reduced from 600 to 200 K at a constant cooling rate, one IMC pair initially forming a hydrogen bond (HB) having a short −OH···OC distance of 10 Å. The other IMC pair, consisting of two wellseparated IMC molecules at the start of cooling, formed a hydrogen bond within 4 ns which remained for the duration of the cooling period. Interestingly, all these changes occurred well above the Tg (391 K) and fluctuations in −OH···OC distance apparently diminished after a cooling time of ∼5.4 ns or at ∼420 K, which was near but ∼30 deg higher than the Tg. The freezing of molecular mobility has been rationalized by a “caging effect” below a critical temperature (Tc) near but somewhat above Tg as predicted by the mode-coupling theory (MCT).53 After glass formation, structural relaxation or physical aging results in a gradual drift of the thermodynamic properties toward equilibrium.11,12 Representative results for the density and internal potential energy of the IMC glass are presented in Figure 5. Over the 100 ns simulation, the density increased by 0.073% while the total potential energy decreased by 0.092%, suggestive of relaxation of the glass toward a more stable
Figure 3. A representative microstructure of an amorphous IMC glass at 298 K (upper panel) and the corresponding water distribution in the same system (lower panel).
Figure 4. Two representative trajectories of the distance (r) between the carboxylic acid proton in one IMC molecule and the carboxylic acid carbonyl oxygen atom in another IMC molecule during the 12 ns cooling process from 600 to 200 K.
105
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
probabilities of 0.36 ± 0.13, 0.28 ± 0.07, and 0.29 ± 0.10, respectively. The carboxylic acid dimers that are dominant in αand γ-crystal forms as supported by the strong band at 1710− 1717 cm−1 in the FTIR spectrum5 were also found in the simulated IMC glass (cf. Figure 6c,d), though at a much lower probability of occurrence. HB patterns that were more easily identified in the simulated amorphous solids are chains (or aggregates) of IMC molecules connected by carboxylic acid− carboxylic acid hydrogen bonds, two representatives of which are shown in Figure 6a,b. The probabilities for finding a hydrogen bond at different HB acceptor sites including water molecules are shown in Figure 8. For the IMC−COOH hydrogen donor, ∼73 ± 6% of the HBs are at the carboxylic acid CO (O4) with the remaining at the benzoyl CO (O1) (17 ± 6%), a motif also found in the α-crystal form, the OCH3 group (O2) (1 ± 1%), or with water molecules (9 ± 2%). In contrast, water molecules acting as hydrogen donors have a similar tendency to form HBs with the benzoyl CO (O1) (39 ± 8%) and the carboxylic acid CO (O4) (37 ± 12%) with the remaining HBs at the OCH3 group (O2) (12 ± 6%) or with other water molecules (12 ± 9%). From the most probable distances in the corresponding radial distribution functions, g(r), shown in Figure 9, the HB length, a measure of HB strength, was found to be 1.540 and 1.605 Å, respectively, for the HBs formed between the carboxylic acid H1 and either O4 or O1. These bond lengths are shorter than those in the αand γ-crystalline forms (1.604 and 1.637 Å),23 suggesting somewhat stronger hydrogen bonding in the simulated amorphous solids. In contrast to the crystal forms, ∼21 ± 3% of the IMC molecules were found to not be involved in a hydrogen bond in the simulated IMC glass. Non-hydrogenbonded carboxylic CO groups have been observed as a shoulder at 1735 cm−1 in the FTIR spectrum for amorphous IMC and have been attributed to the end molecules in hydrogen-bonded IMC chains of varying length.5 A similar observation (i.e., a shoulder at 1759 cm−1) was reported in the FTIR spectrum for acetic acid.5 The formation of hydrogen-bonded IMC chains enhances the local concentration of polar carboxylic acid groups leading to spatial heterogeneity in the distribution of polar and nonpolar residues in the amorphous solid, as illustrated in Figure 6. This may have a significant impact on the location and solvation energy of water molecules in amorphous IMC. Indeed, as shown in Figure 9, a sharp peak at 1.60 Å is evident in the radial distribution function g(r) for the oxygen atoms in water molecules with respect to the IMC carboxylic acid hydrogen atom (H1), suggesting a preferential location for water near this HB group. The distribution of water with respect to the relatively nonpolar hydrogen atoms (H2−H16) in IMC molecules is roughly flat beyond the close-contact distance (2.6 Å), being somewhat less favored at short distances with g(r) ∼ 0.85 at r = 2.6−3.8 Å and approaching a g(r) ∼ 1.0 at r > 5.7 Å. The amide carbonyl atoms C9 and O1 in IMC have large and opposite charges of 0.7459 and −0.5422 (Table 1). This may induce intermolecular dipole−dipole bond formation between the amide carbonyls (CO) of adjacent IMC molecules. Indeed, a short intermolecular C9···O1 distance of 3.19 Å has been observed in the γ-crystalline form.55 However, the radial distribution function between C9 and O1 of different IMC molecules in Figure 9 reveals no significant dipole−dipole bonding in the simulated IMC glasses, reflecting perhaps the more disordered nature of the amorphous material.
Figure 5. Evolution of the density and internal potential energy at 298 K following the formation of the amorphous IMC glass.
structure. However, the changes observed are comparable to the statistical fluctuations in these properties. Intermolecular Interactions in Amorphous IMC. Hydrogen bonding (HB) plays a crucial role in determining thermodynamic (e.g., solubility and miscibility) and kinetic (e.g., nucleation and crystal growth) properties of many condensed materials.54 IMC, as shown in Figure 1, possesses a HB donor carboxylic acid group (−O3−H1) and several HB acceptor oxygen atoms including the amide carbonyl oxygen (O1), the methoxyl (O2), and the carboxylic acid (O4). Because of the more disordered nature and higher mobility in an amorphous solid, the HB patterns in IMC glasses are expected to be more complex.5,8,23 This is consistent with experimental vibrational spectra of amorphous IMC solids which exhibit broader and less intense bands than those found in the crystal forms5,8 and is attributed to a broader range of molecular conformations and intermolecular interaction arrangements in the amorphous solids. The diversity of HB structure in amorphous IMC is highlighted in Figure 6, where some representative HB patterns found in the simulated IMC glass that may resemble those in various IMC crystalline polymorphs are shown. Here a hydrogen bond was defined as one in which the protonacceptor distance is less than 3.5 Å and the angle between the proton-donor bond and the line connecting the donor proton and acceptor atoms is greater than 120 deg. In Figure 7 are shown the probability distributions for the number of HBs in which each IMC and water molecule participates in the simulated amorphous IMC. Most IMC molecules (79 ± 3%) were found to participate in at least one HB with a majority of them involved in two (40 ± 3%) or three HBs (8 ± 2%). More than 90% of the water molecules form 1−3 HBs with 106
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
Figure 6. Some representative hydrogen-bonding patterns found in the simulated IMC glass: (a), a hydrogen-bonded tetramer chain; (b) a hydrogen-bonded trimer chain; (c) a hydrogen-bonded trimer resembling that in the unit cell of the α-crystalline polymorph; (d) a hydrogenbonded dimer resembling that in the unit cell of the γ-polymorph.
Conformation Distribution and Relaxation in Amorphous IMC. Due to hindered rotation of the partial double bond between the N1 and C9 atoms (ϕ1), IMC can exist in different diastereomeric forms having varying anti-inflammatory activity with the Z isomer at ϕ1 ∼ −32° being a much more potent anti-inflammatory agent than the corresponding E isomer at ϕ1 ∼ −123°.56 The distribution of conformations is sensitive to the physical environment as demonstrated by the finding that the Z isomer is preferred in solution while the E isomer is favored in inclusion complexes with β-cyclodextrin.57 Previous quantum mechanical studies have indicated that the torsional energy for different diastereomeric forms (ϕ1) depends on the adjacent torsional angle (ϕ2) for the bond between the amide and phenyl groups, −C9−C10− (cf. Figure 1).24 This conclusion is supported by the present MD study. In MD simulations of a single IMC molecule using existing Amber force-field parameters for the torsional angle ϕ2, virtually no
interconversions were observed between the various stable diastereomeric forms at 298 K over 1000 ns. As shown in Figure 10, only when the torsional barrier for ϕ2 was reduced from 29.0 to 8.0 kcal/mol were frequent transitions between different ϕ 1 states observed at 298 K. Because such interconversions were necessary to ensure equilibration among the different diastereomeric forms, a torsional barrier of 8.0 kcal/mol was assumed for ϕ2 in the MD simulations of amorphous IMC. Because of the 2-fold symmetry of the planar chlorobenzyl group in the present molecular model, only the range of ϕ2 between −90° and 90° is plotted in Figure 10. As shown in Figure 10, the sampling points obtained in both single-molecule and amorphous solid simulations are roughly concentrated around the optimized diastereomeric structures obtained from ab initio calculations.24 The probability distribution for the dihedral angle ϕ1 is shown in Figure 11. As shown in Figure 11, the Z isomer at 107
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
Figure 9. Radial distribution functions, g(r), between the oxygen atoms in water and IMC carboxylic acid hydrogen atoms, H1 (solid circle); relatively nonpolar IMC hydrogen atoms, H2−H16 (open circle); and IMC amide carbonyl atoms C9 and O1 (solid triangle) in the simulated IMC glass at 298 K.
Figure 7. Probability distributions for the average number of hydrogen bonds in which each IMC (black) and water molecule (white) participates in the simulated amorphous IMC at 298 K.
γ- or α-crystalline forms or in amorphous IMC, but it was found to be a locally stable conformer in a recent quantum mechanical calculation by Aubrey-Medendorp et al.24 The evolution of the dihedral angle ϕ1 for a selected IMC molecule as a function of simulation time at 600 and 298 K is illustrated in Figure 12. Multiple transitions are evident for this torsion in the IMC melt, but these events occur rarely in the IMC glass, suggesting much lower local mobility in the amorphous IMC glass. This is unlikely to be due solely to the temperature difference, as fast interconversion on the NMR time scale was observed between different isomers of the methyl ester of IMC in CD2Cl2 solution at an even lower temperature57 and in the present simulation of a single IMC molecule at 298 K (cf. Figure 10). The conformational mobility can be further described by the correlation functions for the rotation of local atomic groups that are largely driven by relevant torsional potentials. Shown in Figure 13 are the rotational correlation functions C(t) for the vector connecting the C9 and Cl1 atoms in the chlorobenzyl amide (left panel) and the vector connecting the C19 and H14 atoms (cf. Figure 1) in the smaller methyl group on the methoxyl side chain (right panel). As shown in the left panel, the slow wobbling motion of the bulky chlorobenzyl amide group is far from complete over the 100 ns time frame of the simulation in the amorphous solid, which makes it difficult to predict an exact relaxation time (τ), though a rough extrapolation would suggest τ is on the order of 1−10 μs. This is in sharp contrast with the much more rapid rotation (τ ∼ 10−100 ps) of the smaller sidechain methyl group. These results are supported by recent dynamics experiments by Carpentier et al.18 Their dielectric relaxation experiments on amorphous IMC prepared via a rapid quench of the IMC liquid gave a relaxation time of 1.8 μs, which they attributed to rotation of the chlorobenzyl group, while their proton NMR experiments gave a relaxation time of 40 ps for the rotation of IMC methyl groups. Molecular Diffusion in Amorphous IMC. Molecular diffusivity is one of the most important dynamic properties determining physical stability (e.g., nucleation and crystal
Figure 8. Probability distribution for different HB acceptor sites to form hydrogen bonds with the HB donors, −COOH in IMC (black) and the H−O groups in water (white) in the simulated IMC glass at 298 K.
−32° and its nearly symmetric counterpart at 51° (M2), both of which have been found in the α-crystalline form, have similar populations in the amorphous IMC glass. The Z isomer is more favorable than the E isomer by a factor of about 5.7 corresponding to an energy difference of 1.0 kcal/mol. Our ab initio calculation for a single optimized IMC molecule at the level of B3LYP/6-311G** indicated that the Z isomer has a similar internal energy to the M2 form (0.03 kcal/mol difference) but is 0.76 kcal/mol lower in internal energy than the E isomer. Interestingly, another isomer at ϕ1 ∼ 136° (M3) was observed in the MD simulations with a population comparable to its nearly symmetric counterpart at ϕ1 ∼ −123°. This isomer has not been found experimentally in either 108
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
Figure 10. Plots of dihedral-angle pair (ϕ1, ϕ2) for a single IMC molecule in a vacuum at 298 K sampled every 2 ns for 1000 ns (left); similar plots generated by sampling all 105 IMC molecules in the simulated amorphous solid every 10 ns. The dihedral-angle pairs for various optimized IMC states obtained from ab initio quantum mechanical calculations24 are also plotted (open circles).
Figure 12. Evolution of the torsional angle ϕ1 for an IMC molecule in the simulated IMC melt at 600 K and in the amorphous glass at 298 K over time.
log plots. Clearly, these plots do not obey a linear model with a slope of 1 as predicted from the Einstein relation. Rather, like diffusion in polymer glasses,59 distinctive regions can be identified in the plots shown in Figure 14, each exhibiting a different diffusivity power law
Figure 11. Probability distributions for the torsional angle ϕ1 in the simulated IMC melt at 600 K and in an amorphous glass at 298 K.
|r(t0 + t ) − r(t0)|2 ∝ t n
growth) in amorphous solids. In fact, recent publications by Mapes19 and Swallen et al.58 suggest that, for single component supercooled liquids, molecular diffusivity is a better predictor of the rate of crystal growth than viscosity and structural relaxation time. According to the Einstein relation, D = lim Dt = lim t −∞
t →∞
|r(t0 + t ) − r(t0)|2 6t
(2)
with n varying from 0.05 to 0.95 and generally increasing with time and temperature. For water diffusion at a temperature (298 K) well below Tg, a fast increase in ⟨r2⟩ up to ∼1.5 Å2 is observed over the first 0.1 ns. This fast rise in MSD can be mostly attributed to the ballistic motion of water within a confining cavity. The MSD increases slowly between 0.1 and 1 ns having a slope of n = 0.20. Supported by representative displacement trajectories for water molecules in the amorphous IMC in Figure 15, water motions in this region (β-like relaxation) are composed mostly of temporary localization of water molecules due to dynamical arrest (caging) as also found in glassy polymers such as 1,4polybutadiene (PBD)30,35,60 with few jump motions of water molecules between cages.36 The latter is mostly diffusive, and
(1)
where the angled brackets represent an ensemble average over all different molecules of the same kind and time origins (t0), molecular diffusivity can be evaluated by monitoring the meansquared displacement (MSD), ⟨r2⟩, over a sufficiently long time (t). In this study, the MSD profiles for IMC and water molecules in the simulated IMC assemblies were calculated at 298 and 400 K. The results are presented in Figure 14 as log− 109
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
Figure 13. Rotational correlation functions C(t) for the vector connecting the C9 and Cl1 atoms in the chlorobenzyl amide (left panel) and the vector connecting the C19 and H14 atoms (cf. Figure 1) in the methyl group (right panel) in the amorphous IMC glass at 298 K.
Figure 15. Representative sections of displacement trajectories for three water molecules in the amorphous IMC at 298 K.
amorphous solid and thereby few molecules to sample for the Einstein equation; and (3) unlike large IMC molecules, whose movements are small and continuous-like, water movements in glassy IMC, as demonstrated in Figure 15, are characterized by entrapment in cavities with occasional discrete jumps between cavities. If the frequency of jumps is relatively scarce, as is usually the case in amorphous solids at low temperature,36 the occurrence of a large (forward or backward) jump could substantially alter the overall MSD profile. At 400 K, which is near but above the Tg (384 K) obtained in the present IMC simulation, n varies from 0.66 within the first ns to 0.95 over a longer period (350 ns), approaching the value of 1.0 predicted by the Einstein relation. The higher n than that at 298 K may result from more frequent jump motions due to higher thermal energy and/or involvement of cooperative motions of neighboring molecules due to the onset of αrelaxation above Tg. Crossover from local β-relaxation to diffusive α-relaxation that is responsible for the glass transition is usually found in amorphous polymers such as 1,4polybutadiene at temperatures somewhat above the Tg.30,63 Nevertheless, the exact contribution of α-relaxation to the diffusion of water at 400 K remains unclear.
Figure 14. Log−log plots of the mean-squared displacements, ⟨r ⟩, for water and IMC in the IMC system at 298 and 400 K versus time. 2
its contribution (relative to local motions in cages) to overall molecular displacement determines the degree (n) of nonEinsteinian behavior. Above 10 ns, the MSD increases more rapidly with n = 0.74, reflecting occurrence of more diffusive jump events over a longer period. Near the end of the simulation (t > 230 ns), a slight falloff in the MSD profile for water molecules at 298 K was noted. Similar phenomena, which are attributed to increasing statistical errors at long time intervals and are usually discarded because of their unreliability in predicting diffusivity, have been observed previously in numerous MD calculations of diffusion coefficients of small solutes in amorphous solids.23,61,62 The statistical errors arise mainly from the following factors: (1) there are fewer sampling origins (t0 in eq 1) at long time; (2) amorphous IMC has a low water content (0−2% w/w), and the value of 0.6% w/w water used in this study provides only 12 water molecules compared to 105 IMC molecules in the 110
dx.doi.org/10.1021/mp3000698 | Mol. Pharmaceutics 2013, 10, 102−114
Molecular Pharmaceutics
Article
Figure 16. Apparent diffusivity Dt for water (left panels) and IMC (right panels) in the IMC glass formed at 298 K (upper panels) and 400 K (lower panels) versus time.
this KWW model to describe the Dt decay profiles in Figure 16, we have
For IMC molecules, as shown in Figure 14, the general trend of changes in non-Einsteinian diffusion is similar to that for the much smaller water molecules, in that there exists a region exhibiting a slow rise in MSDs followed by an upturn after ∼10 ns depending on temperature. At 298 K virtually all of the IMC molecules remained confined within a very small volume during the entire simulation time as indicated by the small MSDs (
