Effects of doping and functionalization on black phosphorus nanosheets for enhanced quercetin delivery

Abstract

Various nanoparticle-based drug delivery platforms have been developed to address the challenges of rapid degradation and immune clearance in cancer treatment while improving tumor-targeting efficiency. Among these, black phosphorus (BPH) has emerged as a promising multifunctional nanodrug platform. This study investigates the behavior of black phosphorus nanosheets doped with nickel, copper, and zinc atoms in the delivery of quercetin through molecular dynamics simulations. Moreover, the combined effects of doping and functionalization on drug adsorption and delivery processes were examined. The results suggested that the tendency for drug adsorption on the nanosheet surface primarily arises from π–π interactions. On the other hand, functionalization of the doped nanosheets significantly improves the adsorption of quercetin molecules. Specifically, we found that the average adsorption energy for pristine black phosphorus nanosheets is −203 kJ mol⁻¹, while functionalized black phosphorus nanosheets doped with nickel atoms show a significantly improved average adsorption energy of −257 kJ mol⁻¹. These findings, along with analysis of the radial distribution function, solvent-accessible surface area, and other molecular dynamics parameters, provide valuable insights, contributing to a deeper understanding of how modifications to black phosphorus nanosheets can optimize drug delivery systems for effective cancer therapy.

---

Introduction

In the last few decades, nanotechnology has played a significant role in medicine and pharmaceuticals and is developing. In this regard, researchers have focused on developing drug-delivery systems based on nanostructures.^1–3 A wide range of nanocarriers, including graphene sheets, carbon nanotubes, dichalcogenides of transition metals, and metal organic frameworks, have been used to load and transport drugs to different parts of the body.^4–8 Among these nanocarriers, two-dimensional materials have received much attention due to their unique structure. For example, graphene and boron nitride have been widely exploited in medicine.^9–12

One of the new two-dimensional nanomaterials discovered in recent years is phosphorene.¹³ It features a layered structure with weak van der Waals interactions between its layers, which distinguishes it from graphene and other two-dimensional materials.^14,15 In addition, its wrinkled shape allows it to offer a higher surface-to-volume ratio compared to graphene and boron nitride.¹⁶ Black phosphorus has excellent biocompatibility and is converted into phosphite and phosphate, which are non-toxic and easily eliminated from the body.^17–19 This nanostructure is specifically designed to pass through biological membranes, selectively interact with patient cells, and enhance the retention time of drugs in the body. These unique properties make it a suitable option for medical research.¹⁹

One effective approach to improving the performance of nanocarriers is through surface modification, which can involve doping with various atoms or functionalizing the surface.^20,21 Doping can significantly enhance the stability of BPHs in biological environments, helping to prevent their rapid degradation. Additionally, this technique can facilitate controlled drug release mechanisms, allowing for the targeted release of therapeutic agents in response to specific stimuli. This approach enhances therapeutic outcomes and overall treatment efficacy.^22,23 Modifying BPHs also improves their ability to target tumor cells specifically. This targeted approach minimizes side effects on healthy tissues, making treatments safer and more effective. By tailoring the properties of these nanocarriers, we can achieve better therapeutic results while reducing unwanted impacts on surrounding healthy cells.²⁴

Zhang et al.²⁵ used black phosphorus nanosheets to deliver doxorubicin in chemotherapy, and the results of their in vitro and in vivo experiments showed that this drug delivery system is biocompatible and has high antitumor activity. In a theoretical study, Yin et al.²⁶ explored the potential of polydopamine-modified BPHs for treating stroke. Their experimental results demonstrated a strong photothermal response, stability, biocompatibility, and significant therapeutic efficacy of this drug delivery system. These findings suggest that surface modification of BPHs holds promise for medical applications. Additionally, Mahboobi et al.¹⁵ examined the impact of doping BPHs with sodium, calcium, and iron ions on the delivery of mercaptopurine. Their density functional theory (DFT) calculations indicated that the doped nanosheets exhibit enhanced drug delivery capabilities. Together, these studies highlight the potential of modified BPHs in improving drug delivery systems for various medical conditions.

In the current study, we developed a black phosphorus-based drug delivery system aimed at delivering the anticancer drug quercetin through molecular dynamics simulations. We investigated the effects of simultaneously doping and functionalizing BPHs with nickel, copper, and zinc atoms, along with a L-lysine functional group, as an innovative approach to drug delivery. To ensure the accuracy of our findings, we employed various analytical methods, including metadynamics calculations, radial distribution functions, and number of contacts assessments. The overall results suggest that the combined strategy of doping and functionalization of 2D nanosheets significantly enhances their potential for effective drug delivery.

Method

Preparing systems

The design of black phosphorus structure, comprising 200 phosphorus atoms, was created using GaussView software.²⁷ To reduce the effects of unsaturation, hydrogen atoms were added to the edges. The three-dimensional structure of the drug quercetin molecule (PubChem CID: 5280343) was obtained from the SwissParam web server.²⁸ Phosphorene nanosheets were individually doped with Ni, Cu, and Zn atoms.^29–32

Ni, Cu, and Zn were selected as dopants because of their known ability to alter the electronic and catalytic properties of two-dimensional (2D) materials. As first-row transition metals, these elements can effectively adjust the band gap, improve charge transport, and enhance the stability of phosphorene.^33–39 A doping level of approximately 0.5% (two metal atoms for every 200 phosphorus atoms) was chosen to represent a low, experimentally relevant concentration. This level allows for the examination of isolated dopant effects without significantly distorting the phosphorene lattice. Consequently, this approach maintains the intrinsic properties of the nanosheet while facilitating meaningful interactions between the dopants and the substrate.

