1. Introduction

Higher-order structures of DNA are the subject of extensive studies for their significance in genome packaging, gene regulation, and enzyme recognition [1,2]. Wrapping, crossover, and braiding are among the most fundamental and essential higher-order structuresof DNA. Since all these higher-order structures are associated with superhelical conformations of DNA, chirality, i.e., right- or left-handedness, of these superhelical motifs should play important roles. This study thereby explores a possible mechanism for the selection of chirality in wrapping, crossover, and braiding of DNA based on a coarse-grained elastic model. Since not only DNA but also many other biopolymers assume helical conformations, roles of helical chirality are of rather universal interest [3,4,5,6,7,8,9].

Wrapping of DNA around core particles is an indispensable step for genome packaging in eukaryotes. As is known, eukaryotic DNA achieves a highly compact and hierarchical packaging in a cell nucleus, which is called chromatin [10,11,12,13,14,15]. At the lowest level of the hierarchy of chromatin, DNA wraps around a protein core particle called histone octamer about 1.75 times in a left-handed manner [1,2]. The uniform left-handedness of this wrapping should be important for the clustering of nucleosomes mediated by charge distributions [16] to form further higher-order structures in chromatin. While there are detailed researches on the wrapping structures of DNA in nucleosomes [16,17,18,19,20,21,22] including researches on the roles of electrostatic effects [23,24] and histone tails [25], the present study tries to develop a minimal model to account for the stability of left-handed wrapping of DNA around a core particle in terms of asymmetric elasticity of DNA. In addition to wrapping of DNA around a core particle, wrapping and unwrapping of DNA around rod-like molecules [26,27], including carbon nanotubes [28,29,30,31,32], is also an interesting subject. The present study also addresses the issue of chirality of DNA adsorbed on a uniform rod.

Crossover and braiding are also important motifs of DNA for genome packaging, gene regulation, and enzyme recognition [33,34,35,36,37,38,39,40,41,42,43,44,45,46]. Recently, Timsit and Várnai [36,37,38] have suggested that right-handed crossover of DNA is more stable than left-handed one in solution due to groove-backbone interaction especially in the presence of divalent cations. Cherstvy [39] argued the chirality of cholesteric twist of DNA in terms of electrostatic and steric effects. Moreover, Kornyshev et al. [40,41,42,43,44] suggested that braiding of DNA molecules could be driven by electrostatic interactions and plays important role in the stabilization of homologous paring. They also suggested that left-handed braiding is preferred over right-handed one due to the helical alignment of charges on DNA. It is getting clearer that chirality of crossover and that of braiding are intimately related: Right-handed crossovers occur in a left-handed braiding of DNA, and left-handed crossovers occur in right-handed braiding of DNA [36,37,38,39,40,41,42,43,44]. While these works have highlighted the electrostatic effects and steric effects on the atomistic-level, our present study explores the roles of mechanical interaction for the selection and stability of left-handed braiding of DNA in terms of the intrinsic asymmetric elasticity, or nonlinear coupling, between bending and writhing of DNA, as the intrinsic property of helical double-stranded polymer.

Indeed, the asymmetric elasticity of DNA has been the subject of extensive studies both theoretically and experimentally for decades. The right-handed helical nature of DNA geometry gives rise to asymmetric coupling among bending, stretching, twisting, and writhing, as the intrinsic mechanical property. For example, an asymmetric and nontrivial coupling between stretching and twisting of DNA has been revealed experimentally in [47,48]. Asymmetric coupling among stretching, twisting, and supercoiling of DNA is also investigated experimentally [49,50,51,52,53]. Coupling between bending and twisting of DNA has been investigated theoretically with implications for higher-order structures [54,55]. Suchfundamental elastic properties of DNA are expected to play significant roles in the formation and maintenance of its higher order structures and is of scientific value to obtain deeper understanding both from experimental and theoretical studies.

Recently, we have explored a possible asymmetric coupling between bending and writhing of a double-stranded helical chain as a prototypical coarse-grained model of B-form DNA [56], where the effect of thermal fluctuations was not taken into consideration. In the present study, we extend our previous model by introducing the effect of thermal fluctuations using a Monte Carlo method. In addition, the present study modifies the DNA model by incorporating the difference between the major groove and the minor groove of DNA. We then study chiral specificity not only in wrapping of DNA around a spherical core particle, but also in wrapping of DNA around a rod and in crossover and braiding of two interacting DNA molecules from a rather common perspective of the asymmetric bend-writhe coupling of DNA: Due to effective DNA-core or DNA-DNA interaction and thermal fluctuations, DNA molecules are generally forced to bend in solution. Even though suchbending has no bias or chirality, it can induce a bias towards left-handed writhing through the intrinsic asymmetric bend-writhe coupling in DNA. As a result, our DNA model provokes the preference for left-handed wrapping and braiding.

Since our present study primarily aims at developing a minimal elastic model of a double-stranded helical chain exhibiting chiral selectivity in wrapping, crossover, and braiding, we simplify the modeling as much as possible. Thus, our DNA model contains only harmonic potentials for bonding and bending of the two mutually intertwining chains. Moreover, we use the simple Morse potential for the effective DNA-core and DNA-DNA interactions without taking into consideration the details of the electrostatic effects such as helical alignments of charges, the effects of image charges, and heterogeneity of dielectricity [57,58]. These can be the drawbacks of our present model as compared to other sophisticated coarse-grained models of DNA [59,60,61,62]. However, it may still be of scientific interest that our simplified model of a double-stranded helical chain can exhibit chiral selectivity in the formation of higher-order structures. In addition, because of the simplicity of the model, the results of the present study may be generalized to the elasticity of other helical biopolymers.

This paper is organized as follows. In Section 2, we introduce the coarse-grained model of DNA. In Section 3, after exploring the fundamental elasticity of the DNA model, we investigate chiral selection in wrapping of the DNA model around a spherical core particle and around a uniform rod using Monte Carlo simulations. We then investigate chiral selection in crossover and braiding of a pair of DNA models. This paper concludes in Section 4 with discussions and implications for future study.

2. Model of Double-Stranded DNA

2.1. Coordinates of the DNA model

The model of DNA in the present study consists of two elastic chains intertwining around a central backbone as presented in Figure 1 for a short segment of DNAwith 5 base pairs. This model is an extension of our previous model introduced in Ref. [56]. For the DNA model with $N$ base pairs, the central backbone consists of $N$ nodal points, whose position vectors are represented by three-dimensional vectors ${r_i}\left( {i = 1, \cdots , N} \right)$, connected with $N - 1$ links. We fix the distance between two adjacent nodal points (i.e., the length of each link) of the central backbone to $b \equiv 0.34 {\text{ nm}}$, which corresponds to the distance between two neighboring base pairs. As shown in Figure 1, $N - 2$ bending angles ${\Theta _i} \left( {i = 2, \cdots , N - 1} \right)$ and $N - 3$ dihedral angles ${\Phi _i}\left( {i = 2, \cdots , N - 2} \right)$ parameterize the conformation of the central backbone as [63,56], ${r_1} = {\left( {0, 0, 0} \right)^T}$, ${r_2} = {\left( {0, 0, b} \right)^T}$, and

${r_i} = {r_{i - 1}} + {\bf{R}}_3^2\left( {\Phi _1}, {\Theta _2} \right){\bf{R}}_4^3\left( {\Phi _2}, {\Theta _3} \right) \cdots {\bf{R}}_i^{i - 1}\left( {\Phi _{i - 2}}, {\Theta _{i - 1}} \right)b, \left( {i = 3, \cdots , N} \right),$

(1)

where $b \equiv {\left( {0, 0, b} \right)^T}$, ${\Phi _1} \equiv 0^\circ$, and

$\mathbf{R}_{i}^{i-1}\left( {{\text{ }\!\!\Phi\!\!\text{ }}_{i-2}},{{\text{ }\!\!\Theta\!\!\text{ }}_{i-1}} \right)=\left( \begin{matrix} \cos {{\text{ }\!\!\Phi\!\!\text{ }}_{i-2}} & -\sin {{\text{ }\!\!\Phi\!\!\text{ }}_{i-2}} & 0 \\ \sin {{\text{ }\!\!\Phi\!\!\text{ }}_{i-2}} & \cos {{\text{ }\!\!\Phi\!\!\text{ }}_{i-2}} & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right)\left( \begin{matrix} 1 & 0 & 0 \\ 0 & \cos {{\text{ }\!\!\Theta\!\!\text{ }}_{i-1}} & -\sin {{\text{ }\!\!\Theta\!\!\text{ }}_{i-1}} \\ 0 & \sin {{\text{ }\!\!\Theta\!\!\text{ }}_{i-1}} & \cos {{\text{ }\!\!\Theta\!\!\text{ }}_{i-1}} \\ \end{matrix} \right).$

