Statistical Behaviors of Semiflexible Polymer Chains Stretched in Rectangular Tubes

We investigated the statistical behaviors of semiflexible polymer chains, which were simultaneously subjected to force stretching and rectangular tube confinement. Based on the wormlike chain model and Odijk deflection theory, we derived a new deflection length, by using which new compact formulas were obtained for the confinement free energy and force–confinement–extension relations. These newly derived formulas were justified by numerical solutions of the eigenvalue problem associated with the Fokker–Planck governing equation and extensive Brownian dynamics simulations based on the so-called generalized bead-rod (GBR) model. We found that, compared to classical deflection theory, these new formulas were valid for a much more extended range of the confinement size/persistence length ratio and had no adjustable fitting parameters for sufficiently long semiflexible chains in the whole deflection regime.


Introduction
Statistical physics properties of single polymer chains can be significantly influenced or even determined by geometrical confinements and applied external forces [1][2][3][4]. A detailed understanding of the behaviors of polymers under such circumstances is still considered to be an unsolved problem in polymer physics after more than half a century. However, even so, advances in the study of geometrically and potentially constrained polymers are playing a role in the development of many existing nanotechnologies of genomics and materials science, etc. [5][6][7].
For polymers under confinement, the effects of constraints have usually been classified into three regimes (weak, moderate, and strong confinements), which are distinguished in terms of the comparison between the polymer's unconfined radius of gyration and the Kuhn length of the typical confinement length scale. In the regime of weak confinement, Casassa [8] has discussed the free energy of ideal chains trapped in pores with different shapes based on the theory of diffusion. Then, de Gennes and his coworkers [9,10] developed the so-called blob model to describe the confinement of a long, flexible polymer. According to this model, the free energy of confinement is given by Here, H is the length scale of confinement, k B is the Boltzmann constant, and T is the absolute temperature. R g ≈ L a N v is the radius of gyration of the unconfined polymer, N is the number of monomers, L a is of the order of the monomer-monomer separation, and v equals 1/2 for ideal polymers and approximately 3/5 for polymers with excluded volume. A widely used model of a polymer with bending energy is the wormlike chain (WLC) model, characterized by the contour length L and persistence length L p , which was first proposed by Kratky and Porod in 1949 [11]. In the strong confinement regime, Odijk [12,13] has argued that polymer behavior can be interpreted in terms of L/λ statistically independent segments with deflection length λ. Based on this picture, he [12][13][14] obtained expressions of the confinement free energy, F, and the average extension of the chain, R || , in terms of λ as where λ fe , λ ext represent the free energy and extension-associated deflection lengths, respectively: λ fe , λ ext ∝ L 1/3 p D 2/3 [12] has been suggested for the confinement of a cylindrical tube with diameter D. For the confinement of a rectangular tube with height and width H h and H w , the deflection length associated with the free energy calculation has been suggested as [15,16] In contrast, this deflection length associated with the average extension is given as We note that Equations (4) and (5) have different expressions and prefactors. The prefactors in Equations (4) and (5) were determined by using various numerical techniques and theoretical derivations, such as the Monte Carlo simulations [17,18] and eigenvalue technique associated with the Fokker-Planck equations [15,19]. Examples of the determined prefactors are illustrated in Table 1. We can see that the prefactors 1/A and α , respectively determined from free energy and extension, have an almost 10 times difference in quantity.  [18] 0.09137 ± 0.00007 [18] 1.1032 ± 0.0001 [19] 0.09143 ± 0.0001 [19] In addition, a slit of separation H can be regarded as a rectangular tube with height H h = H and infinite width H w → ∞ . Statistical properties of polymer chains confined in the slit have been extensively studied [20][21][22][23][24][25] based on Monte Carlo simulations and eigenvalue analysis. The deflection length in a strong confinement regime has been confirmed to follow the Odijk scaling law, Although the deflection length in Equation (5) can be viewed as the combination of that for two slits with heights H h and H w , respectively [15], it can be observed from Equations (4)-(6) that the deflection length in Equation (5) is not consistent with Equation (6) as H w → ∞ .
Beyond the Odijk regime, Chen [24] numerically calculated the confinement free energy by treating the problem of a confined polymer as an eigenvalue problem. He also suggested an interpolating formula that can have very good agreement with that of the numerical calculations for polymers under both strong and weak confinements. In addition, an extended de Gennes regime [20,26,27] has also been identified based on the Monte Carlo simulations.
Interestingly, external forces can pose similar effects to the statistical behaviors of single polymer chains as geometrical confinements. For a polymer chain to be stretched by a sufficiently large force, f s , a deflection length also exists and can be expressed as [3,28] λ f = L p / f , wheref = f s L p /k B T, so that the force-extension relation can be expressed as Polymers in real microenvironments are usually subjected to both geometrical constraints and external forces. Wang and Gao [29] have revealed that the average extension of a strongly tube-confined and force-stretched polymer chain can be equivalent to that of an unconfined chain subjected to an effective stretching force. Li and Wang [30] later confirmed that this equivalence property is still valid for the tube-confined polymers in a much more extended Odijk regime. Therefore, for a semiflexible polymer chain in the deflection confinement regime, one can generally have where one can setf c = L 2 p /λ 2 fe orf c = L 2 p /λ 2 ext as the normalized effective force due to the existence of strong confinement: For the confinements of rectangular tubes, λ fe , λ ext are given by Equations (4) and (5).
However, as shown above, Odijk deflection lengths based on free energy and extension are different. Then a critical question arises. Which deflection length should be used if the polymer chain is under both geometrical confinement and force stretching? Therefore, there are still open questions on how the Odijk length can be uniquely and precisely defined for polymer chains confined in rectangular tubes.
In this study, for semiflexible polymer chains confined in rectangular tubes and slits, we derived a modified deflection length, which was expected to be valid for a more extended range than the classical Odijk deflection length. This extended deflection length was directly used to quantitatively formulate both energy and extension. We will present numerical calculations based on the eigenvalue technique developed by Chen and coworkers [19,31,32] and Brownian dynamics simulations in terms of the generalized bead-rod (GBR) model [30,33,34] to justify our theoretical predictions.