In selecting the doping sites for Ni, Cu, and Zn atoms on the phosphorene nanosheet, the central region of the nanosheet was chosen based on structural stability and electronic effect optimization considerations. The central site of the nanosheet typically provides a more uniform and less strained environment than edge sites, ensuring minimal structural distortion upon doping. This choice also facilitates strong interaction between the dopant atoms and the phosphorene lattice, which can significantly modulate the electronic and catalytic properties while preserving the nanosheet's overall stability. Additionally, doping in the central region allows for a more controlled and reproducible modification since edge sites often have different chemical reactivity and dangling bonds that could introduce variability. This strategy is consistent with computational studies that prioritize central or surface sites for transition metal doping to achieve stable configurations and desirable material properties.^35,40

Magnified representations of the BPH phosphorus atoms adjacent to each doped metal atom are depicted in Fig. S1, highlighting the lattice distortions and variations in P–P and P–metal bond lengths. Furthermore, the atomic partial charges and force field parameters associated with these metal atoms and the L-lysine functional group are summarized in Tables S1 and S2.

The electron configurations of nickel (Ni), copper (Cu), and zinc (Zn) are important factors in their ability to form coordination bonds with quercetin's hydroxyl (OH) and carbonyl (C=O) groups. Ni and Cu both have unpaired d electrons, which theoretically enhance their capacity to form strong coordination bonds with these functional groups on quercetin. Zinc, with a different electron configuration, typically forms weaker interactions. This is why we expect Ni and Cu to have stronger binding affinities, potentially leading to a more stable complex with quercetin, while Zn may show weaker binding, affecting the overall drug–carrier interaction. These interactions are based on well-established coordination chemistry principles, which are consistent with the general understanding of transition metal behavior.

In order to study the extent of drug adsorption, the nanosheets were placed in the center of a simulation box with dimensions of 11 × 11 × 11 nm³, and four drug molecules were placed at a distance of 2 nm from it. In addition, 5 L-lysine chains were covalently linked with doped nanosheets.^41–43 For this purpose, three other boxes were designed and the drug molecules were placed at a certain distance (approximately 2 nm) from the substrate. More details of the studied systems are given in Table S3. Fig. 1 shows the structure of the used carriers and the quercetin molecule.

Fig. 1 Chemical structure of (a) quercetin molecule, (b) L-lysine functional group, (c) pristine BPH nanosheet, and (d) Ni-BPH/L-lysine nanosheet. (The brown phosphorus atoms and the blue atoms in the center of the nanosheet are doped nickel atoms.)

Simulation details

This study's molecular topologies and simulation parameters were derived from the CHARMM36 force field⁴⁴ and relevant previous literature.^45,46 Each simulation box was filled with TIP3P model water molecules⁴⁷ and neutralized using sodium (Na⁺) and chloride (Cl⁻) ions at a concentration of 0.15 M. The energy of the systems was minimized using the steepest descent algorithm.

Subsequently, temperature and pressure were equilibrated using the Nosé–Hoover⁴⁸ and Parrinello–Rahman⁴⁹ algorithms in NVT and NPT ensembles at 310 K and 1 bar, respectively. Broadband electrostatic interactions were calculated using the Particle Mesh Ewald (PME)⁵⁰ method, while nonbonded interactions were managed with a cutoff distance of 1.4 nm. All systems were simulated for 105 ns with a time step of 1.5 fs, utilizing the GROMACS 2022.4 software package.⁵¹ The RMSD profiles (Fig. S2) demonstrate that all studied phosphorene nanosheet systems rapidly achieved equilibrium within the initial 15 ns and maintained stable conformations throughout the 105 ns molecular dynamics simulations, as indicated by steady RMSD plateaus in both panels. This behavior verifies the structural stability of pristine and L-lysine-functionalized, metal-doped phosphorene nanosheets upon quercetin adsorption. The absence of significant post-equilibration fluctuations confirms that each system remains dynamically stable over long timescales and is suitable for subsequent interaction analyses.

Uncertainty calculation

To estimate the uncertainties in the Coulomb and van der Waals (L-J) energies, we used the standard deviation of the energy values observed over the course of the 105 ns simulation period. Since molecular dynamics simulations are time-dependent, energy values fluctuate due to the changing configuration of the system. The uncertainty (σ_E) was computed by calculating the standard deviation of the energy values over time, as follows:^52–54where: std(E) is the standard deviation of the energy values, and N is the total number of data points (simulation time divided by the time step).

Discussion and results

Nanoparticles (NPs) are crucial in drug delivery because of their ability to transport substantial drug doses, which can be affected by various factors. Several factors significantly affect the drug-carrying capacity of nanoparticles (NPs) in drug delivery applications. These factors include the surface area and morphology of the nanosheet, as larger surface areas and more porous structures provide greater adsorption sites for the drug; surface functionalization, where the introduction of functional groups like hydroxyl or amine groups enhances the drug–nanosheet interaction through electrostatic and hydrogen bonding; metal doping, where metals like Ni, Cu, and Zn modify the surface charge and coordination sites, influencing drug binding; solvent conditions, such as pH and ionic strength, which affect the solubility and interaction of the drug with the nanosheet; and the drug's size and molecular structure, where smaller molecules or those with functional groups compatible with the nanosheet surface tend to bind more effectively. Understanding these factors is crucial for optimizing NP design and increasing therapeutic efficacy.^55–58 Hashemzadeh et al.⁵⁹ reported that two-dimensional BPHs could be suitable carriers for the delivery of anticancer drugs and their diffusion into membrane cells. In addition, the results of Tao et al.⁶⁰ showed that the functionalization of BPHs increases the drug loading capacity.

Phosphorene nanosheets have attracted attention as promising drug carriers due to their high surface area, biocompatibility, and ability to interact with various drug molecules. However, pristine phosphorene often exhibits limited interaction strength with drugs, which can reduce loading efficiency. Previous computational studies primarily focused on pristine phosphorene or simple doping with light elements like carbon and nitrogen, resulting in only modest improvements in drug adsorption energy and interfacial interactions.^61–63

