Document Type : Reasearch Paper
Authors
Department of Chemistry, Amirkabir University of Technology (Tehran Polytechnic), P.O.Box:15875-4413, Tehran, Iran.
Abstract
Keywords
INTRODUCTION
Although recently substantial advances have been achieved in developing anticancer drugs, cancer still is one of the leading causes of death worldwide [1, 2]. Factors contributing to ineffectiveness of cancer therapeutics and their drawbacks include lack of tumor specificity, cancer cell heterogeneity, drug resistance and insufficient drug accumulation in tumor tissue [3-5]. Drug resistance can be caused by drug inactivation, change in drug targets and drug efflux [6-8]. Also, some drugs poorly differentiate between normal and cancer cells and this is the main shortcoming of such anticancer drugs due to it can lead to incomplete tumor response, unwanted side effects and ultimately therapeutic failure [9-11]. Thus, it was attempted to develop site specific drug delivery systems (DDSs) and effective drug targeting methods in order to overwhelm the systemic toxicity triggered by cancer drugs [12-17]. In numerous drug delivery platforms, several strategies are commonly used to deliver drugs such as self-assembly, modification by polyethylene glycol (PEG) or PEGylation, stimulus sensitive systems (including pH- and thermo-sensitivity), improved permeability and retention as well as using cell-penetrating moieties and prodrugs [18-24].
Proficient cell penetration is significant for crossing the plasma membrane, decreasing side effects and administering the dose essential for activity [25-29]. Peptides are utilized as cargo carriers in DDSs due to their capability to heighten the cellular uptake of drugs and thereby afford an effectual therapeutic role. The interactions of peptides as well as interactions of drug and carrier have frequently been investigated [8, 30-34]. Most of peptides acting as potential cargo carriers are known as cell penetrating peptides (CPPs) which are also identified as protein transduction domains [35]. CPPs were first discovered in 1988 and have broadly been exploited as drug carriers as a result of their cell penetrability [36].
The CPPs are one of the most prevalent and effective compounds used to accomplish cellular uptake. They are a class of various peptides, usually having 5–30 amino acids, and contrasting most peptides, they are able to cross the cellular membrane [25]. The idea of protein transduction to cells was offered in 1988 by Frankel and Green in parallel. These well-organized transport systems have lower cytotoxicity to diverse cells in comparison to other delivery devices. Hence, such pharmaceutical carriers are employed to enhance efficacy and stability and support intracellular transport but decrease undesirable drugs’ side effects [37].
Zhang and coworkers reported that the RGD peptide could be used for targeted delivery of epirubicin. They showed that RGD-conjugated and EPB-loaded inulin/ibuprofen nanostructures could increase tumor accumulation, thus decrease the systemic toxicity related to EPB in addition to simultaneously improving the effectiveness of the drug delivery [38]. It was indicated that the DDS formulated with RGD-modified, redox-sensitive and prodrug based lipid-polymer nanoparticles was an efficient system for the co-delivery of both cisplatin and paclitaxel for the lung cancer therapy [39].
Literature review indicates that few computational approaches have so far been reported using RGD as the DDS [25]. However, the interactions of RGD with various nanoparticles were examined using the DFT computations. For instance, the interactions of RGD with carbon nanomaterials, comprising defective and pristine graphene sheets as well as graphene oxide were investigated by the DFT calculations and the adsorption energy values indicated that the graphene interacted with all of the RGD functional groups so that the most powerful adsorption was occurred while the RGD was parallel relative to the graphene surface [40]. In another study, the RGD peptide interaction with the rutile surface was examined by DFT computations [41]. It was shown the peptide was bound only through the carboxyl groups of aspartic acid and the bonding had a sigma character. Furthermore, the RGD adsorption was very similar to that of simple carboxylic acids like formic acid adsorption in both geometry and electronics properties [42].
Cyclophosphamide (CP) is an anticancer drug and alkylating agent which is frequently employed in autoimmune disorders [43]. It belongs to the phosphoramide family [44-46]. The mechanism of CP prompted cell damage or killing of all tumor and normal cells is not well defined but it has been indicated that CP is able to generate ROS (reactive oxygen species). All of activated cyclophosphamide metabolites are carried through the blood circulation system to both tumor and healthy organs/tissues to kill the DNA and proteins in cells plus lysosomal and mitochondrial membranes. Additionally, CP is a familiar immunosuppressant in multiple sclerosis, organ transplantation, systemic lupus erythematosus and rheumatoid arthritis [12, 47, 48].
Computational studies are frequently accomplished with the intention of predicting the electronic and structural features of various kinds of systems and matrixes [49-51]. Herein, subsequent to our earlier computational works on DDSs [52-55], the DFT computations were done in ethanol solution using the B3LYP/6-31+G(d,p) method on twenty drug delivery systems composed of the RGD as a promising carrier and the CP anticancer drug. The binding energies, thermodynamics properties, dipole moments, molecular electrostatic potentials, molecular descriptors, density of states spectra and electronic properties of the systems were acquired and evaluated. Besides, properties obtained from the QTAIM data were employed to determine the nature of interactions. Finally, the MD simulations were also conducted in ethanol solvent with the purpose of finding the effects of polyethylene glycol (PEG) polymeric chains added to the CP-RGD system on the drug delivery efficacy.
EXPERIMENTAL
The structures of the CP drug, RGD carrier in addition to their twenty hydrogen bonded structures including CP-RGD-1 to CP-RGD-20 were completely optimized by the DFT calculations using the B3LYP method and 6-31+G(d, p) basis set by means of Gaussian 98 program [56]. The SCRF keyword of the Tomasi’s polarized continuum (PCM) model was used to accomplish the calculations in ethanol solution. The analysis of frontier molecular orbitals were done to obtain the HOMO (Highest Occupied Molecular Orbital), LUMO (Lowest Unoccupied Molecular Orbital) orbitals, band gaps, quantum molecular descriptors, contour maps, electrostatic surface potentials (ESPs), density of states (DOS) as well as partial density of states (PDOS) spectra. The topological analysis of the electron densities for the optimized geometries was conducted by the QTAIM computations [44]. For this purpose, the wave functions of the electron densities were achieved for each structure and applied as input files into the AIM2000 program to calculate the topological features of the electron densities and the Laplacian values. According to the local Virial theory, properties of bond critical points, BCPs, were used to evaluate the nature of interactions such as hydrogen bonds, electrostatic interactions and covalent bonds.
The MD simulations were accomplished using the Materials Studio software [57]. The condensed phase optimized molecular potentials (COMPASS) force field was used in all minimizations and dynamics simulations. Such force field is commonly applied for all kinds of structures in atomistic MD simulations. The COMPASS force field is based upon the Polymer Consistent Force Field (PCFF) that is the first force field parameterized and approved by the condensed phase principles, empirical and ab initio data obtained for different molecules. Along with the non-bonding and bonding interactions utilized in other force fields, also the cross-coupling contacts are measured by the COMPASS. Consequently, it accurately estimates numerous properties of various materials (mainly polymers).
The amorphous cells were generated containing 10 RGD, 3 CP drug and 100 C2H5OH molecules in addition to 0-10 polyethylene glycol (PEG) polymeric chains (containing 20 repeating units) with the purpose of determining the most appropriate PEG amount for the most controlled delivery of the CP anticancer drug. The amorphous cells were initially generated using amorphous cell module with a low density (0.3 g/cm3) to allow them to afford correct equilibrium density values after the MD simulations. All of the configurations used in the amorphous cells were minimized by means of the smart minimizer algorithm combining Newton, steepest descent and conjugate gradient algorithms. The minimizations were done in a cascading method for 4,000,000 iterations till the systems reached the relaxation. Subsequently, the dynamic NVT (constant number of molecules, constant volume and constant temperature) simulations were carried out for 25 ns at 298.15 K and 1 atm. Next, the NPT (constant number of molecules, constant pressure and constant temperature) simulations were accomplished for 25 ns (at 1 atm and 298.15 K). The non-bond interactions such as the van der Waals and the electrostatic potentials were evaluated by the Ewald method. Finally, the NVT simulations were run on the relaxed cells for 5 ns. The last 1 ns trajectories were employed in order to acquire the structural, dynamical as well as energetic properties of all systems.
RESULTS AND DISCUSSION
Optimization, thermodynamics and structural properties
In this work, arginyl-glycyl-aspartic acid (RGD) tripeptide was used as the carrier for the anticancer CP drug which was attached to different RGD active sites to produce twenty DDSs for which their energetic, structural, electronic and spectral properties were investigated with the aim of achieving the most suitable drug vehicle. Table 1 shows the hydrogen bonding data related to the hydrogen bonds generated between the CP drug and the RGD. In this table, the D-H…A generally indicates the hydrogen bonds formed as the Donor-H…Acceptor. Structures 1 to 8 are produced through the hydrogen bonds formed between phosphoryl (P=O) oxygen atom of the CP drug and H atoms of RGD carrier (O-H and N-H functional groups). Compounds 9-12 have been produced through hydrogen bonding formation between the endocyclic N-H group of CP and oxygen atoms of carbonyl groups in the RGD. Also, the hydrogen bonds among Cl atoms of the CP and the H atoms of the RGD create structures 13 to 20. Figs. 1(a-c) and 2(a-c) display the atom numbering formula schemes as well as the optimized structures of the CP drug, RGD carrier and the CP-RGD-7 carrier (with the most negative binding energy), respectively.
The binding energies of drug and carrier were measured by formula ΔEbinding=Emolecule-Σi Ei, i=atom and those of structures 1-20 were computed using formula ΔEbinding=EDrug-Carrier-(EDrug+ECarrier)+BSSE+ΔZPVE. These results and dipole moments of all structures are listed in Table 2. According to the data of this Table, most of the CP-RGD compounds are unstable because of their positive binding energies; thus, it is likely that these compounds will not be formed while structures 1, 6, 7, 8 and 10 with negative ΔEbinding values are stable. Also, it is found that the structure CP-RGD-7 is the most stable compound having a ΔEbinding value of -5.22 kcal/mol. The dipole moment values of both the CP and the RGD molecules are about 8 Debye but it increases to 11.17, 10.63, 14.87, 10.15 and 9.41 in compounds 1, 6, 7, 8 and 10, respectively.
Ethanol is an organic solvent that is most usually used in drug delivery but other alcohols such as isopropyl- and propyl- alcohols have seldom been applied [58]. There are many theoretical and experimental studies that used ethanol in drug delivery systems [59]. As ethanol has a relatively high dipole moment (1.66 Debye) and due to its non-toxicity, it can be a suitable solvent for application in the delivery of the CP-RGD polar structures. The CP and RGD reveal polarizability values of 189.63 and 268.17 Bohr3 but the twenty CP-RGD compounds exhibit values in the range of 458.85-466.43 Bohr3 which indicate the polarizability is highly increased upon the H-bonds formation between the RGD and the CP molecules.
The Gibbs free energy (ΔGinteraction), enthalpy (ΔHinteraction) and entropy (ΔSinteraction) values for the CP-RGD interactions in compound 1-20 are given in Table 3. The ΔG, ΔH as well as ΔS values are determined by the formula ΔGinteraction=GCP-RGD–(GCP+GRGD), ΔHinteraction=HCP-RGD–(HCP+HRGD) as well as ΔG=ΔH-TΔS (T=298.15 K), respectively. It is found that the ΔH and ΔG values are obtained within the range of 0.91-1.90 and 8.71-13.39 kcal/mol, respectively. All of the ΔH and ΔG values are positive for all CP-RGD structures which confirm all of the interactions are endothermic and nonspontaneous. As well, the negative ΔS values demonstrate that the entropies of these systems decrease upon the CP-RGD interactions.
Table S1 illustrates selected bond lengths and bond angles for the CP drug, the RGD carrier and structures 1, 6, 7, 8 and 10. For all compounds, the P=O and the P–O bond lengths are in the range of 1.50–1.51 and 1.61–1.62 Å, respectively, that are comparable with those reported for the P=O (1.45 Å) and the P–O (1.64 Å) bond lengths [60, 61]. The P–N bond lengths are about 1.67 Å which are in accordance with those of the literature reported P–N bond lengths [60, 61]. In all compounds, the C=N and the C–N bond lengths change in the range of 1.29-1.30 and 1.38-1.48 Å while the C=O and the C–O bond lengths are in the range of 1.22-1.24 and 1.32-1.35 Å, respectively; these values are close to those reported elsewhere [62, 63]. The H-bond lengths formed among O1 atom of the CP and the RGD hydrogen atoms in structures 1, 6, 7 and 8 are measured to be 1.83, 1.97, 1.56 and 1.59 Å, respectively. Moreover, in the CP-RGD-10 structure, the hydrogen bond length formed between the hydrogen atom of the endocyclic N-H group in the CP and the O6 atom of the carbonyl group in the RGD is 1.98 Å. As all of these hydrogen bond lengths are below 2 Å, it is confirmed that very strong H-bonds are formed between the functional groups of both the CP drug and the RGD carrier. Particularly, the hydrogen bond length in structure 8 is 1.59 Å that is very much small and indicates a very much strong H-bond has been formed.
The bond angles at around a tetrahedral atom in a molecule are about 109.5°. The six surrounding bond angles around the P atom in the CP drug, i.e. the O1–P–O2, O1–P–N1, O1–P–N2, O2–P–N1, O2–P–N2 and N1–P–N2 angles are 110.53, 113.55, 117.78, 107.09, 103.27 and 103.62°, respectively, indicating an average value of 109.31° that is very close to the value expected for a tetrahedral atom. Furthermore, these six bond angles in the CP-RGD-7 are equal to 110.41, 112.35, 115.03, 108.31, 104.22 and 105.99°, respectively, with an average of 109.39° which represents a tetrahedral configuration for the phosphorus atom.
It is known that the bond angles in a trigonal planar geometry are close to 120°. There are three bond angles of C1–N1–H1, C1–N1–P and H1–N1–P around the N1 atom in the CP which are equal to 114.19, 118.10 and 112.31°, respectively, with an average of 114.87°. Additionally, the bond angles surrounding the N1 atom in the CP-RGD-7 structure and those at around the N2 atom of the CP drug and the CP-RGD-7 exhibit average values equal to 114.95, 118.77 and 116.83°, respectively, demonstrating there are some deviations from 120° (trigonal planar angle). The P…O-H bond angles measured for the structures 1, 6, 7 and 8 as well as the N1-H1…O6 bond angle in the CP-RGD-10 structure are equal to 129.99, 121.45, 130.80, 127.90 and 173.06° respectively.
Quantum molecular descriptors
Table 4 affords the quantum molecular descriptors for the CP drug, the RGD carrier as well as the twenty hydrogen bonded systems including CP-RGD-1 to CP-RGD-20. These descriptors can describe and predict the chemical behaviors and the electronic properties of systems. The global hardness (η) as well as the chemical potential (µ) are two universal descriptions for the structural reactivity and stability which can be calculated by the formula η=(I−A)/٢ and µ=−(I+A)/٢, respectively, in which I (–EHOMO) and A (–ELUMO) exhibit the vertical ionization energy plus the vertical electron affinity, respectively. Also, Table 4 provides the band gap energies that are the differences between the LUMO and the HOMO orbitals. Smaller values of chemical potential, hardness and energy gap can facilitate the charge transfer process and indicate enhanced reactivity of a compound. The electrophilicity index (ω) is measured via equation ω=µ2/2η by means of the chemical hardness and the chemical potential. If a large value is obtained for the electrophilicity index, it exhibits the electrophilic character for the structure [44, 45].
According to the results of Table 2, CP-RGD-6, CP-RGD-7 as well as CP-RGD-8 have the lowest binding energies but since the CP-RGD-6 compound has the smallest values of hardness, chemical potential and energy gap, it has the highest affinity to bind the cancer cells/tissues. Furthermore, the CP-RGD-6 displays the utmost electronegativity and electrophilicity index values that demonstrate its highest binding affinity. Therefore, it can be decided that the CP-RGD-6 is the most suitable DDS.
Fig. 3 shows that the energy difference between the HOMO orbital of the CP and the LUMO orbital of the RGD is smaller compared to that of the HOMO orbital of the RGD and LUMO orbital of the CP; hence, electron transfer occurs from the CP to the RGD. Furthermore, the partial number of electrons transferred (ΔN) is 0.01 that confirms electron transfer happens from the drug to the carrier. It is necessary to mention that ΔN is determined from the formula ΔN=(χRGD-χCP)/2(ηRGD+ηCP) where χ is the electronegativity [44, 45].
QTAIM analysis
In molecular systems, the QTAIM describes various intra- and inter-molecular interactions mainly using the electron density topology [64]. The computed electron density values including ρ(r), Laplacian of electron density, ∇2ρ(r), total electronic energy, H(r), potential energy, V(r), kinetic energy, G(r) along with |V(r)|/G(r) for several bond critical points, BCPs, are listed in Table S2 for the CP anticancer drug, the RGD carrier and the CP-RGD systems 1, 6, 7, 8 and 10. Furthermore, the QTAIM representations of these compounds are depicted in Fig. 4.
Based on this theory, large values for the ρ(r), H(r)<٠ and ∇2ρ(r)<٠ indicate the covalent bonds while small values for the ρ(r), H(r)>٠ and ∇2ρ(r)>٠ denote the non-covalent bonds/interactions [64]. Also, if ρ(r) values are larger than 0.1 au, mainly covalent bonds are formed that are usually along with ∇2ρ(r) values. Nevertheless, for chiefly electrostatic interactions, typically the ρ(r) values are smaller than 0.1 au and ∇2ρ(r) values are positive. The H(r) values of BCPs are used as alternatives for the ∇2ρ(r) values which can be more appropriate indexes to well analyze various interactions. The H(r) and G(r) values are calculated according to the equations H(r)=G(r)+V(r) and (1/4)∇2ρ(r)=٢G(r)+V(r), respectively [65].
Table S2 reveals large values for the ρ(r) at the BCPs of the P=O, P-O, P-N, N-H, C-N, C=N, C-Cl, C-C, C-O and O-H bonds (from 0.165 to 0.413 ea0−3) which confirm their covalent natures. However, the Laplacian of electron density values are positive (from 0.030 to 0.517 ea0−5) at the BCPs of the C-N, O-H, N-H, C-C and C-Cl bonds but they are negative (from -0.004 to -0.371 ea0−5) at the BCPs of the P=O, P-O, P-N, C=N and C-O and bonds. As mentioned above, the H(r) values can be used to more exactly understand the natures of bonds. Here, the negative values of H(r) at the BCPs changing in the range of -0.017 to -0.729 ea0−4 illustrate that these bonds have a covalent character. Moreover, the cyclophosphamide drug and the RGD carrier are held together through hydrogen bonds with small ρ(r) values at the BCPs of the H-bonds that are in the range of 0.021 to 0.059 ea0−3 which corroborate the electrostatic nature of these bonds. The hydrogen bond distances between CP and RGD and the values of ρ(r), G(r), V(r) and H(r) at their bond critical points in structures 1, 6, 7, 8 and 10 are listed in Table 5. Comparison of the hydrogen bond distances shows that structures 7 and 10 have the smallest and the largest values, respectively. This result is in agreement with the data obtained for the ρ(r) values, meaning the highest and lowest electron densities belong to structures 7 and 10, respectively.
It is known that the bond energy is greater for a smaller bond length. Therefore, it is expected that structures 7 and 10 with the shortest and longest hydrogen bond distances will exhibit the highest and lowest bond energies, respectively. It is obvious in Table S2 that the total energy (H) and potential energy (V) at the hydrogen bond critical points formed between the CP and the RGD in these structures, change such that the highest and lowest values belong to compounds 7 and 10, respectively. Moreover, the kinetic energy (G) values for these compounds vary as CP-RGD-10 < CP-RGD-6 < CP-RGD-1 < CP-RGD-8 < CP-RGD-7. This trend is consistent with the results obtained for the H-bond distances, electron densities, potential energies and total energies.
Figs. 4a-4g evidently displays the formation of both inter- and intra-molecular interactions in the CP, RGD, CP-RGD-1, CP-RGD-6, CP-RGD-7, CP-RGD-8 and CP-RGD-10, respectively. Table S2 presents the QTAIM parameters calculated by the analysis of data extracted from the BCPs of such interactions. The CP drug demonstrates one intra-molecular CH…HC interaction while the RGD carrier has one intra-molecular HN…HN interaction. The CP-RGD-1 compound reveals four intra-molecular CH…HC, O…HC, NH…OC and CH..CH interactions and in addition to one inter-molecular NH…OC hydrogen bond. In the CP-RGD-6 compound, five intra-molecular CH…HC, N…HC and NH…OC plus one inter-molecular NH…N interactions are visible. In the CP-RGD-7, there are four intra-molecular CH…HC, N…HC, NH…OC and one inter-molecular NH…O interactions. The CP-RGD-8 reveals two intra-molecular CH…HC and NH…OC along with one inter-molecular NH…O interactions. There are five intra-molecular CH…HC, CH…HN, N…HCl, N…HN, N…HO and two inter-molecular O…HN hydrogen bond interactions in CP-RGD-10. The positive small values of both ρ(r) and ∇2(r) measured at the BCPs for all interactions prove that they have a non-covalent character.
It has been known that the |V(r)|/G(r) ratio can be used as an important factor in order to examine the nature of an interatomic interaction [64]. For this purpose, a non-covalent interaction reveals that H(r)=G(r)+V(r)≥0 or |V(r)|/G(r)≤1, an intermediate interaction depicts that the 1<|V(r)|/G(r)<2 but a covalent/shared interaction shows that |V(r)|/G(r)>2. Table S2 illuminates the |V(r)|/G(r) values and verifies that the C–O, C–N, C–C, C–Cl and N–H bonds have covalent properties, the P=O, P–O and P–N bonds indicate an intermediate nature whereas all intra- and intermolecular hydrogen bonds as well as the intermolecular CH…HC dihydrogen bonds exhibit electrostatic features.
Frontier molecular orbitals
The energies of the HOMO and the LUMO orbitals and their energy differences have major roles in several chemical interactions. The HOMO orbital contains donor electrons and LUMO represents the ability to accept electrons. It should be noted that the HOMO orbital energy is directly corresponds to the ionization potential but the LUMO energy is exactly depended on the electron affinity [66]. Fig. 5 illustrates the HOMO and LUMO orbitals of all structures. It is observed that the HOMO and LUMO orbitals in the CP drug are positioned on the N2 atom, the two chains connected to this atom and oxygen atoms but no LUMO orbital is seen on the oxygen atoms. The HOMO orbital of the RGD is seen onto the arginine’s nitrogen atoms while the LUMO orbital is observed on the glycine carbonyl (C=O) group and aspartic acid atoms. In all drug delivery systems, CP-RGD-1 to CP-RGD-20, the HOMO orbital is placed on the CP while the LUMO orbital distribution is perceived on the RGD and often near the sites it is connected to the drug. These results are consistent with the data obtained from the ΔN value and confirm that the electron transfer happens from the drug to the carrier.
Electronic density of states (DOS) spectra
The changes in band gaps after bonding of drug onto the carrier can be seen in electronic DOS spectra. Fig. 6 indicates the DOS spectra of the CP, RGD, CP-RGD-1, CP-RGD-6, CP-RGD-7, CP-RGD-8 and CP-RGD-10. It is observed that energy of the HOMO orbital of the CP drug in all systems is enhanced upon its attachment onto the tripeptide RGD whereas the energy of the LUMO orbital of the CP is declined. For example, the HOMO and LUMO energies of the CP change from -7.02 and -0.31 eV to -6.58 and -0.69 eV, respectively, after it’s binding onto the carrier in the CP-RGD-7. These spectra obviously illustrate that the band gap energy reduces once the CP drug is attached to the RGD carrier. In addition, Fig. 7 exhibits the PDOS spectra of structures 1, 6, 7, 8 and 10. It is found that the RGD has the most contribution in the PDOS spectra while less intensity is noticed for the CP drug.
Electrostatic surface and contour maps
The electrostatic surface potentials, ESPs, and contour maps depict the distributions of molecular electron densities. Figs. 8 and 9 illustrate the contour maps and ESPs for the CP, RGD and structures 1, 6, 7, 8 and 10. The yellow and orange lines and/or surfaces in these figures specify the distributions of positive and negative charges, respectively. It is seen that the distribution of the positive and negative charges in the CP drug are symmetric but the symmetry is not observed in the RGD. After bonding of the CP to the RGD in structures 1, 6, 7, 8 and 10, the charge distribution becomes asymmetric leading to increasing the dipole moments of these compounds compared to those of the CP and the RGD. For instance, in the CP-RGD-7 structure, that shows the highest dipole moment (14.87 D), the most asymmetric distributions of charges and surface potentials are observed.
MD simulations
Relaxation and equilibrium of cells
The effect of adding PEG polymer chains into the CP-RGD drug delivery system (DDS) on diverse properties of CP-RGD-PEG systems was evaluated in ethanol solvent using the MD simulations. The relaxation of each system was examined by achieving almost constant temperature as well as non-bond, potential and kinetic energies. Figs. 10a-10k demonstrate the relaxed cells of the CP-RGD, CP-RGD-PEG, CP-RGD-2PEG, CP-RGD-3PEG, CP-RGD-4PEG, CP-RGD-5PEG, CP-RGD-6PEG, CP-RGD-7PEG, CP-RGD-8PEG, CP-RGD-9PEG and CP-RGD-10PEG, respectively, attained by the MD simulations. Table 6 provides some simulation results including the density and cell dimension (size) for the CP-RGD-PEG cells. It is observed that both of the cell dimension as well as the cell volume are enhanced once more PEG chains is added to the cell so that the CP-RGD (14724.14 Ǻ3) and CP-RGD-10PEG (27735.58 Ǻ3) reveal the smallest and largest cell volumes, respectively. The CP-RGD has a density of 0.999 g/cm3 while higher density is measured upon adding more PEG chains into the system as the greatest density (1.068 g/cm3) belongs to the CP-RGD-10PEG cell. This result is associated with the enhancement of both cell mass and cell volume when additional PEG chains are introduced into the cell.
Fig. 11 exhibits the changes in the non-bond, potential and kinetic energies in addition to the temperature against simulation time for the CP-RGD cell. It is found that negligible changes occur in all kinds of energies with time confirming the system has fully been relaxed. Analogous data are calculated for all of the CP-RGD-PEG systems. Table S3 presents the average non-bond, potential and kinetic energies achieved for the CP-RGD-PEG cells containing various amounts of PEG chains. For all cells, the kinetic energies reveal positive values and the non-bonded energies are more negative compared to those of the potential energies. A comparison of the kinetic, non-bond and potential energies of different systems proves that increasing the number of PEG leads to more positive energies. For instance, the average non-bond, potential and kinetic energies are -2925, -2909 and 1285 kcal/mol for the CP-RGD cell while those of the CP-RGD-10PEG system are -2554, -2528 and 2547 kcal/mol, respectively.
Surface areas and free volumes
Table 6 gives the surface areas for different CP-RGD-PEG cells demonstrating the CP-RGD has the minimum surface area (3601.67 Å2) whereas the CP-RGD-10PEG indicates the utmost surface area (8525.98 Å2). Also, further addition of the PEG chains to the cells leads to enhancement of the surface area.
The free volume (FV) values for various CP-RGD-PEG systems are collected in Table 6. As well, Figs. 12a-12k illustrate the free volumes with blue color for the CP-RGD, CP-RGD-PEG, CP-RGD-2PEG, CP-RGD-3PEG, CP-RGD-4PEG, CP-RGD-5PEG, CP-RGD-6PEG, CP-RGD-7PEG, CP-RGD-8PEG, CP-RGD-9PEG and CP-RGD-10PEG cells, respectively. It is found that adding more PEG chains to the systems increases the free volume so that the CP-RGD and the CP-RGD-10PEG illuminate the lowest (2903.78 Å3) and the highest (6898.34 Å3) values. Such result is associated with the strongest intermolecular interactions including electrostatic and hydrogen bond attractions happened among different particles inside the CP-RGD-10PEG cell, i.e. the CP drug, the RGD, the PEG chains and the ethanol solvent.
The fractional free volume (FFV) equals the FV divided by the total volume of cell. Among different CP-RGD-PEG systems, the largest FFV (24.87%) is measured for the CP-RGD-10PEG but the smallest FFV (19.72%) is obtained for the CP-RGD. Obviously, like the free volumes, the fractional free volumes are increased upon incorporation of extra PEG chains into the cells. The utmost FFV for the CP-RGD-10PEG cell might be accounted for the simplest CP drug transport in this cell; thus, it affords the maximum drug diffusion coefficient. Nevertheless, to choose the most appropriate drug delivery platform, further properties of systems like mean square displacements and diffusion coefficients need to be estimated.
Mean square displacement (MSD) and diffusion coefficients
The CP drug diffusion inside the CP-RGD-PEG systems was studied through achieving the MSD plots presented in Fig. 13. The approximately linear MSD lines demonstrate that constant CP transports occur within the systems throughout the simulations. Among various cells, the CP-RGD and the CP-RGD-10PEG show the greatest and the lowest MSDs for the CP diffusion, respectively. As an effective drug delivery illustrating the most controlled CP drug carriage should reveal a low MSD, it may be decided that the CP-RGD-10PEG system is the most suitable cell affording the slowest and the most controlled CP diffusion/transport.
Table 7 gives the diffusion coefficient values for the diffusion of the CP molecules in diverse systems measured from the slopes of the MSD curves. The diffusion coefficients are in agreement with the MSD plots indicating among all of the cells, the CP-RGD and the CP-RGD-10PEG exhibit the largest and the smallest diffusion coefficient values of 0.0216×١٠–5 and 0.0067×١٠–5 cm2/s, respectively. Hence, the slowest CP drug diffusion occurs inside the CP-RGD-10PEG system that can result in the most effective and controlled drug delivery. However, a proficient drug carrier system must show a slow drug release rate, suitable FV and FFV values plus a relatively high drug loading. Therefore, the CP-RGD-6PEG is introduced as the most appropriate DDS exhibiting desirable FV (4988.89 Å3) and FFV (22.66%) values, relatively high drug loading as well as small drug diffusion coefficient (0.0114×10–5 cm2/s) indicating low drug release rate.
CONCLUSION
DFT calculations were accomplished on twenty CP-RGD structures formed through strong hydrogen bonds created between the tripeptide RGD and the anticancer drug cyclophosphamide. The QTAIM analysis certified that the C–O, C–N, C–C, C–Cl and N–H bonds had the covalent nature, the P=O, P–O and P–N bonds revealed the intermediate property whereas all intra- and intermolecular hydrogen bonds as well as the intermolecular CH…HC dihydrogen bonds had the electrostatic character. Both of the contour maps and ESPs revealed that after binding of the CP to the RGD in structures 1, 6, 7, 8 and 10, the charge distribution became asymmetric leading to increasing the dipole moments of these compounds compared to those of the CP and the RGD. As the CP-RGD-6 illustrated the smallest values of hardness, band gap energy and chemical potential whereas the utmost electrophilicity index and electronegativity, it had the greatest affinity to attach the cancer cells. Hence, the CP-RGD-6 was chosen as the most suitable DDS for the cancer treatment application. The molecular dynamics simulations confirmed that among diverse CP-RGD-PEG cells, the CP-RGD-6PEG could be introduced as the most promising DDS exhibiting suitable FV (4988.89 Å3) and FFV (22.66%) values, relatively high drug loading as well as small diffusion coefficient (0.0114×10–5 cm2/s) indicating low drug release rate.
ACKNOWLEDGEMENTS
The financial support of this research by the Research Office of Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran, is thankfully appreciated. Furthermore, authors are indebted to High Performance Computing Research Center (HPCRC) of Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran for providing computational facilities (software and hardware) to conduct the MD simulations.
CONFLICTS OF INTEREST
The authors do not have any personal or financial conflicts of interest.
)
