In silico modelling of drug–polymer interactions for pharmaceutical formulations

Samina Ahmad, Blair F. Johnston, Simon P. Mackay, Andreas G. Schatzlein, Paul Gellert, Durba Sengupta, Ijeoma F. Uchegbu


Selecting polymers for drug encapsulation in pharmaceutical formulations is usually made after extensive trial and error experiments. To speed up excipient choice procedures, we have explored coarse-grained computer simulations (dissipative particle dynamics (DPD) and coarse-grained molecular dynamics using the MARTINI force field) of polymer–drug interactions to study the encapsulation of prednisolone (log p = 1.6), paracetamol (log p = 0.3) and isoniazid (log p = −1.1) in poly(l-lactic acid) (PLA) controlled release microspheres, as well as the encapsulation of propofol (log p = 4.1) in bioavailability enhancing quaternary ammonium palmitoyl glycol chitosan (GCPQ) micelles. Simulations have been compared with experimental data. DPD simulations, in good correlation with experimental data, correctly revealed that hydrophobic drugs (prednisolone and paracetamol) could be encapsulated within PLA microspheres and predicted the experimentally observed paracetamol encapsulation levels (5–8% of the initial drug level) in 50 mg ml−1 PLA microspheres, but only when initial paracetamol levels exceeded 5 mg ml−1. However, the mesoscale technique was unable to model the hydrophilic drug (isoniazid) encapsulation (4–9% of the initial drug level) which was observed in experiments. Molecular dynamics simulations using the MARTINI force field indicated that the self-assembly of GCPQ is rapid, with propofol residing at the interface between micellar hydrophobic and hydrophilic groups, and that there is a heterogeneous distribution of propofol within the GCPQ micelle population. GCPQ–propofol experiments also revealed a population of relatively empty and drug-filled GCPQ particles.

1. Introduction

Trial and error experiments dominate selection procedures for drug encapsulating polymeric excipients used in pharmaceutical formulations, with no systematic method available to select an appropriate functional polymer. In this report, we explore the computer simulation of polymer–drug interactions as a possible method for selecting polymers for drug encapsulation. Polymers with their large dimensions pose a considerable challenge for atomistic computer simulations in terms of time and computational power. However, coarse-grained computer simulations, which group a number of atoms or molecules together into single particles (Maiti & McGrother 2004), may be useful in this regard as coarse graining leads to fewer interacting species in a simulation, allowing larger length and time scales to be employed (Frenkel & Smit 2002). Such coarse-grained simulations, which lie between the atomistic and macroscopic scales, lead to the construction of models of 10–1000 nm in size (Groot & Warren 1997; Maiti & McGrother 2004; McGrother et al. 2004). Coarse-grained modelling is useful for studying lipid (Marrink et al. 2009) and block copolymer (Daoulas & Muller 2010) membrane organization, protein–lipid membrane interactions (Venturoli et al. 2005) and polymer properties (Glotzer & Paul 2002). Two coarse-grained computational methods have been employed in this work: a mesoscale approach—dissipative particle dynamics (DPD; Groot & Warren 1997; Maiti & McGrother 2004)—and molecular dynamics simulations using the MARTINI force field in conjunction with the GROMACS (Groningen Machine for Chemical Simulations) package (Marrink et al. 2007).

DPD simulations group atoms and molecules into fluid beads and use bead level interactions to describe the evolution of a system (Maiti & McGrother 2004). For any two beads i and j, the pair wise interaction force (Embedded Image) is the sum of the conservative force (Embedded Image), dissipative force (Embedded Image) and random force (Embedded Image) Embedded Image 1.1 Embedded Image is a soft repulsive force, while Embedded Image is a drag force or frictional force and Embedded Image a random force. Of these three, the conservative force (Embedded Image) best describes the energy of the system (Maiti & McGrother 2004). Groot & Warren (1997) established a connection between DPD beads and a real fluid by defining a relationship between the maximum repulsion between particles ((aij,) which is a function of the conservative force (Embedded Image)) and the Flory–Huggins interaction parameter (χ).