Our study demonstrates that doping phosphorene with transition metals—such as nickel, copper, or zinc—significantly enhances the adsorption of quercetin on the nanosheet. Additionally, functionalizing these metal-doped nanosheets with the L-lysine group further improves drug adsorption by strengthening both electrostatic and van der Waals interactions. These combined modifications lead to increased drug loading capacity and more stable carrier–drug complexes, ultimately enhancing the suitability of phosphorene as an effective drug delivery platform.

This approach addresses the limitations of pristine phosphorene and aligns with recent advances in nanomaterial functionalization strategies aimed at improving drug delivery efficiency.^22,64 The study provides systematic, simulation-based evidence that guides the development of more efficient phosphorene-based nanocarriers for therapeutic molecules. Therefore, in this study, three doping agents and one amino acid were used to investigate the effect of these factors on the drug adsorption process compared to the pristine nanosheet. Initial and final images of the investigated systems are shown in Fig. 2 and S3. It can be seen that the presence of the L-lysine functional group caused better adsorption of the drug molecule quercetin on the nanosheet.

Fig. 2 Initial and final snapshots in the studied systems from top and side views, (a) Ni atom-doped phosphorene nanosheet, and (b) Ni atom-doped phosphorene nanosheet functionalized with L-lysine group.

In additional our simulations revealed the formation of coordination bonds between the Ni-doped nanosheet and quercetin's hydroxyl (–OH) and carbonyl (C=O) groups. These bonds fluctuated over the course of the 105 ns simulation, with coordination bonds being observed in approximately 3–4 frames. This dynamic yet stable interaction suggests that these coordination bonds are significant, though they appear intermittently in different simulation frames. The formation of these bonds was confirmed through the analysis of the simulation trajectory, which is shown in Fig. S3.

The average interaction energy values over the entire 105 ns time interval in different systems are presented in Table 1 and Fig. S4. These energy terms represent the nonbonded interactions between the metal-doped, L-lysine functionalized phosphorene nanosheet (host) and the adsorbed quercetin molecule (guest).

Table 1 The electrostatic (Coul.) and van der Waals (L-J) energy values in the studied systems (kJ mol⁻¹). A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups. (The uncertainty values for the studied systems are presented together with their corresponding energy values.)

System	Coul.	L-J	System	Coul.	L-J
A	−0.023 (0.044)	−203.49 (4.1)	A	−0.023 (0.044)	−203.49 (4.1)
B	−0.049 (0.015)	−226.95 (5.9)	F-B	−0.328 (0.096)	−257.72 (5.1)
C	−0.013 (0.014)	−224.08 (5.7)	F-C	−0.292 (0.12)	−255.44 (10)
D	−0.009 (0.019)	−213.10 (6.5)	F-D	−0.276 (0.14)	−221.65 (9)

---

It is observed that the presence of drug molecules near the substrate and the formation of interactions with the surface are mainly due to the formation of π–π stacks between the nanosheet surface and the drug molecule rings, as a result of which the contribution of van der Waals interactions has increased. A close examination of Table 1 shows that the energy values for functionalized systems have increased compared to the doped systems. This is due to the presence of the L-lysine functional group on the surface, the formation of electrostatic interactions between the two components, and the facilitation of van der Waals interactions. This fact also indicates the greater interaction of the drug molecule for adsorption on the functionalized nanosheet, which has led to the formation of a stable complex. The term “facilitation of van der Waals interactions” refers to the enhancement of these noncovalent forces resulting from surface functionalization with L-lysine. The additional functional groups provide more sites for close contact and favorable geometric alignment with the drug molecule, increasing the contact area and improving the spatial and electronic complementarity between the interacting surfaces. These conditions strengthen the dispersion (induced dipole–induced dipole) interactions between the nanosheet and quercetin, leading to a more stable and energetically favorable complex. Therefore, functionalization not only strengthens the electrostatic component but also promotes more effective van der Waals binding, resulting in greater adsorption stability of the drug molecule on the functionalized nanosheet compared to the doped systems.

Our simulations revealed that the quercetin–Ni complex formed more frequently and exhibited stronger interactions compared to the quercetin–Cu and quercetin–Zn complexes. Specifically, the number of interactions between quercetin and the Ni-doped nanosheet was higher, and the adsorption energy for the quercetin–Ni system was more favorable, indicating stronger binding. In contrast, while π–π interactions were observed between quercetin and both the Cu and Zn nanosheets, the formation of stable complexes in these systems was less frequent, and the interaction energies were lower, reflecting weaker binding compared to the quercetin–Ni system.

As shown in Fig. S4, the van der Waals (vdW) interaction energies between quercetin and phosphorene nanosheets doped with nickel, copper, and zinc were tracked throughout the simulation. All systems reached equilibrium within 20 ns and maintained stable, strong vdW interactions with energies ranging from −200 to −260 kJ mol⁻¹. Following functionalization of the doped nanosheets with L-lysine, these vdW energies decreased further (became more negative), reflecting an enhancement in the attractive interactions between quercetin and the functionalized surfaces. This increase in vdW contribution demonstrates the significant role of L-lysine functionalization in strengthening the adsorption and stability of quercetin on the nanosheets, confirming the effectiveness of functionalization in improving noncovalent binding.

The consistently small electrostatic interaction energies observed for all doped and functionalized phosphorene systems can be attributed to the predominantly neutral or weakly polar nature of the host and guest molecules in these simulations. Phosphorene and L-lysine, under the simulated conditions, possess limited net charge and do not feature strongly ionic sites, leading to minimal coulombic interactions with the quercetin adsorbate. The partial atomic charges assigned by the force field to both the nanosheet and drug molecule further contribute to attenuated electrostatic effects, compared to vdW or coordination-driven interactions. Consequently, the physical meaning of these small energies is that host–guest stabilization predominantly arises through dispersion forces and geometrical coordination, rather than through strong electrostatic attraction. This aligns with trends commonly observed in simulations of organic molecules on 2D nanomaterial surfaces,^63,65 where van der Waals interactions dominate binding unless highly charged species are introduced.

