A new approach to thermal responsivity and hydrodynamics of poly(N-isopropylacrylamide) in aqueous solution: all-atom and mesoscopic simulations

In this study, the impact of chemical and physical characteristics of poly(N-isopropylacrylamide) (PNIPAm) in water on its thermoresponsive behavior and aggregation was investigated via all-atom and mesoscopic simulations. The calculated dimensional and spatial parameters of the PNIPAm proved the existence of a critical chain length consisting of 15 repeating units (1697 g mol−1) which beyond it the PNIPAm conformations depend on the temperature. The influence of chain length on the interfacial behaviors and non-bonded interactions of PNIPAm in water indicated that the thermoresponsive behavior of the PNIPAm chain is independent of the chain length, from interactional aspect. Interestingly, the interactional barriers against the PNIPAm temperature dependency was correlated to the polymer intramolecular interactions and their contribution to heat capacity (C p-intra ). The simulation results indicated the highest barrier energy against the polymer bending at the critical chain length with C p-intra = 1.0189 kcal mol−1. Additionally, the influence of the chain length on its aggregation mechanism in water was evaluated using DPD simulation to correlate the PNIPAm spatial shape to dynamics and frictional solvent forces acting on the solute. According to the simulation results, the PNIPAm has an innate ability for chain aggregation. The higher and less aggregation intensity was occurred for the chain with shorter and longer length than the critical chain length, respectively. The presented studies can be used to design PNIPAm block with a proper molecular weight for preparing an efficient copolymer with chain network formation ability in water.


Introduction
Thermoresponsive polymers are polymers that display a drastic and discontinuous change in their physical and chemical properties versus temperature [1]. These polymers belong to the stratification of stimuli-responsive materials whose physicochemical properties will be changed applying the external stimulants, such as temperature, pH, electric field, and ionic strength [2,3]. Among various mentioned external stimulants, temperature is the most attractive because of its accessibility and low cost [3]. In fact, thermoresponsive polymers almost change all of their behaviors with modifying the environmental conditions, especially in aqueous solution [4][5][6]. Further thought, these polymers also exhibit a miscibility gap in their temperaturecomposition phase diagram via raising the temperature [7,8]. Depending on whether the gap is detected at high or low temperature, an upper or lower critical solution temperature (UCST or LCST) occurs in their aqueous solutions, respectively [9]. Many researchers mainly focus on polymers that represent thermoresponsivity in water for different purposes, i.e. tissue engineering [10], liquid chromatography [11], drug delivery [12], bioseparation [13], and physical network structure [6]. The major research works on aqueous thermoresponsible polymer solutions belong to polymers that illustrate LCST at the predetermined temperature to change their solubility, due to their transition temperature range [12][13][14]. Cao et al [12] studied thermoresponsive peptide-PNIPAm hydrogels to control drug delivery. Their results showed that the PNIPAm molecules collapsed to form condensed globular particles when the temperature increased above its LCST. This issue causing the PNIPAm acts as cross-links to connect peptide nanofibrils together physically and freeze their movements. Bokias et al [13] determined the interaction between the PNIPAm chains and water-soluble polymer with octadecyl short branches using viscometry and turbidimetry methods. Their results showed stronger hydrophobic interactions between the PNIPAm and octadecyl branches above the PNIPAm's LCST leading to increase the viscosity of the polymer systems. The chemical structure of the thermoreversible polymers influence on its efficiency. PNIPAm is one of the most efficient thermoresponsive polymers to prepare thermoreversible aqueous systems [15]. The molecular architecture of PNIPAm consists of a hydrophobic carbon-carbon backbone and a hydrophilic amide group capped by a hydrophobic iso-propyl moiety [16]. This polymer has a coil conformation below its coil-to-globule temperature (T cg ) (miscible state) in water, whereas it collapses into a globule above the temperature (immiscible state) [17]. Researchers fulfilled many investigations to unmask the PNIPAm chain behaviors in aqueous solution from the molecular point of views [15][16][17][18][19][20]. De oliveira et al [15] measured the dimensional behaviors of the single-chain of PNIPAm and poly(N-n-propylacrylamide) (PNnPAm) at different temperatures. They could evaluate the effect of n-propyl and iso-propyl groups on the chain T cg variations. Walter et al [18] studied the conformation change of PNIPAm hydrogels in mixtures of water and methanol solvents. According to their results, the transition in the hydrogel conformation was dependent upon the solution composition (solvents mole fractions) and temperature. Chiessi et al [19] evaluated the influence of chain tacticity on the hydrophobicity of a single PNIPAm chain in water. Their results showed that the affinity of PNIPAm to water molecules decreased with increasing the degree of chain isotacticity. Du et al [20] investigated the effect of salinity on the LCST of PNIPAm in aqueous media. The presence of salt in the aqueous solution shifted its LCST to a lower temperature.
In this work, a comprehensive study on the molecular dynamic behaviors of PNIPAm with various chain lengths in water at different temperatures was investigated via molecular dynamics (MD) and dissipative particle dynamics (DPD) simulation. PNIPAm is a very attractive polymer for industrial and medical applications due to its thermal responsive behavior with a lower critical solution temperature (LCST). Despite many researches showing the PNIPAm behaviors in aqueous solutions, there are still some ambiguities related to changing the PNIPAm physicochemical properties in aqueous media with raising the temperature. It seems that there is an effective factor relating to intra-and intermolecular interactions between the PNIPAm and the solvent molecules, leading to different molecular behavior when the temperature increases. The predominant intra-and intermolecular interactions would influence thermal responsibility and aggregation behavior of the polymer chains in water. Moreover, most of the research works on PNIPAm are related to the chemical aspects of the polymer, while the physical characteristics of the PNIPAm chain can have a significant impact on its thermoresponsive and aggregation behavior. In this regard, investigation of the hydrodynamic behavior and frictional forces acting on the PNIPAm chains and their aggregate structure can clarify the relationship between the chain size and its dynamics, which have been neglected so far. Moreover, the thermoresponsive and aggregation behavior of the polymer chain with various chain lengths has not been studied in coarse-grained level. Generally, PNIPAm has been mainly used to construct a thermoresponsive copolymer with spatial arrangement capability. Therefore, the influence of physical and chemical characteristics of the PNIPAm on its thermal sensitivity and aggregation in an aqueous solution has great importance to tune an efficient thermoresponsive block copolymer. Hence, finding the origin of theses behavioral variations and mechanism of conformational changes of the chain in aqueous media will help molecular engineers to design thermoresponsive block copolymers for desired applications. Accordingly, the impact of physical and chemical characteristics of PNIPAm on its interactional and hydrodynamic behaviors in the media and the frictional forces acting on the PNIPAm chains and its aggregated structures are evacuated from the molecular point of view in atomic and coarse-grained levels. In fact, this study is performed to find the relationship between the physics and chemistry of the PNIPAm in water to reveal the interactional barriers in the system. The achieved results in this work are so useful for polymer scientists to choose a PNIPAm chain with the best molecular weight to tune an effective copolymer with desirable properties.

