Introduction

A continuing challenge in the development of anticancer therapeutics is ensuring that adequate amounts of the drug are delivered to tumors while simultaneously minimizing toxic and adverse effects to healthy tissue [1, 2]. As angiogenesis is crucial for tumor growth, the microvascular endothelial cell and its receptors have become key target modulators in cancer therapy [3]. The αVβ3 integrin and its ligands, vitronectin and fibronectin, are known to promote angiogenesis [4, 5], and tumor cells with non-activated αVβ3 integrins have been shown to exhibit decreased angiogenic behavior [6]. As a result, antagonists of αVβ3 integrins have emerged as a new class of anticancer therapeutics, primarily the Arg-Gly-Asp (RGD) tripeptide and RGD peptidomimetics, which bind to αVβ3 integrins with high affinity to prevent binding of αVβ3 integrins to their natural receptors [7–9].

This study focuses on the preclinical development of a nonpeptide RGD (npRGD)-based drug delivery system: poly(γ-glutamyl-glutamate) paclitaxel tethered to npRGD via poly(ethylene glycol) (PGG-PTX-PEG-npRGD) [10]. Paclitaxel (Taxol®, C47H51NO14) (PTX) is a hydrophobic drug used to treat breast, ovarian, and lung cancers, and poly(γ-glutamyl-glutamate) (PGG) is a hydrophilic polymer conjugated to PTX to improve solubility [11, 12]. To maximize access of npRGD to αVβ3 integrins, PEG is positioned between glutamyl-glutamate and npRGD residues to avoid steric interactions between PGG-PTX and npRGD (see Fig. 1). Conjugation of PEG to therapeutic agents has also been shown to result in improved in vivo stability, protection from proteolytic digestion, increased biological half-life, improved solubility, and decreased toxicity [13].

Fig. 1
figure 1

Chemical structures of poly(γ-glutamyl-glutamate) (GG), GG-paclitaxel (PTX), and GG-poly(ethylene glycol) (PEG)-npRGD [non-peptide residues RGD (Arg-Gly-Asp) peptidomimic]. GGNH3+ and GG-PTX each has a charge of −1; GG and GG-PEG-npRGD each a charge of −2; and GGCOO– a charge of −3. The GG residue differs slightly based on its position on the polymer, e.g., its position at the amino or carboxyl termini. For the GG-PTX residue, paclitaxel is conjugated to a carboxylate group of glutamyl-glutamate via an ester linkage. For the GG-npRGD residue, the nonpeptide RGD is attached to a linear PEG spacer, which is conjugated to the GG residue. The length n of the PEG spacer depends on the molecular weight of the linear PEG molecule: 500 Da; 1,000 Da; and 2,000 Da correspond to approximately 11, 23, and 45 ethylene glycol monomers

Full size image

Preclinical development of anticancer therapeutics usually involves trial-and-error testing of candidate compounds on in vitro and in vivo models. While these methods yield useful results, they are time-consuming, inefficient, and resource-intensive. Given to urgent clinical need for successful therapies there is clear motivation to shorten the course of preclinical development of therapeutic agents. In this study, we performed molecular dynamics (MD) simulations of PGG-PTX to predict their physicochemical properties, such as shape, flexibility, and surface hydrophilicity—factors that may influence how nanoparticles overcome biological barriers [14–18]. It was hypothesized that the density of nonpeptide RGD (4, 8, 12, and 16 ligands) and molecular weight (MW) of PEG spacer (500 Da; 1,000 Da; and 2,000 Da) concomitantly affects the shape, flexibility, and surface hydrophilicity of a PGG-PTX-PEG-npRGD molecule. (Figure 2 shows the patterning schemes for the PGG-PTX-PEG-npRGD molecules.) Theoretical insight from the modeling results may be useful to suggest better molecular designs of PGG-PTX-PEG-npRGD, ultimately expediting its preclinical development.

Fig. 2
figure 2

Abstract representations of the spatial positioning patterns of paclitaxel (PTX) and (PEG) n -nonpeptide RGD molecules on the PGG backbone. Patterning schemes apply to all PEG 500, PEG 1000, and PEG 2000 systems. Each PGG-PTX molecule is composed of 130 poly-γ-glutamyl-glutamate monomers, 26 paclitaxel molecules, and 4, 8, 12, or 16 (PEG) n -nonpeptide RGD molecules. Paclitaxel and (PEG) n -nonpeptide RGD molecules are attached covalently to the PGG backbone in a random fashion. Numbers between residues denote the number of repeating GG residues that are not amino- or carboxyl-termini GG residues. The amino- and carboxyl-termini GG residues are represented by black dots at the ends of each line

Full size image

Prior to running MD simulations, a key challenge was to determine the initial configuration of PGG-PTX-PEG-npRGD [19]. The combination of flexible polymeric components (PGG and PEG) and rigid components (PTX and npRGD) make PGG-PTX-PEG-npRGD a semiflexible molecule. Currently, there is no crystallized conformation representing its experimental form available from the Protein Data Bank. Therefore, circular dichroism (CD) spectroscopy was used to determine the presence of chiral forms of PGG-PTX-PEG-npRGD in aqueous form and obtain a general idea of its initial structure, which was then used to construct all-atom (AA) models.

