Self-assembled xK1 rosette nanotube as a nanocarrier for chemotherapeutic drug molecules: a molecular dynamics study

Abstract

Chemotherapy is widely recognized as an effective treatment for various types of cancers, relying on chemotherapeutic drugs such as chlorambucil (CBL), camptothecin (CPT), doxorubicin (DOX), and flutamide (FLU) to inhibit cancer cell metastasis through their cytotoxic effects. However, these drugs also interact with noncancerous cells, causing significant side effects, including but not limited to nausea, vomiting, hair loss, peripheral neuropathy, lymphedema, and infertility. Developing nano-based drug vehicles is one approach to alleviate such adverse effects. This study introduces a three-ring-xK1-type rosette nanotube (RNT) as a potential nano-based drug-delivery vehicle for these chemotherapeutic agents. Molecular dynamics (MD) simulations aided by the molecular mechanics/three dimensional reference interaction site model (MM/3DRISM) analysis were employed to evaluate the potential of xK1 as a drug carrier. Results revealed that xK1 exhibited structural stability and strong binding affinity with the selected drugs under physiological conditions, with binding energies ranging from −4.66 to −26.10 kcal mol⁻¹. The drugs favorably bind at either end of xK1 due to the prominent aromatic rings facilitating π–π stacking interactions. Furthermore, MD simulations of the release mechanism revealed structural disassembly, validating the responsiveness of the systems under acidic conditions. The findings here offer insights into the viability of xK1 for targeted drug delivery, which warrants further experimental studies.

---

Introduction

The main cancer treatment modality utilized worldwide is chemotherapy, where drugs are intravenously or orally administered to cause cancer cell death. The primary goal of chemotherapy is to use the cytotoxic effects of drugs to halt the spread of cancer cells. However, it also affects non-cancerous cells, causing adverse side effects on patients. Reports indicate that this treatment significantly reduces the quality of life for cancer patients, with physical symptoms—particularly chemotherapeutic side effects—being the strongest factor.¹

Efforts have been made to alleviate the severe side effects of chemotherapeutic drugs by employing strategies, such as combining them with bacterial metabolites,² regulating redox balance using natural products,³ utilizing natural medicines to treat side effects,⁴ and implementing cryotherapy to mitigate the side effects during chemotherapy.⁵ One approach is to develop nano-based drug vehicles, which can alleviate the severe side effects by enabling controlled drug release, protecting the healthy cells from exposure, crossing the biological barriers to target cancer cells directly, and improving the stability and biocompatibility within the bloodstream.^6–8 Numerous studies have investigated the application of nanotechnology in drug vehicles, including carbon-based nanotubes or nanosheets,^9–11 liposome-based¹² and chitosan-based¹³ nanomaterials, carrier proteins,¹⁴ gold nanoparticles,¹⁵ metal-based nanoparticles,¹⁶ quantum nanodots,^17,18 and others.^7,19 A study reported the suitability of chitosan nanoparticles as drug vehicles owing to their physicochemical stability and bioavailability.²⁰ Dendrimers have also gained increasing attention as drug-delivery vehicles owing to their physicochemical features, such as low polydispersity, nano-scale dimensions, high solubility, and the ability to conjugate various molecules through their branched cavities and surface groups.²¹ A significant challenge in the delivery of chemotherapeutic drugs is their limited water solubility. This is commonly addressed by modifying the hydrophobic carriers through hydrophilic attachments or surface functionalization. Among nano-based drug vehicles, carbon-based drug vehicles have been extensively studied. However, these systems can cause kidney injuries during the excretion process, particularly with larger-sized molecules.²² Additionally, studies have shown that they tend to accumulate in organs, such as the lungs, liver, spleen, and kidneys, for prolonged periods, potentially leading to long-term toxicity.^22,23 Thus, there is a need for a nanocarrier that is biocompatible and has a high therapeutic efficacy.

This study explores the potential of rosette nanotubes (RNTs) as a drug vehicle. RNTs exhibit significant properties, such as biocompatibility for biomedical applications, amphiphilicity, low toxicity, and self-assembly behavior, making them potential candidates for the delivery of chemotherapeutic drugs.²⁴ These properties are useful for enhancing the loading, delivery, and solubility of hydrophobic drugs in aqueous environments. RNTs are supramolecular structures composed of self-assembling motifs designed from guanine (G) and cytosine (C) base pairs. These base pairs display complementary hydrogen bonding patterns, namely donor–donor–acceptor for guanine and acceptor–acceptor–donor for cytosine. Through self-assembly, the G∧C motif (Fig. 1a) forms rosette-shaped rings (Fig. 1b), which subsequently stack into stable, hollow nanotubular structures (Fig. 1c) via forces such as hydrophobic interactions, dispersion forces, and stacking interactions.^25,26 A key feature of RNTs is their modifiable surface, which can accommodate a variety of functional groups. This feature enhances their solubility and bioavailability.²⁴ The tubular design of RNTs includes a hydrophobic inner core that can encapsulate hydrophobic drugs, while the functionalized outer surface protects the drugs in physiological environments, enabling controlled drug release.^27,28 In addition, the terminal regions of the RNTs feature prominent aromatic rings, which can interact with drugs containing aromatic moieties. This characteristic further enhances the potential of the system as an efficient drug-delivery vehicle. In this study, the xK1 type of rosette nanotube (RNT) is utilized. The xK1 (Fig. 1c) type of RNT was designed through supramolecular engineering of its G∧C-based motif formed by the same G and C hydrogen bonding arrays as the bicyclic ring, but placed with a pyridine ring in between.²⁵ This was designed to achieve control over the structural parameters of the rosette nanotube (RNT). The incorporation of the pyridine ring allows the fine modulation of π–π stacking and hydrogen-bonding interactions between the rosette subunits, resulting in an increased inner diameter and enhanced hydrophobicity of the internal cavity. This tunability in both geometry and surface chemistry provides a crucial advantage for the accommodation of drug molecules of varying sizes and polarity, making xK1 a more versatile and efficient drug-carrier RNT platform. The defining characteristic of xK1 lies in its enhanced hydrophobicity and amphiphilic tubular structure, which strongly facilitates the formation of assemblies in water.²⁵

Fig. 1 Molecular model of rosette nanotube displaying the (a) motif molecule of xK1, (b) top view of the nanotube, and (c) front view of the nanotube.

The limited biomedical applications of drug nanocarriers are primarily attributed to their poor biocompatibility, renal toxicity, and long-term accumulation in the body, as discussed. In contrast, rosette nanotubes (RNTs), which self-assemble through hydrogen bonding of G∧C motifs, offer a safe and biocompatible alternative. Experimental studies successfully reported the applications of RNTs as drug vehicles of hydrophobic dexamethasone for bone regrowth,²⁸ tamoxifen as an anti-cancer drug,²⁹ and siRNA for gene silencing.³⁰ RNTs have also been utilized for tissue engineering applications as a scaffold for neural regeneration, nerve repair³¹ and cartilage repair.³² These studies suggest that RNTs are highly suitable for biomedical applications, specifically as a drug vehicle for hydrophobic anti-cancer drug molecules.

This study aimed to investigate the structural and energetic stability of xK1 under physiological conditions as a drug vehicle for chemotherapeutic agents, such as chlorambucil (CBL), camptothecin (CPT), doxorubicin (DOX), and flutamide (FLU). Most chemotherapeutic drugs are hydrophobic in nature and known to be cytotoxic. Molecular dynamics (MD) simulations aided by robust energetic analysis were employed in this study to evaluate the potential of xK1 as a targeted drug delivery system for chemotherapeutic molecules. This approach has been employed in investigating the stability and viability of a system consisting of drug-nanocarriers, as demonstrated in a recent study on zwitterion/PAMAM/CQD hybrid nanocarrier for doxorubicin (DOX).³³ The authors investigated the interactions and stability of the hybrid nanocarrier under physiological conditions by employing MD simulations and energetic analysis. Another recent study by Alkriz and Joujeh³⁴ evaluated the potential of RNA aptamers as drug vehicles for chemotherapeutic drugs by employing computational experiments, including molecular dynamics simulations. Thus, the approach in this study is suitable for investigating xK1 as a drug vehicle for chemotherapeutic drugs. Additionally, this approach offers a more efficient and cost-effective alternative to in situ experiments, providing molecular-level insights that are difficult to achieve experimentally. MD simulations allow the comprehensive analysis of stability and molecular interactions between drug molecules and RNTs, which are crucial for understanding the behavior of the drug vehicles under specific conditions. Additionally, energetic analysis using the Molecular Mechanics Poisson Boltzmann Surface Area (MM/PBSA) model and the Three-Dimensional Reference Interaction Site (3DRISM) model provides the system's energies, offering significant support for the overall analysis. MD simulations were also performed under acidic conditions to gain insights into the release mechanisms of the complexes. The results here collectively offer valuable insights into the viability of xK1 as a targeted drug vehicle for chemotherapeutic agents. The findings from this study are essential for developing safe and effective drug-delivery formulations. To the best of our knowledge, this is the first study to investigate xK1 as a drug vehicle for chemotherapeutic drugs, such as CBL, CPT, DOX, and FLU.

