Location of interaction threshold energies? (General question)

[dhollis6]

Hi there,

I have a program which takes the output from output_bonds.py (in analysis tools) and determines whether any stacking or hydrogen bonds are broken throughout the trajectory. I need to define threshold values so that if the energy of an interaction is below this threshold, it is considered bonded.

It looks like the threshold for hydrogen bonding (HB_CUTOFF) is -0.1 as defined in src>Utilities>OrderParameters.h. However I can't find the threshold for normal, coaxial or cross stacking. The stacking observable UnstackedList.h references a value called threshold_fraction, and reads "the nucleotides will be counted as stacked if their stacking energy is at least this fraction of the maximum stacking energy." Where is this defined? And where do I find the maximum stacking energy?

Any guidance would be greatly appreciated. Thank you.

Kind regards Daniel H.

[lorenzo-rovigatti]

Hi,

the maximum stacking energy depends on temperature T, on the model and on whether you are using the sequence-dependent or average-sequence parametrisation. For the average-sequence models:

• average oxDNA: -(1.3448 + 2.6568 * T)
• average oxDNA2: -(1.3523 + 2.6717 * T)
• average oxRNA: -(1.40206 + 2.77 * T)

For the sequence-dependent models the maximum stacking energy also depends on the pair of nucleotides you are looking at:

• seq-dep oxDNA/oxDNA2: -STCK_X_Y (1 - STCK_FACT_EPS + (T 9.0 * STCK_FACT_EPS)
• seq-dep oxRNA: -STCK_X_Y (1.0 + T ST_T_DEP)

where STACK_FACT_EPS and ST_T_DEP and STCK_X_Y (where X and Y are A, C, G, T or U if RNA) can be found in the files that specify the sequence-dependent parametrisations (usually oxDNA/oxDNA1_sequence_dependent_parameters.txt, oxDNA/oxDNA2_sequence_dependent_parameters.txt and oxDNA/rna_sequence_dependent_parameters.txt).

Edit: usually nucleotides are either bonded/stacked or not, meaning that the specific value of the threshold is not super important, as long as you don't pick a value that's too close to the maximum. Using something like -0.2-0.3 should work, but look at the distributions of the energies you get from your simulations before choosing!

[dhollis6]

This is very helpful, thank you! And your last comment confirms a suspicion I had, since I noticed these thresholds are never set to precise values.

I'm assuming these values are only for normal stacking. Would the maximum stacking energies for cross stacking or coaxial stacking be much different? i.e. would the same thresholds work? I'm using the average sequence model.

[lorenzo-rovigatti]

The functional form of the cross- and coaxial stacking is a bit different than the stacking one, and it's a bit harder to obtain the numbers you are after just looking at the code. It's probably faster if you simulate a fully-formed duplex, use output_bonds and look at the energy distributions for the cross- and coaxial-stacking terms.

[dhollis6]

That makes sense. Thanks again!