Analysis of the internal motions of thermoresponsive polymers and single chain nanoparticles

Abstract

Data-driven techniques, such as proper orthogonal decomposition (POD) and uniform manifold approximation & projection (UMAP), are powerful methods for understanding polymer behavior in complex systems that extend beyond ideal conditions. They are based on the principle that low-dimensional behaviors are often embedded within the structure and dynamics of complex systems. Here, the internal motions of a thermoresponsive, LCST polymer are investigated for two cases: (1) the coil-to-globule transition that occurs as the system is heated above its critical temperature and (2) intramolecularly crosslinked, single chain nanoparticles (SCNPs) both above and below the critical temperature (T_C). Our results demonstrate that POD can successfully extract the key features of the dynamics for both polymer globules and SCNPs. In the globular state, our results show that the relaxation modes are distorted relative to the coil state and relaxation times decrease upon chain collapse. After randomly crosslinking a globule to produce a SCNP, we observe a further distortion of the relaxation modes that depends strongly upon the particular set of monomers that are crosslinked. Yet, different sets of crosslinked monomers produce similar relaxation times for the SCNP. We observe that for SCNPs below the critical temperature, the relaxation times decrease with increasing crosslink density while above the critical temperature, they increase as crosslink density increases. Finally, using UMAP we categorize the local structure of SCNPs and examine the influence of the local structure on SCNP relaxation dynamics.

---

1 Introduction

Thermoresponsive polymers have attracted considerable attention in the soft matter community, in large part because of their ability to add functionality to soft materials^1–3 – especially for biomedical applications.^4–6 A commonly exploited behavior is the so-called “coil-to-globule” transition^7–9 that occurs for many lower critical solution temperature (LCST) polymers. Polymers that undergo a coil-to-globule transition collapse above their critical temperature (T_C) due to increasingly unfavorable interactions with their environment as temperature increases. Poly(N-isopropylacrylamide) (PNIPAM) is a widely studied polymer¹⁰ that exhibits a critical temperature in water near T_C ≈ 32 °C, although the exact value of T_C depends on factors that can include the polymer's concentration, terminal group chemistry, architecture, and others.^11–14 Additionally, the presence of a second solvent such as methanol or ethanol can lead to more complex behaviors like co-nonsolvency – a phenomenon in which mixtures of compatible solvents for a polymer result in poor solvation, and which can be utilized to further tune the response of the polymer to its environment.^15,16 Although PNIPAM is widely studied, several alternative thermoresponsive polymers have also received an increasing amount of attention due to their highly tunable critical temperatures, such as with poly[oligo(ethylene glycol) methyl ether methacrylate] (POEGMA),^17,18 or because of a lack of thermal hysteresis, such as for poly(N-cyclopropylacrylamide) (PNCPAM).^19,20

The coil-to-globule transition for PNIPAM has been extensively studied with experimental measurements, computer simulations, and theory.^8,21,22 For instance, Wu and Wang²¹ observed the single chain collapse of PNIPAM by performing light scattering measurements on extremely dilute solutions of high molecular weight chains. They observed a small, monotonic decrease in both the radius of gyration (R_g) and hydrodynamic radius (R_h) as temperature increased. Across a very small temperature window, near T_C, they observed the PNIPAM chains collapse significantly with a shape factor R_g/R_h that decreased from approximately 1.5 (linear chain, good solvent) to 0.56, before recovering to approximately 0.8 (solid sphere/globule). Several distinct states of the PNIPAM chains have been identified during the coil-to-globule transition based on the value of the shape factor. As temperature increases towards T_C, a PNIPAM coil collapses into a crumpled coil near T_C before further collapsing into a molten globule (R_g/R_h = 0.56).²³ The molten globule is thought to be reminiscent of a globule, but with a rough surface that leads to a shape factor below 0.8. Finally, with increasing temperature the PNIPAM chain fully collapses into a globule.

In the collapsed state, individual chains that contain a small amount of di-functional monomers, such as N,N′-methylene bisacrylamide (MBA), can be crosslinked to form single chain nanoparticles (SCNPs),²⁴ which have been suggested as unique materials for drug delivery, imaging, and as nanoreactors among many others.^5,25,26 Neutron scattering techniques have been instrumental for studying the structure and dynamics of SCNPs.^27–29 Small-angle neutron scattering (SANS) measurements by González-Burgos and coworkers examined the structure of poly(methyl methacrylate) (PMMA) SCNPs dissolved in deuterated N,N-dimethyl formamide (DMF) with approximately x_c = 28% of the monomers crosslinked. Upon crosslinking, the polymers became more compact as indicated by a decrease in the radius of gyration (R_g) and increase of the fractal dimension of the particles, where ν is the Flory exponent describing the scaling of R_g with chain length N. However, despite this reduction in the dimensions of the chains upon crosslinking, the authors noted that for many different polymer chemistries (including PNIPAM) and crosslinking densities (x_c ≈ 10 to 30%), SCNPs were not globular in conformation (i.e., ν > 1/3).^27,30 In a similar study, Fisher and coworkers characterized polystyrene and poly(ethyl hexyl methacrylate) SCNPs in deuterated tetrahydrofuran (THF) with SANS, and observed that even at low crosslinking densities, intramolecularly crosslinked polymers develop interfaces making them more akin to particles than polymer chains.²⁸ In fact, González-Burgos et al. observed the appearance of interparticle structure factor peaks in their SANS measurements, indicative of the particle-like character of SCNPs. However, the transition from polymers to particles is gradual as x_c increases. In contrast with measurements from Pomposo et al.,³⁰ Fisher et al. observed a transition from a random coil conformation to a globular conformation as the crosslink density of their SCNPs increased beyond 1% – illustrating potential differences in SCNP structure depending on the synthetic approach and other factors.^25,31,32