Molecular dynamics simulation details
For fundamental study on the chain behaviors of PNIPAm single molecule in aqueous media, all-atom (AA) models of PNIPAm with different chain lengths in water molecules were considered. In conjunction with PNIPAm, Walter et al [17] used various force fields and force field combinations to investigate the PNIPAm molecular behaviors in water. According to their results, the OPLS-AA (Optimized Potentials for Liquid Simulations-All Atom) force field [21] combined with the SPC/E water model [22] was the best combination to determine the PNIPAm chain parameters with the nearest values to experiment results. Accordingly, the OPLS-AA+SPC/E force field combination was used to evaluate the PNIPAm-water solution models at the atomic scale in this research work. The results of Mochizuki et al [23] research work on PNIPAm, eand s Lennard-Jones (LJ) parameters of the non-bonded interactions and also atoms electrostatic charge (q), were used for each atom of NIPAm monomer. The e, s and q values for each atom were separately listed in table S1 (available online at stacks.iop.org/MRX/8/105301/mmedia) (Supplementary Material) according to the atom numbering in figure 1. As seen, the applied parameters for NIPAm ingredients were according to the monomer position because, the inner and end monomers of polymer chain show different behaviors, particularly at oligomers [23]. The LJ cross-interactions parameters were calculated via mix geometric equations, e e e = ij ii jj and s s s = , ij ii jj in which ε and σ are the Lennard-Jones well depths and radii. The interactions of non-bonded dispersive were computed via cut-off distance of 12 Å, and the long-range Coulombic interactions were calculated using Particle Mesh Ewald [24]. Herein, MD simulations for all the PNIPAm aqueous solutions were performed using LAMMPS software package, developed at the Sandia National Laboratories [25]. At the simulation procedure, the temperature and the pressure were controlled via Nose-Hoover thermostat and Berendsen barostat, respectively. The coupling time of 500 fs was selected for both temperature (t T ) and pressure (t P ). In the simulations, the time step of 1 fs was chosen, and the Newton's equations of motion were integrated using the velocity Verlet algorithm. At the end, the configurations of the last 20 ns were used for data sampling.

Model construction
First, the structure of a NIPAm monomer was built and then was geometrically optimized via GAMESS software package with internally stored B3LYP/6-31G (d,p) basis set [26]. Next, the achieved monomer was repeated to generate the PNIPAm chains with three different lengths. The polymer chains were symbolized by pN, in which N is the number of chain repeating units, i.e., p8, p15, and p32 for the chains with 8 (905 g mol −1 ), 15 (1697 g mol −1 ), and 32 (3621 g mol −1 ) repeating units, respectively. Despite the fact that the PNIPAm tacticity has the ability to influence the physicochemical behaviors of the polymer [19], the process was controlled in order to obtain an atactic stereochemistry structure. Since, this type of tacticity is more common at the most of polymerization process. To build the PNIPAm solutions with the mentioned different chain length, the initial pN/water (pN/W) models containing only a single chain were constructed at a dilute regime in 3D structure (periodic boundary condition) to implement the MD simulation. The built 3D boxes' dimensions are large enough to prevent the intermolecular interactions of pN-pN. The dimensions of each initial pN/W dilute solutions boxes and the number of water molecules in the solutions are represented in table 1.
Afterward, the initial pN/W models were energetically optimized to reach a minimum level of energy via using the conjugate gradient (CG) algorithm. Then, after a pre-equilibration process, the simulations boxes were exposed to the NPT ensemble as the production run with the simulation time duration of 50 ns. It seems that the performance of the simulations with total simulation time duration of t tot = 50 ns is adequate, because it is much larger than the structural relaxation time of the PNIPAm chain [15,27], even for the longest chain with N=32. All the solutions were carried out at the constant atmospheric pressure, 1 atm, and the temperatures of 280, 290, 300, 310, 320, and 340 K.

