Debye Screening Effect on Scaling Behavior of Longest Relaxation Time of Biological Polyelectrolyte Chain

Abstract

The scaling relationship of the longest relaxation time of a single chain of semiflexible biological polyelectrolyte has been investigated by performing well-established coarse-grained Brownian dynamics simulations. Two kinds of longest relaxation times were estimated from time-sequences of chain trajectories, and their behaviors were interpreted by applying the scaling law for different molecular weights of polyelectrolyte and Debye lengths. The scaling exponents for longest stress relaxation and rotational relaxation are found in the ranges of 1.67-1.79 and 1.65-1.81, respectively, depending on the physicochemical interaction of electrostatic Debye screening. The scaling exponent increases with decreasing screening effect, which is a special feature of polyelectrolytes differing from neutral polymers. It revealed that the weak screening allows a polyelectrolyte chain to follow the behavior in good solvent due to the strong electrostatic repulsion between beads.

Introduction

In the study of water-soluble polyelectrolyte molecules, as opposed to neutral polymers, some special aspects should be considered. This is because these have the ability to dissociate charges in poor solvents, resulting in charged macroions chains and mobile counter-ions. Biological polyelectrolytes represent an interesting and broad class of soft matter that should be significantly treated, as an additional complexity arising from the Debye screening by physicochemical long-range electrostatic interaction.12 A precise understanding of charged soft matters has become of great interest, since proteins and other biomolecules, such as DNA, polysaccharides, and microtubules, are generally polyelectrolytes. Analyzing their conformational dynamics has become important for several applications, for example, DNA sequencing, microfluidic separations, and diagnosis of diseases.3-9

In polymer physics, the relaxation time of a single molecule represents how quickly it recovers from an applied deformation or stress. Its longest relaxation time is the essential piece of information, offering the natural time constant of the molecule and characterizing how fast a polymer can react in response to an imposed field.1011 Moreover, relaxation information plays an important role in optimum design (e.g., channel dimension and geometry) of microfluidic devices. For instance, the relaxation time of a DNA molecule can qualitatively let us know its conformational change with time when it is stretched in a microchannel to detect its sequence. Bakajin et al.12 examined the relaxation behavior of a T4-DNA in a slit-like channel after it disengaged itself from the obstacles inside a channel. Hsieh et al.13 conducted on the DNA relaxation in nano-slits, establishing a scaling model of rotational relaxation time in terms of channel configuration and molecular weight. Note that there are several types of longest relaxation times for single polymer molecules in dilute solutions, e.g., the longest rotational relaxation time τl,r and the longest stress relaxation time τl,s.1415 For a polymer chain in either theta or good solvent, both Rouse and Zimm models present the fact that τl,s is a half of τl,r.15

The scaling behavior of single polymer molecules with variations of chain properties, such as molecular weight (Mω), is meaningful because polymer solutions can be characterized with the scaling exponent. For instance, the relaxation time is proportional to Mω 2 in the Rouse theory, which is satisfactory for a free-draining Gaussian chain. It is proportional to Mω 3/2 in the Zimm theory, which is adopted for a chain in a theta solvent when the hydrodynamic interaction (HI) is dominant.16 The scaling behavior of neutral polymer chains with a self-avoiding walk has been extensively explored for their structural or relaxation properties. On the contrary, the scaling behavior for polyelectrolyte chains has rarely been reported, in which screening effect should be taken into account.

Various simulation methods have been exploited to study the single chain dynamics. The Brownian dynamics (BD) simulation is a promising method that delivers a more efficient and plausible modeling technique at less computational cost in comparison with others. The most fascinating gist of the BD is perhaps its capability to cope with high-frequency stochastic forces due to a unique characteristic in the time-scale separation. BD simulations for polymer chains have been continually improved to reflect more realistic circumstances. First, BD simulations with a free-draining model were employed to examine DNA molecules as regards their dynamic behavior in shear flows and their transient response under start-up conditions.1718 Coarse-grained models for a polyelectrolyte chain in flows as well as bulk system have also been investigated, incorporating the HI between pairs of beads. The implementation of the HI in BD simulations has reinforced the explicit understanding on conformation and migration of a single DNA molecule in several confined channels and in more complex flows.819-22

