4Center for the Development of Nanoscience and Nanotechnology (CEDENNA), Universidad de Santiago de Chile, Santiago, Chile
Associate Editor: R. Jin
Beilstein J. Nanotechnol. 2018, 9, 1328–1338. doi:10.3762/bjnano.9.126
Received 21 Dec 2017, Accepted 12 Apr 2018, Published 02 May 2018
Nanotherapeutics is a promising field for numerous diseases and represents the forefront of modern medicine. In the present work, full atomistic computer simulations were applied to study poly(lactic acid) (PLA) nanoparticles conjugated with polyethylene glycol (PEG). The formation of this complex system was simulated using the reactive polarizable force field (ReaxFF). A full picture of the morphology, charge and functional group distribution is given. We found that all terminal groups (carboxylic acid, methoxy and amino) are randomly distributed at the surface of the nanoparticles. The surface design of NPs requires that the charged groups must surround the surface region for an optimal functionalization/charge distribution, which is a key factor in determining physicochemical interactions with different biological molecules inside the organism. Another important point that was investigated was the encapsulation of drugs in these nanocarriers and the prediction of the polymer–drug interactions, which provided a better insight into structural features that could affect the effectiveness of drug loading. We employed blind docking to predict NP–drug affinity testing on an antiaggregant compound, cilostazol. The results suggest that the combination of molecular dynamics ReaxFF simulations and blind docking techniques can be used as an explorative tool prior to experiments, which is useful for rational design of new drug delivery systems.
Keywords: drug delivery; PEGylated nanoparticle; PLA; polymeric nanoparticle; reactive force field
In recent years, the use of drug delivery systems based on polymeric nanoparticles (NPs) has generated innovative therapeutic strategies for infection and immune diseases, as well as cancer therapy [1-3]. Polymeric NPs have shown significant advantages compared with many other nanosystems concerning their biodegradation, biocompatibility and highly modifiable physical and mechanical properties [1,4-6]. Polymeric materials also have demonstrated superior efficacy to deliver therapeutic agents at the site of the disease or damaged area and protecting drugs against in vitro and in vivo degradation [7-9].
Poly(lactic acid) (PLA) is one of the most commonly used polymers for the synthesis of NPs. PLA is a synthetic biodegradable, compostable and non-toxic polymer derived from renewable resources [10-12]. Despite its benefits for different formulations in the medical field , the applications of PLA are limited, mainly due to its weak hydrophilicity and low drug loading of polar drugs . However, PLA nanoparticles conjugated with hydrophilic molecules like poly(ethylene glycol) (PEG) are a singular kind of functionalized NPs. Several studies have shown that PEGylation of NPs allows improved blood circulation, clearance, biocompatibility and less cytotoxicity [15-19]. Hydrophilic polymer chains at the surface of NPs act as a steric barrier, reducing the opsonization and the subsequent phagocytosis [20-22]. This amphiphilic nature of the copolymer PLA–PEG supports its high potential for development in drug delivery  because the PLA core acts as a reservoir for hydrophobic drugs while the corona enables a good dispersion in the blood and offers protection against biological attack [23,24].
In addition to the hydrophilicity/hydrophobicity ratio, the size and surface chemistry are crucial factors in NP–cell interactions. Several studies on metallic NPs, carbon nanotubes, and dendrimers have shown that toxicity of different NPs is determined by a combination of these factors [25-29]. In this respect, computational modelling provides useful approaches to help elucidate the behavior at the atomic level of various nanosystems, predict their drug affinity, and at the same time, reduce the number of preliminary experimental tests (or to be used as a complementary approach to experimental work). The structural characterization in silico provides a comprehensive understanding of properties that govern the formation mechanism and drug loading of NPs and help create better designs and promising compositions. Although there are different experimental methods for synthesis and preparation of PLA–PEG NPs [8,11], molecular descriptions using computational methodologies still are not totally addressed.
Wang et al.  carried out dissipative particle dynamics (DPD) simulations to study the PEG–PLA–PEG/paclitaxel delivery system. The results showed a spherical micelle with the drug inside and polymer chains distributed on the surface of the micelle . In a recent investigation , it was found that different computational simulation strategies can successfully predict the experimental drug affinity and drug loading for PLA–PEG NPs.
The aim of this work is to characterize at the atomic level the structural properties of differentially charged polymeric NPs PEGylated with DSPE–PEG(2000) (1,2-distearoyl-sn-glycero-3-phosphoethanolamine-N-[poly(ethylene glycol)-2000]) carboxylic acid-, methoxy- and amino-terminated NPs prepared by the nanoprecipitation method and to study the NP–drug complexation for a model antiaggregant compound, cilostazol [32,33]. PEG chains typically create an interface between the NP core (PLA) and the hydrophilic environment, where the drug encapsulation is largely dependent on the intrinsic affinity between the drug and the PLA core . All-atom molecular dynamics (MD) simulations were performed using the reactive force field (ReaxFF) to prove that the self-assembly process of PLA–PEG NPs is consistent with the experimental method of synthesis , investigate the flexibility of polymer chains, and characterize the ability to protect a potential cargo. Docking calculations were employed in order to determine the spatial distribution of cilostazol in polymeric NPs and to explore its potential use in this kind of drug delivery system.
PLA polymer chains were built in three-dimensional coordinates with a head-to-tail connection and syndiotactic (D,L) configuration (Figure 1A). The molecular weight distribution of the polymer that was used in the experiments was not well characterized. Thus, considering the limitations on size and computational cost of all-atom MD simulations, the length of the PLA polymer was defined as 20 monomers (Figure 1B).
All MD simulations described in this study were performed using the molecular dynamics simulator, large-scale atomic/molecular massively parallel simulator (LAMMPS) [35,36]. The interactions between atomic species were represented with a novel reactive force field (ReaxFF) [36,37]. This force field considers both covalent and electrostatic interactions by employing a bond-order formalism and a polarizable charge description, enabling depiction of chemical bonding (bond formation and rupture) without expensive quantum mechanics calculations [38,39]. In this work, a modified version of the force field developed by Kim and van Duin, optimized for C/O/H/N/P systems , was used.
The convergence criteria for energy minimization, rupture temperature limit and the optimal separation distance for PLA chains to build the PLA core were defined as follows.
Stopping tolerance for energy: For a 20-monomer PLA chain, an energy threshold was defined to perform the energy minimization of the system. The stopping tolerance for energy was set as 1.0e−12 (means an energy tolerance of one part in 1012).
Rupture temperature limit: With the aim of creating an initial maximum temperature value (Ti) to be specified at the start of the equilibration phase, a temperature range from 300 to 600 K for a 20-monomer PLA chain was evaluated. The Ti was set as 600 K.
Optimal separation distance: To ensure an optimal package and minimal energetic repulsion between PLA chains, the optimal separation distance between them was evaluated. The separation distance between PLA chains was established in 7 Å (Figure 2). Thus, 10 PLA chains were packed in a 100 × 100 × 100 Å3 periodic box to assemble the PLA core by a MD simulation in vacuum.
The initially packed configuration was minimized using the Polak–Ribiere version of the conjugate gradient (CG) algorithm  with a convergence criteria of 1.0e−12 followed by an equilibration of 2 ns at 310 K using NVT ensemble (constant temperature and volume) and a timestep of 1 fs. A Nose–Hoover thermostat was used to control the temperature with a damping of 100 timesteps.
To corroborate the correct assembling methodology, several MD simulations of 100 ps were performed for 20, 40 and 80 PLA chains (Figure 1C). Data collected along the trajectories were used to calculate molecular properties such as a radius of gyration and sphericity with the aim of a better characterization of the PLA core shape.
Three PLA formulations were analyzed using DSPE–PEG(2000) carboxylic acid-terminated, methoxy-terminated and amino-terminated molecules (Figure 3) to provide different surface charges, which is a critical factor related to toxicity in nanosystems [25-27,42].
DSPE–PEG structures were built in three-dimensional coordinates. As the starting point, the DSPE–PEG chains were arranged around the PLA core (10 PLA chains) using PACKMOL [43,44], as follows.
Figure 4 shows the built atomistic model and represents the initial conformation for each case: carboxy-, methoxy- and amino-terminated PLA nanoparticle.
Then, this starting configuration was also minimized and equilibrated with the same conditions described above, obtaining the final PLA–DSPE–PEG-coated structure (Figure 5).
Finally, to analyze the stability of this structure, and at the same time reduce the computational cost, Langevin dynamics were performed on these final structures with a friction coefficient of 30 ps−1. In this way, the solvent particles are omitted, and their effects are represented by a combination of random forces and frictional terms. Simulations of 1 ns were carried out. The velocity Verlet-like algorithm was performed at the time of integration. The temperature was controlled at 310 K with a Langevin thermostat  and the pressure was controlled at 1 atm using a Berendsen barostat . The radius of gyration, sphericity and molecular distribution of PLA and DSPE–PEG groups were calculated from the trajectories.
Blind docking calculations were performed using AutoDock 4.2  software. The 3D structure of the drug (Figure 6A), was sketched and preoptimized; partial charges were assigned and rotatable bonds were identified. All files for docking calculations were prepared using AutoDock tools (ADT) . The orientation of polymer chains obtained from the above MD simulations was used to carry out the docking calculations. Due to the fact that the polymer core is the protective portion for drugs in this kind of nanocarrier , only the PLA core was considered for docking calculations.
AutoDock uses a rapid grid-based method for energy evaluation. A grid volume large enough to cover the entire surface of the PLA core was built (126 × 126 × 126 Å3) using a grid spacing of 0.5 Å. The grid parameters were generated using AutoGrid 4.2.6 and the Lamarckian genetic algorithm (LGA) was used to perform a search of the conformational space of the drug. The docking runs were set to 100.
The docking poses were analyzed by examining their binding energy score and the most visited “hot spots” (putative binding sites). The most energetically favorable conformations were selected as the best poses.
The radius of gyration (Rg) is related to the mean distance of atoms in a molecular structure with respect to its center of mass . The radius of gyration is a good indicator of the spatial conformation of the molecule. Thus, a higher value represents a higher spatial disposition. The radius of gyration was calculated for all PLA core models as a function of simulated time, as shown in Figure 7A.
The mean Rg were calculated, including the last 60 ps of simulation, ensuring a steady value for each PLA core conformation. For 10 PLA chains, a value of 15.9 ± 0.2 Å was obtained. Meanwhile, a value of 18.3 ± 0.2 Å, 23.5 ± 0.2 Å and 30.4 ± 0.1 Å was obtained for 20, 40 and 80 PLA chains, respectively.
To identify the grade of spherical shape for the PLA cores of NPs, a modified script developed by Paz et al.  was readjusted, which classifies the geometrical structures of PLA cores into three groups of spheroids: spheres, oblates and prolates . According to the described procedure by Paz et al.  to analyze the geometric deformations over time, it is necessary to define the triaxial parameter:
where a represents the largest axis of symmetry of the spheroid, b is the middle axis and c is the smallest axis. Based on this, κ = 0 indicates a sphere, κ > 0 indicates a prolate and κ < 0 denotes an oblate.
Figure 7B shows the evolution of κ vs time for the PLA core with a different number of PLA chains. The spherical character increased with a major number of PLA chains. It can be observed that the resulting shape is prolate (κ > 0) for 10, 20 and 40 PLA chains.
The sphericity of PLA cores depends directly on the number of chains, independent of their initial configuration. In that context, this behavior is represented by 40 PLA chains (red line) in Figure 7B, in which the initial configuration (at time 0) is more prolate compared with 20 PLA chains (pastel blue line), reaching values of κ = 0.64 and κ = 0.30, respectively. The initial configuration did not affect the final sphericity of cores because after only a few nanoseconds, both lines intersected. The spherical character of PLA cores reached a steady value for 40 and 80 PLA chains after 60 ps of simulation.
The radius of gyration can also characterize the relative DSPE–PEG-coated nanoparticle size and position of the constituent units. The radius of gyration for differentially charged (carboxylic acid-, methoxy- and amino-terminated) DSPE–PEG-coated PLA nanoparticles as a function of time is shown in Figure 8. For carboxylic acid-, methoxy- and amino-terminated molecules, the radius of gyration values were reduced from 62.7 to 20.5, 64.4 to 23.9 and 64.7 to 21.6 Å, respectively.
The different spatial distribution of NPs at the start and at the end of simulation demonstrates the DSPE–PEG-coated process. The initial stretched DSPE–PEG chains were folded around the PLA core, covering the polymer surface (Figure 5).
The mean Rg values were calculated from the last 40 ps of simulation, where the structure size appeared to stabilize. It was observed that the size of all those DSPE–PEG-coated PLA nanoparticles was homogeneous; however, a slightly larger particle size was found with methoxy-terminated chains (24.5 ± 0.6 Å versus 20.7 ± 0.2 Å for carboxylic acid-terminated and 21.8 ± 0.3 Å for amino-terminated NPs). Interestingly, those results were consistent with experimental results measured by dynamic light scattering . Nevertheless, since an all-atom and reactive force-field model was used, which increased computational cost, a smaller model was built maintaining the conformation of the whole system. This is the reason why discrepancies in scale (≈20 Å and ≈200 nm diameter for computational and experimental results, respectively) were observed.
To estimate the charge distribution for each NP, the distance between different terminal groups (COOH for anionic NP, CH3 for neutral and NH2 for cationic NP) and the center of mass of the PLA core was measured. Figure 9 shows the distribution of all terminal groups (ten) of each NP along the simulation.
At the start of the simulation, the terminal groups were observed within 120–160 Å from the center of mass for carboxy-, methoxy- and amino-terminated NPs, which corresponded to the stretched DSPE–PEG chains stage. The distance decreased as simulation time increased in all systems. At the end of the simulations (last 40 ps) the terminal groups were observed at 24.02 ± 5.8, 22.9 ± 6.6 and 24.1 ± 7.1 Å for carboxy-, methoxy- and amino-terminated NPs, respectively. Based on the final values of Rg, it can be suggested that all terminal groups are randomly distributed at the nanoparticle surface. Ensuring a correct PEGylation process is crucial. The surface design of NPs requires that the charged groups must surround the surface region for an optimal functionalization/charge distribution, which is a key factor in determining physicochemical interactions with different biological molecules inside the organism .
The interaction energies between the polymeric phase (PLA) and the lipid portion (DSPE–PEG) (Table 1) were similar in all cases and it is governed by hydrophobic interactions. Despite that van der Waals energy was similar in all formulations (−664,536.34, −661,500.71 and −659,315.18 kcal/mol for carboxy-, methoxy- and amino-terminated NPs, respectively), the slightly lower energetic contribution (total energy (Etotal), potential energy (PE) and kinetic energy (KD), and van der Waals (vdW) energy (EvdW)) obtained in the amino-terminated NPs could lead to significant differences in the experimental synthesis, producing particles with a nonuniform size and/or aggregates.
Table 1: Interaction energy (kcal/mol) between the PLA core and DSPE–PEG lipid portion.
|Nanoparticle||Etotal (PE + KE)||EvdW||EC|
Docking studies showed that cilostazol was trapped in PLA chains and could also be found at the interface of PLA and DSPE–PEG, interacting with both PLA chains and the nonpolar region of the lipids. The high affinity of cilostazol to the PLA core was consistent with the A log P value (octanol/water partition coefficient) of 3.38 calculated using ALOGPS 2.1 software package  and shown in Table 2. Table 2 and Figure 6B show the ten lowest energy conformations of cilostazol and their location within the NP obtained in the docking calculations.
Table 2: Ten lowest computed binding energy values for cilostazol. An A log P value of 3.38 was calculated.
|Binding energy (kcal/mol)|
This NP–drug noncovalent complex was established mainly by hydrophobic interactions, which would allow an easier cargo delivery as the polymer degrades.
It is important to note that this combined computational methodology could be particularly effective in cases where potential drugs have similar hydrophobicity values to achieve a first approximation to guide the experimental phase. It could be expected that a drug with a lower degree of hydrophobicity would be more likely to migrate to the surface of the NP, interacting with the terminal groups of lipid chains (for example through electrostatic interactions) or with the solvent, and thus, the protective function of the nanocarrier would be revoked. However, significant differences in experimental results (drug loading) have been reported when testing drugs with comparable hydrophobicity , which could be explained by NP–drug interactions calculated by these predictive tools.
In the present study, MD simulations provided a comprehensive understanding of the factors that contribute to PLA–PEG nanoparticle formation. To the best of our knowledge, this is the first time that the reactive force field ReaxFF was used to simulate polymer nanoparticle formation.
We have characterized the structure of differentially charged (anionic, neutral and cationic) polymeric DSPE–PEGylated nanoparticles. We detected a self-assembling process in which a core–shell structure is observed with PLA in the core and DSPE–PEG in the shell; this model is consistent with the nanoprecipitation synthesis method previously used . Methoxy-terminated NPs showed a slightly larger size compared to charged particles just as it was observed in the experimental results, and all terminal groups were distributed at the surface of the NPs.
The combination of MD-ReaxFF and blind docking techniques described in this work allowed the investigation of the antiplatelet drug encapsulation ability of polymeric NPs and can be used as a prior explorative tool or as a complementary approach to experiments, which is useful for the rational design of new drug delivery systems.
In future work, this approximation will be correlated with maximum drug loading determined experimentally and the protective role of differentially charged PEG for drug release process will be studied.
M.F.M. acknowledges support from CONICYT-PFCHA/Doctorado Nacional/2014-21140225. M.M.M. thanks the FONCyT PICT-2015-2191, CONICET PIP 11220110100992, Secyt, Universidad Nacional de Cordoba. C.V. acknowledges support from CONICYT under FONDECYT #1161438 and BASAL Grant FB0807, MECESUP PMI-UAB1301, and H2020-MSCA-RISE-2016 #734801 MAGNAMED. The authors thank the High-Performance Computational Center (CCAD UNC) and Escuela de Ingeniería Civil en Bioinformática (Universidad de Talca) for access to supercomputers.