Beyond studies of the structure of thermoresponsive polymers and their associated SCNPs, a handful of studies in the literature have sought to measure their internal motions,^33,34 dynamics,^35–37 diffusion,³⁸ and swelling behavior³⁹ using a combination of theory/simulation and experimental characterizations. González-Burgos et al.²⁷ further characterized PMMA SCNPs using neutron spin echo spectroscopy and observed that internal friction significantly influenced the dynamics of both polymer precursors and SCNPs synthesized from them, with the effects of internal friction becoming stronger with crosslinking. The strong internal friction present in SCNPs has been attributed to non-deformable structural features within the interior of the nanoparticle.²⁷ However, despite the increasing attention being paid to understanding the dynamical behavior of SCNPs and related systems, as of today there are myriad opportunities to better understand their behaviors.

Computer simulations offer an attractive means with which to study the internal motions and relaxation dynamics of thermoresponsive polymers and the SCNPs that they may comprise. In addition, classical theories of polymer dynamics, such as those developed by Rouse and Zimm, provide a means to theoretically understand the internal motions of individual polymers chains in the context of their relaxation modes. The foundation of this approach is an equation of motion for a monomer in a polymer chain:⁴⁰

where D_i is a vector containing the diffusivity of monomer i, ϕ(r_i) is an external potential acting on the monomer, ζ is the friction coefficient for the monomer, and f_i is thermal noise. In the Rouse model, it is assumed that ϕ(r_i) is a harmonic potential and that monomer motion is spatially isotropic (∇D_i = 0). Solutions to the Rouse model produce a series of p sinusoidal modes with relaxation times that scale as τ_p ∼ (N/p)^1+2ν, while the consideration of hydrodynamic interactions (i.e., the Zimm description) produces similar results, but with a scaling τ_p ∼ (N/p)^3ν. Furthermore, the sinusoidal modes that result from an analytical solution of eqn (1), with suitable boundary conditions, can be directly applied to the trajectories of particle based simulations (i.e., “Rouse mode analysis”) to extract relaxation times for polymer systems of interest.^41–44 However, strictly speaking, the set of modes commonly used for Rouse mode analysis are only valid for high molecular weight linear/cyclic polymers under the influence of a harmonic potential. Despite this, the approach works well for polymers under the influence of a finite-extensible nonlinear elastic (FENE) potential, relatively short chains (N ≥ 20), and for certain non-linear architectures (e.g., poly[n]catenanes).⁴⁴ Nevertheless, complementary approaches are required for studying the relaxation dynamics of other classes of systems, such as polymer-grafted nanoparticles or molecular bottlebrushes, with particle-based simulations. Data-driven approaches such as proper orthogonal decomposition (POD)^45–47 can produce quantities such as relaxation modes and relaxation times that are analogous to those which can be derived analytically from classical models or from using Rouse mode analysis on suitable systems. Moreover, they can produce these quantities even for situations in which a governing equation (e.g., the Langevin equation for the Rouse model) is not known, for example, because the potential ϕ(r_i) is unknown or ill-defined. Such approaches rely on snapshots of monomer coordinates from high-fidelity simulations.

Motivated by the need for alternative approaches for studying more complex systems, as well as the opportunities to better understand thermoresponsive polymer/SCNP dynamics with computer simulations, we performed energy-conserving dissipative particle dynamics (eDPD) simulations of a thermoresponsive, LCST polymer under dilute solution conditions. We applied POD to analyze the relaxation dynamics and internal motions of the polymer in two cases: (i) the coil-to-globule transition as the polymer is heated above T_C, and (ii) a SCNP above and below T_C as a function of chain length and crosslink density, and which was created by randomly crosslinking the collapsed globule above T_C. As we will demonstrate below, the modes used for traditional Rouse mode analysis of linear polymers are not applicable to the case of a polymer globule or SCNP. In addition, numerical solutions of the Rouse model in the case of SCNPs are only approximations of the true modes that are present in our simulations due to the presence of strong non-bonded interactions between monomers in the interior of the SCNP. Such numerical solutions are also highly dependent on which pairs of monomers within the SCNP are crosslinked. Our results show that the relaxation times of the polymer decrease as it collapses above T_C, and follow the predictions of the Zimm model. However, distortions of the relaxation modes occur for the polymer globule signifying distinct monomer displacement patterns from a free coil below T_C. Furthermore, larger distortions of the modes occur upon crosslinking the polymer to form a SCNP. Finally, we will demonstrate the interplay between local structure of the SCNP and the variation of the relaxation times of SCNPs with x_c and T by categorizing the internal structure of the SCNPs with uniform manifold approximation & projection (UMAP).

2 Model and methods

2.1 Energy-conserving dissipative particle dynamics

Thermoresponsive polymer solutions were modeled with energy-conserving dissipative particle dynamics (eDPD) simulations using a GPU-accelerated, in-house code (gDPD). GPU acceleration was performed with the NVIDIA CUDA libraries, and simulations were performed on NVIDIA Tesla V100 GPUs.

eDPD can be thought of as an extension to the more commonly used DPD technique, with additional parameters that can account for thermal conduction within the system and its effect on particle interactions.^48,49 We chose to use eDPD instead of DPD because it has been demonstrated that eDPD correctly reproduces the correct temperature dependence of quantities such as the kinematic viscosity and particle diffusivity, whereas DPD does not.⁴⁸ Because the focus of this article is on the dynamic properties of thermoresponsive polymers, it is essential to ensure that such quantities are correctly described. It has been used successfully in previous studies to simulate the collapse of LCST polymers.^49,50 Finally, eDPD naturally describes processes that occur experimentally, such as the variation of interaction strengths as the system is heated from below T_C to above T_C – thus capturing both the effect of faster particle motions as well as increasingly unfavorable interactions as T increases.