The purpose of this study is to accurately predict the longest relaxation time of a single chain of semiflexible polyelectrolyte with submicron size and to clarify a difference to the neutral polymer in view of its scaling behavior. As a biological model polyelectrolyte, we apply the polysaccharide xanthan, which has the double-stranded structure. It has two ionizable carboxyl groups per repeat unit with the side chains imparting a semi-rigid character to the main chain. In this study, its longest rotational and stress relaxation times with relevant scaling have been examined by correlating the Debye screening effect and the molecular weight. Based on the scaling exponent, relaxation property of polyelectrolyte chain can usefully be predicted with various molecular weights.

 

Mesoscale Brownian Dynamics Simulations

For BD simulations, we followed the coarse-grained beadspring model of polyelectrolyte xanthan as well-presented in our previous studies.23-26 Coarse-graining neglects all detailed atomistic interactions, which are replaced by effective interactions. 27 The polyelectrolyte chain is considered to be divided into beads that are connected by springs, each of which fluctuates in length due to thermal excitations. Basic model parameters were determined from experimental rheology data on the purchased xanthan, which had Mω of 1.13 × 106 g/mol (cf. 1,220 monomers). Considering the ratio of contour length (RC) to persistence length (lp), xanthan can be modeled as a mesoscopic discrete wormlike chain with Nb beads connected by Nb ˗ 1 springs. For Nb = 25, a smooth change of lp can be ensured by taking a Kuhn segment length (= 2lp) of about 10 beads.

The dynamics of the polyelectrolyte chain is delineated by the following stochastic equation28 that commonly accounts for the fluctuating bead-bead HI

Here, is the position of the i-th bead, is the total force on the j-th bead at n-th time tn, Etotal is the total potential, is the random displacement due to the solvent, Δt is the time step, and kBT is the Boltzmann thermal energy. Brownian forces are coupled to the velocity perturbations through the fluctuation-dissipation theorem and can be approximated by the White noise process. The Gaussian distribution of yields a quantity with mean (due to the isotropic nature of the collisions) and covariance .

The velocity field generated from the segmental motion of the xanthan molecule is taken into account by implementing the bead-bead HI, which is represented by a chain of point forces acting on the solvent. In the manner of recent studies, the Rotne-Prager diffusion tensor was employed. This is positive-definite for all chain configurations and is expressed as with formulating 3 × 3 block components of 3Nb × 3Nb matrix. Here, η is the solvent viscosity, a is the bead hydrodynamic radius, I is the unit tensor, and Ωij is the HI tensor. HI relates to a velocity perturbation of a bead carried by the surrounding fluid at point ri to a force at point rj.

Total potential Etotal needs to include the spring, bending, Lennard-Jones (LJ) potential, and electrostatic interaction. We describe these briefly here, since details were provided in previous reports.2326 Together with the concept of LJ particles, the xanthan molecule is empirically represented by the finitely extensible nonlinear elastic (FENE) spring model

where l0 and lmax indicate the equilibrium and the maximum allowable bond lengths, respectively. In the limit of large lmax, EFENE reduces to a Hookean potential with force constant kS.

The semiflexible chain is distinguished by the model with the bending potential, compared to the flexible one. The finite rigidity of the chain is modeled by the harmonic bending potential , where kA is the bending force constant, the angle , and is the i-th unit bond vector. The LJ potential describing the dispersion-repulsion interaction between pairs of beads is expressed by introducing the energy parameter (kBT) and the length scale (l0).2930 Incorporating the excluded volume effect, the electrostatic interaction between charged beads is addressed through the screening Coulomb interaction, known as the Debye-Hückel potential, . Here, qb denotes the bead charge, κ is the inverse Debye screening length, and the medium dielectric constant ε is defined in terms of the relative permittivity (= 78.5 for water) and the vacuum permittivity (= 8.854 × 10˗12 C2/J·m).