In previous work on protein simulations, once the initial configuration has been determined, MD simulations were run to equilibrium, at which point a specific structure is attained [20, 21]. At this point, the potential energy level of a system became constant with respect to time. The approach must be different in the present case. Given the semiflexible nature of the PGG-PTX-PEG-npRGD, the molecule most likely exists as a population of related structures, rather than one specific structure. Therefore, the goal is not to attain a static PGG-PTX-PEG-npRGD conformation through MD simulations. Rather, MD simulations are run until minimal structural fluctuations at the molecular level, or a statistical equilibrium [22], is reached, at which point a PGG-PTX-PEG-npRGD molecule is characterized as an ensemble of statistically similar structures.

Given the relatively large sizes of each PGG-PTX-PEG-npRGD molecule (∼90–130 kDa), running MD simulations on AA models in explicit solvent was judged to require an extraordinary amount of CPU time and expense. Therefore, to minimize computational costs and allow access longer time scales, AA models were converted to coarse-grained (CG) models using the mapping scheme in the MARrink’s Toolkit INItiative (MARTINI) force field [23]. CG parameterization of the models was done by using the Boltzmann inversion to determine the equilibrium bonded distances and angles from 100 ns AA MD simulations in implicit solvent.

Experimental and computational methods

Experimental

Sample preparation

A lyophilized sample of PGG-PTX-PEG220-npRGD16 was provided by Nitto Denko Technical Corporation (Oceanside, CA, USA). The sample was weighed and then diluted to 1mg/ml in 1X HyClone modified DPBS buffer (Cat. No. SH30028.03, ThermoScientific, Fremont, CA). The sample was then sonicated in a 37°C water bath for 15 min and allowed to settle at room temperature (25°C) for 10 min. Finally, the sample was filtered using a 0.2 μm filter paper (Part No. 431215, Corning Life Sciences, San Francisco, CA, USA) and a 20G 1½ Precision Glide needle (Becton Dickinson, Franklin Lakes, NJ, USA).

Circular dichroism spectroscopy

CD spectroscopy measurements of the PGG-PTX-PEG220-npRGD16 sample were carried out using an AVIV Model 202 spectrophotometer (AVIV Biomedical, Lakewood, NJ, USA). CD spectra of pure 1X DPBS buffer and FD protein in 50% TFE/50% double distilled H2O/0.1% TFA were also taken as the negative and positive controls, respectively. All measurements were taken at 37°C. The far-UV CD spectra were recorded from 190 to 260 nm using a 1-cm rectangular quartz cuvette, and the CD spectra were collected at every 0.5 nm with 5 sec at each point. For each sample except pure 1X DPBS, the concentration was adjusted and diluted so that the dynode voltage corresponding to the ellipticity signal remained below 500 V throughout the entire spectra collection period. (It has been suggested in Greenfield et al. [24] that the signal-to-noise ratio diminishes greatly once the dynode voltage exceeds 500 V.) The optimum concentrations of each sample leading to the optimal elliptical signal and lowest signal-to-noise ratio were determined to be: 1 mg/ml of PGG-PTX-PEG220-npRGD16 in 1X DPBS and 0.25 mg/ml of FD protein in 50% TFE/50% ddH2O/0.1% TFA.

Computational methods

All-atom modeling

Model construction

The initial structures and input files for the following residues were constructed in the xleap module of AMBER 9.0: GG, GGCOO–, GGNH3+, GG-PTX, GG-PEG500-npRGD, GG-PEG1000-npRGD, and GG-PEG2000-npRGD. The initial structure for paclitaxel (PDB ID: 1JFF) was provided by the Protein Data Bank [25]. The Gaussian (g03) program’s geometry optimization and restrained electrostatic potential fitting (RESP) were used to generate the atom-centered point charges for all residues [26, 27]. For the protonated molecules, hydrogen atoms were added to the carboxyl groups to achieve their ionized states. To obtain the electrostatic potential for these residues, the AM1 geometry scheme and the HF/6-31G* and HF/6-31G** pop = mk iop(6/33 = 2) ab initio level calculations were performed using the Gaussian program. Finally, to derive the equivalent partial atomic charges, RESP fitting was applied on the electrostatic potentials.

Energy minimization and MD simulation

These steps were carried out using the AMBER 9.0 sander module. Each molecule was solvated and minimized in implicit solvent using the modified Generalized-Born model of IGB = 2 and the linear combinations of pairwise overlaps (LPCO) model [28, 29]; 250 steps of steepest descent followed by 1,750 steps of conjugate gradient were carried out. The ionic strength of the implicit solvent was set to 140 mM to mimic the salt concentration of blood plasma. No periodic boundary conditions were applied. A non-bonded electrostatic cutoff of 16.0 Å was used. Trajectory snapshots were saved every 100 steps for further reprocessing.