(2)

The ranges of the angels are $- {180^\circ } ≤ {\Phi _i} < {180^\circ }$ and $- 180^\circ ≤ {\Theta _i} < 180^\circ$ respectively in the present study.

Figure 1(b) shows the two elastic chains intertwining around the central backbone. These two chains, which we call P-chain and Q-chain, represent the two sugar-phosphate chains of DNA. Each of these two chains consists of $N$ nodal points connected with $N - 1$ elastic bonds (harmonic springs) for the model with $N$ base pairs. Position vectors of these nodal points of P-chain and Q-chain are represented by ${P_i}$ and ${Q_i}$ respectively. Distance between $i$-th nodal point of the central backbone and $i$-th nodal points of P-chain and Q-chain is fixed to ${\sigma _0} \equiv 1.0 {\text{ nm}}$, which corresponds to the radius of DNA. Nodal points ${P_i}$ and ${Q_i}$ are linked with a rigid rod representing a hydrogen-bonded base pair. Since we are interested in higher-order structures of DNA and not in the opening of base pairs in the present study, we fix the distance between these two nodal points ${P_i}$ and ${Q_i}$ to ${h_0}$, whose exact value will be specified in Eqs. (7) and (8). Thus, after introducing two kinds of vectors,

$\eqalign{ & {p_i} = {\left( {{\sigma _0}\cos {\phi _i}, {\sigma _0}\sin {\phi _i}, 0} \right)^T}, \quad \quad \quad \left( {i = 1, \cdots , N} \right), \cr & {q_i} = {\left( {{\sigma _0}\cos \left( {{\phi _i} + {\omega _0}} \right), {\sigma _0}\sin \left( {{\phi _i} + {\omega _0}} \right), 0} \right)^T}, \quad \quad \quad \left( {i = 1, \cdots , N} \right), \cr}$

(3)

where ${\phi _i}\left( {i = 1, \cdots , N} \right)$ is the internal twist angle characterizing the relative orientation of $i$-th base pair with respect to -th frame attached to the central backbone (see Figure 1(a)), and ${\omega _0}$ is the common constant angle between the two vectors ${P_i} - {r_i}$ and ${Q_i} - {r_i}$, the positions of nodal points of the P- and Q-chains are expressed as, ${P_1} = {p_1}, {Q_1} = {q_1},$ ${P_2} = {r_2} - {p_2}, {Q_2} = {r_2} - {q_2},$

${P_i} = {r_i} - {\bf{R}}_3^2\left( {\Phi _1}, {\Theta _2} \right){\bf{R}}_4^3\left( {\Phi _2}, {\Theta _3} \right) \cdots {\bf{R}}_i^{i - 1}\left( {\Phi _{i - 2}}, {\Theta _{i - 1}} \right){p_i}, \quad \quad \quad \left( {i = 3, \cdots , N} \right),$

(4)

${Q_i} = {r_i} - {\bf{R}}_3^2\left( {\Phi _1}, {\Theta _2} \right){\bf{R}}_4^3\left( {\Phi _2}, {\Theta _3} \right) \cdots {\bf{R}}_i^{i - 1}\left( {\Phi _{i - 2}}, {\Theta _{i - 1}} \right){q_i}, \quad \quad \quad \left( {i = 3, \cdots , N} \right).$

(5)

Thus, when ${\phi _i}$ changes, the nodal points ${P_i}$ and ${Q_i}$ move together, keeping their mutual distance ${h_0}$, on a circle of radius ${\sigma _0} = 1.0 {\text{ nm}}$, which is centered at $i$-th nodal point of the central backbone and is locally perpendicular to the central backbone.

In the present model of DNA, bending angles ${\Theta _i}\left( {i = 2, \cdots , N - 1} \right)$ and dihedral angles ${\Phi _i}\left( {i = 2, \cdots , N - 2} \right)$ of the central backbone are the key variables that determine global conformation of the DNA model, i.e., bend and writhe of the central backbone of DNA respectively. Positive sign of the dihedral angle ${\Phi _i}$ represents right-handed writhe at $i$-th link of the central backbone, while negative sign of the dihedral angle ${\Phi _i}$ represents left-handed writhe at $i$-th link of the central backbone. The internal twist angles ${\phi _i}\left({i = 1, \cdots , N} \right)$ can also change and characterize the degree of mutual twist of the two sugar-phosphate chains. Since the central backbone is for the coordinatization of the two chains, P-chain and Q-chain, and does not directly influence the elastic energy of the system as we see below, figures hereafter do not show the central backbone upon visualization of the DNA model.

2.2. Equilibrium conformation of the DNA model

The equilibrium conformation of the DNA model is shown in Figure 2. We assume that the central backbone of the DNA model, which is not shown in Figure 2(a), takes a straight conformation without writhe, i.e., ${\Theta _i} = 0^\circ \left( {i = 2, \cdots , N - 1} \right)$ and ${\Phi _i} = 0^\circ \left( {i = 2, \cdots , N - 2} \right)$, at the equilibrium conformation. We then assume that each of the two sugar-phosphate chains, P-chain and Q-chain, of the DNA model completes one helical cycle per every 10 nodal points based on the fact that the double strand of B-form DNA completes one helical cycle per every 10 base pairs approximately [1,2]. Thus, we assume that every two neighboring nodal points in each of P-chain and Q-chain form the constant angle of ${\phi _{\text{0}}} \equiv 36^\circ$ with respect to the central backbone at the equilibrium conformation when vertically projected onto the plane perpendicular to the central backbone (see Figure 2(b)). As a result, the pitch of each helical turn of the double strand is 3.4 nm.

The internal twist angle ${\phi _i}$ for $i$-th nodal point of the P-chain is set to

${\phi _i} = \left( {i - 1} \right) \times {\phi _0}, \quad \quad \quad \left( {i = 1, \cdots , N} \right),$

(6)

at the equilibrium conformation. In order to make the width of major groove 2.2 nm and that of minor groove 1.2 nm at the equilibrium conformation [1,2], we assume that the constant angle ${\omega _0}$ between the two vectors ${P_i} - {r_i}$ and ${Q_i} - {r_i}$ (for all $i$) to be

${\omega _0} = 1.2 \times \frac{{360^\circ }}{{3.4}} = 127.06^\circ .$

(7)

Thus the common length of the rigid rod, ${h_0}$, connecting each pair of nodal points, ${P_i}$ and ${Q_i}$, is

${h_0} = 2{\sigma _0}\sin \left( {\frac{{{\omega _0}}}{2}} \right).$

(8)

Based on the above settings, the equilibrium length ofeach elastic bond connecting the two adjacent nodal points in P- or Q-chain is determined to be

${l_0} = \sqrt {2\sigma _0^2\left( {1 - \cos {\phi _0}} \right) + {b^2}} .$

(9)

Similarly, the equilibrium bending angle at each nodal point of the P- and Q-chain is

${\theta _0} = {\cos ^{ - 1}}\left[{\frac{{2\sigma _0^2\left( {1 - \cos {\phi _0}} \right)\cos {\phi _0} + {b^2}}}{{2\sigma _0^2\left( {1 - \cos {\phi _0}} \right) + {b^2}}}} \right].$

(10)

2.3. Energy functions of the DNA model

We introduce here energy functions into the DNA model. Since we are interested in rather general properties of double-stranded helical chains, we introduce only minimal energy functions necessary for general double-stranded helical chains in the present study. Specifically, our DNA model contains bonding energy of the P- and Q-chains, ${V_{{\text{bond}}}}$, and bending energy of the P- and Q-chains, ${V_{{\text{bend}}}}$. Thus the total energy of the DNA model is ${V_{{\text{DNA}}}} = {V_{{\text{bond}}}} + {V_{{\text{bend}}}}$.

The total energy for bonding in the P- and Q-chains is given by

${V_{{\text{bond}}}} = \mathop \sum \limits_{i = 1}^{N - 1} \frac{1}{2}{k_{{\text{bond}}}}{\left( {\left| {{P_{i + 1}} - {P_i}} \right| - {l_0}} \right)^2} + \mathop \sum \limits_{i = 1}^{N - 1} \frac{1}{2}{k_{{\text{bond}}}}{\left( {\left| {{Q_{i + 1}} - {Q_i}} \right| - {l_0}} \right)^2},$

(11)