Results and discussion
3.1. Role of chain length on PNIPAm thermoresponsivity 3.1.1. Chain dimensions and conformation According to experimental results, there is a coil-to-globule transition temperature (T cg ) for the PNIPAm/water solution close to 305 K [13,17,28]. In this phenomenon, the chain morphology and conformation will be changed significantly with increasing the temperature above T cg [20]. In this work, to investigate the effect of PNIPAm chain length on its coil-to-globule behavior, the dimension variation of the chains with different lengths (N=8, 15 and 32) in water were measured versus temperature. To achieve this goal, the pN radius of gyration (R g ) based on a weight average of the chain atoms distances from their center of mass in the equilibrated solutions was calculated as follows [29]: Where n, m and r represent the total number of the chain atoms, the atom mass and its distance from the chain center of mass, respectively. The subscript i is the chain atom number. In the beginning, the longest polymer chain (p32) in water was simulated at the highest temperature for the extra 50 ns (NPT simulation) to assure the efficiency of the used force field for prediction of the chain dimensional changes. Moreover, the consideration was performed to evaluate the sufficiency of the applied simulation time, as well. The evolution of R g for the p32 chain in water at 340 K, was calculated versus time for the main 50 ns plus extra 50 ns (see figure S1(a)). The chain conformations at the beginning and end of the first 50 ns simulation were also presented in figure S1(a). As seen in the figure, the applied simulation method and the force field could predict the significant variation of the chain R g at the highest temperature (T>T cg ) [13]. Accordingly, the force field combination has an ability to predict the molecular behaviors of PNIPAm [15,16]. For further confident of adequate simulation time, the R g evolution of the p32 versus simulation time at the temperature below T cg , i.e. T=280 K, is shown in figure S1(b). In this temperature region (T<T cg ), the dynamic properties of PNIPAm is blunt because water acts as a good solvent and interacts to the polymer strongly [7]. According to figures S1(a) and (b), the chain dimensions are almost constant after 50 ns production run (main simulation). Hence, the simulation time, 50 ns, is adequate for all the pN/W solution boxes. Oliveira et al [15] also used the simulation time duration of 50 ns for a PNIPAm single chain in water with the same chain length, N=32, and force field combinations. The calculated R g for all the pN chains in the equilibrated solutions at different temperatures were shown in figure 2(a). According to the figure, the dimensions of the p8 and p15 chains did not change with increasing the temperature. In fact, they represented an unexpected behavior with increasing the temperature. The observed behavior can be connected to the chain dimensional and conformational variations. Because, according to the obtained values of R g for the oligomers of p8 and p15, the chains' length was so small that illustrates the coil-toglobule transition behavior with reducing their dimensions. On the other hand, their dimensions were much lower than the chain collapsed state, which is obvious for the longest chain (figure 2(a)). Increasing the pN chain length to N=32 causes the chain achieves folding ability. This issue led to the chain could show the coil-toglobule change with changing its conformation. Also, the end-to-end distance (〈r 0 〉) values of the pN chains with different lengths indicated the aforementioned issuance (see figure 2(b)). As shown, the 〈r 0 〉amounts of the p8 and p15 chains were almost constant without any significant changes, while the 〈r 0 〉value for the p32 chain changed significantly with increasing the temperature.
To determine the coil-to-globule transition temperature, T cg , of the PNIPAm in water, Boltzmann sigmoid function was used only for the p32 chain to fit the R g values (equation (2)). For the p8 and p15 chains, the temperature independency due to their conformational and dimensional change limitations was observed. The proposed Boltzmann's sigmoidal function is based on the logistic sigmoidal equation of = + y e 1 1 .
This function is a relevant starting point for the prediction of transition behavior. The following Boltzmann function involving the required geometric characteristics was used to determine R T g ( )as follows [30]: Where T, R g expanded and R g collapsed are the temperature, the chain radius of gyration below and above the transition temperature, respectively. Factually, R g expanded and R g collapsed indicate the equilibrium amounts of the dependent variable before and after the transition point, respectively. The term of DT is the temperature width that describes the behavior of the slope of the process during the transition. The DT was obtained close to 6 K for the p32 chain which is in a good respect to the value reported elsewhere [15]. To improve the accuracy of the fitting process, the p32/W solution was also simulated at the temperatures of 305 and 315 K, and then their R g values were computed at the equilibrated state. The additional calculations were carried out for a proper curve fitting of equation (2) over the simulation results. Figure 3 indicates the calculated R g amounts and the fitted values versus temperature for the equilibrated p32/W dilute solution. The values are accurately fitted using the Boltzmann sigmoid function with its R-square (R 2 ) value of 0.9904. In fact, the Boltzmann sigmoid function as a mathematical function was used to characterize a sigmoid or S-shaped curve, similar to R g -T curve ( figure 3). Consequently, the T cg of 302.6 K was obtained for the p32 chain in water, which is in good agreement with the experimental result, i.e. 305 K [13,17,28,31] and the calculated value published previously, i.e. 302 K [15]. In fact, the achieved T cg illustrates that the applied simulation method and force field combination in this work was more accurate than the methodology used by other MD researchers.

