## 2 Computational Details  
### Molecular Dynamic Simulations  
In order to study the drug intercalation into DNA, a B-DNA dodecamer, with Protein Data Bank code 1 BNA [23], was used. SAN was manually introduced into the different potential intercalation sites described in the next section. Hydrogens were added to the model system using H++ web-server [24; 25]. In order to obtain SAN parameters, a gas phase geometry optimization was carried out at density functional theory (DFT) level employing the hybrid B3 LYP [26; 27] as exchange and correlation functional coupled with 6-31G(d,p) basis set and using Gaussian 16 package [28]. Atomic charges are derived by fitting the electrostatic potential according to the Merz-Singh-Kollman scheme [29], usingthe RESP fitting procedure. The Lennard-Jones parameters of the GAFF force field were used for SAN [30]. Developed parameters are reported in the Supplementary Materials. Molecular dynamics topology and coordinate files were generated by means of tleap in Amber [31] using the standard DNA.bsc [32] and GAFF force fields [30; 33] together with the newly generated parameters for SAN. A TIP3P solvation model [34; 35; 36] was used to construct an octahedral box of explicit water molecules around the DNA dodecamer with a 14 A buffer distance around the DNA in each direction. Sodium ions were added to neutralize the system. The system was relaxed prior to the molecular dynamics by a 1000 minimization steps with a cutoff distance of 15.0 A and a constant volume of periodic boundaries. DNA and SAN were fixed during this initial minimization by using a force constant of 500 kcal mol\({}^{-1}\). This step was followed by 2500 minimization steps for the whole system including the DNA with the same cutoff distance and constant volume periodic boundaries. The system was then heated from 0 K to 300 K over 10,000 steps for a total of 20 ps. SHAKE algorithm was implemented to constrain bonds involving hydrogen. During the heating, DNA and SAN were restrained by a force constant of 10 kcal mol\({}^{-1}\) while keeping the constant volume periodic boundaries and the same cutoff distance. Equilibration followed by production of molecular dynamics (MD) for 100 ns at 300 K were ran under similar conditions with a 0.002 ps interval with no restraints on DNA, a constrain to bonds involving hydrogen and the cutoff distance of 15.0 angstroms were all maintained. NPT molecular dynamics at 1 atm pressure using the Berendsen barostat and an isotropic position scaling with a relaxation time of 2 ps were used. For evaluating the relative binding free energy between SAN, at different intercalation sites, and DNA, the MM-PBSA method [37], as implemented by MMPBSA.py script in the Amber 16 package, was used considering the whole production trajectory of each intercalation position over 100 ns. Calculation of the interaction energy and solvation free energies for the drug-DNA adducts and DNA and drug alone was followed by results averaging to obtain an estimate of the binding free energy as implemented in the MM-PBSA method in Amber 16. To obtain data about the root mean squared deviation (RMSD), distances between residues, and conformational changes in the DNA in the drug-DNA system, cptrraj [38] of Amber was used. cptrraj was also used to extract the most representative conformations of the SAN-DNA systems, and cluster analysis was performed adopting the algorithm of hierarchical agglomerative (bottom-up approach) kind. RMSD of eight surrounding base pairs was used as parameter for clustering. The clustering ends when the suggested optimal number of clusters, that is 10, remains.