Experimental

System preparation

The xK1 type of RNT used in this study was modelled and prepared using Schrödinger Maestro 2021-1.³⁵ The motif molecule was optimized using the Generalized AMBER Force Field (GAFF) and subsequently multiplied into a rosette ring structure, which was further stacked in a tubular configuration. The modelling and parameterization of the RNT structure were adapted from established studies^24,36 and are detailed in SI Fig. S1.

The molecular structures of the drug molecules (CBL, CPT, DOX, and FLU) were retrieved from PubChem³⁷ and ChemSpider.³⁸ The system topologies were prepared using the GAFF parameters with partial atomic charges assigned via the AM1-BCC method generated using ACPYPE.³⁹ This allows the use of GAFF parameters in GROMACS while maintaining compatibility with AMBER-based tools, such as gmx_MMPBSA, thereby ensuring methodological consistency throughout all simulation and analysis stages. GAFF parameters were also employed for various nanocarriers^40–43 in various applications. In addition, each motif molecule was derived from guanine and cytosine bases, which were developed into a new monomer structure. This resulted in another molecule that required parameterization as a novel small molecule. GAFF, an extension of the AMBER force fields designed for drug-like and organic molecules, was therefore employed.⁴⁴ GAFF covers a broad chemical space, including common atom types, such as C, N, H, F, and Cl. Since CBL and FLU consist of Cl and F atoms, respectively, GAFF parameters are particularly suitable for them.

Molecular docking

Binding regions. Three accessible binding regions are specified for each drug: (1) in the inner channel, (2) at one end with lysine on the right side, and (3) at one end with lysine on the left side. The first configuration (Fig. 2a) is when the drug molecule is in the inner channel. The orientation is described based on the alignment of the longer side of the drug molecule along the x- and y-axes. If the longer side of the drug is aligned along the x-axis in the inner channel of xK1, then it is deemed horizontally positioned. If the longer side is aligned along the y-axis in the inner channel of xK1, then it is deemed vertically positioned. The orientation of each drug molecule along the inner channel of xK1 is based on the results from the molecular docking studies.

Fig. 2 Accessible binding regions of xK1, where the drug molecules are (a) in the inner channel, (b) on the LysR end, and (c) on the LysL end.

The two remaining binding regions are based on the ends of the channel of xK1. These regions are labeled according to the position of lysine (Lys) when viewed from the top in a counterclockwise direction. In one configuration (Fig. 2b), Lys appears on the right side of each motif within the xK1 and is labeled as LysR. In the other configuration (Fig. 2c), Lys appears on the left side of each motif within the RNT and is labeled as LysL.

Docking parameters. Autodock4.2 tool⁴⁵ was used for molecular docking, both for global and local docking. The systems were prepared using AutoDockTools-1.5.7 (ADT) prior to docking. The drugs were treated as flexible systems, while the xK1-RNT was kept rigid. The grid box dimensions and coordinates were adjusted to ensure that all accessible surfaces of the xK1 were enclosed. For local docking, a specific region was defined. Both global and local docking were utilized in this study.

The grid box for both global and local docking was centered at coordinates 50 × 50 × 49.609 Å with a grid spacing of 0.375 Å. A cubic grid box with grid points of 126 × 126 × 126 was used to define the binding sites for global docking. All drugs in the global docking revealed favorable poses at the LysR region (SI Fig. S2). Hence, local docking at LysL and within the inner channel was subsequently employed (SI Fig. S3 and S4). For the region at LysL, the grid box was defined with grid points of 110 × 78 × 26, while for the region within the inner channel, the grid box was defined with grid points of 68 × 48 × 70. The binding poses were analyzed based on the estimated binding free energy per pose. The top-ranked pose within the defined binding regions was selected as the initial state for each complex before MD simulation.

Molecular dynamics

Molecular dynamics (MD) simulations were performed using GROMACS 2024.2 (ref. 46). Complex systems comprising xK1-RNT and one drug molecule (CBL, CPT, DOX, or FLU) were prepared through energy minimization, solvation, and ion addition, employing GAFF parameters. The energy minimization step utilized the steepest descent algorithm, while the solvation process employed the TIP3P water model. A NaCl ion concentration of 0.15 M was added to maintain physiological ionic conditions. Equilibration was carried out in two phases. In the first phase, an NVT ensemble (isothermal–isochoric or canonical) was used to stabilize the system's temperature. This was followed by an NPT ensemble to stabilize the system's pressure. A V-rescale thermostat with a time constant of 0.1 ps was employed for temperature control, while pressure regulation was achieved using the Berendsen barostat with a 1 ps time constant and a compressibility of 4.5 × 10⁻⁵ bar⁻¹. After equilibration, a 500 ns production simulation with an integration time step of 2 fs was performed in the NPT ensemble. All bonds involving hydrogen atoms were constrained using the LINCS algorithm. A force-switch modifier was applied to the van der Waals interactions between 1.0 and 1.2 nm. Long-range electrostatic interactions were calculated with the Particle Mesh Ewald (PME) method. The temperature was maintained at 310 K (37 °C) using the velocity-rescaling thermostat (V-rescale) with a coupling constant of 0.1 ps. Pressure coupling was controlled by the Parrinello–Rahman barostat, with isotropic scaling at 1 bar, a time constant of 2.0 ps, and an isothermal compressibility of 4.5 × 10⁻⁵ bar⁻¹. The temperature was set at 310 K to mimic physiological conditions. Each complex system comprised one selected chemotherapeutic drug (CBL, CPT, DOX, FLU) in three different configurations within xK1 obtained from molecular docking. This provided insights into the stability of each complex system under the designed physiological conditions (310 K, 0.15 M NaCl). Each system with the same initial states was executed in triplicate.

Protonated system. Cytosine-rich fragments are reported to protonate with the change in pH.^47–49 They are reported to fold and unfold in acidic and neutral pH environments. Thus, with the change in pH, cytosine-rich fragments can trigger the release of drug molecules to the target cancer cells. G∧C motifs of xK1, having cytosine as a fundamental component, are hypothesized to be able to release drug molecules into the target cancer cells. The G∧C motif for this system was protonated to mimic an acidic environment. N₂ depicted in SI Fig. S1 for each selected motif molecule was protonated after the MD simulation of a complex system under physiological conditions.

The fraction of protonated cytosines for each RNT was determined using the fraction protonated equation (eqn (1)). The reported pK_a of cytosine is approximately 4.5,⁵⁰ and the tumor microenvironment typically exhibits a pH of 5.1. Based on the calculation, approximately 20% of the motifs within each RNT required protonation. Accordingly, 14 motifs of xK1 were protonated as structurally represented in SI Fig. S5. The protonated rings at the termini of xK1 were selected due to their accessibility within the environment. The structural response of each RNT under these acidic conditions was investigated.

The simulation time for this system was set at 100 ns, and the same simulation parameters described above were consistently applied to ensure consistency across all systems.

Post-MD analysis

In the post-simulation phase, the following analyses were performed: (1) radius of gyration (R_g), (2) root mean square deviation (RMSD), (3) detailed interaction analysis, and (4) radial distribution function (RDF) analysis. R_g was employed to evaluate the compactness and stability of each complex system throughout the simulation. It measures how far the atoms of a molecule are spread out from the center of mass. In this study, the radius of gyration (R_g) was analyzed for both xK1 alone and the xK1–drug complex. A peak in the R_g profile of the complex indicates the movement of the drug molecule away from xK1, suggesting structural instability or a transition of the drug from one binding region to another. RMSD analysis was also employed to evaluate the stability of each complex system throughout the simulation. This method evaluates structural changes from the initial configuration to the final state, with an equilibrating RMSD plot indicating structural stability. RDF analysis measures the radial distribution of the selected atoms at a given distance from a set of reference atoms. In this study, the selected atoms were from water molecules, while the reference atoms were taken from the motif molecules within the binding regions, including the drug molecules. RDF analysis was employed to evaluate the hydrophobicity of the complex systems. R_g, RMSD, and RDF analyses were performed using GROMACS.

Detailed interaction analyses were conducted using Schrödinger's Maestro interaction analysis tool to identify additional interactions contributing to the binding of the drug molecule and xK1. This approach offers a comprehensive understanding of the interaction that stabilizes the complex system.

Energetic analysis