MD simulations were also carried out using the sander module applying the modified version (ff99SB) of the Cornell et al. parm99 and GAFF force fields in the NVT ensemble (in which the number of atoms N, volume V, and temperature T were fixed) at 310K in implicit solvent without periodic boundary conditions [30–32]. Constant temperature scaling was also applied with a time constant of 0.5 ps. Langevin dynamics was used with a collision frequency of 2.0 ps−1. Newton’s equations of motions were integrated with a time step of 2 fs. All bonds involving hydrogen atoms were constrained using the SHAKE algorithm [33]. The rotational and translational degrees of freedom about the center of mass were eliminated. Trajectory snapshots were saved every 100 steps for later reprocessing. Each cycle of the MD simulation was run for 0.1 ns; this step was repeated until a statistical equilibrium was reached at 100 ns.

Coarse-grained modeling

CG mapping

The MARTINI force field was selected for the CG parameterization for its successful application to proteins by experimental validation of their structural and thermodynamic properties as well as the amino acid-based nature of the glutamyl-glutamate backbone [23, 34]. The MARTINI force field dictates that a group of roughly four–five atoms are represented as an interaction center, or bead. The only exception applies to aromatic groups, in which a mapping scheme of two heavy atoms-to-one-bead was applied in order to maintain geometric symmetry, namely, the benzene entities in paclitaxel and nonpeptide RGD. These beads interact through a set of short-range Lennard-Jones potentials to reproduce characteristic properties resulting from AA simulations. Charged groups interact via a Coulombic energy function, and bonded potentials are used to describe the chemical connectivity of the beads. CG mapping of the GG, GGCOO–, GGNH3+, and GG-PTX residues was previously done in Peng et al. [22]. CG mapping of the GG-PEG500-npRGD residue is shown in Fig. 3, and CG mapping of the GG-PEG1000-npRGD and GG-PEG2000-npRGD residues are provided in Fig. S2. Also, each Na + ion is represented as a single Q type bead, and four AA water molecules are mapped to a single P4 type bead (W).

Fig. 3
figure 3

Coarse-grained representation of a GG-PEG500-npRGD residue. A GG-PEG500-npRGD molecule was reduced from 75 heavy atoms to 21 CG beads. (CG representations for the GG-PEG1000-npRGD and GG-PEG2000-npRGD residues are shown in Figs. S1 and S2, respectively)

Full size image

Theory

The dynamics of the beads are described by Newton’s equations of motion. The effective bonded and nonbonded interactions are described by a potential of the form, adapted from Marrink et al. [34]:

$$ V = \sum\limits_m {\,\left[ {V_{{bond}}^m + V_{{angle}}^m + V_{{dihedral}}^m} \right]} + V_{{nonbonded}}^m $$
(2)

where the index m represents a bead from a PGG-PTX-PEG-npRGD molecule, a sodium (Na+) ion, or a water (W) bead. V bond represents the forces between two successively bonded beads, V bond accounts for the forces used to sustain angles between three successively bonded beads, and V dihedral corresponds to the dihedral angle potential for four successively bonded beads. The V nonbonded term accounts for the nonbonded interactions among all of the beads in a system. The bonded interactions are used to describe the chemical connectivity of the beads and are characterized as one of the following (the indices i and j represent two consecutive beads):

$$ V_{{bond}}^m = \sum\limits_{{ij}} {\tfrac{1}{2}k_{{ij}}^m\,{{\left( {r_{{ij}}^m - r_{{0,ij}}^m} \right)}^2}} $$
(3)

where \( r_{{ij}}^m \)is the distance between two consecutive beads i and j, \( r_{{0,ij}}^m \) is the equilibrium bond length, and \( k_{{ij}}^m \) is the distance force constant. The potential for bond angles is:

$$ V_{{angle}}^m = \sum\limits_{{ijk}} {\tfrac{1}{2}k_{{ijk}}^m{{\left[ {\cos \left( {\theta_{{ijk}}^m} \right) - \cos \left( {\theta_{{0,ijk}}^m} \right)} \right]}^2}} $$
(4)

where \( \theta_{{ijk}}^m \) is the angle between beads i, j, and k, \( \theta_{{0,ijk}}^m \) is the equilibrium angle, and \( k_{{ijk}}^m \) is the angle force constant.

The dihedral potential, V dihedral, applies only to backbone protein residues. Since the general structure of PGG-PTX-PEG-npRGD adapts a random coil, V dihedral was set to 0 for all combinations of four successively connected beads:

$$ V_{{dihedral}}^m = \sum\limits_{{ijkl}} {k_{{ijkl}}^m\left[ {1 + \cos \left( {n\chi_{{ijkl}}^m - \delta_{{ijkl}}^m} \right)} \right]} $$
(5)

where n is an integer that represents multiplicity, \( \delta_{{ijkl}}^m \) is the phase shift, \( \chi_{{ijkl}}^m \) is the angle between the plane formed by the first three beads i, j, and k and the plane formed by the last three beads j, k, and l, and \( k_{{ijkl}}^m \)is the dihedral force constant.