Chain solvent accessible surface
The solvent accessible surface (SAS) or solvent accessible surface area (SASA) is a surface area of a polymer chain which is accessible for solvent molecules [32]. In this work, SAS of the pN chains was calculated to evaluate the temperature and the chain length effects on the chain spatial dimensions and its touchable surface for the water molecules. Figure S2 exhibits the temperature dependence of SAS for each pN chain in water. As shown, the SAS amount of the p32 chain decreases as the temperature increases, while its value was independent of the temperature for the shorter chains, i.e., p8 and p15. Factually, the thermoresponsibility of the PNIPAm chain depends on its chain length which originated from its conformational changes with increasing the temperature, from the dimensional aspect. Increasing the temperature led to increasing the chain bonded interactions and reducing the intermolecular interactions between p32 and the water media, which consequently lowered the SAS. Subsequently, p32 tended to diminish its surface area with media and get a globule conformation above T cg . This observation has conformity with the chain conformational and dimensional changes mentioned earlier.

Chain orientation
In order to consider whether the pN chains are folded in water or not, it is needed to find the spatial orientation of each chain in solution after equilibration. This study was accomplished via calculating the order parameter (S) of carbon-carbon bonds of the chain backbone in the equilibrated pN/W solutions using the following equation [33]: Where θ is the angle of the chain backbone bond with the reference direction of x, y, or z coordinate axis. The angular bracket, 〈K〉, means an average over the chain backbone C-C bonds at the data sampling section. For better understanding of the employed theory, the calculation method was schematically presented in figure S3. Depending on the chain backbone orientation, the order parameter can take an amount in the range of −0.5−1.
According to the equation, the values of −0.5 and 1 illustrate a perpendicular and parallel arrangement of the chain backbone relative to the reference direction, respectively. In fact, these two amounts indicate that the chain nearly has a rod conformation and does not have the ability for folding. Moreover, the zero value for S represents a random C−C bond orientation for the chain backbone, and it means that the chain has the ability to fold and collapse. Figure 4 exhibits the calculated S parameter for the pN chains with various lengths at different temperatures. As illustrated, the S values of the p8 and p15 chains were far away from zero at the whole temperature range. Subsequently, these shorter chains had a rod conformation at various T, and the temperature did not affect their orientations. In contrast, the p32 chain could change its conformation to random orientation of its backbone bonds via increasing the temperature above the T cg , due to approaching the S value to almost zero. As seen, the thermoresponsive behavior of the PNIPAm depends on its chain length, from the orientational point of view. Hence, the smaller PNIPAm chains, i.e., p8 and p15, did not show any significant change in their conformational feature with temperature variation. According to the results, the chain with longer length, N>15, can change its conformation and orientation with increasing the temperature. 3.2. Role of chain interactional characteristic on PNIPAm thermoresponsivity 3.2.1. Intermolecular interactions Properties of a polymer aqueous solution depend extremely upon the intermolecular interactions between polymer solute and water molecules [34]. The intermolecular or non-bonded interactions among the components of the solutions can be affected by changing the ambient circumstances, which vary polymer chain behaviors in the solution [35]. Therefore, investigation of the intermolecular interactions and also their contributions amounts at various conditions could reveal important insights into the interactional behaviors of the polymer chain in an aqueous media. In this work, the difference in electrostatic interaction energies (E elc ) at a temperature (T) and the temperature of 280 K, DE T , elc ( ) was calculated as follows: ) are the electrostatic contribution to the non-bonded interactions between the polymer chain and the water molecules at any temperature (T) and the temperature of 280 K, respectively. Similarly, the difference in the Lennard-Jones potential (E LJ ), DE T , LJ ( ) was also calculated according to the following equation.
)are the Lennard-Jones potential contribution to the non-bonded interactions of the pN-W pairs at any temperature (T) and at temperature of 280 K.
These calculations were performed to evaluate the effect of the chain length on the intermolecular interactional characteristics of the PNIPAm in water when the temperature increases. First, the process was carried out via estimating the contributions to the non-bonded interactions, electrostatic (E elc ) and Lennard-Jones (E LJ ) between pN and water in the equilibrated pN/W models as follows [36]: Where r ij is the distance between the element i and j. The calculated E elc and E LJ contributions of the non-bonded interactions of the water molecules with the pN chains and their amide groups per monomer are shown in figure  S4. The DE T elc ( ) and DE T LJ ( )amounts of pN-W and amide groups-W per chain monomer are shown for the p8/W, p15/W and p32/W equilibrated solution models (see figure 5). As seen, increasing the temperature resulted in the D E elc values with higher positive amounts. Actually, a decrease in compatibility between the polymer chain and water was observed with increasing the temperature. This issue occurred for all the models, even for the smaller chains, p8 and p15, without any significant changes in their conformational, dimensional and orientational features (see section 3.1). In fact, all the smaller and larger chains represented a temperature dependency behavior in the solution models, from the interactional aspect because of their intermolecular interactional changes. As shown in figure 5, the intermolecular interactional changes for the p8-water solution model were comparable to other solution models at the coil to globule transition temperature (T cg ). This observation can be related to more solubility of the smaller chains in the media ( Figure S4) while their SAS is constant with increasing the temperature ( Figure S2). On the contrary, the SAS of the longest p32 chain in water was decreased because its conformation changed from coil to collapse state with increasing the temperature. Thus, the thermoresponsive behavior of the p8 chain in water is remarkable in comparison with the other chains, p15 and p32. As shown in figure 5, the lowest change in DE elc values was observed for the p15/W solution model. This is due to lower solubility and less conformational change of the p15 chain in water when compared with smaller p8 and longer p32 chains, respectively. Therefore, there is an interactional obstacle in the p15 chain to reduce its energy level and reach a stable state. In this regard, the consideration of the chains' intramolecular potentials is necessary. In fact, the chain dimensions, chain contact area with the media, chain intra-potential and intermolecular interactions with the media are in competition with one another to control the energy of the solution models at any temperature. Moreover, the amide groups of the pN chains play a crucial role to change DE elc value with the temperature raising, because the amounts of their contribution were almost more than 50% of the DE elc . As shown in figures 5(b) and (d), the DE LJ contribution to the non-bonded interactions was much lower than the DE elc , because of the chains' nature. Hence, the domination intermolecular interactions of the pN chains and their amide groups with the water molecules are the electrostatic contribution, which play the key role to change the pN chain behaviors. Therefore, the electrostatic interactions between the water molecules and the amide moiety govern the pN chains solubility in the water media. The interactional considerations showed that the sensitivity of the PNIPAm chain to temperature at the lowest chain length depends on its interactional characteristics. Figures 2 and 5 clearly exhibit that the PNIPAm has chain length limitation to show the T cg , while its smaller bar-shape oligomers could illustrate partial instability in water with changing the temperature. Their thermoresponsive behavior is originated from the changing the interactional points of view.

