DFT calculationsThe results from the DFT calculations including computation of electronic structure and binding energy are reported in this section. We also discuss the Molecular Electrostatic Potential (MEP), a widely-used concept in computational chemistry for depicting molecular polarization32. The MEP provides valuable insights into the distribution of electron density within a molecule, which is essential for understanding its interactions with other molecules and overall reactivity32. Greater polarization in a molecule refers to a higher separation between its positive and negative charges, which can significantly impact the molecular interactions and binding affinity between PEG and drug molecules. A more polarized molecule generally exhibits stronger dipole moments, leading to enhanced electrostatic interactions with other polarizable species. In the context of PEG-drug complexes, greater polarization in the drug molecule can result in a stronger attraction between the polar regions of the drug and the PEG chain, contributing to more stable complex formation.The MEP results for drug molecules, visualized in Fig. 1, indicate that HINA and HINB display greater polarization than MAYA and CELB, due to more functional groups and oxygen atoms, increasing molecular polarity. High polarization in HINA and HINB could affect interactions with PEG carriers and solvation free energy. The MEP figure for PEG.drug system (as shown in Fig. S1 of the SI) follows the same trend.Figure 1Visualization of MEP obtained from DFT calculations for drug molecules of HINA (a), HINB (b), MAYA (c), and CELB (d). The red color indicates regions of negative electrostatic potential, which are electron-rich areas. The blue color indicates regions of positive electrostatic potential, which are electron-deficient areas.As described earlier in the method sections, the optimized structures of the two lowest-energy conformers for each PEG.drug complex are presented in Fig. 2. We observed that the main interactions between PEG and drug molecules are the vdW interactions between heavy atoms. The heavy atoms primarily involved in these interactions include carbon and oxygen atoms from both the PEG and drug molecules. These interactions are particularly significant in this context because they govern the overall stability of PEG-drug complexes. These interactions contribute to the formation of stable complexes through a network of weak, non-covalent contacts that span across the PEG and drug molecules.Figure 2Optimized structures of unbound PEG.drug complexes obtained from DFT calculations. The PEG molecules are presented in green color for a better visualization.The DFT-computed binding energy and solvation free energy were summarized in Table 1. The values of binding energies (\(\Delta\)E\(_{bind}^{DFT}\)) reveal favorable drug-PEG interactions for all complexes as indicated by negative values between − 69 and − 171 kJ/mol. This suggests that the PEG.drug binding is energetically favorable under the gas-phase conditions.Table 1 also presents the solvation free energies (\(\Delta\)G\(_{solv}^{DFT}\)) of the drug molecules follow the trend: HINB > MAYA > HINA > CELB. This aligns with their MEP results, where HINB shows the highest polarity and CELB displays the lowest polarity (Fig. 1). It is not surprising that PEG exhibits more favourable solvation free energies than the drug molecules due to its hydrophilic behavior1 , which promote stronger interactions with the water solvent.Table 1 Calculated bind energy (kJ/mol) and solvation free energy (kJ/mol) obtained from DFT calculation for unbound PEG.drug complexes.Upon comparing solvation free energies of the drug molecules and their corresponding unbound complexes with PEG carriers, we observed a significant increase in solvation free energy upon complex formation (almost three times higher for the PEG.drug complexes compared to individual drug molecules). However, determining the thermodynamics of these interactions requires considering delta solvation free energy (\(\Delta \Delta\)G\(_{solv}^{DFT}\)), which ranged from 17 to 54 kJ/mol (see Table 1). The positive values of the delta solvation free energy indicate that the solvation process is not favorable for complex formation process.The reaction energies of the PEGylation process for the four bound PEG-drug systems, obtained through DFT calculations, are presented in Table 2. The results reveal that the formation of PEG-HINA exhibits the most favorable \(\Delta\)E\(_{reaction}^{DFT}\) (− 75 kJ/mol), while the formation of PEG-MAYA is unfavorable (80 kJ/mol). As demonstrated in Fig. 3, the MAYA molecule forms a specific bonding position with the PEG fragment, it generally tends to have fewer contacts with the PEG chain compared to the other drug molecules. This reduced interaction may contribute to the less favorable reaction energy observed for the formation of PEG-MAYA.The results indicate that the PEGylation reaction leading to the formation of PEG-drug is favorable for HINA, HINB, and CELB systems, while it is unfavorable for MAYA. The interaction between the drug molecule and the PEG fragment appears to be a crucial factor in determining this reaction.Table 2 Calculated reaction energy (kJ/mol) for the PEGylation process to form covalent bonded PEG-drug molecules.Figure 3Optimized structures of covalent bound PEG-drug molecules obtained from DFT calculations. The PEG fragments are presented in green color for a better visualization.xTB calculationsIn this section, the xTB method was utilized to compute solvation free energies and binding energies of PEG.drug complexes. The xTB results indicate that drug molecules display significantly higher solvation free energies upon forming complexes with PEG carriers, consistent with the DFT findings. The solvation free energy of the complexes is approximately five times greater than that of the individual drug molecule (see Table 3). This enhancement primarily comes from PEG’s favorable solvation. The difference in solvation free energy of the complexed system (\(\Delta \Delta\)G\(_{solv}^{xTB}\)) remains unfavorable (Table 4). The positive \(\Delta \Delta\)G\(_{solv}^{xTB}\) values, ranging from 6 to 34 kJ/mol, again imply an unfavorable solvation energy for the unbound PEG.drug systems.Table 3 Solvation free energy (kJ/mol) obtained from xTB calculations for the unbound PEG.drug systems.Table 4 presents a comprehensive analysis of the binding energies between PEG carriers and various drug molecules, as computed by the xTB method. The table lists the calculated gas-phase interaction energies, denoted as \(\Delta E_{bind}^{xTB}\), which range from − 113 to − 72 kJ/mol. These negative values indicate favorable interactions between the PEG carriers and the drugs in a gaseous state, suggesting that such complexes can form spontaneously under vacuum conditions. Upon including the solvent effects, the originally determined binding energies evolved into binding enthalpies (\(\Delta\)H\(_{bind}^{xTB}\)). The resulting values, as presented in Table 4, show that the binding enthalpy ranges from − 57 to − 79 kJ/mol. This negative values indicate that the interaction between PEG carriers and drug molecules is characterized by a net release of energy, which is typically associated with favorable or exothermic binding processes.However, upon closer examination of the entropy contributions, which were computed through harmonic frequency analysis, it becomes evident that these interactions are not entropically favorable. The entropic effects, quantified by the change in entropy (\(-T\Delta S\)), yield positive values for all complexes. This result indicates that the formation of PEG-drug complexes is associated with a decrease in disorder or entropy at the molecular level. The resulting values for the binding free energies (using Eq. 3), as shown in Table 4, are predominantly positive, with a range of − 8 to 18 kJ/mol. This implies that the favorable interactions observed in the gas phase are largely negated by the entropic penalty and solvation effects upon transitioning to an aqueous environment. Moreover, the positive binding free energies suggest that the formation of PEG-drug complexes is entropy-driven, with water molecules disrupting the ordering of the system upon complexation. This entropic penalty outweighs the enthalpic benefits of the molecular interaction, leading to unfavorable thermodynamics for the stability of the complex in aqueous solution. Consequently, these findings indicate that PEG-drug complexes may be inherently unstable when exposed to water, which has significant implications for their potential application and design in pharmaceutical formulations.Table 4 Calculated binding energy, solvation free energy, enthalpy, entropy and binding free energy obtained from xTB calculations for unbound PEG.drug complexes.The present results highlights the significance of accounting for both thermodynamic and entropic factors when engineering drug delivery systems based on compounds derived from Celastrus hindsii and PEG carriers. While the binding energy between PEG and the drug compound are favorable, the substantial unfavorable entropy change and delta solvation free energy contribute to the overall instability of the complex. The findings from static and implicit water simulations provide valuable insights; however, these approaches might be limited in their ability to accurately capture the intricate interactions between PEG and the drug carrier in a more realistic aqueous environment. To verify and further elucidate these results, a series of MD simulations were conducted. These simulations have incorporated explicit water molecules, allowing for a more detailed investigation into the dynamic behavior and interaction between PEG and the drug carrier in solutions.MD simulationsIn the molecular dynamics study, we first analyzed the separation of PEG carrier and bioactive drug molecules during the simulation process. As depicted in Fig. 4, which shows a typical snapshot of the last MD frame for the unbound system after 100 ns of simulations, it is evident that the PEG carrier and the drug molecule have indeed separated from each other. This observation aligns with findings reported in previous literature on the interaction between PEG and small molecules in the unbound state23. The MD simulations, in conjunction with previous DFT and xTB simulations, consistently demonstrate the separation of the PEG carrier and bioactive drug molecules. This observation shows the reliability and consistency of the computational methods employed in studying the behavior of these systems.Figure 4Representative snapshots of the last MD frame for the unbound systems after 100 ns of simulations, demonstrating the separation of the PEG carrier and the bioactive drug molecules. The green and blue colors represent the PEG and Na, respectively.In contrast to the unbound system, when examining the bound system PEG-drug, we observe that no separation occurs after 100 ns of molecular dynamics simulations (Fig. 5). This striking difference between the unbound and bound systems is attributed to the formation of chemical bonds between small molecules and PEG molecules. The establishment of these strong chemical interactions significantly enhances the overall interaction between the components, effectively preventing their dissociation from one another.Figure 5Representative snapshots of the last MD frame for the bound systems after 100 ns of simulations, demonstrating no separation between the PEG carrier and the bioactive drug molecules due to the formation of chemical bonds. The green and blue colors represent the PEG and Na, respectively.Solvent accessible surface areaThe solvent-accessible surface area (SASA) analysis is one of the most important properties of the PEG and small drug molecules interaction23. We calculated the SASA for each sampled configuration of both unbound and bound systems. By determining the portion of a molecule’s surface area covered by PEG, we gained insights into their dynamic molecular mechanisms. The SASA calculations included assessing total surface areas and interactions with the PEG polymer, considered as a part of the solvent. Combining results from all sampled configurations in histograms (Fig. 6) revealed trends and patterns on how these compounds interact during the simulations.Figure 6Histograms of the fraction of drug surface covered by PEG for unbound (a) and bound (b) systems.In the unbound systems, we observed that there was generally no significant interaction between the PEG carrier and bioactive drug molecules from Celastrus hindsii for most of the simulation time. Consequently, the coverage by PEG over these compounds remained at very low portion (< 0.2) for a majority of the sampling configurations (See Fig. 6a). However, a distinct peak is observed for PEG coverage at a very low frequency. This indicates that during the MD simulation, the PEG carrier and drug molecule only interact in a very small fraction of simulation time. Furthermore, we found that the PEG coverage for the unbound systems involving PEG.HINA and PEG.HINB were approximately 0.65, while those with PEG.MAYA and PEG.CELB exhibited coverage values of around 0.55. This difference in PEG coverage can be attributed to the size disparities between HINA/HINB and MAYA/CELB molecules. Smaller molecules tend to have a higher coverage percentages with the PEG carrier, which may influence their overall behavior within the PEG-based drug delivery system.In constranst, bound systems displayed a significantly different behavior compared to unbound systems, with PEG effectively covering the bioactive drug molecules from Celastrus hindsii most of the time. The absence of low portion peak and a wider PEG coverage range (0.1–0.9) in the histogram data confirmed that chemical bonds favored the interaction between PEG and drug molecules in solutions (See Fig. 6b). In the unbound system, the highest frequency is observed around 0.0 for PEG coverage, indicating that the unbound compounds tend to be separated. However, the bound system exhibits a peak of high frequency at approximately 0.5, suggesting that in the bound state, the compounds are more closely associated, leading to higher PEG coverage values. This observation highlights the potential benefits of PEGylation for enhanced PEG coverage and for protection of the small drug molecules.To further investigate the conformational behaviour of the bound compounds, additional calculations on another conformer of the bound systems, PEG-drug 2, were performed (see SI). The results (Table S1) show that the xTB-calculated relative energy of the PEG-drug 2 system is higher than that of the PEG-drug system. Due to the high computational cost of the DFT method, we did not perform the DFT calculations for the PEG-drug 2 systems. The SASA analysis from the MD simulations of the PEG-drug 2 conformer (See SI) exhibits a similar pattern as observed for the PEG-drug systems. This indicates that the drug molecule in this conformation also experiences increased PEG coverage when bound to the PEG carrier.Radius of gyrationWe also calculated the radius of gyration (Rg) for PEG molecules in both unbound and bound configurations to supplement the SASA analysis. The results demonstrated that there was no significant difference observed between the Rg values for PEG in solution and those encountered within the unbound system simulations (see Fig. 7). These findings were consistent with previous molecular dynamics studies on similar PEG-based systems23, further validating the accuracy of the methodology. The average Rg value obtained for PEG molecules across all unbound systems was approximately 0.85 nm, which in excellent agreement with previous study23.Figure 7x, y, z components of the radius of gyration of PEG molecule in both unbound (PEG.drug) and bound (PEG-drug) systems.In contrast to the findings for unbound systems, we observed a reduction in the Rg values (around 0.77 nm) for PEG carrier when bound with bioactive drug molecules. This observation suggests that an attractive nonbonded interaction may be occurring between the PEG and these compounds, causing the PEG polymer reducing the Rg.PEG and drug interactionIn order to gain further insights into the complex interplay between polyethylene glycool (PEG) and drug molecules during MD simulations, we focused our analysis on characterizing van der Waals (vdW) interactions. These non-covalenent forces are known to be critical in mediating molecular recognition and binding processes within many bio-macromolecular systems. To systematically quantify the vdW interactions between PEG and drug molecules, we adopted a rigorous distance-based approach. Specifically, any atomic pair with an interatomic distance of less than 3.0 Å was deemed to exhibit vdW interactions. This well-defined criterion enabled us to objectively assess the strength and distribution of vdW interactions in both unbound and bound complexes.As depicted in Fig. 8, our analysis reveals a striking disparity in the frequency and distribution of vdW interactions between the two systems. Notably, in the case of unbound PEG-drug pairs, the number of vdW interactions remains remarkably low throughout all simulated trajectories. This indicates that the complexation between PEG and drug molecules is not stable, as they tend to separate in solution.Figure 8Number of vdW interaction between PEG and drug molecules in the unbound systems (a) and bound system (b).In contrast to the unbound systems, our MD simulations revealed a significant difference in vdW interactions for the bound complexes. Specifically, we observed that the number of vdW interactions between PEG and bound drug molecules was nearly four times higher compared to the unbound systems. This substantial increase in vdW interaction frequency suggests that the presence of a chemical bonding between PEG-drug significantly enhances their binding affinity. Such an increase in binding affinity can have profound implications for the stability and efficacy of PEGylated therapeutics, particularly when exposed to physiological environments.Our comprehensive analysis of the PEG-drug complex revealed a striking absence of hydrogen bonds between these two molecules. This finding underscores the unique chemical properties of PEG, including its hydrophilic and flexible nature, which enable it to form strong non-covalent bonds with drug molecules through vdW forces. Our MD simulations demonstrate that van der Waals interactions dominantly mediate PEG-drug recognition and binding, with a substantial increase observed upon formation of a chemical bond.DiscussionIn this section, we will compare and contrast the performance of the xTB and DFT methods in the context of investigating PEG-drug interactions. The xTB approach is known for its computational efficiency and lower accuracy compared to DFT33. Despite its limitations, xTB can serve as a useful tool for providing preliminary insights into the behavior of complex systems. Here, we aim to assess the level of agreement between the results obtained from xTB and DFT calculations for PEG-drug interactions. This comparison will help us understand the reliability of xTB in predicting the binding energy and solvation free energy of these complexes.As depicted in Fig. 9, we observed that the xTB method can capture the trend of the DFT calculations regarding the binding energy between PEG and drug molecules. However, the DFT-calculated binding energies were more negative than the corresponding xTB-calculated binding energies, which suggests a stronger binding interaction between the PEG and drug molecules in the case of DFT calculations. While it would be ideal to compare these values with experimental data, to the best of our knowledge, there are no available experimental data or other DFT data for the binding energy of these specific complexes. To advance the field and address this gap, we propose that future experimental work should focus on measuring the binding free energy for the unbound PEG.drug complex. This would not only confirm the computational trends but also provide a critical benchmark for the accuracy of both DFT and xTB methods in predicting such interactions. Furthermore, extending the scope of future experimental studies to include the energies associated with the PEGylation reactions themselves would be highly beneficial.Figure 9Comparison of calculated values obtained from DFT and xTB calculations for binding energy (a) and for delta solvation free energy (b).Our MD simulations have provided valuable insights into the interplay between PEG and a variety of small drug molecules, comparable to the work established by Li et al.23. The consistency between our MD simulation results and those from the Li et al.23 study is remarkable. Both sets of simulations indicate that PEG exhibits an unfavorable interaction with small drug molecules like paclitaxel and piroxicam. This unfavorability manifests as a propensity for the drug and PEG moieties to separate over time within the simulation environment, suggesting that such small drug and PEG complexation may not be stable in physiological conditions. However, the interaction dynamics change markedly when a larger drug molecule with a porphin ring, such as hematoporphyrin. In contrast to the behavior observed with smaller drugs, MD simulations reveal a strong and stable interaction between PEG and the hematoporphyrin molecule23. This notable difference in interaction strength suggests that the presence of a porphin ring significantly influences the binding affinity between PEG and the drug molecule. The interaction observed with hematoporphyrin can be attributed to favorable lipophilic interactions between the nonpolar –\(\hbox {CH}_{2}\)– groups of the PEG and the porphin ring23.The delta solvation free energy is a critical component of the overall binding free energy for PEG-drug complexes. The analysis showed that the level of agreement between the xTB and DFT methods was notably weaker when examining the delta solvation free energies, as demonstrated in Fig. 9b. The linear regression’s coefficient of determination (R\(^2\)) for this comparison is only 0.13, which indicates a very low correlation between these methods. This discrepancy can be attributed to the less accurate auxiliary basis method employed in xTB and the ALPB model when compared to more sophisticated solvation techniques available in VASP software for DFT-based calculations. This observation underscores the importance of considering both binding and solvation free energies in evaluating computational methods for predicting PEG-drug interactions, as well as highlights the limitations of xTB when it comes to solvation free energy calculations relative to DFT. However, it is important to note that despite its limitations, xTB remains a valuable tool for providing preliminary insights into large and complex systems due to its computational efficiency.While both the xTB and DFT methods, when employed using static calculations, provide valuable insights into the structural properties and binding affinity of PEG.drug complexes, they do not account for the dynamic behavior of these systems at finite temperatures or their explicit interactions with water solvent molecules. This limitation can result in an incomplete picture of the forces governing the formation and stability of these complexes. The MD approach, on the other hand, offers a more dynamic perspective by simulating the system’s behavior at finite temperatures and accounting for the explicit interactions with solvent molecules. This enabled us to observe how the intricate interactions between water and PEG-drug systems influence their overall properties. For instance, we discovered that the majority of the time, there is no interaction between PEG and drug molecules in unbound configurations, which aligns well with the less favorable binding energies calculated using DFT simulations.While the xTB method stands out among semi-empirical quantum chemical methods (such as PM6, AM1, and SM6) for its enhanced accuracy and computational efficiency, it is crucial to recognize the inherent limitations of this approach33. The xTB operates on a framework that incorporates empirical parameters and theoretical approximations to simulate electron correlation effects to strike a balance between predictive power and computational demand. Despite these improvements, xTB is not without its constraints; it inevitably falls short of the precision achievable through experimental observations or more sophisticated ab initio methods based on DFT. The semi-empirical nature of xTB implies that while it can reproduce many properties of molecules and materials with remarkable accuracy, particularly when compared to other methods in its class, it does not capture every nuance of the electron dynamics that would be described by more rigorous theoretical frameworks33. DFT, for instance, relies on accurate exchange-correlation functionals to model the electron density and typically provides results that are in closer agreement with experimental findings than semi-empirical approaches.The MD simulations align well with DFT and xTB results, particularly when considering the separation of unbound PEG-drug complexes. The unfavorable binding free energy, as calculated by both the xTB and DFT methods, indicated a high likelihood that these complexes would remain unbound in solution. This prediction is consistent with the observed MD simulation results, which revealed that drug molecules maintained a considerable distance from PEG carrier molecules even under physiological conditions. These findings collectively support the notion that the unfavorable binding energetics contribute to the separation of these complexes, thereby limiting their formation in solution.In our study, we observed a substantial increase in vdW interactions between PEG-drugs upon the formation of chemical bonds. This enhancement in vdW interactions contributes significantly to their overall binding affinity, ultimately leading to improved PEG coverage and SASA. This indicates that the presence of a covalent bond between PEG and drug molecules plays a crucial role in maintaining close proximity during MD simulations. This stable connection prevents the PEG and drug molecules from separating during simulation time-steps, thereby ensuring consistent vdW interactions and reinforcing their binding affinity. In contrast, when no covalent bond is present, the PEG and drug molecules are more likely to fluctuate and separate throughout the MD simulations, leading to weaker vdW interactions and reduced SASA and PEG coverage.It is important to note that while the MD simulations offer valuable insights into the dynamic behavior of PEG-drug complexes, they also come with their own limitations, such as the need for accurate force fields and the computational cost associated with long simulations21. Therefore, it is essential to carefully consider the appropriate level of theory and simulation methods when studying these systems. In the context of implicit solvent modeling, Wang et al. have recently introduced a novel computational method that aims to refine the prediction of binding energies through a sophisticated approach involving a three-trajectory realization at the energy minimization endpoint34. The potential limitations of the current xTB approach for solvation energy calculations could also be improved by employment of an alternative method grounded in statistical mechanics calculations35. These approaches should be included in the future simulations studies of the interaction between PEG and drug molecules.