Nonbonded interactions between two beads are described by a van der Waals (vdW) potential, characterized by a Lennard-Jones 12-6 potential energy function and a Coulombic potential for charged beads. The following summation is performed over all pairs of nonbonded beads m and n separated by a distance of r mn:

$$ {V_{{nonbonded}}} = \sum\limits_{{m,n}} {4{\varepsilon_{{mn}}}\left[ {{{\left( {\frac{{{\sigma_{{mn}}}}}{{{r_{{mn}}}}}} \right)}^{{12}}} - {{\left( {\frac{{{\sigma_{{mn}}}}}{{{r_{{mn}}}}}} \right)}^6}} \right]} + \sum\limits_{{m,n}} {\frac{{{q_m}{q_n}}}{{4\pi {\varepsilon_0}\varepsilon {r_{{mn}}}}}} $$
(6)

where constants σ mn and ε mn are the vdW parameters for the bead interaction, σ mn refers to the closest distance of approach between two beads m and n, and ε mn denotes the strength of the interactions between beads m and n. The charges q m and q n of the mth and nth beads, respectively, and the vacuum dielectric permittivity ε 0 are included. A relative dielectric constant of ε = 15 is assumed for all electrostatic interactions.

CG parameterization

The bonded distances and angles for the PGG-PTX-PEG-npRGD molecules were determined as follows: (1) MD simulations of the 12 AA PGG-PTX-PEG-npRGD models were run to 100 ns; (2) the bonded distances and values for the assigned beads were extracted from the AA MD simulations using the AMBER 9.0 ptraj module; and (3) the Boltzmann inversion was used to extract the equilibrium bonded distances and angles [35]. These values were averaged for each set of four PGG-PTX-PEG-npRGD molecules with the same PEG MW. (For example, the bonded distances and angles for the four PEG1000 molecules were averaged together; the resulting values were used for the parameterizing the PGG-PTX-PEG1000-npRGD4, PGG-PTX-PEG1000-npRGD8, PGG-PTX-PEG1000-npRGD12, and PGG-PTX-PEG1000-npRGD16 molecules.) CG topologies for the GG, GGCOO–, GGNH3+, and GG-PTX residues were taken from Peng et al. and shown in Table S1 [22]; CG topologies for the GG-PEG500-npRGD residue are provided in Table S2; and CG topologies of the GG-PEG1000-npRGD and GG-PEG2000-npRGD residues are provided in Table S3.

MD simulation

CG models of the 12 PGG-PTX-PEG-npRGD molecules were constructed using the atom2cg and itp-generating scripts provided by the MARTINI tutorials. Each PGG-PTX-PEG-npRGD molecule was solvated in explicit solvent comprising of W beads. Due to the negative charges imparted by the glutamyl-glutamate residues, 234 Na + ions were added in place of water molecules to neutralize the system. The simulations were carried out under NPT (the number of particles N, pressure P, and temperature T are fixed) conditions at a constant temperature of 310 K with a coupling constant of τT = 0.1 ps and a constant pressure of 1 bar with a relaxation time of τP = 0.5 ps. The cutoff length used for the nonbonded interactions was r cut = 1.2 nm. Lennard-Jones and Coulombic forces were considered for r cut   < 0.9 nm and r cut  < 1.2 nm, respectively. Coulombic forces were computed every time step for 1.0 nm and once every 10 time steps for 0.9 nm < r cut < 1.2 nm. The time step in the leap-frog integration scheme was 5 fs. The energies, coordinates and velocities were written every 0.5 ps. MD simulations were run for 1 μs. The timescale was scaled up by a factor of four because CG dynamics are faster than all-atom dynamics, due to the smoothness of the CG interactions. This factor of four is based on the speedup in the diffusional dynamics of CG water as compared to AA water. Therefore, the 1 μs simulation timescale was scaled up by a factor of four, resulting in an effective time of 4 μs [23, 34, 36].

Root-mean-square deviation clustering

Root-mean-square deviation (RMSD) clustering—a quantitative method used to analyze protein structures [37]—was used to determine the PGG-PTX-PEG-npRGD conformations. All molecular conformations accessed throughout an MD simulation were categorized into clusters based on structural similarity. For each cluster, a structure from the MD trajectory closest to the conformational average of all structures in the MD trajectory was regarded as the representative conformation. Using the g_cluster module of the GROMACS 4.0.3 package, the RMSD cutoff length was optimized manually so as to meet the following criteria: (1) there were approximately 40 clusters total, (2) very few clusters contained just one representative conformation, and (3) 90% of the MD trajectory was in clusters fewer than those with one representative conformation. Table 1 provides the optimized RMSD cutoff distances for all 12 PGG-PTX-PEG-npRGD molecules. (Sample optimization of the RMSD cutoff length vs. number of clusters generated for the GG-PTX-PEG1000-npRGD12 molecule is shown in Fig. S2.) After optimization of the RMSD cutoff lengths, the significance of every representative conformation was assessed such that its corresponding frames occupy at least 10% of the entire 4 μs trajectory. Once the representative conformations were judged to be sufficiently significant, each was evaluated for its uniqueness as follows: if the RMSD between the significant representative conformations of each significant cluster were greater than twice the RMSD cutoff length, then those significant representative conformation(s) was/were regarded as unique. (For the PGG-PTX-PEG-npRGD molecules with only one significant representative conformation, that conformation was automatically deemed unique.) To verify the uniqueness among clusters, the RMSD Calculator tool in Visual Molecular Dynamics 1.8.6 was used [38].