where ${l_0}$ is the equilibrium distance, which is common to all the bonds in the P- and Q-chains (see Figure 2(a)), and has been determined in Eq. (9). The parameter ${k_{{\text{bond}}}}$ is the bonding rigidity, which is also common to all bonds in the P- and Q-chains and will be determined below. The total energy for bending of the P- and Q-chains is given by

${V_{{\text{bend}}}} = \mathop \sum \limits_{i = 2}^{N - 1} \frac{1}{2}{k_{{\text{bend}}}}{\left( {\theta _i^{\left( P \right)} - {\theta _0}} \right)^2} + \mathop \sum \limits_{i = 2}^{N - 1} \frac{1}{2}{k_{{\text{bend}}}}{\left( {\theta _i^{\left( Q \right)} - {\theta _0}} \right)^2},$

(12)

where $\theta _i^{\left( P \right)}$ and $\theta _i^{\left( Q \right)}$ are the bending angles in the P- and Q-chains respectively at th nodal points, and ${\theta _0}$ is the common equilibrium angle for all these bending angles (see Figure 2(a)) and has been determined in Eq. (10). The parameter ${k_{{\text{bend}}}}$ is the bending rigidity, which is common to all the bending angles in the P- and Q-chains and will be determined below. Although the bonding and bending energies of Eqs. (11) and (12) are attributed to the P- and Q-chains, it may be possible to assume that these energies take into account the effects of stacking interactions among DNA bases in the present coarse-grained model.

As far as we are aware, there are no experiments that directly measure the spring parameters of the present model. Therefore, we estimate these parameters based on a comparison between simple numerical experiments of the present model and existing experimental data as follows. At a coarse-grained scale, total energy of the DNA model with $N$ base pairs may be approximated as a quadratic function of the bending angles of the central backbone, ${\Theta _i}$, as

$E = \frac{1}{2}K\mathop \sum \limits_{i = 2}^{N - 1} {\Theta _i^2} = \frac{1}{2}K\left( {N - 2} \right){\Theta ^2},$

(13)

where, in the second equality, we have assumed that all bending angles of the central backbone ${\Theta _i}$ take the same value, ${\Theta _i} = {\Theta }\left( {i = 2, \cdots , N - 1} \right)$. The parameter $K$ is the effective bending rigidity for each bending angle of the central backbone of the DNA model.As is known, persistence length of DNA is ${L_{\text{P}}} = 50{\text{ nm}}$, which indicates that the bending rigidity of our model, $K$, should be [58]

$K = \frac{{{L_{\text{P}}}{k_{\text{B}}}T}}{b} = 609 [{\text{pNnm], }}$

(14)

where ${k_{\text{B}}} = 1.3806 \times {10^{ - 23}}{\text{J}}/{\text{K}}$ is the Boltzmann’s constant, and we have assumed that temperature is $T = 300 {\text{ K}}$. We determine the spring constant for bonding ${k_{{\text{bond}}}}$ and that for bending ${k_{{\text{bend}}}}$ so that our DNA model reproduces approximately the same bending rigidity as $K$ in Eq. (14). Based on this consideration, throughout the present article, we set these constants as

${k_{{\text{bond}}}} = 1000 \left[{{\text{pN/nm}}} \right], \quad \quad \quad {k_{{\text{bend}}}} = 7000 {\text{ [pNnm]}}.$

(15)

The solid curve in blue in Figure 3(a) shows the dependence of the total energy of our DNA model ${V_{{\text{DNA}}}} = {V_{{\text{bond}}}} + {V_{{\text{bend}}}}$ with $N = 100$ base pairs on the bending angle ${\Theta _i} = {\Theta }\left( {i = 2, \cdots , N - 1} \right)$ of the central backbone, where all the bending angles are set to the same value ${\Theta }$ and changed together from ${\Theta } = - 3^\circ$ to $+ 3^\circ$. In obtaining the solid curve in blue in Figure 3(a), all dihedral angles of the central backbone are fixed to zero, ${\Phi _i} = 0^\circ \left( {i = 2, \cdots , N - 2} \right)$, and the internal twist angles ${\phi _i}\left( {i = 1, \cdots , N} \right)$ are fixed to the equilibrium values prescribed Eq. (6). In Figure 3(a), the broken curve in red shows the quadratic function of Eq. (13) for $K = 609 {\text{ pNnm}}$. From Figure 3(a), we see that the bending rigidity of the present model is fairly consistent with the known bending rigidity of real DNA for small bending angles, although nonlinearity arisesin the present model of double-stranded DNA. In Figure 3(b), broken curve in red shows the total bonding energy of the present model ${V_{{\text{bond}}}}$ in Eq. (11), and solid curve in green shows the total bending energy of the present model ${V_{{\text{bend}}}}$ in Eq. (12), which are summed to give the solid curve in blue in Figure 3(a). We thus see that the contribution of the bending rigidity in the P- and Q-chains is greater than that of the bonding rigidity of the P- and Q-chains in the overall bending rigidity of the DNA model.

2.4. Monte Carlo algorithm

We use the Metropolis Monte Carlo algorithm [66,67] toincorporatethe effects of thermal and noisy environments. Suppose that we have a “current” conformation of the DNA model with total energy $V$. Then, we randomly choose one of bending angles of the central backbone ${\Theta _i}\left( {i = 2, \cdots , N - 1} \right)$, dihedral angels of the central backbone ${\Theta _i}\left( {i = 2, \cdots , N - 2} \right)$, and internal twist angles ${\phi _i}\left( {i = 1, \cdots , N} \right)$, and perturb it by $+ 1^\circ$ or $- 1^\circ$ randomly to obtain a “trial” conformation. If the total energy of the system for the trial conformation is ${V_{{\text{tr}}}}$ lower than that of the current conformation, the trial conformation ${V_{{\text{tr}}}}$ is accepted. If the total energy of the system for the trial conformation is higher than that of the current conformation, we accept the trial conformation with the probability,

$p = \exp \left( { - \frac{{{V_{{\text{tr}}}} - V}}{{{k_{\text{B}}}T}}} \right).$

(16)

Throughout the present study, temperature is set to .

Note that the Monte Carlo algorithm of the present study uses only the three kinds of angles, ${\Theta _i}$, ${\Phi _i}$, and ${\phi _i}$, as mentioned above for making trial movements. This makes it possible to highlight the roles of the coupling between bending and writhing of the DNA model. On the other hand, however, this algorithm can limit the movements of other degrees of freedom of the DNA model. For example, in this algorithm, the first and the second nodal points of the central backbone of the model DNA, ${r_1}$ and ${r_2}$, do not move throughout the Monte Carlo simulations. This means that one of the endpoints of the DNA model is always fixed to the space in the present Monte Carlo simulations. In the following sections, additional degrees of freedom and additional energy functions will be incorporated into the Mote Carlo simulations depending on the problem settings, which will be detailed in each of the following sections.

3. Results and Discussion

3.1. Asymmetric bend-writhe elasticity of the DNA model

We study here the fundamental elasticity of the DNA model with particular attention tothe coupling between bending and writhing of the central backbone. Figure 4(a)-(c) show the conformations of the DNA model with $N = 100$ base pairs, where all dihedral angels of the central backbone are set to ${\Phi _i} = - 4^\circ {\text{}}\left( {i = 2, \cdots , N - 2} \right)$ and all bending angles of the central backbone are set to (a) ${\Theta _i} = 1^\circ$, (b) ${\Theta _i} = 3^\circ$, (c) ${\Theta _i} = 5^\circ$ for $i = 2, \cdots , N - 1$, respectively. Figure 4(d)-(f) show the conformations of the DNA model, where all dihedral angels of the central backbone are set to ${\Phi _i} = + 4^\circ {\text{}}\left( {i = 2, \cdots , N - 2} \right)$ and all bending angles of the central backbone are set to (d) ${\Theta _i} = 1^\circ$, (e) ${\Theta _i} = 3^\circ$, (f) ${\Theta _i} = 5^\circ$ for $i = 2, \cdots , N - 1$, respectively. From these figures, we confirm that negative dihedral angles ${\Phi _i}$ give rise to left-handed super-helical conformations, while positive dihedral angles ${\Phi _i}$ give rise to right-handed super-helical conformations. We also see that the more bending angles ${\Theta _i}$ increase, the more super-helical pitch of the DNA model decreases and super-helical radius increases, under conditions of fixed dihedral angles ${\Phi _i}{\text{}}\left( {i = 2, \cdots , N - 2} \right)$.