MM/PBSA for stability analysis. Energetic analysis was performed to further evaluate the stability of the drug molecule's association with xK1 throughout the simulation period. The Molecular Mechanics Poisson–Boltzmann Surface Area (MM/PBSA) method was employed to estimate the effective binding free energy of each drug–xK1 complex after system stabilization, as determined from the radius of gyration (R_g) analysis following molecular dynamics (MD) simulations. Binding free energies were calculated for every 10th frame of the trajectory to monitor fluctuations and assess the temporal stability of the interaction. In these calculations, xK1 was defined as the receptor and each drug molecule as the ligand.

MM/PBSA calculations were employed here using the gmx_MMPBSA tool. The tool integrates molecular dynamics (MD) trajectories generated with GROMACS into the AMBER-based MM/PBSA.py workflow.⁵¹ Specifically, AmberTools20 (ref. 52) was employed alongside GROMACS for this computational analysis. The workflow of the gmx_MMPBSA tool is typically divided into three key stages: (1) preparation, (2) calculation, and (3) analysis. During the preparation stage, GROMACS topologies are converted into AMBER-compatible topologies using ParmEd.⁵³

The total free energy of each system can be expressed by eqn (2), where E_MM represents the molecular mechanics energy, G_solvation denotes the solvation free energy, and −TS_solute corresponds to the entropic contribution. The overall binding energy is calculated as the difference between the total energy of the complex and the sum of the energies of the isolated RNT and drug molecule, as shown in eqn (3). Accordingly, the total free energy consists of both enthalpic (ΔE_MM + ΔG_solvation) and entropic energy (eqn (4)). The G_solvation is decomposed as polar and non-polar solvation energies, where the polar energy is estimated using the Poisson–Boltzmann model (eqn (5) and (6)).

G = E_MM + G_solvation − TS_solute (2)

ΔG_bind = G_complex − (G_RNT + G_drug) (3)

ΔG_bind = ΔE_MM + ΔG_solvation − TΔS_solute (4)

ΔG_solvation = ΔG_polar + ΔG_nonpolar (5)

G_solvation = ΔG_PBSA + ΔG_nonpolar (6)

ΔG_bind = ΔE_MM + ΔG_solvation (7)

MM/PBSA calculations mainly account for the enthalpic energy without the entropic contribution (eqn (7)). Thus, binding affinities from MM/PBSA calculations are referred to as effective binding energies.⁵⁴ The entropic term can be estimated through separate methods, such as normal-mode analysis; however, it was excluded in this study due to its high computational cost. Furthermore, the calculation of entropic contributions is computationally demanding, requiring extensive sampling and significant memory resources. Limited snapshot sampling often leads to large fluctuations in entropy estimates, and achieving statistically stable predictions would necessitate a prohibitively large number of frames.^55,56 Since the primary objective of this work was to compare relative binding energies over time, the entropic contribution was not included in this stability analysis.

MM/3DRISM for binding energies. The binding energies of the drug molecules against xK1 were quantified using MM/3DRISM analysis. 3DRISM was combined with molecular mechanics calculations using the gmx_MMPBSA tool again. 3D-RISM is a statistical-mechanical solvation model that treats the solvent as a continuum distribution of molecular interaction sites. It inherently accounts for the entropic contributions of solvation, providing a detailed thermodynamic evaluation of binding energies.

For 3DRISM, water was treated as the solvent, while the complex system was considered the solute. Water's solvent site–site correlation functions were computed using RISM theory coupled with the Kovalenko–Hirata (KH) closure approximation.⁵⁷ The solvent-site distribution functions were then obtained on a three-dimensional grid around the complex system by solving the 3DRISM/KH equations using the site–site correlation and potential functions between the solute and the solvent. These were computed on a uniform 3D rectangular grid of 256 × 256 × 256 (ref. 36) within a supercell of 128 × 128 × 128 Å, providing a sufficient buffer of ≈15 Å on each side to fully solvate the structures. These dimensions ensured an adequate solvation space for reliable energetic estimates. A single trajectory protocol (STP) approach was used, where the bound and unbound states of the receptor and ligand were derived from the same trajectory, thus generating zero-bonded energies (ΔE_bonded = 0). As mentioned above, the total binding energy (eqn (7)) is the sum of its molecular mechanical energy (ΔE_MM) and solvation energy (ΔG_3DRISM). Nonbonded energies, specifically van der Waals energy (ΔE_vdW) and electrostatic energy (ΔE_elec), mainly contribute to the ΔE_MM (eqn (8)). The free solvation energy is calculated from the 3DRISM model (eqn (9)).

ΔE_MM = ΔE_bonded + ΔE_nonbonded = ΔE_vdW + ΔE_elec (8)

ΔG_solvation = ΔG_3DRISM = Δε_solv − TΔS_solv (9)

The 3DRISM model inherently accounts for both the solvation energy (Δε_solv) and the solvent entropic contribution (−TΔS_solv). A negative solvent entropic contribution indicates positive entropy, implying a favorable and more disordered system thermodynamically. The overall binding energy comprises nonbonded energies and the energetic components from the 3DRISM calculations, as described in eqn (10).

ΔG_bind = ΔE_vdW + ΔE_elec + Δε_solv − TΔS_solv (10)

Results and discussion

System stability under physiological conditions

Molecular dynamics (MD) simulations with structural and energetic analyses were performed to study the potential of the xK1 type of RNT as a drug vehicle for chemotherapeutic drugs. The molecular behavior, structural stability, and energetic stability were evaluated under physiological conditions to elucidate the potential of xK1 as the carrier. In addition, the structural changes of xK1 in an acidic environment were examined to provide insights into the release mechanism of the drugs from xK1. Global and local docking were employed to identify the initial states of each complex system subjected to MD simulations. Top-ranked binding poses revealed from global docking were utilized, which revealed that all drugs have a high binding affinity at the LysR (SI Fig. S2). Hence, local docking at the LysL and inner channel was performed subsequently. Accordingly, the top-ranked pose at the LysL and inner channel were utilized (SI Fig. S3 and S4), separately.

Shown in Fig. 3 are the top-ranked poses of each drug molecule within the specified binding regions (LysR, LysL, and inner channel). The drugs at both ends are observed to lie directly on the rings and not within the hollow channel. Each drug molecule exhibits diverse conformations at the ends of xK1 (Fig. 3a–j). On the other hand, each drug within the inner channel occupies a unique position and orientation. CBL (Fig. 3k and l) is horizontally inclined in the middle of the inner channel but closer to the rings than to the middle of the hollow channel. CPT (Fig. 3m) is positioned vertically in the center, closer to the rings than the hollow channel. DOX (Fig. 3n) is closer to the LysL end in a vertical position. Similarly, FLU (Fig. 3o) is closer to the LysL end, but it is slightly inclined vertically. These are the initial states of each complex system subjected to MD simulations under physiological conditions.

Fig. 3 Initial configurations of each complex system obtained from molecular docking. (a) Front view of the CBL positioned at LysR. Top views of the same binding site for (b) CBL, (c) CPT, (d) DOX, and (e) FLU. (f) Front view of CBL positioned at LysL, with corresponding top views for (g) CBL, (h) CPT, (i) DOX, and (j) FLU. (k) Top view of CBL initially positioned within the inner channel. Front views of the same binding site for (l) CBL, (m) CPT, (n) DOX, and (o) FLU.

Visual inspection revealed that CBL, initially at the LysR and LysL regions, remained equilibrated in its respective domains, as evidenced by Fig. 4a, b and SI Videos S1, S2. RMSD analysis (Fig. 4a and b) demonstrated small fluctuations (<0.3 nm), which are attributed to the diverse interactions of FLU at various motifs within the LysR and LysL regions. Despite these small fluctuations, the intactness and stability of CBL within the LysR and LysL, respectively, are preserved. In contrast, CBL transitioned and interacted at the LysL end of xK1, as evidenced in Fig. 4c and SI Video S3. The RMSD plot (Fig. 4c) exhibited a significant increase in peak (from 2 nm to 3 nm) around 100 ns during the transition. However, the equilibration of the plot occurred when CBL interacted at LysL. The equilibration of CBL at both ends of xK1 demonstrated a significant interaction between CBL and the termini compared to the inner channel.

Fig. 4 Clustered structures of the CBL–xK1 complex after MD simulations. The clustered MD trajectories correspond to systems where CBL was initially positioned at (a) LysR, (b) LysL, and (c) the inner channel. The structures were clustered in their last 100 ns to indicate the final states of each complex system. The plots below each structure show the corresponding RMSD values calculated specifically for the localized region where the drug molecule was initially placed.

