3014
Biomacromolecules 2008, 9, 3014–3023
Application of Molecular Dynamics Simulation To Predict the Compatability between Water-Insoluble Drugs and Self-Associating Poly(ethylene oxide)-b-poly(ε-caprolactone) Block Copolymers Sarthak Patel,† Afsaneh Lavasanifar,‡ and Phillip Choi*,† Department of Chemical and Materials Engineering and Faculty of Pharmacy and Pharmaceutical Sciences, University of Alberta, Edmonton, Alberta, Canada Received March 28, 2008; Revised Manuscript Received September 12, 2008
In the present work, molecular dynamics (MD) simulation was applied to study the solubility of two waterinsoluble drugs, fenofibrate and nimodipine, in a series of micelle-forming PEO-b-PCL block copolymers with combinations of blocks having different molecular weights. The solubility predictions based on the MD results were then compared with those obtained from solubility experiments and by the commonly used group contribution method (GCM). The results showed that Flory-Huggins interaction parameters computed by the MD simulations are consistent with the solubility data of the drug/PEO-b-PCL systems, whereas those calculated by the GCM significantly deviate from the experimental observation. We have also accounted for the possibility of drug solubilization in the PEO block of PEO-b-PCL.
Introduction In recent years, polymeric micelles formed by block copolymers have emerged as a novel nanoscopic vehicle for the delivery of water-insoluble drugs in a controlled manner.1-4 Because such block copolymers can be synthesized using a wide variety of chemical moieties, this allows us to develop tailormade carriers. The challenge is to identify the molecular structure of the blocks in the block copolymer that can encapsulate the drug of interest and at the same time provide the desired release properties. It is believed that compatibility between a drug and the hydrophobic block of a block copolymer determines the encapsulation capacity of the micelles.1,5,6 Nonetheless, the use of drug-compatible moieties in the micellar core has been shown to lower the rate of drug release from the carrier.1 To optimize the above performance properties, the current practice is to synthesize the target block copolymer and examine the solubility of drugs in it by using a trial and error approach. This is simply because such drug delivery systems involve complex intermolecular interactions and fittings of molecules of different shapes. Obviously, being able to quantify such interactions between various drugs and block copolymers with different chemical moieties will definitely facilitate the block-copolymer design process. This will in turn allow us to avoid cumbersome and expensive trial and error formulation studies. Molecular simulation seems to be the method of choice. Ideally, one would like to simulate block-copolymer micelles to determine how block-copolymer structure and block lengths affect drug loadings and encapsulation efficiency of the drug of interest. Nevertheless, simulating such systems at the atomistic level on a routine basis is currently totally impossible because computational costs are prohibitively high. Therefore, we have focused on studying binary interactions between the drugs and block copolymers of interest in their liquid state rather * To whom correspondence should be addressed. E-mail: phillip.choi@ ualberta.ca. Fax: (780) 492-2881. † Department of Chemical and Materials Engineering. ‡ Faculty of Pharmacy and Pharmaceutical Sciences.
than in a micelle environment. In particular, the aim of the present work is to test whether molecular dynamics (MD) simulation would be a better method to determine the compatibility of drug/block-copolymer systems than the group contribution method (GCM) that is commonly used in the pharmaceutical industry. As shown in this article, the MD approaches we used seem to outperform the GCM. The concept of the solubility parameter is used to quantify the strength of intermolecular interactions. The two most common solubility parameters used are the Hildebrand (δ) and Hansen (δd, δp, and δh) solubility parameters. The well-known definition of the solubility parameter (δ), as depicted in the following expression, relates the solubility parameter to the interaction energy potential, u(r), and the local arrangement of molecules in the liquid state, g(r) (i.e., the radial distribution function)
δ)
∆EV ) V
2πn2 V2
∞
∫ u(r)g(r)r2dr
(1)
0
where ∆EV/V is the cohesive energy density (CED) and n is the number of molecules. The Hildebrand and Hansen solubility parameters are related by the following equation
δ2 ) (δd2 + δp2 + δh2)
(2)
where the subscripts d, p, and h refer to the contributions due to dispersion forces, polar forces, and hydrogen bonding, respectively. In general, compatible materials tend to have comparable solubility parameters (either δ or δd, δp, and δh). A common method for obtaining Hildebrand and Hansen solubility parameters is the GCM. The GCM is based on the idea that the total intermolecular interactions among the molecules in a liquid is the linear sum of the contributions from various chemical moieties within the molecules. The following equations signify such an idea
10.1021/bm800320z CCC: $40.75 2008 American Chemical Society Published on Web 10/21/2008
Predicting Drug/PEO-b-PCL Compatibility
δd ) δp
∑ Fdi V
( ∑ Fpi2)1/2 )
δh )
V
(∑ Ehi) V
Biomacromolecules, Vol. 9, No. 11, 2008
(3)
Table 1. Reported Experimental Solubilities of Fenofibrate11 and Nimodipine12 in Caprolactone Along with the Solubility Parameters of Both Drugs and PCL Block Calculated by the GCM10 drug
(4)
1/2
(5)
where Fdi, Fpi, and Ehi refer to the contributions from various chemical moieties, and such values can be obtained from the Hoftyzer-Van Krevelen method.7 The molar volume (V) of the material can be determined by the Fedors method.8 The GCM is simple to use and generally provides reasonable compatibility predictions for compounds with simple chemical structures. Nevertheless, the GCM fails to provide a good estimation of the solubility of complex drug molecules and copolymers for the following reasons:7 (a) the method fails to take into consideration the geometry of the molecules involved, (b) it also fails to take into account the excluded volume interactions that are especially prevalent for long-chain copolymers, (c) it would yield the same values of solubility parameters for isomers because of its inability to distinguish between isomers that have identical chemical structures but different constitution and configuration, and (d) it tends to underestimate δp.9 Forrest et al.10 used the GCM to calculate the Hildebrand solubility parameters of a few lipophilic drugs and of the hydrophobic block (i.e., the PCL block) of a commonly used block copolymer poly(ethylene oxide)-b-poly(ε-caprolactone) (PEO-b-PCL, see Figure 1a). Using such δ, they calculated the Flory-Huggins interaction parameters (χ), which quantify the difference in the intermolecular interactions of the components in a binary mixture, and thereby the compatibility, of various pairs of drugs and the PCL block of PEO-b-PCL. According to the Flory-Huggins solution theory, lower χ values indicate better solubility. It is evident in their results that χ values between hydrophobic drugs such as fenofibrate and nimodipine (see Figure 1b,c) and the PCL block of PEO-b-PCL that was predicted using the GCM are inconsistent with their solubility in caprolactone obtained by experimental means. (See Table 1.) In particular, the Hildebrand solubility parameters of fenofibrate and nimodipine computed by the GCM are 22.5 and 21.6 (J/cm3)1/2, respectively, and the corresponding Flory-
3015
PCL fenofibrate nimodipine
δ ((J/cm3)1/2) χdrug/PCL drug/caprolactone (mmol/mol)a 20.20 22.50 21.60
0.539 0.250
120 12
a Reported solubility based on maximum incorporation reported by referenced authors in PEG-b-PCL or PEG-b-PCL-b-PEG nanocarriers.
Huggins interaction parameters (χ) between them and the PCL block are 0.539 and 0.250, respectively. However, the caprolactone solubility of fenofibrate and nimodipine was 120 and 20 mmol/mol, respectively. (See Table 1.) Higher solubility of fenofibrate in PEO-b-PCL was unexpected because their Flory-Huggins interaction parameter is greater than that of the nimodipine/PCL block of the PEO-b-PCL pair. It is worth noting that the experimental solubility data reported by Forrest et al. was obtained from Jette et al.11 and Ge et al.12 It signifies the amount of soluble drug in the micelles formed by the corresponding diblock and triblock block copolymers made up of PEO and PCL blocks for fenofibrate and nimodipine, respectively, but was converted to a basis of mmol drug/mol caprolactone. In the case of nimodipine, no drug-loading data are available for diblock PEO-b-PCL. Considering the inadequacy of the GCM, in the present work, we applied the technique of MD simulation to determine whether it would provide better compatibility predictions for the systems for which the GCM yielded incorrect results. The motivation originates from the work of Choi et al.9,13-15 in which the authors utilized MD simulation to obtain a better estimation of solubility parameters of surfactants as compared with the GCM and showed that such an approach can rectify some of the drawbacks of the GCM. Because one can account not only for the interactions between the drug and the PCL block but also for the interactions between the drug and the PEO block in the MD simulation, we carried out simulations using PEO-b-PCL with different hydrophilic (PEO) and hydrophobic (PCL) block lengths. Recently, an in silico model was proposed by Huynh et al.16 that involves calculations of the solubility parameters of anticancer agent docetaxel and various small molecule excipients by using semiempirical methods and MD simulation. The
Figure 1. Chemical structure of (a) PEO-b-PCL, (b) fenofibrate, and (c) nimodipine.
3016
Biomacromolecules, Vol. 9, No. 11, 2008
Patel et al.
relative in silico solubility of docetaxel in various excipients was in good agreement with the experimental solubility data and hence validated the use of the computational model as a reliable tool for designing drug-excipient mixture formulations. Nevertheless, our study concentrates on the prediction of the thermodynamic compatibility between drug and long-chain block copolymer. Two different MD strategies were used so that effects of drug loading and of hydrophilic and hydrophobic block lengths of PEO-b-PCL on the drug/PEO-b-PCL compatibility could be investigated.
Materials and Methods Software and Force Field. All MD simulations were performed using the Materials Studio software package (MS Modeling version 4.0, Accelrys) run on a Silicon Graphics (SGI) workstation cluster. The initial block-copolymer conformations were generated using the rotational isomeric state (RIS) theory. The COMPASS force field17 was used throughout to describe bonded and nonbonded interactions, as shown in the following equation
E ) Eb + Eθ + Eφ + Eχ + EvdW + EQ
(6)
The first four terms represent bonded interactions that correspond to the energy associated with bond stretching (Eb), bond angle bending (Eθ), torsion angle rotations (Eφ), and Wilson out-of-angle displacements (Eχ). The last two terms represent nonbonded interactions consisting of Lennard-Jones (LJ) 9-6 function for the van der Waals interactions and the Coulombic function for the electrostatic interaction. On the basis of first principle quantum mechanical calculations, the partial atomic charges on the molecules were preset by the COMPASS force field. The electrostatic interaction was calculated using the Ewald summation method because it provides a more effective way of handling long-range interactions.18 Construction of Liquid-State Models. The liquid-state models of the block copolymers, drugs, and their mixtures were built using the methodologies developed by Theodorou and Suter.19 Such structures were constructed using the amorphous builder module that is available in Cerius2. Three PEO-b-PCL block copolymers with different hydrophilic (PEO) and hydrophobic (PCL) block lengths were used. They were PEO(1250)b-PCL(2500), PEO(2500)-b-PCL(2500), and PEO(2500)-bPCL(1250), where the number in the parentheses signifies the molecular weight of the block. Snapshots of liquid-state models of PEO(2500)-b-PCL(2500) block copolymer and a mixture of PEO(2500)-b-PCL(2500) block copolymer and fenofibrate drugs (40% w/w drug/polymer) are shown in Figure 2a,b, respectively. In the case of block copolymers, single-chain conformations in unit cells subjected to 3D periodic boundary conditions were used. Here the periodic boundary conditions provide a means to model bulk liquid state by using only single chains. As singlechain conformations were grown, several constraints were imposed. We avoided hard overlaps by ascribing to each atom a hard-core radius equal to 0.3 times its van der Waals radius. The unit cell was also subjected to density constraint. Because experimental density values of the drugs and their mixtures with various block copolymers at 140 °C (i.e., 413 K) were not available, we carried out an isobaric-isothermal (NPT) ensemble MD simulation (P ) 1 atm, T ) 413 K) using the density values from the GCM as the initial values to determine their density. Here the Nose´ thermostat20 and Andersen barostat21 were used to control the temperature and pressure of the systems. In the NPT MD simulation, the volume of the periodic unit cell (i.e., density) was allowed to fluctuate.
Figure 2. Snapshots of the liquid model of (a) pure PEO(2500)-bPCL(2500) block copolymer and (b) a mixture of PEO(2500)-bPCL(2500) block copolymer and fenofibrate drug molecules (40% w/w drug/block copolymer) subjected to 3D periodic boundary conditions.
The final, and perhaps most important, constraint imposed was the distribution of the torsion angles of the skeletal bonds in the copolymers, which was determined using the RIS theory.22,23 Nevertheless, we note that the RIS theory cannot be used for polymers in their crystalline state. In particular, we determined the torsion angle distribution by applying the Boltzmann weighting factor to the energies of the RIS minima. In this work, the initial distribution of torsion angles was
Predicting Drug/PEO-b-PCL Compatibility Table 2. Unique Torsion Angles for the PEO Block
determined by the following method. We identified a total of seven torsion angles for PEO-b-PCL, and they are described in Tables 2 and 3. We constructed a conformational energy map for each of these torsion angles by rotating the torsion angle through 360° in 10° intervals while simultaneously relaxing all other degrees of freedom to achieve a root-mean-square force of