The MARTINI force field uses, on average, a 4 to 1 mapping of non-hydrogen atoms to interaction centres (although sometimes fewer or more than four atoms are mapped on to an interaction site) and defines the interaction sites into four main types: polar (neutral water soluble atoms), non-polar (mixed groups of polar and apolar atoms), apolar (hydrophobic groups) and charged (groups bearing an ionic charge; Marrink et al. 2004). Non-polar and charged particles are further divided into subtypes that allow a more accurate description of the chemistry of the underlying structure of the representative atoms. Subtypes are distinguished by a letter or number denoting hydrogen-bonding capabilities (d = donor, a = acceptor, da = both donor and acceptor and 0 = no hydrogen-bonding capability). The interactions between the various interaction sites are then described using non-bonded Lennard-Jones potentials (Marrink et al. 2004). Interactions between particles are represented by a number (I–V) where a lower number corresponds to a more attractive polar character (I) and a higher number (V) to a more repulsive (between apolar and polar units) character. Further interactions by charged groups are described using the electrostatic potentials (Marrink et al. 2004) and the interactions by covalently connected interaction sites by standard bond and angle potentials (Marrink et al. 2004).

With both coarse-grained simulations, parametrization of the particles is followed by the computational simulation under defined environmental conditions and over a specified number of time steps. The objective of the current study was to explore the application of these simulation techniques as potential predictive tools for polymer drug encapsulation, by simulating the encapsulation of a range of model drugs within poly(l-lactic acid) (PLA) and quaternary ammonium palmitoyl glycol chitosan (GCPQ) polymer matrices. At the outset, we chose to simulate meaningfully the largest length and time scales possible, as is appropriate for pharmaceutical formulation studies; DPD simulations were thus selected and applied to the homopolymer PLA. Coarse-grained molecular dynamics simulations enable the modelling of molecules with some degree of chemical heterogeneity, since a single molecule may be represented by a collection of chemically distinct beads (including beads with ionic character) and the MARTINI force field was thus selected for the simulation of GCPQ formulations; GCPQ is a polymer with an ionic functional group. PLA is used extensively in pharmaceutical formulations to control the release of drug compounds and thus control bioavailability (Venkatraman et al. 2000; Boussou & van der Walle 2006; Rowe et al. 2006), encapsulating low molecular weight drugs, therapeutic proteins, vaccines and DNA within the interior of microspheres or nanospheres, which are then released by erosion or diffusion from the particle in a controlled manner. GCPQ is a bioavailability enhancing polymer with a proven ability to increase the brain activity of drugs by up to 10-fold (Qu et al. 2006). This amphiphilic polymer self-assembles into polymeric micelles which then encapsulate hydrophobic drugs and control their biodistribution (Qu et al. 2006).

2. Experimental methods

All materials were obtained from Sigma Aldrich, Dorset, UK, unless otherwise stated.

2.1. Polymer microspheres

Microspheres were prepared from PLA (figure 1) polymers with different molecular weights: PLA 1 (Mw = 67 kDa) and PLA 2 (Mw = 102 kDa).

Figure 1.

Quaternary ammonium palmitoyl glycol chitosan (GCPQ) and poly(lactic acid) (PLA).

PLA microspheres encapsulating the hydrophobic drugs prednisolone and paracetamol were prepared using the oil-in-water (o/w) solvent evaporation method (Arshady 1990). Polymers (50 mg ml−1, 1 ml) dissolved in dichloromethane and appropriate volumes of drug solutions (to give the desired initial concentration of drug) in methanol (prednisolone = 10 mg ml−1, paracetamol = 100 mg ml−1) were slowly poured into an aqueous polysorbate 80 dispersion (2 mg ml−1, 5 ml) and the mixture probe-sonicated (1 min with the instrument set at 80% of its maximum output, Soniprep 150, Sanyo, UK). The resulting o/w emulsion was left shaking overnight in a water bath (90 r.p.m., 37°C) to enable evaporation of the organic phase, with the residual organic solvent being removed under a stream of nitrogen at room temperature. The microspheres were recovered by centrifugation (13 000 r.p.m. × 10 min, Z 323 K bench top centrifuge, Hermle, Germany), washed with distilled water (3 ml) and freeze dried.

