Dynamic Structure and Stability of DNA Duplexes Bearing a Dinuclear Hg(II)-Mediated Base Pair

Quantum mechanical (QM) and hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulations of a recently reported dinuclear mercury(II)-mediated base pair were performed aiming to analyse its intramolecular bonding pattern, its stability, and to obtain clues on the mechanism of the incorporation of mercury(II) into the DNA. The dynamic distance constraint was employed to find initial structures, control the dissociation process in an unbiased fashion and to determine the free energy required. A strong influence of the exocyclic carbonyl or amino groups of neighbouring base pairs on both the bonding pattern and the mechanism of incorporation was observed. During the dissociation simulation, an amino group of an adenine moiety of the adjacent base pair acts as a turnstile to rotate the mercury(II) ion out of the DNA core region. The calculations provide an important insight into the mechanism of formation of this dinuclear metal-mediated base pair and indicate that the exact location of a transition metal ion in a metal-mediated base pair may be more ambiguous than derived from simple model building.


Classical Molecular Dynamics
Classical molecular dynamics (MD) simulations were performed using the FIST method of CP2K [73] with the PARMbsc1 force field [70] for the DNA and the Na + counterions, with the modifications for isoguanine described above. Water molecules were described by the TIP3P water model [74,75]. A time step of 0.5 fs was chosen. Equilibration runs were performed in the NPT ensemble using the Nose-Hoover barostat and thermostat with a chain length of 3 and a time constant of 1 ps. The target temperature was set to 300 K in all NVT and NPT runs and the target pressure to 1 bar in all NPT MDs.

Ab Initio Molecular Dynamics and Optimization
DFT based ab initio molecular dynamics simulations were carried out with the CP2K program using the GAPW scheme [76,77]. The PBE functional was chosen together with the TZVP-MOLOPT basis set for all atoms except mercury, for which a DZVP-MOLOPT basis was utilized [78], in conjunction with Goedecker-Teter-Hutter pseudopotentials [79][80][81] and a plane wave cut-off of 400 Ry. Grimme D3 dispersion corrections were employed with Becke-Johnson damping [63,64]. A total of 700 virtual molecular orbitals were calculated for Fermi-Dirac smearing due to the presence of the transition metal mercury. Cholesky inversion was chosen as the solver of the eigenvalue problem, while Broyden mixing was employed for optimization. MD was performed with a time step of 0.5 fs.
The convergence criteria for the geometry optimization were set to 0.003 a.u. for the maximum step size, 0.0015 a.u. for the root-mean-square (RMS) displacement of atomic positions, 0.0004 a.u. for gradients and 0.003 a.u. for the RMS change of the gradients.

QM/MM Simulations
QM/MM molecular dynamics simulations were carried out with the CP2K program, where for the classical region the settings from 2.3 were applied. For the QM region, settings from Section 2.4 were used. The QM region consisted of all atoms of the metal-modified base pair ( Figure 1a) and the base pairs directly above and below. Furthermore, three negatively charged phosphate groups together with the attached backbone structure were incorporated into the QM subsystem to neutralize the threefold positive charge of the metal-modified base pair leading to an overall chargeneutral QM subsystem with multiplicity 1, which proved crucial for obtaining reasonable geometries and stable MD simulations. Coupling was treated with the Gaussian expansion of electrostatic potentials (GEEP) method in CP2K, where the IMOMM method was utilized to treat the links between MM and QM [82].

Dynamic Distance Constraint
Five different model systems termed structure 1-5, which will be presented in Section 3, were investigated by free energy calculations. As a reaction coordinate, we chose the dynamic distance [55].

Dynamic Distance Constraint
Five different model systems termed structure 1-5, which will be presented in Section 3, were investigated by free energy calculations. As a reaction coordinate, we chose the dynamic distance [55].
where NOP are all non-overlapping pairs i, j of atoms between which the constraint is defined, r i and r j are their positions, µ ij = m i m j / m i + m j is their reduced mass and µ * = NOP µ ij is the sum of all reduced masses.
In the case of the full DNA models D = D DNA includes all N···H hydrogen bonds within regular base pairs and the two N···O distances bridged by Hg(II) ions ( Figure 1). In the case of the isolated base-pair models, the two O···N distances according to Figure 1 were selected to enter the dynamic distance D = D ISO .

Thermodynamic Integration
The dynamic distance was incrementally increased and the average Lagrange multiplier at every constraint value, corresponding to the average constraint force <F(D)>, calculated. Integrating over the mean constraint force along the dynamic distance gives the free energy necessary to dissociate the system: Structures 1, 2, and 3 were simulated for 500 fs at the initial dynamic distance before starting the dissociation MD, to allow for further equilibration of the quantum region. The running average of the Lagrange multipliers are shown in Figure S4 for structure 1, Figure S5 for structure 2 and Figure S6 for structure 3. For the calculations involving structure 1 , the dynamic distance was initially set to 3.0 Å and increased/decreased by 0.2 Å at each simulation point, where each point was simulated for at least 1 ps. Overall 41360 simulation steps were performed for structure 1 , equalling 20.58 ps. At the final point of D DNA = 5.8 Å of the trajectory, the modified base pair was found to be dissociated and one Hg(II) ion was in contact with water. Hence, the simulation was not continued further as this would have led to unphysical water-mercury interactions across the QM/MM boundary. Structure 2 was less stable during the dissociation simulation, thus requiring smaller dynamic distance steps of 0.1 Å. The system was simulated from an initial value of D DNA = 3.7 Å up to a final value of D DNA = 10.5 Å, at which point only a change in the bonding pattern occurred instead of dissociation. However, as one mercury ion was already in contact with MM water at this D DNA , the simulation was not continued further. Each D point was simulated for at least 1 ps, amounting to an overall simulation length of 62 ps. Structure 3 was simulated with a dynamic distance step size of 0.1 Å from D DNA = 2.5 Å to D DNA = 3.5 Å and a step size of 0.5 Å for the remaining points up to D DNA = 8.5 Å. The overall simulation length for structure 3 was about 35 ps, after which the modified base pair was sufficiently dissociated (see Figure S3).
In the all-QM simulations, structures 4 and 5 were dissociated by an increase of the dynamic distance between the N and O atoms of the modified base pair (Figure 8), termed D ISO , with settings given in Section 2.4. The system was simulated for 500 fs at the initial value of the dynamic distance constraint, which was subsequently increased in steps of 0.2 Å. An all-QM MD was performed with the CP2k code, with identical settings as before but with all water molecules described by DFT. A total of ca. 6 ps were simulated for structure 4 and ca. 13 ps for structure 5.