Intramolecular interactions of the chains
Despite some studies on PNIPAm in water, less attention has been paid to the impact of the chain length and its solution temperature on the chain intramolecular interactions. This chain characteristic would change the PNIPAm molecular behaviors in aqueous solutions with temperature [37]. Therefore, study of the PNIPAm chain interactions in water can be helpful to disclose why the chain behavior changes at different temperatures, e.g., below and above T cg , and chain lengths. In this study, the difference in total intramolecular interactions of the pN chains per monomer in water (DE intra (T)) at any temperature (E intra (T)) and the temperature of 280 K (E intra (280 K)) were calculated at the equilibrated state (see figure 6). The calculated E intra (T) values versus temperature for the p8, p15 and p32 chains were also presented in figure S5. As shown in figure 6, the DE intra increased monotonously to more positive amounts for all the chains with increasing the temperature.  This increment can be attributed to the chains' tendency to a globule conformation to decrease its contact area with the solvent molecules when the temperature increases capacity [38]. The fitted lines and their R-squared (R 2 ) values are shown in figure 6. In fact, -C p intra is an amount of the intramolecular energy required per a PNIPAm repeating unit to increase the temperature one Kelvin, which is a crucial physicochemical property of the PNIPAm. The -C p intra values of 0.8806, 1.0189 and 0.7884 kcal/mol.K was obtained for the p8, p15 and p32 chain, respectively. According to the simulation results, increasing the chain length deceases the amount of the -C p intra due to increasing the chain folding ability, except for the p15 chain. In fact, the p15 chain has a critical chain length because it has neither sufficient length for chain folding and energy minimization with increasing the temperature nor shorter length for better solubility in water. Consequently, p15 does not move its segments freely. Therefore, the highest -C p intra values were obtained for the p15 when compared to the larger chain, p32, and the smaller one, p8 (figure 6). Although the intramolecular interactions of the PNIPAm chain depend on the temperature, there is a chain length limit to change its conformation.

