Introduction

Despite recent improvements in cancer treatment, most anticancer drugs affect healthy cells due to their non-selective action. Therefore, researchers are investigating ways to destroy cancer cells in a targeted manner. Targeted therapy includes various methods, such as antibody-based targeting, biotin-based targeting, peptide-based targeting, etc1,2. Among the various cellular targets, integrins are a group of overexpressed receptors in all types of cancer cells and tumors that have a crucial role in many biological processes, including interaction with extracellular matrix proteins3. Many studies have shown that these receptors are overexpressed on the cell surface of malignant tumors.

In 1970, the tripeptide RGD (Arg-Gly-Asp), which acts as a cell adhesion ligand, was identified as an integrin binder. RGD is one of the shortest sequences for integrin recognition1,4. To improve the efficiency of cytotoxic agents against cancer cells and reduce side effects on healthy cells, the use of carriers or mediators like RGD can be helpful5. Using the RGD sequence as a carrier, cytotoxic drugs can be selectively delivered to cancer cells because these cells have the highest number of integrin receptors on their surface. This method enables more targeted and selective cancer treatment and helps to reduce damage to healthy cells6.

Conjugated RGD sequence (attachment of cytotoxic drugs to RGD sequence) is clinically significant. Because this sequence can act as a target for cancer cells receptors7,8. As well as, the combination of cytotoxic drugs with RGD sequence may show a synergistic effect in anti-angiogenesis (inhibition of blood vessel growth)9. This means the combination of these two substances may enhance the anti-angiogenic efficacy.

In vitro and in vivo studies, some non-steroidal anti-inflammatory drugs (NSAIDs), such as cyclooxygenase inhibitors, have shown anticancer properties10,11,12. The anticancer properties of NSAIDs appear to occur through inhibition of one of the cyclooxygenase isozymes, but due to side effects, NSAIDs cannot be used as anticancer agents. Hence, by conjugating NSAIDs to the RGD sequence as a tumor-targeting ligand, they can be selectively delivered to cancer cells13,14.

Numerous studies have been conducted to develop a method of targeting the integrin receptor, which is expressed on the surface of cancer cells, by modifying the surface of carriers. In order to effectively deliver these compounds to the targeted site, a suitable carrier is required. Nanostructures or nanoparticles can be attached to the target cell and can destroy the cancer cell. Carbon nanotubes (CNTs) are one of the most widely used nanostructures in biomedicine and medicine for drug delivery and diagnosis15,16,17. Due to their outstanding structure, high surface area, and excellent chemical and thermal stability, these compounds can penetrate the cell membrane and deliver peptides, nucleic acids, and other substances to the cells18.

In an experimental study, Fan et al. investigated the effect of carbon nanotube-based anticancer prodrugs. Their results showed that CNT can be a suitable carrier to transfer certain drugs to the target cells and increase the effectiveness of treatment19.

Maleki et al.. performed molecular dynamics simulations to explore the adsorption and release of doxorubicin from single-walled carbon nanotubes and multi-walled carbon nanotubes. Their findings showed that carbon nanotubes can be effective drug carriers for the transport and release of doxorubicin in cancer tissues20. Mohammadi et al.. reported that conjugation of non-steroidal anti-inflammatory drugs with peptide sequences facilitates targeted delivery and reduces toxicity21.

In the present study, the use of CNT as a carrier for the adsorption and transfer of non-steroidal anti-inflammatory drugs ketoprofen and naproxen in the presence of RGD peptide sequence is investigated with the help of molecular dynamics simulation. Additionally, we will discuss the dynamics and diffusion of the designed systems in the presence of cell membranes, along with the interactions of these complexes with double-layer membranes. This work can be helpful to experimental methods in the treatment of many diseases in the future.

Simulation details

To gain insight into the biological activity of single-walled carbon nanotubes as a nanocarrier, the dynamic adsorption process of peptide-drug conjugates on its surface is investigated by MD simulation.

At first, molecular dynamics simulation was performed to check the efficiency of SWCNT nanocarrier for loading of Ket-RGD and Nap-RGD complexes. Peptide-drug complexes were prepared using GaussView software22. The peptide sequence of RGD was connected to the carboxyl group of the drug molecule through the amine group. Then, to find the best orientation of the complex, it was optimized with the help of Spartan software23. In both systems, the nanocarrier was positioned at the center of simulation boxes with dimensions of 9 × 11 × 9 nm³. The peptide-drug conjugates were then placed at a specific distance (≈ 2 nm) from the nanocarrier. Additionally, the TIP3P water model molecules were added to two simulation boxes.