Model
We first considered a WLC confined in a rectangular tube with width H w and height H h , as shown in Figure 1. We assumed that the tube was small so that the chain's configurations with the so-called "hairpin" structures [14,35] rarely existed. This means statistical behaviors of the chain fell into the deflection regime. In order to obtain a universal deflection length scale, λ, to simultaneously characterize both the free energy and extension of the chain in terms of the ideas of de Gennes [9,10] and Odijk [12], the chain was assumed to behave like L/λ independent free segments, aligning one by one along the tube axis, so that the conformational free energy could still be scaled as Equation (1) and the average extension could be simply the sum of that for each free segment. On the other hand, when considering the average extension of the chain, the quantitative behavior of each segment should be in analogy with a free chain of effective contour length c 1 λ m , in which the parameter c 1 actually reflects the influence of two artificial ends of each such segment. For a free WLC segment of contour length c 1 λ m , projection of the position vector of one end, r(s 2 ), to the tangential vector of the other end, u(s 1 ), can be given by [36] r(s 1 ) · u(s 2 ) = L p (1 − e −c 1 λ m /L p ). Then the average extension of the whole chain can be estimated as in which c 2 is introduced as an unknown dimensionless factor. Equation (10) should reproduce Equation (3) in the deflection regime, which can determine c 1 = ϑ = 8A α and c 2 = 1/ϑ, so that eventually we have For a tightly confined polymer in a channel with a rectangular cross-section, Burkhardt and Yang et al. [16,19] have derived that the confinement free energy of the polymer chain can be scaled by the average length of the tube occupied by the polymer, which is the average extension of the polymer chain, as follows: Assuming that λ m should satisfy Equations (2), (11), and (12), one has or whereĤ h H h /L p andĤ w H w /L p . Equation (14) can be regarded as a new deflection length that fulfills both requirements for the free energy and statistics of geometrical quantities. We can see from Equation (14) that as long as Min(H w , H h )/L p 1, Taylor expansion of this equation yields the result in Equation (4). Inserting Equation (14) into Equation (2), and replacing λ fe by λ m , we can obtain the confinement free energy as follows: For the extension of the chain under both confinement and stretching force, as shown in Figure 1, Wang and Li [37] have suggested the force-extension relation shown in Equation (8), which now can be rewritten as On the other hand, substituting Equation (5) into Equation (3), we can obtain the classical extension relation without stretching in the tight-confinement regime, H h /L p 1, H w /L p 1, as follows: In the case off = 0, Equation (16) in this regime can be reduced to Equations (17) and (18) clearly agree in the special case (H h = H w ) of a tube with a square cross-section. In terms of the aspect ratio β = H h /H w , the right-hand side of Equation (17) is larger than the right-hand side of Equation (18) by the factor The function φ(β) has a single minimum at β = 1, corresponding to φ(1) = 1, and approaches ∞ in the limits β → 0 and β → ∞ . Thus, Equation (18) only agrees with the classical strong confinement limit (17) if the cross-section of the tube is square or close to it.