Figure 4(g) shows the dependence of the total energy of the DNA model with $N = 100$ base pairs, ${V_{{\text{DNA}}}} = {V_{{\text{bond}}}} + {V_{{\text{bend}}}}$ defined by Eqs. (11) and (12), on writhing of the central backbone. All dihedral angels of the central backbone of the DNA model ${\Phi _i}\left( {i = 2, \cdots , N - 2} \right)$ are set to the same value, ${\Phi }$, and changed together from ${\Phi } = - 4^\circ$ to ${\Phi } = + 4^\circ$, while all bending angles of the central backbone ${\Theta _i}\left( {i = 2, \cdots , N - 1} \right)$ are fixed to ${\Theta _i} = 1^\circ$ (solid curve, red), ${\Theta _i} = 3^\circ$ (dashed-dotted curve, green), and ${\Theta _i} = 5^\circ$ (broken curve, blue) respectively. In obtaining the curves in Figure 4(g), the internal twist angles ${\phi _i}\left( {i = 1, \cdots , N} \right)$ are fixed to the equilibrium values prescribed in Eq. (6). We see anevident asymmetry between left-handed writhing $({\Phi } < 0)$ and right-handed writhing $({\Phi } > 0)$. In particular, left-handed writhing gives less elastic stress than right-handed writhing. This asymmetry becomes even more pronounced as the magnitude of bending angles increases. This result clearly indicates the intrinsic asymmetric coupling between bending and writhing of the DNA model. Specifically, the DNA model has a preference for left-handed writhing upon bending. The asymmetric bend-writhe elasticity observed in Figure 4(g) should be the result of the right-handedness of the double-stranded model of DNA since there is no other chirality in the model or in the potential energy functions of the DNA model.

3.2. Chiral selection in wrapping of the DNA model around a spherical core particle

The asymmetric bend-writhe elasticity observed in Sec. 3.1 indicates that DNA may have a preference on the directionality, i.e., chirality, in wrapping around a core particle. Indeed, in a nucleosome, a repeating unit of chromatin, DNA usually wraps around a protein core particle called histone octamer about 1.75 times in a left-handed manner [1,2]. In this subsection, we try to account for the stability of the left-handed wrapping of DNA around a core particle in terms of the asymmetric elasticity of the DNA model.

We perform Monte Carlo simulations for wrapping of the DNA model around a spherical core particle. In a real nucleosome, positively charged core particle and negatively charged DNA attract each other through electrostatic interactions [10]. For simplicity, we use here the Morse potential,

${V_{{\text{core}}}} = \mathop \sum \limits_{i = 1}^N {D_{{\text{core}}}}\left[{\exp \left\{ { - 2{\beta _{{\text{core}}}}\left( {\left| {{r_i} - {r_{{\text{core}}}}} \right| - {\sigma _{{\text{core}}}}} \right)} \right\} - 2\exp \left\{ { - {\beta _{{\text{core}}}}\left( {\left| {{r_i} - {r_{{\text{core}}}}} \right| - {\sigma _{{\text{core}}}}} \right)} \right\}} \right],$

(17)

for the interaction between DNA and the core, where ${r_i}$ is the position of th nodal point of the central backbone of the DNA model, and ${r_{{\text{core}}}}$ is the position of the center of the spherical core, and ${\sigma _{{\bf{core}}}} = 4.5 {\bf{ nm}}$ is the equilibrium distance between the center of the core and each nodal point of the central backbone of the DNA model. The parameter ${\beta _{{\text{core}}}} = 2.0 {\bf{ n}}{{\bf{m}}^{ - 1}}$ determines the width of the Morse potential. The parameter ${D_{{\text{core}}}} = 10 {\text{ pNnm}}$ determines the strength of the interaction between the core and the DNA model.

To avoid overlaps of the DNA model with itself, we introduce a potential for the excluded volume effect. For this sake, we use the repulsive part of the Morse potential,

${V_{{\text{exc}}}} = \mathop \sum \limits_{i = 1}^{N - n} \mathop \sum \limits_{j = i + n}^N {D_{{\text{exc}}}}\exp \left\{ { - 2{\beta _{{\text{exc}}}}\left( {\left| {{r_i} - {r_{\text{j}}}} \right| - {\sigma _{{\text{exc}}}}} \right)} \right\},$

(18)

where we set the parameters as ${\beta _{{\text{exc}}}} = 2.0 {\text{ n}}{{\text{m}}^{ - 1}}$, ${\sigma _{{\text{exc}}}} = 2.1 {\text{ nm}}$, and ${D_{{\text{exc}}}} = 1.0 {\text{ pNnm}}$. The parameter $n$, which we set $n = 7$, removes the repulsions between very close nodal points in comparison to ${\sigma _{{\text{exc}}}}$.

We carried out Monte Carlo simulations for the composite system consisting of the core particle and the DNA model with $N = 200$ base pairs. The total potential energy function $V$ for the Monte Carlo steps in Eq. (16) is

$V = {V_{{\text{bond}}}} + {V_{{\text{bend}}}} + {V_{{\text{core}}}} + {V_{{\text{exc}}}}.$

(19)

In this Monte Carlo simulation, the position of the first and second nodal points of the DNA model, ${r_1}$ and ${r_2}$, and the position of the core particle, ${r_{{\text{core}}}}$, are fixed in the space.

Figure 5 shows the result of Monte Carlo simulations for wrapping of DNA around a core particle. Figure 5(a) shows a typical evolution of the DNA-core system. The DNA model initially assumes the linear equilibrium conformation and the core particle is located right next to the DNA model at ${r_{{\text{core}}}} = {\left( {0, - 4.5, 10.2} \right)^T}$ in units of [nm]. We see that the DNA model gradually wraps around the core particle in a left-handed manner up to about 2 turns in this simulation. Figure 5(b) shows the evolution of the total energy of the system. We see that the total energy tends to decrease and equilibrate as the DNA model wraps around the core, except for the initial abrupt increase.

Figure 5(c) shows the corresponding evolution of a wrapping number [64] defined by

$W = \frac{{b\left( {{N_{{\text{ad}}}} - 1} \right)}}{{2\pi {\sigma _{{\text{core}}}}}},$

(20)

where ${\text{}}{N_{{\text{ad}}}}$ is the number of nodal points of the central backbone of the DNA model adsorbed on the core. In the Monte Carlo simulations of Figure 5, nodal points satisfying $\left| {{r_i} - {r_{{\text{core}}}}} \right| ≤ 4.8 {\text{ nm}}$ are regarded as adsorbed on the core. We see that the wrapping number $W$ increases as DNA wraps around the core steadily and equilibrates at around $W = 1.75$. In order to judge the chirality of wrapping in a systematic manner, we utilize a chirality parameter [64] defined by

$C = {\langle {m_i}\rangle _{{\text{ad}}}} \cdot \left( {{{\langle {r_i}\rangle }_{{\text{tail}}}} - {{\langle {r_i}\rangle }_{{\text{head}}}}} \right),$

(21)

where ${\langle {m_i}\rangle _{{\text{ad}}}}$ is obtained by averaging the vectors ${m_i} \equiv \left( {{r_i} - {r_{i - 1}}} \right) \times \left( {{r_{i + 1}} - {r_i}} \right)$ over all adsorbed nodal points of the central backbone of the DNA model on the surface of the core particle and normalizing it. ${{{\langle {r_i}\rangle }_{{\text{head}}}}}$ represents the average of position vectors of the first half (with smaller subscripts $i$) of the adsorbed nodal points of the central backbone of DNA, whereas ${{{\langle {r_i}\rangle }_{{\text{tail}}}}}$ represents the average of position vectors of the last half (with larger subscripts $i$) of the adsorbed nodal points of the central backbone. Positive sign of $C$ signifies right-handed wrapping while negative sign of $C$ signifies left-handed wrapping. Figure 5(d) shows the evolution of the chirality parameter $C$ for the Monte Carlo simulation of Figure 5(a)-(c), where $C$ takes mostly negative values, characterizing the left-handed wrapping of the DNA model.

Figure 5(e) shows the probability distribution of the chirality parameter $C$ obtained from 10 runs of Monte Carlo simulations, each of which consists of $6 \times {10^5}$ steps. In order to avoid the bias due to the specific initial conditions introduced above, we have discarded the data of initial $3 \times {10^5}$ steps from each Monte Carlo simulation. The clear peak at around $C = - 0.96$ in Figure 5(e) indicates the strong propensity of the DNA model to wrap around the core in the left-handed manner. The probability distribution of Figure 5(e) lying only in the negative range, $C < 0$, indicates that all the 10 runs of Monte Carlo simulation gave rise to the left-handed wrapping. This result suggests that the asymmetric bend-writhe elasticity of DNA can be the origin of the uniform left-handed wrapping of DNA around nucleosome core particles.