Level of water capturing around the amide groups of PNIPAm
In solution, the capturing level of solvents in the hydration shell is obtained via calculating the total number of adjacent solvent molecules, solvent coordination number, around the central solute molecule [15,39]. Reality, this criterion can be used to investigate the structure of solvent molecules arranged close to the solute. In the polymer aqueous solutions, the structural variations in water molecules within the hydration shell play an important role on the spatial conformations of the polymer chain in media. In this study, the change of water coordination number around the pN amide groups (N A-W ) in the equilibrated solution was considered to investigate the effect of chain length and interactional characteristics on PNIPAm thermoresponsive behavior. Hence, the N A-W of the pN chains was calculated per the chain monomer to investigate the hydration level using the following equation [15]: Where l 0 is the range of the first solvation shell or the integration cutoff, which is chosen equal to 3.5 Å [16]. Also, the g(r) A-W is the radial distribution function (RDF) of the center of mass (CM) of the pN amide group with the water CM. The resultant N A-W for the chains at different temperature was shown in figure 7. As seen, the N A-W decreased significantly with increasing the temperature above the transition temperature of the PNIPAm. In fact, the issue illustrates that the coil-to-globule transition happened due to restructuring of the hydration shell caused by the changes of the non-bonded interactions. This behavior occurred for all the pN chains, especially more severely for the p8 chain due to its lower length and constant SAS. Also, more intermolecular interactions were observed for the p8 chain in water when compared to other longer chains. These observations approve the thermoresponsive behavior of the PNIPAm which depends on the chain characteristics. In addition, the thermoresponsive behavior of the PNIPAm with the smallest length in water is more severe, leading to intense aggregation of chains beyond T=T cg .