Model Structures
In total, five different model structures were simulated, three of which were complete DNA models and two were isolated base pairs in solution. Their generation will be described below.

Generation of the Initial DNA Geometry
Five different initial structures are simulated by means of molecular dynamics and geometry optimizations to study the bonding patterns of the novel artificial T:εA base-pair ( Figure 2). As no experimental structures of this mercury(II)-containing DNA duplex are available, the initial geometries were generated by using the 3DNA software [83] for the sequence 3 -d(GAAAGATAGGGAG)-5 /5 -d(CTCCCTATCTTTC)-3 . One DNA strand was then rotated by 180 • and placed back into a helical structure, to generate parallel-stranded DNA as used in the previous experiments ( Figure 3). This was followed by an exchange of the exocyclic O and NH 2 groups of all guanines, transforming them into isoguanines, thus generating the duplex 5 -( i GA i G i G i GATA i GAAA i G)-3 /5 -d(CTCCCTATCTTTC)-3 corresponding to the one used in the original experiments [52] except for the modification of the central A:T base pair.
Molecules 2020, 25, x 5 of 21 experimental structures of this mercury(II)-containing DNA duplex are available, the initial geometries were generated by using the 3DNA software [83] for the sequence 3′d(GAAAGATAGGGAG)-5′/5′-d(CTCCCTATCTTTC)-3′. One DNA strand was then rotated by 180° and placed back into a helical structure, to generate parallel-stranded DNA as used in the previous experiments ( Figure 3). This was followed by an exchange of the exocyclic O and NH2 groups of all guanines, transforming them into isoguanines, thus generating the duplex 5′-( i GA i G i G i GATA i GAAA i G)-3′/5′-d(CTCCCTATCTTTC)-3′ corresponding to the one used in the original experiments [52] except for the modification of the central A:T base pair.