With DPD, atoms are coarse-grained into fluid elements (DPD beads) that may represent several atoms. Each DPD bead interacts with all other neighbors within a cutoff distance r_c = 1.0. This sets the fundamental length scale of the simulations. It has been estimated that r_c ≈ 1 nm for typical values of the DPD parameters described below.^51,52 We assume that all DPD beads have unit mass, and kept the density of DPD beads in the system fixed at ρ = 3r_c⁻³. Interactions between particles i and j are described by a conservative (F^C_ij), random (F^R_ij), and dissipative force (F^D_ij):

F^C_ij = a_ijw^R(r_ij)r̂_ij, (2)

F^R_ij = σ_ij(Δt)^−1/2w^R(r_ij)θ_ijr̂_ij, (3)

and

F^D_ij = −γ_ijw^D(r_ij)(r̂_ij·v_ij)r̂_ij, (4)

where r_ij = r_i − r_j, r_ij = |r_ij|, r̂_ij = r_ij/r_ij, and v_ij = v_i = v_j. In addition, θ_ij is a uniform random number with zero mean and unit variance that is uncorrelated for different pairs of particles and at different times. w^D(r_ij) = [w^R(r_ij)]² are weighting factors for the dissipative and random/conservative forces, respectively. a_ij is the strength of the conservative interaction between particles i and j, and becomes larger as the interaction becomes more unfavorable. Setting a_ij = a₀ for all particles i and j, regardless of their chemical species, produces good solvent conditions in DPD simulations of polymer solutions. The coefficients σ_ij and γ_ij, which are related to the strength of thermal fluctuations and energy dissipation, are discussed in more detail below. The combination of the random and dissipative forces forms a thermostat for the system.

Beyond the three DPD forces, adjacent monomers within a polymer chain are connected to each other through forces derived from the finite-extensible nonlinear elastic (FENE) potential, i.e., F^F_i,i+1 = −∇U_F(r_i,i+1), where

if r_i,i+1 < r_max, and U_F → ∞ otherwise. The bond strength k_F = 50ε₀/r_c², r_max = 2r_c, and r_eq = 0.7r_c where ε₀ is the fundamental energy unit in DPD (= k_BT_r with T_r as the temperature of a surrounding thermal reservoir). The combination of eqn (2)–(5) cannot simulate entanglements, and the addition of slip-springs^53–55 or a segmental repulsive potential (SRP)^56,57 is required to correctly incorporate these behaviors. For that reason, our results do not account for any entanglement effects which could become important within the confined interiors of polymer globules or SCNPs.

To correctly model a thermoresponsive system, in which the system may be heated or cooled while developing corresponding changes in the interaction strengths between the components, it is most convenient to augment the traditional DPD approach to explicitly consider energy flow in the system. Under this approach, each DPD bead i is assigned a temperature T_i, which is governed by an equation of motion:^48,49

where C_v is the heat capacity of the system and Q_i is the net heat flux of particle i, which is the sum of the collisional (Q^C_ij), viscous (Q^V_ij), and random (Q^R_ij) heat fluxes, given by:and

Q^R_ij = υ_ijw^R(r_ij)(Δt)^−1/2ζ_ij, (9)

where ζ_ij is a uniform random number of zero mean and unit variance. In the framework of eDPD, the fluctuation–dissipation theorem requires the random coefficient be related to the dissipative coefficient, σ_ij² = 4γ_ijk_BT_iT_j/(T_i + T_j), meaning that unlike in regular DPD, σ_ij and γ_ij must be calculated at each time step for each pair of interacting particles. Additionally, the collisional and random heat fluxes are coupled through their leading coefficients, υ_ij² = 2k_Bk_ij. The thermal conductivity k_ij = C_v²κ(T_i + T_j)²/4k_B. The heat capacity and heat friction coefficient were chosen to be C_v = 1000ε₀ and κ = 0.001, respectively. The equations of motion for the particles and their heat fluxes were integrated using the velocity-Verlet algorithm with an integration time step Δt = 0.01τ₀. Additional details can be found in a previous publication.⁴⁷

Simulations of single, thermoresponsive polymers with degrees of polymerization N were performed in systems with periodic boundary conditions and sizes of L_x = L_y = L_z = 32r_c, which is at least 4× larger than the radius of gyration of the largest polymer (N = 200). We focused on single chain simulations to eliminate any potential confounding effects caused by aggregation of the polymers above their LCST. The polymers were first generated at a random position in the system, and equilibrated for 50 000 time steps at which point the radius of gyration R_g saturated at a constant value, and the average system temperature saturated at where n_T is the total number of DPD beads in the system. After equilibration, polymers were simulated for a minimum of 4 × 10⁶ time steps to acquire sufficient statistics regarding their size, conformation, and relaxation dynamics. To mitigate any potential for statistical deviations to influence the results, two independent runs were made for each system. We observed only slight differences between the two runs, as described in more detail below.

To simulate the coil-to-globule transition, we took the approach of Karniadakis et al.^48,49 and modeled the interaction strength between a monomer bead and solvent bead, a_MS, according to:

where Λ₀ = −10ε₀/r_c, ΔΛ = 20ε₀/r_c, μ = 300k_B/ε₀, and T_C = 1ε₀/k_B. These parameters produce favorable interactions for T < T_C and unfavorable interactions above T_C, with a sharp transition between the two regimes at T = T_C. We simulated thermoresponsive polymers at two temperatures (T = 0.8ε₀/k_B and T = 1.4ε₀/k_B) to understand their behavior in the coil and globule states, respectively. We kept interactions between same particle types fixed at all temperatures, i.e., a_SS = a_MM = 25ε₀/r_c. From eqn (10), for T = 0.8, a_MS = 15ε₀/r_c while a_MS = 35ε₀/r_c for T = 1.4ε₀/k_B. Thus, at T = 0.8ε₀/k_B, interactions between the polymer and solvent are more favorable than monomer–monomer or solvent–solvent interactions, which we will discuss further in the Results and discussion section. To increase the system temperature from T = 0.8ε₀/k_B to 1.4ε₀/k_B at time t, we set the temperature of the thermal reservoir T_r to the desired temperature, which resulted in an additional thermal flux into the system for each particle: Q^S_i(t) = λC_v[T_r(t) − T_i(t)], with λ = 0.1. This additional heat flux for each particle is added to the right hand side of eqn (6).