3.3. Chirality of the DNA model adsorbed on a uniform rod

Chiral selectivity of the DNA model observed in Sec. 3.2indicates that DNA may also exhibit chiral selectivity in wrapping around a uniform rod. Indeed, conformation of DNA around rod-like molecules, such as carbon nanotubes [28,29,30,31,32], is an interesting issue. While wrapping of single-stranded DNA around carbon nanotubes is frequently studied, where wrapping geometry is known to be very sensitive to the DNA sequence [28,31], we study here the chirality of the double-stranded DNA model adsorbed on a hypothetical uniform rod. This subsection also serves as a preparatory step before considering crossover and braiding of a pair of DNA molecules in Sec. 3.4.

We fix the uniform rod so that the central axis of the rod coincides with the -axis of the space. We then assume an attractive interaction between DNA and the uniform rod using the Morse potential as

${V_{{\text{rod}}}} = \mathop \sum \limits_{i = 1}^N {D_{{\text{rod}}}}\left[{\exp \left\{ { - 2{\beta _{{\text{rod}}}}\left( {\sqrt {x_i^2 + y_i^2} - {\sigma _{{\text{rod}}}}} \right)} \right\} - 2\exp \left\{ { - {\beta _{{\text{rod}}}}\left( {\sqrt {x_i^2 + y_i^2} - {\sigma _{{\text{rod}}}}} \right)} \right\}} \right],$

(22)

where ${x_i}$ and ${y_i}$ represent $x$- and y-components of the position vector of the central backbone of the DNA model, ${r_i}$. The parameter ${\sigma _{{\text{rod}}}}$ determines the effective radius of the rod including the radius of DNA. In the following, we examine three different values of ${\sigma _{{\text{rod}}}} = 2.1 {\text{ nm}}$, ${\sigma _{{\text{rod}}}} = 3.1 {\text{ nm}}$, and ${\sigma _{{\text{rod}}}} = 4.1 {\text{ nm}}$. The parameter ${D_{{\text{rod}}}}$ determines the strength of the interaction between the rod and DNA, which we set ${D_{{\text{rod}}}} = 20 {\text{ pNnm}}$. The parameter ${\beta _{{\text{rod}}}} = 2.0 {\text{ n}}{{\text{m}}^{ - 1}}$ determines the width of the Morse potential. Here, we use the DNA model with $N = 100$ base pairs. We also introduce the excluded volume effect in the DNA model using exactly the same potential function ${V_{{\text{exc}}}}$ as in Eq. (18).

We perform Monte Carlo simulations for the system of the rod and the DNA model molecule. Here, the total potential energy function for the Monte Carlo steps in Eq. (16) is

$V = {V_{{\text{bond}}}} + {V_{{\text{bend}}}} + {V_{{\text{rod}}}} + {V_{{\text{exc}}}}.$

(23)

The first nodal point of the central backbone of the DNA model is fixed to the space at ${r_1} = {\left( {{\sigma _{{\text{rod}}}}, 0, 0} \right)^T}$ throughout the Monte Carlo simulation in this subsection. In this Monte Carlo simulation, the DNA model initially assumes thestraight equilibrium conformation and is alignedparallel to the uniform rod in the direction of positive $z$-axis. Initial separation distance between the central backbone of the DNA model and that of the rod is ${\sigma _{{\text{rod}}}}$. This initial configuration of the system corresponds to the minima of the DNA-rod interaction potential, ${V_{{\text{rod}}}}$ in Eq. (22). Therefore, there is no external torque that drives the DNA model to wrap around the uniform rod in the Monte Carlo simulation. However, because of the thermal fluctuations, the DNA model can still deform from the initial conformation to take superhelical or undulating conformations around the rod.

In order to characterize chirality and the number of wrapping turns of the DNA model around a uniform rod, we introduce a wrapping number as

${W_{{\text{rod}}}} = \frac{1}{{2\pi }}\mathop \sum \limits_{i = 2}^{N - 1} {\text{sgn}}\left[{\hat z \cdot {n_i}} \right]{\cos ^{ - 1}}\left( {\frac{{{x_i}{x_{i + 1}} + {y_i}{y_{i + 1}}}}{{\sqrt {x_i^2 + y_i^2} \sqrt {x_{i + 1}^2 + y_{i + 1}^2} }}} \right),$

(24)

where $\hat z \equiv {\left( {0, 0, 1} \right)^T}$ is a unit vector along $z$-axis, ${n_i}$ is a unit vector perpendicular to the central backbone of the DNA model atth nodal point defined by ${n_i} \equiv \left( {{r_i} - {r_{i - 1}}} \right) \times \left( {{r_{i + 1}} - {r_i}} \right)/\left| {\left( {{r_i} - {r_{i - 1}}} \right) \times \left( {{r_{i + 1}} - {r_i}} \right)} \right|$, and ${\text{sgn}}\left[\cdot \right]$ represents the sign of ${W_{{\text{rod}}}}$ its argument. Positive sign of signifies right-handed wrapping of the DNA model around the rod while negative sign of ${W_{{\text{rod}}}}$ signifies left-handed wrapping. Absolute value of ${W_{{\text{rod}}}}$ can be regarded as a measure of the number of wrapping turns of the DNA model around the rod.

Figure 6 shows the probability distributions of the wrapping number ${W_{{\text{rod}}}}$, each of which is obtained from 20 runs of Monte Carlo simulations of the DNA model adsorbed on a uniform rod. Each Monte Carlo run consists of $2 \times {10^6}$ steps, where the data of the first $2 \times {10^5}$ steps are discarded to obtain the probability distributions in Figure 6 in order to avoid the effects of the initial conditions. The rod radius parameter is ${\sigma _{{\text{rod}}}}$ (a) ${\sigma _{{\text{rod}}}} = 2.1 {\text{ nm}}$, (b) ${\sigma _{{\text{rod}}}} = 3.1 {\text{ nm}}$, and (c) ${\sigma _{{\text{rod}}}} = 4.1 {\text{ nm}}$, respectively. From Figure 6, we see that the probability distributions of the wrapping number ${W_{{\text{rod}}}}$ are mostly symmetric for all the three values of the rod radius parameter ${\sigma _{{\text{rod}}}}$. This indicates that the conformation of the DNA model fluctuates almost equally for left and right directions around the straight equilibrium conformation on the surface of the rod. We see that the absolute value of the wrapping number $\left| {{W_{{\text{rod}}}}} \right|$ tends to increase as the radius of the rod ${\sigma _{{\text{rod}}}}$ decreases. Since the absolute value of the wrapping number is small (mostly $\left| {{W_{{\text{rod}}}}} \right| < 0.1$) for all the three values of the rod radius parameter ${\sigma _{{\text{rod}}}}$, we see that the DNA model seldom achieves a clear wrapping around a rod. Above each graph of Figure 6, typical configuration of the model system corresponding to the peak of the probability distribution is shown with the corresponding wrapping number ${W_{{\text{rod}}}}$.

The interaction potential introduced in Eq. (22) may correspond to the situation where the attraction between DNA and the rod is “local”, limited only to the surfaces of DNA and the rod. In this situation, the parallel juxtaposition of DNA and the rod corresponds to the minimum of the total potential energy. This is why the DNA model seldom achieves a clear wrapping around the rod, and the probability distributions in Figure 6 are mostly symmetric around ${W_{{\text{rod}}}} = 0$. It should be noted however that another parameter that we examined, other than the wrapping number ${W_{{\text{rod}}}}$, has indicated a bias towards left-handed wrapping for the same Monte Carlo simulations as in Figure 6. This indicates that the DNA model still exhibits a weak tendency for left-handed wrapping around the uniform rod in thermal fluctuations, but this tendency may be too weak to be detected by the wrapping number ${W_{{\text{rod}}}}$.

On the other hand, when the attraction between DNA and a rod is “global”, in the sense that a part of the rod attracts many different parts of DNA and a part of DNA attracts many different parts the rod, wrapped states of DNA around the rod could be stabilized at the expense of the bending energy of DNA. Based on this consideration, we also examine a global interaction between the DNA model and a uniform rod by using Morse potential as

${V_{{\text{rod}}}} = \mathop \sum \limits_{i = 1}^N \mathop \sum \limits_{j = 1}^N {D_{{\text{rod}}}}\left[{\exp \left\{ { - 2{\beta _{{\text{rod}}}}\left( {\left| {{r_i} - r_j^{{\text{rod}}}} \right| - {\sigma _{{\text{rod}}}}} \right)} \right\} - 2\exp \left\{ { - {\beta _{{\text{rod}}}}\left( {\left| {{r_i} - r_j^{{\text{rod}}}} \right| - {\sigma _{{\text{rod}}}}} \right)} \right\}} \right],$