Table 1 Summary of data from root-mean-square deviation (RMSD) clustering method. The RMSD cutoff lengths optimized for each PGG-PTX-PEG-npRGD molecule and the corresponding number of clusters are listed, as well as the unique and significant central members and the corresponding RMSD between each central member
Full size table

Density analysis

Typically in MD simulations, the radial distribution function, g(r), is used to characterize the structure of the biomolecules. The g AB(r) between species A and B is defined as follows [39]:

$$ {g_{{AB}}}(r) = \frac{{{\rho_{{AB}}}(r)}}{{{\rho_{{B,local}}}}} $$
(7)

where ρ AB(r) is the density of species B at a distance r from species A and ρ B,local is the local density of species B at maximum distance r max from species A. For this study, species B refers to a PGG-PTX-PEG-npRGD molecule and species A its center-of-mass (COM). It is important to emphasize that the COM does not necessarily correspond to the geometric center, particularly for irregularly shaped molecules, and the COM is always towards the more massive end.

Since the local densities of the PGG-PTX-PEG-npRGD molecules are different, their g AB(r) values would not properly describe how the structure of PGG-PTX-PEG-npRGD is affected by PEG MW and npRGD density. Therefore, the density of PGG-PTX-PEG-npRGD from its COM, ρ AB(r), was used to describe structure. The local density of each PGG-PTX-PEG-npRGD molecule, ρ B,local, was determined from its molecular mass and simulation box size (see Table 2). The g_rdf module of GROMACS 4.0.3 was then used to extract the g AB(r) values of each PGG-PTX-PEG-npRGD molecule throughout the 4 μs trajectory. Using the definition of g AB(r), the density of PGG-PTX-PEG-npRGD from its spatial center was determined by multiplying g AB(r) and ρ B,local: ρ AB(r) = g AB(r) × ρ B,local.

Table 2 Determination of local densities of PGG-PTX-PEG-npRGD molecules. The local density, ρ B,local, of each PGG-PTX-PEG-npRGD molecule was determined by dividing the mass of each molecule by the simulation box volume using the following unit conversions: 6.022 × 1023 Da = 1 gram and 1 nm3 = 1.0 × 10−21 ml
Full size table

Solvent-accessible surface area analysis

The solvent-accessible surface area (SASA) indicates the area of a particular species that is exposed to the solvent. The SASA of the PGG-PTX-PEG-npRGD molecules was analyzed in VMD 1.8.6 using a solvent probe radius of 2.4 nm. For each PGG-PTX-PEG-npRGD molecule, the SASA of the entire molecule was calculated, along with the SASA of PGG, PEG, PTX, and npRGD entities of each PGG-PTX-PEG-npRGD molecule. The contribution of each entity was then determined by taking the fraction of the SASA for a particular entity with respect to the SASA of the whole PGG-PTX-PEG-npRGD molecule, or the percentage (%) SASA. Table 3 shows the % SASA for the PGG, PEG, PTX, and npRGD entities for all PGG-PTX-PEG-npRGD molecules. The hydrophilic SASA (% SASAphil) was determined from the sum of the % SASAPGG and % SASAPEG residues.

Table 3 Percentage solvent-accessible surface area (SASA) of PGG-PTX-PEG-npRGD molecules. Lists the percentage SASA values for the PGG, PEG, PTX, and npRGD entities with respect to the entire SASA of each PGG-PTX-PEG-npRGD molecule
Full size table

Results and discussion

Circular dichroism spectroscopy

Circular dichroism spectroscopy was used to determine the secondary structure of a PGG-PTX-PEG220-npRGD16 molecule. This method of determining the initial structure was originally applied to PGG-PTX in Peng et al. [22]. Although PGG-PTX-PEG220-npRGD16 was not one of the PGG-PTX-PEG-npRGD molecules pursued in this study, its chemical structure is similar enough to provide an accurate representation of the PGG-PTX-PEG-npRGD molecules being studied.

Figure 4 shows the CD spectra of 1mg/ml PGG-PTX-PEG220-npRGD16 dissolved in 1X DPBS. The curve is most representative of a random coil at physiological temperature (37°C). The negative control of 1X DPBS at 37°C indicates that there is no signal interference that may have influenced the conditions while the CD spectrum of PGG-PTX-PEG220-npRGD16 was being taken. For the positive control, the shape of the alpha-helical FD protein, a bZIP transcription factor in a floral pathway [40], in 50% TFE/50% ddH2O/0.1% TFA buffer does indeed correspond to an alpha-helix, as compared with the CD spectrum of the alpha-helical poly(γ-tyrosine) [41]. The lowest part of the curve is ∼215 nm and the highest point of the curve is ∼200 nm, collectively indicating that the FD protein is alpha-helical. Since PGG-PTX-PEG220-npRGD16 exists as a random coil, this data suggest that the PGG-PTX-PEG-npRGD molecules being modeled most likely exist as random coils. Therefore, their existence as random coils imposes no stipulations on the construction of AA models of PGG-PTX-PEG-npRGD. Also, the CD spectroscopy results further support the theory that the conformation of PGG-PTX-PEG-npRGD is not limited to a specific structure, and structural analysis must account for multiple conformations of PGG-PTX-PEG-npRGD.

