Reconstructing the Free Energy Profiles Describing the Switching Mechanism of a pH-Dependent DNA Nanodevice from ABMD Simulations

The pH-responsive behavior of six triple-helix DNA nanoswitches, differing in the number of protonation centers (two or four) and in the length of the linker (5, 15 or 25 bases), connecting the double-helical region to the single-strand triplex-forming region, was characterized at the atomistic level through Adaptively Biased Molecular Dynamics simulations. The reconstruction of the free energy profiles of triplex-forming oligonucleotide unbinding from the double helix identified a different minimum energy path for the three diprotic nanoswitches, depending on the length of the connecting linker and leading to a different per-base unbinding profile. The same analyses carried out on the tetraprotic switches indicated that, in the presence of four protonation centers, the unbinding process occurs independently of the linker length. The simulation data provide an atomistic explanation for previously published experimental results showing, only in the diprotic switch, a two unit increase in the pKa switching mechanism decreasing the linker length from 25 to 5 bases, endorsing the validity of computational methods for the design and refinement of functional DNA nanodevices.


Introduction
The possibility to synthesize DNA sequences of any length and scale, the specificity of Watson-Crick base pairing and the possibility to introduce chemical modifications has led to the current use of DNA as an efficient nanoscale building material [1][2][3]. Indeed, DNA was exploited to build static three-dimensional structures, based on polyhedral geometries [4][5][6][7], as well as objects with complex shapes, mainly based on the DNA origami assembly technique [8][9][10][11]. DNA nanotechnology is now applied to solve real-world problems, developing functional and dynamic structures such as nanodevices and nanomachines [9,[12][13][14][15][16]. One important application concerns biological sensing, where nanosensors have the potential to make simple, cheap, and fast the detection of specific biological materials [17][18][19][20]. DNA triplexes have drawn a lot of attention for their pH-dependent conformational changes, which can be integrated into complex nanomachines [21,22]. Using DNA triplexes as pH-sensitive locks, a selective, reversible and enantioselective control of reconfigurable chiral plasmonic metamolecules assembled by DNA origami was implemented [23]. Similarly, the Linko group designed a DNA origami nanocapsule equipped with DNA triple helix locks, which can be loaded with various types of molecular cargo and is able to switch between open and closed conformational states upon a pH change, and thus expose or protect the encapsulated molecules [24]. Recently, a O6-methyl-guanine 2 of 12 (O6-MeG) modified triple helix-based nanoswitch was designed, with the ability to monitor the activity of DNA repair enzymes. [25] In this context, the full understanding of the DNA triplex-based devices' atomistic behavior represents an important aspect for a fine prediction of the rules governing their structure/dynamics/function relationship. In the last few years, we demonstrated the importance of Molecular Dynamics (MD) simulations coupled to experiments in describing and predicting the atomistic behavior of DNA nanoswitches [26,27] integrated into complex nanostructures [28]. Recently, it was experimentally shown that the pH-responsive behavior of a nucleic acid nanoswitch, which can form an intramolecular triplex structure through hydrogen bonds (Hoogsteen interactions) between a hairpin double helix (DH) and a single-strand triplex-forming oligo (TFO) with two protonation centers, strongly depends on the length of the linker connecting the two domains ( Figure 1A) [29]. Furthermore, it has also been demonstrated that the linker-dependent modulation can be dissipated by the introduction of four protonation centers in the single stranded triplex-forming portion [29]. In this work, we apply Adaptively Biased Molecular Dynamics (ABMD) simulations to describe the effect of varying the length of the linker from 5 to 25 bases in regulating the pHdependent conformational unbinding of the TFO from the DH, simulating three diprotic and three tetraprotic DNA nanoswitches. The results provide an atomistic depiction of the mechanism of unbinding, and they indicate that enhanced sampling techniques can be a valuable tool to rationally design the function of these devices, allowing a correct prediction of their real-world applications.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 2 of 12 states upon a pH change, and thus expose or protect the encapsulated molecules [24]. Recently, a O6-methyl-guanine (O6-MeG) modified triple helix-based nanoswitch was designed, with the ability to monitor the activity of DNA repair enzymes. [25] In this context, the full understanding of the DNA triplex-based devices' atomistic behavior represents an important aspect for a fine prediction of the rules governing their structure/dynamics/function relationship. In the last few years, we demonstrated the importance of Molecular Dynamics (MD) simulations coupled to experiments in describing and predicting the atomistic behavior of DNA nanoswitches [26,27] integrated into complex nanostructures [28]. Recently, it was experimentally shown that the pH-responsive behavior of a nucleic acid nanoswitch, which can form an intramolecular triplex structure through hydrogen bonds (Hoogsteen interactions) between a hairpin double helix (DH) and a single-strand triplex-forming oligo (TFO) with two protonation centers, strongly depends on the length of the linker connecting the two domains ( Figure 1A) [29]. Furthermore, it has also been demonstrated that the linker-dependent modulation can be dissipated by the introduction of four protonation centers in the single stranded triplex-forming portion [29]. In this work, we apply Adaptively Biased Molecular Dynamics (ABMD) simulations to describe the effect of varying the length of the linker from 5 to 25 bases in regulating the pH-dependent conformational unbinding of the TFO from the DH, simulating three diprotic and three tetraprotic DNA nanoswitches. The results provide an atomistic depiction of the mechanism of unbinding, and they indicate that enhanced sampling techniques can be a valuable tool to rationally design the function of these devices, allowing a correct prediction of their real-world applications.

DNA Nanoswitch Modelling
The scheme of the six simulated diprotic or tetraprotic triplex-based nanoswitches is shown in Figure 1A and in Figure S1A of Supplementary Materials. The diprotic systems share the same three-dimensional structure, except for the loop connecting the DH to the TFO, which is composed of 5, 15 or 25 bases for the DIPRO5 ( Figure 1B), DIPRO15 ( Figure  1C) and DIPRO25 ( Figure 1D) nanoswitches, respectively. The tetraprotic systems share the same structural features but contain four protonation centers in the TETRA5,

DNA Nanoswitch Modelling
The scheme of the six simulated diprotic or tetraprotic triplex-based nanoswitches is shown in Figure 1A and in Figure S1A of Supplementary Materials. The diprotic systems share the same three-dimensional structure, except for the loop connecting the DH to the TFO, which is composed of 5, 15 or 25 bases for the DIPRO5 ( Figure 1B), DIPRO15 ( Figure  1C) and DIPRO25 ( Figure 1D) nanoswitches, respectively. The tetraprotic systems share the same structural features but contain four protonation centers in the TETRA5, TETRA15 and TETRA25 nanoswitches ( Figure S1B-D of Supplementary Materials, respectively). The sequences used for the model building are identical to those reported in the experimental work [29], and are displayed here:  Figure S1B-D of Supplementary Materials). The three-dimensional structures of the six nanoswitches were built using the fiber and mutate_bases modules of the X3DNA package [31], which were used to independently produce a 10-base triple helix and the single strands representing the TFO and the loops. The PyMol sculpting module (PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC; https://pymol.org, accessed on 3 March 2021) was used to connect the DNA fragments through phosphodiester bonds.

MD and ABMD Simulations
The topology and coordinate files of the six nanoswitches were generated using the tLeap module of the AmberTools program [32], parametrizing the structures through the AMBER parmbsc1 [33] force field, mimicking high pH values (i.e., deprotonated cytosines). The structures were placed into a rectangular box, solvated with TIP3P water molecules [34] and neutralized by adding Mg 2+ ions, forcing a minimum distance between the structure and the box sides of 16 Å. For each structure, a minimization run of 500 steps using the steepest descent algorithm followed by 1500 steps of a conjugate gradient was performed in order to remove any unfavorable interaction. The systems were gradually heated from 0 to 300 K in the NVT ensemble over a period of 500 ps using the Langevin thermostat [35], applying a restraint of 0.5 kcal·mol −1 ·Å −2 on each nucleotide atom to relax the solvent. Through 1.0 nslong equilibration runs, the restraint forces were gradually decreased to 0.1 kcal·mol −1 ·Å −2 . The systems were simulated using an isobaric-isothermal (NPT) ensemble for 1.0 ns, setting the temperature to 300 K and the pressure to 1.0 atm using the Langevin barostat [36]. The SHAKE algorithm [37,38] was used to constrain the covalent bonds involving hydrogen atoms. The systems were then subjected to a 10.0 ns equilibration run before starting the ABMD simulations, with a time step of 2.0 fs, using the PME method [39] for long-range interactions and a cut-off of 9.0 Å for the short-range interactions.
The unbinding of the triplex-forming region of the equilibrated systems was simulated by means of 50.0 ns-long WT-ABMD simulations [40,41] using the NFE-toolkit of the Amber suite. ABMD enhances the sampling of the high-energy regions of the free energy surface (FES) by adding a biasing history-dependent term to the total energy of the system, obtained by the sum of the Gaussian hills laying on the subspace described by a set of user-defined collective variables (CVs). The most visited regions of conformational space (low free energy) are slowly filled and the biasing potential boosts the exploration of the less probable (high free energy) regions. The FES can be reconstructed at the end of the simulation as the sum of the added Gaussian hills using the nfe-umbrella-slice tool of the Amber suite. In the simulation, two CVs were used to describe the triplex-forming oligo unbinding from the duplex region: (1) the distance between the center of mass of the C2 atoms of the nucleotides belonging to the duplex and that of the C2 atoms belonging to the triplex-forming region; (2) the number of hydrogen bonds between the triplex-forming region and the nucleotides belonging to the duplex, using the coordination number (CN) CV as implemented in AMBER: where r i and r j are the coordinates of the atoms involved in the hydrogen bond interactions. The value of the parameter r 0 was set to 3.2 Å, as calculated from the 10.0 ns MD equilibration. The data here obtained cannot be related to the TFO dissociation constant because a full sampling of the unbound state is required for to calculate K d . However, the evaluated free energy difference between the transition state (Ts) and the bound state (B) represents an estimate of the 'activation energy' for the unbinding process and, therefore, is related to the different pH-dependent behavior of the switches.

Trajectory Analyses
RMSF, principal component analyses (PCA), distance, and hydrogen bond percentages were calculated over the entire 50-ns ABMD trajectories using the GROMACS 2020.3 tools [42]. The hydrogen bond number was evaluated through the gmx hbond module, using an angle cut-off (hydrogen-donor-acceptor) of 180 ± 30 • and a maximum donoracceptor distance of 3.5 Å. Cluster analysis was performed through the gmx cluster module of GROMACS, applying the gromos algorithm [43] on all of the saved configurations. The DNA geometrical parameters were estimated using the CURVES+ program [44]. Contour plots in Figures 2 and 3 were generated with the Matplotlib Python library [45]; the plots in Figures 4-6 and Figure S2 were produced using R and the ggplot2 package [46], while the structures represented in Figures 1-3 and Figure S1 were drawn with UCSF Chimera [30]. The analyses were performed using one node, for a total of 48 CPUs, on the CRESCO6 partition of the ENEA HPC cluster [47].

Unbinding of the TFO from the Double Helix
The TFO unbinding profile from the double helix was analyzed in terms of its free energy surface (FES) and plotted as a function of the number of hydrogen bonds (HBs) established between the TFO and the DH, and of the distance between their centers of mass, as reported in Figures 2 and 3 for the diprotic and tetraprotic switches, respectively. For all of the six systems, the FES shows the presence of three minima corresponding to the bound state (B) and two intermediate (I 0 and I 1 ) states, separated by two transition states (Ts and Ts unb ). In particular, the Ts unb state represents the final step of the ABMD simulations ( For the DIPRO5 switch, the B state is stable, being characterized by the presence of 13-14 HBs and a distance between TFO and DH centers of mass close to 2 Å (Figure 2A). The I 0 state has a lower number of HBs and an average TFO-DH distance of 6 Å (Figure 2A) due to the unbinding of three nucleotides at the 3 end of the TFO. In the I 1 state, the number of HB interactions is reduced to four, and in the final Ts unb state only two to three HBs are present. In fact, in the Ts unb state, the TFO is not fully detached because the first three bases at its 5 end are still interacting with the DH. The per-base unbinding profile of the DIPRO5 switch, reported in Figure 4A, provides a clear view of the detaching mechanism described by the FES, showing the progressive loss of HBs from the 3 to the 5 end of the TFO. The DIPRO15 system shows the presence of less stable minima than those observed in the DIPRO5 one ( Figure 2B). In detail, the B state shows the presence of 8-10 HBs, although the TFO is still quite close to the DH. The I 0 state is similar to that of the DIPRO5 switch, showing the unbinding of the first three nucleotides at the TFO Appl. Sci. 2021, 11, 4052 5 of 12 3 end. The I 1 state samples several loose structures, showing a concurrent loss of three base pair interactions at both the 3 and 5 ends of the TFO ( Figure 2B). Starting from this configuration, the Ts unb state is reached by an extensive sampling of a high-energy conformational basin, characterized by a TFO-DH distance of about 20 Å and the loss of all HBs (Figure 2B), indicative of a full unbinding of the TFO from the DH. The per-base unbinding profile of the DIPRO15 confirms the removal of the base-pair interactions at both the 3 and 5 ends of the TFO, followed by the loss of the internal bases' interactions, leading to the complete detachment of the TFO from the DH ( Figure 4B). The DIPRO25 switch is the less stable one, being characterized by high-energy minima, with the B state easily evolving towards the I 0 state ( Figure 2C). In this case, the unbinding process of the TFO from the DH, as observed through the structures sampled in the I 0 and I 1 states and from the per-base unbinding profile in Figure 4C, follows an opposite path when compared to the DIPRO5 system because the detachment starts from the TFO 5 end. The Ts unb state, identified by structures with a TFO-DH distance of 24-25 Å and characterized by a complete loss of HBs interactions, is reached through the sampling of a high-energy conformational basin ( Figure 2C).
HBs, although the TFO is still quite close to the DH. The I0 state is similar to that of the DIPRO5 switch, showing the unbinding of the first three nucleotides at the TFO 3′ end. The I1 state samples several loose structures, showing a concurrent loss of three base pair interactions at both the 3′ and 5′ ends of the TFO ( Figure 2B). Starting from this configuration, the Tsunb state is reached by an extensive sampling of a high-energy conformational basin, characterized by a TFO-DH distance of about 20 Å and the loss of all HBs ( Figure  2B), indicative of a full unbinding of the TFO from the DH. The per-base unbinding profile of the DIPRO15 confirms the removal of the base-pair interactions at both the 3′and 5′ ends of the TFO, followed by the loss of the internal bases' interactions, leading to the complete detachment of the TFO from the DH ( Figure 4B). The DIPRO25 switch is the less stable one, being characterized by high-energy minima, with the B state easily evolving towards the I0 state ( Figure 2C). In this case, the unbinding process of the TFO from the DH, as observed through the structures sampled in the I0 and I1 states and from the per-base unbinding profile in Figure 4C, follows an opposite path when compared to the DIPRO5 system because the detachment starts from the TFO 5′ end. The Tsunb state, identified by structures with a TFO-DH distance of 24-25 Å and characterized by a complete loss of HBs interactions, is reached through the sampling of a high-energy conformational basin ( Figure 2C). For all the tetraprotic switches, the B state is characterized by metastable high energy minima, which evolve towards the I0 and I1 states ( Figure 3A-C) by sampling several conformations characterized by a rapid loss of HBs. The I1 states of these switches are characterized by low energy minima in which the bases at the 5′ end of the TFO are still interacting with the duplex. The Tsunb state, which is characterized by structures with a TFO-DH distance of 23-24 Å , and by a complete loss of HBs interactions, is reached through the sampling of a high-energy conformational basin, as observed for the diprotic switch (Figure 3A-C). For all the tetraprotic switches, the B state is characterized by metastable high energy minima, which evolve towards the I 0 and I 1 states ( Figure 3A-C) by sampling several conformations characterized by a rapid loss of HBs. The I 1 states of these switches are characterized by low energy minima in which the bases at the 5 end of the TFO are still interacting with the duplex. The Ts unb state, which is characterized by structures with a TFO-DH distance of 23-24 Å, and by a complete loss of HBs interactions, is reached through the sampling of a high-energy conformational basin, as observed for the diprotic switch ( Figure 3A-C).
These observations are confirmed by the per-base unbinding profiles that are similar for all of the three tetraprotic switches, which are characterized by the removal of the basepair interactions starting from the 3 to the 5 end of the TFO ( Figure 4D-F), independently of the linker length.  These observations are confirmed by the per-base unbinding profiles that are similar for all of the three tetraprotic switches, which are characterized by the removal of the basepair interactions starting from the 3′ to the 5′ end of the TFO ( Figure 4D-F), independently of the linker length.
A comparative free energy profile along a minimum energy path, as a function of a generalized reaction coordinate, is plotted in Figure 5A and B for the DIPRO and TETRA switches, respectively. The energy value of the B state for the DIPRO5 switch ( Figure 5A red line) is two times lower than those observed for the DIPRO15 and DIPRO25 systems ( Figure 5A green and blue lines, respectively). Moreover, the DIPRO5 B and I1 states are separated by a relatively high activation energy ( Figure 5A, red line), again almost twice than that required for the TFO unbinding in the DIPRO15 and DIPRO25 systems ( Figure  5A, green and blue lines, respectively). In the case of the TETRA5 ( Figure 5B, red line), TETRA15 ( Figure 5B, blue line) and TETRA25 switches ( Figure 5B, green line), the energy values observed for the B state are half of those observed for the diprotic switches, with the I0 and I1 states separated by very low activation energies, which favour the transition towards the unbound states. A comparative free energy profile along a minimum energy path, as a function of a generalized reaction coordinate, is plotted in Figure 5A,B for the DIPRO and TETRA switches, respectively. The energy value of the B state for the DIPRO5 switch ( Figure 5A red line) is two times lower than those observed for the DIPRO15 and DIPRO25 systems ( Figure 5A green and blue lines, respectively). Moreover, the DIPRO5 B and I 1 states are separated by a relatively high activation energy ( Figure 5A, red line), again almost twice than that required for the TFO unbinding in the DIPRO15 and DIPRO25 systems ( Figure 5A, green and blue lines, respectively). In the case of the TETRA5 ( Figure 5B, red line), TETRA15 ( Figure 5B, blue line) and TETRA25 switches ( Figure 5B, green line), the energy values observed for the B state are half of those observed for the diprotic switches, with the I 0 and I 1 states separated by very low activation energies, which favour the transition towards the unbound states.
These data indicate that the length of the linker connecting DH with TFO has a crucial role in stabilizing the DH-TFO interaction only in the presence of two protonation sites. In this case, the DH-TFO interaction energy strongly decreases upon increasing the loop length, indicating that the entropic contribution related to the loop length has a direct effect on the TFO unfolding energy. This effect is not observed for the tetraprotic switches. It is interesting to note that the experimental results indicate a two unit increase in the pKa switching mechanism only for the diprotic system, when the linker length is reduced from 25 to 5 bases [29], in line with the simulative findings presented here.   These data indicate that the length of the linker connecting DH with TFO has a crucial role in stabilizing the DH-TFO interaction only in the presence of two protonation sites. In this case, the DH-TFO interaction energy strongly decreases upon increasing the loop length, indicating that the entropic contribution related to the loop length has a direct effect on the TFO unfolding energy. This effect is not observed for the tetraprotic switches. It is interesting to note that the experimental results indicate a two unit increase in the pKa switching mechanism only for the diprotic system, when the linker length is reduced from 25 to 5 bases [29], in line with the simulative findings presented here.

Conformational Variability of the Linker and Stability of the Double Helix
The conformational variability of the linker region connecting the DH to TFO was evaluated through a PCA of the motion, carried out on the C2′ atom trajectories of the six systems. The projection of the loop C2′ atom trajectories on the plane defined by the first and second eigenvectors ( Figure 6A,B for the diprotic and tetraprotic switches, respectively) indicates that the loop in the six systems samples a very different conformational space. In the case of the DIPRO5 and TETRA5 systems (red circles), the sampled landscape is highly restricted, while it increases for the DIPRO15 and TETRA15 systems (green circles) and attains the widest space for the DIPRO25 and TETRA25 switches (blue circles). The reduced sampling in the DIPRO5 system is in agreement with the need to provide more energy to this system in order to observe the TFO unbinding, and it is also in line  It is interesting to note that the different conformational variability of the loop does not affect the stability of the double helix portion of the switch in addition to its flexibility.

Conformational Variability of the Linker and Stability of the Double Helix
The conformational variability of the linker region connecting the DH to TFO was evaluated through a PCA of the motion, carried out on the C2 atom trajectories of the six systems. The projection of the loop C2 atom trajectories on the plane defined by the first and second eigenvectors ( Figure 6A,B for the diprotic and tetraprotic switches, respectively) indicates that the loop in the six systems samples a very different conformational space. In the case of the DIPRO5 and TETRA5 systems (red circles), the sampled landscape is highly restricted, while it increases for the DIPRO15 and TETRA15 systems (green circles) and attains the widest space for the DIPRO25 and TETRA25 switches (blue circles). The reduced sampling in the DIPRO5 system is in agreement with the need to provide more energy to this system in order to observe the TFO unbinding, and it is also in line with previous results obtained by simulating the unfolding of similar clamp-triplex nanoswitches with different linker lengths [27].
It is interesting to note that the different conformational variability of the loop does not affect the stability of the double helix portion of the switch in addition to its flexibility. In fact, the DH RMSF are identical for both the diprotic and tetraprotic switches ( Figure S2A,B of Supplementary Materials, respectively), with the exception of the 5 end of the DIPRO25 switch (pink line, Figure S2A of Supplementary Materials), demonstrating that the conformational variability of the loop does not have a destabilizing effect on the double helix. In line with these observations, the averaged geometrical parameters of the 10 base pairs forming the DH in the diprotic and tetraprotic switches, shown in Tables S1 and S2 of Supplementary Materials, are close to those of the standard B-DNA DH, confirming that the duplex geometry is maintained over all of the simulation time, independently of the linker length and protonation centers.

Conclusions
In this manuscript, we analyzed the effect of increasing the length of the variable loop integrated on three diprotic and three tetraprotic triple helix nanoswitches, using the same sequences experimentally characterized in a recent work [29]. The atomistic investigation, carried out through ABMD simulations, permitted us to finely describe the unbinding process of the TFO from the double helical portion of the nanoswitch. In detail, the analysis of the six trajectories identified a different TFO unbinding path for the diprotic switch, depending on the length of the variable loop, both in terms of the energy required to switch from the bound to unbound states (Figures 2 and 4) and in terms of the per-base unbinding profile ( Figure 5A). In fact, the DIPRO5 switch requires, for the TFO unbinding, more energy than the DIPRO15 and DIPRO25 ones ( Figure 5A). Moreover, for the DIPRO5, the unbinding is incomplete because the bases directly connected to the variable loop continue to interact with the DH even in the final state ( Figure 4A). Finally, the per-base unbinding profile is different following a 3 to 5 direction for the DIPRO5 and a 5 to 3 for the DIPRO25 one, while it concomitantly occurs at both the 3 and 5 ends of the TFO for the DIPRO15, demonstrating the importance to the loop length in modulating the mechanism. Analyzing the results obtained for the tetraprotic switches, we did not find differences in the per-base unbinding profiles ( Figure 4B), and we observed only slight energy differences in the free energy profiles ( Figure 5B), confirming that linker-dependent modulation is dissipated upon increasing the number of protonation centers in the TFO.
These results provide an atomistic explanation to the experimental results showing, for the TFO unbinding in the diprotic switches, a decrease of more than two units of the pK a increasing the linker length from 5 to 25 bases, whilst the pH-dependent switching mechanism-dependent on the linker length-is lost when the protonation centers are four. Interestingly, our data indicate that the presence of a longer variable loop does not destabilize the DH, as can be observed in Figure S2A,B and in Tables S1 and S2 of Supplementary Materials. This is an important point to ensure a correct functioning of the switch as a pH-sensing device, because a loss of the DH stability would not permit the reversibility of the pH-dependent switching process.
Computational approaches, such as enhanced sampling MD simulations, proved to be of fundamental importance for the design and tuning of the functional properties of DNA nanoswitches. Software like mfold [48] or NUPACK [49], although allowing the prediction of the thermodynamics of Watson-Crick interactions, are unable to predict the behaviour of DNA nanodevices characterized by the presence of non-canonical interactions or sequence-independent effects. Therefore, the use of complex and advanced simulation techniques has become an essential tool in the field of DNA nanotechnology, which is also thanks to the reduced computational time required to perform MD simulations due to the improvement of the available workstations.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/app11094052/s1. Figure S1: Schematic and cartoon representations of the simulated tetraprotic DNA nanoswitches; Figure S2: Per-nucleotide RMSF values calculated for the DH of the six DNA nanoswitches; Table S1: DNA DH geometrical parameters calculated for the diprotic switches; Table  S2: DNA DH geometrical parameters calculated for the tetraprotic switches.