On the other hand, it is observed that the adsorption process in both pristine and functionalized complexes follows the following trend: Ni > Cu > Zn.

Nickel doping on phosphorene creates local coordination sites that enhance interactions with the functional groups of quercetin. The nickel atom, thanks to its electronic structure and spatial characteristics, serves as an effective center that promotes stronger electrostatic and van der Waals interactions with the hydroxyl (–OH) and carbonyl (C=O) groups of quercetin. Molecular dynamics simulations reveal these heightened interactions through changes in adsorption geometry and increased binding energy, resulting in a more stable Ni-BPH–quercetin complex compared to systems doped with copper or zinc. This increased stability is primarily due to the geometric arrangement and distribution of partial charges. The coordination effects induced by nickel are observed in both biological and materials modeling contexts, confirming the superior formation and stability of the complex.^66–68

Most drug molecules are adsorbed on the surface through their aromatic rings, and the presence of functional groups also causes hydrogen bonding between the surface and the drug molecule.

In this system, H-bonds were identified using a donor–H⋯acceptor distance ≤ 3.5 Å and donor–H–acceptor angle ≥ 150°. In our systems, the principal donors are quercetin phenolic O–H groups and lysine –NH₃⁺/–NH– groups (when present), with water also acting as a donor/acceptor in bridging events. The main acceptors are quercetin carbonyl oxygen (C=O) and phenolic oxygen (when deprotonated or engaged as acceptor), surface phosphorus atoms of BPH (lone-pair bearing P acting as weak acceptors), and water oxygen. Metal centers (Ni, Cu, Zn) were not counted as H-bond acceptors; their interactions with quercetin are treated separately as coordination.

Fig. 3 shows the number of hydrogen bonds between the drug molecules and the lysine functionalized substrate. It is noted that the hydrogen bonds formed correspond to the electrostatic energy values.

Fig. 3 Number of hydrogen bond diagram for the studied systems A pristine phosphorene, F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups.

To assess the probability of drug molecules being present around the substrate, we conducted an analysis using radial distribution function (RDF) calculations. Fig. 4 presents the RDF diagrams for the systems under investigation.

Fig. 4 RDF diagram for the studied systems A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups.

In systems doped with various atoms (Fig. 4, panel (a)), the probability of drug molecules being present around the substrate is higher compared to pure phosphorene. This finding indicates that doping the system with different atoms can enhance its adsorption capacity. While there are minor variations in the trends regarding drug presence around the three substrates doped with different atoms, it is noteworthy that the phosphorene nanosheet doped with nickel exhibits the highest probability of drug molecule presence. This finding highlights the potential of selective doping to optimize drug adsorption in phosphorene-based materials.

Fig. 4, panel (b), presents the radial distribution function for the functionalized systems. This diagram reveals a significant increase in the probability of drug molecules being present around the substrate compared to both pure phosphorene and phosphorene doped with various atoms. The major peaks in these plots occur at a distance of approximately 0.6 nm, indicating π–π interactions between the drug molecule rings and the nanosheet surface. Additionally, an examination of Fig. S5 confirms that the drug molecules interact with the nanosheet surface through π–π interactions and are oriented parallel to the surface.

During the simulation, drug molecules come into contact with the nanosheet surface over time. The frequency of these contacts is calculated using eqn (2).⁶⁹

Here, N_BPH and N_QES represent the total number of BPH and QES atoms, respectively. The variable r_j denotes the distance from the j-th atom of QES to the i-th atom of BPH. The δ function is utilized to count all drug atoms that are within 0.6 nm of BPH atoms.

Fig. 5 illustrates the number of contacts between the drug molecule and the nanosheet. This figure shows that, for all the systems studied, the drug molecule does not make contact with the surface at the beginning of the process. In the doped systems with Ni and Zn (Fig. 5a), contact between the drug molecule and the surface begins at approximately 7 ns. In the pristine phosphorene system and the copper atom-doped nanosheet, a peak in contacts is observed between 10 and 15 ns. This peak is attributed to the orientation of the drug molecules on the surface as they seek the optimal configuration. Eventually, in all systems, the number of collisions reaches a maximum value and then fluctuates around a constant level (2100).

Fig. 5 Diagram of the number of contacts for the studied systems A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups.

Fig. 5b illustrates the number of contacts between the drug molecule and the functionalized systems. The data indicate that in systems functionalized with L-lysine groups, the number of contacts increased compared to the pristine nanosheet. In systems F-B and F-D, the number of contacts peaks at approximately 7 ns and then fluctuates around this value for the remainder of the simulation. In contrast, system F-C shows a decrease in contacts between 30 and 60 ns, which is attributed to the repeated reorientation of the drug molecule on the surface (see Fig. S6). After this period, the number of contacts increases as the drug molecule adopts a more favorable position. These observations align with the findings from the energy analysis and the results of the radial distribution function (RDF).

Additionally, the number of interactions between the doped metal atoms on the nanosheet and the drug molecules was calculated, as shown in Fig. S7. Generally, the unfunctionalized doped systems exhibit greater fluctuations in contact number, reflecting intermittent periods of strong interaction between the drug and the doped nanosheets. After functionalization, peak contact values tend to decrease, and the interaction profiles become more stable and sustained compared to the unfunctionalized systems. Functionalization with an amino acid typically reduces the frequency of contact peaks, suggesting that surface modification can modulate and stabilize the nanosheet–drug interactions, possibly through steric hindrance or altered surface energy. These interaction profiles offer valuable insights into how nanosheet surface chemistry influences drug dynamics and can guide the design of drug delivery platforms with controllable interaction characteristics.