2.2 Proper orthogonal decomposition

The relaxation modes and relaxation times of the polymers were computed by proper orthogonal decomposition (POD), also referred to as principal component analysis (PCA) depending on the context.⁵⁸ In this approach, snapshots of the monomer coordinates are stored every Δs = 500 time steps, separately for each value of T_r, during the course of the simulation, and post-processed to determine the distance between each monomer and the center of mass of the molecule as a function of time, r̃_i(t). Three covariance matrices M^q are constructed for each component q = x, y and z of r̃_ij:where s is the index of the snapshot and S is the total number of snapshots. The time at which each snapshot is taken t = s(Δs)Δt. M^q is then factored by singular value decomposition (SVD), i.e., M^q = UΣV. Each column of the left singular matrix U contains the eigenfunction ψ_p(n) that describes the fundamental displacement patterns of the monomers (i.e., the modes), ordered left-to-right from the lowest mode (p = 1) to the highest (p = N). The singular value matrix Σ contains the eigenvalues of the factorization along the diagonal, appearing in descending order of their magnitude. The singular values, normalized by Tr(Σ), describe the relative amount of variance in the monomer motions each mode captures. The right singular matrix, V, describes the temporal evolution of the monomer coordinates and is not used for the analysis of the polymer dynamics. However, because the dynamics of a free polymer are expected to be isotropic in the x, y, and z-directions, the eigenfunctions from the three covariance matrices are averaged together to reduce statistical fluctuations.

To extract relaxation times of the polymer, normal coordinates X_p(t) for each mode p are computed by projecting the coordinates of each monomer on to the basis describing the motions of that mode, i.e., X_p(t) = ψ_p(n) × [x(t)y(t)z(t)]. The temporal autocorrelation of the normal coordinates C_p(t) is fit to a stretched exponential function:

Typically, β ≈ 1, although if the dynamics of the polymer are sufficiently heterogeneous β < 1. A detailed description of the POD approach to analyzing polymer dynamics is provided elsewhere.^45–47

3 Results and discussion

3.1 The coil-to-globule transition

We first modeled the collapse of a single, thermoresponsive chain as the temperature of the external reservoir was increased from T_r = 0.8ε₀/k_B to T_r = 1.4ε₀/k_B for chain lengths N = 20, 50, 100, and 200. As T_r increased, a net heat flux developed that added energy to the system such that 〈T〉 ≈ T_r. At the same time, eqn (10) produced an increase of the polymer–solvent (PS) interaction parameter from a_MS = 15ε₀/r_c (T_r = 0.8ε₀/k_B) to a_MS = 35ε₀/r_c (T_r = 1.4ε₀/k_B), relative to solvent–solvent (SS)/monomer–monomer (MM) interaction strength of a_SS = a_MM = 25ε₀/r_c. As the temperature was increased above T_C, the polymer coils collapsed into globules within 5–10τ₀.

The average radius of gyration of the polymer across 4 × 10⁶ time steps (40 000τ₀), and averaged between two independent simulations, is plotted in Fig. 1 for T < T_C (T_r = 0.8ε₀/k_B, black) and T > T_C (T_r = 1.4ε₀/k_B, red). Representative snapshots of the polymer conformation are shown in the figure for a chain of N = 100 monomers. The two independent simulations yielded almost identical values for R_g, as indicated by vanishing size of the error bars at each data point. The dashed lines in the figure denote the scaling exponent of R_g, i.e. R_g ∼ N^ν, where ν = 0.63 for T < T_C and ν = 0.30 for T > T_C – indicating that the expected scaling relationships between R_g and N for a coil (ν ≈ 3/5) and a globule (ν ≈ 1/3) are essentially reproduced by the eDPD simulations.

Fig. 1 Average radius of gyration, R_g, plotted as a function of degree of polymerization, N. Circles correspond to coil conformations below T_C while squares correspond to globular conformations above T_C. The rendered polymers depict representative coil (top) and globule (bottom) conformations for N = 100. Error bars, corresponding to one standard deviation in the value of R_g are smaller than the size of the points.

We next focused on analyzing the relaxation dynamics of the coil and globule conformations of the polymer using POD. Shown in Fig. 2 are the first three relaxation modes ψ_p(n), determined from POD, at temperatures (A) above and (B) below the critical temperature for a chain of N = 200 monomers. The modes for other chain lengths (not shown) were similar. In this figure, the modes extracted by POD in the x-, y-, and z-directions were averaged into a single mode. Shown as dashed lines in the figure are the Rouse modes calculated from the numerical solution of eqn (1) for a harmonic potential with ∇D_i = 0. The POD modes represent fundamental monomer displacement patterns that underlie the monomer motions in the eDPD simulations. In particular, p = 1 (black) describes the relaxation of the chain ends relative to each other (i.e., motions about the central monomer of the chain). As p increases, the modes describe the motions/relaxations of increasingly smaller regions in the interior of the chain.^59,60

Fig. 2 Displacement patterns (i.e., eigenfunctions) of the first three relaxation modes of a single, free polymer chain (N = 200) determined from proper orthogonal decomposition (POD, points) compared to the analytical result from the Rouse Model (dashed lines) for (A) T < T_C and (B) T > T_C.