(25)

where ${r_i}$ represents the position of th nodal point of the central backbone of the model DNA. Here, we have introduced $N = 100$ nodal points on the central axis of the rod ($z$-axis), $r_j^{{\text{rod}}}\left( {j = 1, \cdots , N} \right)$, with equal separation of $b = 0.34 {\text{nm}}$, which is the same as the separation between two adjacent nodal points of the central backbone of the DNA model. Thus the rod is regarded to have the same length as the DNA model, $b\left( {N - 1} \right)$, and to lie over the range, $0 ≤ z ≤ b\left( {N - 1} \right)$. As we see from Eq. (25), each nodal point of the DNA model interacts with all nodal points of the rod and vice versa. For this “global” potential of Eq. (25), we set the parameters as ${D_{{\text{rod}}}} = 5 {\text{ pNnm}}$ and ${\beta _{{\text{rod}}}} = 2.0 {\text{ n}}{{\text{m}}^{ - 1}}$.

We perform Monte Carlo simulations in the same manner as in Figure 6, after replacing the interaction potential between the rod and the DNA model in Eq. (22) with the “global” potential in Eq. (25). Figure 7 shows the results of probability distributions of the wrapping number ${W_{{\text{rod}}}}$ for the “global” potential in Eq. (25), each of which is obtained from 20 runs of Monte Carlo simulations. Each Monte Carlo run consists of $2 \times {10^6}$ steps, where the data of the first $2 \times {10^5}$ steps are discarded.The rod radius parameter ${\sigma _{{\text{rod}}}}$ is (a) ${\sigma _{{\text{rod}}}} = 2.1 {\text{ nm}}$, (b) ${\sigma _{{\text{rod}}}} = 3.1 {\text{ nm}}$, and (c) ${\sigma _{{\text{rod}}}} = 4.1 {\text{ nm}}$, respectively. In each panel of Figure 7, the probability distribution shows two major peaks; one is in the negative range, ${W_{{\text{rod}}}} < 0$, corresponding to left-handed wrapping, and the other is in the positive range, ${W_{{\text{rod}}}} > 0$, corresponding to right-handed wrapping. Above each graph of Figure 7, typical configurations of the model system corresponding to the two major peaks are shown with the corresponding wrapping number ${W_{{\text{rod}}}}$. Existence of these two peaks indicates that wrapped states are stabilized at the expense of bending energy of the model DNA for the case of the “global” potential of Eq. (25). Moreover, absolute value of wrapping number $\left| {{W_{{\text{rod}}}}} \right|$ for the “global” potential in Figure 7 tends to be large as compared to that for “local” potential in Figure 6. These are significant differences between Figure 6 and Figure 7. Importantly, the peaks in the range ${W_{{\text{rod}}}} < 0$ are always higher than those in the range ${W_{{\text{rod}}}} > 0$ for all the three values of rod radius parameter in Figure 7, indicating the clear preference for left-handed wrapping. In the case of “global” potential ${\sigma _{{\text{rod}}}}$ in Figure 7, the DNA model is forced to bend more largely than in Figure 6. Such strong bending in turn provokes the intrinsic propensity of the DNA model to writhe in the left direction through the asymmetric coupling between bending and writhing as revealed in Figure 4. In this manner, the DNA model tends to select left-handed wrapping for the case of “global” interaction even though there is no external torque. Based on this result, we will investigate crossover and braiding of two DNA molecules interacting through a “global” potential in the next subsection.

3.4. Chiral selection in crossover and braiding of two juxtaposed DNA molecules

The results of Sec. 3.2 and Sec. 3.3 indicate that DNA may have a preferred chirality for crossover and braiding due to its asymmetric elasticity. In fact, chirality in crossover and braiding of juxtaposed DNA molecules has been an important issue. As noted in the Introduction, recent works of Timsit and Várnai [36,37,38] and those of Kornyshev et al. [40,41,42,43,44] indicate that right-handed crossover and left-handed braiding are preferable for two juxtaposed DNA molecules in solution in the presence of cations. It should be noted here that right-handed crossover and left-handed braiding are intimately related (or similar) phenomena since right-handed crossovers occur in a left-handed braiding of DNA molecules [36,37,38,39,40,41,42,43,44]. While these works have highlighted the electrostatic effects and steric effects, we investigate here the roles of asymmetric elasticity of DNA in braiding and crossover based on our coarse-grained model.

DNA molecules can attract each other effectively in solution depending on salt condition. Suppose there are two DNA molecules, A and B, attracting each other. We model theattractive interaction and the excluded volume effect between these DNA molecules (DNA-DNA interaction) at the same time using the Morse potential for simplicity as

$V_{{\text{ext}}}^{{\text{AB}}} = \mathop \sum \limits_{i = 1}^N \mathop \sum \limits_{j = 1}^N {D_{{\text{DNA}}}}\left[{\exp \left\{ { - 2{\beta _{{\text{DNA}}}}\left( {\left| {r_i^A - r_j^{\text{B}}} \right| - {\sigma _{{\text{DNA}}}}} \right)} \right\} - 2\exp \left\{ { - {\beta _{{\text{DNA}}}}\left( {\left| {r_i^{\text{A}} - r_j^{\text{B}}} \right| - {\sigma _{{\text{DNA}}}}} \right)} \right\}} \right],$

(26)

where $r_i^{\text{A}}$ and $r_i^{\text{B}}$ represent the positions of i-th and j-th nodal points of the central backbones of the two DNA models, A and B. The parameter ${\sigma _{{\text{DNA}}}}$ represents the equilibrium distance between the nodal points of the central backbones of the two DNA molecules. We set ${\sigma _{{\text{DNA}}}} = 2.1 {\text{ nm}}$, which is roughly the sum of the radii of the two DNA molecules. The parameter ${\beta _{{\text{DNA}}}} = 0.8 {\text{ n}}{{\text{m}}^{ - 1}}$ determines the width of the Morse potential. The parameter ${D_{{\text{DNA}}}}$ determines the strength of the DNA-DNA interaction per each pair of nodal points and serves as the control parameter in the present study. Later, we change the value of ${D_{{\text{DNA}}}}$ from ${D_{{\text{DNA}}}} = 0.5 {\text{ pNnm}}$ to ${D_{{\text{DNA}}}} = 1.5 {\text{ pNnm}}$. It should be noted that the interaction potential in Eq. (26) is “global” similarly to Eq. (25) in the sense that a part of one DNA model interacts with many different parts of the other DNA model and vice versa, as is expected for the DNA-DNA interaction in solution mediated by salt.

When two DNA molecules attract each other due to salt condition, internal attraction should also exist inside each DNA molecule. Therefore we model thisattractive interaction as well as the excluded volume effect in each DNA molecule using the Morse potential as

$V_{{\text{int}}}^{\text{X}} = \mathop \sum \limits_{i = 1}^{N - n} \mathop \sum \limits_{j = i + n}^N {D_{{\text{DNA}}}}\left[{\exp \left\{ { - 2{\beta _{{\text{DNA}}}}\left( {\left| {r_i^X - r_j^X} \right| - {\sigma _{{\text{DNA}}}}} \right)} \right\} - 2\exp \left\{ { - {\beta _{{\text{DNA}}}}\left( {\left| {r_i^X - r_j^X} \right| - {\sigma _{{\text{DNA}}}}} \right)} \right\}} \right],$

(27)

where the superscript X distinguishes the two DNA models, i.e., $X = A {\text{or}} B$, and $r_i^X$ represents the position of i-th nodal point of central backbone of the DNA model X. We use the same values for the parameters ${\sigma _{{\text{DNA}}}}, {\beta _{{\text{DNA}}}}, {\text{and}} {D_{{\text{DNA}}}}$ as in the DNA-DNA interaction potential in Eq. (26). We set $n = 7$ to remove the repulsions between very close nodal points in comparison to ${\sigma _{{\text{DNA}}}}$. Note that the interaction potentials in Eqs. (26) and (27) do not break the symmetry between left- and right-handed crossovers and braidings of the two DNA models.

In order to characterize chirality and degree of braiding of two DNA molecules, we introduce a braiding number as