CPT, when initially placed at LysR and LysL, remained in its position throughout the simulation, as evidenced by its MD trajectories (Fig. 5a, b and SI Videos S4, S5). RMSD analysis revealed the stability of the system throughout the simulation as it maintained its initial state configuration until the end of the simulation. Occasional fluctuations (<0.1 nm) are observed, primarily due to a brief transition of CPT from one motif to the next within the initial state. Similar to CBL, CPT transitioned and equilibrated at the LysL end of xK1, as depicted in Fig. 5c and SI Video S6. Nevertheless, the transition was gradual, as evidenced by its RMSD analysis, which indicated a sudden surge in peak (≈1 nm) at approximately 350 ns, and the onset of equilibration was observed afterwards. The equilibration of CPT at the termini of xK1 suggests a stronger interaction at the ends compared to the inner channel.

Fig. 5 Clustered structures of the CPT–xK1 complex after MD simulations. The clustered MD trajectories correspond to systems where CPT was initially positioned at (a) LysR, (b) LysL, and (c) the inner channel. The plots below each structure show the corresponding RMSD values calculated specifically for the localized region where the drug molecule was initially placed.

DOX, when initially placed at LysR, remained in its position throughout the simulation, as evidenced by its MD trajectories (Fig. 6a and SI Videos S7). Similar to CPT, it maintained its initial state configuration throughout the simulation. The RMSD analysis revealed the plot's equilibration throughout the simulation, suggesting a stable binding of DOX to LysR while maintaining its configuration. DOX, initially bound to LysL (Fig. 6b and SI Video S8), remained in the binding site. However, it exhibited varying interactions across different motifs within the binding site, resulting in minimal fluctuations (<0.2 nm) in its RMSD plot. However, the minimal fluctuations suggest the intactness of DOX at the LysL. Analogous to CBL and CPT, DOX transitioned and equilibrated to one end of xK1 (Fig. 6c and SI Video S9). DOX equilibrated at the LysR end of xK1, resulting in an increased peak at 350 ns (≈1 nm), but an onset of equilibration was observed afterward. These findings suggest that DOX is favoured at the ends of the channel rather than in the hollow region of xK1.

Fig. 6 Clustered structures of the DOX–xK1 complex after MD simulations. The clustered MD trajectories correspond to systems where DOX was initially positioned at (a) LysR, (b) LysL, and (c) the inner channel. The plots below each structure show the corresponding RMSD values calculated specifically for the localized region where the drug molecule was initially placed.

The FLU drug molecule remained in the LysR end throughout the simulation (Fig. 7a and SI Video S10). Minor fluctuations are observed in the RMSD analysis, which are due to the shifts of the FLU from one motif to another motif of xK1 within the LysR end. FLU remained at LysL (Fig. 7b and SI Video S11), but it abruptly interacted with the bulk solution for a very short time around 90 ns, which resulted in a very abrupt spike in its RMSD analysis. Frequent small fluctuations (<0.3 nm) are also observed in its RMSD due to frequent shifts of FLU from one motif to another in the LysL region throughout the simulation. Meanwhile, FLU swiftly transitioned from the inner channel to the LysL of xK1, resulting in a sudden peak in the RMSD at approximately 175 ns (Fig. 7c and SI Video S12). Notably, frequent transitions of FLU between different motifs are also observed after this transition. These transitions suggest that FLU possesses multiple metastable states at the termini of xK1. However, the fact that they did not further distance themselves from xK1 indicates strong binding of FLU within the termini of xK1.

Fig. 7 Clustered structures of the FLU–xK1 complex after MD simulations. The clustered MD trajectories correspond to systems where FLU was initially positioned at (a) LysR, (b) LysL, and (c) the inner channel. The plots below each structure show the corresponding RMSD values calculated specifically for the localized region where the drug molecule was initially placed.

Following each MD simulation, comprehensive analyses were conducted, encompassing stability, structural, and energetic evaluations. The radius of gyration (R_g) was utilized to quantify the stability of each complex system. It measures the overall compactness or intactness of the complex system over time.

Additionally, root mean square deviation (RMSD) measurements were performed for the isolated drug, the isolated RNT, and the localized binding region, as previously discussed, to augment the structural analysis. The R_g plots (Fig. 8) reported herein show a comparison of the RNT alone (red line) and the complex system consisting of the drug and the RNT together (gray line). Thus, when the R_g for the complex system exhibits a sudden rise in its peak, it indicates that the drug has moved farther away from the RNT, suggesting the repulsion of the drug and, consequently, instability in the system.

Fig. 8 Radius of gyration of each complex system consisting of CBL (a–c), CPT (d–f), DOX (g–i), and FLU (j–l), initially at the LysR, LysL, and the inner channel of xK1, respectively. The initial states of each complex were obtained from the molecular docking study.

Each drug molecule remained intact in its initial position at LysR throughout the 500 ns simulation (Fig. 8a, d, g and j). Similarly, CBL, CPT, and DOX remained stable at LysL throughout the simulation (Fig. 8a, e and h). FLU also exhibited overall intactness; however, the abrupt elevation observed in its plot (Fig. 8k) can be attributed to a brief interaction of FLU with the bulk solution around 100 ns. This phenomenon was transient, as FLU shortly re-associated with LysL. For systems where the drugs were initially positioned within the inner channel, the trajectories indicated a transition toward the terminal regions of xK1, resulting in a temporary rise in the corresponding R_g values. Despite this transient increase, the R_g later stabilized, indicating that the drugs achieved stable interactions at the terminal regions of xK1 following their initial confinement within the inner channel (Fig. 8c, f, i and l). Overall, the R_g analysis confirmed the stability of the drug molecules at the terminal regions (LysR and LysL). These findings correspond with the observed molecular behavior manifested in their MD trajectories.

The MD simulations were conducted in triplicate. The triplicate runs of the complex system comprising CPT and DOX demonstrated the consistency of their molecular behaviors, suggesting their stability from their initial states, particularly at the LysR and LysL ends, as observed in their R_g plots (SI Fig. S7 and S8). Similarly, these drugs, when initially positioned within the inner channel, consistently exhibited stabilization at the channel ends after transition, further indicating a strong interaction between CPT and DOX at the termini of xK1. This also suggests that once CPT and DOX attain their most stable configurations at the terminal regions of xK1, they are likely to remain bound there under physiological conditions.

In the case of CBL and FLU (SI Fig. S6 and S9), distinct molecular behaviors were observed across different simulation runs. In some trajectories, even when initially positioned at one end of xK1, both drugs eventually shifted toward the opposite end. Before shifting, each drug briefly interacted with the bulk solution but rapidly re-associated with the terminal regions, indicating a stronger affinity for the ends of xK1 compared to the surrounding solvent. This behavior suggests that CBL and FLU possess multiple metastable states at the terminal regions of xK1, leading to frequent transitions between motifs. Nevertheless, these observations still imply strong interactions, as no permanent dissociation or repulsion toward the bulk solvent was detected. Under realistic conditions, the xK1 structures exist as multiple connected fragments during drug delivery, so such movements could represent transitions of the drugs between carriers. The discussed MD trajectories in this paper were focused on the system that exhibited prolonged stability in their R_g and RMSD analyses. Moreover, RMSD analysis for each drug alone and RNT alone was conducted. The RMSD analyses consistently highlighted their structural stability throughout the simulations (SI Fig. S10 and S11). This indicates that RNT, as a carrier, is stable under physiological conditions.

Plots of the binding energies by the MM/PBSA model were also obtained here for a robust stability analysis. The binding energies of each complex when they were stabilized were plotted over time. The binding energy profiles of CBL at LysR and after equilibration at LysL (from the inner channel) revealed occasional occurrences of positive binding energies. At LysR (Fig. 9a), transient spikes were observed at 128 ns (4.61 kcal mol⁻¹), 334 ns (2.94 kcal mol⁻¹), and 435 ns (1.27 kcal mol⁻¹). Similarly, when CBL equilibrated at LysL, minor positive values appeared at 358 ns (0.65 kcal mol⁻¹) and 475 ns (0.87 kcal mol⁻¹). These brief positive fluctuations may correspond to the diverse interactions of CBL within the LysR and LysL regions, where the molecule transitions among motifs. Nevertheless, the predominance of negative binding energies throughout the trajectories indicates that the overall complexes remain stable, reflecting strong interactions across multiple configurations within the binding sites. This observation supports the earlier conclusion that CBL likely exists in multiple metastable binding states. In addition, (Fig. 9a, b and c) the distribution of energy in the plots for CBL was highly concentrated around the mean of each complex (−14.19, −12.6, and −12.16 kcal mol⁻¹), respectively. This indicates that CBL has a favorable binding affinity with the LysR and LysL of xK1.

Fig. 9 MM/PBSA binding energy plots of complexes in their equilibrated phase. The dotted line above the data points denotes the mean binding energy of each complex. Illustrated are the binding energies over time when CBL (a–c) and CPT (d–f) were initially at LysR, LysL, and the inner channel, respectively.