PLA microspheres encapsulating the hydrophilic model drug isoniazid were prepared using a double emulsion technique. An appropriate volume of an aqueous solution of isoniazid (100 mg ml−1) to give the desired initial concentration of drug was added to a solution of the polymer (50 mg ml−1, 1 ml) in dichloromethane and emulsified using probe sonication for 1 min. This primary emulsion was then poured into an aqueous solution of polysorbate 80 (2 mg ml−1, 5 ml) and the double emulsion probe-sonicated for 1 min. The resulting water in oil in water (w/o/w) emulsion was left in a shaking water bath overnight to enable the evaporation of the organic solvent (90 r.p.m., 37°C). The residual solvent was dried under a stream of nitrogen and the isoniazid microspheres isolated as detailed above for prednisolone and paracetamol microspheres.

Drug encapsulation in the microspheres was analysed by dissolving each formulation of dried microspheres (5 mg) in dichloromethane (1 ml) followed by the addition of methanol (2 ml) to precipitate the polymer. The samples were then filtered (0.45 µm), an aliquot (100 µL) dried under a stream of nitrogen and the residue dissolved in the mobile phase (64% v/v acetontrile for prednisolone, 10% methanol for paracetamol and 5% isopropanol in ammonium formate (0.08 M) for isoniazid) containing an appropriate internal standard (6-methyl prednisolone for prednisolone, nicotinamide for paracetamol and paracetamol for isoniazid). Samples were then chromatographed over a reverse phase HPLC column (C18 symmetry—75 × 4.6 mm (Waters, UK) for prednisolone and paracetamol and a Cyanol (CN) 250 × 4.6 mm (Phenomenex, UK) column for isoniazid) using a Waters 717 autosampler and a Waters 515 isocratic pump. Samples were detected using a Waters 486 variable wavelength UV detector (λ = 243 nm for prednisolone, λ = 254 nm for paracetamol and λ = 262 nm for isoniazid). Data were analysed using Empower software (Waters, UK). The encapsulation efficiency was expressed as a percentage of the initial drug encapsulated in the microspheres.

To study the release of drug compounds from polymer matrices, drug loaded microspheres (5 mg) were suspended in phosphate buffered saline (PBS, pH = 7.4, 10 ml) and incubated at 37°C in a shaking water bath (30 r.p.m.). At regular time intervals an aliquot of the samples (0.5 ml) was withdrawn, centrifuged (13 000 r.p.m. × 30 min) and the supernatant analysed for drug content using the HPLC methods outlined above. Fresh buffer (0.5 ml) was added to replace the withdrawn sample.

2.2. Quaternary ammonium palmitoyl glycol chitosan micelles

Glycol chitosan (GC; 2 g) was degraded by heating at 50°C in hydrochloric acid (4 M, 150 ml) for 48 h (Wang et al. 2001). GCPQ was prepared from the degraded GC as previously described (Qu et al. 2006) by reacting GC (500 mg) in sodium bicarbonate (376 mg) dissolved in water (50 ml) with a solution of palmitic acid-N-hydroxysuccinimide in ethanol (5.3 mg ml−1, 200 ml) to form palmitoyl glycol chitosan (PGC). PGC was then reacted in N-methyl-2-pyrrolidone (63 ml) with methyl iodide (1.1 ml) in the presence of sodium hydroxide (100 mg) and sodium iodide (100 mg). Yield = 334 mg.

The molecular weight of GC was determined by gel permeation chromatography and laser light scattering and the structure of GCPQ was confirmed using 1H NMR.

Propofol was encapsulated within GCPQ micelles by adding GCPQ and propofol (20 mg) to water (1 ml), followed by probe sonication for 1 min and filtration (0.45 µm). To analyse the level of propofol encapsulated, formulations were diluted in the mobile phase (80% v/v methanol) and samples chromatographed over a 250 × 4.6 mm ODS2 column (Waters UK) using a Waters 717 autosampler, a Waters 515 isocratic pump and a Waters 486 variable wavelength UV detector (λ = 229 nm).

