Effects of Solvent Stabilization on Pharmaceutical Crystallization

May 17, 2019 - The thermodynamics of solvent inclusion in the Form II·MeOH ... Department of Chemistry and Biochemistry, Boise State University, 1910...
1 downloads 0 Views 3MB Size
Download PDF
Article Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX

pubs.acs.org/JPCA

Effects of Solvent Stabilization on Pharmaceutical Crystallization: Investigating Conformational Polymorphism of Probucol Using Combined Solid-State Density Functional Theory, Molecular Dynamics, and Terahertz Spectroscopy Neilson R. Rexrode,† Jordan Orien,‡ and Matthew D. King*,†,‡ Micron School of Materials Science and Engineering and ‡Department of Chemistry and Biochemistry, Boise State University, 1910 University Drive, Boise, Idaho 83725, United States

Downloaded via UNIV OF ROCHESTER on May 19, 2019 at 10:00:36 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: Solid-state density functional theory (DFT), molecular dynamics (MD), and terahertz (THz) spectroscopy were used to study the formation of enantiotropically related conformational Form I and Form II polymorphs of the pharmaceutical compound, probucol. DFT calculations were performed on the crystal systems to compare relative lattice energies and the solvent stabilization of the metastable Form II structure. The thermodynamics of solvent inclusion in the Form II·MeOH crystal system were determined from MD simulations, as was the favored conformation of molecular probucol in methanol and ethanol solutions. The findings from both solid-state DFT and MD calculations suggest that the preferred molecular orientations of the probucol molecule in solution and the probable inclusion of methanol in the crystal lattice during the crystallization process lead to the solvent selectivity of the probucol polymorph formation. The additional stabilization energy provided by the crystallization solvent facilitates the nucleation and growth of the Form II polymorph under conditions that favor this metastable crystal form over the thermodynamically stable Form I, despite the higher energy molecular and crystalline configurations of probucol Form II. We demonstrate the influence of solvent on the formation of pharmaceutical polymorphs and provide a molecular-level view of complex interactions leading to polymorphism using a combination of computational methods and THz spectral data.

■

different carrier solutions.9 MD has also been employed to provide insight into the process of crystallization of different polymorphs.10−12 DFT and MD can be utilized as a powerful combination of computational tools for the study and prediction of crystalline structures and properties. This study investigates two known polymorphic forms of a pharmaceutical compound, probucol, to uncover the underlying physical processes influencing solvent selectivity in polymorph formation.13 Probucol, or 4,4′(isopropylideneithio)bis(2,6-di-tert-butylphenol), was developed for the treatment of coronary heart disease and remains the subject of research into its physical and pharmacological properties.14−18 Two primary crystal forms of probucol (referred to as Form I and Form II, herein) have been identified and investigated in detail in the literature.13,19 Form I is the thermodynamically stable polymorph with a melting point of 125 °C, compared to 116 °C for the metastable Form II polymorph. Form II will also convert to Form I over time at

INTRODUCTION Polymorphism and solvatomorphism are two major concerns in today’s pharmaceutical industry.1 Different crystalline structures of the same medicinal compound can exhibit significantly varying properties, such as stability, solubility, dissolution rate, hygroscopicity, bioavailability, and therapeutic index.2,3 Such properties can have a profound influence on the efficacy of a drug, making it essential to develop effective methods for the characterization and prediction of crystalline structures and properties of an active pharmaceutical ingredient (API). Computational modeling approaches using density functional theory (DFT) and molecular dynamics (MD) are often implemented in the study of polymorphism and solvatomorphism to provide a better understanding of the structure and properties of such systems. It has been demonstrated that DFT is capable of reliably predicting the structures and vibrational terahertz (THz) spectra of pharmaceutically relevant polymorphs, such as those formed by ibuprofen,4 diclofenac acid,5 aripiprazole,6 salicylic acid,7 and irbesartan.8 Likewise, it has been shown that MD may be used to predict the solubility of drug molecules, such as indomethacin, in © XXXX American Chemical Society

Received: January 25, 2019 Revised: May 6, 2019

A

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A room temperature.19 The probucol crystal Forms I and II are considered constitutional polymorphs and are enantiotropically related as a transition state is present between the two polymorphic crystal forms.20 As the physical properties of pharmaceutical polymorphs are governed by their molecular packing arrangements, it is imperative to understanding the relationships among the molecular structure, crystal structure, and molecule−solvent interactions to aid the drug development process. The probucol molecular structure and the unit cells of the Form I and II polymorphs are shown in Figure 1. Both

frequency THz spectra indicates that our applied computational methodology for computing the potential energy surface is suitable for determining the energetics and crystal properties. The results of our DFT and MD calculations reveal the role of solvent molecules in the formation of metastable polymorphs via molecular conformational preferences because of solvent interactions and through the energetic favorability of solvent inclusions within the crystal lattice.

■

MATERIALS AND METHODS Solid-State DFT. The structures and THz vibrational spectra of probucol Forms I and II were calculated by solidstate DFT using the Crystal14 software package.21,22 Published X-ray crystal structures of probucol available from the Cambridge Crystallographic Database were used for the initial atomic positions and lattice dimensions for the computational model (HAXHET and HAXHET01 for Forms I and II, respectively).13,23 Initial DFT structural calculations utilized the generalized gradient approximation (GGA) density functional Perdew−Burke−Ernzerhof (PBE),24 the hybrid meta-GGA functional M06,25 and the global hybrid functional B3LYP,26,27 with all calculations using an atom-centered 631G(d,p) basis set.28,29 On the basis of the resulting qualities of structural reproductions, only PBE/6-31G(d,p) was used in subsequent frequency and structural calculations. A shrinking factor of 5 for both polymorphs was used in defining the kpoint sampling of the density matrix and commensurate grid in reciprocal space to achieve a self-consistent field convergence precision of 10−8 hartree for geometry optimizations and 10−11 hartree for normal mode calculations.30,31 Convergence on the root-mean-square (rms) of the energy gradient of 10−4 and 10−5 hartree were applied for geometry optimizations and frequency calculations, respectively. A pruned (75 974) integration grid was used to define the radial and angular distribution of grid points. A pairwise empirical energy correction for London-type dispersion forces as incorporated in the Crystal14 software package was included in all DFT calculations. The global scaling factors, s6, used for the dispersion correction term were 0.75, 0.25, and 1.05 for PBE, M06, and B3LYP, respectively.25,32−34 Additional atomic dispersion parameters, including C6 coefficients, van der Waals radii, and cutoff distance, were taken from the parameterization by Grimme et al.32,33 Geometry optimizations of atomic coordinates and lattice dimensions were performed within the constraints of space group symmetry. Energy corrections for the basis set superposition error (BSSE) were performed using the counterpoise method.35 Normal mode frequencies were calculated within a harmonic approximation by numerical differentiation of the analytical gradient of the potential energy with respect to the atomic positions.36 The infrared intensities for normal modes were calculated analytically using a first-order coupled perturbed Hartree−Fock approach for calculating Born charges.37−39 MD Simulations. For MD simulations, the GROMACS v5.0.4 software package40−42 was used with an optimized potentials for liquid simulations (OPLS)-AA force field parameter set.43 Developing a suitable molecular mechanical description of organic molecules in the crystal phase remains a challenge, and commonly used force fields have been employed for studies of such systems.44−48 The OPLS-AA force field was chosen for this study because of its initial parameterization to reproduce small organic molecule properties and for the availability of the existing models for