The binding of CPT to both LysR and LysL was found to be stable, as evidenced by the consistently negative binding energy values throughout the simulation after stabilization. In contrast to CBL, CPT exhibits smaller fluctuations, which aligns with the observed behavior of CPT wherein it remained localized at one to two motif molecules throughout the simulation (SI Videos S4 and S5). For the CPT at the LysR complex (Fig. 9d), the distribution of data points was largely concentrated around the mean of −17.33 kcal mol⁻¹, suggesting a stable interaction between CPT and the LysR binding site. Similarly, as shown in Fig. 9e, CPT at LysL exhibited consistently negative binding energies throughout the simulation, indicating a stable interaction. Likewise, the energy distribution was highly concentrated around the mean of −18.40, further supporting the stability of the CPT at LysL. Meanwhile, Fig. 9f illustrates that CPT, after equilibrating at LysR around 380 ns following its transition from the inner channel, also maintained negative binding energies, signifying a strong binding affinity. Overall, the consistent negative binding energies in the plots indicate strong and stable binding of CPT at LysR and LysL.

Among the drugs, DOX possesses the largest molecular size; consequently, its binding energy plots exhibited relatively broad intervals between fluctuations. Despite this, DOX consistently showed negative binding energies across all binding regions (Fig. 10a–c), indicating favorable and stable interactions at both LysR and LysL. Similar to CBL and CPT, the energy distributions were highly concentrated around their respective mean values (−26.41, −18.99, and −20.50 kcal mol⁻¹), indicating the stability of the DOX at LysR and LysL of xK1. In contrast, FLU exhibited a few instances of positive binding energies (0.33 and 2.32 kcal mol⁻¹), as shown in Fig. 10d and e. Nonetheless, the predominance of negative binding energies, like that for CBL, throughout the simulations suggests strong and stable associations of FLU at both LysR and LysL. The energy distributions were also concentrated around the mean values (−12.19, −11.28, and −11.03 kcal mol⁻¹), further confirming complex stability. Overall, all drug molecules demonstrated strong and stable interactions at the LysR and LysL binding regions following equilibration. These findings highlight the suitability of xK1 as a potential drug delivery vehicle for CBL, CPT, DOX, and FLU, given its structural stability and sustained binding behavior under physiological conditions.

Fig. 10 MM/PBSA binding energy plots of complexes in their equilibrated phase. The dotted line above the data points denotes the mean binding energy of each complex. Illustrated are the binding energies over time when DOX (a–c) and FLU (d–f) were initially at LysR, LysL, and the inner channel, respectively.

Further structural interaction analysis was conducted to identify additional interaction forces that contributed to the stability of the complex at the LysR (Fig. 11) and LysL (Fig. 12) regions. The π–π stacking orientations were observed for all complexes, as illustrated in Fig. 11. The π–π stacking interaction is anticipated due to the aromatic moieties of the drugs and the aromatic moieties at the end of xK1. They predominantly exhibited a parallel displaced arrangement. This orientation is known to contribute more substantially to system stability compared to other stacking configurations.⁵⁸ DOX, aside from its π–π stacking interaction, forms a hydrogen bond with the terminal regions (LysR and LysL) of xK1. FLU was revealed to have a π–cation stacking of its nitrogen atom with the aromatic compound of xK1, both at the LysR and LysL regions. The π–cation interactions are recognized for their stronger stabilizing effect compared to π–π stacking interactions.⁵⁹ Thus, the π–cation interaction played a critical role in enhancing the stability of the FLU–xK1 complex.

Fig. 11 Structural interaction analysis of each drug with xK1 when they remained at LysR.

Fig. 12 Structural interaction analysis of each drug with xK1 when they remained at LysL.

Binding energies by MM/3DRISM

Energetic analysis using the Molecular Mechanics/Three-Dimensional Reference Interaction Site Model (MM/3DRISM) was employed to determine the binding affinity of the drug within the xK1 in their stable period. In the MM/3DRISM analysis for each system, the receptor is the RNT and the ligand is the drug molecule. The MM/3DRISM analysis was performed for the system in the equilibrated state. The equilibrated state based on the R_g analysis was clustered as one structure, and it was used as the reference structure for the analysis. The equilibrated state was based on the R_g analysis of each system. Systems with more negative binding energy values indicate stronger binding energy and a more stable complex.

It is observed that the solvation energy (Δε_solv) for each complex is positive due to the nonpolar nature of both RNT and drug molecules. This indicates that each complex system has significant resistance to water. The main forces driving the interaction between xK1 and each drug are nonbonded and nonpolar solvation energies. Hydrophobic interactions play a crucial role, drawing the nonpolar drug molecule toward the nonpolar xK1 in the polar aqueous environment. Once the drug molecule is close enough, nonbonded interactions like van der Waals and electrostatic forces take over, as these are strongly dependent on the proximity of the interacting atoms. This is true and consistent for all complex systems. The entropic contribution (−TΔS_solv) represents the entropic forces arising from the solvent, primarily due to the displacement of water molecules upon the interaction of the drug with the binding site of xK1. A negative entropic energy (−TΔS_solv) corresponds to a positive entropy change, indicating a more thermodynamically favorable and disordered system. Therefore, a more negative entropic contribution reflects greater water displacement upon drug binding, signifying stronger binding affinity, enhanced stability, and a favorable contribution to the overall binding energy.

The MM/3DRISM analysis of all complex systems (Table 1) revealed strong binding energies of each drug molecule to the LysR and LysL sites. This aligns with the earlier observed molecular behavior and the obtained structural analysis. Most binding energies are more negative than the reported binding energy of a small drug molecule to xK1 (−5.81 kcal mol⁻¹).³⁶ When CBL and FLU transitioned from the inner channel to the LysL site, their resulting negative binding energies (−5.21 and −4.66 kcal mol⁻¹) were weaker than those reported in a previous study. Nevertheless, since these energies are not significantly beyond the threshold, the negative values still indicate strong binding affinity. DOX consistently showed stronger binding energies (−24.10 kcal mol⁻¹ and −16.25 kcal mol⁻¹) at LysR and LysL than the rest of the drugs, which is likely due to its larger size than the rest of the molecules. The prominent aromatic rings at the ends of xK1 and each drug molecule, which facilitate π–π stacking interactions, favorably contribute to the van der Waals energy, which significantly contributes to the overall binding energy of each complex. Additionally, the entropic forces exhibited favorable contributions to the overall binding energy. The overall results suggest that these drugs have favorable binding energies with xK1, particularly at the ends.

Table 1 MM/3DRISM analysis of each complex system in kcal mol⁻¹ units^a

	Initial state	Final state	ΔEvdW	ΔEelec	Δεsolv	−TΔSsolv	ΔGbind
a The selected system for the calculations is based on the stability observed from the Rg plots. The most stable plot in Rg was selected for further MM/3DRISM calculations, while the rest were calculated from their first simulation run.b Calculations applied to the system from run 3 among the triplicate simulation runs.c Calculations employed in the system from run 2 among the triplicate simulation runs.
CBL	LysR	LysR^b	−26.61	−23.81	47.97	−6.48	−8.94
	LysL	LysL^b	−26.52	12.87	9.22	−7.13	−11.56
	Inner	LysL	−16.93	−16.06	37.19	−9.40	−5.21
CPT	LysR	LysR	−32.77	7.18	25.14	−8.39	−8.85
	LysL	LysL	−30.03	−15.28	40.54	−7.96	−12.73
	Inner	LysR	−32.40	−8.33	41.32	−10.34	−9.74
DOX	LysR	LysR	−41.32	−19.28	48.51	−12.01	−24.10
	LysL	LysL	−38.50	−17.63	47.68	−7.80	−16.25
	Inner	LysL	−40.37	−17.35	57.67	−10.11	−10.16
FLU	LysR	LysR^c	−21.48	−36.16	56.70	−6.20	−7.15
	LysL	LysL^b	−21.54	6.74	14.00	−5.12	−5.92
	Inner	LysL^c	−20.17	−23.61	41.57	−2.45	−4.66

---

As discussed, hydrophobic interactions significantly contribute to the binding of the drug molecules to xK1, due to their nonpolar nature. Radial distribution function (RDF) analysis was employed to quantify the hydrophobicity of the system, wherein the radial distribution of water molecules from the complex system was analyzed.

Fig. 13 illustrates the radial distribution of water molecules at a distance from the complex system. Each line in the plot represents a binding region of xK1, such as LysR, LysL, and the inner channel. In the analysis, the drug molecules together with the motifs from each binding region were selected as the reference atoms. The water molecules, on the other hand, were selected as the atoms to be evaluated. Fig. 13 demonstrates that water molecules consistently have a higher distribution at a farther distance of around 4 nm, indicating the repulsion of water from the complex. This quantifies and validates the hydrophobicity of each complex.

Fig. 13 RDF analysis of water molecules at a distance from the complex system.

