Efficient simulation of semiflexible polymers

Using a recently developed bead-spring model for semiflexible polymers that takes into account their natural extensibility, we report an efficient algorithm to simulate the dynamics for polymers like double-stranded DNA (dsDNA) in the absence of hydrodynamic interactions. The dsDNA is modeled with one bead-spring element per base pair, and the polymer dynamics is described by the Langevin equation. The key to efficiency is that we describe the equations of motion for the polymer in terms of the amplitudes of the polymer's fluctuation modes, as opposed to the use of the physical positions of the beads. We show that, within an accuracy tolerance level of $5\%$ of several key observables, the model allows for single Langevin time steps of $\approx 1.6$, 8, 16, and 16 ps for a dsDNA model chain consisting of 64, 128, 256, and 512 base pairs (i.e., chains of 0.55, 1.11, 2.24, and 4.48 persistence lengths), respectively. Correspondingly, in 1 h, a standard desktop computer can simulate 0.23, 0.56, 0.56, and 0.26 ms of these dsDNA chains, respectively. We compare our results to those obtained from other methods, in particular, the (inextensible discretized) wormlike chain (WLC) model. Importantly, we demonstrate that at the same level of discretization, i.e., when each discretization element is one base pair long, our algorithm gains about five to six orders of magnitude in the size of time steps over the inextensible WLC model. Further, we show that our model can be mapped one on one to a discretized version of the extensible WLC model, implying that the speed-up we achieve in our model must hold equally well for the latter. We also demonstrate the use of the method by simulating efficiently the tumbling behavior of a dsDNA segment in a shear flow.

Figure 1

The ratio ${[\langle H_p^2\rangle /\langle \mathrm{\Delta }{H}_p^2\rangle ]}^{1/2}$ for a dsDNA chain of length $N=63$.

Figure 2

Determination of $\mathrm{\Delta }{\tau }_{\text{max}}$ from the equilibrium values $Q\left(0\right)$. The autocorrelation functions in time of the end-to-end vector are shown in black circles, the middle bond is shown in blue squares, and the mean-square displacement (msd) of the middle bead with respect to the position of the center of mass of the chain is shown in red diamonds. The data are rescaled by their MC values: (a) $N=63$, (b) $N=127$, (c) $N=255$, and (d) $N=511$. The error bars for the autocorrelation functions in time of the middle bond are not shown since they are smaller than the symbol size. The yellow band represents $5\%$ validity thresholds. See text for details.

Figure 3

Determination of $\mathrm{\Delta }{\tau }_{\text{max}}$ from $Q\left(t\right)$. The autocorrelation functions in time of the end-to-end vector are shown in black circles, the middle bond is shown in blue squares, and the mean-square displacement (msd) of the middle bead with respect to the position of the center of mass of the chain is shown in red diamonds. The data are rescaled by $\tilde{Q}\left(\mathrm{\Delta }{\tau }_{\text{min}}\right)$. (a) $N=63$ ($\mathrm{\Delta }{\tau }_{\text{min}}=0.1$): $\tau =100$ (solid line), 1000 (long-dashed), 10 000 (dash-dot-dashed), 100 000 (dash-dot-dot-dashed). (b) $N=127$, ($\mathrm{\Delta }{\tau }_{\text{min}}=1$): $\tau =1000$ (solid line), 10 000 (long-dashed), 100 000 (dash-dot-dashed), $1\times {10}^6$ (dash-dot-dot-dashed). (c) $N=255$ ($\mathrm{\Delta }{\tau }_{\text{min}}=1$): $\tau =10000$ (solid line), 100 000 (long-dashed), $1\times {10}^6$ (dashed-dot-dashed), $1\times {10}^7$ (dash-dot-dot-dashed). (d) $N=511$ ($\mathrm{\Delta }{\tau }_{\text{min}}=5$): $\tau =100000$ (solid line), $1\times {10}^6$ (long-dashed), $1\times {10}^7$ (dash-dot-dashed), $1\times {10}^8$ (dash-dot-dot-dashed). The yellow band represents $10\%$ validity thresholds. See text for details.

Figure 4

The $Q\left(t\right)$ curves for $\mathrm{\Delta }\tau =\mathrm{\Delta }{\tau }_{\text{min}}$ to $\mathrm{\Delta }{\tau }_{\text{max}}$ for $N=255$: the normalized autocorrelation function of (a) the end-to-end vector $\mathbf{L}$ and (b) the middle bond vector ${\mathbf{u}}_m$ and (c) the mean-square displacement $\langle \mathrm{\Delta }{\tilde{r}}_m^2\left(t\right)\rangle$ of the middle bead measured with respect to the center of mass of the chain. The data for all values of $\mathrm{\Delta }\tau$ coincide, as they should. The figure thus demonstrates the validity of our procedure.

Figure 5

The end-to-end vector data for a chain of 256 base pairs and that of the corresponding chain of 128 dimerized base pairs.