Fig. 4
figure 4

Circular dichroism spectra of samples. The spectra for 16 nonpeptide RGD-targeted PGG-PTX (red), PGG-PTX (blue), FD protein (green), and pure 1X DPBS (black) are shown

Full size image

RMSD values

Figure 5 shows the RMSD time evolutions for all PGG-PTX-PEG-npRGD molecules. Most time evolutions increased up to ∼6-7 nm and reached a statistical equilibrium by 4 μs. Given the relatively large sizes of these systems, these values are reasonable. The RMSD values for the PEG2000 molecules are slightly higher than those of the PEG500 and PEG1000 molecules (by ∼1 nm), whereas the RMSD values for the PEG500 and PEG1000 molecules are comparable. This 1 nm difference is in accordance with the size of the molecular systems, as the PEG2000 molecules are the largest. In addition, the time it takes for PGG-PTX-PEG-npRGD molecules to attain statistical equilibrium depends on the PEG MW. The PEG1000 molecules take the longest time to reach equilibrium, as some molecules do not reached a plateau until 3 μs. Conversely, the PEG2000 molecules take the least time to reach equilibrium; by 1 μs, all molecules have attained a statistical equilibrium. The time it takes for the PEG500 molecules to level off is ∼2 μs, which is intermediate compared to the time it takes for the PEG1000 and PEG2000 molecules to reach equilibrium. This behavior suggests that the PEG1000 molecules fluctuate the most, accessing a wide range of molecular conformations throughout the 4 μs MD trajectory.

Fig. 5
figure 5

RMSD time evolutions of PGG-PTX-PEG-npRGD molecules. Data for molecules with 4 (PEG) n -nonpeptide RGD molecules (blue), 8 (PEG) n -nonpeptide RGD molecules (green), 12 (PEG) n -nonpeptide RGD molecules (red), and 16 (PEG) n -nonpeptide RGD molecules (black) are shown

Full size image

Size of a PGG-PTX-PEG-npRGD molecule

Figure 6 shows the radius of gyration (R gyr) trajectories for all PGG-PTX-PEG-npRGD molecules over the 4 μs time scale. For all molecules, the R gyr values begin at ∼7–8 nm and decrease to ∼3 nm, nearly half their original size. The PEG2000 and PEG1000 systems take the shortest and longest to reach their final R gyr values, respectively. These values correspond to the RMSD values: the PEG2000 and PEG1000 systems also take the shortest and longest to reach their final RMSD values, respectively. Furthermore, as shown in Fig. 7, the R gyr of each PGG-PTX-PEG-npRGD representative conformation is approximately 2-4 nm. Figure 8 further supports the notion that neither the PEG MW nor the npRGD density significantly influences the size of a PGG-PTX-PEG-npRGD molecule.

Fig. 6
figure 6

Radius of gyration trajectories of PGG-PTX-PEG-npRGD PTX molecules. Data for molecules with 4 (PEG)n-nonpeptide RGD molecules (blue), 8 (PEG)n-nonpeptide RGD molecules (green), 12 (PEG)n-nonpeptide RGD molecules (red), and 16 (PEG)n-nonpeptide RGD (black) molecules are shown

Full size image
Fig. 7
figure 7

Representative conformations of PGG-PTX-PEG-npRGD molecules. Each structure shows GG (grey), PTX (yellow), PEG (blue), and npRGD (red) residues. Explicit W beads and Na + ions are not shown. Representing an actual frame in the trajectory, the central member embodies the configuration that is the most similar to the average of all configurations in that cluster. The percentage indicates the population of frames, or % trajectory occupancy, corresponding to that particular cluster for the full 4 μs MD trajectory. Radius of gyration values for each PGG-PTX-PEG-npRGD molecule is also shown below the % trajectory occupancy

Full size image

Structure of a PGG-PTX-PEG-npRGD molecule

