Coil-globule transition of a single semiflexible chain in slitlike confinement

Single polymer chains undergo a phase transition from coiled conformations to globular conformations as the effective attraction between monomers becomes strong enough. In this work, we investigated the coil-globule transition of a semiflexible chain confined between two parallel plates, i.e. a slit, using the lattice model and Pruned-enriched Rosenbluth method (PERM) algorithm. We find that as the slit height decreases, the critical attraction for the coil-globule transition changes non-monotonically due to the competition of the confinement free energies of the coiled and globular states. In wide (narrow) slits, the coiled state experiences more (less) confinement free energy, and hence the transition becomes easier (more difficult). In addition, we find that the transition becomes less sharp with the decreasing slit height. Here, the sharpness refers to the sensitivity of thermodynamic quantities when varying the attraction around the critical value. The relevant experiments can be performed for DNA condensation in microfluidic devices.

because proteins usually experience spatial confinement in cells. Confinement has also been shown to facilitate compaction of DNA by crowding agents 33 . In addition, Das and Chakraborty have performed Brownian dynamic simulations for a collapsing flexible chain in slitlike confinement 34 . However, these simulations focus on the collapse kinetics and provide limited information regarding the critical attraction for the coil-globule transition. The effect of confinement on the critical attraction of the transition has been studied for flexible chains using simulations 35,36 .
In this work, we use simulations to calculate the densities of states of confined semiflexible chains, and the densities of states are used to determine the critical attraction and free energy landscape of coil-globule transitions.

Simulation and theory
The polymer model. The semiflexible chain is modelled as a series of N m monomers in a cubic lattice of spacing a. The model includes four interactions: bending energy along the chain, self-avoiding interaction between monomers, attraction between monomers, and repulsion between monomers and the slit walls. The bending energy is defined as where θ i is the angle formed by two adjacent monomers. The bending angle θ = π is forbidden due to self-avoidance. κ bend , as well as all other energies in this manuscript, are presented in units of k T B . The bending energy leads a correlation in the orientations of segments as well as the persistence length (See supporting information): The self-avoiding interaction between two monomers is described as A pairwise attraction is applied to the contact pair, which is defined as a pair of two non-bonded monomers with distance equal the lattice spacing a. For two monomers with positions  x i and  x j , the attraction interaction follows The attraction strength is ε in the unit of k T B . The total attraction energy for the entire chain is εN c , where N c is the total number of contact pairs. The hardcore repulsions with slit walls are described as where z i is the z-coordinate of i-th monomer, and H is the slit height. Then, the total interaction energy for an allowed conformation is In simulations, we vary the bending rigidity, the contour length, and the slit height to obtain a systematic picture of confinement effect on the coil-globule transition. Due to the expensive computational cost, we only investigate three contour lengths: = N 256 m , 512 and 1024, and three bending rigidities: κ = , . 0 Effective density of states. For every allowed configuration (no overlap with self or slit walls), the total energy is Then, the partition function for the semiflexible chain is is the inverse of thermal energy, and the summation is over all allowed configurations. For the convenience of investigating ε -dependence, the partition function can be written as where the summation is over all allowed configurations for a specific N c . Note that in the case of flexible chains with = E 0 bend , the effective density of states is simply the number of states with N c contact pairs, a commonly calculated quantity in previous studies of the coil-globule transition of flexible chains 16,37,38 . The effective density of states ( ) g N eff c is a quantity independent of ε , and hence we can apply ( ) g N eff c for any given value of ε in the calculation of other thermodynamic quantities, which facilitates the investigation of coil-globule transition while varying ε.
Simulation algorithm. We used the flat version 39 of Pruned-enriched Rosenbluth method (PERM) 40 , called flatPERM, to efficiently sample chain configurations. Prellberg and Krawczyk have previously provided comprehensive description of the flatPERM algorithm 39 for simulations of flexible chains in bulk. We adapt the flatPERM algorithm for semiflexible chains in confinement. The detail of the algorithm is presented in the supporting information.

