Accurate prediction of DNA-Intercalator binding energies: Ensemble of short or long molecular dynamics simulations?

Despite the wide use of molecular dynamics (MD) simulations for binding energy predictions in biomolecular systems, results from single MD simulations are non-reproducible and often deviate from experimental values, even when longer simulations are used. This study addresses these limitations using ensemble MD simulations for the formation of DNA-intercalator complexes. Twenty-five replicas of short (10 ns) and long (100 ns) MD simulations were performed on different intercalators binding into DNA. The MM/PBSA and MM/GBSA binding energies of the Doxorubicin intercalating into DNA, including entropy and deformation energy corrections, are -7.3 ± 2.0 kcal/mol and -8.9 ± 1.6 kcal/mol, using 25 replicas of 100 ns. These values were closely reproduced even with shorter simulations of 10 ns, where the energies, averaged over 25 replicas, are -7.6 ± 2.4 kcal/mol (MM/PBSA) and -8.3 ± 2.9 kcal/mol (MM/GBSA). In both cases, the energies align well with the experimental range of -7.7 ± 0.3 to -9.9 ± 0.1 kcal/mol. This shows that reproducibility and accuracy of the binding energies depend more on the number of replicas than on the simulation length. The study was repeated for the DNA-Proflavine system, where the corrected MM/PBSA and MM/GBSA binding energies, averaged over 25 replicas of 10 ns each, are -5.6 ± 1.4 and -5.3 ± 2.3 kcal/mol, respectively. These are congruent with the experimental range of -5.9 to -7.1 kcal/mol. Bootstrap analyses revealed that 6 replicas of 100 ns or 8 replicas of 10 ns provide a good balance between computational efficiency and accuracy within 1.0 kcal/mol from experimental values.