Figure 9 shows that each PGG-PTX-PEG-npRGD molecule is most densely packed at its COM. Most PGG-PTX-PEG-npRGD molecules exhibit a peak density of ∼1 g/ml just outside its center-of-mass at r ∼ 1 nm, 1.5 nm, and 2 nm for the PEG500, PEG1000, and PEG2000 molecules. This increase in distance of peak density is correlated with increasing PEG MW. Figure 8 shows how the density of each residue type (PGG, PTX, PEG, and npRGD) varies from the COM of a PGG-PTX-PEG-npRGD molecule. The PEG MW impacts how the density of each residue changes with distance from COM. For instance, PTX and npRGD residues exhibit peak densities at greater distances as PEG MW increases, suggesting that higher PEG MW causes PTX and npRGD residues to approach closer to the surface of the PGG-PTX-PEG-npRGD. The peak densities of PTX and npRGD in Fig. 8 agree with the peak densities of PGG-PTX-PEG-npRGD in Fig. 9, suggesting that the presence of PTX and npRGD residues is mostly responsible for this spatial shift in the peak density of the overall PGG-PTX-PEG-npRGD molecule. On the contrary, for PEG, the distance of peak density tends to decrease as the PEG MW increases, which suggests that increasing the PEG spacer length causes PEG to shift away from the external surface and move internally towards the COM of a PGG-PTX-PEG-npRGD molecule. In general, the difference in peak densities of PEG and peak densities of PTX and npRGD becomes more pronounced as the PEG MW increases, which signifies that the PGG-PTX-PEG-npRGD molecule tends to become more spatially divided (by residue type) at higher PEG MW. Finally, taking the R gyr values into account, these results suggest that while the PEG MW does not affect the size of a PGG-PTX-PEG-npRGD molecule, the PEG MW does impact structure. Figure 8 also suggests that the npRGD density also does not significantly affect the structure of the PGG-PTX-PEG-npRGD molecule.

Fig. 8
figure 8

Density of PGG-PTX-PEG-npRGD molecules from COM. Data for molecules with 4 (PEG)n-nonpeptide RGD molecules (blue), 8 (PEG)n-nonpeptide RGD molecules (green), 12 (PEG)n-nonpeptide RGD molecules (red), and 16 (PEG)n-nonpeptide RGD molecules (black) are shown

Full size image

Implications of shape of a PGG-PTX-PEG-npRGD molecule

Discher and colleagues have demonstrated that filamentous particles retains a longer circulation half-life in the bloodstream, thus decreasing levels of clearance via the reticoendothelial system (RES) and improving the chances of the nanoparticle to reach the target tumor [16]. Also, nonspherical nanoparticles have been shown to adhere more firmly to walls of tumor endothelia [17]. Figure 7 shows the representative conformations of all PGG-PTX-PEG-npRGD molecules. The geometry of each molecule is characterized as a filament (wormlike and threadlike) or globule (round). In general, most molecules are globules, which is most likely attributed to the hydrophobic driving force among PTX molecules that are uniformly distributed along the PGG backbone. The uniform distribution of PTX causes the entire PGG-PTX-PEG-npRGD to self-assemble into a globule in which PTX molecules are located in the core, where it is the furthest from the solvent and surrounded by hydrophilic PGG. Other molecules, such as PGG-PTX-PEG1000-npRGD4 and PGG-PTX-PEG1000-npRGD12, exhibit both filamentous and globular behavior. The shape of the PGG-PTX-PEG-npRGD molecule is not significantly affected by the npRGD density but by the PEG MW. The PEG1000 molecules exhibit the most filamentous behavior, as each one of its four molecules is characterized by a filamentous conformation. The higher prominence of filaments in the PEG500 and PEG1000 molecules is possibly a result of the sterical interactions among the rigid, conformationally constrained nonpeptide RGD residues that hinder the formation of a more compact globule. For the PEG2000 molecules, however, the prominence of the globular shapes (∼90%) suggests that 2,000 Da is the PEG MW most likely to prevent sterical interactions among nonpeptide RGD residues, thus allowing the PGG-PTX-PEG-npRGD molecules to form a globule. Overall, the results suggest that PGG-PTX-PEG-npRGD molecules with PEG500 and PEG1000 spacers are more efficacious than the PGG-PTX-PEG-npRGD molecules with PEG2000 spacers. While the molecular weight of the PEG spacer influences the shape of a PGG-PTX-PEG-npRGD molecule, the npRGD density has only minimal impact.

Fig. 9
figure 9

Density of PGG-PTX-PEG-npRGD molecules from COM. Data for molecules with 4 (PEG)n-nonpeptide RGD molecules (blue), 8 (PEG)n-nonpeptide RGD molecules (green), 12 (PEG)n-nonpeptide RGD molecules (red), and 16 (PEG)n-nonpeptide RGD molecules (black) are shown

Full size image

Implications of flexibility of a PGG-PTX-PEG-npRGD molecule