As demonstrated in Fig. 2(A), there is good quantitative agreement between the Rouse modes and the POD modes for coils below the critical temperature, without the need for any adjustable parameters. However, the modes in Fig. 2(B) display small distortions above the critical temperature, implying that the motions of the monomers in the globular state are distinct from those of a random coil. These distortions likely arise in a couple of ways due to the local environment inside the globule. First, as compared to the Rouse model, the potential ϕ(r_i) experienced by monomers in the globule consists of both the bonded and non-bonded interactions, arising from the FENE and conservative DPD potential, respectively. For the largest chain length simulated (N = 200), R_g ≈ 2r_c, implying that a large number of monomer pairs exist within a cutoff distance r_c, producing non-negligible contributions to the dynamics (e.g., in eqn (1)). Second, the collapse of the coil into a globular state may introduce physical restrictions on the motions of the monomers such that underlying monomer displacement patterns differ from those of a free coil below T_C. In particular, above T_C, the three modes show a common behavior of displacements moving towards the interior of the chain, as indicated by the movement of the maxima of ψ_p(n) away from the chain ends (i.e., n = 1 and n = 200). Similar restrictions on internal motions of polymer globules were observed experimentally by Dai et al.^33,34 using laser light scattering.

In Fig. 3 the fraction of variance in the monomer motions captured by each mode, i.e., where λ_p is the singular value of mode p, is plotted for the first 6 modes below (black) and above (red) T_C. While the first 4 modes capture the majority of the polymer dynamics below T_C, we observe many more modes are required above T_C to fully describe the motions of monomers in the globular state. For instance, while the lowest p = 1 mode captures more than 70% of the dynamics below T_C, the lowest mode only captures approximately 20% of the dynamics above T_C – necessitating a larger number of modes to fully describe the dynamics. Nevertheless, the relaxation times at both temperatures can still be extracted from fitting C_p(t) to eqn (12) with β = 1 (solid lines, Fig. 4). Shown in the inset of Fig. 4 is the scaling of the relaxation time of mode p = 1 (τ₁) as a function of chain length. According to the Zimm description of dilute polymer dynamics in the presence of hydrodynamic interactions, τ₁ ∼ N^3ν. The scaling exponents extracted from the eDPD simulations are shown as dashed lines in the figure. While the simulations slightly under-predict this value for T < T_C (i.e., 3ν = 1.6 rather than 1.8), for T > T_C the result exactly matches the Zimm prediction (3ν = 0.9 for ν = 0.3). Although solvent is not expected within the interior of the globule, hydrodynamic effects at the globule's surface may lead to similar dynamics as a coil in dilute solution. Thus, although the motions of the monomers are slightly distorted from the ideal Rouse modes, the relaxation times are in accord with the classical description.

Fig. 3 Fraction of the variance in the monomer motions captured by modes p = 1 to 6 for a single polymer (N = 200) in the coil state (T < T_C, black) and globule state (T > T_C, red).

Fig. 4 Normal mode temporal autocorrelation function C_p(t) plotted for the lowest mode p = 1 as a function of time above T_C. Red curves are fits to eqn (12) with β = 1 for (left to right) N = 20, 50, 100, and 200. The inset presents the scaling of the longest relaxation time τ₁ with degree of polymerization N. The dashed curves correspond to the scaling exponent for T < T_C (black) and T > T_C (red).

3.2 Single chain nanoparticles (SCNPs)

After determining the behavior of a thermoresponsive coil both above and below its critical temperature, we focused on understanding single chain nanoparticles (SCNPs) that were produced by randomly crosslinking x_c% monomers within a globule at T > T_C. This process is demonstrated in Fig. 5 where a globule is randomly crosslinked by connecting random pairs of monomers in close proximity using the FENE potential with identical parameters to those described above in Section 2. SCNPs were then simulated at T_r = 1.4ε₀/k_B and T_r = 0.8ε₀/k_B. As will be discussed in more detail below, the random crosslinking procedure resulted in modes ψ_p(n) that differed between different particles/independent simulations.

Fig. 5 Simulation snapshots depicting the process of forming a single-chain nanoparticle (SCNP) from a collapsed globule (N = 200) at T > T_C by crosslinking x_c = 10% random, non-adjacent monomers, followed by an optional re-swelling of the particle at T < T_C.

Representative displacement patterns of the first mode, ψ₁(n) are shown in Fig. 6 for a SCNP with N = 100 and a crosslink density of x_c = 15% (red line), compared to ψ₁(n) for a free chain/globule (black line) and a numerical solution of eqn (1) (blue dashed line) for an identical set of crosslinks. Crosslink points are superimposed on ψ₁(n) as filled symbols. Note that different pairs of crosslinked monomers will produce different values for ψ_p(n), but nearly identical values of C_p(t), as shown in Fig. S1 and S2 in the ESI.† The difference between the numerical solution and the POD modes demonstrates the effect that local interactions, for example from the conservative force, have on the monomer motions. For instance, although the numerical Rouse mode is very similar to the POD mode in terms of the locations of maxima and minima of the function, the amplitudes of local regions of the modes can differ between the two approaches. At the top of Fig. 6 are representative snapshots of the SCNP along different planes with crosslinked monomers rendered in red and non-crosslinked monomers as translucent beads. Upon crosslinking, the relaxation modes become more distorted than those for a non-crosslinked globule. In addition, we observe that regions of ψ₁(n) with many adjacent crosslinks (e.g., near n = 30 and n = 50 in the figure) show less variation than those regions between crosslink points, implying these regions may be relatively spatially static compared to regions of the polymer farther along the contour (e.g., near n = 80). This can be visualized with the snapshots in Fig. 6 where we observe regions of crosslinked monomers in close proximity to one another, and is reminiscent of the immobile domains/clusters proposed by González-Burgos et al.²⁷ on the basis of their neutron spin echo measurements. Although ψ_p(n) varies with the specific monomers that are crosslinked within a SCNP, independent simulations found that both C₁(t) and τ₁ were very similar between different sets of crosslinked monomers. However, certain sets of crosslinked monomers resulted in autocorrelation functions that did not fully decay to zero at intermediate times, the origin of which will be discussed in more detail below.

Fig. 6 Monomer displacement patterns for a SCNP (N = 100, x_c = 15%) above T_C. The snapshots show the SCNP along several different orientations, with crosslinked monomers shaded in red and non-crosslinked monomers as translucent beads. The plot compares ψ₁(n) for a globule (x_c = 0%, black) to ψ₁(n) computed for a SCNP (x_c = 15%, red) using POD. The result of numerically solving the Rouse model for an ideal SCNP with identical crosslinks is shown in blue. Crosslink locations are superimposed on ψ₁(n) as discrete points.