The reported structural and energetic analyses demonstrate that xK1 is a potential vehicle for chemotherapeutic drug molecules. All drugs, when initially at LysR and LysL, remain in their respective places. All systems also revealed that, starting from the inner channel, each drug molecule transitioned and eventually equilibrated at either end (LysR or LysL). This indicates that they are more successful on either end than in the inner channel of xK1. The equilibration of the drug molecules at LysR and LysL further supports findings from previous research, which suggest that xK1, in particular, lacks a significant potential energy barrier at either end.³⁶ The potential energetic barrier is the energy that the molecule must overcome or pass through in order to be released from the inner channel of the RNT. This explains why the drugs were drawn by strong binding forces to the ends rather than remaining in the inner channel. In addition, most chemotherapeutic drugs contain aromatic rings, a fundamental component of RNTs that is more prominent at the ends. This characteristic contributes to their localization at either end due to π–π stacking interactions. Some papers reported that π–π stacking interactions increase the stability of the chemotherapeutic drugs with the carrier due to their aromatic compounds.^60–62

It was discussed in a previous study that there are two potential loading mechanisms for the three types of RNT (K1, xK1, and iEt-xK1): (1) drugs will attach to either end after the stabilization of RNT, and (2) drugs will be encapsulated upon the self-assembly process.³⁶ A similar loading mechanism is proposed in this paper, with the former mechanism being potentially more applicable for loading chemotherapeutic drug molecules in xK1-type RNTs. This is likely due to the prominent aromatic rings in xK1. As to the need to protect drug molecules from environmental exposure, particularly in noncancerous environments (e.g., normal and healthy cells), the attachment of the drug vehicle at either end does not compromise its functionality in delivering the drug to the target cells. The strong binding energies and the observed structural stability of xK1 enhance its ability to shield the drug molecules from noncancerous cells, due to their strong associations (−4.66 kcal mol⁻¹ to −24.10 kcal mol⁻¹) with the vehicle. Several studies have reported the attachment of chemotherapeutic drugs to the surfaces of drug carriers, such as carbon nanotubes, graphene oxide, quantum dots, and other vehicles, which exhibit stability under physiological conditions.^9,62,63 Furthermore, rosette nanotubes (RNTs) are typically longer in reality (≈100 nm) than in the simulations, suggesting that when multiple drug molecules are confined within the inner channel, they may remain entrapped for longer durations before transitioning to the ends due to the extended length of the nanotube.

Releasing mechanism in acidic environments

Regarding drug release, several studies have explored the release behavior of nanocarriers under acidic conditions through molecular simulations.^64–69 A common approach involves protonating the nanocarrier to examine its structural behavior at low pH levels. A carbon nanotube conjugated with DNA fragments mainly composed of cytosine was employed in one study.⁴⁹ The researchers studied the structural behavior of the nanocarrier in a protonated environment to better understand the release mechanism at low pH levels.⁴⁹ They aimed to determine the behavior of the nanocarrier system when exposed to acidic environments. In this study, the G∧C motifs within each RNT were selectively protonated on the cytosine residues, following the approach described by the mentioned study.

Earlier, each drug exhibited stability when it remained at the LysR site of xK1. In this study, stable complex systems, where drugs remained at LysR, were protonated to elucidate the release mechanism. From the visual inspection of the MD trajectories (Fig. 14a–d and SI representative Video S13), xK1 began to disassemble within the first 10 ns, as evidenced by the increasing R_g (Fig. 15a–d) during this period. This behavior indicates the successful release of the drug from the RNT structure. Although some drug molecules appeared to interact with residual aggregates from the disrupted RNTs, the overall breakdown still reflects a release mechanism.

Fig. 14 Clustered MD trajectories of the complex systems at pH 5.2, consisting of (a) CBL, (b) CPT, (c) DOX, and (d) FLU at the LysR site. The structures shown correspond to the clusters obtained from the last 10 ns of the simulation, with the first representative structure from each cluster displayed. The initial states used here were derived from the stable conformations of each complex under physiological conditions.

Fig. 15 Radius of gyration analysis of the MD simulations of (a) CBL, (b) CPT, (c) DOX, and (d) FLU under acidic conditions.

The inherent biocompatibility of the RNT motifs supports their nontoxic persistence within biological environments, as demonstrated in prior biomedical applications. For instance, RNTs were utilized for tissue engineering applications as a scaffold for neural regeneration and nerve repair.³¹ One study showed that RNT motifs are internalized via temperature-dependent endocytosis, confirming that endocytic pathways are primarily responsible for their cellular uptake.³⁰ This suggests that, following RNT disassembly, target cells can efficiently absorb both the drug and the motif molecules through natural cellular mechanisms. Thus, an endocytosis mechanism is anticipated to occur within the target cancer cells.

Conclusions

This study comprehensively investigated the structural and energetic stability of the xK1 type of RNT as a potential delivery vehicle for chemotherapeutic drugs under physiological conditions, with the purpose of reducing the adverse effects associated with these drugs. This study also investigated the structural response of RNT in acidic environments to elucidate the release mechanism. Molecular dynamics simulations, aided by robust energetic analysis, were conducted to assess the viability of xK1 as a drug carrier and to elucidate factors contributing to its stability in physiological environments. Results from R_g, RMSD, and MM/3DRISM analyses revealed that the xK1 type of RNT exhibits both structural and energetic stability with CBL, CPT, DOX, and FLU as the drug molecules, under physiological conditions, with binding energies ranging from −4.66 to −24.10 kcal mol⁻¹. Visual inspection further demonstrated that these drugs preferentially localize at the terminal regions of the xK1 structure, attributed to the prominent aromatic rings at both ends, which facilitate π–π stacking and π–cation interactions. On the other hand, to elucidate the mechanism of drug release within target cells, molecular dynamics simulations under acidic conditions were performed. Given that the tumor microenvironment is often characterized by acidic pH, drug release was investigated by protonating the stable RNT (pH = 5.1). Protonating the RNT at low pH levels confirmed the successful release of each drug following structural disassembly within the first 10 ns, supporting its potential as a responsive drug-delivery platform. These results offer a coherent basis for the experimental validation and future optimization of RNT-based nanocarriers.

Author contributions

The manuscript was written through the contributions from all the authors. All authors have given approval to the final version of the manuscript. Hanah Nasifa M. Ali: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, resources, validation, visualization, writing-original draft, and writing-review & editing. Arthur A. Gonzales III: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, software, project administration, resources, software, supervision, validation, visualization, writing-original draft, and writing-review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data presented in this study are available in the article and the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5ra04657b.

Acknowledgements

The authors acknowledge the support provided by the Department of Science and Technology – Engineering Research and Development for Technology (DOST-ERDT), HPC Cluster Team of the College of Engineering, the University of the Philippines Diliman, the Office of the Vice-Chancellor for Research and Development (UPD-OVCRD), and the UPD – Engineering Research and Development Foundation, Inc. (UPD-ERDFI), which have been instrumental in the completion of this research.

References

