Radial of distribution function (RDF) and Kirkwood-Buff Integrals (KBI)The radial distribution function (RDF), g(r), describes the probability of finding a particle at a specific distance r from a reference particle and is given by:$$g\left( r \right) = \frac{V}{{4N\pi r^{2} \Delta r}}\mathop \sum \limits_{i} \left\langle {\mathop \sum \limits_{ij}^{N} \delta \left( {r - r_{ij} } \right)} \right\rangle$$(3)where N is the total number of particles, and r is the distance between pairs of particles23. Figure 2 (panels a, and b) illustrates the probability of finding a j particle around an i particle at a distance r. As illustrated in Fig. 2, both RDF values converge toward unity, perfectly mirroring the behavior anticipated in an ideal gas. In an aqueous solution of the drug, the function g(r) starts at zero because it does not account for the reference particle. It then rises to a peak at the distance that defines the first shell of particles surrounding the reference particle, known as the first solvation shell. As the distance increases, g(r) approaches a value of one, represented by the dotted lines. RDFs between TPTD and water molecules were calculated to examine the spatial distribution of solvent around TPTD across the simulated systems. As shown in Fig. 2 (panel a), a distinct solvation shell is observed at approximately 1.85 Å in all systems. The relatively low peak height indicates weak hydration, suggesting the hydrophobic nature of TPTD in water.Fig. 2(a) Radial distribution function (RDF) of water molecules around TPTD for systems containing 1, 5, and 10 TPTD molecules, illustrating the solvation structure and the formation of the first hydration shell. (b) RDF of TPTD molecules around other TPTD molecules for systems with 5 and 10 TPTD molecules, highlighting the intermolecular interactions and aggregation behavior as TPTD concentration increases.Full size imageRDFs between TPTD molecules were also computed to explore drug-drug interactions, as shown in Fig. 2 (panel b). Systems with 5 and 10 TPTD molecules exhibited sharp peaks in the RDFs, suggesting significant intermolecular interactions. The first peak, located at ~1.05 Å, corresponds to close contact between TPTD molecules. With increasing drug concentration, the intensity of this peak also increased, indicating enhanced aggregation. Additionally, a secondary peak at ~2.25 Å reflects intramolecular hydrogen bonding within TPTD. Notably, in the 10-TPTD system, the TPTD–TPTD RDF peak is higher than that for TPTD–water, demonstrating a preference for self-association over solvation. This suggests that TPTD molecules tend to cluster together, displacing water molecules and potentially aligning their hydrophilic regions inward and hydrophobic regions outward.The number of neighboring molecules encompassed within the first solvation shell, surrounding the central particle, is defined as the coordination number, CN (rmin). This parameter can be determined using the following equation:$$CN = 4\pi \rho_{0} \mathop \int \limits_{0}^{{r_{min} }} g\left( r \right) r^{2} dr,$$(4)In this equation, the limits of integration are from 0 to rmin, where rmin represents the first minimum after the primary peak24. As a result, this integral calculates the number of particles within a spherical shell for a specific range. Owing to treating the radial distribution function as a normalized density distribution function, we multiplied it by the total number of water molecules, as well as by the total number of TPTD molecules, to determine the number of water molecule neighbours around TPTD and the number of TPTD neighbours surrounding the central TPTD molecule, respectively. The coordination number for TPTD–TPTD interactions further supports this trend. The coordination number, representing the average number of solvent molecules within the first solvation shell, was calculated by integrating the RDF up to the first minimum. The values were 2070, 1736, and 1452 for systems containing 1, 5, and 10 TPTD molecules, respectively. This trend indicates that as TPTD concentration increases, water molecules are less able to solvate the drug, likely due to increased TPTD–TPTD interactions.For the 5-TPTD and 10-TPTD systems, the coordination numbers were 3 and 6, respectively (Table 1), indicating that the number of neighbouring TPTD molecules doubles as the drug concentration increases. These RDF analyses collectively provide valuable insights into the interplay between solvation and aggregation behavior in TPTD systems under physiological conditions.Table 1 Coordination number (n(r′)) of TPTD molecules around a reference TPTD molecule at the first peak position (r) in the radial distribution function (g(r)) for 5-TPTD and 10-TPTD systems in aqueous solution.Full size tableThe thermodynamic properties of aqueous solutions offer critical insight into molecular interactions. Among the prominent theoretical frameworks, Kirkwood–Buff (KB) theory stands out for linking microscopic molecular distributions to macroscopic thermodynamic behavior in solutions. This theory provides a foundation for analyzing deviations from ideal solution behavior, particularly in aqueous environments25. The KBI, first introduced by Kirkwood and Buff, have become essential tools in the study of solution thermodynamics. For instance, Shimizu employed KB theory to explore the spatial distribution of water and co-solvent molecules around proteins26.KB theory relies on the integration of radial distribution functions (RDFs) within the grand canonical ensemble (μVT). KB integrals can be obtained from experimental data, such as small-angle scattering, or through computational approaches including molecular dynamics (MD) and Monte Carlo (MC) simulations. While these simulations are often conducted in various ensembles—canonical (NVT), isothermal-isobaric (NPT), and microcanonical (NVE)—the RDFs derived from the NPT ensemble closely approximate those from the μVT ensemble at a defined cutoff distance (Rc). A known limitation in calculating KB integrals through simulations is the finite size of the simulation box, which can prevent accurate convergence to the thermodynamic limit. In large systems, RDFs tend to converge more reliably, but inconsistencies can still arise. To correct for finite-size effects, Vegt and colleagues proposed a method that adjusts RDFs based on the bulk density of component j around component i. This method accounts for the excess or depletion of particles within a sphere of radius r centred on a reference particle. The corrected RDF is expressed as:$$g_{ij}^{correct} \left( r \right) = g_{ij} \left( r \right)\frac{{N_{j} \left( {1 - \frac{{\left( \frac{4}{3} \right)\pi r^{3} }}{V}} \right)}}{{N_{j} \left( {1 - \frac{{\left( \frac{4}{3} \right)\pi r^{3} }}{V}} \right) - \Delta N_{ij} \left( r \right) - \delta_{ij} }},{ }$$(5)Here, \({N}_{j}\) is the number of particles of type j, V is the system volume, \(\Delta {N}_{ij}\left(r\right)\) is the excess number of j particles within a sphere of radius r, and \({\delta }_{ij}\) is the Kronecker delta, which equals 1 if i = j, and 0 otherwise.The general expression for the Kirkwood-Buff integral is:$$G_{ij} \approx 4\pi \mathop \int \limits_{0}^{{R_{c} }} \left[ {g_{ij}^{correct - NPT} \left( r \right) - 1} \right]r^{2} dr,$$(6)In this equation, \({G}_{ij}\) denotes the corrected KB integral under NPT conditions, and \({g}_{ij}^{correct}(r)\) is the radial distribution function between components i and j. Notably, the KB integrals are symmetric, such that \({G}_{ij}\) = Gij27.Table 2 summarizes the KB integrals computed via MD simulations for all studied systems, based on RDFs integrated up to 20 Å and varying TPTD concentrations. The values of GTPTD–TPTD and GWater–Water generally increase with the number of TPTD molecules in the system. The positive values of GTPTD–TPTD in systems with 5 and 10 TPTD molecules indicate favourable, attractive interactions among drug molecules. Similarly, the positive GWater–Water values confirm strong water–water interactions.Table 2 Kirkwood–Buff integrals (KBI), Gij and Gii, calculated from MD simulations for the interactions between TPTD and water (GTPTD–Water), TPTD and TPTD (GTPTD–TPTD), and water and water (G Water–Water) in the 1-TPTD, 5-TPTD, and 10-TPTD systems. KBI values are reported in cm3/mol, with corresponding standard deviations.Full size tableIn contrast, GTPTD–Water exhibits negative values across all systems, indicating a depletion of water molecules in the vicinity of TPTD and suggesting repulsive interactions between drug and water molecules. As the concentration of TPTD increases, the GTPTD–Water values become more negative, reflecting a reduced affinity between the drug and solvent molecules. This trend emphasizes that the self-association of TPTD molecules becomes more favourable than their interaction with water at higher concentrations.The strong positive values of GTPTD–TPTD in the 10-TPTD system further highlight the pronounced aggregation behavior of TPTD. These results suggest that TPTD molecules preferentially associate with each other, excluding water and forming a more tightly bound drug network, a critical insight into their aggregation and solubility behavior in aqueous environments.Root mean square deviation (RMSD)The root mean square deviation (RMSD) analysis of the Cα atoms in TPTD was performed to assess the structural stability of the drug molecules over time and to investigate how varying concentrations of TPTD influence this stability. The RMSD was calculated using Eq. (7):$$RMSD = \sqrt {\frac{1}{{N_{atoms} }}\mathop \sum \limits_{i}^{N} \left( {r_{i} \left( {t_{x} } \right) - r_{i} \left( {t_{ref} } \right)} \right)^{2} } { },$$(7)Here, \({N}_{atoms}\) is the total number of atoms, \({r}_{i}\left({t}_{x}\right)\) is the position of atom i at time t, and \({r}_{i}({t}_{ref})\) is its position in the reference structure, which was taken from the initial frame of the molecular dynamics (MD) trajectory28. In RMSD calculations, the atom positions used for evaluation are determined by the selected frame. In this study, the reference structure was the NPT-equilibrated structure of TPTD, specifically the structure at the first frame of the trajectory simulation, which corresponds to time step 0. This selection provides a clear point of comparison for assessing molecular movements and transformations29.Initially, RMSD values were high in all systems due to structural adjustments at the start of the simulation. However, fluctuations decreased over time and stabilized after 90 ns. The higher RMSD observed in the 1-TPTD system suggests increased conformational flexibility due to extensive interactions with surrounding water molecules. Conversely, the reduced RMSD in the 5-TPTD and 10-TPTD systems indicates increased structural stability, likely due to drug–drug aggregation limiting molecular motion.These findings suggest that TPTD aggregation enhances structural rigidity, while fewer drug molecules allow TPTD to interact more freely with water, leading to greater flexibility. Aggregation among drug molecules may decrease conformational entropy and contribute to a more stable molecular configuration during the simulation.Figure 3 presents the RMSD profiles of TPTD molecules over a 100 ns imulation period for all studied systems. Equilibration was achieved around 90 ns across all simulations. The average RMSD for a single TPTD molecule in the 1-TPTD system between 90 and 100 ns was 8.31 Å. In comparison, the RMSD values for systems with 5 and 10 TPTD molecules were lower, at 5.87 Å and 2.50 Å, respectively.Fig. 3(a) Root mean square deviation (RMSD) of the Cα atoms of TPTD molecules over a 100 ns molecular dynamics simulation for systems containing 1, 5, and 10 TPTD molecules in water, indicating structural deviations over time. (b) Average RMSD values for each system, highlighting the decreasing trend in structural fluctuations with increasing TPTD concentrationFull size imageRoot mean square fluctuation (RMSF)The root mean square fluctuation (RMSF) is a valuable metric for assessing the flexibility and average positional displacement of individual atoms or groups of atoms—such as amino acid residues—within a molecular structure over time. Figure 4 illustrates the flexibility profile of TPTD residues across different drug concentrations, with peak heights corresponding to the magnitude of fluctuation for each amino acid.Fig. 4(a) Root mean square fluctuation (RMSF) analysis of TPTD residues during molecular dynamics simulations. (b) Average RMSF values with standard deviations for all systems, showing a decreasing trend in residue flexibility with increasing TPTD concentration, suggesting enhanced structural stability due to drug aggregation.Full size imageAs shown in Fig. 4 (panel a), the RMSF values in the 1-TPTD system range from 25 to 33 Å, with Valine 31 and Histidine 32 exhibiting significantly higher fluctuations compared to other residues. These elevated values suggest that these residues are more solvent-exposed and possess higher conformational mobility in the aqueous environment.In the systems containing higher concentrations of TPTD, residue flexibility is notably reduced. Specifically, RMSF values for the 5-TPTD and 10-TPTD systems fall within the ranges of 7.26–21.47 Å and 5.68–14.59 Å, respectively. The residues exhibiting the lowest fluctuations in the 5-TPTD system include Asparagine 78, Glutamic Acid 90, Tryptophan 91, Leucine 92, Arginine 93, and Lysine 94. Similarly, the most stable residues in the 10-TPTD system are Glutamine 74, Leucine 75, Asparagine 78, Leucine 219, Asparagine 220, and Lysine 319. These residues are likely involved in stabilizing intermolecular interactions during the aggregation process.Panel b of Fig. 4 presents the average RMSF values across all residues for each system. The trend reveals a decrease in overall residue flexibility with an increasing number of drug molecules. This suggests that spatial hindrance caused by drug–drug interactions reduces the freedom of motion for individual residues, thereby enhancing structural stability in more concentrated systems.The water residence time at the TPTD surface and dipole momentsThe residence time of water refers to the duration that a water molecule remains on the surface of the TPTD during a single stay. To calculate the residence time of water on the drug molecule, we utilize the time correlation function method, described as follows:$${\text{p}}\left( {\uptau } \right) = \mathop \int \limits_{0}^{\infty } \delta (N\left( t \right) - \left( {N\left( {t + \tau } \right)} \right)dt,$$(8)In this equation, the function δ(x–y) equals 1 when x = y (with the condition that x ≠ 0 and y ≠ 0) and equals 0 when either x ≠ y or both x and y equal 0, indicating that the site is unoccupied. In this context, N(t) represents the index of the water molecule located at the hydration site at time t. The site is considered occupied if a water molecule is found within a spherical volume of radius 3.5 Å centered at the coordinates of the site. The resulting time correlation function is fitted using a double exponential model. Often, a Kohlrausch–Williams–Watts stretched exponential is employed, which is expressed as:$${\text{p}}\left( {\uptau } \right) = W_{0} \left[ {\left( {1 - w} \right)\exp \left( { - k_{1} \tau } \right) + w {\text{exp}}\left( { - k_{2} \tau } \right)} \right],$$(9)In this context, \({W}_{0}\) represents the average site occupancy, defined as the fraction of time that any water molecule occupies a site. The parameter of \({k}_{1}\) and \({k}_{2 }\) are the rate constants, while \(\left(1-w\right)\), and \(w\) denote the weights of the fast and the slow components, respectively. The short and the long residence times are defined as \({\tau }_{1}=1/ {k}_{1}\) and \({\tau }_{2}=1/ {k}_{2}\)30.Figure 5 (panel a) shows the water molecules present at the surface of TPTD in all the studied systems. Notably, an increase in TPTD concentrations leads to a decrease in the number of water molecules on the drug’s surface over time. The average residence times are t1 = 0.002 ps, t1 = 0.0006 ps, and t1 = 0.0003 ps. These calculations indicate that water molecules near the hydrophilic regions of the TPTD molecule have significantly longer residence times compared to those near the hydrophobic region. This suggests the formation of more stable water structures near the hydrophilic regions. This phenomenon is likely due to water molecules at the 1-TPTD surface forming stronger hydrogen bonds with the amino acids of the TPTD molecule, thereby establishing a stable network of water molecules near the hydrophilic region.Fig. 5(a) The number of water molecules at distribution around TPTD surface. (b) Spatial distribution of water molecules within 5 Å of TPTD molecules in the 5-TPTD system after 100 ns of simulation. (c) Distribution of water molecules around TPTD in the 10-TPTD system after 100 ns. Magenta arrows represent the orientation of water dipole moments, suggesting increasing repulsion and reduced hydration as TPTD concentration increases.Full size imageTo further investigate the orientation and interaction tendencies of TPTD, dipole moments of the drug molecules were calculated using the Dipole Watcher plugin in VMD21 over the last 10 ns of simulation. The normalized dipole moments were 202.62 D ± 0.031, 58.36 D ± 0.015, and 32.04 D ± 0.086 Debye for the 1-TPTD, 5-TPTD, and 10-TPTD systems, respectively. In the 1-TPTD system, water remains closely associated with the drug throughout the entire 100 ns simulation, suggesting a preference for drug–solvent interaction.However, as TPTD concentration increases, aggregation occurs rapidly due to enhanced intermolecular dipole–dipole interactions. These interactions contribute to the stabilization of TPTD aggregates, pushing water molecules away from the drug core. Visualizations (Fig. 5, panels b and c) indicate that the dipole moment vectors of water molecules orient outward from the drug clusters, supporting the notion that higher dipole interactions among TPTD molecules drive aggregation by repelling solvent molecules.Normalized radius of gyrationThe radius of gyration (\({R}_{g}\)) is a key structural parameter that reflects the compactness of a molecule. It represents the mass-weighted root mean square distance of a system’s atoms from its center of mass and is calculated using the following equation:$$R_{g} = \left( {\frac{{\mathop \sum \nolimits_{i} r_{i}^{2} m_{i} }}{{\mathop \sum \nolimits_{i} m_{i} }}} \right)^{\frac{1}{2}} ,{ }$$(10)where \({m}_{i}\) is the mass of atom i, \({r}_{i}\) is the distance of atom i from the center of mass, and Rg denotes the radius of gyration of the molecule31.Figure 6 (panel a) illustrates the time evolution of the Rg values for TPTD molecules in systems containing 1, 5, and 10 drug molecules. The corresponding initial and final Rg values are presented in Table 3. Structural snapshots from the final frame of the simulations (t = 100 ns) for the 5-TPTD and 10-TPTD systems are shown in Fig. 6 (panels b and c).Fig. 6(a) Radius of gyration (Rg) values over the 100 ns molecular dynamics simulation for systems containing 1, 5, and 10 TPTD molecules in water, illustrating the effect of increasing drug concentration on molecular compactness. (b,c) Representative snapshots at 100 ns of the 5-TPTD (b) and 10-TPTD (c) systems, showing the spatial distribution and compactness of drug molecules in the aqueous environment.Full size imageTable 3 Calculated radius of gyration (Rg) values for TPTD molecules in different systems (1-TPTD, 5-TPTD, and 10-TPTD) at the initial and final time points of the molecular dynamics simulation.Full size tableFor the 1-TPTD system, the Rg decreases from 12.35 Å at the beginning of the simulation to 10.22 Å at the end, indicating a slight compaction of the structure over time. In the 5-TPTD system, the Rg increases from 5.46 Å initially to 8.32 Å at the end, suggesting structural reorganization or expansion during aggregation. For the 10-TPTD system, Rg starts at 3.55 Å and gradually increases to 4.37 Å, remaining relatively stable, which indicates that the system maintains a more compact conformation throughout the simulation.Among the three systems, the 1-TPTD system consistently exhibits the highest Rg values, while the 10-TPTD system shows the lowest, reflecting greater molecular compactness due to stronger aggregation. The smaller difference between initial and final Rg values in the 10-TPTD system further supports the idea of a stable, aggregated structure.The mean and standard deviation of the Rg values for 1-TPTD, 5-TPTD, and 10-TPTD are 10.66 Å ± 0.58, 7.07 Å ± 0.98, and 4.33 Å ± 0.067, respectively. The highest average Rg and greatest fluctuations are observed in the 1-TPTD system, suggesting more pronounced structural flexibility. This behavior likely results from the single drug molecule extending toward surrounding water molecules to maximize hydrogen bonding interactions. Conversely, in the systems containing multiple TPTD molecules, the average Rg values decrease, indicating enhanced compactness due to aggregation and reduced flexibility in the aqueous environment.Hydrogen bond (H-bond)Hydrogen bonding plays a critical role in stabilizing molecular aggregates, as it facilitates strong intermolecular interactions between drug molecules. When two molecules with suitable donor and acceptor groups come within close proximity, a hydrogen bond may form—typically when the donor–acceptor distance is less than 3.5 Å and the donor–hydrogen–acceptor angle exceeds 150°32.The stability of protein structures is often correlated with the number of hydrogen bonds; a greater number generally indicates higher structural stability. In this study, the number of intermolecular hydrogen bonds among TPTD molecules was analyzed over the final 10 ns of the molecular dynamics simulation using the VMD software. The average number of hydrogen bonds formed was modeled using a Gaussian distribution function:$$F\left( x \right) = \frac{a}{{\sigma \sqrt {2\pi } }}exp\frac{{ - (x - \overline{x})^{2} }}{{2\sigma^{2} }},{ }$$(11)Here, \(\overline{x }\) is the average number of hydrogen bonds, \(\sigma\) is the standard deviation, and \(a\) is an adjustable amplitude parameter. The fitting was performed using the GNUPLOT software33.Figure 7 (panel d) presents a summary of the average number of hydrogen bonds observed across all systems studied. The data reveal that as the concentration of TPTD molecules in aqueous solution increases, there is a corresponding rise in the number of hydrogen bonds formed between water and TPTD. Specifically, the average number of hydrogen bonds in the 5-TPTD system is recorded at 938.80 ± 0.135, while the 10-TPTD system shows a significantly higher average of 2585.28 ± 0.156. Additionally, the analysis indicates that a greater concentration of TPTD in water leads to a notable decrease in the normalized number of hydrogen bonds. This reduction from the 1-TPTD to the 10-TPTD system implies that enhanced intermolecular interactions among TPTD molecules arise when the hydrophobic regions of TPTD interact more insignificantly with the surrounding solvent molecules. Figure 7 (panel a) shows the Gaussian distribution of normalized intermolecular hydrogen bonds with a slight reduction in bond formation in the 10-TPTD system, indicating a decreased tendency to form hydrogen bonds for any TPTD in this system. As noted in the radial distribution function (RDF) analysis (Fig. 2), the number of solvent molecules in the first solvation shell around TPTD decreases with increasing drug concentration, while the intensity of the RDF peaks between TPTD molecules increases. This shift reflects reduced drug–solvent interactions and enhanced drug–drug aggregation, ultimately reducing solubility.Fig. 7(a) Gaussian distribution of the normalized number of hydrogen bonds formed between TPTD molecules in the 5-TPTD and 10-TPTD systems. (b) Occupancy percentage of hydrogen bonds between TPTD and surrounding water molecules. (c) Occupancy percentage of hydrogen bonds between TPTD molecules. (d) Average number of hydrogen bonds formed between TPTD and water molecules (yellow bars), and between TPTD molecules (blue bars) across 1-, 5-, and 10-TPTD systems.Full size imageHydrogen bond occupancy—defined as the percentage of simulation time during which a particular hydrogen bond exists—was also calculated to assess the stability of these interactions. As shown in Fig. 7 (panel b), the 1-TPTD system exhibited the highest occupancy values, confirming strong and persistent hydrogen bonding with water. In contrast, the 5-TPTD and 10-TPTD systems showed lower occupancy values, beginning at only 10% and 30% of frames, respectively. These results reinforce that higher drug concentrations reduce solvent interaction and hydrogen bond stability.The occupancy of hydrogen bonds between drug molecules themselves is presented in Fig. 7 (panel c). The 10-TPTD system shows the highest occupancy, indicating more stable drug–drug hydrogen bonding and a more compact aggregate structure. These intramolecular hydrogen bonds contribute significantly to the stabilization of TPTD aggregates in aqueous solution. Specific residues and atoms involved in TPTD–TPTD hydrogen bonding are listed in Table S2 (Supporting Information).Solvent accessible surface areaThe solvent-accessible surface area (SASA) of proteins is a key factor in understanding their structural behavior in a solvent environment. SASA refers to the surface area of a molecule that is exposed to a solvent34, which is typically calculated using the ‘rolling ball’ algorithm. This algorithm involves a sphere with a radius equivalent to the size of a solvent molecule (1.4 Å for water molecules) to define the molecule’s surface.Figure 8 (panel a) illustrates the SASA values throughout the simulation time, while panel b shows the average SASA values for TPTD molecules. A comparison of the SASA profiles reveals that the 1-TPTD system exhibits higher SASA values compared to higher concentration solutions. This observation can be attributed to the molecular structure and size of TPTD. A single TPTD molecule has a relatively low molecular mass compared to systems with 5 or 10 TPTD molecules, leading to greater solvation and, thus, higher SASA values for the 1-TPTD system (3494.82 Å2 ± 0.153). In contrast, the average SASA values per TPTD molecule for the 5-TPTD and 10-TPTD systems are lower (3042.53 Å2 ± 0.039 and 2980.20 Å2 ± 0.095, respectively), reflecting the increased molecular mass in these higher concentration systems.Fig. 8(a) Time evolution of the solvent-accessible surface area (SASA) for 1-, 5-, and 10-TPTD systems per TPTD molecule in aqueous solution at 310 K. Data are color-coded as: blue (1-TPTD), red (5-TPTD), and green (10-TPTD). (b) The SASA values per TPTD molecule with standard deviations for each system calculated over the last 10 ns of the molecular dynamics simulation.Full size imageOur findings indicate a general decrease in SASA values by the end of the simulation across all systems. This decrease is likely due to a slight disruption in the hydrogen bonding network among water molecules when TPTD molecules are introduced.Specifically, the more significant reduction in SASA values in the 5- and 10-TPTD systems, compared to the 1-TPTD system, is linked to the decreased polarity of these systems as additional drug molecules are added. As discussed in the dipole section, 1-TPTD has the highest average dipole moment, reinforcing its polar character.There is a strong correlation between hydrogen bonding and SASA. Changes in the number of hydrogen bonds lead to variations in the accessible surface area of TPTD molecules in water. The increase in SASA for the 1-TPTD system is associated with a higher number of hydrogen bonds, which allows the drug to attract more water molecules due to its polar nature. However, as the concentration of TPTD molecules increases (5 and 10 TPTDs), this ability to attract water molecules diminishes, resulting in lower SASA values. Consequently, the decrease in SASA values (Fig. 9, panel b) suggests that the addition of more TPTD molecules reduces the drug’s hydrophilicity and causes a more compact structure.Fig. 9Electrostatic (Elec), van der Waals (vdW), and total non-bonded interaction energies between TPTD and water molecules (1-, 5-, and 10-TPTD systems), and between drug molecules (5-TPTD and 10-TPTD systems). Energies were calculated using the NAMD Energy plugin in VMD. Error bars represent standard deviations from the last 10 ns of the molecular dynamics simulations.Full size imageInteraction energyThe non-bonded interaction energy between water and TPTD molecules during the final 10 ns of the simulation was calculated using the “NAMD Energy” plugin35 within the VMD package. This analysis included electrostatic (elec) and van der Waals (vdW) components, which together constitute the total non-bonded interaction energy. These energies are defined by the following expression:$$E_{non{\text{-}}bonded} { } = E_{elec} + E_{vdw} = \mathop \sum \limits_{i > j} \frac{{q_{i} q_{j} }}{{\varepsilon r_{ij} }} + 4\varepsilon { }\left[ {\left( {\frac{\sigma }{r}} \right)^{12} - \left( {\frac{\sigma }{r}} \right)^{6} } \right],$$(12)In this equation, \({q}_{i}\) and \({q}_{j}\) are the partial charges on atoms i and j, \({r}_{ij}\) is the distance between them, and \(\epsilon\) is the dielectric constant. For the vdW term, \(\varepsilon\) represents the depth of the potential well, and \(\sigma\) is the interatomic distance at which the potential is zero.Gaussian distributions were fitted to the interaction energy data using the GNUPLOT program for all simulated systems. These fitted distributions are displayed in Fig. S2 (supplementary information), while average interaction energy values (electrostatic, vdW, and total non-bonded) are summarized in Table 4 and further visualized in Fig. 9. These results highlight the dominant role of electrostatic interactions in the overall interaction between water and TPTD molecules.Table 4 Average electrostatic (Elec), van der Waals (VDW), and total interaction energies (in kcal/mol) between TPTD and water, as well as TPTD–TPTD interactions, during the final 10 ns of molecular dynamics simulations for the 1-TPTD, 5-TPTD, and 10-TPTD systems. Values are presented as mean ± standard deviation.Full size tableAs shown in Table 4, electrostatic interactions contribute more significantly than vdW forces across all systems. Notably, the total non-bonded interaction energy between TPTD and water decreases as the concentration of TPTD increases. This reduction is primarily attributed to a decline in electrostatic interactions. In the 1-TPTD system, the interaction between the drug and water is stronger, but as more drug molecules are introduced, the COM distances between water and TPTD molecules increase (as discussed in the COM distance section). This spatial separation leads to weakened electrostatic interactions and a corresponding drop in non-bonded energy, indicating that water molecules are repelled or displaced as TPTD concentration rises.Further insights are provided in Fig. S2 (panels a and b), which show that electrostatic interactions between TPTD molecules—due to their charged amino acid residues—contribute substantially to the total interaction energy. In fact, electrostatic energy is higher in the 10-TPTD system compared to the 5-TPTD system. Similarly, vdW interaction energy reaches its maximum value in the 10-TPTD system, correlating with increased TPTD concentration in the aqueous environment. As the number of drug molecules increases, they tend to cluster together, resulting in enhanced charge distribution and stronger intermolecular interactions.The short-range vdW interactions also influence the compactness of the TPTD systems. This is supported by the low radius of gyration observed in the 10-TPTD system (Fig. 6), indicating a more compact molecular structure due to enhanced vdW attractions. These interaction energies facilitate the formation of intermolecular bonds between TPTD molecules, emphasizing the critical role of concentration in modulating molecular behavior. At higher concentrations, TPTD molecules exhibit stronger mutual interactions, underscoring the impact of drug aggregation in aqueous environments.Secondary structures of TPTDsThe three-dimensional structure of proteins provides essential insights into their conformational behavior and potential for aggregation. In this study, structural changes in TPTD molecules were assessed in relation to their secondary structure as drug concentration increased. The influence of drug addition on TPTD’s secondary structure was systematically analyzed across all simulated systems, and the probabilities of each secondary structure type were calculated. These results are summarized in Table 5 and visualized in Fig. 10.Table 5 Average secondary structure content (310-helix, α-helix, coil, turn, and isolated-bridge) of TPTD in systems containing 1, 5, and 10 molecules, calculated from the last 10 ns of MD trajectories.Full size tableFig. 10Percentage distribution of secondary structure elements (turn, α-helix, 3₁₀-helix, coil, and isolated-bridge) in TPTD molecules over the last 10 ns of MD simulations for systems containing 1, 5, and 10 TPTD molecules in water. Structural types are represented in red (turn), green (α-helix), pink (3₁₀-helix), blue (coil), and violet (isolated-bridge).Full size imageThe coil content in the 1-TPTD system was found to be 18.38%. However, as more drug molecules were introduced, the coil content decreased to 18.15% in the 5-TPTD system and further to 14.51% in the 10-TPTD system (see Table 5). This reduction suggests that intermolecular interactions between drug molecules reduce the flexibility of the protein structure, causing portions of the coil to transition into other structures such as turns, 310-helix, and isolated bridges. Additionally, reduced coil content may result from limited solvent accessibility as TPTD molecules become more compact in higher concentrations.Notably, the α-helix content increased with rising TPTD concentrations. In the 10-TPTD system, α-helices formed rigid, stable regions that likely contribute to the binding and stabilization of TPTD molecules. An increase in α-helix content correlates with a greater probability of intermolecular contact, supporting the notion of structural compactness and potential aggregation.Turn structures also increased as the number of TPTD molecules rose. This is particularly evident in the 10-TPTD system, where the highest turn content was observed. Turns facilitate directional changes in the peptide chain, enabling TPTD to fold upon itself and adopt a more compact, aggregation-prone conformation. These findings are consistent with the results from the radius of gyration analysis, which also indicated increased structural compression with higher peptide concentrations.Interestingly, the 310-helix content decreased in the 5-TPTD system, likely due to changes in the orientation and packing of the drug molecules, rather than their proximity. As the system becomes more compact, some secondary structure elements may be suppressed or restructured.The content of isolated bridges remained low across all systems but showed a slight increase with concentration—from negligible in the 1-TPTD system to 0.22% and 0.37% in the 5- and 10-TPTD systems, respectively. This increase reflects the presence of non-covalent interactions between TPTD molecules, particularly in the highly concentrated 10-TPTD system.Hydrogen bonding plays a crucial role in these structural transformations. As detailed in Fig. 7 (panels a and d), the number of hydrogen bonds between TPTD molecules increases significantly with concentration. These interactions are central to the formation of stable α-helices and turn structures, promoting a more compact and potentially aggregative state.In the 10-TPTD system, the TPTD chains undergo notable conformational changes. The two ends of each molecule are often linked through α-helical and turn structures, forming a continuous surface that is structurally favorable for aggregation. These findings underscore the concentration-dependent conformational behavior of TPTD and its implications for aggregation in aqueous environments.Potential of mean forceGibbs free energy plays a vital role in understanding protein–protein interactions. The free-energy landscape, representing the minimum energy configuration on a protein’s surface, is referred to as the potential of mean force (PMF). PMF is a function of molecular coordinates and typically exhibits a global minimum or several local minima, especially in complex molecular systems. The free energy landscape is defined as:$$\Delta G\left( {x,y} \right) = - k_{B} Tln\left( {\rho \left( {x,y} \right)} \right),{ }$$(13)where x and y correspond to the root mean square deviation (RMSD) and radius of gyration (Rg), respectively; T is the temperature, \(k_{B}\) is Boltzmann’s constant, and ρ(x,y) is the probability distribution of states throughout the simulations36.PMF plots were generated using MolAICAL software37, mapping RMSD versus Rg for the centers of mass of TPTD molecules. The corresponding free-energy landscapes were visualized using Origin 2016 (OriginLab)38. For comparison, Kumar et al. previously demonstrated the stability of amyloid-β-enzyme complexes through a similar approach39.Figure 11 presents both 2D contour maps and 3D surface plots of the Gibbs free-energy landscape for systems containing 1, 5, and 10 TPTD molecules. The free energy ranges from 0 to 4 kcal/mol for 1-TPTD, 0 to 3 kcal/mol for 5-TPTD, and 0 to 4 kcal/mol for 10-TPTD. Notably, the blue regions in the plots (particularly in Fig. 11e) indicate the lowest energy states. These regions become smaller and more centralized in the 10-TPTD system, suggesting increased structural stability.Fig. 11Two-dimensional (2D) contour plots and three-dimensional (3D) surface plots of the Gibbs free-energy landscapes as a function of root mean square deviation (RMSD) and radius of gyration (RoG) for TPTD–water systems: (a,b) 1-TPTD, (c,d) 5-TPTD, and (e,f) 10-TPTD. The energy landscapes were generated from the last 10 ns of molecular dynamics simulations using MolAICAL and visualized with OriginPro 2016 (OriginLab Inc., Northampton, MA, USA). The color scale represents the free energy values in kcal/mol, with blue regions indicating the lowest energy states and thus the most stable conformations. RMSD and RoG serve as reaction coordinates to depict conformational stability and structural compaction of TPTD across different aggregation states.Full size imageAmong the systems studied, the 1-TPTD configuration exhibits the most centralized and deepest energy minima, indicating a highly stable state. During the simulations, TPTD undergoes conformational changes, with the lowest energy configurations corresponding to its most stable folded states. The narrow and deep energy funnels seen in the 3D plots further emphasize the presence of these stable conformations.Overall, the energy landscapes support a more stable folding and aggregation behavior for the 10-TPTD system compared to 5-TPTD. These findings align well with the radius of gyration results, reinforcing the hypothesis that the 10-TPTD system forms a more stable molecular assembly in aqueous environments.Cluster analysis of the molecular dynamics simulation trajectoriesThe size aggregation of the TPTD molecules can be examined by tracking the protein size using the gmx clustsize tool in the GROMACS 2024.440 program. The number of clusters formed and the maximum number of TPTD molecules in each cluster are shown in Fig. 12 (panels a, and b). In the early stages of the simulations, the number of clusters varied between 1 and 3 for the 5-TPTD system and between 2 and 3 for the 10-TPTD system in the first 15 ns. After 15 ns, the molecules began to move apart, which increased the number of clusters. However, after 80 ns, a decrease was observed in the number of clusters, indicating that smaller clusters were merging into larger formations. Figure 12 (panel a) shows that in the 5-TPTD system, the number of clusters ranges from approximately 50–85 ns, while in the 10-TPTD system, it ranges from about 15–60 ns on a larger scale. It is crucial to note that the largest TPTD clusters, particularly those containing 5 and 10 TPTD molecules in the 5-TPTD and 10-TPTD systems, show instability. These clusters have a fleeting existence; they either break apart spontaneously into smaller aggregates or lose one or two TPTD molecules, resulting in abrupt fluctuations in both the quantity and size of the clusters observed throughout the trajectory.Fig. 12Time evolution of the (a) number of clusters, (b) maximum number of TPTD molecules within the clusters.Full size imageFigure 12 (panel b) shows the number of TPTD molecules in the largest clusters for the 5-TPTD and 10-TPTD systems. In the 5-TPTD system, the largest cluster at the end of the simulation consisted of 2 TPTD molecules, indicating that the remains in solution were more dispersed. Conversely, the 10-TPTD system demonstrated a much stronger clustering effect, with the largest clusters containing 7 TPTD molecules, and an additional 2 molecules remaining dispersed in the aqueous solution. These results highlight the significant role of drug concentration in influencing molecular aggregation behavior.