$B = \frac{1}{{2\pi }}\mathop \sum \limits_{i = 1}^{N - 1} {\text{sgn}}\left[{\left( {{C_{i + 1}} - {C_i}} \right) \cdot \left( {{B_i} \times {B_{i + 1}}} \right)} \right]{\sin ^{ - 1}}\left| {{B_i} \times {B_{i + 1}}} \right|,$

(28)

where

${C_i} = \frac{{r_i^{\text{A}} + r_i^{\text{B}}}}{2}, \quad {\text{}}{B_i} = {\text{}}\frac{{r_i^{\text{A}} - r_i^{\text{B}}}}{{\left| {r_i^{\text{A}} - r_i^{\text{B}}} \right|}},$

(29)

and ${\text{sgn}}\left[\cdot \right]$ represents the sign of its argument. Because of the definition of Eq. (28), positive values of the braiding number ($B > 0$) signify right-handed braiding between two DNA molecules and negative values of the braiding number ($B < 0$) signify left-handed braiding. Absolute value of the braiding number $\left| B \right|$ is a measure of the number of braiding turns of two DNA molecules.

We perform Monte Carlo simulations for two juxtaposed DNA model molecules, each of which consists of $N = 100$ base pairs. The total potential energy for the Monte Carlo steps in Eq. (16) is

$V = {V_{{\text{bond}}}} + {V_{{\text{bend}}}} + V_{{\text{ext}}}^{{\text{AB}}} + V_{{\text{int}}}^{\text{A}} + V_{{\text{int}}}^{\text{B}}.$

(30)

In this Monte Carlo simulation, the two DNA molecules initially assume thestraight equilibrium conformation and are aligned parallel to the $z$-direction of the space with the separation distance of . In this Monte Carlo simulation, we randomly choose one of the bending angles ${\Theta _i}$, dihedral angels ${\Phi _i}$, and internal twist angles ${\phi _i}$ of one of the two DNA model molecules, A and B, and perturb it by $+ 1^\circ$ or $- 1^\circ$ randomly to obtain a trial conformation. The first nodal points of the central backbones of the two DNA models are fixed to the space at $r_1^{\text{A}} = {\left( {0, 0, 0} \right)^T}$ and $r_1^{\text{B}} = {\left( {4, 0, 0} \right)^T}$ in units of [nm] throughout the Monte Carlo simulations.

Figure 8(a) shows an evolution of the braiding number $B$ defined in Eq. (28) in a Monte Carlo simulation of the two juxtaposed DNA models, where the strength of DNA-DNA interaction is set as ${D_{{\text{DNA}}}} = 0.7 {\text{ pNnm}}$. The upper pictures of a pair of DNA molecules in Figure 8(a) show the overall configuration of the system at the indicated Monte Carlo steps. From these pictures, we see that the two DNA model molecules tend to braid each other spontaneously even in the absence of external torque. This spontaneous braiding indicates that the braided conformations are stabilized because of the global attractive interaction between the two DNA molecules. The graph of Figure 8(a) confirms that the braiding number properly characterizes the chirality of braiding. In the simulation of Figure 8(a), the braiding number mostly takes negative values (see the conformation of the system at 1152000th step), indicating that the left-handed braiding is preferred in this simulation. In this simulation, the braiding number $B$ approaches to zero at around 380000th step, whose corresponding configuration is also shown in Figure 8(a).

Because of the global attractions between and within the DNA model molecules, the two juxtaposed DNA model molecules can sometimes collapse and form a globule, which complicates the analysis of braiding using the braiding number $B$. Therefore, we also performed Monte Carlo simulations under tension in the $z$-direction of the space. In order to incorporate the tension into the Monte Carlo simulation, we introduce the following potential functions for tension in the $z$ direction as

$V_{{\text{tension}}}^A = - Fz_N^A, \qquad V_{{\text{tension}}}^B = - Fz_N^B,$

(31)

into the right-hand side of Eq. (30), where $F$ is the strength of tension, and $z_N^A$ and $z_N^B$ are the $z$-coordinates of the $N$-th nodal points (terminal nodal points) of the central backbones of the two DNA models, A and B, respectively. Figure 8(b) and (c) show the results of Monte Carlo simulations under tensions of $F = 60 {\text{ pN}}$ and $F = 100 {\text{ pN}}$ respectively. From the results of Figure 8(b) and (c), we see that the absolute values of the braiding number tend to decrease as the tension increases. This indicates that the two DNA model molecules do not braid largely but just crossoveras shown in the figures. However, we can still see that the braiding number $B$ tends to take negative values, indicating that the two DNA model molecules tend to prefer left-handed braiding. In the present study, configurations with the absolute value of the braiding number greater than 1.0 may be regarded as braiding ($\left| B \right| > 1$), while configurations with the absolute value of the braiding number less than 1.0 may be regarded as crossover ($\left| B \right| < 1$). In the terminology of crossover [36], left-handed braiding corresponds to right-handed crossover. Thus, the results of Figure 8(b) and (c) indicate that the two DNA molecules show the tendency to right-handed crossover.

To characterize chiral selectivity in the braiding of a pair of DNA model molecules statistically, Figure 9 shows the probability distributions for the braiding number of two juxtaposed DNA molecules, $B$, for different magnitudes of tension, $F$. Each probability distribution in Figure 9 is obtained from 20 runs of Monte Carlo simulations, each of which consists of $2 \times {10^6}$ steps. In order to avoid bias due to the specific initial conditions, we have discarded the data of initial $2 \times {10^5}$ steps from each Monte Carlo simulation. Integrated probabilities for left-handed ($B < 0$) and right-handed ($B > 0$) braiding are shown inside each panel of Figure 9. Above the graphs of Figure 9, typical configurations of a pair of DNA molecules are shown for the braiding numbers that roughly correspond to the prominent peaks of the probability distributions. Figure 9(a) shows the results for no tension, where we see that the probability is distributed over a wide range of the braiding number $B$. Above the graph of Figure 9(a), two typical configurations of DNA molecules are shown for the braiding numbers, $B = - 1.2$ and $B = + 1.38$. From Figure 9(a), we clearly observe that the braiding number $B$ shows higher probability in the negative range than in the positive range. This indicates that the two DNA molecules have a tendency to braid in a left-handed manner.

Figure 9(b) shows the probability distribution of the braiding number $B$ under tension of $F = 60 {\text{ pN}}$. The probability distribution of the braiding number in Figure 9(b) is narrower than that of Figure 9(a). This indicates that the tendency for braiding is reduced because of the tension. However, we can still see that the distribution in Figure 9(b) is biased towards negative braiding numbers centered around $B = - 0.11$. Figure 9(c) shows the probability distribution of the braiding number for even stronger tension, $F = 100 {\text{ pN}}$. We see that the probability distribution of the braiding number becomes even narrower and symmetric as compared to Figure 9(b) because of the even stronger tension. However, we can still see that the distribution in Figure 9(c) is biased towards to negative braiding numbers centered around $B = - 0.04$.

Figure 10 summarizes the dependence of the probabilities for left-handed braiding and that for right-handed braiding on the magnitude of tension, $F$. The probabilities for left-handed and right-handed braiding plotted in Figure 10 are defined as the integrated areas of the probability distributions as in Figure 9 over the ranges $B < 0$ and $B > 0$ respectively. We see that the chiral selectivity is the highest when there is no tension, $F = 0 {\text{ pN}}$, and it decreases as the tension increases, which is consistent with the observations in Figure 9. It should be noted that the collapse of the DNA molecules occurred twice in the 20 runs of Monte Carlo simulation with no tension, $F = 0 {\text{ pN}}$. Therefore, the data of braiding numbers of such collapsed conformations are included in the probabilities in Figure 9(a) and Figure 10 for $= 0 {\text{ pN}}$. For other values of tension, $F \ne 0{\text{ pN}}$, we did not observe the collapse.

3.5. Dependence of braiding chirality on the strength of DNA-DNA interaction

Dependence of braiding behaviors on the strength of the interaction between DNA molecules is an interesting issue. This is because the effective interaction forces between DNA molecules can change significantly in solution depending on salt condition. We thus study here the dependence of chiral selectivity in DNA braiding on the Morse parameter ${D_{{\text{DNA}}}}$ in Eqs. (26) and (27), which determines the strength of the interaction between and inside the two DNA model molecules. The procedures for the Monte Carlo simulations in this subsection are the same as in the previous subsection, Sec. 3.4. To avoid collapse of the DNA model molecules, we introduce a tension of $F = 40{\text{ pN}}$ in the $z$-direction of the space.