1. M. F. Manalo, S. Ng and S. Ozdemir, et al., Quality of life and psychological distress of patients with advanced cancer in the Philippines, Qual. Life Res., 2023, 32(8), 2271–2279,  DOI:10.1007/s11136-023-03389-y.
2. A. Chatterjee, R. Khan and T. Mukherjee, et al., Harnessing bacterial metabolites for enhanced cancer chemotherapy: unveiling unique therapeutic potentials, Arch. Microbiol., 2024, 206(11), 449,  DOI:10.1007/s00203-024-04179-x.
3. Y. Sun, Q. Li and Y. Huang, et al., Natural products for enhancing the sensitivity or decreasing the adverse effects of anticancer drugs through regulating the redox balance, Chin. Med., 2024, 19(1), 110,  DOI:10.1186/s13020-024-00982-2.
4. F. Yu, Y. Li and J. Zou, et al., The Chinese herb Xiaoaiping protects against breast cancer chemotherapy-induced alopecia and other side effects: a randomized controlled trial, J. Int. Med. Res., 2019, 47(6), 2607–2614,  DOI:10.1177/0300060519842781.
5. F. Rosenbaek, H. S. Holm, J. v. B. Hjelmborg, M. Ewertz and J. D. Jensen, Effect of cryotherapy on dose of adjuvant paclitaxel in early-stage breast cancer, Support. Care Cancer, 2020, 28(8), 3763–3769,  DOI:10.1007/s00520-019-05196-z.
6. H. Ma, Z. Pan, B. Lai, C. Zan and H. Liu, Recent Research Advances in Nano-Based Drug Delivery Systems for Local Anesthetics, Drug Des., Dev. Ther., 2023, 17, 2639–2655,  DOI:10.2147/DDDT.S417051.
7. S. Prakash, Nano-based drug delivery system for therapeutics: a comprehensive review, Biomedical Physics & Engineering Express, 2023, 9(5), 052002,  DOI:10.1088/2057-1976/acedb2.
8. P. K. Gupta, Nano-Based Drug Delivery Systems, in Nanotoxicology in Nanobiomedicine, ed. P. Gupta, Springer International Publishing, 2023, pp. 91–110,  DOI:10.1007/978-3-031-24287-8_6.
9. M. I. Khan, I. Nadeem, A. Majid and M. Shakil, Adsorption mechanism of Palbociclib anticancer drug on two different functionalized nanotubes as a drug delivery vehicle: A first principle's study, Appl. Surf. Sci., 2021, 546, 149129,  DOI:10.1016/j.apsusc.2021.149129.
10. H. Ijaz, A. Mahmood and M. M. Abdel-Daim, et al., Review on carbon nanotubes (CNTs) and their chemical and physical characteristics, with particular emphasis on potential applications in biomedicine, Inorg. Chem. Commun., 2023, 155, 111020,  DOI:10.1016/j.inoche.2023.111020.
11. A. Kordzadeh, M. Zarif and S. Amjad-Iranagh, Molecular dynamics insight of interaction between the functionalized-carbon nanotube and cancerous cell membrane in doxorubicin delivery, Computer Methods and Programs in Biomedicine, 2023, 230, 107332,  DOI:10.1016/j.cmpb.2022.107332.
12. O. E. Glukhova, Liposome Drug Delivery System across Endothelial Plasma Membrane: Role of Distance between Endothelial Cells and Blood Flow Rate, Molecules, 2020, 25(8), 1875,  DOI:10.3390/molecules25081875.
13. J. Kurczewska, Chitosan-Based Nanoparticles with Optimized Parameters for Targeted Delivery of a Specific Anticancer Drug—A Comprehensive Review, Pharmaceutics, 2023, 15(2), 503.
14. D. Agudelo, G. Bérubé and H. A. Tajmir-Riahi, An overview on the delivery of antitumor drug doxorubicin by carrier proteins, Int. J. Biol. Macromol., 2016, 88, 354–360,  DOI:10.1016/j.ijbiomac.2016.03.060.
15. N. Farhadian, M. S. Kazemi, F. Moosavi Baigi and M. Khalaj, Molecular dynamics simulation of drug delivery across the cell membrane by applying gold nanoparticle carrier: Flutamide as hydrophobic and glutathione as hydrophilic drugs as the case studies, J. Mol. Graphics Modell., 2022, 116, 108271,  DOI:10.1016/j.jmgm.2022.108271.
16. D. Lombardo, M. A. Kiselev and M. T. Caccamo, Smart Nanoparticles for Drug Delivery Application: Development of Versatile Nanocarrier Platforms in Biotechnology and Nanomedicine, J. Nanomater., 2019, 2019, e3702518,  DOI:10.1155/2019/3702518.
17. S. A. Dehdast and O. Pourdakan, Quantum Dots Nanocomposites as Drug Delivery Vehicle, in Quantum Dots Based Nanocomposites: Design, Fabrication and Emerging Applications, ed. S. Thomas, P. Das and S. Ganguly, Springer, Nature Switzerland, 2024, pp. 323–336,  DOI:10.1007/978-3-031-54779-9_16.
18. A. Mikhraliieva, O. Tkachenko and R. Freire, et al., Carbon Nanodots with Solvatochromic Photoluminescence for the Electrochemical Determination of Estrogenic Steroids, ACS Appl. Nano Mater., 2022, 5(8), 10962–10972,  DOI:10.1021/acsanm.2c02219.
19. B. M. Sahoo, D. Shree, B. K. Banik, P. Borah, A. Jain and M. K. Mahapatra, Chapter 2 - Nano-based drug delivery of anticancer agents, in Green Approaches in Medicinal Chemistry for Sustainable Drug Design, ed. B. K. Banik, Advances Green and Sustainable Chemistry, Elsevier, 2nd edn, 2024, vol. 2, pp. 43–60,  DOI:10.1016/B978-0-443-16164-3.00002-9.
20. S. Tanreh, M. Rezvani and M. Darvish Ganji, Molecular simulation investigations on interaction properties of the teriflunomide–chitosan complex in aqueous solution, J. Phys. Chem. Solids, 2023, 174, 111171,  DOI:10.1016/j.jpcs.2022.111171.
21. M. R. Fadaei, M. S. Fadaei and A. E. Kheirieh, et al., Overview of dendrimers as promising drug delivery systems with insight into anticancer and anti-microbial applications, Int. J. Pharm.:X, 2025, 10, 100390,  DOI:10.1016/j.ijpx.2025.100390.
22. W. Chen, B. Wang and S. Liang, et al., Renal clearance of graphene oxide: glomerular filtration or tubular secretion and selective kidney injury association with its lateral dimension, J. Nanobiotechnol., 2023, 21(1), 51,  DOI:10.1186/s12951-023-01781-x.
23. S. Liang, B. Wang and X. Li, et al., In vivo pharmacokinetics, transfer and clearance study of graphene oxide by La/Ce dual elemental labelling method, NanoImpact, 2020, 17, 100213,  DOI:10.1016/j.impact.2020.100213.
24. H. Fenniri, P. Mathivanan and K. L. Vidale, et al., Helical Rosette Nanotubes: Design, Self-Assembly, and Characterization, J. Am. Chem. Soc., 2001, 123(16), 3854–3855,  DOI:10.1021/ja005886l.
25. R. L. Beingessner, Y. Fan and H. Fenniri, Molecular and supramolecular chemistry of rosette nanotubes, RSC Adv., 2016, 6(79), 75820–75838,  10.1039/C6RA16315G.
26. H. Fenniri, B. L. Deng, A. E. Ribbe, K. Hallenga, J. Jacob and P. Thiyagarajan, Entropically driven self-assembly of multichannel rosette nanotubes, Proc. Natl. Acad. Sci. U. S. A., 2002, 99(suppl 2), 6487–6492,  DOI:10.1073/pnas.032527099.
27. T. Yamazaki, H. Fenniri and A. Kovalenko, Structural Water Drives Self-assembly of Organic Rosette Nanotubes and Holds Host Atoms in the Channel, ChemPhysChem, 2010, 11(2), 361–367,  DOI:10.1002/cphc.200900324.
28. Y. Chen, S. Song, Z. Yan, H. Fenirri and T. J. Webster, Self-assembled rosette nanotubes encapsulate and slowly release dexamethasone, Int. J. Nanomed., 2011, 6, 1035–1044,  DOI:10.2147/IJN.S18755.
29. S. Song, Z. Yan, H. Fenniri and T. J. Webster, Self-assembled rosette nanotubes for incorporating hydrophobic drugs in physiological environments, Int. J. Nanomed., 2011, 6, 101–107,  DOI:10.2147/IJN.S11957.
30. U. Ho, M. El-Bakkari and A. Alshamsan, et al., Delivery of siRNA using cationic rosette nanotubes for gene silencing, Biomater. Sci., 2023, 11(21), 7169–7178,  10.1039/D3BM01115A.
31. M. L. Puzan, B. Legesse, R. A. Koppes, H. Fenniri and A. N. Koppes, Bioactive Organic Rosette Nanotubes Support Sensory Neurite Outgrowth, ACS Biomater. Sci. Eng., 2018, 4(5), 1630–1640,  DOI:10.1021/acsbiomaterials.8b00326.
32. Y. Chen, B. Bilgen and R. Pareta, et al., Self-Assembled Rosette Nanotube/Hydrogel Composites for Cartilage Tissue Engineering, Tissue Eng., Part C, 2010, 16, 1233–1243,  DOI:10.1089/ten.TEC.2009.0400.
33. P. Wolski and T. Panczyk, Insight Into Interfacial Behaviors between Doxorubicin and Zwitterion/PAMAM/CQD Hybrid Nanocarrier. A Molecular Dynamics Simulations Study, J. Phys. Chem. B, 2024, 128(37), 8946–8955,  DOI:10.1021/acs.jpcb.4c05572.
34. M. M. Alkriz and D. Joujeh, In silico Design of a Quadruplex RNA Aptamer Against Carcinoembryonic Antigen, and Evaluating its Potential as a Drug Delivery Vehicle, Asian Pac. J. Cancer Biol., 2024, 9(3), 307–314,  DOI:10.31557/apjcb.2024.9.3.307-314.
35. Schrödinger Release 2021-1, Maestro, Schrödinger, LLC, New York, NY, 2021.
36. T. Yamazaki and H. Fenniri, Encapsulation of ferrocene by self-assembled rosette nanotubes: An investigation using statistical mechanical theory of molecular liquids, J. Mol. Liq., 2016, 217, 70–74,  DOI:10.1016/j.molliq.2015.11.055.
37. PubChem, Structure Search, https://pubchem.ncbi.nlm.nih.gov/docs/structure-search, accessed 20 December 2024.
38. H. E. Pence and A. Williams, ChemSpider: An Online Chemical Information Resource, J. Chem. Educ., 2010, 87(11), 1123–1124,  DOI:10.1021/ed100697w.
39. A. W. Sousa da Silva and W. F. Vranken, ACPYPE - AnteChamber PYthon Parser interfacE, BMC Res. Notes, 2012, 5(1), 367,  DOI:10.1186/1756-0500-5-367.
40. M. Vögele, J. Köfinger and G. Hummer, Molecular dynamics simulations of carbon nanotube porins in lipid bilayers, Faraday Discuss., 2018, 209, 341–358,  10.1039/C8FD00011E.
41. R. Esfandiarpour, F. Badalkhani-Khamseh and N. L. Hadipour, Exploration of phosphorene as doxorubicin nanocarrier: An atomistic view from DFT calculations and MD simulations, Colloids Surf., B, 2022, 215, 112513,  DOI:10.1016/j.colsurfb.2022.112513.
42. V. Jain, P. K. Maiti and P. V. Bharatam, Atomic level insights into realistic molecular models of dendrimer-drug complexes through MD simulations, J. Chem. Phys., 2016, 145(12), 124902,  DOI:10.1063/1.4962582.
43. H. Lei, X. Wang and C. Wu, Early stage intercalation of doxorubicin to DNA fragments observed in molecular dynamics binding simulations, J. Mol. Graphics Modell., 2012, 38, 279–289,  DOI:10.1016/j.jmgm.2012.05.006.
44. J. šponer, X. Cang and T. E. Cheatham, Molecular dynamics simulations of G-DNA and perspectives on the simulation of nucleic acid structures, Methods, 2012, 57(1), 25–39,  DOI:10.1016/j.ymeth.2012.04.005.
45. G. M. Morris, R. Huey and W. Lindstrom, et al., AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility, J. Comput. Chem., 2009, 30(16), 2785–2791.
46. M. Abraham, A. Alekseenko and V. Basov, et al., GROMACS 2024.2 Manual, 2024,  DOI:10.5281/zenodo.11148638.
47. P. Wolski, K. Nieszporek and T. Panczyk, Pegylated and folic acid functionalized carbon nanotubes as pH controlled carriers of doxorubicin. Molecular dynamics analysis of the stability and drug release mechanism, Phys. Chem. Chem. Phys., 2017, 19(13), 9300–9312,  10.1039/C7CP00702G.
48. P. Wolski, K. Nieszporek and T. Panczyk, Cytosine-Rich DNA Fragments Covalently Bound to Carbon Nanotube as Factors Triggering Doxorubicin Release at Acidic pH. A Molecular Dynamics Study, Int. J. Mol. Sci., 2021, 22(16), 8466,  DOI:10.3390/ijms22168466.
49. P. Wolski, K. Nieszporek and T. Panczyk, Carbon Nanotubes and Short Cytosine-Rich Telomeric DNA Oligomeres as Platforms for Controlled Release of Doxorubicin—A Molecular Dynamics Study, Int. J. Mol. Sci., 2020, 21(10), 3619,  DOI:10.3390/ijms21103619.
50. A. Furukawa, E. Walinda, K. Arita and K. Sugase, Structural dynamics of double-stranded DNA with epigenome modification, Nucleic Acids Res., 2021, 49(2), 1152–1162,  DOI:10.1093/nar/gkaa1210.
51. B. R. Miller, T. D. McGee, J. M. Swails, N. Homeyer, H. Gohlke and A. E. Roitberg, MMPBSA.py: An Efficient Program for End-State Free Energy Calculations, J. Chem. Theory Comput., 2012, 8(9), 3314–3321,  DOI:10.1021/ct300418h.
52. D. A. Case, H. M. Aktulga and K. Belfon, et al., Amber 2022, University of California, San Francisco, 2022.
53. J. Swails, C. Hernandez, D. L. Mobley, H. Nguyen, L. P. Wang and P. Janowski, ParmEd, 2010, https://github.com/ParmEd/ParmEd.
54. N. Homeyer and H. Gohlke, Free Energy Calculations by the Molecular Mechanics Poisson−Boltzmann Surface Area Method, Mol. Inf., 2012, 31(2), 114–122,  DOI:10.1002/minf.201100135.
55. S. Genheden and U. Ryde, The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities, Expert Opin. Drug Discovery, 2015, 10(5), 449–461,  DOI:10.1517/17460441.2015.1032936.
56. T. Hou, J. Wang, Y. Li and W. Wang, Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations, J. Chem. Inf. Model., 2011, 51(1), 69–82,  DOI:10.1021/ci100275a.
57. A. Kovalenko and F. Hirata, Three-dimensional density profiles of water in contact with a solute of arbitrary shape: a RISM approach, Chem. Phys. Lett., 1998, 290(1), 237–244,  DOI:10.1016/S0009-2614(98)00471-0.
58. G. B. McGaughey, M. Gagné and A. K. Rappé, pi-Stacking interactions. Alive and well in proteins, J. Biol. Chem., 1998, 273(25), 15458–15463,  DOI:10.1074/jbc.273.25.15458.
59. A. S. Mahadevi and G. N. Sastry, Cation−π Interaction: Its Role and Relevance in Chemistry, Biology, and Material Science, Chem. Rev., 2013, 113(3), 2100–2138,  DOI:10.1021/cr300222d.
60. K. Bhushan A, D. Ammunje, S. Kunjiappan, J. Anbu, M. Devi and P. Pavadai, An in-Silico approach to evaluate the binding efficacy and stability profile of MWCNT entangled rutin for breast cancer treatment, Nanomed. J., 2023, 10(3), 210–215,  DOI:10.22038/NMJ.2023.70663.1756.
61. W. Dai, P. Jin and X. Li, et al., A carrier-free nano-drug assembled via π–π stacking interaction for the treatment of osteoarthritis, Biomed. Pharmacother., 2023, 164, 114881,  DOI:10.1016/j.biopha.2023.114881.
62. X. Tang, S. Zhou, X. Tao, J. Wang, F. Wang and Y. Liang, Targeted delivery of docetaxel via Pi-Pi stacking stabilized dendritic polymeric micelles for enhanced therapy of liver cancer, Mater. Sci. Eng., C, 2017, 75, 1042–1048,  DOI:10.1016/j.msec.2017.02.098.
63. A. Zaboli, H. Raissi, H. Hashemzadeh and F. Farzad, Graphene Oxide Hosting a pH-Sensitive Prodrug: An In Silico Investigation of Graphene Oxide-Based Nanovehicle toward Cancer Therapy, ACS Appl. Bio Mater., 2023, 6(7), 2826–2836,  DOI:10.1021/acsabm.3c00276.
64. N. S. Kulkarni, Y. Guererro, N. Gupta, A. Muth and V. Gupta, Exploring potential of quantum dots as dual modality for cancer therapy and diagnosis, J. Drug Delivery Sci. Technol., 2019, 49, 352–364,  DOI:10.1016/j.jddst.2018.12.010.
65. J. Lasham, A. Djurabekova, V. Zickermann, J. Vonck and V. Sharma, Role of Protonation States in the Stability of Molecular Dynamics Simulations of High-Resolution Membrane Protein Structures, J. Phys. Chem. B, 2024, 128(10), 2304–2316,  DOI:10.1021/acs.jpcb.3c07421.
66. J. Li, S. Ying and H. Ren, et al., Molecular dynamics study on the encapsulation and release of anti-cancer drug doxorubicin by chitosan, Int. J. Pharm., 2020, 580, 119241,  DOI:10.1016/j.ijpharm.2020.119241.
67. R. Maleki, H. H. Afrouzi, M. Hosseini, D. Toghraie, A. Piranfar and S. Rostami, pH-sensitive loading/releasing of doxorubicin using single-walled carbon nanotube and multi-walled carbon nanotube: A molecular dynamics study, Computer Methods and Programs in Biomedicine, 2020, 186, 105210,  DOI:10.1016/j.cmpb.2019.105210.
68. R. Maleki, H. H. Afrouzi, M. Hosseini, D. Toghraie and S. Rostami, Molecular dynamics simulation of Doxorubicin loading with N-isopropyl acrylamide carbon nanotube in a drug delivery system, Computer Methods and Programs in Biomedicine, 2020, 184, 105303,  DOI:10.1016/j.cmpb.2019.105303.
69. G. Zhang, X. Zhang, K. Yue and W. Zhong, Mechanistic study of enhanced drug release in mixed pH-responsive peptide-loaded liposomes, J. Biomol. Struct. Dyn., 2025, 1–15,  DOI:10.1080/07391102.2025.2481581.

---