Zigzag configuration of SWCNT carbon nanotube consisting of 936 carbon atoms and 36 hydrogen atoms, with chirality index (10, 0)24, tube diameter of 1.4 nm, and length of 5.38 nm as a suitable substrate using the Nanotube Builder plugin package VMD was prepared. The edge carbon atoms were saturated with hydrogen to avoid boundary effects. All CNT carbon atoms were considered as uncharged Lennard-Jones particles. Two drugs, Ketoprofen and Naproxen, were each complexed with RGD (arginine-glycine-aspartic acid) peptide sequence, and a structural search was performed using the Spartan package to find optimal structures (Fig. 1).

Fig. 1
figure 1

The structure of drugs (a) ketoprofen, (b) naproxen and (c) peptide sequence of RGD.

Full size image

Molecular structures of drugs Ket and Nap are provided from the PubChem database25. Force field parameters of drug molecules and topology were obtained from the Swiss Param server26.

Considering that the drug delivery performance of CNT depends on its ability to diffuse into the cell membrane, the interaction of Ket-RGD/SWCNT and Nap-RGD/SWCNT complexes with lipid bilayers has been investigated. To investigate this process, the final orientation of Ket-RGD and Nap-RGD molecules on the carrier were extracted from the loading system. Next, these assemblies were placed in a parallel orientation above the surface of the membrane at a distance of ≈ 2 nm. The dimensions of designed simulation boxes are 11.44 × 11.44 × 20 nm3. The selected membrane in this study contains 400 lipid palmitoyloleoyl phosphatidylcholine (POPC) and 100 cholesterol (CHL), about 20% of the total membrane27. The membrane was constructed using the protocol described on the CHARMM-GUI28 website and equilibrated for 30 ns. The GROMACS software package (version 2019.2)29 was used to perform all the MD simulations in this study. CHARMM36 force field30 is chosen to obtain carrier and membrane parameters. Figure 2 shows the chemical structure of the CNT nanotube, Ket-RGD molecule, Nap-RGD, and bilayer membrane.

To imitate the physiological conditions, water molecules were explicitly used with the TIP3P model, and 0.15 M sodium chloride (NaCl) was added in to the simulation boxes to neutralize the system. After the steepest initial energy reduction, the systems in the conventional set (NVT) are equilibrated at 310 K using a Nose-Hoover thermostat31. After that, they reach equilibrium in NPT set by the Parrinello-Rahman algorithm32 at 1 atmosphere.

Fig. 2
figure 2

The Structure of a, b) Ket-RGD and Nap-RGD complexes, c) carbon nanotube, d) two layers of POPC respectively.

Full size image

MD simulations for loading and membrane systems under boundary conditions are run in three directions for 100 and 300 ns, respectively. In addition, the integration time step is set to 1.5 fs. In these systems, the LINCS algorithm33 is used to maintain all links in their equilibrium position. The van der Waals (vdW.) interactions are evaluated at a cut-off distance of 1.4 nm, and the Particle Mesh Ewald (PME)34 method is used to calculate the electrostatic forces (Ele). In order to visualize the simulation results, VMD software has been used.

Free energy calculation can accurately evaluate the DDS diffusion process in the membrane. Metadynamic simulation is a technique used to calculate the free energy level. Gromacs 2019.2 software patched with PLUMED plugin version 2.5.235 was used in this work. Metadynamics uses Gaussian biases to overcome high energy barriers and is smoothed by imposing age-long limitations on the bias potential. To determine the free energy level, one can use Eq. 136.