Figure 1. Molecular structure, conformation, and crystal structures of probucol. (A) Molecular structure of probucol. (B) Form I polymorph is the thermodynamically favored crystal form and contains the lowest energy molecular conformer. (C) Form II polymorph consists of the higher energy molecular conformer in a metastable crystal form.

polymorphs exist in a monoclinic crystal setting with four molecules per unit cell. There is a marked difference in molecular conformation between the two forms, with a lower energy conformer present in Form I. Both forms are held exclusively by weak dispersion forces despite the presence of hydroxyl groups in the probucol molecular structure. Preparation of the individual forms can be achieved using a variety of solvent systems through evaporation, cooling, or precipitation. Generally, crystals of the stable Form I are obtained by slow evaporation in most common solvents (e.g., acetone, ethanol, and chloroform), allowing the probucol molecules to reach their thermodynamically favored molecular and crystalline conformations.13 Form II can be readily crystallized by rapid cooling of a saturated probucol solution or by precipitation with the addition of an antisolvent such as water.19 The rapid crystallization process yields the apparent kinetically favored Form II polymorph. However, Form II can also be obtained by a slow evaporation of methanol, indicating thermodynamic favorability for the Form II molecular conformation in this solvent system. The polymorphic forms of probucol were investigated using a combination of DFT, MD, and THz spectroscopy to provide answers as to the formation of crystalline products as a result of different dissolution processes, with a particular interest in the energetics of the Form II molecular conformation and crystal system. The acquisition of THz spectra allows for the identification of various polymorphs of a pharmaceutical solid and also provides valuable experimental data to compare and validate the computational results. The ability to reproduce the experimental crystal structures and lowB

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. Percent Change of Lattice Parameters between the Calculated and Experimental X-ray Crystal Structuresa Form I

% change

parameter

exp.b

PBE

M06

B3LYP

PBE

M06

B3LYP

a (Å) b (Å) c (Å) β (deg) V (Å3) Form II

16.972 10.534 19.030 113.66 3116.0

16.5315 10.1347 18.4259 109.529 2909.5

16.4601 9.8444 18.0603 113.067 2692.5

16.6680 10.0477 18.0207 113.585 2765.9

−2.60 −3.79 −3.17 −3.36 −6.6

−3.02 −6.55 −5.10 −0.52 −13.6 % change

−1.79 −4.62 −5.30 −0.07 −11.2

parameter

exp.a

PBE

M06

B3LYP

PBE

M06

B3LYP

a (Å) b (Å) c (Å) β (deg) V (Å3)

11.226 15.981 18.800 104.04 3272.0

11.0263 15.4038 18.2482 102.431 3026.7

10.8897 15.1758 17.8994 102.359 2889.5

10.9568 15.2372 17.9595 102.130 2931.4

−1.78 −3.61 −2.94 −1.55 −7.5

−3.00 −5.04 −4.79 −1.62 −11.7

−2.40 −4.65 −4.47 −1.84 −10.4

a Values for the lattice dimensions (Å), β angle (°), and volume (Å3) calculated using the PBE, M06, and B3LYP functionals and compared to the experimental values determined via published X-ray crystal structures. bExperimental parameters from ref 13.

developing small-molecule topologies.49 The molecular topologies for probucol, methanol, and ethanol were initially generated using the LigParGen server.49 Atomic charges were calculated from Mulliken population analyses at the DFT/ PBE/6-31G(d,p) level. 50 Topology files used in MD simulations for probucol, methanol, and ethanol are available in the Supporting Information. A 3 × 3 × 3 supercell of the Form II·MeOH solvate was constructed for the MD simulations from the available Form II crystal structure. The system was prepared for simulations by performing a conjugate gradient energy minimization, followed by equilibrations in NVT (constant number of particles, volume, and temperature) and NPT (constant number of particles, pressure, and temperature) ensembles. A 100 ps equilibration was conducted under the NVT ensemble using a v-rescale thermostat at 300 K with a coupling constant of 0.1 ps.51 The subsequent NPT equilibration was performed for 100 ps using an anisotropic Parrinello−Rahman barostat with a time constant of 1.0 ps, external pressure of 1.01325 bar, and isothermal compressibility of 1.2 × 10−5 bar−1.52 The integration time steps for both equilibration steps were 1.0 fs. An initial 10 ns MD simulation was executed using 2.0 fs time steps. All bonds were constrained to equilibrium values using the LINCS algorithm.53 The Verlet cutoff scheme was used in the calculation of short-range Coulomb and van der Waals forces, and particle mesh Ewald was used for long-range electrostatics.54−56 The free energy of methanol inclusion in the Form II· MeOH crystal was determined by a perturbative method calculating the energy differences as the interaction energies of a single bound methanol molecule were iteratively turned off in two steps involving three states: A, B, and C. State A represents the methanol molecule with full interaction potentials, state B corresponds to an intermediate conformation with the partial atomic charges of the methanol molecule set to zero, and state C represents the methanol molecule with all Coulomb and van der Waals interactions set to zero. In the first step, the partial charges were decreased by varying the coupling parameter, λ, from λ = 0 to λ = 1 through 20 evenly spaced intermediate λ values. Next, the van der Waals radii were decreased, again by varying the coupling parameter from λ = 0 to λ = 1 through 20 intermediate steps. At each λ value, simulations were performed for 1 ns at 1.0 fs time steps. ∂H/∂λ was saved