The normal coordinate autocorrelation function C_p(t) is plotted in Fig. 7(A) at temperatures above (red squares) and below (black circles) the critical temperature for a fixed chain length of N = 200 and varied crosslink density x_c. As indicated by the directions of the arrows, curves from left to right correspond to increasing x_c for T > T_C, while for T < T_C, C_p(t) shifts to the left as x_c increases. However, in all cases C_p(t) is well-described by eqn (12) with β = 1 regardless of temperature or crosslink density. The variation of C_p(t) with x_c implies that relaxation times increase in the crosslinked globular state as x_c increases, whereas they decrease in the crosslinked coil state as x_c increases. The relaxation times determined by eqn (12) are plotted in Fig. 7(B), where we also observe that τ₁ appears to decrease as x_c increases below the critical temperature, while τ₁ appears to increase and saturate at a constant value as x_c increases. For the lowest crosslink density of x_c = 2.5%, the relaxation times at the two temperatures differ by an order of magnitude and differ by approximately 500τ₀ for x_c = 20%. Although we did not simulate SCNPs for x_c > 20%, we expect that the two branches in Fig. 7B should converge as x_c → 100%.

Fig. 7 (A) Temporal autocorrelation function C_p(t) (p = 1) for SCNPs (N = 200) as a function of x_c for T < T_C (black circles) and T > T_C (red squares). The lines are fits to eqn (12) with β = 1. The arrows denote the variation of x_c. (B) Relaxation times of the first mode, obtained from the data in (A), as a function of x_c, plotted for SCNPs below T_C (black circles) and above T_C (red squares). The lines are guides for the reader.

An explanation for the different variation in τ₁ above and below the critical temperature can be formulated on the basis of the snapshots of a SCNP (N = 200, x_c = 2.5%) in Fig. 8, and those shown in Fig. 5 for N = 200 and x_c = 10%. First, the SCNP undergoes an additional collapse for T > T_C due to the poor solvent condition as shown in Fig. 8, resulting in shorter relaxation times above the critical temperature. Second, the snapshots in Fig. 8A demonstrate how distant crosslink points along the polymer contour connect distant points together, leading to a more compact structure than if the chain was not crosslinked. As a result, as more crosslinks are added and the chain becomes increasingly compact, the relaxation times decrease by a factor of 2 to 3 as x_c increases from 2.5% to 20%. In this case, both the conformation of the polymer and constraints on the monomer motions vary with x_c. In contrast, for T > T_C, the size of the globule and the SCNP are similar to each other (cf., Fig. 5), and increasing x_c leads almost exclusively to additional constraints on monomer motions – leading to correlations in X_p(t) over longer spans of time and a small increase in τ₁.

Fig. 8 Representative snapshots of a SCNP (N = 200, x_c = 2.5%) in the (left to right) XZ-plane, XY-plane, and YZ-plane for (A) T < T_C and (B) T > T_C. All snapshots are at the same magnification with identical pairs of crosslinked monomers. For reference, the radius of the SCNP in (B) is approximately 2r_c. Crosslinked monomers are shaded in red, while non-crosslinked monomers appear translucent.

Finally, although C_p(t) decayed exponentially for all polymer coils and globules, we observed that for certain SCNPs, particular sets of crosslinked monomers showed residual amounts of correlations in C_p(t) at intermediate times. This behavior was observed almost exclusively in short chain SCNPs with N = 50, and did not appear to depend on x_c. However, most systems showed full exponential decay of C_p(t). For those systems with residual correlations, the degree of correlation varied between the x-, y-, and z-directions within the SCNP, implying that the residual correlations may be related to the interior structure of the SCNP. An example of this behavior is shown in Fig. 9A for a SCNP with N = 50 and x_c = 10%, where the black circles represent the average C_p(t) function, and the colored points show C_p(t) in the x-, y-, and z-directions. r_0x, r_0y, and r_0z denote the amount of residual correlations in each of the three directions, where for full exponential decay of C_p(t), r₀ = r_0x = r_0y = r_0z = 0. Consideration of quantities such the distribution of crosslinks or the distances between pairs of crosslinked monomers along the polymer contour did not produce any meaningful connection between the SCNP structure and the behavior of C_p(t). Instead, we adopted the data-driven approach of Reinhart^61,62 for constructing sequence-morphology phase maps of sequence-defined polymers to characterize the local structure of our SCNPs using the structure of triads of DPD particles {i, j, k} In this approach, uniform manifold approximation & projection (UMAP) is performed on structural histograms to embed systems in a “structure space” to quantitatively determine whether SCNPs with residual correlations in C_p(t) were structurally similar (i.e., in close proximity in structure space). Full details of implementing this approach are provided in the ESI.† The embedding of SCNP local structure in structure space Z = (Z₀, Z₁, Z₂) using UMAP is shown in Fig. 9B in the Z₁ and Z₂ directions, with red crosses representing systems that show residual correlations and black circles representing systems that fully decay exponentially. This projection in Z groups SCNPs on the basis of the distance between particles j and k (d_jk) within the SCNP, and the sum of the bond lengths between all three monomers (i.e., l_ijk = d_ij + d_ik = |r_j − r_i| + |r_k − r_i|). We observe that systems with r₀ > 0 formed a single group in Z that was separate from those systems with r₀ = 0, further supporting the hypothesis that residual correlations in C_p(t) were the result of the local structure of the SCNP that resulted from random crosslinking. Projections of the UMAP embedding along the Z₀ direction are provided in the ESI,† but did not show as clear of a grouping between the two classes of systems.