Propofol–GCPQ formulations and GCPQ dispersions were imaged using negative stained transmission electron microscopy (TEM). The specimens were mounted on a carbon-coated plastic grid, allowed to settle and then covered with a drop of trehalose followed by a drop of uranyl acetate (1% w/v) and imaged using a Zeiss 912 AB energy filtering transmission electron microscope (EFTEM) operating at 120 kV.

2.3. In silico methods

2.3.1. DPD simulation

Mesoscale simulation required identifying and defining the chemically distinct components in the system, defining the coarse-graining parameters and defining the interaction parameters between the various chemical species. Drugs, polymer, surfactant and solvent were each represented as distinct bead types. The bead number/chain length of polymers was determined using Embedded Image 2.1 where Mp is the polymer molecular weight, Mm the monomer molecular weight and C the polymer characteristic ratio.

Bead–bead interactions were entered as DPD repulsion units (aij): Embedded Image 2.2 where χ is the Flory–Huggins interaction parameter (van Krevelen 1997), given by Embedded Image 2.3 where Vref is the arithmetic mean of the molar volume of interacting species, δ the solubility parameter of each specie, R the real gas constant (8.314 J K−1 mol−1) and T the temperature in K. The solubility parameters δ were determined using the additive group contribution technique (van Krevelen 1997; table 1) and the DPD repulsion parameters also calculated (table 2).

View this table:
Table 1.

Mesoscale topology and solubility parameters of the various components.

View this table:
Table 2.

DPD repulsion parameters (aij) between species pairs.

DPD simulations of drug–polymer interactions were performed using the Materials Studio modelling software (v. 3.0, Accelrys Inc). The size of the simulation cell used was 30 × 30 × 30 DPD units. A bead density of 3 was used. The temperature of the simulation was taken to be the default value of 300 K and the duration of the simulation was 10 000 DPD time steps (24.4 ns). Encapsulation efficiency in simulations was determined by counting the number of drug beads entrapped by the polymer in relation to the total number of drug beads in the system. Drug, polymer ratios/concentrations were simulated by using comparative drug, polymer and water molar ratios.

Simulations of in vitro release were performed by forcing the model drugs into the polymer matrix using artificial χ values and then restarting the simulations with real χ values. The χ values used for forcing the drug into the polymer were based on the hypothesis that negative χ value corresponds to a favourable interaction (Groot & Warren 1997). A smaller simulation cell was used (15 × 15 × 15 DPD units) in order to facilitate drug entrapment within reasonable time. Simulation durations for the simulated release experiments ranged from 24.4 to 244.3 ns.

2.3.2. Molecular dynamics simulations

The MARTINI force field was used to model the interaction of GCPQ with propofol in aqueous medium (Marrink et al. 2007) and simulations were performed using the GROMACS simulation package v. 3.0 (Lindahl et al. 2001). The GCPQ monomer was modelled using 16 CG sites (figure 2). The palmitoyl chain (representing 15 methylene groups) was modelled based on the CG model for hexadecane (Marrink et al. 2004). A 4 to 1 mapping was used for the palmitoyl chain and a 2 to 1 mapping for the chitosan ring. This 2 to 1 mapping approach for the ring was first validated by coarse-graining glucose molecules and simulating their interaction with water; the coarse-grained model of glucose was based on the structural comparison to atomistic models. The quaternary ammonium ion group was modelled as a positively charged Q0 particle type (subscript 0 is used to denote the absence of hydrogen bonding capability; Marrink et al. 2004). The carbonyl group (linking the chitosan ring and palmitoyl chain) was modelled as a non-polar particle type Na (the subtype a was chosen because of the hydrogen bond acceptor capability of the carbonyl oxygen).

Figure 2.

(a) Mapping in the coarse-grained parametrization of GCPQ; (b) a coarse-grained model of GCPQ using 16 coarse grain interaction sites: four palmitoyl group interaction sites (green), one carbonyl group interaction site (brown), six chitosan ring interaction sites (orange), four glycol side chain interaction sites (red) and one quaternary ammonium group interaction site (blue). (c) Mapping in the coarse-grained parametrization of propofol; (d) a coarse-grained model of propofol consisting of five interaction sites, two sites for the ring (green), two sites for the methyl groups (green) and one site for the hydroxyl group (purple).