Figure 11 shows typical evolutions of the system in the Monte Carlo simulations, where the strengths of the DNA-DNA interaction are (a) ${D_{{\text{DNA}}}} = 0.5 {\text{ pNnm}}$, (b) ${D_{{\text{DNA}}}} = 1.0 {\text{ pNnm}}$, and (c) ${D_{{\text{DNA}}}} = 1.5 {\text{ pNnm}}$, respectively. When the interaction between DNA molecules is weaker as in Figure 11(a), the braiding number $B$ switches between positive and negative values during a single Monte Carlo simulation, indicating that the two DNA model molecules switch between left-handed and right-handed braiding.

When the interaction between DNA molecules becomes stronger as in Figure 11(b), frequency of switching between left-handed and right-handed braiding is reduced, but the absolute value of the braiding number can become larger. In the simulation of Figure 11(b), we also see the tendency that the braiding number tends to take negative values more frequently than positive values, indicating the preference for left-handed braiding. When the attractive interaction between DNA molecules becomes even stronger as in Figure 11(c), the switching between left-handed and right-handed braiding hardly occurs during the single Monte Carlo simulation. In the simulation of Figure 11(c), the two DNA molecules appear to have settled into left-handed braiding.

To study the dependence of chiral selectivity in the braiding of two DNA molecules on the strength of DNA-DNA interaction, we show in Figure 12 the probability distributions of the braiding number $B$ fordifferentstrengths of the DNA-DNA interaction; (a) ${D_{{\text{DNA}}}} = 0.5 {\text{ pNnm}}$, (b) ${D_{{\text{DNA}}}} = 1.0 {\text{ pNnm}}$, and (c) ${D_{{\text{DNA}}}} = 1.5 {\text{ pNnm}}$. Each probability distribution in Figure 12 is obtained from 20 runs of Monte Carlo simulation, each of which consists of $2 \times {10^6}$ steps. In order to avoid bias due to the specific initial conditions, we have discarded the data of initial $2 \times {10^5}$ steps from each Monte Carlo simulation. Above each graph of Figure 12, typical configurations of the system are shown, which have the braiding numbers roughly corresponding to the prominent peaks of the probability distributions. When the attractive interactions between and inside DNA models are weak, as in Figure 12(a), the distribution of the braiding number is narrow and single-peaked, indicating that the DNA models do not braid largely but they just crossover. However, we still see from Figure 12(a) that the distribution of the braiding number is biased towards the negative direction centered around $B = - 0.09$. This indicates that the DNA models provoke the tendency to left-handed braiding.

When the attractive interaction becomes stronger as in Figure 12(b), the probability distribution of the braiding number becomes wider and is still biased towards the negative range possessing a faint second peak at around $B = - 0.55$ in addition to the primary peak at around $B = - 0.15$. The second peak at around $B = - 0.55$ corresponds to a deeper crossover configuration. The emergence of such state of deeper crossover is due to the tendency of the two DNA model molecules to increase the contacts at the expense of bending energy, which is in turn due to the strong attraction. When the attractive interaction becomes even stronger as in Figure 12(c), the probability distribution of the braiding number becomes even wider in both positive and negative regions. Such widening of distribution indicates that the two DNA molecules can crossover even deeper and finally braid each other. The distribution of Figure 12(c) possesses several peaks, indicating a stepwise braiding. In the negative range, $B < 0$, the probability distribution is biased towards the peak that has larger absolute value of the braiding number $B$, while in the positive range $B > 0$, the probability distribution is biased towards the peak that has smaller absolute value of the braiding number $B$. This indicates that left-handed braiding tends to have larger number of braiding turns, while right-handed braiding tends to have smaller number of braiding turns. In this manner, the probability distribution of the braiding number is asymmetric. However, the total integrated probability in the negative range $B < 0$ and that in the positive range $B > 0$ are similar.

Figure 13 shows the dependence of the probability for left-handed braiding and that for right-handed braiding on the strength of the DNA-DNA interaction, ${D_{{\text{DNA}}}}$. Here the probabilities for left-handed and right-handed braidings are essentially the integrated areas of the probability distributions as in Figure 12 over the ranges $B < 0$ and $B > 0$ respectively. We see that chiral selectivity is the most pronounced at around intermediate strength of the DNA-DNA interaction, ${D_{{\text{DNA}}}} = 1.1 {\text{ pNnm}}$. When the DNA-DNA interaction becomes weaker than ${D_{{\text{DNA}}}} = 1.1{\text{ pNnm}}$, chiral selectivity becomes slightly weaker. This may be because of the dominance of the tension applied on the two DNA models over the DNA-DNA interaction. When the DNA-DNA interaction becomes stronger than ${D_{{\text{DNA}}}} = 1.1{\text{ pNnm}}$, probabilities for left-handed and right-handed braidings approach to a similar value ($50\%$). In this sense, chiral selectivity disappears for the strong DNA-DNA interaction. However, it should be noted that, even for such strong DNA-DNA interaction, probability distribution of the braiding number is still asymmetric in terms of left-handed and right-handed braiding, as we have seen in Figure 12(c).

4. Conclusion

We have investigated the elastic mechanisms for the selection of chirality in wrapping, crossover, and braiding of DNA based on a coarse-grained model. The DNA model has consisted of two elastic chains that mutually intertwine in a right-handed manner forming a double-stranded helix with the distinction between major and minor grooves. Regardless of the scheme of modeling of DNA that local potential energy functions have no chirality, the model has exhibited an asymmetric propensity as a whole to writhe in the left direction upon bending due to the right-handed helical geometry. The result of Monte Carlo simulations of this model has indicated that DNA has a preference for the left-handed wrapping around a spherical core particle and also around a uniform rod due to the asymmetric elastic coupling between bending and writhing. This result has suggested an elastic origin of the uniform left-handed wrapping of DNA in nucleosomes and also has provided implications on the wrapping of double-stranded DNA around rod-like molecules. Monte Carlo simulations of the DNA model have also suggested that two juxtaposed DNA molecules can braid each other spontaneously under moderate attractive interactions with the preference for left-handed braiding due to the asymmetric coupling between bending and writhing. This result thus indicates the importance of asymmetric elasticity to cause the spontaneous left-handed braiding of DNA molecules.

The present study has highlighted almost exclusively the roles of asymmetric bend-writhe elasticity of DNA in the formation of higher-order structures. Therefore, the present DNA model has disregarded many other important characteristics of real DNA, which may be taken into consideration in the future studies for reality. For one thing, the Monte Carlo simulations of the present study have allowed the movements of only limited degrees of freedom of the coarse-grained DNA model, such as the bending angles and dihedral angels of the central backbone and the internal twist anglesIn reality, however, many other internal degrees of freedom of DNA could also couple each other giving rise to more complex behaviors. In addition, in the present model of DNA, we have determined the spring constants for bonding ${k_{{\text{bond}}}}$ and for bending ${k_{{\text{bend}}}}$ so that our DNA model reproduces approximately the same bending persistence length as real DNA. Use of additional data such as stretching and twisting rigidities of real DNA would improve the evaluation of the spring constants of the present model. However, we note that the asymmetric elasticity of the DNA model to writhe in the left direction upon bending has been almost independent of the setting of these spring constants qualitatively.

For another thing, the present study has used the Morse potential frequently to avoid the complications of electrostatic effects. However, since electrostatic effects play essential roles in the interactions of DNA with nucleosome core particles [23,24], carbon nanotubes [28,29,30,31,32], and other DNA molecules [36,37,38,39,40,41,42,43,44], it will be important to take into consideration electrostatic effects in more rigorous manners in the future study. Indeed, the asymmetric bend-writhe elasticity suggested by our present DNA model may play a complementary role to the reported electrostatic effects due to the helical alignments of charges and image charges [44,57,58] in the selection of chirality in wrapping, crossover and braiding of DNA. In nature, combination of the elastic effects and the electrostatics and steric effects based on atomistic level of interactions [35,36,37,38,39,40,41,42,43,44,45,46] would be important for proper chiral selection of DNA.

It is also an interesting issue to extend the implications of the present study to other biopolymers. Since helical geometry is ubiquitous among biopolymers such as actin filaments [3,4], microtubules [5,6], and collagen [7,8], etc., the asymmetric elastic properties of our simplified model may be applicable to a wide class of biopolymers, especially to grasp the essentials of giant macromolecules through coarse-graining.

Acknowledgments

We thank Prof. A. A. Kornyshev for stimulating discussions. Graphics of the DNA model in the present paper have been created by VMD [68]. This study has been partially supported by JSPS Grant-in-Aid (No. 26800207), and Grants-in-Aid, KAKENHI (25103012; 15H02121).

Conflict of Interest

The authors certify that there is no conflict of interest.