Fig. 9 (A) Example autocorrelation function for a SCNP (N = 50, x_c = 10%) that shows residual correlations r₀ (dashed lines) at intermediate times (t ≈ 10² to 10³τ₀). The black circles are an average of the autocorrelation functions in the x- (red squares), y- (blue diamonds), and z-directions (green triangles). (B) Slice of the manifold obtained by UMAP in the Z₁ and Z₂ directions, showing the grouping/structural similarities between those points with residual correlations (red crosses) and those which fully decay (black circles).

4 Conclusions

In summary, we have utilized two data-driven analysis techniques – proper orthogonal decomposition (POD) and uniform manifold approximation & projection (UMAP) – to study the relaxation dynamics of a thermoresponsive, LCST polymer during (i) the coil-to-globule transition and (ii) when crosslinked to form a single-chain nanoparticle (SCNP). Below the critical temperature T_C, individual polymer coils display relaxation modes that exactly match those predicted by analysis of a Rouse chain (cf.eqn (1)). The relaxation times in this regime roughly scale with N as predicted by the Zimm model for dilute polymer solutions. However, above the critical temperature we observe a distortion of these modes and the relaxation modes describe processes that occur towards the interior of the polymer contour. Nevertheless, when these modes are used to compute the normal coordinates X_p(t) and their temporal autocorrelation C_p(t), we obtain relaxation times τ that exactly scale with polymer chain length as predicted by Zimm: τ₁ ∼ N^0.9.

In contrast to the behavior of the polymers during the coil-to-globule transition, when globules are randomly crosslinked to form SCNPs, the modes present in the nanoparticle depend strongly on the set of monomers that are crosslinked. Below T_C, relaxation times decrease with increasing crosslink density x_c while above T_C, they increase. This behavior was attributed to the interplay between the conformation of the SCNP and restriction of monomer motions by the crosslinks. In particular, while the polymer conformation and constraints on motions vary with x_c below T_C, only the constraints on the monomer motions vary with x_c above T_C. The local, internal structure of the SCNP, categorized using UMAP, strongly influenced the decay of C_p(t) for short chain SCNPs.

Looking to the future, techniques such as POD and UMAP can enable investigations of the structure and dynamics of more complex systems and polymer architectures – with an eye towards predicting the macroscopic properties of soft materials. For instance, these techniques may be directly applicable to more structured polymer architectures, such as those which possess secondary structures (e.g., folded proteins), sequence-defined polymers, and micro/nanoplastic contaminants. At the same time, experimental characterizations of the dynamics of thermoresponsive polymers, for example by neutron spectroscopy, are well-warranted to compare the experimental behaviors of thermoresponsive polymers (e.g., PNIPAM) with the behaviors predicted by eDPD simulations.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work made use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Case Western Reserve University. The author thanks Prof. Wesley Reinhart (Pennsylvania State University) for helpful discussions regarding the application of UMAP to interpreting local structural information.

References