Solutions to the Fokker-Planck Equation
In order to verify the derived free energy expression, we considered the solutions to the Fokker-Planck equation, which can be used to describe the statistical behaviors of confined polymer chains. We first introduce q(r, u, s) to represent the probability that a polymer chain at arc length s has an end position vector r and an end unit tangential vector u. Then, we can have the partition function of the chain with contour length L as Z = drduq(r, u, L) and the Fokker-Planck equation [31,32] where is the potential energy per unit length due to the confinement of a rectangular tube. As suggested by Chen [19,32], the solution of Equation (20) can be expanded into a series of eigenfunctions associated with negative exponential terms of eigenvalues. By noting that the chain is sufficiently long and that the high-order eigenvalues are large enough, then the solution can be approximated by the ground state eigenfunction Ψ 0 (r, u) and eigenvalue µ 0 as follows: Then the free energy can be written as Comparing Equation (2) to Equation (23) gives Chen and his coworkers [19,31,32] have proposed an iteration method to numerically determine the ground state eigenvalue and eigenfunction. In this study, we adopted this method to calculate the confinement free energy. As examples, we considered polymer chains confined in slits with different heights H. When using Chen's method to calculate µ 0 , we set the tolerance error as 10 −4 . Figure 2 shows the comparison of the confinement free energy as a function of H/L p obtained by numerical solutions of the ground state eigenvalue, Equation (15), and Equation (2) in terms of the classical Odijk length. It can be seen from Figure 2 that free energy based on the modified deflection length had a better agreement with the numerical results than that based on the classical one.