In the palmitoyl chain, an angle potential with an equilibrium bond angle of 180° and a force constant of kangle = 25 kJ mol−1 rad−2 was used; this angle potential reproduces the properties of aliphatic chains. An angle potential with a smaller equilibrium bond angle of 120° and a force constant of kangle = 40 kJ mol−1 rad−2 was used to model the chitosan backbone. The linkage in the chitosan backbone was defined between particle SP1 and SN0 in terms of an equilibrium distance (Rbond = 0.28 nm) and a force constant (kbond = 5000 kJ mol−1 nm−2). Improper dihedral angle potentials, Vid(θ), were introduced in the chitosan rings to keep the rings planar: Embedded Image 2.4 where θ = the angle between the planes constituted between atoms i, j, k and j, k, l, θid = the equilibrium angle and kid = the force constant of the improper dihedral.

The model drug propofol was modelled using five CG sites, three beads in a triangular configuration representing the phenol ring and two beads representing the isopropyl side chains (figure 2). The ring had particle types SP2 (representing ROH) and the other two particles were both SC2 types (representing C = C; Marrink et al. 2004). The isopropyl side chains were modelled by particle type C1 (subscript 1 denotes least polar character; Marrink et al. 2004). In the phenol ring, the distances between the three particles were constrained in order to prevent bond vibrations. The force field included the set of bonds, angles and dihedral potentials used to keep the ring structure rigid.

The force fields for GCPQ and propofol were generated and simulations were performed at a temperature of 300 K. A time step of 10 fs was used and the total simulation duration was 3 ns (300 000 steps). The system size was 7.740 × 7.898 × 7.340 nm (4128 CG particles in total) and an appropriate number of chloride counter ions (counter to the quaternary ammonium groups) were added to the simulation cell in order to preserve overall charge neutrality. Two different systems were modelled: system I was used to model the GCPQ–water interaction and system II was used to simulate the interaction of model drug propofol with the GCPQ–water system. System I consisted of eight GCPQ polymer fragments, each being a string of eight monomers (64 GCPQ monomers in total) in aqueous medium. System II consisted of 20 propofol molecules and eight GCPQ polymer fragments (n = 8) in aqueous medium. Water was modelled by coarse graining four molecules into one bead; the properties of bulk water have previously been reproduced using this CG approach (Marrink et al. 2004).

The output trajectories were analysed to determine some of the properties of the system, namely the energy of the system, system density (drugs and micelles) and the radius of gyration of the micelles (both in the presence and absence of drug).

3. Results and discussion

3.1. DPD simulations

The encapsulation of prednisolone within PLA microspheres was studied at an initial prednisolone, polymer concentration of 5 : 50 mg ml−1 and computational data revealed an increase in prednisolone encapsulation with an increase in molecular weight; a trend in agreement with experimental findings (figure 3). However, there was an underestimation of the actual level of prednisolone encapsulated within the polymer matrices in the computational data (figure 3a). Higher initial levels of prednisolone could not be studied experimentally as the unencapsulated prednisolone precipitate could not be easily separated from the microspheres. The encapsulation of paracetamol within PLA microspheres was studied at a fixed polymer concentration of 50 mg ml−1 and at varying initial drug concentrations (5, 9 and 12 mg ml−1). At initial paracetamol levels in excess of 5 mg ml−1, the simulations returned drug encapsulation levels that were in remarkable agreement with the experimental data (figure 4). It is possible that the coarse-graining technique, by its very nature, hampers the study of low drug concentrations effectively. However, useful data may be gleaned from the studies carried out at a drug concentration of 5 mg ml−1. For example, with the less hydrophobic drug, paracetamol, experimental and simulated paracetamol encapsulation results were lower than with the more hydrophobic prednisolone (figures 3a and 4a). During preparation of the microspheres, hydrophobic drugs will bind to hydrophobic polymer matrices as both are unable to hydrogen bond with aqueous media, and to varying degrees, dependent on drug physical chemistry, such hydrophobic drugs could become entrapped within precipitating polymer microspheres. These computational outputs indicate that the relative encapsulation efficiencies of a polymer for a variety of drugs may be successfully modelled. In experimental studies, there was no change in the per cent paracetamol encapsulated within the microspheres as the initial concentration of the drug increased, although the amount encapsulated ranged from 0.44 to 0.74 mg ml−1 (figure 4a).