The spring parameters (e.g., ks, l0, lmax, kA) were previously determined from preliminary BD simulations.2324 Dimensional properties of the chain were also verified by experimental values of RC and lp at different screening lengths κ˗1. For a semiflexible chain, the optimum stretching constant was chosen along with l0 = 10 nm, and lmax = 150 nm. The optimal bead charge qb of ˗35e for the xanthan chain corresponds to ~35% of the maximum value (i.e., ˗98e) for a bead composed of 49 monomers, leading to the contour distance per elementary charge lq (= RC/(Nbzb)) of 0.66 nm. The optimized kA determined as 4kBT yields lp = 4b, where b satisfies the relation lp = bkA/(kBT). We set the time step Δt = 2 ns, satisfying m0D0/kBT << Δt << a2/D0 (cf., D0 is the self-diffusion of a single bead) for the bead mass m0. Here, the motions of interest occur on a time scale much longer than the momentum relaxation time. Input values for each parameter were already found to be satisfactorily suitable to verify the dynamics of single xanthan chain under relevant conditions.2426

 

Longest Relaxation Times of Polyelectrolyte Chain

Additional algorithm to calculate the longest relaxation time was implemented in the BD simulations. The longest stress relaxation time, τl,s, can be determined from the exponential decay rate of the stretched molecule.31 Here, the transition of the extended chain (cf., initial set by 90% of full extension) to its coiled state is recorded along the time. Then, τl,s is calculated by fitting the transient relaxation curve from t = t to t = ∞, . Here, xE is the stretching position along the x-axis at any time, and L is the fully extended length of the chain.

We admit the in-plane two-dimensional (2D) rotation from 3D simulations in view of proper experimental comparisons. The longest rotational relaxation time, τl,r, is determined from the exponential decay rate of the autocorrelation function for the orientation angle of a xanthan chain.13 Its orientation angle is related with the radius of gyration (Rg) tensor

where subscripts x and y represent in-plane coordinates. It offers structural information about the size, shape, and orientation of the molecule. Its each component is given as

Here, is the center of mass. Next, Rg is the square root of the trace of G, such that

As mentioned before, the longest rotational relaxation time, τl,r, is extracted from the rate of change of the xanthan orientation. The orientation angle θ between the principal eigenvector and the x-axis should be first calculated as

where λ(t) is the eigenvalue associated with the principal eigenvector.

Further, the time autocorrelation function Cr of the orientation angle is defined by

where < > denotes an ensemble average. The autocorrelation function describes the interrelationship between physical properties at any given time (t) and lagged time (t+δt). Its increase implies that the orientation angles are highly correlated between the specified time and the lagged time. Here, Cr of the xanthan chain decreases with increasing time lag. Considering its exponential decay rate, τl,r can be calculated from Eq. (7).

 

Results and Discussion

Time-sequence configurations of xanthan chain were calculated from the position vector of each bead. We obtained the average over multiple simulations for 40 trajectories of the xanthan molecule. Figure 1 displays the mean square stretch of the xanthan molecule versus the simulation time for different Nb from 10 to 40 corresponding to 4.5 × 105 g/ mol ≤ Mω ≤ 1.8 × 106 g/mol. Full range of the screening is considered as highest screening (κ˗1 = 0.97 nm, lp ≅ 120 nm, and RC ≅ 580 nm at bulk ionic concentration I = 100 mM) and lowest screening (κ˗1 = 970 nm, lp ≅ 530 nm, and RC ≅ 1.1 μm at I = 10˗4 mM). Here, the Debye screening length, is determined from Avogadro’s number NA, the elementary charge e, the medium ionic strength I, and the valence of i ion Λi. The xanthan molecule extended to 90% of its maximum length (i.e., the initial dimensionless mean square stretch was 0.81 with RC = 3.24 μm) was initially placed in no bulk flow regime. The dimensionless mean square stretch decreases as time elapse. The xanthan molecules with different Nb exhibit different exponential decay patterns of mean square stretch , as shown in Figure 1. As the molecular weight (i.e., Nb) of the xanthan molecule increases, the decay rate gets slower caused by the increase of molecule size. Thus, it takes a longer time for a large molecule to relax.