Figure 2.
Originally proposed [52] bonding pattern of the dinuclear mercury(II)-mediated base pair between 1,N 6 -ethenoadenine (εA) and thymine (T), including the atom numbering according to IUPAC recommendations in blue. Only standard bases are available in the 3DNA program; thus, a thymine:adenine (T:A) base pair was generated at the site of the modified base pair. The atomic vacuum structure of the modified thymine:1,N 6 -ethenoadenine (T:εA) pair was optimized using the settings from Section 2.1. The electric charge of the base pair is +3, its spin multiplicity 1. The resulting structure is shown in Figure  4. This structure closely resembles that reported in reference 53 for this base pair.  [52] bonding pattern of the dinuclear mercury(II)-mediated base pair between 1,N 6 -ethenoadenine (εA) and thymine (T), including the atom numbering according to IUPAC recommendations in blue.
Molecules 2020, 25, x 5 of 21 experimental structures of this mercury(II)-containing DNA duplex are available, the initial geometries were generated by using the 3DNA software [83] for the sequence 3′d(GAAAGATAGGGAG)-5′/5′-d(CTCCCTATCTTTC)-3′. One DNA strand was then rotated by 180° and placed back into a helical structure, to generate parallel-stranded DNA as used in the previous experiments ( Figure 3). This was followed by an exchange of the exocyclic O and NH2 groups of all guanines, transforming them into isoguanines, thus generating the duplex 5′-( i GA i G i G i GATA i GAAA i G)-3′/5′-d(CTCCCTATCTTTC)-3′ corresponding to the one used in the original experiments [52] except for the modification of the central A:T base pair.  Only standard bases are available in the 3DNA program; thus, a thymine:adenine (T:A) base pair was generated at the site of the modified base pair. The atomic vacuum structure of the modified thymine:1,N 6 -ethenoadenine (T:εA) pair was optimized using the settings from Section 2.1. The electric charge of the base pair is +3, its spin multiplicity 1. The resulting structure is shown in Figure  4. This structure closely resembles that reported in reference 53 for this base pair. Only standard bases are available in the 3DNA program; thus, a thymine:adenine (T:A) base pair was generated at the site of the modified base pair. The atomic vacuum structure of the modified thymine:1,N 6 -ethenoadenine (T:εA) pair was optimized using the settings from Section 2.1. The electric charge of the base pair is +3, its spin multiplicity 1. The resulting structure is shown in Figure 4. This structure closely resembles that reported in reference 53 for this base pair.
To assess the influence of the neighbouring base pairs onto the equilibrium geometry of the artificial base pair, further DFT geometry optimizations including the base pairs above and below T: εA have been performed with identical functionals and basis set as the metal-mediated base pair alone. In this model, adjacent base pairs are connected via a phosphate backbone. Due to the four phosphate moieties connecting the three base pairs, the charge of the system is now −1 and the multiplicity 1. Figure 5 displays the resulting structure. This structure is more similar to the originally proposed base pair geometry from reference 52.  To assess the influence of the neighbouring base pairs onto the equilibrium geometry of the artificial base pair, further DFT geometry optimizations including the base pairs above and below T: εA have been performed with identical functionals and basis set as the metal-mediated base pair alone. In this model, adjacent base pairs are connected via a phosphate backbone. Due to the four phosphate moieties connecting the three base pairs, the charge of the system is now -1 and the multiplicity 1. Figure 5 displays the resulting structure. This structure is more similar to the originally proposed base pair geometry from reference 52. Two different bonding patterns emerged during the geometry optimization. Within the isolated base pair, the mercury ion Hg2 is located in between the εA-N6 and T-O2 atoms, while Hg1 coordinates to the εA-N7 and T-O4 atoms (NO bonding pattern, Figure 4). When the adjacent pairs are also taken into account, Hg2 remains located between the εA-N6 and T-O2 atoms, while Hg1 changes its bonding pattern and is now coordinated to εA-N7 and T-N3 (NN bonding pattern, Figure  5). These two structures correspond to two minima on the potential energy surface and exist for both the isolated and the stacked base pair. In the isolated case, the NO bonding pattern represents the global minimum structure. This has been confirmed for a number of different exchange-correlation functionals given in Section 2.1; an overview of the results is tabulated in Figure S7. In the three base pair stack, the NN bonding pattern is favoured ( Figure 5), where the neighbouring N and O atoms from the nucleobases above and below attract the mercury ions, if they are situated directly on top or below them.  To assess the influence of the neighbouring base pairs onto the equilibrium geometry of the artificial base pair, further DFT geometry optimizations including the base pairs above and below T: εA have been performed with identical functionals and basis set as the metal-mediated base pair alone. In this model, adjacent base pairs are connected via a phosphate backbone. Due to the four phosphate moieties connecting the three base pairs, the charge of the system is now -1 and the multiplicity 1. Figure 5 displays the resulting structure. This structure is more similar to the originally proposed base pair geometry from reference 52. Two different bonding patterns emerged during the geometry optimization. Within the isolated base pair, the mercury ion Hg2 is located in between the εA-N6 and T-O2 atoms, while Hg1 coordinates to the εA-N7 and T-O4 atoms (NO bonding pattern, Figure 4). When the adjacent pairs are also taken into account, Hg2 remains located between the εA-N6 and T-O2 atoms, while Hg1 changes its bonding pattern and is now coordinated to εA-N7 and T-N3 (NN bonding pattern, Figure  5). These two structures correspond to two minima on the potential energy surface and exist for both the isolated and the stacked base pair. In the isolated case, the NO bonding pattern represents the global minimum structure. This has been confirmed for a number of different exchange-correlation functionals given in Section 2.1; an overview of the results is tabulated in Figure S7. In the three base pair stack, the NN bonding pattern is favoured ( Figure 5), where the neighbouring N and O atoms from the nucleobases above and below attract the mercury ions, if they are situated directly on top or below them. Two different bonding patterns emerged during the geometry optimization. Within the isolated base pair, the mercury ion Hg2 is located in between the εA-N6 and T-O2 atoms, while Hg1 coordinates to the εA-N7 and T-O4 atoms (NO bonding pattern, Figure 4). When the adjacent pairs are also taken into account, Hg2 remains located between the εA-N6 and T-O2 atoms, while Hg1 changes its bonding pattern and is now coordinated to εA-N7 and T-N3 (NN bonding pattern, Figure 5). These two structures correspond to two minima on the potential energy surface and exist for both the isolated and the stacked base pair. In the isolated case, the NO bonding pattern represents the global minimum structure. This has been confirmed for a number of different exchange-correlation functionals given in Section 2.1; an overview of the results is tabulated in Figure S7. In the three base pair stack, the NN bonding pattern is favoured ( Figure 5), where the neighbouring N and O atoms from the nucleobases above and below attract the mercury ions, if they are situated directly on top or below them.
As the PBE functional correctly predicts the NO pattern to be the more stable configuration for the isolated base pair, it was chosen as the exchange-correlation functional for the QM/MM and the all-QM MD runs in view of its computational advantage compared with hybrid functionals. Relativistic corrections did not change the observed bonding pattern (see Figure S7) and were thus not included in the MD runs. Because the TZVP basis led to a further stabilization of the NO compared to the NO pattern, the TZVP-MOLOPT basis was chosen for use during the MD.
Replacing the nucleobases of the T:A base pair in the centre of the DNA duplex by the modified base pair T:εA within the geometry from Figure 4 led to the desired initial structure Relativistic corrections did not change the observed bonding pattern (see Figure S7) and were thus not included in the MD runs. Because the TZVP basis led to a further stabilization of the NO compared to the NO pattern, the TZVP-MOLOPT basis was chosen for use during the MD.
Replacing the nucleobases of the T:A base pair in the centre of the DNA duplex by the modified base pair T:εA within the geometry from Figure 4 led to the desired initial structure 5′-( i GA i G i G i GATA i GAAA i G)-3′ / 5-d'(CTCCCTεATCTTTC)-3′, which is depicted in Figure 6 along with the nucleobase numbering scheme. The topology of the DNA for classical MD was generated as described in Section 2.2.

Zipping up the DNA by a Dynamic Distance Constraint
The initially prepared structure was still in a very high-energy, non-equilibrated state. To relax the geometry, classical MD runs and geometry optimizations were performed using the dynamic distance constraint (Section 2.6.1). This enables fixing the overall hydrogen-bond structure of the DNA during equilibration, while allowing individual hydrogen bond distances to change, which is needed to equilibrate the double helix without it breaking apart into two strands. Because no parametrization of the N-Hg(II) and O-Hg(II) bonds within the modified bases was available, all atoms of the modified base pair (except for the deoxyribose moieties and the phosphate backbone) were fixed at the equilibrium geometry from the optimization from Section 3.1, Figure 4 during purely classical MD, with the computational details given in Section 2.3. With all hydrogen bonds between base pairs selected as a dynamic distance constraint DEQUI, this dynamic distance was gradually reduced from an initial value of DEQUI = 1.88 Å to DEQUI = 1.56 Å during an NPT MD simulation with the settings given in Section 2.3 over 10000 steps, thus "zipping up" the DNA hydrogen bond structure ( Figure 7). The topology of the DNA for classical MD was generated as described in Section 2.2.