The distance between the nanosheet and the drug molecule as a function of simulation time is illustrated in Fig. 6. It can be observed that for systems doped with different atoms (Fig. 6a), the distance between the centers of mass of the two components decreases from approximately 1.2 nm to about 0.35 nm. Fig. 6b shows that with the functionalization of the nanosheet surface, the interaction distance between the substrate and the drug is reduced to 0.2 nm. This finding confirms that the presence of functional groups improves the adsorption process. On the other hand, this is consistent with the results obtained from the energy and number of interatomic contacts.

Fig. 6 Diagram of the COM distance between quercetin molecule and nanosheet for the studied systems A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups.

The results presented in Table 2 show the D_i values, to estimate the mobility of molecules in the systems under study. System D exhibits the highest diffusion coefficient in both its doped and functionalized states, indicating a weaker interaction between the drug molecule and the substrate. Conversely, the lower D_i values in the functionalized systems suggest a stronger interaction between the drug and the functionalized nanosheet. Additionally, the diffusion coefficient results are supported by the slope of the Mean Squared Displacement (MSD) diagram shown in Fig. 7, which indicates that system B has the highest diffusion rate around the substrates.

Table 2 Diffusion coefficient values for the studied systems A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups

System	Di (×10^−5 cm^2 s^−1)	System	Di (×10^−5 cm^2 s^−1)
A	0.0264 (±0.0005)	A	0.0264 (±0.0005)
B	0.0079 (±0.0003)	F-B	0.0011 (±0.0007)
C	0.0138 (±0.0016)	F-C	0.0022 (±0.0003)
D	0.0386 (±0.0008)	F-D	0.0063 (±0.0008)

---

Fig. 7 MSD diagram for the studied systems A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups.

For the interpretation of Fig. 7, we select the final 25 ns because this segment shows quasi-steady, near-linear MSD growth suitable for extracting D_i earlier times are dominated by approach, re-orientation, and adsorption transients that are not representative of steady surface diffusion. Fig. 7 reports the mean-squared displacement, MSD, for the final 25 ns of the 105 ns trajectories (80–105 ns), a time window in which MSD grows approximately linearly and thus permits estimation of diffusion coefficients via the Einstein relation.

One of the key parameters for predicting how drug molecules interact with nanosheets, as well as for understanding molecular interactions and behaviors in biological systems, is the solvent-accessible surface area (SASA). Introduced by Lee and Richards in 1971,⁷⁰ SASA quantifies the surface area of a molecule that is accessible to solvent molecules. This measurement effectively indicates the extent to which a molecule can interact with its environment, providing valuable insights into its potential biological activity and interactions.

The solvent-accessible surface area (SASA) for the drug molecule is illustrated in Fig. 8. In doped systems, the SASA remains relatively consistent across all configurations, fluctuating around 19.7 nm² (see Fig. 8a). This consistency suggests that the adsorption of drug molecules in these systems is nearly uniform.

Fig. 8 Hydrophobic SASA for the studied systems A pristine phosphorene, B, C, and D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms, and F-B, F-C, and F-D, respectively, phosphorene nanosheets doped with Ni, Cu, and Zn atoms functionalized with L-lysine groups.

In contrast, the functionalized systems demonstrate a notable reduction in SASA, particularly in the F-B system (Fig. 8b). This reduction indicates that drug molecules in this system are well-oriented on the substrate surface, enhancing their interaction with it and consequently decreasing their interaction with water. A comparison of SASA values between functionalized and doped systems highlights the significant impact of the L-lysine functional group on the drug adsorption process. The average SASA values over the 0–105 ns simulation period reveal minor variations among the systems, consistent with the nearly overlapping curves observed in Fig. 8, and the corresponding data are summarized in Table S4. This finding underscores the importance of functionalization in optimizing drug–nanosheet interactions.

Metadynamics

Parrinello⁷¹ and Laio⁷² utilized metadynamic simulation techniques to explore the free energy landscape of molecular systems. This approach is specifically designed to address the limitations of traditional molecular dynamics simulations, particularly in the context of rare events and high-energy barriers.

Metadynamics relies on a predefined set of collective variables (CVs) that characterize the system's configuration along a reaction coordinate (CV (collective variable): a predefined set of parameters or coordinates used in metadynamics simulations to describe the system's configuration along a reaction coordinate. In this case, the CV represents the distance between the centers of mass of the quercetin molecule and the BPH nanosheet surface, which is monitored throughout the simulation to evaluate the adsorption process and the free energy profile). This method, implemented through the PLUMED plugin,⁷³ operates based on the distances between the centers of mass of two components. Hence, a thorough metadynamic simulation for the pristine BPH system and Ni-BPH/L-lysine system was conducted throughout 105 ns.

The estimate of the free energy level calculated using the sum_hills program is depicted in Fig. 9. A detailed analysis shows that as the drug molecule transitions from the aqueous phase to the adsorbent, the free energy level decreases in both systems.

Fig. 9 The free energy profile for the adsorption of (a) quercetin on the pristine phosphorene nanosheet and (b) quercetin on the Ni atom-doped phosphorene nanosheet functionalized with L-lysine group as a function of the center of mass.

In the pristine phosphorene system, the drug molecule overcomes two energy barriers of 5 kJ mol⁻¹ and 13 kJ mol⁻¹ at distances of approximately 3.5 nm and 2.5 nm, respectively. It then reaches its lowest free energy value of −186 kJ mol⁻¹ at a distance of about 1 nm, after overcoming an additional energy barrier of 10 kJ mol⁻¹ (as shown in Fig. 9a). In contrast, in the Ni-BPH/L-lysine system, the drug molecule initially surmounts an energy barrier of 28 kJ mol⁻¹ at a distance of approximately 6 nm. This increased barrier may be attributed to the functionalization of the phosphorene substrate surface. Furthermore, as the drug molecule approaches a distance of around 3.4 to 2.5 nm, it reaches its most negative free energy value after crossing a significant energy barrier of approximately 48 kJ mol⁻¹. The global minimum free energy value for this system is observed to be −248 kJ mol⁻¹ (Fig. 9b).

