Diffusion dynamics of star-shaped macromolecules in dilute solutions

Polymeric fluids are a unique class of complex fluids that show a plethora of fascinating non-Newtonian behaviours¹, which can be understood with the transport and rheological properties of the fluid. One key challenge of the complex fluid community is understanding how the macroscopic properties of the polymeric fluids arise from microscopic interactions of the macromolecules². The polymer chain dissolved in a solvent can take a multitude of conformations due to its large internal degrees of freedom. The average shape and size are used to characterize a polymer chain geometrically: the radius of gyration is widely used to describe the size and the eigenvalues of the gyration tensor for the shape calibration³. Advances in controlled polymerization have led to synthesizing complex polymer structures like star-shaped, comb-shaped, H-shaped, ring, and many more.⁴. The diffusive dynamics of such complex polymer chains are of fundamental interest in the biophysics community and are ubiquitous in numerous engineering applications. Star-shaped polymers are used for the controlled delivery of drugs in biomedical applications⁵ and as viscosity index modifiers in oil industries⁶. Polyethylene glycol stars are used for protein delivery⁷. Understanding macromolecular diffusion in a biological cell is important for its various functions, such as the movement of plasmids and transport of amino acids ^{8 9}. Hence, the influence of the shape and size of the polymer chain on its diffusive dynamics in solution is essential.

Most studies on the diffusion of complex polymer chains have been done by keeping the length of the polymer chain constant. Using fluorescence microscopy, Robertson et al.¹⁰ have shown a lower radius of gyration and a higher diffusion coefficient for circular DNA molecules than linear ones for the same chain length. Using Brownian dynamics simulations with hydrodynamic interaction, Kanaeda and Deguchi¹¹ have reported higher diffusion coefficients for the ring polymers than for the linear polymers of the same chain length. Hegde et al.¹² have reported similar findings for the ring polymers in comparison to the linear chains for the same chain length by using three different simulation techniques. Singh et al.¹³ have shown, using Multi-particle Collision Dynamics(MPCD), that the center-of-mass diffusion coefficient of the star polymer chains decreases with an increase in their radius of gyration. Hence, when it comes to size, it is clear that the higher the size of the polymer chain, the lower the center-of-mass diffusion coefficient. However, it is difficult to comment on the influence of the shape of the polymer chain on its diffusion from the same polymer chain length study as both shape and size are distinct for different polymer chain architectures¹⁴. The diffusion study of complex polymer chains for the same size case is scanty. Hegde et al.¹² have reported a higher diffusion coefficient for the linear chains than the ring chains for the same radius of gyration using Molecular Dynamics, MPCD, and the Lattice Boltzmann method and also noted that size could not be the only factor that influences the diffusion of the chain. Therefore, the effect of the shape parameter on the center-of-mass diffusion of the polymer chains still remains an open question.

In this work, the effect of the shape parameter on the center-of-mass diffusion of the star-shaped polymer chains in solution is studied in the limit of infinite dilution using a mesoscopic coarse-grained simulation method, namely MPCD. For simulating the Brownian motion of the complex polymer chains in a solution, MPCD is frequently used as it incorporates both thermal fluctuation and long-range hydrodynamic interactions¹⁵. At first, the shape and size of star-shaped polymer chains with different functionality are analyzed using the gyration tensor and compared with linear polymer chains at the same chain length. Subsequently, the translational diffusion of six different types of chains (one linear and five star-shaped chains) with the same radius of gyration is studied using the center-of-mass mean square displacement, followed by their rotational diffusion using the reorientation correlation function. Finally, the diffusion study is correlated to the shape characterization study in order to find the effect of the shape parameter on the center-of-mass diffusion of the star-shaped polymer chains in a solution.

The coarse-grained bead-spring model represents the polymer chains dissolved in the solution. To replicate good solvent conditions, the excluded volume interactions between the monomer beads are modeled using the repulsive part of the 12-6 Lennard-Jones (LJ) potential, also known as Weeks-Chandler-Andersen potential¹⁶ ($U_{WCA}$), defined as:

where $\sigma_{p}$ is the diameter of a monomer bead, $r$ is the distance between two beads and $\varepsilon=k_{B}T$ is the strength of interaction, $k_{B}$ is the Boltzmann’s constant and $T$ is temperature. The neighboring monomers of the polymer chain are connected with springs, the potential of which is given by Finitely Extensible Nonlinear Elastic (FENE)¹⁷, defined as:

where $k$ is the spring constant, and $r_{0}$ is the maximum length of the extension. The values of $\kappa$ and $r_{0}$ are 30 $\frac{k_{B}T}{\sigma_{p}^{2}}$ and $1.5\sigma_{p}$ respectively, as recommended by Kremer and Grest¹⁸. The bead spring model, FENE potential, and Kreme and Grest parameters are widely used in the coarse-grained modeling of the polymer chains¹⁹. The star-shaped polymer chains of varying degrees of functionality (number of arms) have been modeled by connecting different linear arms at their ends instead of connecting them to a single central monomer, ensuring equal flexibility of the arms for all the functionalities.

The solvent is modeled explicitly as an ensemble of non-interacting point particles of finite mass ($m$) using a mesoscopic coarse-grained simulation technique, MPCD²⁰. MPCD consists of alternating streaming and collision steps. In the streaming step, the MPCD particles with velocity $\boldsymbol{v}_{i}$ undergo ballistic motion and their positions ($\boldsymbol{r}_{i}$) are updated as:

In the collision step, the simulation box is divided into cubic cells of equal size($a$), and all the particles within a cell undergo stochastic collision. The collision of the MPCD particles is modeled using a momentum-conserving version of the Andersen Thermostat, also known as MPCD-AT²¹, in which the particle velocities ($\boldsymbol{v}_{i}$) are updated as:

where $\textbf{v}_{cm}$ is the center-of-mass velocity of the collision cell, $\boldsymbol{v}_{i}^{ran}$ is a random velocity selected from a Maxwell-Boltzmann distribution, and $\Delta\textbf{v}_{cm}^{ran}$ is the change in center-of-mass velocity of the collision cell due to the addition of $\boldsymbol{v}_{i}^{ran}$. During the streaming interval of MPCD, the positions, and velocities of the monomer beads evolve by the velocity-Verlet algorithm²² with a time step $\delta t_{MD}$. During the collision step, the monomers are considered MPCD particles and undergo stochastic collisions. The three components of $\boldsymbol{v}_{i}^{ran}$ are selected from a Gaussian distribution with variance $k_{B}T/m$ for the solvent particles and $k_{B}T/M$ for the monomer beads, where $M$ is the mass of a monomer. This way of considering the monomers just like other MPCD particles in the collision step for modeling solvent-monomer interaction is often used in recent studies^1223242526 due to its advantage of avoiding spurious depletion forces²⁷ which could lead to breakage of FENE bonds. Galilean invariance is ensured by randomly shifting the cells before each collision step by a vector with the three components randomly chosen from $[-a/2,a/2]$²⁸. All the simulations have been performed using the MPCD-AT routines in LAMMPS²⁹(Chen et al.^{26 30}).

The size of the collision cells is taken to be the same as the size of the monomer beads, $a=\sigma_{p}$. The average density of the MPCD solvent equals 5$\frac{m}{\sigma_{p}^{3}}$. The mass of a monomer($M$) is taken as 5$m$ to achieve neutral buoyancy. The MD time step ($\delta t_{MD}$) equals $0.002\tau$. The MPCD collision time step ($\delta t$) is $0.09\tau$, where $\tau$ is the intrinsic unit of time equals $\sqrt{m\sigma_{p}^{2}/k_{B}T}$. The resulting viscosity and Schmidt number($Sc$) of the MPCD fluid are $4\frac{k_{B}T}{\sigma_{p}^{3}}$ and 12, respectively. The size of the cubic simulation box is increased with polymer chain length to avoid the finite size effects following the previous studies.¹²²⁵ Box size($L$) equals $32\sigma_{p}$ for the set of chain lengths $\{24\sigma_{p},36\sigma_{p},48\sigma_{p}\}$ , $48\sigma_{p}$ for the set of chain lengths $\{60\sigma_{p},84\sigma_{p}\}$, and $64\sigma_{p}$ for the set of chain lengths $\{108\sigma_{p},192\sigma_{p}\}$. The equilibration simulation run is performed for $2\times 10^{6}$ MD time steps. The results are time averaged over $5\times 10^{8}$ MD time steps and ensemble-averaged over five system replicas, each with a unique set of random velocities at starting of the simulation and during the stochastic collision, both taken from Maxwell-Boltzmann distribution. The measured parameters will be expressed in reduced units using the energy scale $k_{B}T$, length scale $\sigma_{p}$, and mass scale $m$. Periodic boundary conditions are implemented in all directions. The snapshots of the simulations are shown in Figure 1.