Figure 3.

(a) The encapsulation of prednisolone within various PLA microspheres: initial polymer, prednisolone levels were 50 mg ml−1 and 5 mg ml−1, respectively. Numbers in parentheses are the concentration of prednisolone (mg ml−1) in the final experimental formulation; the asterisk represents statistically significant difference between experimental encapsulation efficiencies (p < 0.05); filled bar, experimental; open bar, simulation. (b) Simulated distribution of prednisolone (green) in PLA (red) microsphere formulations: initial polymer, prednisolone levels were 50 mg ml−1 and 5 mg ml−1, respectively; equilibrium encapsulation efficiency was 11.1% and 33.3% for PLA 1 and PLA 2 microspheres, respectively. Polysorbate 80 is represented by a pair of blue and purple beads, while water is represented as small light blue dots.

Figure 4.

(a) The encapsulation of paracetamol within PLA 1 microspheres; initial polymer levels were 50 mg ml−1, experimental encapsulation efficiencies are not statistically significantly different (p > 0.05; open square, PLA 1 simulation; filled square, PLA 1 experiment). (b) Simulated distribution of paracetamol (green) in PLA 1 (red) microspheres at various initial drug loads; initial polymer levels were 50 mg ml−1 and equilibrium encapsulation efficiencies were 22%, 5.9% and 8.3% when initial drug levels of 5 mg ml−1, 9 mg ml−1 and 12 mg ml−1 were used, respectively. Polysorbate 80 is represented by a pair of blue and purple beads, while water is represented as small light blue dots.

From the foregoing, it is apparent that DPD simulations may be used to predict the encapsulation within PLA polymers of hydrophobic drugs such as prednisolone and paracetamol and even predict the amount of drug that could be encapsulated in the case of paracetamol as well as the relative encapsulation efficiencies for a range of drugs. Such computational methods may be used as an initial screening tool for polymers prior to experimental encapsulation procedures.

The encapsulation of hydrophilic drugs within hydrophobic polyester microspheres is problematic (Odonnell & McGinity 1997), as hydrophilic drugs remain in the aqueous phase during processing and hydrogen bonding with water prevents co-precipitation in the polymer matrix or association with the surface of the polymer microsphere. With the hydrophilic drug isoniazid, encapsulation levels were low (1–3%) apart from when the initial isoniazid concentration was increased to 12 mg ml−1 (figure 5). With an initial drug level of 12 mg ml−1, experimental encapsulation efficiencies were between 4 and 9 per cent, depending on the polymer type. The limitation of the DPD simulations method was seen when attempting to model PLA interactions with the hydrophilic drug isoniazid, in that simulations predicted zero drug encapsulation (figure 5). The solubility parameters of water (δ = 44.75) and isoniazid, calculated using the additive group contribution technique (van Krevelen 1997; δ = 47.05), were similar, ensuring a preferred association between isoniazid and the aqueous media during the simulation run. It could be argued that low levels of encapsulation may be difficult to observe using mesoscale simulations due to the coarse graining of multiple molecules to form one bead; however, encapsulation efficiencies of 6 per cent were observed in the case of paracetamol simulations (figure 4), but the simulation could not reproduce the experimental data of 9 per cent isoniazid encapsulation efficiency when the highest initial isoniazid levels were used (figure 5). As the simulation could not predict the encapsulation of low levels of isoniazid, when isoniazid beads were added to PLA beads in the presence of water beads, we sought to confirm the inability to encapsulate the hydrophilic drug by investigating if the drug would be released from the polymer matrix if encapsulated prior to the simulation's onset. Isoniazid was forced to associate with the microspheres (using a false parametrization, which promoted association of isoniazid with the polymer), followed by the input of real isoniazid–water–PLA 1 parameters (figure 6). As expected simulated release was rapid and complete within 96 ns (figure 6d), confirming the inability of hydrophilic drugs to be encapsulated within PLA 1 matrices in DPD simulations. In experiments 87 per cent of isoniazid was released within the first 1 h (figure 6a).

Figure 5.