Conclusion

This study investigates the behavior of the drug quercetin in its adsorption onto phosphorene nanosheets. The results of molecular dynamics simulation calculations indicate that doping phosphorene nanosheets can play a significant role in the delivery of anticancer drugs. Moreover, the combined effects of doping and functionalization of the nanosheets lead to promising outcomes for drug adsorption. Functionalization improves drug adsorption by enhancing electrostatic and van der Waals interactions, which increases drug retention and reduces its diffusion into the surrounding environment. The presence of doped atoms, around which drug molecules adsorb, along with the formation of π–π interactions, significantly influences the overall drug adsorption process. The correlation between metadynamic calculations, distance, solvent-accessible surface area (SASA), and diffusion coefficients supports our findings. The insights gained from this study highlight the potential of doped and functionalized phosphorene nanosheets as effective carriers for anticancer drugs, paving the way for advancements in targeted drug delivery systems.

Author contributions

Fatemeh Talaei: devised the computational protocol and prepared the model systems, performed all calculations, analyzed the data, and wrote and edited the original and the revised manuscript. Farzaneh Farzad: supervision. Reviewing – editing, edited the original and the revised version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The codes that support the findings of this study are available from the corresponding authors upon reasonable request. The authors confirm that the simulation data are presented in the form of graphs and tables and are analyzed within the text of the article. If you require the original data in specialized files for molecular dynamics simulation due to their large volume, you can send an email to E-mail: fatemeh_talaei96@birjand.ac.ir, and the data will be made available to you.

Supplementary information is available. See DOI: https://doi.org/10.1039/d5ra06137g.

References