Brownian Dynamics Simulations (Appendix A)
We used the technique of statistical dynamics simulations to verify the derived force-extension relation on polymer chains subjected to both confinement of rectangular tubes and stretching of external forces. We performed the simulations by using our GBR model for Brownian dynamics of semiflexible polymer chains in confinements [33]. This model has been successfully applied to the quantitative analysis of statistical behaviors of polymers confined on spherical surfaces [34] and in cylindrical tubes [29,30] and subjected to stretching forces [29,30]. In this GBR model, we considered the polymer chain as a discrete WLC with N identical virtual beads of radius a at different positions r k (t) = {x k (t), y k (t), z k (t)} }, where k = 1, 2, . . . , N, linked by N-1 rods with inextensible length b. The virtual beads are used to feel the hydrodynamic interactions, and angle changes of the adjacent rods are used to account for the bending deformation. As long as the position vectors of all N beads at the nth time step, denoted as r (n) = {r 1,(n) , r 2,(n) , . . . , r N,(n) }', are obtained, the new position vector at the (n + 1)th time step, r (n+1) , can be calculated from [33,34] where ∆t is the time step, δ nn is the Kronecker delta symbol, F (n) is the collective vector of internal and external forces, I − T (n) B (n) is a projection matrix (which together with T (n) d sets the inextensible constraints), χ wall (n) is the penalty displacement vector for the tube/slit walls, D (n) is the translational diffusion matrix determined through hydrodynamic interactions between beads, and ξ (n) is the vector of random force generated at each time step from a Gaussian distribution with zero mean and variance equal to We performed Brownian dynamics simulations for WLCs confined in square tubes, rectangular tubes, and narrow slits of different sizes and subjected to various stretching forces. In all simulations, the chains were initially set in a straight configuration. Confinements and constant tensile forces were then applied during chain relaxation. At the nth time step, we recorded the end-to-end distance along the z axis, z N,(n) − z 0,(n) . For each simulation, we needed to keep the steady extension states lasting a sufficiently long enough time to generate enough numbers of different equilibrium configurations. For each case, average extension of the WLC was obtained by first averaging over time and then averaging over a large number of independent trajectories with different random seeds, which was then denoted as R || . For the simulation parameters, we should note that the contour length should be larger than (at least) two times the persistence length and much larger than the deflection length scale λ m . As we were only interested in the equilibrium properties of the polymer chains, therefore specific values of the bead radius, the hydrodynamic interaction between beads, and time steps were not the key factors as long as sufficiently large numbers of different configurations of the polymer chain could be generated for averaging. In addition, the bond length should be selected to be much smaller than the deflection length scale λ m and the persistence length. For all these Brownian dynamics simulations, we chose persistence length of the chain as L p = 50 nm, the viscosity of water as η 0 = 1.005725 × 10 −4 Pa·s, and the absolute temperature as T = 293 K. Simulation parameters on bead radius a, bond length b, time step ∆t, contour length L, total simulation time, and total number of different trajectories for the ensemble average are listed in Tables 2-6 for different chains in different confinements.   Figure 3 shows convergence of the simulations for the evolution of the ensemble average of the extension, R || , for slit-and square tube-confined WLCs under stretching. It can be seen from Figure 3 that the equilibrium state could last a sufficiently long enough time to guarantee the effectiveness of time averaging.   Table 2. It can be seen from Figure 4 that the prediction in terms of the modified deflection length agreed well with the simulation results, and could be reduced to that of the classical one at a tight confinement limit. For the case of a large H w /L p , the prediction in terms of the classical deflection length exhibited a large discrepancy with the simulation results.   (8), associated with the classical deflection lengths λ fe in Equation (4), λ ext in Equation (5), and that based on Equation (17) associated with the present new deflection length λ m in Equation (16), for the normalized average extension of the WLCs stretched by different forces and confined in square tubes, rectangular tubes, and slits of different sizes, respectively.    It can be seen from Figures 5-7 that results based on the newly derived formula on the average extension of the confined WLC agreed with the simulation results very well, and those based on the classical deflection length showed an apparent discrepancy with the simulation results when the tube diameter became large or the stretching force was small.

Conclusions
Based on WLC theory and existing results on statistical properties of strongly confined polymers, we theoretically and numerically studied confinement free energy and force-confinement-extension relations of rectangular tube-confined semiflexible polymer chains under stretching in a deflection regime. We derived a modified deflection length without any adjustable parameters, which was valid for quantitative formulations of both free energy and geometrical extension. By using this deflection length scale, we obtained compact formulas on the confinement free energy and force-extension relation without any fitting parameters. Numerical analysis based on the eigenvalue problem of the governing Fokker-Planck equations and the GBR Brownian dynamics simulations justified these theoretical predictions to be valid for a much more extended range of the confinement/persistence length ratio than that based on the classical deflection length.
Author Contributions: J.W. and K.L. conceived and designed the study; K.L. performed the simulations; J.W. and K.L. analyzed the data; J.W. contributed reagents/materials/analysis tools; J.W. and K.L. wrote the paper.