every 10 fs for postprocessing and free-energy calculations using Bennett’s acceptance ratio perturbation method.57 THz Spectroscopy. THz spectra were recorded using a custom time-domain pulsed THz spectrometer. A diodepumped solid-state Nd:YVO4 laser (Coherent, Verdi-5) was used to pump a Ti:sapphire oscillator (Coherent, Mira-SPO). The 800 nm output pulse seeded a Ti:sapphire regenerative amplifier (Coherent, Legend-HE-USP) pumped by a Qswitched solid-state Nd:YLF laser (Coherent, Evolution-30). The output beam from the regenerative amplifier consisted of ∼800 nm pulses with a pulse duration of 40 fs and energy of 2.5 mJ/pulse. The near-infrared (NIR) pulses were focused on a 1 mm thick 110-cut ZnTe crystal to generate ∼200 fs THz pulses with a bandwidth of 0.1−2.8 THz by optical rectification.58 The THz pulses were directed onto the sample by a pair of off-axis parabolic mirrors. The samples and blanks for measurement were held in a cryostat (Janis Research, VPF ST-100) equipped with 3 mm thick polymethylpentene windows. The transmitted THz radiation was directed onto a ZnTe detector crystal by a second pair of off-axis parabolic mirrors and was detected by free-space electro-optic sampling.59 The polarization of the reflected NIR gate pulse was monitored by a balanced photoreceiver (New Focus, 2307) connected to a lock-in amplifier (Stanford Research Systems, SR-830). The amplifier was locked to an optical chopper (Digirad, C-980) synchronized to a subharmonic of the regenerative amplifier repetition rate. A linear delay stage (Newport Corporation, ILS150CCHA) was used to vary the arrival time of the THz pulses relative to the NIR pulses at the detector crystal, yielding the time-domain THz waveform. The samples for THz measurements were mixed with polytetrafluoroethylene (PTFE) powder at concentrations of approximately 8% by mass. The samples were pulverized using a stainless steel amalgamator (Dentsply Rinn 3110-3A) to minimize the particle size and scattering of THz radiation. Approximately 0.55 g of the sample mixtures were pressed into pellets under a measured pressure of 2000 psi using a hydraulic press (ICL EZ-Press 12) equipped with a 13 mm stainless steel die. The PTFE pellets for use as reference blanks were prepared in the same manner. Data were acquired at 293 and 78 K with the samples under vacuum. The samples and blanks were scanned 32 times, and the data were averaged for each individual data set. A 32 ps C

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 2. rmsds for Molecular Parameters between the Calculated and Experimental X-ray Crystal Structuresa Form I bond lengths (Å) bond angles (deg) torsion angles (deg)

Form II

PBE

M06

B3LYP

PBE

M06

B3LYP

0.01629 0.6909 1.442

0.01049 0.8945 3.925

0.01297 0.7538 1.935

0.01819 0.7538 1.398

0.01094 0.8734 44.31

0.01439 0.8301 2.336

a

rmsds for bond lengths (Å), bond angles (deg), and dihedral angles (deg) for probucol molecules in the Form I and II crystal structures produced from calculations with the PBE, M06, and B3LYP functionals compared against the published X-ray crystal structures.

0.670 Å, respectively. Likewise, the Form II structure was best reproduced by the PBE functional, having an rmsd of 0.475 Å. The resulting rmsds from the M06 and B3LYP optimizations were 0.724 and 0.666 Å, respectively. The quality of the optimized probucol structures was also assessed in terms of rmsds in bond lengths, angles, and dihedrals containing heavy atoms (nonhydrogen atoms) in comparison to the experimental X-ray structures (Table 2). Complete lists of the experimental and calculated bond lengths, bond angles, and dihedral angles can be found in the Supporting Information. With regard to bond lengths, each functional resulted in a low rmsds for both polymorphs. M06 performed the best with the rmsds of 0.01049 and 0.01094 Å, respectively, for Forms I and II. The PBE functional exhibited the largest rmsds of 0.01629 and 0.01819 Å for Forms I and II, respectively. The opposite trend was obtained for bond angles, with PBE being the best performer followed by B3LYP and M06 for both polymorphs. Again, the rmsds for bond angles resulting from each functional were of high quality. Much greater differences in structural reproduction between functionals were observed in the evaluation of dihedral angle rmsds. For both polymorphs, the PBE calculations resulted in the lowest rmsds of 1.442° for Form I and 1.398° for Form II. The structures calculated with B3LYP had slightly higher rmsds. However, the M06 functional performed quite poor in terms of reproduction of dihedral angles. For Form I, rmsd was a reasonable 3.925°, whereas for Form II the resulting dihedral angle rmsd was 44.31°. The largest structural deviations predominantly involved the rotation of the isobutyl moieties. Most notable was the >10° rotation of the isobutyl group of the central C10 atom (see Supporting Information for atom labeling). It was concluded from the results of the structural optimizations that the PBE functional was best suited for modeling the probucol polymorphic crystal systems. This density functional has been likewise shown to perform well with the structure and vibrational frequency calculations of smaller pharmaceutical systems.4,5,61 As such, the PBE functional was used exclusively for all subsequent structure, energy, and normal mode calculations. It is important to note that GGAs, such as PBE, not only are computationally less demanding than hybrid and meta-GGA density functionals, but they also commonly outperform these more intensive functionals in terms of accuracy in the frequencies of optical phonon modes, even when used with the relatively modest double-zeta 6-31G(d,p) basis set.62−65 Given the significant computational resource requirements for solid-state DFT calculations of large crystal systems, the ability to apply a favorable density functional/basis set combination for the analysis of such systems is noteworthy. It is imperative that the applied computational model be able to accurately predict the lattice vibrational mode frequencies as low-frequency phonon modes contribute most

scan window consisting of 3200 data points was used to capture the THz waveform, which was then symmetrically zero-padded to a total of 6000 data points for data transforms. The effective instrument resolution was approximately 1.0 cm−1. The ratio of the power spectra obtained from the Fourier-transformed data sets of the sample and blank produced the THz absorption spectrum. Each THz spectrum presented in this work is the average of four individual THz spectra, each representing a complete set of sample and blank measurements.

■

RESULTS AND DISCUSSION Although solid-state DFT has been widely used in recent years in the computational analyses of biopharmaceutical molecules and small APIs, there exists little literature precedent for the application of atom-centered solid-state DFT calculations on systems as large (83 atoms per molecule) as the crystalline probucol polymorphs.8,60,61 This is primarily because of the considerable computational resources required to perform such calculations. As an example, a single calculation of vibrational frequencies for one of the probucol systems required an excess of 75 000 cpu h (using the PBE/631G(d,p) density functional/basis set combination). This otherwise intractable resource demand, when considering the multiple structural optimizations and frequency calculations, can be alleviated in a highly parallel computing environment and software capable of efficiently scaling large jobs across multiple computing nodes. For the calculations presented in this study, individual jobs typically utilized 224 Intel Xeon 2.4 GHz computing cores across eight compute nodes. Structural optimizations were initially performed using the PBE, M06, and B3LYP density functionals to evaluate their performance in the structural reproduction of the probucol crystal systems. Comparisons of the calculated unit cell dimensions with the experimental X-ray structure are provided in Table 1. For the DFT-optimized structures, contraction along all lattice vectors was observed with each density functional. This was an expected result for two reasons: (1) the experimental crystal dimensions were acquired at room temperature, whereas the calculated structures are effectively computed at 0 K, thus eliminating any thermal expansion; (2) the usage of a finite 6-31G(d,p) basis set results in a small artificial increase in the cohesive energy of molecules within the unit cell because of BSSE. The amount of BSSE was calculated at 16.6 and 17.3% of the total lattice energy for Forms I and II, respectively. However, the resulting unit cell structures exhibited reasonable reductions in volume for the applied computational parameters, with Form I undergoing a 6.6−13.6% decrease in volume and Form II with 7.5−11.5% decrease. The best reproduction of the Form I crystal structure was obtained with the PBE functional, producing an rms deviation (rmsd) in lattice vector lengths of 0.489 Å, compared to the M06 and B3LYP functionals with rmsds of 0.748 and D

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