Zipping up the DNA by a Dynamic Distance Constraint
The initially prepared structure was still in a very high-energy, non-equilibrated state. To relax the geometry, classical MD runs and geometry optimizations were performed using the dynamic distance constraint (Section 2.6.1). This enables fixing the overall hydrogen-bond structure of the DNA during equilibration, while allowing individual hydrogen bond distances to change, which is needed to equilibrate the double helix without it breaking apart into two strands. Because no parametrization of the N-Hg(II) and O-Hg(II) bonds within the modified bases was available, all atoms of the modified base pair (except for the deoxyribose moieties and the phosphate backbone) were fixed at the equilibrium geometry from the optimization from Section 3.1, Figure 4 during purely classical MD, with the computational details given in Section 2.3. With all hydrogen bonds between base pairs selected as a dynamic distance constraint D EQUI , this dynamic distance was gradually reduced from an initial value of D EQUI = 1.88 Å to D EQUI = 1.56 Å during an NPT MD simulation with the settings given in Section 2.3 over 10,000 steps, thus "zipping up" the DNA hydrogen bond structure ( Figure 7).

Relaxing the DNA Backbone by Geometry Optimization
Simulated annealing starting from the NPT-equilibrated temperature of 300 K was applied during the subsequent 30,000 steps with an annealing factor of 0.99 to slowly relax the DNA duplex. This was followed by a QM/MM geometry optimization [84,85] using CP2K with all constraints released. The computational details are given in Section 2.5.

Relaxing the DNA Backbone by Geometry Optimization
Simulated annealing starting from the NPT-equilibrated temperature of 300 K was applied during the subsequent 30,000 steps with an annealing factor of 0.99 to slowly relax the DNA duplex. This was followed by a QM/MM geometry optimization [84,85] using CP2K with all constraints released. The computational details are given in Section 2.5. Subsequent to the geometry optimization, a classical MD simulation was performed in the NPT ensemble for pressure equilibration for 1.7 ns and settings from Section 2.3, keeping all atoms of the modified base pair (except for the sugar and phosphate backbone) fixed due to the absence of suitable force field parameters, however releasing the dynamic distance constraint. As the hydrogen bond structure persisted during optimization and MD without the use of the dynamic distance constraint, we considered the DNA duplex to be in a chemical reasonable structure, which we chose as a basis for our further investigations. This resulted in a unit cell of size 43.29 Å × 44.53 Å × 63.33 Å, corresponding to a density of 985 kg/m 3 . The resulting structure is referred to as structure 1 (Figure 8).

Structure 2
The dynamic distance method was subsequently employed again as a reaction coordinate to simulate the dissociation of DNA and to determine the corresponding free energy profile. In this case, the dynamic distance DDNA comprises all N···H hydrogen bonds of all regular base pairs together with the two N···O distances bridged by the two Hg ions in the metal-mediated base pair (Figure 1a). Thereby all the A-N1···T-H3 and G-H1···C-N3 hydrogen bonds plus the εA-N7···T-O4 and the εA-N6···T-O2 distances in the modified base pair are subject to the dynamic distance constraint. These are the nonoverlapping pairs (NOP) specified in the formula for D given in Section 2.6.1.
This collective variable has the advantage of not biasing the dissociation process, i.e., not favouring a particular pathway. At the initial point of DDNA = 3.7 Å, a QM/MM molecular dynamics simulation was performed for 4000 steps, resulting in a different structure with Hg1 positioned inbetween εA20, T7 and T8, and a proton being transferred from T19 to A6 (Figure 9b). Hg2 then cross-

Structure 2
The dynamic distance method was subsequently employed again as a reaction coordinate to simulate the dissociation of DNA and to determine the corresponding free energy profile. In this case, the dynamic distance D DNA comprises all N···H hydrogen bonds of all regular base pairs together with the two N···O distances bridged by the two Hg ions in the metal-mediated base pair (Figure 1a). Thereby all the A-N1···T-H3 and G-H1···C-N3 hydrogen bonds plus the εA-N7···T-O4 and the εA-N6···T-O2 distances in the modified base pair are subject to the dynamic distance constraint. These are the nonoverlapping pairs (NOP) specified in the formula for D given in Section 2.6.1.
This collective variable has the advantage of not biasing the dissociation process, i.e., not favouring a particular pathway. At the initial point of D DNA = 3.7 Å, a QM/MM molecular dynamics simulation was performed for 4000 steps, resulting in a different structure with Hg1 positioned in-between εA20, T7 and T8, and a proton being transferred from T19 to A6 (Figure 9b). Hg2 then cross-links the bases T7 and T19 originally involved in neighbouring base pairs (Figure 9c). The time evolution of the N3···H and N3···Hg2 distances during deprotonation of T19-N3 are depicted in Figure 10, showing that as the proton leaves, its binding position approaches the Hg2 ion. The resulting structure is used as the initial structure 2 for the second dissociation path of the DNA.  three phosphate groups neutralizing the subsystem of the constrained QM/MM molecular dynamics starting from structure 1 at DDNA = 3.7 Å. The N3H proton of T19, which is transferred to A6-N1, is marked with an arrow and coloured in purple. (a) t = 0 fs, the mercury ions (yellow) are initially located within the modified base pair, while Hg2 attracts the electron cloud from T19, thus favouring deprotonation. (b) t = 140 fs, proton transfer transition state. (c) t = 200 fs, the proton has been transferred to A6 and a T19-Hg(II)-T7 base pair is formed.