Mesoscopic study on aggregation behavior of PNIPAm
3.3.1. DPD simulation 3.3.1.1. Theoretical background Dissipative particle dynamics (DPD) is an efficient simulation method to study polymer multi-chain aggregation phenomenon in polymer solutions in the mesoscale level [40]. In this methodology, a group of atoms are coarsegrained into a particle as a DPD bead [41]. Each bead indicates the center of mass of a cluster of atoms, which occupy continuous position and velocity of r i and v i , respectively. The beads in DPD simulation obey Newton's equations of motion (equation (8)), which are used to consider the time (t) evolution of the beads [42].
Where i is the bead number, and m illustrates the bead mass. In the simulation, the non-bonded force acting on the i th bead is presented as f i having given sum over all the beads j that lie within a fixed cut-off distance [42].
Where the first, second and third term in ∑ illustrates a conservative, a dissipative and a random force, respectively. Conservative contribution of f i gives a chemical identity, while dissipative and random forces together create a thermostat that can keep the system at a constant temperature. Also, F i B and F i A exhibit the conservative forces because of the bonded interactions (B) and angles (A) between the bonded beads, respectively. F , ij C F ij D and F ij R pairwise additive forces are given as follows [43]: ii jj B / c ij and ρ are Flory−Huggins parameter and the number density, respectively [43]. g and s D are the friction coefficient and noise amplitude, respectively (s ). In equations (11) and (12), w D (r) and w R (r) are dimensionless weight functions for the dissipative and random forces ( , which is represented as follows: z ij is a delta-correlated Gaussian random variable, which its details were presented in Groot et al [42].
Additionally, F i B and F i A are defined by harmonic potentials for polymeric materials as follows [44]: Where, C B and C A are the modulus of bond and angel potentials, respectively. r ij 0 and q 0 as well as indicate an equilibrium length and angel, respectively.
The evolution of a DPD system is followed by integration of the equations of motion of each bead using the algorithm of velocity-Verlet as follows [45]:

Simulation details
In the current work, three solution systems of p8/W, p15/W and p32/W, with the same concentration of 3.006 g lit −1 were considered using DPD simulation to study the PNIPAm chain length effect on its aggregation intensity and hydrodynamic behaviors in water when the temperature increases. The simulations were carried out at the temperature of 300 and 310 K, the temperatures at below and above the polymer transition temperature, respectively. According to figure S6, one PNIPAm repeating unit and three molecules of water were selected to build their DPD beads. The former bead is shown as P and the latter bead is designated as W. MD simulation was used to calculate the c ij for PW pair, c , PW and the repulsion parameter between the beads, i.e., a PW . The all-atom simulations were performed on the atomistic models of the P and W beads with a box volume of 125 nm 3 . The amount of c PW at the temperature of 300 and 310 K was calculated based on the last 20 ns of atomistic simulations run (the simulation methodology presented in section 2.1) as follows [7,46]: Where E T , coh ( ) V m and AE i are the cohesive energy at temperature T, molar volume, and volume fraction of the component i, respectively. Consequently, the amounts of 26.60 and 29.95 were obtained for the a PW parameter at the temperature of 300 and 310 K, respectively. These results are in a good agreement with calculated values published previously [47].
To implement DPD simulation on the DPD models of pN/W (DpN/W), a set of required parameters, in DPD and real units [48], were listed in table 2. Moreover, the parameters of conservative forces responsible for the bonded interactions, PP bond and angle of PPP, presented in equations (14) and (15), were applied according to Li et al [44] research work. The boxes' volume is 8×10 3 nm 3 and all the DPD simulations were performed for 2×10 6 DPD steps.

Aggregation behavior
The equilibrated solution models of DpN/W at the end of DPD simulations were shown in figure 8. As seen, the DpN chains exhibit miscibility in water at 300 K<T cg . Increasing the temperature to 310 K led to the chains aggregation and reduction of the contact area between the chains and media. It is necessary to mention that all the aggregated structure of DpN chains (ADpN) were stable at the end of the simulation procedure. This behavior is due to the decrement of chains' tendency to water (see figure 5) [49]. For further investigation, the relative density (r r ) of the DpN chains in water beads was calculated in three x, y, and z directions and at the temperature of 300 and 310 K (see figure 9). According to figure 9, the r r of the DpN chains shows small peaks in three directions at 300 K, which means the chains' miscibility in water media. However, the r r -distance curve illustrates a strong and single peak in three directions of x, y and z indicating the spherical aggregation of DpN chains at the temperature of 310 K. This behavior was observed for all the DpN chains with different lengths in ii ij ji water beads. As seen in figure 9, the uniformity of ADpN structures was the highest for Dp32, Dp8 and Dp15 chains, due to the chain folding ability, small length, and having interactional barriers (see section 3.2.2), respectively. The results have a good conformity with the dimensional, conformational, orientational and interactional considerations of PNIPAm at different chain length, in the atomistic level presented in sections 3.1 and 3.2. In fact, the immiscibility and followed by aggregation of the chains in water, at the temperature above T cg , depend on the chain molecular characteristics. Additionally, to study the aggregation intensity of the DpN chains in the DpN/W models, inter-and intramolecular radial distribution function of the chains, i.e. g inter (r) and g intra (r), respectively, were calculated at the temperature above T cg , 310 K (figure 10). As seen in figure 10(a), the highest and lowest amounts of the g inter (r), implying the level of chain aggregation intensity, belonged to the DpN chains with N=8 and 32, respectively. These obtained results are in respect to the full-atom intermolecular interactions calculated for the pN-water models (see figure S4). Accordingly, the highest change in intermolecular interactions between the chain and water with increasing the temperature, was achieved for the smallest chain, p8. This remarkable change led to intensive thermoresponsive behavior and followed by more chain aggregation of the Dp8 chain when compared to the larger chains, Dp15 and Dp32. Moreover, less steric hindrance of the smallest chain compared to the longer chain also influenced the chain aggregation.
According to the g intra (r) results for the DpN chains ( figure 10(b)), increasing the chain length enhances further the chain folding ability and facilitates the coil-to-globule transition. As shown, the highest peaks on the g intra -r curves were observed for the Dp32 chains while the lowest one belongs to the Dp8 chains. Thus, the Dp32 chains considerably collapsed in comparison with the smaller chains. This observed behavior is in good agreement with the atomistic results (see figure 2).

Hydrodynamic behavior of immersed and aggregated chains
Hydrodynamic measurements are very important to characterize the dynamic behaviors of macromolecules in solution [50]. This consideration can provide deep understanding into the relationship between the dimensions and dynamics of a polymer, and the forces acting on the molecule by media [51] . This issue has crucial importance in micellization, aggregation and many other assembly phenomena in an aqueous media [52]. Because, the stability, dynamics and interactions of the polymeric structure with the media control the solution properties, particularly in drug delivery [53]. In this work, the hydrodynamic behaviors of the DpN chains and its aggregated structure, ADpN, before and after the aggregation, respectively were investigated. This evaluation was implemented to create a deep insight into the impact of chain length on the physical behaviors of DpN chains, below and above the polymer transition temperature in the water media.
In the beginning, the dimensions of DpN chains and ADpN structures were estimated via calculation of their gyration and hydrodynamics radii as follows [54,55]: Where k is the total number of beads, and the subscript CM represents the center of mass of DpN chain (in 300 K) and ADpN structure (in 310 K). r ij also indicates the distance between the beads i and j, and 〈K〉 also illustrates the ensemble average. The R h is the hydrodynamic radius of the solute (based on the Kirkwood definition [56]), the corresponding radius of a hard sphere which diffuses with the same rate as that solute. The calculated dimensional parameters were listed in table 3, for DpN chains below T cg , and ADpN structures above T cg . As seen in the table, the R g and R h of DpN chains were increased with raising the chain length at 300 K, and this is obvious [57]. Indeed, the important dimensional parameter is the ratio of R h /R g , which exhibits the spatial shape and conformation of the chains and an assembly structure in solvent media. According to the Mansfield et al [58,59] research works, the conformation of a polymeric structure in solution, with the shape of a smooth sphere, a random walk and an infinite long rod, gives the R h /R g value of 1.29, 0.79 and 0, respectively [60]. The calculated dimensional ratios for the DpN chains at 300 K show a random coil conformation for the chains. Nonetheless, Dp8 and particularly Dp15 tended to have a less compact conformation compared to Dp32, due to their lower amount of R h /R g [58]. In addition, the obtained dimensional parameters for ADpN structures at T>T cg indicate a decrease in R g and R h values with increasing N, due to enhancing the chain length and its ability for folding, particularly Dp32. The achieved amounts of R h /R g for aggregated structures, illustrate clearly that the ADpN had a spherical structure because of R h /R g @ 1.29 [58][59][60], and this issue was more severe for ADp32, due to the chain collapse. The achievements have a good conformity with the atomistic considerations of PNIPAm at different lengths, currently.
To study the dynamics behavior of the chains and its aggregated structure, mean square displacement (MSD) and velocity (V ) of the solutes were calculated as follows [61,62]: Since, the number of steps devoted for the diffusion of DpN chains and ADpNs is sufficiently long, subsequently, the diffusion regime will be Fickian [62]. The step evolution of MSD and V were shown in figures 11(a) and (b), respectively, for the solutes at 300 and 310 K. The obtained D amounts according to slopes of the best linear fits on the MSD-St curves were listed in table 3, with the fitting R 2 values. As seen in figure 11 and table 3, the chain mobility was decreased with increasing the chain length at T=300 K, due to enhancing the chain friction with the media because of the increment of chain size [8]. In fact, the frictional force acting on the DpN chains reduced their movement, due to the existing friction between the solute and the media. In addition, the behavior led to lowering the diffusivity and velocity of ADp8 and ADp15 compared to ADp32, as well. The schematic of this phenomenon was shown in figure S7 for a ADp32 structure. To quantitate this property, friction coefficient between the solute and the water media ( f ) and as well as frictional force or drag force (F D ) acting on the solute by the media were calculated according to the Stockes law as follows [8,63]: The calculated frictional parameters were presented in table 3, as well. As shown in the table, the increment of friction between the media and DpN chains at 300 K caused the enhancement of the frictional force and followed by the blunt dynamics of the solute. In 310 K, the frictional force acting on ADp15 was the highest compared to ADp32, even more than ADp8 despite the nearly same size. The obtained dimensional ratios clearly indicate that ADp15 had lower compression than other ADpN structures, due to the least amount of R h /R g . Moreover, a little deviation of ADp15 structure from the sphere, compared to the other aggregated structures, lowered the solute mobility in the water media (see figure 9). As seen, ADp32 had the lowest frictional interaction with the media, due to its pellet-like structure.

Conclusions
In this research work, MD and DPD simulations were implemented to study the effect of PNIPAm chemical and physical characteristics on its thermoresponsive, hydrodynamic and aggregation behaviors in water at the temperature below and above T cg . Atomistic simulations were performed at a dilute concentration regime and atmospheric pressure. First, the influence of PNIPAm chain length on its coil-to-globule transition behavior was considered. This issue was accomplished via calculating R g and 〈r 0 〉of the chains and followed by estimating the PNIPAm coil-to-globule transition temperature (T cg ). The simulation results showed that the p8 and p15 chains tended to have a rod conformation in water without any changes in their dimensions with increasing the temperature. The transition temperature of 302.6 K was calculated via fitting the Boltzmann sigmoid function on the achieved radius of gyration of the p32 chain, which had a good agreement with the experimental and calculated results published elsewhere. In addition, the computed solvent accessible surface area of the chain, the contact area to interact with the water molecules, was also supported by the aforementioned results. The obtained order parameters for the PNIPAm chains were confirmed the conformational and dimensional consideration. To investigate the thermoresponsivity intensity of PNIPAm, the non-bonded interactions contributions and differences of the interactions at temperature T compared to 280 K were determined. Despite the consideration of the chains' intramolecular interactions and their difference at each T compared to 280 K, the contribution of the intramolecular interactions to the thermal capacity was determined based on the total intramolecular interactions values. The interactional results indicate that the strength of intermolecular interactions between the PNIPAm chain and water decreases with increasing the temperature, particularly at the shortest chain. This change caused that the chains' intramolecular potential increased with increasing the temperature. The interactional variations occurred for all the chains with N=8, 15 and 32, although there was a critical chain length for the coil-to-globule transition, i.e. N=15. In addition, a decrement in the N A-W amount with the temperature was observed for all the pN chains in the hydration shell. Finally, the PNIPAm chain aggregation and hydrodynamic behaviors in water with increasing the temperature was evaluated via DPD simulation. According to the resultant data, the PNIPAm chains with different lengths can aggregate together, while the aggregated structures by Dp15 chains had the slowest dynamic and the most friction with the media, due to the chains interactional barrier at this length. Moreover, a little deviation of ADp15 from the spherical shape was also effective. Factually, all the considerations in this work clearly showed that the thermoresponsive and aggregation behavior of the PNIPAm occurs at any chain length from the interactional aspects. Nevertheless, the PNIPAm exhibited the chain length limitation to dedicate the conformational change, particularly the coil-to-globule transition, with increasing the temperature. The deep insight into the effects of chemical and physical features of the PNIPAm on its thermal responsivity and aggregation phenomenon provides the ability to design efficient thermoresponsive copolymers to exhibit special characteristics in aqueous solutions. In fact, understanding the effect of temperature and chain length of thermoreversible polymers on the polymer interfacial and hydrodynamic behaviors provides the capability to tune an effective thermoresponsive copolymer with specified features. Designing a block copolymer consisting of an inner water-soluble block and PNIPAm end blocks with a thermal transition temperature and selfassembly behavior is one of our future goal.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).