systems, but rather artifacts arising from data collection and signal processing. The low-frequency vibrational modes and IR intensities for both probucol polymorphs were calculated by solid-state DFT. High-quality calculated THz spectra for probucol Forms I and II were obtained in good agreement with the experimental THz vibrational spectra (Figure 3). Nearly all experimental

appreciably to the zero-point vibrational energies of crystal systems. To properly evaluate crystal energies, the computational data must utilize the experimental vibrational spectra for validation of the reliability of vibrational frequency calculations. In this study, THz spectroscopy was used to acquire low-frequency vibrational dynamics of the crystalline probucol polymorphs for comparison with the predicted phonon mode frequencies. The THz spectra of probucol Forms I and II obtained at 78 K are shown in Figure 2. The exclusive

Figure 2. THz vibrational spectra of probucol Forms I and II. The THz spectra from 10 to 90 cm−1 of probucol obtained at 78 K display unique vibrational profiles between the polymorphic forms as a result of differences in crystalline packing.

vibrational mode frequencies and the absence of common peaks between the two spectra indicate the high purity of the individual probucol polymorphs grown under different solvent conditions. There are 12 discernible spectral features in the Form I spectrum centered at 31.4, 36.1, 41.0, 46.5, 50.4, 53.2, 58.6, 62.2, 66.6, 69.9, 79.3, and 83.1 cm−1. Nine vibrational mode locations can be readily resolved. The feature at 46.5 cm−1 appears as a shoulder convolved into the peak at 50.4 cm−1. The modes centered at 31.4 and 41.0 cm−1 have very low absorption strengths but can be effectively identified as real features above the baseline. Additionally, the offset in the baseline adjacent the 36.1 cm−1 peak indicates absorption contributions from the underlying vibrational modes. Peak centers were determined through fitting of individual baselinecorrected peaks with Lorentzian functions. The average width of Form I spectral absorptions at full width at half-maximum (fwhm) was 3.4 cm−1. The probucol Form II spectrum contains 11 features, although not as clearly visible as in the Form I spectrum (Figure 2). Slightly broader absorption line shapes of 5.0 cm−1 at fwhm and the close proximity of some features result in most of the peak locations able to be determined only through peak fitting. The three absorptions at 29.8, 74.3, and 82.8 cm−1 are clearly resolved, with the remaining peaks at 26.3, 37.3, 44.4, 53.4, 58.0, 62.4, 67.3, and 77.0 cm−1 identified through a combined multiple peak fitting of the full THz spectrum, truncated at the base of the peak centered at 82.8 cm−1, and the strong absorption extending beyond the instrument detection bandwidth. It should also be noted that the apparent absorptions below 20 cm−1 in both the Form I and Form II spectra are not real spectral features of the probucol crystal

Figure 3. Comparison of the probucol experimental THz vibrational spectra and those calculated by solid-state DFT. Overlay of the experimental and calculated THz spectrum for (A) Form I polymorph and (B) Form II polymorph. Vertical bars indicate the calculated vibrational modes and normalized IR intensities. Lorentzian line shapes were convolved into calculated modes using a fwhm of 3.4 and 5.0 cm−1 for Forms I and II, respectively, determined by fitting experimental data.

absorption peaks in the complex vibrational spectra could be assigned to the calculated normal modes of vibration (Table 3). Mode assignments to the experimental spectral features were made via inspection of peak location and intensity and provided as a direct comparison of experimental and calculated results. For the graphical representation in Figure 3, the calculated Form I mode frequencies were scaled by 0.88 to better align the experimental and calculated spectra. A slight overestimation of phonon vibrational frequencies is commonly observed in solid-state DFT calculations.64,66−68 In addition, the full unit cell optimizations yielding contracted lattice vectors typically result in the vibrational modes being shifted to high frequency because of the steepening of vibrational potentials. The calculated Form II spectra displayed in Figure 3B were not scaled accordingly. E

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 3. Experimental and Calculated Vibrational Frequencies of Probucol Forms I and II Phonon Modes Form I a

Form II

exp.

PBE/6-31G(d,p)

freq.b

freq.

int.c

22.34 28.48 32.09 33.21 35.06 37.03 39.31 43.20 48.58 50.08 55.26 55.43 61.80 62.72 63.01 67.80 71.16 74.99 80.36 82.04 85.75 85.93 91.08 93.29 98.18 101.87 107.49 107.86

0.68 0.43 1.13 0.69 0.04 0.73 0.40 0.70 2.02 0.49 2.36 2.72 3.59 2.31 1.06 0.02 0.82 1.03 2.50 1.43 2.55 0.59 1.94 7.65 3.11 4.82 0.07 5.77

31.4 36.1

41.0 46.5 50.4

53.2

58.6 62.2 66.6 69.9 79.3 83.1

exp.

PBE/6-31G(d,p)

freq.

freq.

int.

14.68 26.85 29.65 31.10 38.87 41.04 43.64 45.39 49.26 51.50 52.41 56.20 59.59 60.05 60.98 61.85 65.96 70.68 70.92 72.01 72.80 80.38 84.70 84.71 85.87 85.92 91.55 96.77 96.79 109.17 109.75

0.48 2.13 0.03 0.79 0.86 0.03 1.08 0.01 0.17 0.22 1.74 0.91 1.86 1.66 0.04 0.01 0.63 0.09 1.00 1.81 6.21 0.16 0.01 0.51 1.42 3.78 2.48 1.98 3.98 4.15 0.76

29.8

37.3

44.4 53.4 58.0

62.4 67.3

74.3 77.0 82.8

a Experimental. bFrequency (cm−1) (calculated frequencies for Form I are scaled). cIntensity (km mol−1).

Upon fitting of the experimental spectra, the widths of the Lorentzian fitting functions were determined to have average fwhm values of 3.4 and 5.0 cm−1 for Forms I and II, respectively. It was unexpected for similar crystal systems to exhibit substantial variations in spectral line widths. This discrepancy may be attributed to either anharmonicity, which is unlikely because of the symmetry of the experimental line shapes, or possibly to solvent molecules entrapped in the crystal in a random, or nonperiodic, manner leading to inhomogeneous broadening of spectral features. However, the goodness-of-fit of the calculated spectrum of Form II with the THz vibrational spectrum indicates that the crystal symmetry remains unbroken despite the possible presence of disperse methanol molecules contained within the crystal. To begin investigating the influence of the crystallization solvent on polymorph selection, the differences in crystal packing between Forms I and II were examined. The packing efficiency and crystal void spaces were mapped for the DFToptimized unit cells of probucol Forms I and II using CrystalExplorer17 (Figure 4).69 For these measurements, the unit cells were expanded by plus/minus 5 Å on each lattice vector. A probe sphere was rolled over a fused-sphere representation of the molecular crystal at an isovalue of

