Structural evolution analysis of simulated systems of NIPAM molecules with different chain lengthsTo minimize the effect of factors other than chain length, all five simulation systems constructed used the same simulation box, in which four curcumin molecules are randomly placed, and their relative positions are the same in all five boxes. Subsequently, 60 4-mer NIPAM molecules (NIPAM4) are randomly inserted into one of the boxes. Each of the remaining boxes is subjected to a similar operation, differing, in the number of molecules inserted and the length of the molecular chains, which are, respectively, 30 8-mer NIPAM molecules (NIPAM8), 20 12-mer NIPAM molecules (NIPAM12), 15 16-mer NIPAM molecules (NIPAM16), and 12 20-mer NIPAM molecules (NIPAM20). These numbers of molecules are set to ensure that the molecular mass (concentration) of the solutes in each box is similar. To facilitate differentiation and description, these five systems are abbreviated as CUR-NIPAM4, CUR-NIPAM8, CUR-NIPAM12, CUR-NIPAM16, and CUR-NIPAM20, respectively. The above systems all undergo the same 60 ns time simulation. The root mean square deviation (RMSD) is an important reference quantity for the change and stability of the structure of a substance during the simulation process44. RMSD can be used to characterize the net change in conformation of the NIPAM4, NIPAM8, NIPAM12, NIPAM16, and NIPAM20 molecular systems after overall rotations or translations relative to the equilibrium coordinates, which are calculated and displayed in Fig. 1. From the figure, it can be seen that the five systems with different chain lengths of NIPAM molecules show not exactly the same RMSD variations over time, but the changing trends and the time nodes are almost the same, i.e., they all show drastic changes before 20 ns and then fluctuate around a certain value, and these fluctuations become relatively smooth in the last 10 ns, which suggests that the configurations of the five systems have gradually reached a stable state with the simulation time. In detail, the overall mean value of the RMSD of NIPAM4 over time is relatively high, which is a reflection of the greater structural changes in this molecular system. As the chain length of NIPAM molecules increases, there is a gradual decrease in RMSD. However, during the last 10 ns, this trend becomes less pronounced.Fig. 1: RMSD (root mean square deviation) of the NIPAM (N-isopropylacrylamide) molecular groups in the simulation systems.Variation of RMSD with simulation time for the NIPAM4 (4-mer N-isopropylacrylamide), NIPAM8 (8-mer N-isopropylacrylamide), NIPAM12 (12-mer N-isopropylacrylamide), NIPAM16 (16-mer N-isopropylacrylamide), and NIPAM20 (20-mer N-isopropylacrylamide) groups corresponding to the CUR-NIPAM4 (a system consisting of curcumin molecules and NIPAM4 molecules), CUR-NIPAM8 (a system consisting of curcumin molecules and NIPAM8 molecules), CUR-NIPAM12 (a system consisting of curcumin molecules and NIPAM12 molecules), CUR-NIPAM16 (a system consisting of curcumin molecules and NIPAM16 molecules), and CUR-NIPAM20 (a system consisting of curcumin molecules and NIPAM20 molecules) systems, respectively.Figure 2a depicts the solvent-accessible surface area (SASA) of the NIPAM molecular system during the simulation process. As the simulations are run, the SASA values for each system gradually decrease and remain at a relatively smooth and low level during the last 10 ns of time. The decrease in SASA values serves as compelling evidence for the transition of NIPAM molecules from dispersion to aggregation. The initial values of SASA for different molecular systems in Fig. 2a exhibit notable discrepancies, which can be attributed to differences in their molecular shapes and quantities. Examining the total amount of change in SASA (or the total rate of change over this period of time) from 0 ns to 60 ns basically satisfies the relationship of NIPAM4 > NIPAM8 > NIPAM12 ~ NIPAM16 ~ NIPAM20, which exhibits the extent to which NIPAM systems with different molecular chain lengths change in conformation from a molecular aggregation point of view, which is consistent with the inference of Fig. 1.Fig. 2: SASA (solvent accessible surface area) and Rg (radius of gyration) of the NIPAM (N-isopropylacrylamide) molecules group in the simulated system.a SASA versus simulation time for the NIPAM4 (4-mer N-isopropylacrylamide), NIPAM8 (8-mer N-isopropylacrylamide), NIPAM12 (12-mer N-isopropylacrylamide), NIPAM16 (16-mer N-isopropylacrylamide), and NIPAM20 (20-mer N-isopropylacrylamide) groups corresponding to the CUR-NIPAM4 (a system consisting of curcumin molecules and NIPAM4 molecules), CUR-NIPAM8 (a system consisting of curcumin molecules and NIPAM8 molecules), CUR-NIPAM12 (a system consisting of curcumin molecules and NIPAM12 molecules), CUR-NIPAM16 (a system consisting of curcumin molecules and NIPAM16 molecules), and CUR-NIPAM20 (a system consisting of curcumin molecules and NIPAM20 molecules) systems, respectively. b Data on Rg for the NIPAM molecules group described in (a).The radius of gyration (Rg) is another informative quantity that characterizes the evolution of the structure of a molecular system, and it is calculated by the following formula (defining equation)45:$${R}_{g}={\left(\frac{{\sum }_{i}{\left|{r}_{i}\right|}^{2}{m}_{i}}{\sum _{i}{m}_{i}}\right)}^{\frac{1}{2}}$$
(1)
where \({m}_{i}\) is the mass of atom i and \({r}_{i}\) is the position of atom i relative to the center of mass of the system. The size of Rg is an important response to the bulkiness of the system45. The calculated results of Rg for the NIPAM4, NIPAM8, and NIPAM12 molecular systems can be found in Fig. 2. The values of all three are almost equal in the initial state, which is consistent with the fact that these systems are constructed in simulation boxes with the same dimensions and have a nearly comparable number of NIPAM monomers dispersed in them. As the simulation time increases, these Rg show different degrees of large fluctuations, but the overall trend is all gradually decreasing. This reveals that running the simulations prompts a gradual decrease in the fluffiness (and a gradual increase in the densification) of the NIPAM molecular system, apparently as a result of the gradual coalescence of those randomly dispersed molecules into clusters. One detail is that the NIPAM4 system (and also the NIPAM8 system) has a relatively low Rg, implying its higher degree of compactness, and hence inferring that its molecules undergo a more intense aggregation behavior (more pronounced structural changes), which is consistent with the conclusions drawn from the larger RMSD of NIPAM4 in Fig. 1. However, when continuing to increase the chain length of NIPAM, the Rg of the NIPAM16 system and the NIPAM20 system do not show an overall decrease during the large fluctuations. On the contrary, they gradually increase during the last 15 ns. Accordingly, it can be speculated that although the long-chain NIPAM molecular systems also undergo the process of self-assembling into clusters during the simulation as the short-chain molecular systems do, they have a relatively large number of clusters, and the shapes of these clusters and the relative positions of the clusters are unstable, and these clusters are over-dispersed in the box (each dispersed in different corners of the box), which results in the overall Rg of the whole system to be increased instead. That is, the NIPAM molecular systems tend not to aggregate into a single cluster but rather often form multiple dynamically changing clusters. As illustrated in Supplementary Fig. 2, the evolution of the configuration of the CUR-NIPAM4 system over time aptly demonstrates this process. It is worth noting that the data presented in Fig. 1, Fig. 2a, b are all derived from considering the entire simulated system rather than focusing solely on the formation process of a single cluster. This explains why these data exhibit significant fluctuations (attributed to multiple clusters undergoing dynamical changes) and why different types of data (thus different physical meanings and calculations) do not show the same changing trends (e.g., all are gradually increasing, or all decreases gradually) and the same timing of transitions (e.g., all transitioning at x ns).Supports for the above deductions from Fig. 2 can be found visually in Fig. 3. Figure 3a–e shows snapshots of the configurations of the CUR-NIPAM4, CUR-NIPAM8, CUR-NIPAM12, CUR-NIPAM16, and CUR-NIPAM20 systems at the 0 ns moment (before simulation) and the 60 ns moment (after simulation), respectively, during the simulation. As can be seen from the figure, the CUR and NIPAM molecules in all systems are dispersed in the simulation box before running the simulation, and at the end of the simulation, these molecules are highly aggregated and formed into clusters, which realize the encapsulation of the CUR molecules. This is expected to be seen since such results demonstrate that NIPAM molecules with different chain lengths in the dispersed-in-solvent state all possess the kinetic conditions for self-organizing to form CUR-loaded DDS, with the difference that they exhibit different potentials for self-assembly. Specifically, in Fig. 3a, almost all of the NIPAM4 molecules end up aggregating at the same site, forming the only cluster affixed to the periphery of the four CUR molecules, demonstrating the relatively highest degree of molecular aggregation, whereas in Fig. 3b, c, both the NIPAM8 and NIPAM12 molecular systems exhibit two aggregation sites, and the curcumin molecules are situated at both sites and surrounded by the NIPAM molecules. In Fig. 3d, e, the NIPAM16 and NIPAM20 molecular systems form multiple clusters and are dispersed at different corners of the box.Fig. 3: Configuration evolution of the simulated systems.a Configurations of the CUR-NIPAM4 (a system consisting of curcumin molecules and 4-mer N-isopropylacrylamide molecules) system before and after running simulations for 60 ns time, b–e correspond to the situations of the CUR-NIPAM8 (a system consisting of curcumin molecules and 8-mer N-isopropylacrylamide molecules), CUR-NIPAM12 (a system consisting of curcumin molecules and 12-mer N-isopropylacrylamide molecules), CUR-NIPAM16 (a system consisting of curcumin molecules and 16-mer N-isopropylacrylamide molecules), and CUR-NIPAM20 (a system consisting of curcumin molecules and 20-mer N-isopropylacrylamide molecules) systems. For clarity, different geometries are used in the figure to show NIPAM (N-isopropylacrylamide), CUR (curcumin), and SOL (solvent).Analysis of the interactions between NIPAM and CUR moleculesIn the initial stages of a chemical reaction or interaction, molecules approach each other through electrostatic attraction. The molecular surface electrostatic potential (ESP) distribution can be utilized to predict or explain the relative orientation of molecules, binding strength, receptor-ligand binding modes, adsorption of molecules, etc. in complexes. Its basis is that molecules are easily in contact with each other in a complementary ESP manner, that is, to minimize the overall energy, the positive areas of the molecular surface ESP tend to contact the negative regions, and the larger the positive and negative values, the stronger this tendency46. The surface ESP of a molecule is calculated by the following formula47:$$V({\boldsymbol{r}})=\mathop{\sum }\limits_{A} \frac{{z}_{A}}{\left|{{\boldsymbol{R}}}_{A}-{\boldsymbol{r}}\right|}-\int \frac{\rho ({{\boldsymbol{r}}}^{{\prime} })}{\left|{{\boldsymbol{r}}}^{{\prime} }-{\boldsymbol{r}}\right|}d{{\boldsymbol{r}}}^{{\prime} }$$
(2)
where \({\bf{r}}\) is the positional coordinate of a point on the surface or outside of the molecule, \({Z}_{A}\) is the nuclear charge of atom A, \({{\bf{R}}}_{{\rm{A}}}\) is the positional coordinate of the nucleus, and \(\rho ({{\bf{r}}}^{{\prime} })\) is the electron density within the volume element \(d{{\bf{r}}}^{{\prime} }\). To facilitate the analysis, based on the calculated van der Waals surface ESP of the NIPAM4 and CUR molecules, an electron-density isosurface coloring method is employed to visualize the distribution of these values on the molecular surfaces, and the results are displayed in Fig. 4.Fig. 4: ESP (electrostatic potential) mapped the van der Waals surface of molecules studied.a NIPAM4 (4-mer N-isopropylacrylamide) molecule. b CUR (curcumin) molecule. Significant surface local minima and maxima of ESP are represented as red and blue. The unit is in a.u. (atomic unit).As shown in Fig. 4a, the side of oxygen atoms of the amide group (-CONH2) in NIPAM4 exhibits negative ESP (red), while the other side and almost all of the isopropyl (–CH(CH3)2) edges reveal positive ESP (blue). Regarding the surface of the CUR molecule shown in Fig. 4b, most of the regions of the two phenolic ring planes and their connected main chains display negative values, while the edges of the molecule display local positive values. Based on the above results, the interaction sites of NIPAM with CUR can be qualitatively classified into two cases, one of which is the binding of the oxygen atoms side of the amide groups of NIPAM to the side edges of the trailing ends of the two phenolic rings of CUR, and the other isopropyl groups of NIPAM to the frontal faces of the two phenolic rings of CUR and their connected backbone regions, and it is clear that the latter has a higher probability of being the predominant mode of binding due to the possession of more contact areas.According to the AIM theory, \({sign}({\lambda }_{2})\rho\) can be defined as a function to differentiate the type and strength of interactions in the system, where \(\rho\) is the actual electron density of the current system, \({\lambda }_{2}\) is the second largest eigenvalue of the Hessian matrix for electron density, and sign() stands for taking the sign of \({\lambda }_{2}\)48. The correspondence between the \({sign}({\lambda }_{2})\rho\) values and the interactions can be visualized in Fig. 5a, where the types and strengths of the interactions are differentiated by colors49. Furthermore, Lefebvre et al.50 define a function \({\delta }_{g}\) for the difference of the electron density gradient based on their proposed independent gradient model, which is denoted as follows:$${\delta }_{g({\boldsymbol{r}})}=\left|\sum _{i}{abs}\left[\nabla {\rho }_{i}({\boldsymbol{r}})\right]\right|-\left|\sum _{i}\nabla {\rho }_{i}({\boldsymbol{r}})\right|$$
(3)
where \(i\) is the atomic number, \(\nabla \rho\) is a vector of the electron density gradient, and \({abs}(\nabla \rho )\) means that absolute values are taken for each component of the \(\nabla \rho\) vector. The physical meaning of \({\delta }_{g}\) is that the stronger the interaction between the systems, the larger the value of \({\delta }_{g}\) in the interaction region will be.Fig. 5: The visual representation of the types and regions of intermolecular interactions using the part circled by the red dashed line in Fig. 4a as the studied object.a Correspondence between the values of \({sign}({\lambda }_{2})\rho\) and the types and strengths of the interactions. b Scatter plots of \({\delta }_{g}^{{inter}}\) vs \({sign}({\lambda }_{2})\rho\) (red scatter points) and \({\delta }_{g}^{{intra}}\) vs \({sign}({\lambda }_{2})\rho\) (black scatter points) based on the independent gradient model for the system consisting of NIPAM4 (4-mer N-isopropylacrylamide) molecules and CUR (curcumin) molecule defined as two fragments, respectively, where \({\delta }_{g}^{{inter}}\) and \({\delta }_{g}^{{intra}}\) represent the difference function of the electron density gradients between the fragments and that within the fragments, respectively. c The system composed of NIPAM4 molecules and the CUR molecule is defined as two fragments, and the isosurface map of the calculated \({\delta }_{g}^{{inter}}\) is plotted to show the types and regions of interactions between the fragments. d Each NIPAM4 molecule is defined as a fragment separately, and the corresponding isosurface map of \({\delta }_{g}^{{inter}}\) is plotted to visualize the types and regions of interactions between NIPAM4 molecules.In this paper, to reveal the essence of the intermolecular interactions between NIPAM and CUR, the portion circled by the red dashed line in the configuration obtained after simulation of the CUR-NIPAM4 system shown in Fig. 3a is extracted as the study region of interest. The region is defined as two fragments, one consisting of NIPAM4 molecules and the other being CUR. The \({\delta }_{g}\) in the interaction regions between the two fragments mentioned above can be expressed as \({\delta }_{g}^{{inter}}\), and they are calculated and used to plot into scatter plots of \({\delta }_{g}^{{inter}}\) vs \({sign}({\lambda }_{2})\rho\) as shown by the red scatter points in Fig. 5b. The \({\delta }_{g}\) of the interaction regions between the atoms within the fragments is denoted as \({\delta }_{g}^{{intra}}\), and they are plotted as the black scatter points in Fig. 5b. As shown in Fig. 5b, there is a peak consisting of red scatter points at the position where the value of \({sign}({\lambda }_{2})\rho\) is about −0.04. Since the electron density at this position is not very large and the peak position is close to 0, it can be inferred that there is a small amount of hydrogen bonds between the two fragments (the NIPAM4 molecules and the CUR molecule). There are a large number of black scatter points in the region where \({sign}({\lambda }_{2})\rho\) is significantly greater than 0, indicating the presence of site-blocking in the system, and there are large areas of scatter points at larger negative values of \({sign}({\lambda }_{2})\rho\), which can be considered as points corresponding to the chemical bonding region due to the high electron density at these positions. Plotting the isosurfaces of the calculated \({\delta }_{g}^{{inter}}\) and projecting the corresponding values of \({sign}({\lambda }_{2})\rho\) onto these isosurfaces in different colors (color scale shown in Fig. 5a) allows a clear examination of the regions of interactions existing between the fragments and the type and strength of the interactions, the results of which are shown in Fig. 5c. In Fig. 5c, the interaction regions are mainly green in color, indicating that the interaction between the NIPAM4 molecular system and the CUR molecule is dominated by van der Waals forces. The centers of some green areas show a conspicuous blue color, corresponding to hydrogen bonding, which is consistent with the inference of Fig. 5b. For details, these interaction regions are mainly spread along the planes of the two phenolic rings of CUR, and these interactions are more contributed by the atoms on the planes of the two phenolic rings and on the isopropyl group of NIPAM4, which is generally in accordance with the conclusion obtained from Fig. 4. Besides, defining each NIPAM4 molecule in the study region as a fragment separately and plotting the isosurface map by using the corresponding calculated \({\delta }_{g}^{{inter}}\) can demonstrate the interactions between NIPAM4 molecules, such results are displayed in Fig. 5d. The various NIPAM4 molecules in Fig. 5d are colored differently to facilitate the investigation of their interactions. The results show that van der Waals interactions and hydrogen bonding are the main driving forces for their aggregation into clusters.A concern that may arise is whether NIPAM can produce a similar kinetic behavior for other drug molecules that also have multiple ring structures as it does for CUR. For this reason, doxorubicin (DOX), and paclitaxel (PTX) are selected as the comparison drugs, and the DOX-NIPAM4 and PTX-NIPAM4 systems are constructed, which possess almost identical initial conditions as CUR-NIPAM4, and they are also subjected to run simulations of 60 ns time length. The architectures obtained at the completion of the simulations are shown in Supplementary Fig. 3. From the figure, it can be observed that similar to CUR-NIPAM4, the self-assembled NIPAM4 clusters in the DOX-NIPAM4 and PTX-NIPAM4 systems are capable of adsorbing (encapsulating) drug molecules. In detail, the cyclic structures of the three-drug molecules can all be inserted into the NIPAM4 system to varying degrees (or, the cyclic structures can attract NIPAM4 molecules), with the characteristic distance between the cyclic structure and the nearest atoms on its two sides being the same (measured to be 3–8 Å). This indicates that weak interactions serve as the cohesive force between them. Other differences regarding the relative positions and changes in force between NIPAM4 clusters and drug molecules can be attributed to structural and property differences in these drug molecules themselves, and their discussion is not the central goal of this paper. Supplementary Fig. 3 illustrates the universal applicability of NIPAM in DDS, at least in producing favorable drug loading patterns for various drugs with cyclic structures.The influence of chain length of NIPAM polymers on their efficacy in loading CURThe number of atomic contacts between molecules is an important parameter reflecting the degree of intimacy between molecules. If the distance between two atoms is defined as mutual contact when it is less than or equal to 0.6 nm, the number of atomic contacts between the systems51 can be counted using the following equation:$${N}_{C}\left(t\right)=\mathop{\sum }\limits_{i=1}^{{N}_{A}} \mathop{\sum }\limits_{j=1}^{{N}_{B}} \int ^{{r}{i}+0.6{nm}}_{{r}_{i}}\delta \left(r\left(t\right)-{r}_{j}(t)\right){dr}$$
(4)
where \({N}_{A}\) and \({N}_{B}\) represent the total number of atoms in system A and system B, respectively, and \({r}_{i}\) is the distance from the jth atom in B to the ith atom in A. Figure 6a presents the number of contacts between the CUR and NIPAM molecular systems. As shown in the figure, the number of contacts first rises rapidly after the beginning of the simulation and then fluctuates around a certain value during the last ~15 ns. This confirms the process of dispersed NIPAM molecules gradually aggregating and covering the CUR to form stable complexes as illustrated in Fig. 3. Examining the behavior of these contacts over the final time frame in which equilibrium has been reached it is easy to conclude that the NIPAM4 molecules are in closer contact with the CUR molecules and that this closeness significantly decreases as the chain length of the NIPAM molecules increases from NIPAM4 to NIPAM12. However, further increases in the chain length to NIPAM16 and NIPAM20 do not show a continued decrease in contact; instead, the values are comparable to those of NIPAM12.Fig. 6: Distribution of relative positions between CUR (curcumin) and NIPAM (N-isopropylacrylamide) molecules in the simulated systems.a Variation of the number of contacts between CUR and NIPAM molecules with simulation time in the simulated systems CUR-NIPAM4 (a system consisting of curcumin molecules and 4-mer N-isopropylacrylamide molecules), CUR-NIPAM8 (a system consisting of curcumin molecules and 8-mer N-isopropylacrylamide molecules), CUR-NIPAM12 (a system consisting of curcumin molecules and 12-mer N-isopropylacrylamide molecules), CUR-NIPAM16 (a system consisting of curcumin molecules and 16-mer N-isopropylacrylamide molecules), and CUR-NIPAM20 (a system consisting of curcumin molecules and 20-mer N-isopropylacrylamide molecules). b Calculated radial distribution functions of NIPAM molecules around the CUR in the CUR-NIPAM4, CUR-NIPAM8, CUR-NIPAM12, CUR-NIPAM16, and CUR-NIPAM20 simulated systems.The radial distribution function (RDF) expresses the probability of finding an atom/particle in the \({dr}\) shell layer at a distance r from an atom/particle as a reference point. Its formula (defined equation) can be expressed as follow45:$${g}_{A-B}\left(r\right)=\frac{\left\langle {\rho }_{B}(r)\right\rangle }{{\left\langle {\rho }_{B}\right\rangle }_{\rm{local}}}\frac{1}{{\left\langle {\rho }_{B}\right\rangle }_{\rm{local}}}\frac{1}{{N}_{A}}{\int_{i\in A}^{{N}_{A}}} {\int_{j\in B}^{{N}_{B}}}\frac{\delta ({r}_{{ij}}-r)}{4\pi {r}^{2}}$$
(5)
Where \(\langle {\rho }_{B}(r)\rangle\) is the partial density of component B at a distance r from A, and \({\langle {\rho }_{B}\rangle }_{\rm{local}}\) is the partial density of component B in all spheres of radius r around component A. To further characterize the distribution of NIPAM molecules around the CUR, the RDF of the studied molecular systems is calculated and shown in Fig. 6b. It is shown that the peak positions of the RDF curves of the five simulated systems are very similar, and the values of the peak positions exhibited imply that the main distance at which the interactions between the CUR and NIPAM molecules occur is about 0.53 nm. The main difference between these RDF curves lies in their intensities, among which CUR-NIPAM4 exhibits the highest one, indicating that CUR has the highest adsorption capacity for NIPAM4, and this adsorption capacity gradually decreases as the chain length of the NIPAM molecules increases from NIPAM4 to NIPAM12. This trend does not continue with further increases in the NIPAM chain length (NIPAM16 and NIPAM20); instead, the values are very similar to those of NIPAM12, and such a result well supports the inference obtained from Fig. 6a.Typically, hydrogen bond can be determined using a combination of a cutoff value for the hydrogen donor-acceptor angle (30°) and a cutoff value for the donor–acceptor distance (0.35 nm), where the hydrogen donors are OH and NH groups and the acceptors are by default O and N52. Figure 7 displays the variation of the number of hydrogen bonds with simulation time for the simulated systems with five different chain lengths of polymerized molecules studied.Fig. 7: The number of hydrogen bonds in the simulated systems.a Variation of the number of hydrogen bonds between NIPAM4 (4-mer N-isopropylacrylamide) and solvent molecules (NIPAM4-SOL), between NIPAM4 and NIPAM4 molecules (NIPAM4-NIPAM4), and between NIPAM4 and CUR (curcumin) molecules (NIPAM4-CUR) in the CUR-NIPAM4 (a system consisting of curcumin molecules and 4-mer N-isopropylacrylamide molecules) with simulation time. b–e correspond to the situations of CUR-NIPAM8 (a system consisting of curcumin molecules and 8-mer N-isopropylacrylamide molecules), CUR-NIPAM12 (a system consisting of curcumin molecules and 12-mer N-isopropylacrylamide molecules), CUR-NIPAM16 (a system consisting of curcumin molecules and 16-mer N-isopropylacrylamide molecules), and CUR-NIPAM20 (a system consisting of curcumin molecules and 20-mer N-isopropylacrylamide molecules) systems, respectively. f The number of hydrogen bonds between NIPAM and CUR molecules in CUR-NIPAM4, CUR-NIPAM8, CUR-NIPAM12, CUR-NIPAM16, and CUR-NIPAM20 systems is displayed in the same coordinate system for comparison.As shown in Fig. 7a, the number of hydrogen bonds between NIPAM4 and water molecules (NIPAM4-SOL) has a large value initially, which is caused by the fact that at this time, the NIPAM4 molecules are dispersed in the water in a stretched manner, which can interact with a significant number of water molecules. With the running of the simulation, the hydrogen bonds of NIPAM4-SOL decrease rapidly, while those of NIPAM4-NIPAM4 and NIPAM4-CUR increase, which is attributed to the fact that the randomly dispersed NIPAM4 molecules gradually gather into clusters and wrap around the CURs, so that the contacts between NIPAM4 and water molecules reduce, while that between NIPAM4 molecules and that between NIPAM4 and CUR molecules rise, and the increase in the number of contacts creates the conditions for the formation of more hydrogen bonds between these molecules. It should be pointed out that the numerical values and increasing trends of the NIPAM4-CUR curve are relatively small, and their rationality can be guaranteed by the fact that there are only four CUR molecules in the simulated system. The number of hydrogen bonds between different molecules in Fig. 7a all maintain small fluctuations for nearly 15 ns near the end of the simulation, which is a reflection of the simulated system reaching a steady state. The same curve trends apply to the NIPAM8 simulation system (Fig. 7b), the NIPAM12 simulation system (Fig. 7c), the NIPAM16 simulation system (Fig. 7d), and the NIPAM20 simulation system (Fig. 7e). Overall, Fig. 7a–e present one form of evidence for the evolution of NIPAM molecules from dispersion to aggregation and reveal the important role that hydrogen bonds play in such a process, which also confirms the reliability of the visualized hydrogen bonding interactions demonstrated in Fig. 5c, d. In addition, when the NIPAM-polymer cluster is used as a DDS, its interactions with the CURs may be of greatest interest, and for this reason, the number of hydrogen bonds of NIPAM-CUR for the five systems is displayed in the same coordinates for comparison in Fig. 7f. While the difference between the five curves is not very significant, it basically conforms to the numerical magnitude relationship of NIPAM4 > NIPAM8 > NIPAM12 > NIPAM16 ~ NIPAM20, which suggests that the increase in chain length is unfavorable for the formation of hydrogen bonding interactions between the polymer molecules and the CURs.Taking the previous data analysis into account, it can be concluded that when the chain length of the NIPAM molecules is increased to NIPAM16 and NIPAM20, the simulated systems do not differ sufficiently to allow conclusions to be drawn about the apparent regularity associated with chain length changes. Therefore, when discussing the change in the interaction energy between molecules over simulation time in Fig. 8, the data for NIPAM16 and NIPAM20 can be considered similar to NIPAM12 and are not presented in the figure in order to facilitate comparative analyses between the curves. As shown in Fig. 8a, at the initial state, the electrostatic interaction energies (Elec.) among NIPAM molecules satisfy the relationship NIPAM4 > NIPAM8 > NIPAM12, as does the van der Waals interaction energies (Vdw.). This numerical difference is mainly caused by the differences in the position and structure (chain length) of randomly inserted molecules in the initial system. Elec. are some large positive values indicating the presence of electrostatic repulsion between some molecules in the systems, and their contribution to the Elec. is higher than the electrostatic attraction between other molecules. The values of Elec. decrease sharply in the first 5 ns of simulation time and then flatten out, suggesting that the interspaces between some of the repulsive molecules in the systems expand rapidly at the beginning of the simulation and that this process is more pronounced than the reduction of intermolecular distances due to electrostatic attraction in the other part of the molecules. The negative values of Vdw. imply the presence of significant van der Waals forces of attraction in the systems. The Vdw. decrease gradually with simulation time until they remain almost constant in the last 15 ns, which corresponds to the process by which mutually attracted molecules in the systems progressively shorten their molecular spacing and finally stabilize. Overall, the synergistic combination of attractive and repulsive forces drives and maintains changes in the systems. Obviously, molecular aggregations guided by attractive forces (including van der Waals attraction and electrostatic attraction) are dominant, which can be found as intuitive evidence from the evolution of the overall conformation of the NIPAM molecules illustrated in Fig. 3.Fig. 8: Interaction energies in the simulated systems.a Elec. (Electrostatic interaction energy) and Vdw. (van der Waals interaction energy) between NIPAM (N-isopropylacrylamide) molecules (NIPAM4 (4-mer N-isopropylacrylamide)-NIPAM4, NIPAM8 (8-mer N-isopropylacrylamide)-NIPAM8, and NIPAM12 (12-mer N-isopropylacrylamide)-NIPAM12) corresponding to the CUR-NIPAM4 (a system consisting of curcumin molecules and NIPAM4 molecules), CUR-NIPAM8 (a system consisting of curcumin molecules and NIPAM8 molecules), and CUR-NIPAM12 (a system consisting of curcumin molecules and NIPAM12 molecules) systems, respectively. b The Elec. and Vdw. between NIPAM molecules and CURs (curcumin) (NIPAM4-CUR, NIPAM8-CUR, and NIPAM12-CUR) in the CUR-NIPAM4, CUR-NIPAM8, and CUR-NIPAM12 systems, respectively.Figure 8b compares the interaction energies between the NIPAM molecules with different chain lengths and the CUR molecules. Initially, all Elec. and Vdw. are close to zero, which is consistent with the fact that the four CUR molecules are dispersed in the systems and thus have very large spacings from the vast majority of NIPAM molecules. All Elec. and Vdw. reduce significantly with the prolongation of the simulation time, indicating that the electrostatic and van der Waals attractions are gradually narrowing the distances between the NIPAM molecules and the CURs, and the magnitude of the larger reduction in Vdw. suggests that the van der Waals forces are dominant in this process. By comparison, it can be seen that both Elec. and Vdw., NIPAM4-CUR show the relatively smallest values (i.e., the largest amount of change), indicating the strongest adhesion between NIPAM4s and CURs. Correspondingly, NIPAM12-CUR exhibits the relatively highest values (the smallest amount of change), corresponding to the weakest interaction. This quantitatively confirms the conclusion that an increase in chain length is detrimental to the formation of a stabilizing cover for CUR by NIPAM molecules. In the simulation system, CUR and NIPAM4 molecules continuously adjust their positions and orientations within the simulation box over time. Within a very short time frame, these adjustments are necessary for the molecular system to reach a relatively lower energy state (an instantaneous stable state). This is the reason why the curves in Figs. 7 and 8 exhibit noticeable fluctuations over short time intervals. As the above process continues, it can be seen in Supplementary Fig. 2 that the relative positions of the CUR and NIPAM4 molecules are dynamically changing with a general tendency to become closer (more molecular contacts), which corresponds to the overall curvilinear trend of an increased number of hydrogen bonds and larger interaction energies (absolute values) exhibited in Figs. 7 and 8.Temperature-sensitive effects on loading CUR by NIPAM molecular clustersFor exploring the influence of the temperature-sensitive effect of NIPAMs on their efficacy in loading CUR, the simulated data of the CUR-NIPAM4 system at temperatures of 300 K (conventional ambient temperature) and 310 K (human body’s characteristic temperature) are extracted for comparative analyses, as shown in Fig. 9. Since the differences exhibited by the systems compared to each other are of most interest after their simulations have been run to near stability, Fig. 9 shows only the data for the systems during the last 15 ns of simulation time.Fig. 9: Effect of temperature on the configurations and intermolecular interactions of simulated systems.a Demonstration of the difference in the intermolecular Vdw. (van der Waals interaction energy) of the CUR-NIPAM4 (a system consisting of curcumin molecules and 4-mer N-isopropylacrylamide molecules) system in the last 15 ns time range under simulated conditions at 300 K and 310 K. Correspondingly, the b corresponds to the Elec. (Electrostatic interaction energy), the c corresponds to the Rg (radius of gyration) of the NIPAM4 (4-mer N-isopropylacrylamide) molecular system, and the d corresponds to the variation in the number of hydrogen bonds between the solute molecular system (NIPAM4&CUR (curcumin)) and the SOL (solvent).Comparing the black and red curves in Fig. 9a, Vdw. increases after elevating the temperature, which indicates that the van der Waals interactions between the NIPAM4 molecular cluster and the CUR are weakened. The values of the black and red curves over time in Fig. 9b do not show a significant difference in magnitude, implying that the raised temperature does not have a noticeable effect on the electrostatic interactions between NIPAM4 and CUR. Taken together, the elevated temperature cut down the adhesion between NIPAM and CUR molecules. Such a conclusion will bring important implications as it points out that the NIPAM molecular cluster can release CUR by elevating the temperature, which is particularly important for researchers to design anti-tumor DDSs with temperature as the excitation factor based on the fact that the temperature of the tumor site is higher than that of normal tissues in vivo.Comparing the gray and pink curves in Fig. 9a, as the temperature increases, Vdw. decreases significantly, indicating that the distances between NIPAM4 molecules are further reduced and their van der Waals interactions are dramatically intensified. However, in Fig. 9b (gray and pink curves), the elevation of temperature makes the positive Elec. decrease, suggesting that the repulsive forces between NIPAM molecules are further released, which corresponds to the process of the expanding of the NIPAM4 molecular interspaces. It is clear that the process described in (a) is dominant since the amount of change in Vdw. is significantly higher than that in Elec. when the temperature is elevated. Such a conclusion can be supported by Fig. 9c, in which the higher-temperature NIPAM4 molecular system has significantly smaller Rg, implying that this system has higher compactness as a result of further contraction of the NIPAM4 molecules. In addition, solubility in water is an important reference for bioavailability, which can be reflected by the number of hydrogen bonds between the target substance and water molecules. In Fig. 9d, the number of hydrogen bonds between PNIPAM4-CUR and water molecules markedly decreases after the temperature increase, demonstrating that the elevation of temperature weakens the water solubility of the solute system.To further investigate the effect that elevated temperatures have on the relative positional distributions between different molecules, the RDF of the molecules of interest at different temperatures are calculated and shown in Fig. 10. Figure 10a shows the distribution of NIPAM4 molecules around CUR. In this figure, the shape and peak positions of the RDF curves at different temperatures are similar, with the main difference being the significant increase in peak intensity after heating. This can be considered as the further aggregation of NIPAM4 molecules into smaller molecular clusters in the system, which increases the molecular distribution density around CUR. Figure 10b, c shows the RDFs of water molecules around NIPAM4 and CUR molecules, respectively. The shapes and peak positions of the RDF curves in the figure are similar at different temperatures, but the peak intensities decrease after heating, indicating that increasing temperature leads to a decrease in the number of water molecules around NIPAM4 and CUR. This is consistent with the results of the decrease in solubility of the NIPAM4-CUR solute system shown in Fig. 9d.Fig. 10: Effect of temperature on the distribution of relative positions between molecules in the simulated system.RDF (radial distribution function) of NIPAM4 (4-mer N-isopropylacrylamide) molecules around CUR (curcumin) (a), RDF of SOL (solvent molecules) around NIPAM4 (b), and RDF of SOL around CUR (c) for the CUR-NIPAM4 (a system consisting of curcumin molecules and NIPAM4 molecules) system under simulated conditions at 300 K and 310 K.The temperature sensitivity of NIPAM can be attributed to its unique polymeric structure53. The structure of NIPAM contains two distinct parts, the hydrophobic isopropyl group and the hydrophilic amide group. When the temperature is low, the polymerization chain of poly-NIPAM is mainly in the foaming state, and the different groups on the chain are in a relatively stretched state due to the weak interactions. An increase in temperature will cause strong aggregation between the isopropyl groups of poly-NIPAM chain segments and even form hydrophobic microdomains. Such behavior of further coalescence between molecular groups and between molecules occupies the spaces for the water molecules in the clusters, causing them to be expelled. At the same time, the reduction in the volume of the molecular clusters (and thus the surface areas) can lead to a reduction in the probability of water molecules contacting the hydrophilic groups and the whole system. This is the reason for the results in Fig. 10b, c, which also explains why the water solubility of the NIPAM4-CUR molecular cluster in Fig. 9d decreases with increasing temperature. Going back to the previous discussion on Figs. 4 and 5, it is concluded that the position of the isopropyl groups in the NIPAM4 structure and the CUR phenolic ring plane are the main interaction regions. After the temperature increases, the condensation and contraction of isopropyl groups within and between molecules in the NIPAM4 molecular system change the layout of alkyl chains on these groups that are originally dispersed around the CUR. Based on the distribution characteristics of CURs embedded at the edge of the NIPAM4 molecular system rather than embedded in its center (see Fig. 3), the behavior of the alkyl chains on the isopropyl groups of the NIPAM4 molecules to coalesce into clusters throughout the inner part of the molecular system would reduce their contact with the CURs to some extent, i.e., the elevated temperature reduces the probability of contact between the effective interaction regions in NIPAM4-CUR. This is the reason for the weakened adhesion between NIPAM4 clusters and CURs. Overall, our simulation results not only corroborate existing experimental findings that NIPAM-related polymers can be stimulated to release drugs upon heating33,34,35,37, but more importantly, they elucidate the microscopic mechanisms underlying this process at the molecular level. These mechanistic insights, combined with previous findings, can provide effective guidance for the design of DDS. Specifically, based on the simulation results and analysis of systems with different NIPAM chain lengths, it is contemplated that when preparing DDS either from pure NIPAM or from NIPAM crosslinked with polymers having similar chain-like structures, molecular systems with shorter chains are more effective in encapsulating and entrapping CUR. Starting from the research conclusion that electrostatic and van der Waals interactions are the primary driving forces behind the self-assembly process of NIPAM and CUR, it is suggested that during the preparation of DDS, optimizing the solvent environment and ion strength can be considered to enhance these interactions, thereby improving drug loading or release efficiency. The simulation results of the CUR-NIPAM4 system at different temperatures reveal that, in order to achieve stable delivery and release of CUR drugs, it is possible to consider preparing NIPAM-based DDSs at room temperature or even lower temperatures, and stimulate drug release at the tumor or inflammation site with the help of its temperature higher than the body temperature (or the local temperature raised by the stimulation of the in vitro device). This requires the determination of the appropriate target temperature based on the specific DDS structure to ensure optimal drug release characteristics. From the analysis of molecular surface ESP distribution and electron density difference function \({\delta }_{g}\), the screening of the system of DDS drugs can be carried out by considering a sufficient number of isopropyl groups or its similar structures in the DDS matrix as well as a sufficient number of cyclic structures in the drug molecules, which are likely to result in stable interactions and sensitive temperature response to enhance the drug delivery performance.
Molecular understanding of the self-assembly of an N-isopropylacrylamide delivery system for the loading and temperature-dependent release of curcumin