Structure 3
For comparison, an analogous DNA structure without any Hg ions (termed initial structure 3) was generated by deleting the Hg ions from initial structure 1, protonating the central thymine residue T7 (Figure 1b) and neutralizing the system by the addition of three counter ions. Then a geometry optimization was performed, followed by 500 ps of NPT-MD with settings given in Section 2.3. The resulting structure 3 is shown in Figure S8.

Structure 4
To have access to an all-QM calculation, the isolated modified base pair without Hg (Figure 1b) was solvated by 168 water molecules ( Figure 11). The computational details are given in Section 2.3. The force field allowed for a 150 ps NPT equilibration simulation of the solvated base pair. This resulted in a unit cell of size 11.2 Å × 24.2 Å × 20.4 Å, corresponding to a density of 1037 kg/m 3 . This geometry will be referred to as structure 4. For comparison, an analogous DNA structure without any Hg ions (termed initial structure 3) was generated by deleting the Hg ions from initial structure 1, protonating the central thymine residue T7 (Figure 1b) and neutralizing the system by the addition of three counter ions. Then a geometry optimization was performed, followed by 500 ps of NPT-MD with settings given in Section 2.3. The resulting structure 3 is shown in Figure S8.

Structure 4
To have access to an all-QM calculation, the isolated modified base pair without Hg (Figure 1b) was solvated by 168 water molecules ( Figure 11). The computational details are given in Section 2.3. The force field allowed for a 150 ps NPT equilibration simulation of the solvated base pair. This resulted in a unit cell of size 11.2 Å × 24.2 Å × 20.4 Å, corresponding to a density of 1037 kg / m 3 . This geometry will be referred to as structure 4. Structure 4 was used as the initial geometry for an all-QM dissociation simulation of the isolated base pair without Hg ions in-between the bases.

Structure 5
Adding two Hg ions in-between the bases while removing the proton from T-N3 in structure 4 led to the initial structure for the dissociation simulation of the isolated mercurated base pair Structure 4 was used as the initial geometry for an all-QM dissociation simulation of the isolated base pair without Hg ions in-between the bases.

Structure 5
Adding two Hg ions in-between the bases while removing the proton from T-N3 in structure 4 led to the initial structure for the dissociation simulation of the isolated mercurated base pair (structure 5).

Free Energy Simulations
Having generated these five different model systems, free energy simulations were performed to analyze the dissociation path and thus, using the principle of micro-reversibility, the formation path of the DNA duplex. QM/MM MD was performed with the settings outlined in Section 2.5 and the dynamic distance as defined in Section 2.6.1.

DNA Dissociation
The calculated free energy profiles for DNA dissociation are presented in Figure 12. At the lowest energy dynamic distance of D DNA = 2.9 Å of structure 1, the average distance between the carbon atoms C1 of the deoxyribose attached to the nucleobases of the modified base pair is 10.4 Å. This is somewhat smaller than the experimental value of 11.4 Å reported for a parallel-stranded duplex composed of A:T base pair in the reversed Watson-Crick geometry [86]. Within structure 2, where T19 was deprotonated and the Hg(II) ions were bound across different base layers, the free energy minimum is found at D DNA = 3.6 Å with an average C1 ···C1 distance of 11.8 Å, which is closer to the experimental value. During the simulation of structure 2, both the NO and the NN bonding pattern are observed, with the system oscillating between the two minima (Figures 13 and 14). The first free energy minimum observed around D DNA = 3.6 Å arises because the amino groups of adenine residues A6-N6H 2 and A8-N6H 2 involved in one of the base pairs adjacent to the metal-mediated one stabilize the Hg2 ion in-between them.

