|Home | About | Journals | Submit | Contact Us | Français|
A combination of Molecular Dynamics (MD) simulations and docking calculations was employed to model and predict polymer-drug interactions in self-assembled nanoparticles consisting of ABA-type triblock copolymers, where A-blocks are poly(ethylene glycol) units and B-blocks are low molecular weight tyrosine-derived polyarylates. This new computational approach was tested on three representative model compounds: nutraceutical curcumin, anti-cancer drug paclitaxel and pre-hormone vitamin D3. Based on this methodology, the calculated binding energies of polymer-drug complexes can be correlated with maximum drug loading determined experimentally. Furthermore, the modeling results provide an enhanced understanding of polymer-drug interactions, revealing subtle structural features that can significantly affect the effectiveness of drug loading (as demonstrated for a fourth tested compound, anticancer drug camptothecin). The present study suggests that computational calculations of polymer-drug pairs hold the potential of becoming a powerful prescreening tool in the process of discovery, development and optimization of new drug delivery systems, reducing both the time and the cost of the process.
Many drugs, including a number of leading anti-tumor agents, anti-depressants, statins, corticosteroids and vitamins, are hydrophobic and therefore require a solubilization process to enable their administration. While dispersing agents are customarily used to overcome this problem, they can contribute significant side effects beyond the inherent toxicity of the drugs.1 Recent progress in nanoparticle technology holds promise to improve the delivery of pharmaceuticals with poor bioavailability by improving their stability, circulation times in the body and permeability through cell membranes. Currently, the potential of nanoparticles for targeted drug delivery is being investigated for several life-threatening disorders, including cancer,2-5 HIV,6,7 and diabetes.8 One of the most promising technologies is the use of synthetic amphiphilic block copolymers that self-assemble into stable and safe (biodegradable and non-cytotoxic) nanoparticles,9,10 in which the core provides the reservoir for hydrophobic drugs11 while the corona enables aqueous stable dispersion and offers protection from protein adsorption and biological attack.12
In general, the development of such drug delivery systems is a wholly empirical process that begins by pairing a drug with a delivery vehicle (nanoparticle) based solely on their hydrophobic compatibility. However, it has been shown that stable incorporation of drugs into nanoparticles is also governed by other physical and structural factors, such as rigidity, conformation and configuration of the drug and the hydrophobic inner core-forming polymer block.13-15 Experimental capabilities remain limited, and often the failed experimental outcome (e.g. unstable formulation or low drug binding) is not fully understood and does not reveal ways to improve the drug/carrier system. Despite significant efforts to optimize nanoparticulate drug delivery systems, including new nanoparticle designs, materials, and drug encapsulation processes, fully functioning therapies based on nanoparticulate delivery systems are nascent at best,16 and the trial and error approach remains ineffective, time consuming and costly. A computational model that could predict, quantify and compare drug-particle interactions across different polymeric carriers a priori (relying solely on the chemical structures of the drug and polymer) would vastly reduce the costs and time associated with development of effective drug-delivery systems. Additionally, better understanding of the fundamentals governing drug-polymer interactions and binding could help explain the experimental results that, based on hydrophobic compatibility alone, are currently counter intuitive. Such insights could also provide further understanding of drug stability and release from nanoparticle-drug complexes.
In spite of this great potential, only a few groups have explored the capabilities of computational models for drug delivery.15,17-21 This at first seems surprising, since similar computational approaches are commonly used in the pharmaceutical field and are recognized as important tools in rational drug design, discovery and development. As such, various docking algorithms are used to predict small molecule interactions with proteins22 and large molecular databases are rapidly screened against certain targets (usually a protein).23-25 However, unlike proteins, polymers do not have a secondary/tertiary structure or specific binding sites for drugs, making the translation of knowledge from ligand-protein science more challenging.
Based on the challenges and limitations listed above, the aims of this study are to develop a computational method that provides (a) a better understanding of structural and physico-chemical factors that contribute to polymer-drug interactions and (b) ab-initio estimation of binding energies of polymer—drug complexes. This work represents the first step towards developing a scoring function, based entirely on ab initio computation, for rank ordering of drug-polymer pairs in terms of drug loading efficiency of the drug into polymeric nanospheres.
For the development of this computational model, we used our previously reported ABA-type amphiphilic triblock copolymers that self-assemble into nanospheres at low critical aggregation concentrations.9 The A-blocks (corona) of these nanospheres are composed of poly(ethylene glycol) PEG (5 KDa), and the B-blocks (core) are oligomers of desaminotyrosyl-tyrosine octyl esters (DTO) and nontoxic suberic diacids (number average molecular weight: ~ 19 KDa). The three model compounds, nutraceutical curcumin, anti-cancer drug paclitaxel, and pro-hormone vitamin D3, were selected based on their hydrophobicity and vast difference in the chemical structure. In this paper, for brevity, all three are denoted simply as “drugs”.
Molecular Dynamics (MD) simulations were combined with docking calculations to simulate, predict and better understand the interactions between each drug and the polymeric carrier. The entire surface of the simulated nanosphere core was scanned by blind docking, allowing for (a) qualitative evaluation of the physical interactions observed (i.e. types of interactions) and (b) quantitative assessment of the number of interactions and possible putative sites for each drug within polymer. For the three model drugs, the free energy of drug binding obtained computationally was correlated with the maximum drug loading as measured experimentally. Additionally, by applying the methodology to the anti-cancer drug camptothecin, which was previously reported to have a very low binding with tyrosine-derived nanospheres despite its relative hydrophobicity,9 a possible explanation to this counter intuitive behavior was revealed. The combination of MD and blind docking described in this paper has the potential of turning into a powerful screening tool that can accelerate the selection of suitable polymer-drug pairs from large virtual polymer libraries and drugs of interest, thereby speeding the development of drug-delivery stratagems for a variety of critical needs.
Drug-loaded nanospheres were prepared as described previously9 and their maximum stable drug loading was determined by increasing the initial drug feed concentration while keeping all other parameters constant (e.g. polymer concentration and encapsulation conditions). This experimentally determined drug loading (eq.1) is defined as percent ratio of the mass of strongly bound drug compared to the mass of the nanospheres.9
The first aim of this study was to provide a better understanding of the factors that contribute to polymer-drug interactions in nanoparticle drug delivery systems.
The DTO-SA repeat unit of the B-block of the ABA-type triblock copolymer (Scheme 1, DTO-SA/5K in 2D) was built in three dimensional (3D) coordinates using MOE's (Molecular Operating Environment software package, Chemical Computing Group, Canada) builder tool.26 The polymer chain was built with the head-to tail connection, also available in MOE, and the oligo(DTO-SA) block was terminated on both of sides with SA unit. The size of the polymer chain was scaled by 50% compared to the experimentally measured (based on the number average degree of polymerization as measured by Gel Permeation Chromatography)9 due to the software limitation regarding the number of atoms it can handle. Finally, the 5K methoxy-terminated PEG units (also scaled by 50%) were attached at both ends.
The 3D structure of the polymer chain was pre-optimized before running Molecular Dynamics (MD) simulations using the all atom MMFF94× force field with no constraints, which is parameterized for organic molecules. This particular force field was selected in order to avoid the inherent biases of other force fields commonly used in MOE (such as AMBER, CHARMM, Engh-Huber, etc.) that are calibrated with experimental data for proteins.
The self-assembly of the tyrosine-derived triblock copolymer was demonstrated via MD simulations performed on the DTO-SA/5K triblock copolymer (prepared as described at bullet (a)) in both implicit and explicit water environment, and using a canonical ensemble (NVT) with a target temperature of 300 K and an integration time step of 2 fs in a total of 5 ns simulation time. Given that the trajectory of time vs. potential energy converged significantly earlier than 5 ns, this time point was considered sufficient to generate a reliable conformation. In the case of explicit (fully) hydrated models the chains of the polymers looked slightly more packed and it was used to ensure that the solvent effect (e.g. hydrogen bonding with water molecules) was not omitted. The conformations obtained in fully hydrated conditions were used in subsequent studies. In order to account for the effect of the nanosphere size on the free energy of binding, four types of the polymers chains assemblies were used: one, two, four, and eight strands of DTO-SA/5K. These strands were initially placed in close proximity and parallel to each other, and then minimized (multiple chain MD was performed in all cases except for the one strand nanosphere). After the minimization step, the terminal PEG blocks were removed in order to reduce the computational times for the docking calculations. Since the nanosphere core acts as a sink for hydrophobic drugs, only the aggregated B-block units were considered for docking calculations. This assumption was confirmed by MD simulations, where it was observed that the A-blocks accumulated separately from B-blocks in water environment (supporting information Fig. S1). The partial atomic charges were automatically assigned by MOE based on a bond-charge increments method. Right before docking, they were re-calculated using Autodock Tools (ADT)27 as a requirement for docking calculations. ADT assigns Gasteiger partial charges. The orientation of the polymer strands used for docking calculations was exactly as obtained from the MD step described above, seemingly entangled chains with an overall spherical shape.
The 3D structures of the three drugs, curcumin (CUR), paclitaxel (PTX), and vitamin D3 (VD3) (Figure 1 A, B, and C, respectively) were pre-optimized using the same MMFF94× force field. Then, each drug was docked into the minimized, hydrated polymer structures using Autodock4 software package from The Scripps Research Institute.28 The files for docking calculations as well as the figures for this article were prepared using AutoDock Tools (ADT).27
Autodock uses a rapid grid-based method for energy evaluation. An initial grid refinement study was performed to insure convergence of the computed expected values (i.e. most probable value) of the binding energy. The size between the grid points in the grid space was varied from 1 to 0.125 Å and the energy values converged when the grid space was set between 0.125-0.5 Å29 (supporting information Fig. S2). Thus, for all docking calculations the size of the grids was varied between 64 × 64 × 64 Å and 126 × 126 × 126 Å using grid spaces between 0.375 - 0.5 Å for each core assemble. For the docking protocol, three built-in algorithms of Autodock4, Lamarkian Genetic Algorithm (LGA),30,31 Genetic Algorithm (GA),32 and Simulated Annealing (SA)33 were probed in order to find the most favorable drug-polymer complex geometry or, in other words, the most favorable interactions (supporting information Table S1). The LGA algorithm predicted the strongest drug-polymer affinity in all cases, and therefore, throughout this study, the LGA-computed expected binding energies were used. For all the simulations, the docking runs were set to 100. For GA and LGA searching algorithms the number of energy evaluation was set to 25,000,000 while the population size was set to 150. The other docking parameters were set to the default values for all three algorithms.
Our computational approach simulates the experimental study on a “minimal” scale (the actual nanosphere is orders of magnitude larger that the simulated one). The four assembles (one, two, four and eight polymer strands) were used in docking calculations. The estimated binding energy values converged in the assemblies made from four and eight strands, and thus to reduce the computational time/cost the four strands assembly was used as an optimal size to simulate the nanosphere (supporting information Table S1). All the binding energy results presented in this manuscript are generated based on this nanosphere size. Since the binding site(s) within the polymer are not defined, blind docking was applied to the entire nanosphere core represented by four strands. The Autodock program utilizes a three dimensional grid (in our case encompassing the entire surface of the assembly of four strands) to calculate the energy between a receptor (the core assembly) and a ligand (tested compound). An additional advantage of using blind docking is that it allows not only to quantify the possible physical interactions between the drug and strands assembly, but also to identify the number of putative binding sites for each drug. In contrast with proteins, which have unique binding sites due to different amino acids sequences, our nanosphere have the same repeat unit in the core, DTO-SA, and as a result there are several identical, or very similar, “hot spots” (putative binding sites) where the drug binds preferably (similar to the binding sites of proteins).
As mentioned above, in some cases the nanosphere drug loading is believed to parallel the hydrophobic compatibility between the drug and the polymer.3,14 The AlogPs values (octanol/water partition coefficient) for each of the three drugs, representing the quantitative measurement of the hydrophobicity of the drug, were calculated using ALOGPS 2.1 software package34,35 and are shown in Table 1, along with other molecular properties. As expected, the most hydrophobic drug, vitamin D3, showed the highest affinity to this particular oligo(DTO-SA) assembly, as evidenced by the highest drug loading determined experimentally (Table 1). For curcumin and paclitaxel, the hydrophobicity values are similar, and one would expect that they would have similar drug loading with DTO-SA/5K nanospheres. However, experimental results showed that the maximum amount of curcumin encapsulated is more than double that for paclitaxel. This experimental outcome can be explained if we consider the possible polymer-drug interactions, as revealed from our docking studies (Figure 2 (A-B)).
Figure 2 (A-D) shows the lowest energy conformation of each drug within core assembly after the in silico docking calculations. Curcumin (Figures 1A and and2A)2A) forms two H-bonds and four π-π interactions with the core. Curcumin's structure allows for almost an ideal alignment with the DTO unit of the polymer strands resulting in more inside tight binding in the nanosphere core. In contrast, paclitaxel (Figures 1B and and2B)2B) has fewer interactions: only one hydrogen bond and two π-π interactions with the polymer. In addition, paclitaxel is much bulkier (17 Å in diameter as opposed to only 5 Å for curcumin) and is more rigid due to the internal H-bonding and polycyclic structure and, therefore could mainly bind at the interface of the nanosphere core and corona, with a lower degree of inside binding. Interestingly, although vitamin D3 (Figures 1C and and2C)2C) has only one possible hydrogen donor/acceptor and no π-π interactions, its flexibility and extremely high hydrophobicity allowed for many hydrophobic interactions within the hydrophobic nanosphere core. The same complex is shown in Figure 2D, but in a surface view, where it is easier to observe the highly hydrophobic vitamin D3 buried in the polymer strands of the core assembly, confirming high incorporation into the hydrophobic part of the nanosphere. Overall, we found that the binding affinity of drugs to polymers depends not only on their hydrophobic compatibility, but also on physical interactions, such as hydrogen bonding, π-π stacking, or hydrophobic interactions. Additionally, structural features such as size and flexibility/rigidity can have a significant effect on drug loading into nanospheres (as seen it in case of paclitaxel).
The second aim of this study was to find a correlation between the experimentally measured drug loading of nanospheres, and the free energy of binding calculated by docking. The free energy of binding is defined as the energy released when two molecules are brought together to form a complex, and takes into account van der Waals, hydrogen bonding, electrostatic, entropy penalty, and desolvation component energies.30,31 In this study it was assumed that the drug is a small molecule compared to the size of the entire nanosphere, and thus it binds independently of any other drug molecule present. Furthermore, it was assumed that the binding energy field (as a function of position relative to the nanospheres and conformation of the drug) is a continuous function, and therefore low binding energies exist in the vicinity of the location defined as the lowest binding energy based on the AutoDock calculations. These two assumptions imply that the binding energy correlates with the mass of the drug per mass of the polymer when the maximum mass of drug per mass of polymer is reached. Based on above, the appropriate comparison for rank ordering of the drug binding to nanospheres (or polymer strands) is binding energy (BE) versus drug loading at maximum loading (see equation1).
During docking, the search algorithm identifies the global minimum in energy, or the most energetically favorable geometry of the complex. Therefore, for a fair comparison with experimental values, these values must represent the optimum experimental condition for a maximum of stable drug incorporation. This condition becomes crucial since experimentally it was observed that the polymer to drug ratio in the initial feed affects the LE significantly, as depicted in Figure 3. If all the drugs are tested experimentally at the same fixed feed ratio (6% w/w), curcumin has the highest drug loading (e.g. highest incorporation) (Figure 3-Grey). The polymer to drug ratio in the feed is one of the critical parameters.
However, the rank order of experimentally determined LE is changed when the optimized feed ratios were used and loading efficeiencies of of 9, 15, and 31% w/w were obtained for PTX, CUR, and VD3, respectively, (Figure 3-Black). Using these optimized conditions for each drug, a strong correlation (R2 = 0.93) is now observed in the rank ordering of computationally calculated binding energies (BE) and experimentally determined maximum drug loading (Figure 4).
The computational model correlates a drug's free energy of binding with drug loading in polymer strands assembly at maximum loading. The search algorithm used to calculate the free energy of binding identifies the lowest possible energy of the polymer-drug complex under optimized experimental conditions. When drug binding experiments are performed without carefully optimizing the experimental parameters, the correlation between the drug's free energy of binding and drug loading in nanospheres breaks down.
The strong correlation obtained in our study case (one polymer composition and three drugs) needs to be validated on a larger experimental set and possibly further refined. Nonetheless, our present results are encouraging and warrant further attention because they open the possibility of a rapid screening of large virtual polymer libraries against specific targets (small and large therapeutic molecules) in terms of binding affinity to a nanoparticulate vehicle. Thus, this approach will limit the number of polymers that need to be synthesized and characterized and can facilitate the optimization of the formulation and drug loading processes.
As stated above, a largely empirical method that selects drugs/polymer complexes based only on drug hydrophobicity, could not predict which drugs would be more effectively entrapped by the polymeric nanospheres, especially when the values of AlogP were similar. Our computational approach showed that, along with the physical interactions between drug-polymer pairs, the rigidity and conformation of the structures of both the drug and the polymer have a significant influence on the binding affinity and efficient incorporation of the drug into nanospheres. An example is camptothecin, an FDA approved anti-cancer drug, which is a small and mildly hydrophobic molecule (calculated ALogP =1.91, ALOGPS 2.1 program). It has a very rigid conformation with five associated rings and contains one possible hydrogen bond donor (Figure 5A). Previous findings showed poor encapsulation regardless of the polymer:drug initial feed ratio (experimental loading efficiencies of ~1 to 3 %).9 In contrast, our docking calculations predict a binding energy of -9.27 kcal/mol, which is indicative of a good binding affinity and would rank this drug between curcumin and vitamin D3.
As discussed above, during “blind docking” calculations, the search space consists of the entire volume occupied by the polymer strands assembly (representing the nanosphere's core). In the case of camptothecin, and unlike the other drugs tested, only an extremelly small part of the search space was preferred by the drug (supporting information Fig. S3). As presented in Figure 5B, in about 70% of the docking runs, camptothecin bound preferably in a single, certain spot, with the drug having same alignment. The consequence is that very few putative binding sites, or “hot spots”, within the DTO-SA-based nanosphere core are available for this molecule. To explain the scarcity of “hot spots” for camptothecin, we need to consider the drug-polymer pair not only in the context of their relative interactions and drug's rigidity, but also in the context of the kinetic competition between drug encapsulation and drug stability. Therefore, even though theoretically camptothecin can form one hydrogen bond with the substrate and multiple π-π interactions, its rigidity allows being accomodated only by certain polymer conformation. We suggest that, once this “hot spot” is occupied, it becomes much more difficult for camptothecin to bind because polymer rearrangement is a slower process and the rest of unbound molecules are loosly associated, if at all, in the interface of the nanosphere hydrophobic core and PEG corona. Since the nanospheres are dispersed in the aqueous media and camptothecin is not water soluble, the loosly bound molecules will disassociate from nanopshere complex. This explanation is confirmed by the experimental visual observations of yellow precipitates of camptothecin constantly formed in all camptothecin feed ratios. What was experimentally observed could be explained as follow. Hence, while favored thermodynamically, camptothecin's binding is unfavored kinetically due to the drug's structural rigidity. This observation can be the basis of a further refinement of our method, where the number of putative binding sites needs to be considered when predicting the drug loading of the drug.
The camptothecin study suggests that a simple correlation between the experimentally determined drug loading and the calculated binding energy is not sufficient to rank order outlier drugs. Thus, the development of a scoring function is needed to improve the method of searching for the effective drug delivery system. In this scope, more polymer systems (not only tyrosine-derived triblock copolymers) will need to be virtually screened against drugs of interest to help in refining the rank ordering of the drugs.
A computational model to predict the affinity between the drug and its nanoparticulate carrier, in this case DTO-SA/5K nanospheres, was developed by combining MD simulations with docking calculations. The model successfully rank ordered three drugs in terms of experimentally observed drug loading at maximum into the nanospheres, by comparing it with computationally calculated binding energies (BE). Moreover, the model was able to provide an understanding of additional factors that affect the drug binding and encapsulation. Thus far, this methodology was tested only on one polymer and four drugs with successful prediction outcome or explanation of the failing experimental outcome demonstrating the feasibility of this method. This modeling approach is not intended to substitute the experimental work, but to complement it to faster achieve the established aims.
Future work with additional drugs and polymer families will be undertaken to further improve the scope of this methodology, and in developing a scoring function to improve the method for search for the most efficient delivery system for a desired drug. The proposed method has the potential of becoming an important step in the discovery and development of new drug delivery vehicles for target drugs by offering the tools for a rational design of the of nanoparticles drug complexes that are based on fundamentals of physico-chemical properties of drugs and polymeric materials. If successful, this model will make it possible at the same, relatively short time (which is impossible to perform experimentally) to analyze a variety of different materials and potential drugs for addressing the needs of a variety of medical applications.
We thank Prof. Robert Latour and Dr. Xianfeng Li of Clemson University, Clemson, SC for fruitful discussions and suggestion for this current research. This work was supported by RESBIO (Integrated Technology Resource for Polymeric Biomaterials) funded by National Institutes of Health (NIBIB and NCMHD) under grant P41 EB001046. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH, NIBIB or NCMHD).
Supporting Information Available: The three dimensional (3D) structure of the ABA-type triblock copolymer, Grid refinement studies, Estimated binding energy values using three searching algorithms built in Autodock4, and Main putative binding sites of vitamin D3, curcumin, and paclitaxel docked into four relaxed strand of oligo(DTO-SA) core. This material is available free of charge via the Internet at http://pubs.acs.org.