Polymer-looping-under-tension-WLC

Code and data used for publication in Physical Review Research: "Theoretical models for tension-dependent DNA looping time", by Wout Laeremans and Wouter G. Ellenbroek.

LAMMPS

The main LAMMPS file is main.lmp, which contains the code for a semiflexible chain polymer under tension and outputs the timestep at which a loop is formed (i.e., when the two outermost beads come within a capture distance $r_\mathrm{c}$). In the uploaded example, the polymer consists of 50 beads (so $N = 49$ is the number of bonds), and the corresponding data file is polymer50.data. The bond length $b$ is set to 1.0 and can be changed in line 15. The force $f$ is set to 1.0 as well and can be changed in lines 31-32. This file can be easily adapted to simulate a different number of beads, as well as other parameters used in the example.

The files rxn_pre.data_template, rxn_post.data_template, and rxn.map are included to determine when a loop is formed. For this, we use the REACTER package and define a reaction when the two outermost beads come within the capture distance $r_\mathrm{c}$. This reaction is implemented in lines 50-52 of main.lmp. In this example, the capture distance is set to $r_\mathrm{c} = 1.0$. To change it, modify the 1.0 on line 51 to the desired value. Note that a larger capture distance may require small changes to both the main script and data file. These are explained in: https://github.com/WLaeremans/Polymer-looping-under-tension-FJC-

Additionally, this reaction setup requires that the two outermost beads on each end of the polymer are of type 2 and 3 on the upper side, and type 3 and 2 on the lower side, as specified in the example file polymer50.data.

To enable output of the exact timestep when a loop forms, we slightly modified the LAMMPS source code. Specifically, we changed fix_bond_react.cpp as shown in this commit: https://github.com/wouterel/lammps/commit/314e9ca4e797eabb391b2d5b575cb6f5e259d4e0. With this modification, main.lmp prints the exact timestep at which a loop is formed.

Finally, the file main_sampling.lmp is provided to output the distance between the end beads every 1000 timesteps. This data was used to numerically extract the free energy profiles, as presented in the appendix of the publication.

MATLAB

For the MATLAB code, MATLAB R2023a was used. The first script comparison_experiment.m reproduces Fig. 2 in the article.

The second script 3_theories.m reproduces Fig. 3 in the article and requires the data for the timesteps of loop formation in the zip file data_WLC_3_theories.zip. In order for the code to run, the data should be extracted into a folder located in the same directory as the code and added to the main MATLAB path. The data files are named as follows: Nx_fy.txt, where x represents the value of $N$ (with $N + 1$ beads), and y is the value of the force. Each file contains 100 timesteps, corresponding to 100 simulations with different random seeds for each combination of $N$ and $f$.

To reproduce Fig. 4 of the article, the script finite_size.m can be used. This requires the zip files data_finite_size.zip and data_WLC_3_theories.zip. The former contains loop formation timesteps in files named N22_fy_fse.txt, where y again denotes the force. For $N = 11$, the bond length was $b = 1$, while for $N = 22$, the bond length was $b = 0.5$. Additionally, the bending stiffness $\kappa$ was doubled (LAMMPS main.lmp line 18, 5 was replaced by 10) for $N = 22$ so that both cases represent the same physical polymer. As before, to run finite_size.m, the corresponding data folders must be added to the MATLAB path.

The script zero_force_scaling.m reproduces Fig. 5 in the article and requires the data in data_zero_force_scaling.zip.

To reproduce Fig. 6, run the script diff_rc.m. The required data is in the folder diff_rc_data.zip, which should be extracted and added to the MATLAB path. This folder includes the same data for $N = 11$ and $N = 23$ as in data_WLC_3_theories.zip, along with additional data for $r_\mathrm{c} = 0.5$. For example, the file N23_f0_6_rc05.txt corresponds to $N = 23$, $f = 0.6$, and $r_\mathrm{c} = 0.5$.

The final figure in the manuscript (Fig. 7) can be generated using the script corrected_free_energy.m. Before running, load the following .mat files into the MATLAB workspace: final_F_r_N11_f0_with_fit.mat, final_F_r_N23_f0_3_with_fit.mat, final_F_r_N23_f0_6_with_fit.mat and final_F_r_N23_f0_with_fit.mat. The last three MAT-file objects contain the free energy data extracted from simulation and fitted curves, necessary to apply Eq. 23 in the manuscript and generate subfigure $a$. The first MAT-file object contains the fit for the zero-force case when $N = 11$, which is used in the numerical integration for the barrier escape calculation shown in subfigure $b$. These free energies were generated by running the LAMMPS script main_sampling.lmp, and subsequently binning the output distances to obtain a probability distribution, as explained in the article.