Acknowledgments:
The authors wish to thank Jeff Z. Y. Chen (Waterloo, Ontario, Canada) for help on solving the eigenvalue problem, and the anonymous reviewers whose comments have greatly improved this manuscript. This research is supported by grants from the National Natural Science Foundation of China (11472119) and the 111 Project (B14044). The National Supercomputing Center in Guangzhou provided the computational resources.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. The GBR Model for Brownian Dynamics Simulations
In the GBR model, a semiflexible chain is described as N identical virtual beads of radius a, which is connected by N−1 inextensible rods of length b. Each bead has a different position in a Cartesian coordinate system, r k (t) = {x k (t), y k (t), z k (t)} , k = 1, 2, . . . , N, which is introduced for accounting hydrodynamic interactions between different chain segments. These segments are treated as hard rods, whose geometrical constraints are realized by the so-called linear constraint solver (LINCS) [33,38]. Evolution of chain configurations is characterized by a collective Brownian motion of N identical beads in solution. By denoting the time step, ∆t, and position vectors of all the beads at the current n ∆t, r (n) = {r 1,(n) , r 2,(n) , . . . , r N,(n) }', then the position vectors r (n+1) at (n + 1) ∆t can be determined from the governing Langevin equation in terms of the first-order difference method as where D (n) is the translational diffusion matrix consisting of 3 × 3 sub-blocks D jk (n) , j, k=1, 2, . . . , N, for which the Rotne-Prager tensor [39] is used to represent the hydrodynamic interactions between beads j and k. In a case where only equilibrium properties are interested, the sub-blocks D jk (n) can be simply set as Such a simplification will not affect various average results in the equilibrium state as long as a sufficiently large number of different configurations can be generated. In Equation (A11), F (n) is the collective vector, including internal bending forces and external forces; and B (n) is the associated gradient matrix related to the N − 1 time-independent constraint equations for the inextensible rods g i,(n) (r) = r i,(n) − r i+1,(n) − d i = 0 (A3) and B (n) = ∂g i,(n) (r) ∂r , i = 1, 2, · · · , N − 1 , where d i is the length of the ith rod, and the dimensions of gradient matrix B (n) are (N − 1) × 3N. Vector λ (n) consists of N − 1 Lagrange multipliers, so that B T (n) λ (n) can be regarded as collective constraint forces keeping the rod length constant. The vector ξ (n) represents random forces generated at each time step from a Gaussian distribution with zero mean and variance ξ (n) ξ (n ) = 2D (n) ∆tδ nn , where δ nn is the Kronecker delta symbol. As B (n) is time-independent, the time derivative of Equation (A3) yields Based on Equations (A1) and (A6), we can determine the Lagrange multipliers as Inserting Equation (A7) into Equation (A1) and adding an additional term −T (n) (B (n) r (n) − d) to suppress the accumulation of numerical errors, then the new position vector can be further expressed as r (n+1) = (I − T (n) B (n) )(r (n) + ∆t k B T D (n) F (n) + ξ (n) ) + T (n) d, where T (n) = D (n) B T (n) (B (n) D (n) B T (n) ) −1 , and I − T (n) B (n) is a projection matrix that sets the constraints.
For the confinement of a rectangular tube with width and height H w and H h , we used a semi-analytical technique proposed by Peters and Barenbrug [40] to treat the collective stochastic motion of N beads near the reflecting tube walls. Considering the jth bead with current position r j,(n) , we defined S j,(n) as the current distance of the jth bead from a wall. The bead can be viewed as close enough to the tube wall as long as S j,(n) ≤ 5∆tk B T/6πηa, j = 1, 2, · · · , N. (A9) If the chain is stretched by external forces, then the force vector in the above equation can be expressed as where F b (n) represents the internal forces on all beads due to bending deformation [33], and F t (n) represents the tensile forces.