DNA Dissociation
The calculated free energy profiles for DNA dissociation are presented in Figure 12. At th west energy dynamic distance of DDNA = 2.9 Å of structure 1, the average distance between th rbon atoms C1′ of the deoxyribose attached to the nucleobases of the modified base pair is 10.4 is is somewhat smaller than the experimental value of 11.4 Å reported for a parallel-strande plex composed of A:T base pair in the reversed Watson-Crick geometry [86]. Within structure here T19 was deprotonated and the Hg(II) ions were bound across different base layers, the fr ergy minimum is found at DDNA = 3.6 Å with an average C1′···C1′ distance of 11.8 Å, which is clos the experimental value. During the simulation of structure 2, both the NO and the NN bondin ttern are observed, with the system oscillating between the two minima (Figures 13 and 14). Th st free energy minimum observed around DDNA = 3.6 Å arises because the amino groups of adeni sidues A6-N6H2 and A8-N6H2 involved in one of the base pairs adjacent to the metal-mediated on abilize the Hg2 ion in-between them.    1 and 2 contain the dinuclear Hg(II)-mediated base pair, whereas structure 3 is the analogous mercury-free duplex. Structure 1 starts from a QMMM-optimized geometry, whereas structure 2 starts after QMMM-MD, and structure 3 is simulated by both classical and ab initio MD. In all simulations, the modified base pair proved to be the weak spot of the DNA duplex, opening before the other base pairs. The bonding was found to be weakest if no mercury was incorporated into the DNA as can be seen from the free energy curves in Figure 12. As the dissociation simulation of structure 3 led to a breaking of the inner-pair hydrogen bonds of the artificial base pair (see Figure S3), the artificial base pair in case of structure 3 can be considered dissociated. For structures 1 and 2, only a change of bonding patterns occurred; thus, the DNA with Hg(II) can be considered more stable. To check for a possible influence of the classical force field parameters, the dissociation MD of structure 3 was repeated by a purely classical force-field approach. This was possible because no mercury ions are present in structure 3. Ten points along the trajectory of structure 3 were sampled for 55 ps, each, by classical MD. The NN pattern can be identified after 0.6 ps, corresponding to the situation in Figure 13a). At 1.5 ps the system is in the NO pattern, depicted in Figure 13b).
Structure 1, where the calculations started from the optimized QM subsystem, emerges as the most stable structure, exhibiting the biggest energy gradient and thus the largest required force to promote the system out of its initial bonding pattern. It is followed by structure 2, which exhibits a second energy minimum along the dissociation path at DDNA = 5.5 Å. At this point, Hg1 is still located in the modified base pair, whereas Hg2 forms a diagonal T7-Hg(II)-T19 base pair involving nucleobases that were originally located in neighbouring base pairs. Such an unusual metal-mediated base pair was recently observed in a crystal structure, too [87]. This bonding pattern was enabled by deprotonation of T19, which is favoured by the presence of Hg(II) in the adjacent modified base pair. Structure 3 shows the weakest binding in the classical description. Its QM description led to lower binding energies than observed in the metal-mediated DNA duplexes. At DDNA= 8.5 Å, structure 3 opened a cavity at the modified base pair, allowing water to interact with the base pairs of the DNA duplex.
All computed free energy differences lie in the range of 100-200 kJ/mol, which are typical binding energies of DNA base pairs [88]. These calculated energies correspond to the dissociation process depicted in Figures S1-S3 for the simulations of the whole DNA and to the binding energy of the isolated base pair in water from Section 4.2. As the simulated partial dissociation in the whole DNA mainly affects the modified base pair, the calculated energies are indicative of the amount of energy required to induce the observed change of bonding patterns within the modified base pair.
Within structures 1 and 2, the Hg(II) ions did not leave the DNA during the dissociation simulation, but changed their bonding pattern. However, structure 1 exhibited a dissociation path that allowed Hg2 to get into contact with water. This dissociation path can be analysed as follows. From the initial geometry of being situated within the modified base pair, the final geometries for structure 1 at DDNA = 5.8 Å showed Hg2 in an A6-Hg(II)-εA20 environment (Figure 15a). The median distances from Hg2 to the coordinated nitrogen atoms are 2.15 Å to εA20-N6 and 2.17 Å to A6-N1. Hg1 is involved in a T7-Hg(II)-A8 pair dangling on the side of only one strand with a median distance from the coordinated nitrogen atoms to Hg1 of 2.17 Å (T7-N3) and 2.12 Å (A8-N6) ( Figure  15b). These values suggest the formation of coordinate bonds to the Hg ions. The NH2 group of A8 changes its hybridization from a sp 2 -like planar state to a sp 3 -like tetrahedral state, with Hg1 acting as the fourth binding partner. Such a hybridization shift of the exocyclic amino group of adenine has been reported before in structurally characterized metal complexes [89]. In all simulations, the modified base pair proved to be the weak spot of the DNA duplex, opening before the other base pairs. The bonding was found to be weakest if no mercury was incorporated into the DNA as can be seen from the free energy curves in Figure 12. As the dissociation simulation of structure 3 led to a breaking of the inner-pair hydrogen bonds of the artificial base pair (see Figure S3), the artificial base pair in case of structure 3 can be considered dissociated. For structures 1 and 2, only a change of bonding patterns occurred; thus, the DNA with Hg(II) can be considered more stable.
To check for a possible influence of the classical force field parameters, the dissociation MD of structure 3 was repeated by a purely classical force-field approach. This was possible because no mercury ions are present in structure 3. Ten points along the trajectory of structure 3 were sampled for 55 ps, each, by classical MD.
Structure 1, where the calculations started from the optimized QM subsystem, emerges as the most stable structure, exhibiting the biggest energy gradient and thus the largest required force to promote the system out of its initial bonding pattern. It is followed by structure 2, which exhibits a second energy minimum along the dissociation path at D DNA = 5.5 Å. At this point, Hg1 is still located in the modified base pair, whereas Hg2 forms a diagonal T7-Hg(II)-T19 base pair involving nucleobases that were originally located in neighbouring base pairs. Such an unusual metal-mediated base pair was recently observed in a crystal structure, too [87]. This bonding pattern was enabled by deprotonation of T19, which is favoured by the presence of Hg(II) in the adjacent modified base pair. Structure 3 shows the weakest binding in the classical description. Its QM description led to lower binding energies than observed in the metal-mediated DNA duplexes. At D DNA = 8.5 Å, structure 3 opened a cavity at the modified base pair, allowing water to interact with the base pairs of the DNA duplex.
All computed free energy differences lie in the range of 100-200 kJ/mol, which are typical binding energies of DNA base pairs [88]. These calculated energies correspond to the dissociation process depicted in Figures S1-S3 for the simulations of the whole DNA and to the binding energy of the isolated base pair in water from Section 4.2. As the simulated partial dissociation in the whole DNA mainly affects the modified base pair, the calculated energies are indicative of the amount of energy required to induce the observed change of bonding patterns within the modified base pair.
Within structures 1 and 2, the Hg(II) ions did not leave the DNA during the dissociation simulation, but changed their bonding pattern. However, structure 1 exhibited a dissociation path that allowed Hg2 to get into contact with water. This dissociation path can be analysed as follows. From the initial geometry of being situated within the modified base pair, the final geometries for structure 1 at D DNA = 5.8 Å showed Hg2 in an A6-Hg(II)-εA20 environment (Figure 15a). The median distances from Hg2 to the coordinated nitrogen atoms are 2.15 Å to εA20-N6 and 2.17 Å to A6-N1. Hg1 is involved in a T7-Hg(II)-A8 pair dangling on the side of only one strand with a median distance from the coordinated nitrogen atoms to Hg1 of 2.17 Å (T7-N3) and 2.12 Å (A8-N6) (Figure 15b). These values suggest the formation of coordinate bonds to the Hg ions. The NH 2 group of A8 changes its hybridization from a sp 2 -like planar state to a sp 3 -like tetrahedral state, with Hg1 acting as the fourth binding partner. Such a hybridization shift of the exocyclic amino group of adenine has been reported before in structurally characterized metal complexes [89]. The dissociation MD of structure 2 proceeded to a final value of DDNA = 10.4 Å, when the quantum-mechanically treated Hg1 came into contact with classically treated water. At this point, Hg2 is involved in a T7-Hg(II)-T19 pair, with equal Hg-N3 distances of 2.15 Å on either side ( Figure  16a). In this structure, Hg1 is contained in an εA20-Hg(II)-A8 pair with median distances from the nitrogen atoms coordinated to Hg(II) of 2.20 Å (A-N1) and 2.18 Å (εA-N3), which are close to the distances found in the εA20-Hg(II)-A6 base pair in structure 1. All distances suggest coordinate bonding of the Hg(II) ions (Figure 16b). The only migration of an Hg ion from the DNA inner core to the water solvent occurred in structure 1. The inverse of the pathway out of the DNA duplex can be considered a possible pathway into the duplex, using the principle of microscopic reversibility. The transition out of the DNA duplex is mediated by the exocyclic amino group of A8-N6H2 from the layer below the metal-mediated base pair. Hg1 approaches the NH2 group, thus leading to an sp 3 -like hybridization of its nitrogen atom. Being attached to this moiety, which can rotate due to its single bond to the adenine, Hg1 can rotate out of the DNA core, while remaining attached to this group (Figure 17). The dissociation MD of structure 2 proceeded to a final value of D DNA = 10.4 Å, when the quantum-mechanically treated Hg1 came into contact with classically treated water. At this point, Hg2 is involved in a T7-Hg(II)-T19 pair, with equal Hg-N3 distances of 2.15 Å on either side (Figure 16a). In this structure, Hg1 is contained in an εA20-Hg(II)-A8 pair with median distances from the nitrogen atoms coordinated to Hg(II) of 2.20 Å (A-N1) and 2.18 Å (εA-N3), which are close to the distances found in the εA20-Hg(II)-A6 base pair in structure 1. All distances suggest coordinate bonding of the Hg(II) ions (Figure 16b). The dissociation MD of structure 2 proceeded to a final value of DDNA = 10.4 Å, when the quantum-mechanically treated Hg1 came into contact with classically treated water. At this point, Hg2 is involved in a T7-Hg(II)-T19 pair, with equal Hg-N3 distances of 2.15 Å on either side ( Figure  16a). In this structure, Hg1 is contained in an εA20-Hg(II)-A8 pair with median distances from the nitrogen atoms coordinated to Hg(II) of 2.20 Å (A-N1) and 2.18 Å (εA-N3), which are close to the distances found in the εA20-Hg(II)-A6 base pair in structure 1. All distances suggest coordinate bonding of the Hg(II) ions (Figure 16b). The only migration of an Hg ion from the DNA inner core to the water solvent occurred in structure 1. The inverse of the pathway out of the DNA duplex can be considered a possible pathway into the duplex, using the principle of microscopic reversibility. The transition out of the DNA duplex is mediated by the exocyclic amino group of A8-N6H2 from the layer below the metal-mediated base pair. Hg1 approaches the NH2 group, thus leading to an sp 3 -like hybridization of its nitrogen atom. Being attached to this moiety, which can rotate due to its single bond to the adenine, Hg1 can rotate out of the DNA core, while remaining attached to this group ( Figure 17). The only migration of an Hg ion from the DNA inner core to the water solvent occurred in structure 1. The inverse of the pathway out of the DNA duplex can be considered a possible pathway into the duplex, using the principle of microscopic reversibility. The transition out of the DNA duplex is mediated by the exocyclic amino group of A8-N6H 2 from the layer below the metal-mediated base pair. Hg1 approaches the NH 2 group, thus leading to an sp 3 -like hybridization of its nitrogen atom. Being attached to this moiety, which can rotate due to its single bond to the adenine, Hg1 can rotate out of the DNA core, while remaining attached to this group ( Figure 17). Figure 16. Final geometry of structure 2 at DDNA = 10.4 Å with no water displayed. (a) Hg2 is bound between T19-N3 from the base pair originally adjacent the modified base and T7-N3 from the modified base pair. T19 transferred its N3-H proton to A6-N1 during the equilibration MD due to the presence of Hg(II). (b) Final geometry of structure 2 after thermodynamic integration at DDNA = 10.4 Å with water not displayed. Hg1 is bound by A8-N1 from the base pair originally below the modified base and εA20-N7.
The only migration of an Hg ion from the DNA inner core to the water solvent occurred in structure 1. The inverse of the pathway out of the DNA duplex can be considered a possible pathway into the duplex, using the principle of microscopic reversibility. The transition out of the DNA duplex is mediated by the exocyclic amino group of A8-N6H2 from the layer below the metal-mediated base pair. Hg1 approaches the NH2 group, thus leading to an sp 3 -like hybridization of its nitrogen atom. Being attached to this moiety, which can rotate due to its single bond to the adenine, Hg1 can rotate out of the DNA core, while remaining attached to this group ( Figure 17). The distance from Hg1 to A8-N6 is shown in Figure S9. During the rotation, the distance oscillates around at a mean value of 2.2 Å, while after the rotation, the Hg ion is coordinated by T7-O4 and εA20-N7 and the distance from Hg1 to A8-N6 reaches 2.8 Å, suggesting that the bond to the NH 2 group is broken.
It is tempting to speculate that the relative location of the neighbouring adenine residue with respect to the εA:T mispair is important for the incorporation of the second Hg(II) ion into the dinuclear base pair. This would be in good agreement with the observation that only a mononuclear Hg(II)-mediated base pair is formed between εA and T when the complementary sequences are arranged in an antiparallel-stranded manner, i.e., when the relative orientation of the neighbouring base pairs is different (see Supporting Information for details).