The encapsulation of isoniazid within various polymers at various initial isoniazid drug levels; initial polymer levels were 50 mg ml−1 (filled square, PLA 1; filled triangle, PLA 2). Inset: simulated distribution of isoniazid (green) within PLA 2 (red) microspheres; initial isoniazid level = 12 mg ml−1, initial polymer level = 50 mg ml−1 and 0% encapsulation was recorded at equilibrium. Polysorbate 80 is represented by a pair of blue and purple beads, while water is represented as small light blue dots.

Figure 6.

(a) The release of isoniazid from PLA 2 microspheres; 1.05 mg ml−1 of isoniazid was encapsulated by 50 mg ml−1 PLA 2 microspheres at the start of the release experiment. (b) DPD simulation of isoniazid (green) release from PLA 2 (red) microspheres after 72 ns; (c) DPD simulation of isoniazid (green) release from PLA 2 (red) microspheres after 96 ns; (d) DPD simulation of isoniazid (green) release from PLA 2 (red) microspheres after 96.002 ns, isoniazid released as a burst at just over 96 ns.

It appears that DPD simulations are deficient in predicting the polymer microencapsulation of water soluble drugs such as isoniazid: drugs with negative log p-values.

It should be noted that DPD modelling of poly(lactic acid-co-glycolic acid)–drug interactions was not well correlated to experimental data (data not shown). Furthermore attempts to use DPD simulations to study GCPQ–drug interactions revealed a difficulty in modelling the properties of an ionic functional group.

3.2. MARTINI simulations

In the coarse-grained simulation, eight GCPQ polymer fragments, each with a degree of polymerization (DP) of 8, self-assembled to form micelles on addition to aqueous medium (figure 7b,c). The hydrophobic palmitoyl chains comprised the micelle core, while the ethylene glycol and quaternary ammonium units constituted the surface of the micelle (figure 7c). The aggregation process (involving 1024 GCPQ beads and 3104 water beads) was very rapid and equilibrium was reached within 3 ns. Two micelles were formed, one micelle comprising six GCPQ polymer chains (figure 7b,c) and one micelle comprising two GCPQ polymer chains (data not shown), with the simulation suggesting that the micelle population could be heterogeneous in aggregation number, although longer equilibration times may be required before firm conclusions may be drawn.

Figure 7.

(a) An unfiltered dispersion of GCPQ micelles in water (8 mg ml−1; scale bar, 200 nm). (b) Snapshots taken at the end of the simulations showing the molecular arrangement in GCPQ micelles composed of six polymer fragments (DP = 8), either alone or in the presence of propofol (yellow and purple); in the three-dimensional view, the chitosan backbone (dark yellow) along with the hydrophilic glycol (red) and quaternary ammonium ion (blue) groups form the micelle surface while the palmitoyl chains (green) comprise the micelle hydrophobic core, and chloride counter ions (dark brown) and water molecules (blue) are seen in the background. (c) A slice through the micelle, key as in (b). (d) A snapshot of a propofol containing GCPQ micelle showing very little incorporation of propofol in the micelles. (e) A snapshot of a propofol loaded GCPQ micelle showing an abundance of propofol molecules loaded into the micelle; propofol is located mainly at the interface between the hydrophilic and hydrophobic groups. (f) Density plots showing the average propofol density within the micellar systems, micelle 1 (the larger micelle) possesses a lower drug density (propofol 1) when compared with micelle 2 (the smaller micelle, propofol 2; black line, micelle 2; red line, propofol 2; blue line, micelle 1; green line, propofol 1). (g) An unfiltered formulation of GCPQ (8 mg ml−1) and propofol (5.4 mg ml−1) nanoparticles (arrow 2) in water, relatively empty GCPQ micelles are also present (arrow 1; scale bar, 200 nm). (h) The experimental encapsulation of propofol within GCPQ micelles.