Figure 4. Measurement of crystal void volumes. Calculated isosurfaces for the optimized, computed unit cells of (A) Form I polymorph and (B) Form II polymorph. (C) Position of methanol molecules placed in regions of large voids in the Form II crystal structure, as indicated by the green circles.

0.002 e au3. This approach locates the isosurface of the estimated procrystal electron density and allows for the visualization and calculation of crystal void volumes. The thermodynamically favored Form I unit cell was found to have a packing efficiency of 83.7%, compared to the Form II structure having a lower packing efficiency of 78.0% (Table 4). The Form II crystal structure, however, had a lower packing efficiency of 78.0%. The additional void space was highly localized and associated with the positioning of each symmetry-equivalent probucol molecule of the unit cell (Figure 4B). Each local void volume was determined to be sizable enough to support the inclusion of a methanol F

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 4. Void Volumes and Lattice Parameters for Probucol Form I, Form II, and Form II·MeOH Solvate 3

total volume (Å ) void volume (Å3) packing efficiency a (Å) b (Å) c (Å) β (deg)

Form I

Form II

Form II·MeOH

2909.5 473.2 83.7% 16.5315 10.1347 18.4259 109.529

3026.7 665.9 78.0% 11.02632 15.40380 18.24815 102.4309

3006.8 460.0 84.7% 11.02718 15.36955 18.16789 102.4441

Figure 5. Root-mean-squared fluctuation of molecular center of masses for probucol and methanol molecules in the MD simulation cell. The red line indicates a 10-point running average to guide the eye.

molecule, which has an approximate molecular volume of 25.8 Å3. On the basis of the measured void volumes of the Form II polymorph, four symmetrically equivalent methanol molecules (Z′ = 1) were inserted into the Form II crystal structure in the locations indicated in Figure 4C. The generated Form II· MeOH solvate crystal structure was then optimized using solid-state DFT. The void volumes of the resulting Form II· MeOH solvate unit cell were measured as previously described (Table 4). A slight decrease of 19.9 Å3 in total cell volume was observed with the addition of methanol to the Form II structure, as well as a significant improvement in packing efficiency at 84.7%. This suggests an increase in cohesive interactions between the contents of the unit cell, which was verified by comparing the relative lattice energies for the Form II and Form II·MeOH systems (Table 5). As Form I is known to be the most stable polymorph, the resulting 72.7 kJ/mol difference in energy between Form I and Form II was expected. However, it was found that the Form II·MeOH crystal system has a calculated lattice energy of 97.94 kJ/mol lower than that of Form I, indicating an enhanced crystal stability upon the inclusion of methanol solvates. This, however, only offers a static picture of the electronic energy, neglecting thermal and kinetic energy contributions because of nuclear motion, but does support the increased stability of the crystal system from a purely structural standpoint. MD simulations were employed to further investigate the solvent interactions and energy contributions to the crystallization of the probucol Form II polymorph. An initial 10 ns MD simulation of the Form II·MeOH system was performed, and the center-of-mass motion for probucol and methanol molecules was evaluated. Figure 5 shows the rms fluctuation (rmsf) of each molecule contained in the simulation cell. It is evident from the simulation that the relative motion of the methanol molecules is notably greater than that of probucol, with rmsfs of 0.572 (±0.037) and 0.429 (±0.048) Å, respectively. This is to be expected given the disparity in molecular weight between the two molecule types. The simulation also shows that although the methanol molecules are confined to a single void space within the lattice structure, they are able to flip orientation by 180°, which would lead to

solvent disorder or noncrystallinity of methanol molecules within the solvated crystal. The estimation from MD simulations of the rate at which methanol molecules are able to flip orientation is approximately 1.5 flips/ns per methanol molecule (Figure S2, Supporting Information). The large fluctuations in position also increase the probability of solvent migration through the network of crystal voids shown in Figure 4B. Solvent migration would lead to the rapid conversion of the Form II·MeOH solvate to the pure Form II polymorph. Applying the MD free-energy perturbation simulations to the Form II·MeOH solvate system, the ΔG value for inclusion of methanol in the crystal structure was calculated to be −2.74 kJ/mol (Figure 6). A favorable ΔG suggests a net increase in

Figure 6. Free-energy changes over MD perturbative decoupling of methanol from the Form II·MeOH solvate crystal structure, yielding a ΔG value of 2.74 ± 0.16 kJ mol−1.

crystal stability with methanol incorporated into the crystal structure despite the more substantial molecular motion observed of the methanol molecules. The value of −2.74 (±0.16) kJ/mol (the reverse process yielding a positive ΔG for decoupling shown in Figure 6) is unexpectedly low considering the magnitude of change in lattice energy observed from the

Table 5. Calculated Energies (kJ mol−1) for Form I, Form II, and Form II·MeOH Solvate Form I Form II Form II·MeOH

Eunit cell

Esingle moleculea

Elatticeb

Erelc

−2.2636206 × 107 −2.2636107 × 107 −2.3850533 × 107

−5.6587613 × 106 −5.6587549 × 106 −5.9623187 × 106

−1.1603814 × 103 −1.0877576 × 103 −1.2583209 × 103

97.9 170.6 0.0

a

Energy of a single probucol molecule in the unit cells of Forms I and II and the sum of the individual probucol and methanol molecules in Form II·MeOH. bDifference between Eunit cell and Esingle molecule. cEnergies relative to minimum value. G

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A DFT calculations with the inclusion of methanol. This difference is attributed to the thermal and kinetic energy contributions, and these findings do support the enhanced energetic stability of the Form II·MeOH system. However, the low relative free-energy contribution along with the large center-of-mass motion would suggest a metastable crystal phase, resulting in the probable loss of methanol from the crystal lattice, yielding the Form II structure. Subsequent MD simulations were performed on a single probucol molecule contained within a solvation box of either pure methanol or ethanol. Over the simulation time, the preferred molecular conformation of the probucol molecule about the central S−C−S−C dihedral angle was monitored. Prior to the simulations, the torsion energy was parameterized to provide an accurate representation of the true torsion potential. A potential energy scan of the dihedral was performed using DFT/PBE/6-31G(d,p) (Figure 7A). The resulting potential energy curve was fit to the Ryckaert− Bellemans potential, and the resulting coefficients were applied to the OPLS-AA force field and the probucol topology describing bonding potentials. A simulation box was constructed around a single probucol molecule with dimensions of 35.8 Å × 35.8 Å × 25.3 Å, which should be large enough to allow the molecule to move freely without interaction with the molecules in a neighboring cell (as simulations are performed in periodic boundary conditions). The void space was filled with methanol or ethanol molecules by random insertion as to not overlap the van der Waals radii of the central probucol molecule. For the methanol system, the number of solvent molecules added was 571, and for the ethanol system, this number was 308. Figure 7B shows the S− C−S−C dihedral angle distributions for the methanol and ethanol systems for the initial simulation time period of 0−100 ps from an initial molecular conformation having a 0° dihedral angle, such as that found in the Form II crystal structure (subsequently referred to as conformer II). During this time frame, the probucol molecule in both solvent systems maintained its starting conformation. Figure 7C shows the S−C−S−C dihedral angle distribution over the last 1 ns of the 10 ns MD simulations. The probucol molecule in ethanol solution adopted the stable conformation expected from the curvature of the torsional potential energy surface, which is the molecular conformation of the probucol molecule in the Form I crystal (referred to as conformer I). However, the molecule in methanol solution showed a bimodal dihedral angle distribution of both the conformer II and I at an approximate average ratio of 2:1. This ratio was found to be consistent in subsequent MD simulations, regardless of the initial conformation of the probucol molecule in solution. Figure S3 provides a graphical representation of the dihedral angle distribution frequency at time points along a 10 ns MD simulation. The conformational preference of the probucol molecule between solvents is likely because of the combination of the solvent molecule size, the amphiphilic nature of both methanol and ethanol, and the solvent-accessible surface area (SASA) of each probucol conformer. The SASA of both conformers extracted from the DFT-optimized crystal structures was measured using the double cubic lattice method with the solvent probe radii of 1.84 and 2.14 Å, representing the hardsphere radii for methanol and ethanol, respectively. The measured SASA was 89.52 Å2 for conformer I and 94.66 Å2 for conformer II, using the methanol probe, and 95.60 and 101.20