Dissociation of the Isolated Base Pair
The isolated εA-Hg(II) 2 -T base pair with Hg(II) ions incorporated (structure 5) exhibited a binding free energy of about 160 kJ/mol, while the isolated non-metal-modified base pair εA:T (structure 4) separated without any additional energy in a water environment ( Figure 18). The latter can be explained by the fact that water molecules approaching from above or below have the propensity to form strong hydrogen bonding with the individual bases of the εA-:T pair, thus favouring dissociation (Figure 19b). Within the DNA helix, this effect is compensated by the hydrophobic inner DNA core, which blocks access of water from above or below, thus leading to a substantial binding energy of about 100 kJ/mol (MM) and 125 kJ/mol (QM/MM) for the non-metal-mediated DNA duplex ( Figure 12).
Along the dissociation path of structure 5, the dinuclear base pair splits up such that each base keeps an Hg(II) ion attached to it. Hg1 bonded to T-N3 coordinates an additional OH − formed from a nearby water molecule. The former binding site of Hg1 at εA is filled up by a hydrogen-bonded water molecule, pointing its oxygen atom towards Hg2. At DISO = 8.6 Å, the attractive Coulomb interactions between the bases are screened by water molecules, facilitating further dissociation ( Figure 20).

Conclusions
Free energy QM and QMMM molecular dynamics simulations have been performed of an experimentally established parallel-stranded DNA duplex bearing a T-Hg(II)2-εA base pair. Two main bonding patterns have been identified for the mercury(II)-mediated base pair. For the most stable structure 2, the C1′···C1′ distance of the carbon atoms involved in the glycosidic bonds is close to the experimental value derived from a related experimental structure [86]. While the isolated nucleobases T and εA do not form a base pair in the absence of Hg(II), the inclusion of Hg(II) leads to a stable dinuclear metal-mediated base pair. If neighbouring base pairs are considered in the calculations, the amino groups of their adenine nucleobases influence the bonding pattern by binding to the Hg(II) from above and below. This unexpected structural insight shows that the precise location