1. J. Kim, J. A. Hanna, M. Byun, C. D. Santangelo and R. C. Hayward, Science, 2012, 335, 1201–1205.
2. A. S. Kuenstler, M. Lahikainen, H. Zhou, W. Xu, A. Priimagi and R. C. Hayward, ACS Macro Lett., 2020, 9, 1172–1177.
3. S.-J. Jeon and R. C. Hayward, Soft Matter, 2020, 16, 688–694.
4. N. A. Cortez-Lemus and A. Licea-Claverie, Prog. Polym. Sci., 2016, 53, 1–51.
5. A. Bordat, T. Boissenot, J. Nicolas and N. Tsapis, Adv. Drug Delivery Rev., 2019, 138, 167–192.
6. F. Doberenz, K. Zeng, C. Willems, K. Zhang and T. Groth, J. Mater. Chem. B, 2020, 8, 607–628.
7. J. Yu, Z. Wang and B. Chu, Macromolecules, 1992, 25, 1618–1620.
8. X. Wang, X. Qiu and C. Wu, Macromolecules, 1998, 31, 2972–2976.
9. D. S. Simmons and I. C. Sanchez, Macromolecules, 2008, 41, 5885–5889.
10. A. Halperin, M. Kröger and F. M. Winnik, Angew. Chem., Int. Ed., 2015, 54, 15342–15367.
11. Y. Pei, J. Chen, L. Yang, L. Shi, Q. Tao, B. Hui and J. Li, J. Biomater. Sci., Polym. Ed., 2004, 15, 585–594.
12. K. N. Plunkett, X. Zhu, J. S. Moore and D. E. Leckband, Langmuir, 2006, 22, 4259–4266.
13. Y. Zhang, S. Furyk, L. B. Sagle, Y. Cho, D. E. Bergbreiter and P. S. Cremer, J. Phys. Chem. C, 2007, 111, 8916–8924.
14. T. E. de Oliveira, D. Mukherji, K. Kremer and P. A. Netz, J. Chem. Phys., 2017, 146, 034904.
15. H. G. Schild, M. Muthukumar and D. A. Tirrell, Macromolecules, 1991, 24, 948–952.
16. Y.-Y. Zhang, X.-M. Jia, R. Shi, S.-J. Li, H. Zhao, H.-J. Qian and Z.-Y. Lu, Macromol. Rapid Commun., 2020, 41, 1900655.
17. J.-F. Lutz and A. Hoth, Macromolecules, 2006, 39, 893–896.
18. J.-F. Lutz, Ö. Akdemir and A. Hoth, J. Am. Chem. Soc., 2006, 128, 13046–13047.
19. D. Roy, W. L. Brooks and B. S. Sumerlin, Chem. Soc. Rev., 2013, 42, 7214–7243.
20. X. Lang, A. D. Patrick, B. Hammouda and M. J. Hore, Polymer, 2018, 145, 137–147.
21. C. Wu and X. Wang, Phys. Rev. Lett., 1998, 80, 4092.
22. L. Tavagnacco, E. Zaccarelli and E. Chiessi, Phys. Chem. Chem. Phys., 2018, 20, 9997–10010.
23. C. Wu and S. Zhou, Phys. Rev. Lett., 1996, 77, 3053.
24. X. Qiu and C. Wu, Macromolecules, 1997, 30, 7921–7926.
25. R. Chen and E. B. Berda, ACS Macro Lett., 2020, 9, 1836–1843.
26. A. Nitti, R. Carfora, G. Assanelli, M. Notari and D. Pasini, ACS Appl. Nano Mater., 2022, 5, 13985–13997.
27. M. González-Burgos, I. Asenjo-Sanz, J. A. Pomposo, A. Radulescu, O. Ivanova, S. Pasini, A. Arbe and J. Colmenero, Macromolecules, 2020, 53, 8068–8082.
28. J. Fischer, L. Han, T. Saito and M. Dadmun, Nanoscale Adv., 2022, 4, 5164–5177.
29. B. Robles-Hernández, P. Malo de Molina, I. Asenjo-Sanz, M. Gonzalez-Burgos, S. Pasini, J. A. Pomposo, A. Arbe and J. Colmenero, Macromolecules, 2024, 57, 4706–4716.
30. J. A. Pomposo, I. Perez-Baena, F. Lo Verso, A. J. Moreno, A. Arbe and J. Colmenero, ACS Macro Lett., 2014, 3, 767–772.
31. S. Liao, L. Wei, L. A. Abriata and F. Stellacci, Macromolecules, 2021, 54, 11459–11467.
32. T. P. Le, L. Cavalcanti, J. P. Tellam and P. Malo de Molina, Biomacromolecules, 2024, 25, 6602–6610.
33. Z. Dai and C. Wu, Macromolecules, 2010, 43, 10064–10070.
34. Z. Dai, T. Ngai and C. Wu, Soft Matter, 2011, 7, 4111–4121.
35. A. Arbe, J. Rubio-Cervilla, A. Alegra, A. J. Moreno, J. A. Pomposo, B. Robles-Hernández, P. Malo de Molina, P. Fouquet, F. Juranyi and J. Colmenero, Macromolecules, 2019, 52, 6935–6942.
36. E. Verde-Sesto, A. Arbe, A. J. Moreno, D. Cangialosi, A. Alegra, J. Colmenero and J. A. Pomposo, Mater. Horiz., 2020, 7, 2292–2313.
37. Y. Ruan, Q. Zou, H. Zhang, Y. Zhang, H. Zhang, W. Wu, H. Liu, J. Yan and G. Liu, Macromolecules, 2024, 57, 6583–6592.
38. B. Mei, A. J. Moreno and K. S. Schweizer, ACS Nano, 2024, 15529–15544.
39. H. Rabbel, P. Breier and J.-U. Sommer, Macromolecules, 2017, 50, 7410–7418.
40. J. T. Padding, Advanced Courses in Macroscopic Physical Chemistry (Han-Sur-Lesse Winterschool 2005), 2005, pp. 21–28.
41. J. T. Kalathi, S. K. Kumar, M. Rubinstein and G. S. Grest, Macromolecules, 2014, 47, 6925–6931.
42. J. T. Kalathi, S. K. Kumar, M. Rubinstein and G. S. Grest, Soft Matter, 2015, 11, 4123–4132.
43. H.-P. Hsu and K. Kremer, Eur. Phys. J.: Spec. Top., 2017, 226, 693–703.
44. P. M. Rauscher, S. J. Rowan and J. J. de Pablo, ACS Macro Lett., 2018, 7, 938–943.
45. C. P. J. Wong and P. Choi, Macromol. Theory Simul., 2019, 28, 1800072.
46. C. A. Miller and M. J. Hore, ACS Polym. Au, 2021, 2, 157–168.
47. R. S. Mukkamala and M. J. Hore, Macromolecules, 2024, 57, 445–455.
48. Z. Li, Y.-H. Tang, H. Lei, B. Caswell and G. E. Karniadakis, J. Comput. Phys., 2014, 265, 113–127.
49. Z. Li, Y.-H. Tang, X. Li and G. E. Karniadakis, Chem. Commun., 2015, 51, 11038–11040.
50. E. Cudjoe, S. Khani, A. E. Way, M. J. Hore, J. Maia and S. J. Rowan, ACS Cent. Sci., 2017, 3, 886–894.
51. F. Lahmar and B. Rousseau, Polymer, 2007, 48, 3584–3592.
52. J. Cheng, A. Vishnyakov and A. V. Neimark, Langmuir, 2014, 30, 12932–12940.
53. M. Langeloth, Y. Masubuchi, M. C. Böhm and F. Müller-Plathe, J. Chem. Phys., 2013, 138, 104907.
54. A. Ramrez-Hernández, F. A. Detcheverry, B. L. Peters, V. C. Chappa, K. S. Schweizer, M. Müller and J. J. De Pablo, Macromolecules, 2013, 46, 6287–6299.
55. J. Schneider, A. Z. Panagiotopoulos and F. Müller-Plathe, J. Phys. Chem. C, 2017, 121, 27664–27673.
56. P. Nikunen, I. Vattulainen and M. Karttunen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 036713.
57. T. W. Sirk, Y. R. Slizoberg, J. K. Brennan, M. Lisal and J. W. Andzelm, J. Chem. Phys., 2012, 136, 134903.
58. S. L. Brunton and J. N. Kutz, Data-driven science and engineering: Machine learning, dynamical systems, and control, Cambridge University Press, 2022.
59. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, 1988, vol. 73.
60. G. R. Strobl, The Physics of Polymers, Springer, 1997, vol. 2.
61. W. F. Reinhart, Comput. Mater. Sci., 2021, 196, 110511.
62. A. Statt, D. C. Kleeblatt and W. F. Reinhart, Soft Matter, 2021, 17, 7697–7707.

---

Footnote

† Electronic supplementary information (ESI) available: Details and results of UMAP procedure. See DOI: https://doi.org/10.1039/d4sm01308e

---