Once stable micelles had been formed, 20 molecules of propofol were introduced into the GCPQ–water system. As the simulation proceeded, drug molecules were seen diffusing into the micelle until the simulations converged (as determined by the total kinetic and potential energy achieving minimum values). Snapshots of the final structures obtained in the propofol–GCPQ–water simulations are shown in figure 7d,e. Propofol was seen to be associated with the interface between the hydrophilic and the hydrophobic parts of the micelle, presumably due to hydrogen bonding between the phenolic hydroxyl group and the chitosan monomers. Another interesting observation was that the smaller micelle (figure 7e) encapsulated more propofol compared with the larger micelle (figure 7d). The density profiles (figure 7f) were used to derive the average density of the micelles and density of the drug within the micelles. The density plot confirmed that the larger micelle (micelle 1) contained less drug when compared with the smaller micelle (micelle 2; figure 7f). The radius of gyration of the micelles (Rg) also increases on encapsulation of propofol (table 3).

View this table:
Table 3.

The radius of gyration (Rg) of GCPQ micelles.

In experimental studies GCPQ was synthesized with a level of palmitoylation of 30 mol% and a level of quaternization of 8.9 mol%, The levels of palmitoylation and quaternization were determined by 1H NMR analyses (Qu et al. 2006). The DP was 32 as determined from GPC-MALLS analysis on the parent polymer GC. This polymer is known to form 20–30 nm micelles in aqueous media (Qu et al. 2006; figure 7a). Electron microscopy images show that on encapsulation of propofol within GCPQ self assemblies, larger propofol nanoparticles are seen to coexist with smaller relatively empty micelles (figure 7g). Thus, the GROMACS simulation provided information that is coherent with experimental data (both the self-assembly of GCPQ into micelles and the heterogeneity of drug distribution within the micelle population) and information that is not readily available from experimental data (the location of the drug at the boundary between the hydrophobic core and hydrophilic surface of the micelles). The simulation experiments also confirmed the high level of propofol encapsulation by individual GCPQ micelles (figure 7e); in experiments one mole of GCPQ encapsulates 35 moles of propofol (8 mg ml−1 of GCPQ with an estimated molecular weight of 8920 Da encapsulates 5.4 mg ml−1 of propofol; figure 7h). GCPQ–propofol formulations are of pharmaceutical interest as GCPQ is able to increase the activity of propofol by 10-fold, with high levels of drug encapsulation playing a role in the mechanism of action of this system (Qu et al. 2006).

4. Concluding remarks

In this work, two types of coarse-grained approaches to modelling polymer–drug interactions have been explored: mesoscale modelling and coarse-grained molecular dynamics simulations. Mesoscale modelling, which constrains multiple molecules into beads (DPD), is useful for predicting the interaction of hydrophobic drugs with PLA polymer matrices, correctly revealing that the drug–polymer associations increase with polymer molecular weight and with the initial amount of hydrophobic drug used, provided that sufficient drug is used during the encapsulation process. With hydrophilic drugs, the drug model proves to be too hydrophilic to associate with the polymer matrices and we conclude that DPD simulations are not useful for modelling the interaction of hydrophilic drugs with PLA matrices. The DPD approach was useful for modelling homopolymers such as PLA (mesoscopic modelling of PLGA–drug interactions were not well correlated with experimental data—data not shown) but was not appropriate for the polyelectrolyte copolymer GCPQ, and not just because GCPQ is not a homopolymer but chiefly because DPD modelling does not take electrostatic interactions (such as those arising from the quaternary ammonium ion) into account. We thus used coarse-grained molecular dynamics simulations (the MARTINI force field) to model GCPQ. This method, which constrains atoms into a series of chemically distinct interaction sites, has been used to accurately model the self-assembly of GCPQ in aqueous media. The model indicates that propofol resides at the interface between the hydrophobic and hydrophilic portions of the micelle and that the drug propofol distributes heterogeneously within GCPQ micelles; the latter also evident from electron microscopy studies following drug encapsulation experiments.


The authors gratefully acknowledge the financial assistance from Astra Zeneca and the University of Strathclyde. The authors would also like to acknowledge Siewart Jan Marrink of the University of Groningen for his help with the coarse-grained modelling work involving GROMACS.


  • One contribution to a Theme Supplement ‘Scaling the heights—challenges in medical materials: an issue in honour of William Bonfield, Part I. Particles and drug delivery’.

  • Received March 30, 2010.
  • Accepted May 10, 2010.


View Abstract