ARTICLE pubs.acs.org/JPCB
Prediction of Solubility Parameters and Miscibility of Pharmaceutical Compounds by Molecular Dynamics Simulations Jasmine Gupta,*,†,‡ Cletus Nunes,‡ Shyam Vyas,# § and Sriramakamal Jonnalagadda*,† ,
†
Department of Pharmaceutics, Philadelphia College of Pharmacy, University of the Sciences in Philadelphia, Philadelphia, Pennsylvania 19104, United States ‡ Biopharmaceutics R&D, Bristol-Myers Squibb Co., New Brunswick, New Jersey 08903, United States § Accelrys, Inc., San Diego, California 92121, United States ABSTRACT: The objectives of this study were (i) to develop a computational model based on molecular dynamics technique to predict the miscibility of indomethacin in carriers (polyethylene oxide, glucose, and sucrose) and (ii) to experimentally verify the in silico predictions by characterizing the drug-carrier mixtures using thermoanalytical techniques. Molecular dynamics (MD) simulations were performed using the COMPASS force field, and the cohesive energy density and the solubility parameters were determined for the model compounds. The magnitude of difference in the solubility parameters of drug and carrier is indicative of their miscibility. The MD simulations predicted indomethacin to be miscible with polyethylene oxide and to be borderline miscible with sucrose and immiscible with glucose. The solubility parameter values obtained using the MD simulations values were in reasonable agreement with those calculated using group contribution methods. Differential scanning calorimetry showed melting point depression of polyethylene oxide with increasing levels of indomethacin accompanied by peak broadening, confirming miscibility. In contrast, thermal analysis of blends of indomethacin with sucrose and glucose verified general immiscibility. The findings demonstrate that molecular modeling is a powerful technique for determining the solubility parameters and predicting miscibility of pharmaceutical compounds.
1. INTRODUCTION The development of amorphous drugs is becoming increasingly important in the pharmaceutical industry. Amorphous materials have a higher free energy than their crystalline counterparts and, therefore, often exhibit higher apparent solubility and faster dissolution rates.1,2 Formulations containing a poorly soluble active (BCS Class II/IV; Biopharmaceutics Classification Systems based on aqueous solubility and permeability) in the amorphous state can thus provide increased bioavailability. The key challenge is to prevent the crystallization of the amorphous active during the shelf life of the drug product. A widely used approach in the pharmaceutical industry is to formulate an “amorphous dispersion” wherein the drug is dispersed in a carrier matrix (polymeric or nonpolymeric).3-6 Various methods, such as spray-drying, precipitation, and hot-melt extrusion, can be used to prepare the amorphous dispersions. These can then be readily converted to the desired solid dosage forms, such as capsules and tablets.7 This approach offers several advantages over solubilized liquid or semisolid formulations in terms of a simplified manufacturing process, a potential for high drug loading, and feasibility of fixed-dose combination products.7 Factors such as the glass transition temperature (Tg) of the active pharmaceutical ingredient (API), the Tg/Tm ratio for the API, log P, and Tg of the API-carrier composite play an important role in the manufacturability and stability of amorphous systems.8-11 The key consideration for developing a physically stable dispersion is the miscibility of the drug in the carrier matrix.12,13 The cohesive energy and the solubility parameter are important r 2011 American Chemical Society
descriptors in this context. The cohesive energy represents the total attractive forces within a condensed state resulting from intermolecular interactions and consists of electrostatic interactions, van der Waals interactions, and hydrogen bonds interactions.14 It is equivalent to the amount of energy required to separate the constituent atoms/molecules to an infinite distance where it approaches zero potential energy. Knowledge of cohesive energy is vital toward understanding the physical properties of a material, such as solubility, melting point, and predicting the interactions between two materials (e.g., drug and carrier). Cohesive energy density (CED) refers to the energy required to vaporize a mole of liquid per unit volume and is mathematically defined by eq 1,14,15 ΔEv 1=2 ΔHv - RT 1=2 1=2 δ ¼ ðCEDÞ ¼ ¼ ð1Þ Vm Vm where ΔEv = internal energy change of vaporization, ΔHv = enthalpy of vaporization, Vm = molar volume of the liquid at the temperature of vaporization, R = gas constant, and T = absolute temperature. The solubility parameter (δ), as defined by Hildebrand and Scott,15 is the square root of the cohesive energy density. The concept was initially developed for liquids on the basis of regular Received: September 7, 2010 Revised: December 8, 2010 Published: February 09, 2011 2014
dx.doi.org/10.1021/jp108540n | J. Phys. Chem. B 2011, 115, 2014–2023
The Journal of Physical Chemistry B
ARTICLE
solution theory and has been extended to solids with appropriate assumptions.14 According to Hildebrand and Scott, the enthalpy change of mixing ΔHm is expressed as eq 2,14,15 ΔHm ¼ V 0 ðδ1 - δ2 Þ2 φ1 φ2
ð2Þ
where V0 is the volume of mixture; φ is the volume fraction of components 1 and 2. ΔGm ¼ ΔHm - TΔSm
ð3Þ
The Gibbs free energy of mixing at constant pressure (ΔGm) governs the process of dissolution between two constituents and is related to the changes in the enthalpy (ΔHm) and the entropy (ΔSm) of the system (eq 3). For mixing to occur between two components, ΔGm must be negative. The ΔSm term is usually positive, and hence, ΔHm is the determining factor.16 As per eq 2, ΔHm becomes zero when δ values of the two components are equal, in which case mutual solubility is achieved due to a negative entropy term. In other words, mixtures of materials with similar solubility parameter values are thermodynamically compatible. Experimentally, the solubility parameter can be estimated by direct and indirect methods, as reported in the literature. The direct methods include measuring the heat of vaporization by calorimetry or measuring the solubility/miscibility of compounds of interest in solvents of a known solubility parameter.17,18 The former method cannot be used for pharmaceutical compounds with thermal instabilities and low vapor pressures. The latter method, though widely used in the field of polymer science, is very subjective and time-consuming. An indirect method utilizing inverse gas chromatography has been used for determining solubility parameters.19-21 This method is based on retention time of gases of known δ values; however, the technique suffers from the limitation of being a nonequilibrium method and is sensitive to surface heterogeneities on the stationary phase. Even though solubility parameter is an elegant tool to predict the compatibility of materials, it is very difficult to experimentally determine these parameters. Group contribution methods for estimating the solubility parameter are based on the knowledge of structural fragments within the molecule and are by far the most popular way of calculating δ, despite the known shortcomings.14,22 Such theoretical models are typically applicable for simple molecular structures, wherein van der Waals forces predominate. The method has limitations when dealing with complex molecules containing highly directional interactions (e.g., hydrogen bonding) or longrange interactions (e.g., electrostatic).22 Recent advances in the computational field are enabling mechanistic understanding of complex molecules and polymeric systems via techniques such as molecular mechanics and dynamics.23-26 For example, molecular modeling allows calculation of the time-dependent behavior of a molecular system with structural and conformational changes such as chirality and tacticity.27 Choi et al.28 and Kavassalis et al.29 have demonstrated the use of molecular modeling in estimating the solubility parameters of nonionic surfactants. In a recent study published by Simperler et al.,30 molecular dynamics simulations were used to determine the glass transition temperatures of carbohydrates. Our first objective in the current study was to develop a molecular dynamics technique to understand the miscibility of a model drug in three carrier matrices. The second goal was to corroborate the computational predictions using experimental tools; namely, differential scanning calorimetry and hot stage microscopy. Indomethacin (IND) was used as the model drug;
Figure 1. Molecular structures of (a) indomethacin, (b) polyethylene oxide, (c) glucose, and (d) sucrose.
polyethylene oxide (PEO), glucose (GLU), and sucrose (SUC) were selected as the polymeric and nonpolymeric carriers, respectively. Using full atomistic molecular dynamics (MD) simulations, we have constructed the atomic models and performed simulations under thermodynamic constraints to calculate the properties of interest, thus allowing miscibility prediction. We have also examined the nature and strength of various intermolecular interactions contributing to the miscibility.
2. MATERIALS AND METHODS 2.1. Materials. Crystalline indomethacin, sucrose, and glucose were obtained from Sigma Aldrich (St. Louis, Missouri, USA). Polyethylene oxide (Polyox WSR-N80, MW 200 000) was obtained from Dow Chemical Company (Midland, WI, USA). The materials were passed through a 60-mesh (250 μm) screen for research experiments. 2.2. Methods. 2.2.1. Molecular Dynamics Simulations. The MD simulations were performed using the Materials Studio 4.0 (Accelrys, San Diego, CA). The molecular structures of indomethacin, sucrose, and glucose were obtained from the Cambridge Structural Database (Figure 1). In the case of PEO, a homopolymer of PEO containing 30 repeat units of ethylene oxide was constructed using the Materials Visualizer module within Materials Studio. The ab initio COMPASS (condensedphase optimized molecular potentials for atomistic simulation studies) force field was utilized for MD simulations.31 Ten in silico amorphous models (“amorphous cells”) each for indomethacin (60 molecules), sucrose (72 molecules), glucose (144 molecules), and PEO (MW 18509, 14 chains each containing 30 monomer units of ethylene oxide), were generated at 298 K. The dimensions of the amorphous cells were 3 nm 3 nm 3 nm, with periodic boundary conditions. The lowest energy configuration of the 10 amorphous cells for each compound was selected for energy minimization. The energies of the amorphous cells were minimized using the steepest decent and FletcherReeves conjugate gradient procedures to remove unfavorable interactions and attain the lowest energy state. 2015
dx.doi.org/10.1021/jp108540n |J. Phys. Chem. B 2011, 115, 2014–2023
The Journal of Physical Chemistry B Following minimization, the MD simulations were carried out in two phases: the equilibration phase and the production phase. In the equilibration phase, the amorphous cells were allowed to relax for 2 ns under isothermal (NVT) or isobaric-isothermal conditions (NPT-NVT) at 298 K. In the production phase, the equilibrated structure was processed via the NVT ensemble for 200 ps at 298 K with a time step of 1 fs. For PEO, the MD simulation time was increased to 500 ps to enable equilibration of the polymeric chains for a longer time. The Anderson thermostat and barostat was used to maintain the temperature and pressure; respectively. The nonbonded van der Waals and electrostatic interactions were truncated using group based cutoff distance of 1.25 nm. Trajectory frames were captured during the production run and, the data from the final 50 ps was used for computing the cohesive energy density and solubility parameter for the compound. The procedure was repeated to generate three independent data sets for predicting the average CED and δ. 2.2.2. Theoretical Solubility Parameters. The Hansen solubility parameters of the compounds were calculated from the chemical structures using Hoftyzer/Van Krevelen and Hoy methods.32,33 The details are presented in the Results and Discussion section. For polyethylene oxide, determination of the solubility parameter was based on the average molecular weight of the polymer. 2.2.3. Differential Scanning Calorimetry (DSC). A differential scanning calorimeter (Q1000, TA Instruments, New Castle, DE, USA) equipped with a refrigerated cooling accessory was used. The instrument was calibrated with indium. Physical mixtures of indomethacin and excipients were prepared at four mass ratios (10:90, 25:75, 50:50, and 75:25) by geometrically mixing. Samples (∼8-10 mg) were packed in aluminum pans and crimped with lids having pinholes. The DSC sample chamber was purged with dry nitrogen (50 mL/min) during the experiments. The samples were then subjected to a heat-cool-heat cycle as follows to determine the melt onset, Tg, and the miscibility of indomethacin with the excipients. Cycle 1: Samples were initially heated from 25 °C to the target temperature (above or below the melting point) at a rate of 1 °C/min. The samples were then held isothermally at the target temperature for 5 min. Cycle 2: Samples were subsequently cooled at a rate of 20 °C/min to a temperature of -60 °C and held at -60 °C for 10 min. Cycle 3: Finally, the samples were heated again at a rate of 1 °C/min to the desired temperature. 2.2.4. Hot Stage Microscopy (HSM). Binary mixtures of indomethacin and excipients in various ratios (10:90, 50:50, 75:25) were observed by hot stage microscopy using an Olympus BX51 polarized microscope (Olympus America Inc., Center Valley, PA, USA) coupled with Instec’s HCS400 hot stage and STC200 Standalone Temperature Controller (Instec Research Instrumentation Technology, Boulder, CO, USA). The mixtures were heated at a rate of 5 °C to the melting point of the constituents. A polarized light source was utilized, and images were captured at regular intervals to observe miscibility.
3. RESULTS AND DISCUSSION 3.1. Calculation of the Solubility Parameter by Molecular Modeling. In molecular dynamics (MD) simulations, the mi-
crostates of a system are specified by the position and momentum of the constituent atoms and are obtained by solving the Newtonian equation of motion (F = ma) for a set of interacting
ARTICLE
Figure 2. Simulated amorphous cells of (a) IND with 60 molecules, (b) PEO (30-mer) with 14 chains, (c) GLU with 144 molecules, and (d) SUC with 72 molecules.
particles. The equation is transformed to a time independent Hamiltonian function which represents the total energy expressed as sum of the kinetic and potential energy of the system. The future positions and velocities of the particles are calculated thereof using a numerical algorithm based on the initial positions and velocities. In this work, MD simulations were used to compute the solubility parameters of the compounds of interest. The simulations were performed using the ab initio COMPASS force field. In this force field, the total potential energy Epot is defined as the sum of several terms arising from the parametrization of the bonding and nonbonding interactions and is expressed as follows31 Epot ¼ Eb þ Eθ þ Eφ þ Eχ þ Ecross þ ECou þ EvdW
ð4Þ
where Eb = bond stretching energy, Eθ = valence angle bending energy, EΦ = dihedral torsion energy, Eχ = out-of-plane energy, Ecross = cross-term interaction energy, ECou = Coulombic interaction energy, and EvdW = van der Waals interaction energy. The first five terms parametrize the short-range bonded interactions, whereas the last two terms parametrize the inter- and intramolecular nonbonded interactions. The electrostatic or Coulombic interactions stem from the attraction or repulsion between the two charged atoms (i, j) and is inversely proportional to the distance separating the two charges as defined in eq 5.34 qi qj Cðrij Þ ¼ ð5Þ 4πε0 rij Here, rij represents the distance between atoms i and j, q is the atom charge, and ε0 is the dielectric permittivity of the vacuum. The van der Waals interactions include attractive and repulsive forces arising from the interactions between molecules polarizing each other and from the van der Waals volume overlap, respectively. This is often defined using a Lennard-Jones potential function. 2016
dx.doi.org/10.1021/jp108540n |J. Phys. Chem. B 2011, 115, 2014–2023
The Journal of Physical Chemistry B
ARTICLE
Figure 3. Energy-minimized structure of indomethacin relative to its initial configuration at 298 K.
The total energy for an isolated molecule and for the bulk system with periodic boundary conditions is computed using eq 4. The difference in the energy between the two is an indicator of the energy required to vaporize the molecules and is related to CED via eq 1. Mathematically, it is expressed as ΔEvap ¼ EtotðbulkÞ - EtotðmoleculeÞ
ð6Þ
Bulk amorphous cell models of IND, SUC, GLU, and PEO were constructed (Figure 2) using periodic boundary conditions, as described in the Methods section. The constructed amorphous cells were subjected to a series of energy minimization to attain a state of minimal potential energy and to allow efficient packing of molecules in the cell. This was achieved by using the steepest descent initially followed by conjugate gradient methods.35 The energy-minimized structure of the unit cell of IND illustrating the efficiency in cell packing relative to its initial configuration is shown in Figure 3. The energy-minimized cells were used as a start point for the MD simulations (i.e., the equilibration stage). In the MD equilibration stage, NVT followed by the NPT ensemble was employed to allow each cell to attain the equilibrium density for each structure at 298 K. The equilibration time depends on the number of atoms in a specific system, and therefore, the simulations were carried out until the total energy of each system was stabilized. The time evolution energy profiles during MD simulations of indomethacin (nonpolymeric) and PEO (polymeric) are shown in Figure 4. The results show that the duration of 2 ns was adequate to attain the equilibration without any significant changes in the density of the structures. Further, to avoid entrapment of simulated systems in a metastable state of local high-energy minima, MD simulations were carried out in two stages (NVT followed by NPT) to provide thermal energies to cross-energy barriers between local minima. The resultant density values are reported in Table 1 and are in good agreement with the experimental values reported in the literature.23,36-38 The close agreement between the two density values signifies the effective optimization of intermolecular interactions achieved between the molecules and, thus, its packing in reality within the structures. In the MD production stage, each structure was subjected to the NVT ensemble for 200-500 ps to calculate the physical properties of interest (e.g., CED and δ). A snapshot of the final trajectory
Figure 4. Plot of the total potential energy, nonbonded energy, and temperature of (a) IND and (b) PEO as a function of MD simulation time.
Table 1. MD Simulated and Experimental Density (g/cm3) at 298 K compounds
experimental density
simulated density
IND
1.32
1.30 ( 0.03
PEO
1.11
1.12 ( 0.00
GLU
1.50
1.52 ( 0.00
SUC
1.43
1.44 ( 0.02
frame of amorphous cell of IND resulting from the NVT production run is shown in Figure 5. Furthermore, to determine minimum 2017
dx.doi.org/10.1021/jp108540n |J. Phys. Chem. B 2011, 115, 2014–2023
The Journal of Physical Chemistry B
ARTICLE
The solubility parameter as per the Hoy method is given by the following equation, ð9Þ δt ¼ ðFt þ BÞ = V
Figure 5. Final trajectory frame of IND cell resulting from the NVT production run at 298 K.
Table 2. Solubility Parameters from MD Simulations and Group Contribution Methods (MPa0.5) compounds
δVan Krevelen
δHoy
δMD simulations
IND
23.3
21.9
23.9 ( 0.3
PEO GLU
22.9 39.0
20.0 37.8
22.2 ( 0.2 34.8 ( 0.2
SUC
36.0
33.5
29.9 ( 0.5
chain length that represents the real polymer chain, the solubility parameter of PEO was examined as a function of chain length (20, 30, 40, 55 monomeric units). It is well established that the solubility parameter increases with increasing chain length and reaches a point where no further increase is observed.39 A stable value of δ indicates the number of repeating units is sufficient for simulations. For PEO, it was observed that the δ values remained stable beyond 30 monomer units, suggesting it to be a reasonable chain length for performing further simulations. The solubility parameter values for IND and excipients computed from the MD simulations (production stage) at 298 K are summarized in Table 2. The values were compared with those calculated using the Hoftyzer-Van Krevelen and Hoy group contribution methods as described below. As per the Hoftyzer-Van Krevelen method, P 2 1=2 P P 1=2 Fpi Fhi Ehi δh ¼ δp ¼ δh ¼ ð7Þ V V V and δt 2 ¼ δd 2 þ δp 2 þ δh 2
ð8Þ
where Fdi and Fpi are the molar attraction constants due to dispersion and polar components, Ehi is hydrogen bonding energy and V is the molar volume. δt is the total solubility parameter (also referred to as the Hansen solubility parameter). δd, δp, and δh represent the dispersive, electrostatic-polar, and hydrogen bonding forces, respectively.
where Ft = molar attraction constant, B = base constant, and V = molar volume. The molar volume of the compounds was calculated on the basis of its density and molecular weight. The values of Fdi, Fpi, Ehi, Ft, and B for various groups were calculated from the literature.32,33,40 As seen in Table 2, there was a reasonable agreement between the solubility parameter values computed from MD simulations and calculated using the group contribution methods for each compound. This was particularly true in the case of IND and PEO, for which the two methods yielded comparable values. In contrast, the solubility parameter values for GLU and SUC obtained from the MD simulations were somewhat lower than those calculated from the group contribution methods. The difference can be attributed to the predominance of the hydrogen bonding in GLU and SUC, which cannot be adequately described using the group contribution methods. The H-bond network was visualized using the final MD trajectory frames for amorphous indomethacin, glucose, and sucrose. As evident from Figure 6, sucrose and glucose show a significantly greater extent of H-bond network as compared with indomethacin. The advantage of MD is that it takes into account such directional dependent interactions during computations. In contrast, the group contribution methods generally do not account for directional interactions and, hence, may overestimate the CED.28 The study demonstrates that the MD simulations can be used as an effective tool for screening and determining solubility parameters of pharmaceutical compounds. Moreover, the computational models can afford a molecular level insight into the nature of interactions within complex systems. Another aspect is that the in silico analysis can be performed at various temperatures pertinent to amorphous dispersions (e.g., hot-melt extrusion). This can potentially overcome the limitations of the experimental methods that are typically restricted to room temperature.28 Langer et al.22 have also reported that the calculation of solubility parameters solely via group contribution approaches can yield different values because different contributions for one structural unit can be obtained using group contributions either for the complete functional group or for the sum of its single components. This discrepancy is pronounced for systems containing tertiary, secondary, and primary amide groups or aromatic systems. 3.2. In Silico Prediction of Miscibility of Drug Carrier. It has been well established in the literature that compounds with similar values of δ are thermodynamically miscible. This is because the nonbonded interaction energies are balanced during mixing. A classification system based on the solubility parameter differences (Δδ) in the molten state has been cited in literature for predicting the miscibility of drug and excipients. For instance, Greenhalgh et al.41 have demonstrated that compounds with Δδ < 7.0 MPa0.5 are likely to be miscible, whereas Δδ > 10.0 MPa0.5 is indicative of an immiscible system. Forster et al.13 have shown that compounds with Δδ < 2.0 MPa0.5 are miscible, as exemplified by the lacidipine-PEG 8 system. The solubility parameter differences, as computed from MD simulations, between IND and the carriers have been summarized in Table 3. It is evident from the data that IND/PEO is likely to be miscible owing to a small Δδ (