1. T. Sahu, Y. K. Ratre, S. Chauhan, L. Bhaskar, M. P. Nair and H. K. Verma, Nanotechnology based drug delivery system: Current strategies and emerging therapeutic potential for medical science, J. Drug Delivery Sci. Technol., 2021, 63, 102487.
2. M. I. Khan, et al., Recent progress in nanostructured smart drug delivery systems for cancer therapy: a review, ACS Appl. Bio Mater., 2022, 5(3), 971–1012.
3. G.-A. Croitoru, et al., Nanostructured drug delivery systems in immunotherapy: an updated overview of nanotechnology-based therapeutic innovations, Appl. Sci., 2024, 14(19), 8948.
4. S. Riaz, K. Jaffar, M. Perveen, A. Riaz, S. Nazir and J. Iqbal, Computational study of therapeutic potential of phosphorene as a nano-carrier for drug delivery of nebivolol for the prohibition of cardiovascular diseases: a DFT study, J. Mol. Model., 2021, 27, 1–15.
5. L. Razavi, H. Raissi and F. Farzad, Graphene oxide and silicene as 2D platforms for complexation and intracellular delivery of siRNA, J. Drug Delivery Sci. Technol., 2024, 95, 105514.
6. F. Talaei, F. Farzad and A. Yaghobi, Molecular insights into functionalized carbon nanotubes for the adsorption of therapeutic peptides, Results Mater., 2025, 26, 100704.
7. B. Chen, Y. Dai, S. Yang, C. Chen and L. Wang, Recent progress on transition metal dichalcogenide-based composites for cancer therapy, Nanoscale, 2025, 17, 7552–7573.
8. A. Haqyar, H. Raissi, F. Farzad and H. Hashemzadeh, A strategy toward therapeutic improvement of electric field-sensitive gemcitabine prodrugs in 2D metal–organic frameworks in view of their structure and interactions, Inorg. Chem. Commun., 2022, 138, 109281.
9. T. Hu, X. Mei, Y. Wang, X. Weng, R. Liang and M. Wei, Two-dimensional nanomaterials: fascinating materials in biomedical field, Sci. Bull., 2019, 64(22), 1707–1727.
10. J. Santos, M. Moschetta, J. Rodrigues, P. Alpuim and A. Capasso, Interactions between 2D materials and living matter: a review on graphene and hexagonal boron nitride coatings, Front. Bioeng. Biotechnol., 2021, 9, 612669.
11. G. Peng and B. Fadeel, Understanding the bidirectional interactions between two-dimensional materials, microorganisms, and the immune system, Adv. Drug Delivery Rev., 2022, 188, 114422.
12. F. Najafi, F. Farzad and S. Pasban, Interactions of boron nitride nanosheet with amino acids of differential polarity, Sci. Rep., 2022, 12(1), 11156.
13. C. M. Ledur, I. Zanella and S. B. Fagan, Ab initio simulations of black and blue phosphorene functionalised with chemical groups for biomolecule anchoring, J. Mol. Model., 2021, 27, 1–11.
14. T. Mahboobi, M. R. Zardoost and M. R. Toosi, A comparative DFT study of the effect of doping atoms of groups III, IV, and V on the electronic properties of phosphorene, Struct. Chem., 2022, 33(1), 131–145.
15. T. Mahboobi, M. R. Zardoost and M. R. Toosi, Phosphorene and Na-, Ca-, and Fe-doped phosphorene as candidates for delivery of mercaptopurine and fluorouracil anticancer drugs, J. Mol. Model., 2020, 26(10), 269.
16. M. Fu, et al., Studying adsorption and cellular toxicity of boron nitride nanostructure versus melphalan anti-ovarian cancer drug, Curr. Mol. Med., 2021, 21(8), 698–705.
17. H. Karim, et al., A DFT investigation on theranostic potential of alkaline earth metal doped phosphorenes for ifosfamide anti-cancer drug, Appl. Surf. Sci., 2022, 596, 153618.
18. M. Qiu, et al., Novel concept of the smart NIR-light–controlled drug release of black phosphorus nanostructure for cancer therapy, Proc. Natl. Acad. Sci. U. S. A., 2018, 115(3), 501–506.
19. A. Tariq, S. Nazir, A. W. Arshad, F. Nawaz, K. Ayub and J. Iqbal, DFT study of the therapeutic potential of phosphorene as a new drug-delivery system to treat cancer, RSC Adv., 2019, 9(42), 24325–24332.
20. A. Thakur, A. Kumar, S. Kaya, R. Marzouki, F. Zhang and L. Guo, Recent advancements in surface modification, characterization and functionalization for enhancing the biocompatibility and corrosion resistance of biomedical implants, Coatings, 2022, 12(10), 1459.
21. F. Ahmad, et al., Unique properties of surface-functionalized nanoparticles for bio-application: functionalization mechanisms and importance in application, Nanomaterials, 2022, 12(8), 1333.
22. A. Bigham, et al., Black phosphorus-based nanoplatforms for cancer therapy: chemistry, design, biological and therapeutic behaviors, Chem. Soc. Rev., 2025, 54(2), 827–897.
23. X. Xuan, Z. Zhang and W. Guo, Doping-stabilized two-dimensional black phosphorus, Nanoscale, 2018, 10(17), 7898–7904.
24. X. Yang, et al., Functional black phosphorus nanosheets for mitochondria-targeting photothermal/photodynamic synergistic cancer therapy, Chem. Sci., 2019, 10(13), 3779–3785.
25. Z. Zhang, et al., Construction and evaluation of a folate-mediated adriamycin phosphorene thermosensitive hydrogel delivery system against osteosarcoma, Appl. Mater. Today, 2024, 40, 102386.
26. S. Yin, et al., Polydopamine-modified black phosphorus nanosheet drug delivery system for the treatment of ischemic stroke, Regener. Biomater., 2024, 11, rbae046.
27. R. Dennington, T. Keith and J. Millam, et al., GaussView, Version 5, 2009.
28. V. Zoete, M. A. Cuendet, A. Grosdidier and O. Michielin, SwissParam: a fast force field generation tool for small organic molecules, J. Comput. Chem., 2011, 32(11), 2359–2368.
29. M. Ghashghaee and M. Ghambarian, Defect engineering and zinc oxide doping of black phosphorene for nitrogen dioxide capture and detection: quantum-chemical calculations, Appl. Surf. Sci., 2020, 523, 146527.
30. L. M. Moreau, D.-H. Ha, C. R. Bealing, H. Zhang, R. G. Hennig and R. D. Robinson, Unintended phosphorus doping of nickel nanoparticles during synthesis with TOP: a discovery through structural analysis, Nano Lett., 2012, 12(9), 4530–4539.
31. A. Kumar, et al., Black Phosphorus n-type doping by Cu: a microscopic surface investigation, J. Phys. Chem. C, 2021, 125(24), 13477–13484.
32. M. Mishra, S. Sushama, S. K. Pandey and S. Chakrabarti, Phosphorus doping of ZnO using spin-on dopant process: A better choice than costly and destructive ion-implantation technique, J. Lumin., 2021, 233, 117921.
33. J. Calvario-Coyotl, S. J. Gutierrez-Ojeda, P. Rubio-Pereda, S. Esteban-Gómez, L. M. de la Garza and G. H. Cocoletzi, Doping effects in the alanine-black phosphorene interactions: first principles studies, Chem. Phys., 2024, 579, 112214.
34. H. Sun, Y. Shang, Y. Yang and M. Guo, Realization of N-Type Semiconducting of Phosphorene through Surface Metal Doping and Work Function Study, J. Nanomater., 2018, 2018(1), 6863890.
35. S. Gazzari, K. Wrighton-Araneda and D. Cortes-Arriagada, A first-principles description of the stability of transition-metal doped phosphorene nanosheets, Surf. Interfaces, 2020, 21, 100786.
36. P. Srivastava, K. Hembram, H. Mizuseki, K.-R. Lee, S. S. Han and S. Kim, Tuning the electronic and magnetic properties of phosphorene by vacancies and adatoms, J. Phys. Chem. C, 2015, 119(12), 6530–6538.
37. S. Gazzari, K. Wrighton-Araneda and D. Cortés-Arriagada, Theoretical description of the physical/chemical contributions determining the stability of transition-metal doped phosphorene nanosheets, ChemRxiv, 2020, preprint,  DOI:10.26434/chemrxiv.12686012.v1.
38. A. Marjani, M. Ghambarian and M. Ghashghaee, Alkali metal doping of black phosphorus monolayer for ultrasensitive capture and detection of nitrogen dioxide, Sci. Rep., 2021, 11(1), 842.
39. A. S. Nair, R. Ahuja and B. Pathak, Unraveling the single-atom electrocatalytic activity of transition metal-doped phosphorene, Nanoscale Adv., 2020, 2(6), 2410–2421.
40. R. Babar and M. Kabir, Transition metal and vacancy defect complexes in phosphorene: a spintronic perspective, J. Phys. Chem. C, 2016, 120(27), 14991–15000.
41. M. Zheng, et al., Poly(α-l-lysine)-based nanomaterials for versatile biomedical applications: Current advances and perspectives, Bioact. Mater., 2021, 6(7), 1878–1909.
42. A. Li, S. Wang, Z. Zhang, N. Xu, G. Ling and P. Zhang, Poly-l-lysine derivative-coated black phosphorus as a nanoplatform for photothermal chemotherapy to enhance anti-tumor efficiency, J. Mater. Chem. B, 2022, 10(27), 5191–5202.
43. Y. Zhao, Y.-H. Zhang, Z. Zhuge, Y.-H. Tang, J.-W. Tao and Y. Chen, Synthesis of a poly-L-lysine/black phosphorus hybrid for biosensors, Anal. Chem., 2018, 90(5), 3149–3155.
44. J. Huang and A. D. MacKerell Jr, CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data, J. Comput. Chem., 2013, 34(25), 2135–2145.
45. P. P. Putra, N. Fitria, A. Asnawi and A. Djamaan, Metalloprotein parameters in molecular dynamics simulation for AMBER, CHARMM, GROMACS, and NAMD: A systematic review, Int. J. Appl. Pharm., 2024, 16(5), 90–98.
46. A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III and W. M. Skiff, UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations, J. Am. Chem. Soc., 1992, 114(25), 10024–10035.
47. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79(2), 926–935.
48. D. J. Evans and B. L. Holian, The Nose–Hoover thermostat, J. Chem. Phys., 1985, 83(8), 4069–4074.
49. H. Saito, H. Nagao, K. Nishikawa and K. Kinugawa, Molecular collective dynamics in solid para-hydrogen and ortho-deuterium: The Parrinello–Rahman-type path integral centroid molecular dynamics approach, J. Chem. Phys., 2003, 119(2), 953–963.
50. H. G. Petersen, Accuracy and efficiency of the particle mesh Ewald method, J. Chem. Phys., 1995, 103(9), 3668–3679.
51. M. J. Abraham, et al., GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1, 19–25.
52. S. Wan, R. C. Sinclair and P. V Coveney, Uncertainty quantification in classical molecular dynamics, Philos. Trans. R. Soc., A, 2021, 379(2197), 20200082.
53. V. Zaverkin, et al., Uncertainty-biased molecular dynamics for learning uniformly accurate interatomic potentials, npj Comput. Mater., 2024, 10(1), 83.
54. A. P. Bhati, S. Wan, Y. Hu, B. Sherborne and P. V Coveney, Uncertainty quantification in alchemical free energy methods, J. Chem. Theory Comput., 2018, 14(6), 2867–2880.
55. Z. Chenxi, A. Hemmat, N. H. Thi and M. Afrand, Nanoparticle-enhanced drug delivery systems: An up-to-date review, J. Mol. Liq., 2025, 126999.
56. M. J. Mitchell, M. M. Billingsley, R. M. Haley, M. E. Wechsler, N. A. Peppas and R. Langer, Engineering precision nanoparticles for drug delivery, Nat. Rev. Drug Discovery, 2021, 20(2), 101–124.
57. K. Elumalai, S. Srinivasan and A. Shanmugam, Review of the efficacy of nanoparticle-based drug delivery systems for cancer treatment, Biomed. Technol., 2024, 5, 109–122.
58. D. M. Cahyani, et al., Nanoparticle tools for maximizing oral drug delivery, Braz. J. Med. Biol. Res., 2025, 58, e14459.
59. H. Hashemzadeh and H. Raissi, Design of new drug delivery platform based on surface functionalization of black phosphorus nanosheet with a smart polymer for enhancing the efficiency of doxorubicin in the treatment of cancer, J. Biomed. Mater. Res., Part A, 2021, 109(10), 1912–1921.
60. W. Tao, et al., Black phosphorus nanosheets as a robust delivery platform for cancer theranostics, Adv. Mater., 2016, 29(1), 1603276.
61. W. Chen, et al., Black phosphorus nanosheet-based drug delivery system for synergistic photodynamic/photothermal/chemotherapy of cancer, Adv. Mater., 2016, 29(5), 10.1002/adma.201603864.
62. W. Liang and X. Luo, Theoretical studies of MoS₂ and phosphorene drug delivery for antituberculosis drugs, J. Phys. Chem. C, 2020, 124(15), 8279–8287.
63. R. Esfandiarpour, F. Badalkhani-Khamseh and N. L. Hadipour, Theoretical studies of phosphorene as a drug delivery nanocarrier for fluorouracil, RSC Adv., 2023, 13(26), 18058–18069.
64. M. S. Islam, M. T. Islam and M. R. Hossain, Phosphorene: A novel nanomaterial revolutionizing biomedicine, JCIS Open, 2024, 16, 100124.
65. S. Gao, Y. Wang, J. Sun and Z. Zhang, Recent development in black phosphorus nanomaterials for anticancer applications, Front. Biomater. Sci., 2023, 2, 1172524.
66. A. van Teijlingen, H. W. A. Swanson, K. H. A. Lau and T. Tuttle, Constant pH coarse-grained molecular dynamics with stochastic charge neutralization, J. Phys. Chem. Lett., 2022, 13(18), 4046–4051.
67. A. P. Latham, E. S. Levy, B. D. Sellers and D. H. Leung, Utilizing molecular simulations to examine nanosuspension stability, Pharmaceutics, 2023, 16(1), 50.
68. G. A. Ross, A. S. Rustenburg, P. B. Grinaway, J. Fass and J. D. Chodera, Biomolecular simulations under realistic macroscopic salt conditions, J. Phys. Chem. B, 2018, 122(21), 5466–5486.
69. H. Hashemzadeh and H. Raissi, The functionalization of carbon nanotubes to enhance the efficacy of the anticancer drug paclitaxel: a molecular dynamics simulation study, J. Mol. Model., 2017, 23, 1–10.
70. B. Lee and F. M. Richards, The interpretation of protein structures: estimation of static accessibility, J. Mol. Biol., 1971, 55(3), 379–400.
71. E. Gandini, F. Dapiaggi, F. Oliva, S. Pieraccini and M. Sironi, Well-Tempered MetaDynamics based method to evaluate universal peptidomimetics, Chem. Phys. Lett., 2018, 706, 729–735.
72. G. Bussi and A. Laio, Using metadynamics to explore complex free-energy landscapes, Nat. Rev. Phys., 2020, 2(4), 200–212.
73. M. Bonomi, et al., PLUMED: A portable plugin for free-energy calculations with molecular dynamics, Comput. Phys. Commun., 2009, 180(10), 1961–1972.

---