Figure 1.Polyelectrolyte chain relaxation as a function of time for different number of beads. (a) highest screening (κ ˗1 = 0.97 nm) and (b) lowest screening (κ ˗1 = 970 nm).

In Figure 1, when the simulation time elapsed sufficiently, the mean square stretch of xanthan molecules with different Nb values approached their own equilibriums. The final square stretch value was dependent on the Debye length of the xanthan molecule, such that /L2 ~ 0.012 and RC = 575 nm in highest screening and /L2 ~ 0.10 and RC = 1.08 μm in lowest screening. This means that the conformation of the xanthan chain is significantly influenced by the screening effect according to the change in the medium ionic concentration. For the weak screening case with a lower salt concentration (i.e., a larger Debye length), the equilibrium mean square stretch of the xanthan molecule becomes larger, which is attributed to the increased repulsive force between beads. Figure 2 evidently portrays the different evolution of chain conformation according to the passage of time for Debye length κ˗1 of 0.97 nm (left) and 970 nm (right) with Nb = 25. In the equilibrium state (at 10 ms), the final stretch in the case of small κ˗1 is less than that in the case of large κ˗1.

The longest stress relaxation times with different Nb and screening effects were calculated by curve fitting of the mean square stretch (Fig. 3). In the case of small κ˗1, repulsive electrostatic interaction between beads is weak, making the transit rate from the extended state to the coiled one slow and thus the relaxation time of xanthan increases. In order to ascertain their scaling behaviors, we try to plot the longest relaxation time against the molecular weight (i.e., Nb). It is found that a xanthan with higher Nb has a larger relaxation time because of the larger molecule size. The slope exponent α obtained for τl,s ~ from the best fit of simulation results is in the range 1.67-1.79, depending on the screening effect. The slope increases with increasing κ˗1. Note that the exponent α is 2.0 for the Rouse theory, 1.5 for the theta solvent of Zimm theory, and 1.8 for the good solvent of Zimm theory (self-avoiding walk).16 Considering the HI in xanthan polymer, the slope of the chains can be applicable to the regime of Zimm theory as the screening becomes weaker. Moreover, the property of xanthan solution follows that of the good solvent in the Zimm theory, revealing that the excluded volume effect of beads becomes stronger due to the strong electrostatic repulsion. Good solvents imply that polymer-solvent contacts are favorable, so the polymer chain expands to contact as many solvent molecules as possible. For theta solvents, polymer chains collapse oppositely due to unfavorable polymer-solvent contacts.

Figure 2.Typical time-sequence snapshots of polyelectrolyte xanthan chain under different screening conditions: The projection is chosen to display the longest axis in the paper plane.

Figure 3.The longest stress relaxation time versus molecular weight (Mω) of xanthan chain with different screening effects.

Unlike the longest stress relaxation time, the longest rotational relaxation time, which indicates how fast its orientation change, can be evaluated from the autocorrelation function for the orientation angle of xanthan molecule (cf., Eq. (7)). Figure 4 shows the temporal changes in the orientation angle of the molecule with Nb = 25 at lowest screening. Owing to the Brownian motion of beads within a single molecule, albeit its equilibrium state, the orientation angles unceasingly fluctuate with time progress. Even in the same simulation case, the orientation angle of the molecule determined from the gyration tensor shows different trajectories according to their Brownian forces. Figure 5 shows the average gyration tensor components (Gxx, Gxy, Gyy) along the orientation angle provided in Figure 4. When the molecule is placed on the x-axis (i.e., θ = 0), the averaged Gxx and Gyy become maximum and minimum, respectively. On the other hand, when the molecule lies on the y-axis (i.e., θ = ˗π/2 or π/2), the averaged values of Gxx and Gyy change to opposite. Also, the absolute value of Gxy is the highest at θ = ˗π/4 or π/4. It is evident that patterns of gyration tensor components along the orientation angle would be physically unchanged under different screening conditions.