Results and Discussions
Simulation results and analysis. Figure 1 shows randomly chosen conformations with the contact number  bundle with a number of local defects. Note that in the off-lattice model 14,19 as well as the experiments 19,20 , the globular state can assume not only a bundle structure but also a toroidal structure.
Next, we proceed to more quantitative results. Figure  . The peaks were also observed in the previous simulations of flexible chains 16,37,38 . The cause of these peaks is that sufficiently long chains are likely to form a certain contact number randomly. Beyond the peaks, (  . Based on the numerical values of ( ) g N eff c , we performed a series of analysis for the coil-globule transition as presented in the following figures. At first, we calculate the fluctuation in the number of contact number for any given attractive strength ε using the following equation: Figure 3 shows the fluctuation as a function of the attractive strength for κ = 3 bend k T B and = N 1024 m . In both bulk and in slits, the fluctuation exhibits a peak, corresponding to the critical attractive strength ε ⁎ of the coil-globule transition. With the decreasing slit height, the peak monotonically broadens and reduces in magnitude. The critical attraction ε ⁎ , i.e. the peak location, is plotted in Fig. 4 for each slit height. The critical attractions for the chains with different bending rigidities and contour lengths are included for comparison. It is worthwhile to mention that the coil-globule transition is a first-order transition only when the bending rigidity is larger than a critical value κ ⁎ bend . The critical bending rigidity in free space is found to be κ ≈  All curves in Fig. 4 exhibit the same trend: with decreasing slit height the critical attraction first decreases to a minimum value for a critical slit height ⁎ H and then increases. The critical slit heights corresponding to the minimum ε ⁎ are plotted in Fig. 5. The critical slit height generally becomes larger for the longer and more flexible chains. For κ = bend 3 k T B , the critical slit height is always = ⁎ H a 2 for the three contour lengths. However, if we look into the shapes of the three curves for κ = bend 3 k T B , the critical slit height appears to slightly increase with the increasing contour length. For example, if we assume the critical slit height locating at the cross point of two segments: the connection of two data points = H a and = H a 2 , and the connection of two data points = H a 3 and = H a 4 , we will have = .    Fig. 4c, such non-monotonic behaviour may be related to the special feature of = H a. Specifically, = H a corresponds to a purely two-dimensional case. The change of dimensionality from 3-D or quasi-2D to pure 2-D may lead to some fundamental differences that result in the non-monotonic behavior. For example, the projection of the chain onto a slit wall is non-self-crossing for = H a, while in all other slits, the projected chains can cross itself. Whether such special feature of = H a is responsible for the non-monotonic behavior can be clarified by simulation of very long chains.
The critical attraction as a function of the slit height has been calculated for the flexible chains in the previous study by Hsu and Grassberger 35 , and is also plotted in the thick green line of Fig. 4(a). The difference between our results and the green line is probably because Hsu and Grassberger 35 35 . In addition to the critical attraction, the curves in Fig. 3 also provide the transition widths or the sharpness of the coil-globule transition. To quantify the sharpness of the transition, we plot the full width at half maximum as a function of the slit height in Fig. 6. With the decreasing slit height, the coil-globule transition becomes less sharp. Here, the sharpness refers to the sensitivity of thermodynamic quantities, e.g. the size of chain, when varying the attraction around the critical value.
Then, we turn to the free energy landscape for the coil-globule transition. For a given attractive strength, the probability of the chain containing N c contact pairs can be calculated using the equation The probability can be converted to the free energy using c c Using the critical attractions in Fig. 4, the free energy landscapes at phase equilibriums for = N 1024 m and κ = kT 3 bend are calculated and shown in Fig. 7. All curves are shifted to offset the local minimum in the range ≤ / ≤ .
In both bulk and slits, the free energy landscapes exhibit two local minima, corresponding to the coiled and globular states. The locations of two local minima become closer with the decreasing slit height. Free energy barriers appear between the two states. To quantify the barriers, we define the free energy barrier as the difference between the free energy maximum in the range . ≤ / ≤ . . Figure 8 shows the free energy barrier as a function of the slit height. The barrier monotonically decreases with the decreasing slit height. It is important to mention that the free energy barriers shown in Figs 7 and 8 cannot be directly used for kinetics of coil-globule transition. This is because the reaction coordinate in Fig. 7 is the contact number, and the motion of the chain is not in the space of contact number but in the space of XYZ coordinates.
To understand why confinement induces a non-monotonic change in coil-globule transition, we analyzed the confinement free energy as a function of the contact number, while the confinement free energy is defined as  Figure 9 shows the confinement free energy using the data in Fig. 2. There are two competing effects of confinement on the free energy. First, in the most cases, the confinement free energy decreases with the increasing contact number because larger contact number corresponds to more compact configurations, which are less likely to be restricted by the geometrical confinement. Second, in very strong confinement with = H a, the confinement free energy increases with the increasing contact number. This is because the strong confinement deforms the bundle structure and increases the number of surface monomers. The surface monomers have less contact partners, and hence have higher energy than the inner monomers. The difference in energy between the surface monomers and inner monomers is commonly referred to as the surface energy. When the slit height decreases from = H a 2 to = H a, the conformation switches from a double-layer structure to a single-layer structure, which leads to the large increase of surface energy and is responsible for the increase of ε ⁎ from = H a 2 to = H a. With the information in Fig. 9, we can understand the non-monotonic behaviors of ε ⁎ versus H in a simple way, as shown in Fig. 10. In wide slits with   R H R globule bulk coil bulk , only the coiled state experiences a significant confinement free energy. Here, R globule bulk and R coil bulk refer to the conformation size of the globular and coiled states in free space, respectively. In narrow slits with  H R globule bulk , the confinement free energy of the globular state increases rapidly with the decreasing slit height, mainly because the deformation of the globular state exposes more monomers to the surface and increases the surface energy. The critical slit height between the two trends is expected to  be close to R globue bulk . For a spherical shape, R globue bulk is relatively easy to define. For a bundle conformation of a semiflexible chain, R globue bulk should refer to the smallest size among three dimensions, because the bundle conformation can rotate to avoid deformation. Due to the bundle conformation, R globue bulk of a semiflexible chain is smaller than the one of a flexible chain for a fixed N m . As a result, the critical slit height decreases with the increasing bending rigidity in Fig. 5.

Implications in experiments.
DNA condensation is an example of the coil-globule transition of nearly homogenous chains. Extensive experiments of DNA condensation have been performed in bulk using various approaches to induce the transitions, such as adding multivalent ions, ethanol or crowding agents 13 . The transitions in some experiments have been confirmed to be a first-order transition 13 . The chain parameters used in our simulations are not far away from parameters in common DNA experiments. There are two dimensionless parameters to describe a semiflexible chain: / L L p and / w L p , where L is the contour length, and w is the effective chain width. In our simulations, the effective chain width can be approximated as a because the closest distance between two monomers is a. Then, the dimensionless parameters are / = / . ≈ L L 1024 5 51 186 p and / = / . ≈ . w L 1 5 51 0 18 p . Applying these two dimensionless parameters to DNA and using a typical DNA persistence length of 50 nm 41 , we can estimate that the chains in our simulations correspond to a contour length of approximately 9.3 μ m and a chain width of approximately 9 nm. These two values are close to the parameters of λ -DNA, widely used in DNA condensation experiments. The contour length of λ -DNA is approximately 16 μ m and 22 μ m respectively, before and after YOYO-1 labelling 42 . For the coiled state, the effective chain width typically ranges from 2 nm to 20 nm, depending on the ionic strength of solution 43 .
These values indicate that direct experimental testing of our results in microfluidic confinement is possible. Our simulation results suggest the critical attraction can be reduced by confinement, i.e. the transition requires a lower concentration of multivalent ions or ethanol.

Conclusions
Our simulations reveal that with decreasing slit height the critical attraction for the coil-globule transition first decreases and then increases. The non-monotonic behavior is due to the competing of the confinement free energies of the coiled state and the globular states. Our results suggest that the critical slit height corresponding to the minimum critical attraction mainly depends on the smallest dimension of the globule transition. We expect the critical slit height increases with the contour length for stiff chains, but we do not directly observe the trend in our simulations due to insufficiently long chains. Such trend needs to be validated in the future simulations of very long semiflexible chains. In addition, our analysis suggests the coil-globule transition becomes less sharp with the decreasing slit height.