Article pubs.acs.org/Macromolecules
Equilibrium Dynamics and Shear Rheology of Semiflexible Polymers in Solution Arash Nikoubashman*,† and Michael P. Howard‡ †
Institute of Physics, Johannes Gutenberg University Mainz, Staudingerweg 7, 55128 Mainz, Germany Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States
‡
ABSTRACT: We study the structure and dynamics of semidilute solutions of semiflexible polymers at rest and under shear using hybrid molecular dynamics simulations that take hydrodynamic interactions into account. We show that the polymer center-of-mass diffusion coefficient significantly decreases with increasing chain stiffness at fixed monomer density. The zeroshear viscosity shows a corresponding increase due to the intermolecular interactions of stiffer chains. We apply steady shear flow to the polymer solutions and show that at high shear rates the flow properties become almost independent of polymer stiffness. We characterize the polymer conformations under shear and find that in this regime polymers are elongated and aligned along the flow direction, but semiflexible polymers surprisingly form an increased number of hairpin-like conformations due to increased tumbling.
1. INTRODUCTION
In the past decades, several studies have focused on the conformational and rheological properties of single semiflexible polymers under flow using theory,25,26 simulations,27−31 and experiments.32−35 Here, shear alignment of the polymer and shear thinning of the solution were observed at large shear rates. Further, it was found that semiflexible polymers behaved increasingly like flexible polymers as the shear rate was increased. Previous theoretical considerations25,26 neglected excluded volume interactions, which is a reasonable approximation in the limit of infinite dilution and for sufficiently high stiffness. However, polymer concentrations are often close or above the overlap concentration ρ* (see section 3 for definition) in many practical situations, such as viscoelastic focusing in microfluidic devices36−38 or the interior of living cells.39 In these scenarios, excluded volume can play an important role, which makes a purely theoretical study of this problem considerably more difficult. In this work, we performed molecular dynamics (MD) simulations of a microscopic polymer model to study the equilibrium dynamics and shear rheology of semiflexible polymers in the semidilute concentration regime. MD simulations are an ideal tool for investigating the dynamic properties of such complex liquids because they provide a detailed picture on the nanoscale and allow for a systematic exploration of the effects of polymer stiffness and concentration on the solution properties. In particular, we studied a generic bead−spring polymer model with a tunable local bending stiffness spanning the flexible and semiflexible regimes at dilute and semidilute concentrations. Hydrodynamic interactions
Dilute and semidilute polymer solutions play an integral role for many technological applications, for example, as viscosity modifiers in the food and cosmetic industry1 or as additives in enhanced oil recovery.2,3 To date, extensive experimental,4−6 theoretical,7−9 and computational works10−14 have been conducted to study the rheological properties of these complex liquids and to identify general physical trends that are independent of the specific polymer chemistry. Of particular interest is how the microscopic properties of the macromolecules, such as the degree of polymerization or polymer topology, influence the macroscopic properties of the solution, such as the zero-shear viscosity. The majority of previous research has focused on fully flexible polymers, which are characterized, both in solution and in the melt, by a self-similar fractal structure from the scale of a monomer up to the coil size.7,15−17 This peculiar property allows for the development of elegant and simple scaling arguments for the dynamic properties of polymer solutions.7 However, many biological and synthetic macromolecules are not fully flexible and instead have a local bending rigidity, usually characterized in terms of their persistence length lp.16 If the persistence length is much greater than the chain contour length L, the polymer behaves like a rigid rod, whereas fully flexible chains have lp ≪ L. Semiflexible polymers occupy an intermediate regime of stiffness where lp ∼ L, which can be realized for e.g. double-stranded DNA (lp ≈ 50 nm)18 and actin filaments (lp ≈ 10 μm).19 For such semiflexible macromolecules, there are several crossover length scales relevant for describing the static structure,20,21 and it is even more challenging to characterize their dynamics due to the large number and disparity of the relevant time scales.22−24 © XXXX American Chemical Society
Received: August 31, 2017 Revised: September 29, 2017
A
DOI: 10.1021/acs.macromol.7b01876 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
stochastic collision steps. During a streaming step of time Δt, the solvent particle positions were simply updated according to
were incorporated into the model using the multiparticle collision dynamics (MPCD) technique.40−43 We show that the zero-shear viscosity of the polymer solution increases substantially as the persistence length is increased at fixed monomer density due to intermolecular excluded volume effects. For high shear rates, the flow properties of the polymer solutions become almost independent of polymer stiffness. In this regime, the polymers are increasingly elongated and aligned along the flow direction, and at the same time, an increasing number of the semiflexible polymers exhibit hairpin-like formations. The rest of this article is organized as follows. In section 2 we provide a detailed description of the employed model and simulation methods. In section 3 we present our results on the dynamic properties of semiflexible chains under quiescent conditions. In section 4.1 we study the rheological properties of polymer solutions under shear and then draw connections to the microscopic structural chain properties in section 4.2. We summarize our main findings in section 5.
ri(t + Δt ) = ri(t ) + vi(t )Δt
with ri and vi being the position and velocity of solvent particle i, respectively. During this interval, the solute particle positions and velocities evolved according to Newton’s equations of motion, which were integrated using the velocity-Verlet algorithm47 with time step ΔtMD. In the collision step, both solute and solvent particles were sorted into cubic cells with edge length a. The cells were shifted by a random three-dimensional vector with components drawn uniformly on [−a/2, +a/2] to ensure Galilean invariance.48 The particle velocities were then updated through an Andersen thermostat collision scheme42 vi(t + Δt ) = u(t ) + Δv iran − Δu ran
(1)
τsim = ma 2 /(kBT ) is the intrinsic unit of time. In the rest of this article, we will report all simulation parameters and results in implicitly reduced units using the energy scale kBT, the bead diameter a, and the solvent mass m, unless explicitly stated otherwise. All simulations were performed in a cubic simulation box with edge length Lbox = 50, and periodic boundary conditions were applied in all directions. Note that we chose Lbox > L to prevent unphysical self-interactions in the case of a fully extended polymer. Initial configurations were generated by randomly placing the chains into the simulation box at monomer number densities 0.025 ≤ ρ ≤ 0.1. Starting configurations were equilibrated by running the simulations until the radius of gyration, Rg (see eq 15), reached a plateau for a given κ. During this step, we turned off the MPCD algorithm to expedite the simulations and used an Andersen thermostat to maintain constant temperature.49 We then simulated the systems for t = 2 × 105 to 5 × 105 in order to sample the systems at rest and under shear. Shear was applied using the reverse perturbation method.50,51 This approach differs from other nonequilibrium methods for generating shear flow, such as Lees−Edwards boundary conditions52 and the SLLOD algorithm,53 because the external stress σ is imposed by generating a momentum flux, and the flow field emerges rather than being explicitly prescribed. In practice, the momentum flux is created by swapping solvent particle velocities in the following way. First, the periodic simulation box is subdivided into equally sized slabs with thickness a along the gradient direction (x) of the flow (z). Then, the particle i in the x = 0 slab with the largest positive z
where ε = kBT controls the strength of the interaction and r is the distance between beads (kB is Boltzmann’s constant and T is the temperature). Bonds were represented through the finitely extensible nonlinear elastic (FENE) potential45 ⎧ ⎡ ⎛ r ⎞2 ⎤ 1 2 ⎢ ⎪ ⎪− kr0 ln 1 − ⎜ ⎟ ⎥ , r ≤ r0 ⎢⎣ UFENE(r ) = ⎨ 2 ⎝ r0 ⎠ ⎥⎦ ⎪ ⎪∞ , r > r0 ⎩
(2)
We employed the standard Kremer−Grest parameters46 for the spring constant (k = 30kBT/a2) and maximum spring extension (r0 = 1.5a). With these parameters, the average bond length was b = 0.97a, which precluded unphysical bond crossing. The contour length of a polymer was then L = (N − 1)b ≈ 45.6a. Local bending rigidity was incorporated into the model using Ubend(Θijk ) = κ(1 − cos Θijk )
(5)
where u is the center-of-mass (CM) velocity of the cell to which particle i belongs, and Δvran is a random velocity with i components drawn from a Gaussian distribution with variance kBT /m for the solvent particles and kBT /M for the solute particles. The last term Δuran is the change in the cell CM velocity from adding Δvran i , and subtracting it enforces local (and hence also global) conservation of linear momentum. Note that the collision rule eq 5 acts as an implicit thermostat for the fluid. The dynamic properties of the solvent are controlled by the solvent density ρs and the time step Δt. The solvent density was set to ρs = 5a−3, and the mass of the polymer beads was chosen as M = 5m to match the average solvent mass in an MPCD collision cell. The time between MPCD collisions was set to Δt = 0.1τsim, and the MD time step was ΔtMD = 0.01τsim, where
2. MODEL AND SIMULATION METHODS We studied a generic semiflexible bead−spring polymer model. Each polymer chain consisted of N = 48 beads of diameter a and mass M linearly connected by springs. The excludedvolume interactions between the beads were modeled by the purely repulsive Weeks−Chandler−Andersen (WCA) potential44 ⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ a a 1/6 ⎪ ⎪ 4ε⎢⎝⎜ ⎠⎟ − ⎝⎜ ⎠⎟ ⎥ + ε , r ≤ 2 a r r ⎣ ⎦ ⎨ UWCA(r ) = ⎪ ⎪ 0, r > 21/6a ⎩
(4)
(3)
where Θijk is the angle between subsequent bonds connecting beads i, j, and k. (An angle of Θijk = 0° corresponds to three beads in a line.) The parameter κ controls the bending stiffness. The persistence length of the polymers is then defined as lp = −b/ln⟨cos Θijk⟩, with lp/b ≈ κ/(kBT) for κ/(kBT) ≳ 2.24 We employed a hybrid simulation approach that combines MD simulations for the polymers (solute) with the MPCD technique40−43 to account for hydrodynamic interactions mediated by the solvent. In MPCD, the solvent is modeled explicitly as point particles with mass m, and the particle motion is governed by alternating ballistic streaming and B
DOI: 10.1021/acs.macromol.7b01876 Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
observed a slight shrinking of the polymers with Rh = 3.79 ± 0.01 and Rg = 4.53 ± 0.05. At this density, the system is still slightly below the overlap concentration of the solution, ρ* ≡ N/(4πRg3/3) ≈ 0.11, which explains the small change in coil size. For the semiflexible polymers, both Rg and Rh increased with increasing persistence length, as expected. Interestingly, Rg exhibited a more pronounced increase compared to Rh as shown in Figure 1. Further, we did not find any appreciable ρ
component of the velocity and the particle j in the x = Lx/2 slab with the largest negative z component of the velocity are identified, and their velocities are swapped. This swapping procedure artificially generates a momentum flux, which gives rise to a physical flow. If both particles have the same mass, as is the case in our model, swapping conserves both linear momentum and kinetic energy. The imposed shear stress can be controlled by the number and frequency of momentum swaps. For our geometry, the average shear stress can be computed as
⟨σ ⟩ =
⟨Δpz ⟩ 2ΔtL box 2
(6)
where ⟨Δpz⟩ is the average total z component of momentum exchanged during one time step. The shear rate can then be extracted from the measured flow profile in the linear response limit. In nonequilibrium MPCD simulations, it is recommended to keep the solvent velocities well below the speed of sound in the medium cs to avoid artificial compressibility effects.54,55 Therefore, we restricted our simulations to Mach numbers Ma = vmax/cs < 0.5, with cs = 1 for an MPCD liquid under isothermal conditions. Further, note that the velocity profile generated by this method has a triangular shape with maxima in the exchange regions (see Figure 2 in ref 50), leading to a change in the direction of shear there. Although this nonuniformity in the shear field is negligible for small particles, it can potentially cause complications for objects that span multiple streamlines in those exchange regions. In fact, we found that the polymer conformations were nonuniform along the gradient direction of flow in the vicinity of the exchange regions. To exclude such finite size effects from our simulations, we excluded all polymers with their centers of mass within 5a of these exchange regions from our data analysis.
Figure 1. Radius of gyration, Rg, and hydrodynamic radius, Rh, of a semiflexible chain in ultradilute solution as a function of κ. The dashed line indicates the ratio Rh/Rg.
dependence of Rg, Rh, and lp for κ > 2 in the investigated density regime. Note that the overlap concentration ρ* decreases substantially with increasing persistence length due to the accompanied increase of Rg and is, for example, only ρ* ≈ 6.0 × 10−3 for κ = 50. To determine D via eq 7, we also need to know the dynamic viscosity of the solution, η0, as well as the monomer diffusion coefficient, D0. For sufficiently low polymer concentrations ρ ≪ ρ*, the solution viscosity can be approximated by the solvent viscosity, ηs, for which an analytic expression has been derived.42,59 For the employed MPCD parameters, we find ηs = 3.71, which is in excellent agreement with the value we measured via shear simulations of the pure solvent. An analytic expression has also been developed for the selfdiffusion coefficient D0,s of an MPCD solvent particle under the molecular-chaos assumption.59,60 However, it has been shown in ref 59 that this expression underpredicts the actual diffusion coefficient when the solvent mean free path is too short. It was argued that this mismatch stems from hydrodynamic backflow, which is not captured in the analytic expression for D0,s. We tested the accuracy of the prediction by computing D0,s from the mean-squared displacement of the solvent particles in a pure MPCD liquid and found D0,s = 0.093 ± 0.001, which is indeed significantly larger than the predicted value of D0,s = 0.075. When hydrodynamics were disabled by randomly swapping solvent velocities after each collision step,61,62 we measured a value of D0,s = 0.075 ± 0.001 for the solvent particles, which is in excellent agreement with the theoretical prediction. Hence, the diffusion coefficient of the monomers cannot be calculated from theory but needs to be determined from simulation. We measured D0 for free monomers at all investigated densities ρ and found D0 = 0.051 for the ultradilute case and D0 = 0.046 for ρ = 0.1. The measurement uncertainty was