Figure 7. Parameterization of the S−C−S−C torsion potential and dihedral angle distributions of the probucol molecules in solution from MD simulations. (A) Parameterization of the S−C−S−C torsion potential for MD simulations. (B) Calculated S−C−S−C dihedral angle distribution over MD simulation times of 0−0.1 ns for a probucol molecule solvated in pure methanol and pure ethanol. (C) Calculated S−C−S−C dihedral angle distribution over MD simulation times of 9−10 ns for a probucol molecule solvated in pure methanol and pure ethanol.

Å2, respectively, using the ethanol probe.70,71 Comparison of the SASA surfaces mapped with electrostatic potentials, computed at the PBE/6-31G(d,p) level, for the two conformers shows that conformer II has a considerably larger solvent accessibility in regions of high electrostatic potential (Figure 8).72 This conformation would lend itself to increased favorable interactions with smaller, higher polarity molecules, such as methanol. Electrostatic interactions between methanol and the probucol sulfur atoms could overcome the 8.2 kJ mol−1 difference between the energetically favored conformer I H

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

molecular-level view of complex interactions leading to polymorphism using a combination of computational methods and THz spectral data. Understanding the complexity of the crystallization process of pharmaceutical compounds and how specific solvent interactions yield various crystal forms will aid in the design and creation of effective polymorphs, solvatomorphs, and cocrystals with enhanced physical properties with respect to therapeutic performance.

Figure 8. Solvent-accessible and electrostatic potential surfaces for probucol conformers. SASA map of probucol generated using a 1.84 Å probe sphere for (A) conformer I and (B) conformer II. Blue surface indicates the surface inaccessible to the solvent, and red indicates high solvent accessibility. Electrostatic potential mapped on the SASA of probucol in (C) conformer I and (D) conformer II.

■

CONCLUSIONS

■

ASSOCIATED CONTENT

The conformational polymorphs of probucol were studied by a combinatorial approach using solid-state DFT, MD, and THz spectroscopy, and the solvent influence on the selective formation of the enantiotropically related Form I and Form II polymorphs is presented. The results of the solid-state DFT and MD calculations support that the preferred molecular orientations of the probucol molecule in solution and the possible inclusion of methanol in the Form II crystal lattice during the crystallization process provide the additional stabilization energy required to reach the metastable Form II crystal state over the thermodynamically favored Form I polymorph. This molecular-level study of pharmaceutical crystal formation demonstrates that the nucleation and growth of metastable polymorphs are aided significantly through subtle differences in solvent−molecule interactions and not a result of crystal packing forces. This is particularly noteworthy in that the probucol crystal structures are held only by weak dispersion interactions where solvent selection based on distinct electrostatic forces is not achievable. Solvent crystallization of pharmaceutical compounds for which clear intermolecular interactions are absent makes the question of polymorph selection quite difficult to address. As such, most screening procedures in search of pharmaceutical polymorphs rely on a trial-and-error approach. Gaining better understanding on the complexity of such crystallization processes will guide the discovery of new and more effective pharmaceutical formulations with enhanced therapeutic efficacy.

