Base pair number not staying constant in simulations

[piia600]

Hi all,

I have managed to briefly analyze my simulations with dnaMD Python module. I managed to analyze the curvature of the unbound DNA but I ran into issues while analyzing the bound DNA.

The 3DNA calculations finish fine, but they identify different number of base pairs in different frames. There are 40 base pairs in the starting structure, but in some frames this drops to 35.

I can load the helical axis with step range defined:

complexdna.set_helical_axis('complex2/HelAxis_complex2_test.dat',step_range=True, step=[10,25])

But I can't continue the curvature analysis, because the command

complexdna.generate_smooth_axis(step_range=True, step=[10,25], smooth=500, spline=2, fill_point=6, cut_off_angle=180)

Leads to an error (which happens around the time when the base pair number goes very low):

Fitting spline curve on helical axis of frame 6000 out of 80000 frames
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[24], line 2
      1 # Generate global axis by interpolation (smoothening)
----> 2 complexdna.generate_smooth_axis(step_range=True, step=[10,25], smooth=500, spline=2, fill_point=6, cut_off_angle=180)
      4 # Calculate curvature and tangent along global helical axis
      5 #freedna.calculate_curvature_tangent(store_tangent=True)

File ~/.local/lib/python3.9/site-packages/dnaMD/dnaMD.py:1322, in DNA.generate_smooth_axis(self, step_range, step, smooth, spline, fill_point, cut_off_angle)
   1318     sys.stdout.write("\rFitting spline curve on helical axis of frame %d out of %d frames" % (
   1319         frame_number, nframes))
   1320     sys.stdout.flush()
-> 1322 xsmooth, ysmooth, zsmooth, mask = fit_axis(bp_idx, frame_number, RawX, RawY, RawZ, smooth, spline, fill_point, cut_off_angle)
   1323 maskArray = mask
   1325 smoothX.append(xsmooth)

File ~/.local/lib/python3.9/site-packages/dnaMD/dnaMD.py:2621, in fit_axis(bp_idx, nframe, RawX, RawY, RawZ, smooth, spline, fill_point, cut_off_angle)
   2618 points = fill_point * len(orig_x)
   2620 nest = -1
-> 2621 tckp, u = splprep([orig_x, orig_y, orig_z], s=smooth, k=spline, nest=nest)
   2623 xnew, ynew, znew = splev(np.linspace(0, 1, points), tckp)
   2625 new_axis = np.array([xnew, ynew, znew]).T

File /usr/local/lib/python3.9/site-packages/scipy/interpolate/_fitpack_py.py:155, in splprep(x, w, u, ub, ue, k, task, s, t, full_output, nest, per, quiet)
     13 def splprep(x, w=None, u=None, ub=None, ue=None, k=3, task=0, s=None, t=None,
     14             full_output=0, nest=None, per=0, quiet=1):
     15     """
     16     Find the B-spline representation of an N-D curve.
     17
   (...)
    152
    153     """
--> 155     res = _impl.splprep(x, w, u, ub, ue, k, task, s, t, full_output, nest, per,
    156                         quiet)
    157     return res

File /usr/local/lib/python3.9/site-packages/scipy/interpolate/_fitpack_impl.py:175, in splprep(x, w, u, ub, ue, k, task, s, t, full_output, nest, per, quiet)
    173 wrk = _parcur_cache['wrk']
    174 iwrk = _parcur_cache['iwrk']
--> 175 t, c, o = _fitpack._parcur(ravel(transpose(x)), w, u, ub, ue, k,
    176                            task, ipar, s, t, nest, wrk, iwrk, per)
    177 _parcur_cache['u'] = o['u']
    178 _parcur_cache['ub'] = o['ub']

ValueError: Invalid inputs.

Also the 3DNA HelAxis*.dat files contained ---- for some base pairs which were not understood by dnaMD. I needed to replace them with 0.00 to get the file read.

Could you help me get the curvature analyzed also for the bound DNA? I only used the DNA atoms in the analyzed bound trajectory to get a comparison with the unbound DNA.

If needed, I can share some files for testing.

[xiangjun]

Thanks for the detailed report. While @rkumar should be the primary contact for dnaMD related issues, I'd like to address the general question of why the base pair number not staying constant in simulations.

In 3DNA, the find_pair program is used to identify base pairs in an input structure. For MD trajectories, when each frame is processed with auto-detected base pairs, the numbers can fluctuate due to the dynamic nature of the system. The 3DNA suite includes the Ruby script x3dna_ensemble, and the beginning portion of the "x3dna_ensemble analyze -h" command is as follows. Basically, it requires a template base-pair input file, possibly generated with ‘find_pair’ and manually edited as necessary.

------------------------------------------------------------------------
Analyze a MODEL/ENDMDL delineated ensemble of NMR structures or MD
trajectories. All models must correspond to different conformations
of the same molecule. For the analysis of duplexes (default), a template
base-pair input file, generated with 'find_pair' and manually edited
as necessary, must be provided.

Usage:
        x3dna_ensemble analyze options
Examples:
        x3dna_ensemble analyze -b bpfile.dat -e sample_md0.pdb

In DSSR, the --nmr (or --md) option can be used with --pair-list-input to analyze MD trajectories with a customized set of base pairs of interest. See the DSSR User Manual for more details in Sections "3.13 The --nmr option" and "3.9 The --pair-list options".

Best regards,

Xiang-Jun

[rkumar]

Hi,

Apologies for the late reply. Please use  -ref option to keep the base-pairs constant. However, any parameters calculated for broken base-pairs will be of extreme values. If base-pairs are broken at the terminals, then those could be ignored during analysis by dnaMD.

Thanks,
Rajendra