Figure 4.Orientation angle of xanthan chain along the time for Nb = 25 and κ˗1 = 970 nm.

Figure 5.Average gyration tensor components (Gxx, Gyy, Gxy) of the xanthan molecule along the orientation angle.

Each orientation angle of the xanthan chain with different screening effects along the elapsed time is depicted in Figure 6. Strong screening (κ˗1 = 0.97 nm, 9.7 nm) results in highly fluctuations in the orientation angle, compared to the case of weak screening (κ˗1 = 97 nm, 970 nm). As exhibited in Figure 7, this result can be represented as the temporal configuration of the single molecule, illustrating the conformation changes with screening effect. Molecules remain in the coiled state under strong screening, however, they evolve transition to the extended one under weak screening. When the xanthan molecule is in the coiled state, little movement is necessary to change the orientation angle. However, in the extended state, a large movement should be required to change the orientation angle. Based on this fact, the angular change in the coiled state is higher, resulting in a larger relaxation time.

Figure 6.Screening effect on the orientation angle of xanthan chain.

Figure 7.Temporal images of xanthan chain with different screening effects.

Figure 8.Autocorrelation function for the orientation angle of xanthan chain with Nb = 25.

Figure 8 shows the autocorrelation function with the lag time in the case that Nb is 25 and Debye lengths are 0.97 and 970 nm. From the exponentially decay rate, the longest rotational relaxation time can be calculated by Eq. (7). The dashed line is the single-exponential best fit that extracts the longest rotational relaxation time from the autocorrelation function. Figure 9 depicts the longest rotational relaxation times of xanthan molecules with several Debye lengths along Nb. The scaling pattern of the longest rotational relaxation time is actually analogous to that of the longest stress relaxation time. In Figure 9, the slope exponent through the best fit of the simulation data ranges from 1.65 to 1.81, demonstrating the doubled relaxation time in comparison with the longest stress relaxation time under the similar slope range. These results correctly follow the aforementioned theoretical expectation on relaxation times.

Figure 9.The longest rotational relaxation time versus molecular weight (Mω) of xanthan chain with different screening effects.

Based on the scaling behavior of the longest stress and rotational relaxation times of xanthan chain, the difference between water-soluble biological polyelectrolytes and neutral polymers can be clarified. For the polyelectrolyte chain, it should be emphasized that the electrostatic screening has a substantial effect on the solvent property. Since the HI force between beads was considered, the range of the scaling exponent exists within the regime of the Zimm theory. It is also worthwhile mentioning that the scaling behavior of the longest relaxation time of xanthan chain with screening effect is similar to the case of radius of gyration but opposed to the case of diffusivity.2426

 

Conclusion

In this study, the scaling behavior of the longest relaxation time in a single chain of semiflexible polyelectrolyte was examined by employing coarse-grained BD simulations on the model biopolymer xanthan. First, the longest stress relaxation time was determined by calculating transient configurations of the xanthan molecule: i.e., the transition rate from extended to equilibrium states. It implies how quickly the chain recovers from its extended state after the external deformation. Second, the longest rotational relaxation time was evaluated from the exponentially decay rate of the autocorrelation function regarding the orientation angle of the xanthan molecule. Dynamic properties of the polyelectrolyte chain and screening effects were elucidated, by analyzing the scaling plot of the relaxation time with respect to the molecular weight. The scaling exponents of τl,s and τl,r were in the ranges of 1.67-1.79 and 1.65-1.81, respectively, according to a change of the Debye screening. It turns out that the scaling exponent increases as the Debye screening becomes lower. Furthermore, the strong repulsive electrostatic interaction between the beads allows a polyelectroyte chain to match apparently the behavior in good solvent.