and conformer II in stabilizing probucol in the higher energy molecular conformation (Figure 7A). In addition, the smaller size of methanol, as compared to ethanol, allows for a greater contact between the extended conformer II because of the increased SASA. Evaluating short-range Coulomb and Lennard-Jones energies between probucol and the solvent from MD simulations showed conformer II to have an average Coulombic contribution of −65.1 (±3.5) kJ mol−1 in methanol compared to −48.6 (±3.8) kJ mol−1 in ethanol solution. This interaction energy decreased considerably more with conformer I in ethanol, which was determined to be −36.4 (±4.0) kJ mol−1. This supports that the increased stabilization of conformer II is achieved because of the greater accessibility of methanol to regions of higher electrostatic potential. In ethanol solution, the Lennard-Jones interactions between probucol and the solvent were dominant because of the significant relative increase in solvent hydrophobicity of ethanol as compared to methanol. Simulations in ethanol exhibited a −136.9 (±2.9) kJ mol−1 contribution from the probucol/solvent Lennard-Jones interactions, whereas this value was only −71.0 (±2.4) kJ mol−1 in methanol solution. The findings from both solid-state DFT and MD calculations suggest that this preferred molecular orientation of the probucol molecule in methanol solution, combined with the likely nonperiodic methanol inclusion stabilizing the crystal lattice during the crystallization process, allows the Form II polymorph to materialize under conditions that favor this metastable crystal form over the thermodynamically stable Form I. The preferred molecular conformation of probucol in methanol supports energetically advantageous nucleation and growth of the Form II system, even though this system consists of both a higher energy molecular and crystalline configuration. Initial formation of the metastable Form II is aided by additional stabilizing interactions as a disordered methanol solvate. The collective set of weak dispersion interactions holding the probucol and methanol molecules in the Form II· MeOH crystal is not networked in such a manner to provide stability to the solvated system; thus, the methanol molecules are able to diffuse from the crystal system through the significant void spaces, yielding the experimentally observed pure probucol Form II polymorph. The results of this study shed light on the crystallization process of probucol and the stabilization of higher energy polymorphs and molecular conformations by the crystallization solvent. We have demonstrated the influence of the solvent on the formation of pharmaceutical polymorphs and provide a

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.9b00792. Calculated bond lengths, angles, and dihedrals from solid-state DFT optimizations; MD analysis of methanol dynamics in the Form II·MeOH crystal; dihedral angle frequency distribution for the probucol molecule in methanol; and MD topologies in Gromacs format for probucol, methanol, and ethanol (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: (208) 4261019. Fax: (208) 426-1311. ORCID

Matthew D. King: 0000-0003-1011-3108 Notes

The authors declare no competing financial interest. I

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

■

synaptic impairments induced by amyloid β peptide in mice. Exp. Neurol. 2012, 233, 767−775. (18) Favari, E.; Zanotti, I.; Zimetti, F.; Ronda, N.; Bernini, F.; Rothblat, G. H. Probucol Inhibits Abca1-Mediated Cellular Lipid Efflux. Arterioscler., Thromb., Vasc. Biol. 2004, 24, 2345−2350. (19) Kawakami, K.; Ohba, C. Crystallization of Probucol from Solution and the Glassy State. Int. J. Pharm. 2017, 517, 322−328. (20) Burger, A.; Ramberger, R. On the Polymorphism of Pharmaceuticals and Other Molecular Crystals. Ii. Microchim. Acta 1979, 72, 273−316. (21) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; De La Pierre, M.; D’Arco, P.; et al. CRYSTAL14: A program for theab initioinvestigation of crystalline solids. Int. J. Quantum Chem. 2014, 114, 1287−1317. (22) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J., et al. Crystal14 User’s Manual; University of Torino: Torino, 2014. (23) Groom, C. R.; Bruno, I. J.; Lightfoot, M. P.; Ward, S. C. The Cambridge Structural Database. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2016, 72, 171−179. (24) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (25) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (26) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (27) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (28) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. Xii. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257−2261. (29) Hariharan, P. C.; Pople, J. A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213−222. (30) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. (31) Gilat, G. Analysis of Methods for Calculating Spectral Properties in Solids. J. Comput. Phys. 1972, 10, 432−465. (32) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (33) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (34) Karton, A.; Gruzman, D.; Martin, J. M. L. Benchmark Thermochemistry of the CnH2n+2 Alkane Isomers (n = 2−8) and Performance of Dft and Composite Ab Initio Methods for DispersionDriven Isomeric Equilibria. J. Phys. Chem. A 2009, 113, 8434−8447. (35) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (36) Pascale, F.; Zicovich-Wilson, C. M.; López Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. The Calculation of the Vibrational Frequencies of Crystalline Compounds and Its Implementation in the Crystal Code. J. Comput. Chem. 2004, 25, 888−897. (37) Ferrero, M.; Rérat, M.; Kirtman, B.; Dovesi, R. Calculation of First and Second Static Hyperpolarizabilities of One- to ThreeDimensional Periodic Compounds. Implementation in the Crystal Code. J. Chem. Phys. 2008, 129, 244110. (38) Ferrero, M.; Rérat, M.; Orlando, R.; Dovesi, R. Coupled Perturbed Hartree-Fock for Periodic Systems: The Role of Symmetry

ACKNOWLEDGMENTS The authors acknowledge the high-performance computing support of the R2 computed cluster (DOI: 10.18122/ B2S41H) provided by Boise State University’s Research Computing Department. This material is also based upon the work supported by the National Science Foundation under grant no. 1229709. The authors thank Boise State University for continued support.

■

REFERENCES

(1) Brittain, H. G. Polymorphism and Solvatomorphism 2010. J. Pharm. Sci. 2012, 101, 464−484. (2) Karpinski, P. H. Polymorphism of Active Pharmaceutical Ingredients. Chem. Eng. Technol. 2006, 29, 233−237. (3) Yu, L. X. Pharmaceutical Quality by Design: Product and Process Development, Understanding, and Control. Pharm. Res. 2008, 25, 781−791. (4) King, M. D.; Buchanan, W. D.; Korter, T. M. Understanding the Terahertz Spectra of Crystalline Pharmaceuticals: Terahertz Spectroscopy and Solid-State Density Functional Theory Study of (S)(+)-Ibuprofen and (Rs)-Ibuprofen. J. Pharm. Sci. 2011, 100, 1116− 1129. (5) King, M. D.; Buchanan, W. D.; Korter, T. M. Identification and Quantification of Polymorphism in the Pharmaceutical Compound Diclofenac Acid by Terahertz Spectroscopy and Solid-State Density Functional Theory. Anal. Chem. 2011, 83, 3786−3792. (6) Delaney, S. P.; Pan, D.; Yin, S. X.; Smith, T. M.; Korter, T. M. Evaluating the Roles of Conformational Strain and Cohesive Binding in Crystalline Polymorphs of Aripiprazole. Cryst. Growth Des. 2013, 13, 2943−2952. (7) Saito, S.; Inerbaev, T. M.; Mizuseki, H.; Igarashi, N.; Kawazoe, Y. Terahertz Vibrational Modes of Crystalline Salicylic Acid by Numerical Model Using Periodic Density Functional Theory. Jpn. J. Appl. Phys. 2006, 45, 4170−4175. (8) Delaney, S. P.; Pan, D.; Galella, M.; Yin, S. X.; Korter, T. M. Understanding the Origins of Conformational Disorder in the Crystalline Polymorphs of Irbesartan. Cryst. Growth Des. 2012, 12, 5017−5024. (9) Gupta, J.; Nunes, C.; Vyas, S.; Jonnalagadda, S. Prediction of Solubility Parameters and Miscibility of Pharmaceutical Compounds by Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 2014−2023. (10) Xiang, T.-X.; Anderson, B. D. Molecular Dynamics Simulation of Amorphous Indomethacin. Mol. Pharm. 2013, 10, 102−114. (11) Hamad, S.; Moon, C.; Catlow, C. R. A.; Hulme, A. T.; Price, S. L. Kinetic Insights into the Role of the Solvent in the Polymorphism of 5-Fluorouracil from Molecular Dynamics Simulations. J. Phys. Chem. B 2006, 110, 3323−3329. (12) Gerges, J.; Affouard, F. Predictive Calculation of the Crystallization Tendency of Model Pharmaceuticals in the Supercooled State from Molecular Dynamics Simulations. J. Phys. Chem. B 2015, 119, 10768−10783. (13) Gerber, J. J.; Caira, M. R.; Lötter, A. P. Structures of Two Conformational Polymorphs of the Cholesterol-Lowering Drug Probucol. J. Crystallogr. Spectrosc. Res. 1993, 23, 863−869. (14) Kim, J. H.; Hong, K. W.; Bae, S. S.; Shin, Y.-I.; Choi, B. T.; Shin, H. K. Probucol plus cilostazol attenuate hypercholesterolemiainduced exacerbation in ischemic brain injury via anti-inflammatory effects. Int. J. Mol. Med. 2014, 34, 687−694. (15) Yamashita, S.; Matsuzawa, Y. Where Are We with Probucol: A New Life for an Old Drug? Atherosclerosis 2009, 207, 16−23. (16) Al-Majed, A. A. Probucol Attenuates Oxidative Stress, Energy Starvation, and Nitric Acid Production Following Transient Forebrain Ischemia in the Rat Hippocampus. Oxid. Med. Cell. Longevity 2011, 2011, 471590. (17) Santos, D. B.; Peres, K. C.; Ribeiro, R. P.; Colle, D.; Santos, A. A. d.; Moreira, E. L. G.; Souza, D. O. G.; Figueiredo, C. P.; Farina, M. Probucol, a lipid-lowering drug, prevents cognitive and hippocampal J

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A and Related Computational Aspects. J. Chem. Phys. 2008, 128, 014110. (39) Ferrero, M.; Rérat, M.; Orlando, R.; Dovesi, R. The Calculation of Static Polarizabilities of 1-3D Periodic Compounds. The Implementation in the Crystal Code. J. Comput. Chem. 2008, 29, 1450−1459. (40) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1-2, 19−25. (41) Abraham, M. J.; Van Der Spoel, D.; Lindahl, E.; Hess, B. Gromacs User Manual Version 2018. 2018. (42) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. Gromacs: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (43) Jorgensen, W. L.; Tirado-Rives, J. The OPLS [Optimized Potentials for Liquid Simulations] Potential Functions for Proteins, Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (44) Simões, R. G.; Bernardes, C. E. S.; Joseph, A.; Piedade, M. F.; Emmerling, W.; Diogo, H. P.; Minas da Piedade, M. E.; Minas da Piedade, M. E. Polymorphism in Simvastatin: Twinning, Disorder, and Enantiotropic Phase Transitions. Mol. Pharm. 2018, 15, 5349− 5360. (45) Nemkevich, A.; Bürgi, H.-B.; Spackman, M. A.; Corry, B. Molecular Dynamics Simulations of Structure and Dynamics of Organic Molecular Crystals. Phys. Chem. Chem. Phys. 2010, 12, 14916−14929. (46) Larsen, A. S.; Ruggiero, M. T.; Johansson, K. E.; Zeitler, J. A.; Rantanen, J. Tracking Dehydration Mechanisms in Crystalline Hydrates with Molecular Dynamics Simulations. Cryst. Growth Des. 2017, 17, 5017−5022. (47) Gao, Y.; Olsen, K. W. Molecular Dynamics of Drug Crystal Dissolution: Simulation of Acetaminophen Form I in Water. Mol. Pharm. 2013, 10, 905−917. (48) Agrawal, P. M.; Rice, B. M.; Zheng, L.; Velardez, G. F.; Thompson, D. L. Molecular Dynamics Simulations of the Melting of 1,3,3-Trinitroazetidine. J. Phys. Chem. B 2006, 110, 5721−5726. (49) Dodda, L. S.; Cabeza de Vaca, I.; Tirado-Rives, J.; Jorgensen, W. L. Ligpargen Web Server: An Automatic OPLS-AA Parameter Generator for Organic Ligands. Nucleic Acids Res. 2017, 45, W331− W336. (50) Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. I. J. Chem. Phys. 1955, 23, 1833−1840. (51) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (52) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (53) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (54) Verlet, L. Computer ″Experiments″ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98−103. (55) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N· log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (56) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (57) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22, 245−268. (58) Rice, A.; Jin, Y.; Ma, X. F.; Zhang, X. C.; Bliss, D.; Larkin, J.; Alexander, M. Terahertz optical rectification from ⟨110⟩ zinc-blende crystals. Appl. Phys. Lett. 1994, 64, 1324−1326. (59) Wu, Q.; Litz, M.; Zhang, X. C. Broadband detection capability of ZnTe electro-optic field detectors. Appl. Phys. Lett. 1996, 68, 2924−2926.

(60) Alvarez-Ros, M.; Palafox, M. Conformational Analysis, Molecular Structure and Solid State Simulation of the Antiviral Drug Acyclovir (Zovirax) Using Density Functional Theory Methods. Pharmaceuticals 2014, 7, 695−722. (61) King, M. D.; Buchanan, W. D.; Korter, T. M. Application of London-Type Dispersion Corrections to the Solid-State Density Functional Theory Simulation of the Terahertz Spectra of Crystalline Pharmaceuticals. Phys. Chem. Chem. Phys. 2011, 13, 4250−4259. (62) Hakey, P. M.; Allis, D. G.; Hudson, M. R.; Ouellette, W.; Korter, T. M. Terahertz Spectroscopic Investigation ofS-(+)-Ketamine Hydrochloride and Vibrational Assignment by Density Functional Theory. J. Phys. Chem. A 2010, 114, 4364−4374. (63) King, M. D.; Hakey, P. M.; Korter, T. M. Discrimination of Chiral Solids: A Terahertz Spectroscopic Investigation ofl- anddlSerine. J. Phys. Chem. A 2010, 114, 2945−2953. (64) King, M. D.; Buchanan, W. D.; Korter, T. M. Investigating the Anharmonicity of Lattice Vibrations in Water-Containing Molecular Crystals through the Terahertz Spectroscopy ofl-Serine Monohydrate. J. Phys. Chem. A 2010, 114, 9570−9578. (65) Quertinmont, J.; Carletta, A.; Tumanov, N. A.; Leyssens, T.; Wouters, J.; Champagne, B. Assessing Density Functional Theory Approaches for Predicting the Structure and Relative Energy of Salicylideneaniline Molecular Switches in the Solid State. J. Phys. Chem. C 2017, 121, 6898−6908. (66) Ruggiero, M. T.; Zeitler, J. A.; Erba, A. Intermolecular Anharmonicity in Molecular Crystals: Interplay Between Experimental Low-Frequency Dynamics and Quantum Quasi-Harmonic Simulations of Solid Purine. Chem. Commun. 2017, 53, 3781−3784. (67) Ruggiero, M. T.; Zeitler, J. A. Resolving the Origins of Crystalline Anharmonicity Using Terahertz Time-Domain Spectroscopy and Ab Initio Simulations. J. Phys. Chem. B 2016, 120, 11733− 11739. (68) Zhang, F.; Wang, H.-W.; Tominaga, K.; Hayashi, M. Characteristics of Low-Frequency Molecular Phonon Modes Studied by Thz Spectroscopy and Solid-State Ab Initio Theory: Polymorphs I and III of Diflunisal. J. Phys. Chem. B 2016, 120, 1698−1710. (69) Turner, M. J.; McKinnon, J. J.; Wolff, S. K.; Grimwood, D. J.; Spackman, P. R.; Jayatilaka, D.; Spackman, M. A. Crystalexplorer17; University of Western Australia, 2017. (70) Eisenhaber, F.; Lijnzaad, P.; Argos, P.; Sander, C.; Scharf, M. The Double Cubic Lattice Method: Efficient Approaches to Numerical Integration of Surface Area and Volume and to Dot Surface Contouring of Molecular Assemblies. J. Comput. Chem. 1995, 16, 273−284. (71) Gogonea, V.; Bãleanu-Gogonea, C.; Osawa, E. Solvent Hard Sphere Diameter from van der Waals Volume a Statistical Analysis of Computed and Solubility Determined Solvent Diameters. J. Mol. Struct.: THEOCHEM 1998, 432, 177−189. (72) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera?A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605−1612.

K

DOI: 10.1021/acs.jpca.9b00792 J. Phys. Chem. A XXXX, XXX, XXX−XXX