To validate the MPCD routines, the Brownian motion of 250 colloidal particles of the same size as the monomer beads are modeled in the MPCD solvent. The variation of their average mean square displacement (MSD) with lag time ($\Delta t$) is plotted in Figure 2(a). Typically, power law describes the variation of MSD with lag time: MSD $\propto\Delta t^{b}$. The dynamics of the solutes can be diffusive ($b=1$), sub-diffusive ($b<1$), or super-diffusive ($b>1$). From Figure 2(a), we note that we obtain b = 1 as expected for the colloids. Hence, the dynamics of the colloids are captured well by the simulation. Further, the radius of gyration($R_{g}$) vs. chain length ($N$) is plotted for the linear chain in Figure 2(b). A power law behavior can be observed with a scaling exponent of 0.623 for the linear chain. The value of the power law exponent reported by Chen et al.²⁶ using similar simulation parameters is 0.61. Linear chains have been studied widely, and their corresponding scaling exponent values for good solvent conditions are summarized in Table 1. The exponent value calculated in the present work agrees well with earlier studies. In addition, the diffusion coefficient($D$) is calculated from the center-of-mass mean square displacement(MSD) vs. lag time plot for the linear chains of different lengths. The variation of $D$ with $N$ is shown in Figure 2(b). It also follows a power law $D\sim N^{-\nu_{d}}$, where $\nu_{d}=0.622$. The equality of $\nu$ with $\nu_{d}$ confirms the Zimm theory for the diffusion of a polymer chain with intra-chain hydrodynamic interactions, which predicts $D\sim\frac{1}{R_{g}}$.

The difference between a good solvent and a poor solvent can be demonstrated by introducing attraction using the standard 12-6 LJ potential with a cutoff of $2.5\sigma_{p}$ for the pairwise interaction between monomer beads instead of WCA potential as explained by Peng et al.³⁸. In a good solvent, the polymer chain forms a coil, whereas it becomes a globule in the case of a poor solvent. This coil-to-globule transition is shown in Figure 3 and can be observed by measuring the radius of gyration, which is approximately $5\sigma_{p}$ in a good solvent for the chain length of $56\sigma_{p}$, contrary to $2\sigma_{p}$ for the poor solvent case. This reduction in the size of the polymer chain is also visible in the values of the resulting diffusion coefficients, which are $0.0061\sigma_{p}^{2}/\tau$ and $0.0032\sigma_{p}^{2}/\tau$ for poor and good solvents, respectively.

In this work, the Brownian diffusion of the linear and star-shaped polymer chains of different functionalities is simulated using MPCD. It is shown that the radius of gyration of the star-shaped polymer chains follows a functionality-independent scaling law with chain length, in which the scaling exponent $\nu\sim 0.627$. The linear chain is shown to be more anisotropic than the star-shaped chains, and for star-shaped chains, the value of relative shape anisotropy decreases with an increase in functionality. For the same radius of gyration, the linear chain diffuses at a faster rate along the translational direction and has a slower rate of rotational diffusion than the star-shaped chains. Among star-shaped chains with the same radius of gyration, higher functionality leads to a higher value of rotational diffusion coefficient and a slower rate of diffusion along the translational direction. In terms of the shape parameter, we conclude that the star-shaped chains with higher values of relative shape anisotropy have a slower rate of rotational diffusion and therefore diffuse at a faster rate along the translational direction. Hence, shape anisotropy leads to faster center-of-mass diffusion of star-shaped polymer chains in a solution.

G.T. acknowledges partial support from the Department of Science and Technology National Supercomputing Mission HPC system in the Supercomputing Education and Research Center-Indian Institute of Science. A.K. acknowledges partial support from SERB CRG/2022/005381. P.K.P. acknowledges partial support from the Ministry of Education, Government of India.