Conclusions
Free energy QM and QMMM molecular dynamics simulations have been performed of an experimentally established parallel-stranded DNA duplex bearing a T-Hg(II)2-εA base pair. Two main bonding patterns have been identified for the mercury(II)-mediated base pair. For the most

Conclusions
Free energy QM and QMMM molecular dynamics simulations have been performed of an experimentally established parallel-stranded DNA duplex bearing a T-Hg(II) 2 -εA base pair. Two main bonding patterns have been identified for the mercury(II)-mediated base pair. For the most stable structure 2, the C1 ···C1 distance of the carbon atoms involved in the glycosidic bonds is close to the experimental value derived from a related experimental structure [86]. While the isolated nucleobases T and εA do not form a base pair in the absence of Hg(II), the inclusion of Hg(II) leads to a stable dinuclear metal-mediated base pair. If neighbouring base pairs are considered in the calculations, the amino groups of their adenine nucleobases influence the bonding pattern by binding to the Hg(II) from above and below. This unexpected structural insight shows that the precise location of a metal ion within a metal-mediated base pair is not necessarily identical to that in the ideal geometry of the isolated base pair but may also depend on the identity of the adjacent canonical base pairs.
In the complete DNA duplex, the artificial base pair proved to be the weak point in both classical and quantum mechanics simulations if no Hg(II) is incorporated. Upon to the incorporation of Hg(II), the base pair is strengthened. In the dissociation simulations of the Hg(II)-bound DNA duplex, a shift in bonding patterns is observed, with the Hg(II) ions engaging in the original T-Hg(II) 2 -εA base pair, in an inter-planar T-Hg(II)-T base pair, or in an unprecedented A-Hg(II)-εA base pair. Irrespective of the precise coordination environment of the Hg(II) ions, the Hg(II)-containing DNA is predicted to be more stable than its non-metalated counterpart.
During one simulation, a mercury ion was transferred from the inner DNA region to the surrounding water, where the amino group of an adjacent adenine moiety acted as a turnstile to transport the mercury(II) ion out of the DNA. Due to the principle of microscopic reversibility, this adenine residue can be postulated to act as a gateway for Hg(II) during the formation of the metal-mediated base pair. This indicates another important contribution of the oligonucleotide sequence on the mechanism of metal-mediated base pair formation, which will need to be investigated in detail in future experiments. This is in good agreement with the experimental observation that in a different sequence context, only one Hg(II) is incorporated into the εA:T pair. The present results may thus prove useful in guiding the future design of Hg-modified DNA.
As Hg(II) has a linear coordination environment, the findings of this paper concerning the structural arrangement can only be transferred to metal ions which coordinate linearly, which is not the case for the other metal ions mentioned in the introduction, except for Ag(I). Similarly, they are not applicable to Ag(I), because as a monovalent cation it is less capable of deprotonating aqua ligands to hydroxido ligands. However, the observed turnstile mechanism might be relevant to the incorporation of other metal cations into DNA, as they can bind to the lone electron pair of the adenine amino group.
Aside from the chemical insights gained, it is worth noting that the dynamic distance constraint proved to be a useful method to investigate bonding patterns in systems with competing (hydrogen) bonds-not only for quantitative free energy simulations, but also for constructing DNA geometries with predefined hydrogen bonding patterns. Since this is a common issue in the investigation of solvated biomolecules, the current approach should also be helpful for many other systems.
Supplementary Materials: The following are available online, Figure S1: Final geometry of structure 1, Figure S2: Final geometry of structure 2, Figure S3: Final geometry of structure 3, Figure S4: Lagrange multipliers for path 1, Figure S5: Lagrange multipliers for path 2, Figure S6: Lagrange multipliers for path 3, Figure S7: Energies and end-to-end distances of optimized geometries for different functionals, Figure S8: DNA structure without Hg, Figure S9