Stretching Wormlike Chains in Narrow Tubes of Arbitrary Cross-Sections.

We considered the stretching of semiflexible polymer chains confined in narrow tubes with arbitrary cross-sections. Based on the wormlike chain model and technique of normal mode decomposition in statistical physics, we derived a compact analytical expression on the force-confinement-extension relation of the chains. This single formula was generalized to be valid for tube confinements with arbitrary cross-sections. In addition, we extended the generalized bead-rod model for Brownian dynamics simulations of confined polymer chains subjected to force stretching, so that the confinement effects to the chains applied by the tubes with arbitrary cross-sections can be quantitatively taken into account through numerical simulations. Extensive simulation examples on the wormlike chains confined in tubes of various shapes quantitatively justified the theoretically derived generalized formula on the force-confinement-extension relation of the chains.

1. Introduction

The statistical behaviors of semiflexible polymers confined in nano- and micro-tubes are fundamental problems in polymer physics and have been investigated both experimentally and theoretically for decades [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. A thorough understanding of these problems is very important in the development of various application techniques that exploit the effects of confinement and stretching on polymers, including genome mapping [15,16], DNA sorting [17], and DNA denaturation mapping [18,19], etc.

The conformational behavior of a polymer in a tube is determined by the competition of three interactions: bending, excluded volume interacting, and confining [1]. Several distinct conformational regimes have been identified, each with its scaling properties for both configuration and free energy. For sufficiently narrow tubes, the polymer is highly extended and lies in the Odijk regime [20,21,22]. At the opposite extreme of large tubes, the polymer falls into the classic de Gennes regime [23]. Between these limits, rich physical regimes have been gradually revealed, which include the extended de Gennes [24,25,26] and the backfolded Odijk regimes, in which the latter was first predicted by Odijk [27,28] and later characterized by Muralidhar et al. [29,30] using techniques of computer simulations.

Although the effects of geometric constraint on the semiflexible chains are relatively well described, chain statistics under the simultaneous implementation of tube-like geometric confinement and force stretching has been much less understood [31]. However, there are still a few exceptions, for example, quantitative and compact formulae on the force-extension relation of a wormlike chain (WLC) confined in cylindrical and rectangular tubes [32,33,34,35] have been derived and numerically verified by using the Brownian dynamics simulations in terms of the Generalized-Bead-Rod model (GBR) [36]. Recently, these results have been successfully applied to characterize the experimental measurements on the thermal fluctuation of internal segments of DNA confined in a nanochannel [37], the theoretical description on entropically driven motion of polymers in nonuniform nanochannels [38], and the quantitative influence on how domain topology, stability, and translation speed may determine mechanical force generation on the ribosome [39], and so on. In addition, for slit-like confinements, Haan et al. [40] have derived an interpolation formula on the force-extension relation of DNA chains by introducing an effective dimensionality, which is considered as a generalization of the Marko-Siggia force-extension relation for WLCs [41] or a more general one given by Rosa et al. [42] valid for not only WLCs, but also freely jointed chains (FJCs); Also, Taloni et al. furnished a scaling framework to identify three distinct regimes (linear, nonlinear, bulk-like) attained by a slit-confined polymer subjected to an external pulling force [43].

Nevertheless, all these studies were aimed to understand the statics and dynamics of polymers in nanoslits and nanotubes with regular square or circular cross sections. It is still not clear how the more complex shapes of a confining tube may quantitatively influence the statistics of the confined polymer chains when a stretching force is simultaneously applied. This lack of knowledge contrasts with the increasing use of nanochannels with complex shapes in nanotechnologies. Most promisingly, theoretical research works have emerged, and techniques based on tubes with triangle cross-sections have been applied as polymeric devices for biomanipulation [44,45,46,47,48,49]. For example, Reinhart et al. [44] investigated the extension of DNA chains in isosceles triangular nanochannels by using Monte Carlo simulations of a touching-bead model. They found that the channel corners may enhance the extension, relative to a circular nanochannel of the same effective size. In addition, Manneschi et al. [45] also studied the conformations of DNA chains in the triangular, rectangular, and square channels by performing coarse-grained Monte Carlo simulations. Meanwhile, recent progresses in experimental technology have been able to fabricate nanofluidic channel systems with various cross-sections [47].

However, despite these progresses, it is still very urgent and necessary to develop a generalized theoretical model to quantitatively describe the statistical behavior of a polymer chain confined in channels with arbitrary cross-sections and simultaneously subjected to external stretching forces. In the present contribution, we theoretically investigated the conformations of WLCs confined in narrow tubes with arbitrary cross-sections and established a unique force-confinement-extension relation which is quantitatively applicable to nanotubes with any cross-sections. In order to validate the derived model, based on the semi-analytical method in the treatment of particle-wall interaction developed by Peters and Barenbrug [50], we extended our GBR model [36] to realize the Brownian dynamics simulations of WLCs confined in channels with arbitrary shapes.

2. Materials and Methods

2.1. Model

Figure 1 shows the schematic of a WLC confined in a tube with an arbitrary convex cross-section and subjected to stretching by a tensile force. We assume that the geometric center of each cross-section of the tube forms the tube axis along a straight line. A set of Cartesian coordinates are introduced so that the axis is along the axis of the tube, and the profile of the cross-section of the tube at position can be described by the equation, , where . Then defines the set of all internal points enclosed by each cross-section boundary. For instance, the cross-section of a cylindrical tube of radius can be described by the equation , and represents all the spatial points inside this cylindrical tube.

Figure 1

Schematic of a wormlike chain (WLC) confined in a tube with complex shape and subjected to force stretching.

In addition, we assume that the confined chain is stretched by an applied tensile force, where is the unit vector along the -axis. The position vector at the arc length, , of the chain can be expressed as where vector is perpendicular to the -axis at position , and . The derivatives define the tangential vector and its components. As shown in [32,51], in the case of strong confinement, undulation of the chain due to thermal fluctuation is small so that . The inextensibility condition of the WLCs, [52,53,54], together with Equation (2) leads to

Based on Equation (3), one can obtain the extension of the chain along the tube axis [51]

We can further derive from Equations (1)–(3) that [32]

According to references [20,22,41], the Hamiltonian of the confined WLC without stretching force can be described as the summation of bending and potential energy as in which is the confinement potential per unit length due to the tube wall, and is the persistence length of the WLC. However, it can be very difficult to obtain analytical solutions of Equation (6) when hard wall boundary, Equation (7), is considered [55,56]. Instead, Burkhardt [55,56], Wang and Gao [32], and Wang and Li [35] found that the hard wall potential for confinements of cylindrical and rectangular tubes can be well approximated by effective harmonic potentials.

Following these studies, for the confinements of tubes with arbitrary cross-sections, we try to approximate the effect of hard wall boundaries by the harmonic potential where we assume that the parameter, , only depends on the shape/size of the confinement and the chain’s persistence length. In order to determine the parameter, , by considering Equations (2) and (5)–(8), we can have [32]

Following earlier investigations [32,41], we define the Fourier transform of

Then we have

Based on the Parseval theorem in mathematics, Equation (9) can be equivalently rewritten as

The average energy contributed by each mode [32,41] is written as

According to the Equipartition theorem in statistical physics [41,57], is equal to for two degrees of freedom, then Equation (13) becomes

Using the Parseval theorem once again, we have

According to Equations (4) and (15), we can obtain the average extension of the confined chain without stretching force as

Eventually, we can determine from Equation (16) as

Equation (17) shows the relation between the “spring” constant of the effective harmonic potential, , and the extension of the chain due to confinement without stretching force.

We further assume that Equation (17) is still valid when a stretching force, , is applied to the confined chain. In this case, the Hamiltonian becomes [20,41]

By noting that and discarding the constant term, Equation (18) can be approximately expressed as [32]

Similarly, we can obtain [32],

Further, we have [32]

Insertion of Equation (17) into Equation (22) yields

Equation (23) shows the force-confinement-extension relation for the WLC confined in a tube with arbitrary cross-section and simultaneously subjected to force stretching. From Equation (23), we can see that, as , we have , and as , Equation (23) can be reduced to which is exactly the classic force-extension relation for large forces. This fact implies that Equation (23) is exact when the tensile force is either very small or very large. For these two extremes, the normalized tube size only influences how small or large the tensile force in which Equation (23) can give an accurate prediction.

On the other hand, as the tube size is infinitely small, we have . Then the tensile forces become irrelevant to the extension of the confined chains due to the inextensibility of the WLCs. As , Taylor expansion of Equation (23) in terms of yields by omitting higher order terms.

2.2. Brownian Dynamics Simulations

In order to verify the derived force-confinement-extension relation, we consider the GBR model established by Wang and Gao [36] for the Brownian dynamics simulations of confined WLCs under stretching. This model has been successfully applied to the quantitative analysis of statistical properties of polymers confined on spherical surfaces [58] and in cylindrical and rectangular tubes [32,33,35]. In the GBR model, a WLC is depicted as identical virtual beads of radius, , connected by inextensible rods of length with the unit tangent vectors, (). The contour length of the polymer chain is . The virtual beads with coordinates, (), are introduced for simulating hydrodynamic interactions between different sections.

Motion equation of the beads can be expressed as follows in terms of the second Newton law where is a diagonal matrix containing the mass of each bead, is the position vector at time , is the configuration-dependent friction vector, is the configuration-dependent collective vector of internal and external forces including bending, stretching, and geometrical constraining [59,60,61], and is the randomly fluctuating force exerted on the beads by the surrounding fluid.

If the inertia effect is negligible, we have the overdamped equation of motion in the form

Equation (27) can be readily solved by using the Euler method where is the time step, subscript “(n)” represents the nth time step, is the translational diffusion matrix (see Ref. [36] for the detail), represents the collective random displacements at the nth time step generated from a Gaussian distribution with zero mean and variance where is the Kronecker delta symbol. As a simulation model of discrete WLCs, a linear constraint solver is applied to realize the rod length constraint [61]

If we define , and add the constraints to the motion equation through the Lagrange multiplier vector, eventually we obtain the following equation to determine the new position vector of the beads where , , and as a projection matrix which sets the constraints.

In the GBR simulation model proposed by Wang and Gao [36], a semi-analytical method on the random motion of a bead near a reflecting wall established by Peters and Barenrug [50] was taken into account for the cylindrical tube confinements. In this study, as tubes with arbitrary cross-sections are considered, therefore extension to the original GBR model is necessary to realize the complex spatial confinement. To this end, we assume that profile of the tube can be expressed as where represents the position vector of a point on the tube surface. Considering the jth bead in the chain with current position, , the distance from this bead to the point on the tube surface can be given by

In the following, we derive the minimum value of as a function of position vector, . To do this, we introduce the Lagrange function [62] as follows where is the Lagrange multiplier. The gradient of the function, , can be given by

Setting , so that

Solution of Equation (36) can give the position on the tube surface, , so that reach its minimum value. In the case that is piecewise continuous, we simply neglect the effect of undifferentiable connecting lines on the surface.

For a one-dimensional bead moving in its coordinate, , Peters and Barenrug [50] studied the effect of a flat reflecting wall at position, , based on the solution of the corresponding Fokker–Planck equation for the probability distribution of the bead position. They derived the stochastic movement of the bead after one time-step as where represents the diffusion coefficient of the bead, and the random variable with first moment and second moment . The functions and have the forms of [50]

We can imagine that as long as the time step is sufficiently small, the local curvature of the tube surface will be no longer significant to the reflecting effect of the tube wall. Therefore, similar to the case of the flat reflecting wall, when the bead is sufficiently close to the tube wall, say , we may set the next stochastic movement of the jth bead due to the reflecting tube surface, , as [36,50]

Taking into account the confinement effect in terms of Equation (40), then the positions of N beads at the next time step can be determined by

In the following, we show several examples of tubes with different profiles.

Polygonal tubes

For a straight tube with polygon cross-sections, we assume that the polygon is n-sided and locates in the plane at the position, , coordinate of its jth vertex is denoted as , . Then, the whole profile of such a tube can be expressed as where .

Elliptic tubes

For straight tubes with elliptic cross-sections, their profiles can be described as where and are positive real constants.

Wavy tubes

For periodic wavy tubes with effective wave length , assuming that profile of the first wavelength can be defined by , then we can express its overall profile as

3. Results and Discussion

Based on the extended GBR model [32,36], Brownian dynamics simulations were performed for WLCs confined in tubes of different shapes. In all the simulations, the chains were initially set in straight configurations. Tube confinements and stretching forces were then applied during the chains’ relaxation. We recorded the average extension of the chain along the -axis, , at a fixed time increment. For each simulation, the average extension of the corresponding WLC was obtained by first averaging over time and then over numerous independent trajectories obtained by using different random seeds. In all the simulations, we set the persistence length of the WLCs as Lp = 50 nm, the viscosity of the solvent, , the absolute temperature, , bead radius, , bond length, 2.0 nm, time step, 6.0 ps, contour length, 100–120 nm and total simulation time, . For each simulation, ensemble average was performed over 14–200 different trajectories with different random seeds. As we have shown in Appendix A, the extended GBR model can realize the spatial confinement of tubes with complex shapes and the equilibrium states of confined and stretched WLCs can be reached within a few tens microseconds.

To numerically verify the derived force-confinement-extension relation in Equation (23) for the confined and stretched WLCs, in the following, we consider the Brownian dynamics simulations of the chains confined in tubes of various shapes. The simulation results are represented by hollow circles, squares, rhombuses, and triangles, and corresponding theoretical predictions based on Equations (23) and (25) by solid, dashed and dash-dotted lines. Figure 2 illustrates the relation between the applied normalized stretching force, , and the normalized extension of the confined WLCs, where relevant geometrical and simulation parameters are listed in Table 1. From Figure 2, we can see that theoretical predictions based on Equation (23) are in good agreement with all the simulation results for a large range of tensile forces.

Figure 2

Comparison of simulation results on the force-confinement-extension relations of confined WLCs and their corresponding theoretical predictions based on Equation (23) under confinements of tubes with various cross-sections.

Table 1

Tube geometries and simulation parameters.

Tube Geometries		Simulation Parameters
	S(r)=x^2+y^2lc2−1	Figure 2al_c=2.5 nm,5.0 nm,7.5 nmL=100 nm,110 nm,120 nmSimulationtime=90 μs
	S(r)=y−(x−x_m)y_m+1−y_mx_m+1−x_m−y_m,x_m<x<x_m+1,m=1,2,3,4.(x_1,y_1)=−(x_3,y_3)=(−l_r1/2,−l_r2/2),(x_2,y_2)=−(x_4,y_4)=(l_r1/2,−l_r2/2).	Figure 2bl_r1=l_r2=5 nm,7.5 nm,10 nmL=100 nm,110 nm,120 nmSimulationtime=90 μs
		Figure 2cl_r1=2.5 nm,5 nm,7.5 nml_r2=5 nm,7.5 nm,10 nmL=100 nm,110 nm,120 nmSimulationtime=90 μs
	S(r)=y−(x−x_m)y_m+1−y_mx_m+1−x_m−y_m,x_m<x<x_m+1,m=1,2,3.(x_1,y_1)=(−l_t/2,−3l_t/6),(x_2,y_2)=(l_t/2,−3l_t/6),(x_3,y_3)=(0,3l_t/3).	Figure 2dl_t=10 nm,15 nm,20 nmL=100 nm,110 nm,120 nmSimulationtime=90 μs
	S(r)=y−(x−x_m)y_m+1−y_mx_m+1−x_m−y_m,x_m<x<x_m+1,m=1,2,3,4.(x_1,y_1)=−(x_3,y_3)=(−l_rsinθ,0),(x_2,y_2)=−(x_4,y_4)=(0,−l_rcosθ).	Figure 2eθ=0.83,L=100nml_r=5nm,10nm,15nmSimulationtime=120 μs
		Figure 2fθ=1.09,L=100 nml_r=5 nm,10 nm,15 nmSimulationtime=90 μs
	S(r)=y−(x−x_m)y_m+1−y_mx_m+1−x_m−y_m,x_m<x<x_m+1,m=1,2,3,4,5,6.(x_1,y_1)=−(x_4,y_4)=(−l_h,0),(x_2,y_2)=−(x_5,y_5)=(−l_h/2,−3l_h/2),(x_3,y_3)=−(x_6,y_6)=(l_h/2,−3l_h/2).	Figure 2gl_h=2.5nm,4nm,5.5nmL=100 nm,110 nm,120 nmSimulationtime=60 μs
	S(r)=S_0(x,y,z−nl),n∈ZS_0(r)={x^2+y^2+z^2−Rs2,z∈(−d,d)x^2+y^2−Rc2,   z∈(d,l−d)whered=Rs2−Rc2andR_s>R_c.	Figure 2hR_c=2.5 nm,3.75 nm,5 nmR_s=3.75 nm,5nm,6.25 nml=20.6 nm,21.6 nm,22.5 nmL=100 nm,110 nm,120 nmSimulationtime=60 μs
	S(r)=S_0(x,y,z−nl),n∈ZS_0(r)={x^2+y^2−Rc12,z∈(0,l_1)x^2+y^2−Rc22,z∈(l_1,l−l_1)	Figure 2iR_c1=2.5 nm,3.75 nm,5 nmR_c2=3.75 nm,5 nm,6.25 nml_1=30 nm,l=60 nm,L=100 nmSimulationtime=60 μs

In order to investigate how the tube size or the unstretched extension, , may influence the force-confinement-extension relation, we considered the confinements of tubes with circular and triangular cross-sections. We performed extensive Brownian dynamics simulations based on the GBR model to investigate the extension of the chains confined in tubes of various sizes and stretched by different tensile forces. Figure 3 shows the relation between the normalized extension, , under tensile force, , and the unstretched extension, , that decreases monotonically with the increase of tube size. Solid and dashed lines in Figure 3 are theoretical predictions based on Equations (23) and (25) under different tensile forces. In Figure 3a, hollow squares, circles and rhombuses represent the simulation results with simulation parameters of contour lengths, L = 100 nm, 150 nm, and 200 nm, total simulation times, 90 µs, 90 µs, and 150 µs, and total simulation trajectories, 69, 112, 69, respectively. In Figure 3b, the simulation results are still represented by the hollow squares, circles and rhombuses, for the contour lengths, L = 100 nm, 150 nm, and 200 nm, total simulation times, 90 µs, 90 µs, and 300 µs, and total simulation trajectories, 84, 50, and 140, respectively. We can observe from Figure 3 that numerical simulations agree very well with the theoretical predictions for wide ranges of tensile forces and tube sizes, though the derivation of Equation (23) assumes strong confinement.

Figure 3

Comparison of simulation results on the force-confinement-extension relations of the confined WLCs and their corresponding theoretical predictions based on Equations (23) and (25) under confinements of tubes with cross-sections of (a) circular, and (b) triangular shapes. Note that the simulation data for show very large horizontal error bars, which are not plotted in the figure for clarity.

4. Conclusions

We studied the statistical mechanics behavior of WLCs confined in narrow tubes with arbitrary cross-sections and derived a generally applicable compact expression on the force-extension relation of the confined chains. We also extended our generalized bead-rod (GBR) model for Brownian dynamics simulations of confined polymer chains subjected to force stretching, so that the confinement effects arising by tubes with complex shapes to the chains can be quantitatively accounted for during the numerical simulations. Sufficiently large numbers of simulations on the WLCs confined in tubes of various shapes were performed, which quantitatively justified the theoretically derived generalized force-confinement-extension relation.