Takeoka has speculated that higher molecular flexibility of nanoparticles promotes more surface interactions to tumor cells [18]. The flexibility of each PGG-PTX-PEG-npRGD molecule was assessed by the number of representative conformations of each molecule and, for those systems with multiple representative conformations, the similarity among representative conformations. A higher number of representative conformations indicates a stronger tendency for that molecule to fluctuate and change shape. The PEG1000 molecules have the most representative conformations; the PGG-PTX-PEG1000-npRGD8 and PGG-PTX-PEG1000-npRGD12 molecules are characterized by three representative conformations, although with low significance (∼10–20%). On the contrary, each PEG2000 molecule has only one unique shape with very high significance (∼90%). Moreover, the degree of similarity among the representative conformations is also an indicator of flexibility: a PGG-PTX-PEG-npRGD molecule with two globular representative conformations is less flexible than a molecule characterized by globular and filamentous conformations. For instance, PGG-PTX-PEG500-npRGD4 and PGG-PTX-PEG1000-npRGD4 each has two representative conformations, but a larger structural difference exists between the globular and filamentous shapes of the PGG-PTX-PEG500-npRGD4 molecule, as compared to the two filamentous conformations of the PGG-PTX-PEG1000-npRGD4 molecule. In all, these observations show that the PEG2000 molecules are the least flexible, whereas the PEG500 and PEG1000 molecules are more flexible and would therefore bind more effectively to tumor cells. Finally, the npRGD ligand density does not have a significant effect on the flexibility of a PGG-PTX-PEG-npRGD molecule.

Implications of surface hydrophilicity of a PGG-PTX-PEG-npRGD molecule

Nanoparticles with higher degrees of surface hydrophilicity have been shown to possess longer circulation half-lives in the bloodstream, suggesting a higher preference for accumulating in leaky vasculature of tumors [14]. The degree of surface hydrophilicity of each PGG-PTX-PEG-npRGD molecule was assessed by the % SASA of the hydrophilic entities: PGG and PEG (see Table 3). These % SASA values are related to the density plots in Fig. 9. In general, the high % SASAPGG values can be attributed to the high density values at r ∼ 3 nm, at the external surface of the PGG-PTX-PEG-npRGD molecule. Likewise, the % SASAPEG values are the lowest, which is a result of the high concentration of PEG near the center-of-mass and not the surface. Figure 8 also shows that the densities of PTX and npRGD residues are relatively high near the surface of a PGG-PTX-PEG-npRGD molecule, but their low molecular weight ratio with respect to the entire molecule results in lower % SASAPTX and % SASAnpRGD values than % SASAPGG. As expected, the surface hydrophilicity increases with PEG MW, and an increase in npRGD density results in higher % SASAnpRGD. A higher npRGD density also tends to cause a decrease in the % SASAphil. While nonpeptide RGD is used to promote active targeting of PGG-PTX to tumor cells, too many nonpeptide RGD ligands may decrease surface hydrophilicity, which may unfavorably decrease to the circulation half-life. PGG-PTX-PEG-npRGD must remain in the bloodstream long enough to diffuse convectively into leaky tumor vasculature. At the same time, adequate amounts of nonpeptide RGD ligands must be exposed to the solvent to achieve effective binding to αVβ3 integrins. Therefore, optimization of npRGD density is necessary to attain optimal efficacy of PGG-PTX-PEG-npRGD. Taking these factors into consideration, the PGG-PTX-PEG-npRGD molecules with 8 and 12 npRGD and PEG1000 and PEG2000 MW are most likely to confer the longest circulation half-lives and effective targeting to αVβ3 integrins on tumor cells.

Recommendations

This study was founded on the hypothesis that the degree of PEGylation (500 Da; 1,000 Da; and 2,000 Da) and npRGD density (4, 8, 12, and 16 nonpeptide RGD per PGG-PTX-PEG-npRGD molecule) collectively affect the size, structure, shape, flexibility, and surface hydrophilicity of PGG-PTX-PEG-npRGD. While the size of a PGG-PTX-PEG-npRGD molecule is consistent at all PEG MW and npRGD densities, the molecular density tends to be is higher at the surface at higher PEG MW. The results also show that the PEG500 and PEG1000 molecules exhibited more filamentous behavior and higher flexibility, whereas the PEG2000 molecules are globular and more rigid. An increase in the nonpeptide RGD ligand density may result in lower surface hydrophilicity, yet an insufficient number of nonpeptide RGD ligands may prevent PGG-PTX-PEG-npRGD from achieving effective targeting to tumor cells. All factors considered, the ideal PGG-PTX-PEG-npRGD candidate would possess the following properties: filamentous shape, high flexibility, and high surface hydrophilicity. Based on the results, the PGG-PTX-PEG1000-npRGD4 and PGG-PTX-PEG1000-npRGD8 molecules most satisfy these criteria and are the most promising candidates.

Conclusions

The roles of shape, flexibility, and surface hydrophilicity in influencing the delivery of therapeutic nanoparticles to tumors have attracted increasing attention. These factors can serve as a basis for the design of anticancer therapeutics, such as PGG-PTX-PEG-npRGD. While it may be technically feasible to determine these properties using experimental methods, the effort required to synthesize all candidate compounds and subsequently perform experimental testing would be very time- and resource-intensive, especially since the testing is usually conducted in a trial-and-error manner. Therefore, MD simulations were used to explore the size, shape, flexibility, and surface hydrophilicity of candidate PGG-PTX-PEG-npRGD compounds in a more efficient and controlled manner. The theoretical results can be used to further guide experiments, potentially shortening the preclinical development of PGG-PTX-PEG-npRGD. This study also introduces a new structural analysis method of characterizing semiflexible molecules without a specific, crystallized conformation.