$$\:{V}_{bias}\left(s'\:t\to\:\infty\:\right)=-\frac{\varDelta\:T}{T+\varDelta\:T}\:F\left(s\right)+C$$
(1)

In this relationship, s is collective variables, t is time, F(S) is free energy and ∆T is a parameter with temperature dimension that adjusts the speed of FES discovery.

Results and discussion

In this section, we investigate the adsorption behavior of Ket-RGD and Nap-RGD complexes on CNT and their interaction with cell membranes through molecular dynamics simulation using a step-by-step approach.

Adsorption of Ket-RGD and Nap-RGD complexes on CNT

This study initially focuses on investigating the adsorption mechanism of Ket-RGD (KRG) and Nap-RGD (NRG) complexes with CNT to understand their interaction. Figure 3 shows images of the studied systems. Our results show that both complexes have a great tendency to adsorb on the CNT surface.

Fig. 3
figure 3

Initial and final snapshots of the investigated systems, (a) ketoprofen and (b) naproxen conjugated with RGD sequence in the presence of CNT.

Full size image

Upon examining the snapshots, it becomes clear that as the simulation progresses, the distance between the drug molecules and the CNT surface decreases, ultimately resulting in the adsorption of the molecules on the carrier surface. Furthermore, it is apparent that due to the presence of aromatic rings in the structure of drug molecules and the CNT surface (see Figure S2), the adsorbed drug molecules interact with the carrier surface via π-π stacking37,38, and the RGD sequence is physisorbed entirely to the carrier surface. Also, the distance between the aromatic rings of the drug and CNT is in the range of the distance reported in the studies in which π-π stacking exist39,40. In order to interpret the adsorption process, the Lennard-Jones interaction energy values during the MD process are calculated and shown in Fig. 4.

Fig. 4
figure 4

The interaction energy between KRG، NRG complexes and carrier as a function of time.

Full size image

According to Fig. 4, at the beginning of the simulation, drug molecules are located farther than the cut-off radius and have no interaction with the substrate, and the Lennard-Jones energy is zero. For the Ket-RGD system, the energy has a downward trend until about 50 ns, and after that, due to the increase of interactions between the surface of the drug and the CNT, the energy has reached its lowest value.

In the case of the Nap-RGD system, as the drug molecules approached the surface of the carrier, the energy decreased until about 60 ns with a gentle slope, and then due to the change in the orientation of the drug on the surface of the carrier, there was a decrease in the amount of energy, then in time 80 ns has reached its lowest value.

In order to gather additional information about the adsorption process, we have calculated radial distribution functions (RDF). These functions display the likelihood of finding drug molecules at a given distance from the CNT. Figure 5 shows the RDF diagram for the examined systems.

According to the shape of the main peak of the RDF diagram for both compounds, which appeared at a distance of 0.5 nm, it is obvious that the peak intensity for KRG is greater than that of NRG, which confirms more presence of KRG around the nanotube. The existance of more π-π and van der Waals interactions between KRG and CNT has caused the peak intensity of this drug.

Fig. 5
figure 5

Radial distribution function (RDF) for two KRG, and NRG complexes in the presence of CNT.

Full size image

Figure S1 shows the drug density diagram around the substrate. It can be seen that the drug molecules are distributed in the range of 1.5 to 0.5 nm around the substrate. It is clear that the drug density of KRG is higher than the other system, which is consistent with the energy results.

To investigate the drug release behavior around the CNT, MSD diagrams as a function of simulation time are calculated and depicted in Fig. 6. Furthermore, presents the self-diffusion coefficient values of drug molecules compared to CNT using Einstein’s relation (Eq. 2)41.

Fig. 6
figure 6

Mean square displacement (MSD) of studied complexes against simulation time with corresponding Di values.

Full size image
$$\:\mathop {{\text{lim}}}\limits_{t \to \:\infty \:} \langle ||{r_i}\left( t \right) - {r_i}\left( 0 \right)|{|^2}\rangle > = 6Dt$$
(2)

where the expression in the bracket <…> indicates the MSD of the drug and ri(0), ri(t) are the positions of the particle at time 0 and t, respectively.

According to the Fig. 6, it is clear that the KRG complex has a lower diffusion speed than the NRG complex around the CNT. For both systems, it can be seen that initially, the graphs show an increasing trend over time, but the average diffusion coefficient in the time range studied in CNT-KRG is lower than in the CNT-NRG system (see Fig. 6). This fact confirms that the mobility of the adsorbed KRG complex on the CNT surface is less than that of the other system.

In order to provide further insight into the extent of adsorption, the number of contacts between the KRG and NRG complexes with the nanotube is calculated and results shown in Fig. 7. A detailed analysis of this figure shows that the complexes consistently approach and contacts with the CNT surface, with these collisions gradually increasing over time. It can be seen the KRG system moves towards the substrate from the beginning of the simulation (~ 9 ns) and starts interacting with it. Over time, the KRG complex has approached the CNT surface, and as a result, the number of collisions has increased.

Fig. 7
figure 7

The number of contacts between the studied complexes and CNT during the simulation time.

Full size image

After approximately 50 ns, the structure is aligned parallel to the surface of the nanotube and most of its atoms have been involved in the interaction, leading to an increase in the number of collisions. In the case of the NRG system, there was no interaction with the CNT surface until about 20 ns. After that, the number of collisions increased with a gentle slope with the complex’s adsorption on the substrate’s surface, and from approximately 95 ns to the end of the simulation time, there was the highest number of contacts between the drug and the substrate. There is a good agreement between the number of collision values and other obtained results.

The permeability of the drug delivery cell in the double-layer membrane

The final orientation of KRG and NRG complexes on CNT was extracted from the last moments of the path of the loading system and used to investigate its interaction with the bilayer membrane. According to Fig. 8, this set is parallel to the membrane cell above it (0 ns). Both systems of the DDS move from the water phase to the membrane during the simulation. The CNT-KRG complex remains parallel to the membrane’s surface and is physically adsorbed on its surface, and in approximately 190 ns, the complex penetrates the membrane (See Figure S3). But about the CNT-NRG complex was observed, after 15 ns passed from the start of the simulation, the DDS changed its direction and was placed almost vertically on the membrane. This sudden change in direction may be due to the accumulation of drug molecules on one side of the carrier. Also, owing to the competition between water and membrane interacting with DDS, the system remains in this situation for 30 ns. Eventually, the DDS fully penetrated the bilayer membrane after 50 ns (Figure S4).

Fig. 8
figure 8

Snapshots of the interaction of DDS with the cell membrane at the beginning and end of the simulation: (a) CNT-KRG, (b) CNT-NRG.

Full size image

In order to study how DDSs interact with the membrane, we calculated the van der Waals energy values and the number of contacts (See Fig. 9). The results of energy calculations indicate that the diffusion process through the lipid bilayer occurred spontaneously. A close inspection of Fig. 9 panel a reveals that the energy changes are relatively small in the first 20 ns due to less interaction of the DDS with the membrane. But around 20–35 ns for the CNT-KRG system, the energy decrease becomes steeper, which indicates the movement of DDS towards the membrane and the increase of interaction with the membrane.

Fig. 9
figure 9

(a) Van der Waals energy and (b) number of contacts between the two studied complexes and the cell membrane as a function of simulation time.

Full size image
Table 1 Average electrostatic interaction energies (Ele) and Van Der Waals (L-j) in the studied systems.
Full size table

The interaction of DDS with the nitrogen and oxygen groups of POPC increases the electrostatic interactions and stability of this complex on the membrane. It can also be seen that at approximately 195 to 230 ns, we have a sharp decrease in the amount of van der Waals energy, which can be caused by the penetration of DDS into the membrane and the increase of its interaction with POPC groups. Inspection of the CNT-NRG system shows a sudden decrease in the energy level in the time interval of 20 to 50 ns (see Fig. 9 panel a). This reduction is due to the high tendency of DDS to interact with the carrier, where van der Waals and electrostatic interactions contribute to this interaction, but as it is evident in Table 1, van der Waals interactions play a decisive role. Finally, the energy reaches a constant level and fluctuates around an average value (-4200 kJ/mol). In addition, the number of contacts between the membrane and DDSs supports this conclusion (see Fig. 9 panel b). It can be seen that with the passage of simulation time, the number of collisions has increased in both systems. Therefore, it can be said that the use of a drug-conjugated peptide sequence (RGD) in the presence of CNT is suitable for transferring the drug to the membrane surface.

Free energy calculations can be used as a practical method in computational chemistry to study different biological systems. Metadynamics is a technique used to calculate the energy of long-term processes, enabling the system to escape deep energy wells and regenerate the free energy level. The calculation of FES, which is a function of the distance between the center of mass of the membrane and DDS, is done using the sum_hills program. Figure 10 provides a visual representation of this calculation.

In addition, In Fig. 10, density diagrams along the z-axis are presented to provide a deeper understanding of the studied systems. It is necessary to mention that the “gmx density” module was applied to calculate density and the POPC membrane was selected as a reference. The analysis of the density plots indicates that, at first the CNT-KRG and CNT-NRG complexes are approximately 3–4 nm away from the membrane. The obtained results show that both complexes have a great tendency to adsorb on the surface of the membrane and even penetrate it. A close detection of Fig. 10 also shows that both DDS move from the water phase to the membrane. In addition, the drug and CNT density graphs overlapped until the end of the simulation.

Fig. 10
figure 10

Density profiles of POPC, water, drug and carrier at different times and free energy profile for DDS penetration into cell membrane, (a) CNT-KRG, (b) CNT-NRG.

Full size image

Figure 10 panel a shows that as the complex moves towards the membrane, the energy decreases and oscillates around a constant value near the membrane surface until reaching a global minimum (about − 175 kJ/mol). The minimum point for the CNT-NRG system is approximately 1 nm from the center of the membrane (about − 183 kJ/mol) (see Fig. 10 panel b). It is worth noting that the drug molecules (KRG and NRG) and CNT both played a role in the adsorption of DDS on the surface of the membrane and did not face any energy barrier to reach the minimum point of these systems.

Conclusion

In this study, the ability to transport Ketoprofen and Naproxen drug molecules conjugated with RGD peptide sequence and using carbon nanotubes as their carrier was studied by MD simulation. The results of the simulation process show that in the most stable state of the system, the drug molecules choose an orientation parallel to the CNT surface and form van der Waals and π-π interactions. In addition, the diffusion process of CNT-KRG and CNT-NRG complexes to the cell membrane has been investigated. Also, the results indicate that the main reason for the ability to adsorb and penetrate the membrane surface of these complexes is the formation of π-π stacks. The results of metadynamic calculations show that as the complexes approach the surface of the membrane, the free energy values decrease, and its global minimum is around − 183 kJ/mol. In this way, the use of carbon nanotubes in the presence of a peptide sequence (RGD) provides new insight into transferring non-selective drugs to the membrane surface.