All-atom Molecular Dynamics simulations of spin labelled double and single-strand DNA for EPR studies

Please note that technical editing may introduce minor changes to the text and/or graphics, which may alter content. The journal’s standard Terms & Conditions and the ethical guidelines, outlined in our author and reviewer resource centre, still apply. In no event shall the Royal Society of Chemistry be held responsible for any errors or omissions in this Accepted Manuscript or any consequences arising from the use of any information it contains. Accepted Manuscript


Introduction
The internal dynamics and conformational variability of DNA are known to be crucial to its biological functions including interaction with proteins and expression of the genetic code. 1,2 Because of the semi-flexible nature of DNA its dynamic regime consists of multiple contributions from breathing, bending and twisting modes, as well as groove fluctuations of the helix. Analysis of such motions is challenging and has been the subject of extensive studies. Various experimental techniques including nuclear magnetic resonance (NMR), 3 fluorescence, 4 dynamic light scattering 5 and Fourier transform infrared difference spectroscopy 6 amongst others have been utilised to probe different motional components. Electron Paramagnetic Resonance (EPR) with introduced spin labels is a particularly suitable spectroscopic technique to study the dynamics and conformational changes in complex bio-molecular systems such as DNA. EPR, being a 'fast' magnetic technique, is able to directly resolve molecular dynamics and structural changes on the sub-nanosecond time scale. [7][8][9][10][11] Furthermore, due to recent advances in spin labelling methodologies (e.g. the use of click chemistry) an increasingly wide range of spin labels are being developed for DNA studies by EPR, allowing greater sequence selectivity. [12][13][14][15][16][17] However, the need to introduce a synthetic spin label with its own internal dynamics, combined with the potential to influence the local DNA structure might significantly complicate the analysis of EPR data. Additionally, global tumbling motions of relatively short DNA fragments employed in many model experimental studies are significant on the EPR timescale, thus requiring sophisticated spectral modelling to interpret the results. 7,9 In the case of nitroxide labels attached via flexible tethers this typically requires the fitting of EPR spectra with multiple adjustable parameters using the slowly relaxing local structure (SRLS) model of Freed and co-workers 7 with global tumbling rates based on hydrodynamic theory. 18,19 Spin labels such as the quinolonyl derived Q (Fig. 1) 10 or alkyne tethered labels such as the cytosine derived C* (Fig. 1) 20 (or the analogous thymine derived T* 21 ) were the first ones developed by Robinson and co-workers for DNA studies using EPR. [8][9][10][11]20,22 These labels were assumed to report accurately on the motions of the host DNA, with the spectra typically interpreted through the application of either diffusion of a rigid cylinder 10 or weakly bending rods 9 models. In particular the Q spin label is assumed to exhibit negligible internal motion, as well as being highly thermally stable. However it has no natural analogue and requires a synthetic 2-aminopurinyl complementary base pair, P ( Fig. 1), with the effect on the local dynamics and structure unclear from EPR alone. 10 C* represents an alkyne-linked spin label covalently attached to conventional base pairs and thus is assumed to be less structurally perturbing. 10,11 Although such labels have been demonstrated to be highly sensitive to the dynamics of DNA duplexes, they exhibit internal rotation about the alkyne linkage that potentially compromises the structural and dynamic information obtained from the analysis of EPR spectra. 11 Over the past decade novel approaches have been developed allowing for the prediction of motional Continuous wave (CW) EPR spectra from molecular dynamics (MD) simulations by either indirect [23][24][25] or direct propagation [25][26][27][28][29] calculation methods. MD-EPR prediction techniques can greatly simplify the interpretation and analysis of EPR experimental spectra, and hence provide unambiguous conclusions about molecular order and motions. They have been successfully applied to spin labelled proteins 24,[27][28][29][30][31] and soft matter systems such as liquid crystals, both thermotropic and lyotropic. [32][33][34][35] Fully atomistic MD simulations have already provided significant insights into the sequence-dependent flexibility of DNA. [36][37][38][39][40] In particular it is becoming increasingly apparent that DNA sequences can adopt a wide variety of conformations, depending on the chemical environment, and that generally the structure of DNA should be considered in terms of conformational ensembles. 40,41 In the past modelling studies on DNA have been hindered by the lack of generating sufficiently long MD trajectories required for representing adequately the helix properties on different timescales. 41 Recent advances in computing power and refined force fields, such as the recently parameterised parmbsc1, 42 allow the accurate description of conformational behaviour with trajectories of up to 10 ms currently achievable at least for relatively short DNA chains. 40 The purpose of this paper was to achieve the following. Firstly, given that the EPR spectra are highly sensitive to the motions and order of the spin probes, simulation of EPR line shapes from the results of MD provides an ultimate test bed for the force fields currently employed to model DNA and also RNA structures. Additionally we have analysed the perturbing effects from the presence of spin probes on the DNA local structure. Secondly, the use of fully atomistic MD combined with MD-EPR simulation methodology allowed us to provide detailed quantitative analysis of different motional contributions, both internal and global associated with the spin probe's and DNA motions, respectively, into the dynamics of the Q and C* labels in both duplex and single-stranded DNA fragments. Finally, the application of fully atomistic MD simulations for modelling single strand and duplex labelled DNA complexes allowed us to test the validity of previously employed simplified models for the simulation and interpretation of EPR spectra that are based on the application of so-called Model-Free (M-F) approach. The M-F approach assumes that the global and local motions are independent. In addition, in many simulation strategies employed previously the local dynamics of the label was assumed to be in the so-called fast motional regime on the EPR timescale that justifies the partial averaging of the magnetic tensors A and g. 22,43 The advantage of performing MD is that the statistically averaged parameters employed in such models can be readily calculated from MD trajectories and used in the simulation of the EPR spectra. In our case both global (DNA) and local (nitroxide spin probe) motional contributions have an impact on CW X-band EPR spectra thus making the reported DNA systems an ideal test bed for simplified models.

Molecular dynamics modelling
Initial 14-mer DNA configurations were constructed using the analysis module of the w3DNA web server developed by Olson and co-workers. 44 Conventional base pairs were modelled using the parmbsc1 force field. 42 Parameters for both spin labels have been generated in analogy with other spin probes developed by us previously, initially in the General AMBER force field (GAFF) 45 and subsequently adapted to the parmbsc1 format. Quantum chemical calculations of Q and C* were performed with the Gaussian 09 software package. 46 Force field parameters for the new atom types of the nitroxide moieties (the unsaturated carbon atoms of the nitroxide ring, the saturated carbon atoms of the nitroxide ring, the nitrogen and the oxygen) were taken from a combination of geometry optimization calculations in the gas phase and previous calculations. Equilibrium bond lengths and angles were taken directly from minimized energy structures. Force constants were interpolated using the reference values in the AMBER99 force field 47 and the quantum mechanical calculations of Barone and co-workers. 48,49 In the case of the C* label two dummy atoms in the triple bond (see Fig. S1 of ESI †) were introduced in order to define the torsional rotation between the cytosine and nitroxide rings. The parameters for the dihedral angle were determined by fitting to the QM potential energy scan. Partial charges for spin labels were calculated using a multi-conformational Restrained Electrostatic Potential (RESP) fit 47 at the HF/6-31G* level of theory. The calculated force field parameters and partial charges for the Q and C* probes are included in the ESI † ( Fig. S1-S3 and Tables S1-S5). The SPC/E water model 50 was used with the Smith-Dang ion parameters for sodium counter ions 51 with the systems compositions given in Table 1. All MD calculations were performed using the AMBER 14 52 software package.
An NPT ensemble was maintained at a pressure of 1 atm using the Berendsen algorithm 53 with a coupling constant of 5 ps. The SHAKE algorithm 54 was used to maintain the hydrogen bond lengths. The centre of mass motion was removed every 20 ps to limit the build-up of translational kinetic energy, allowing for a time step of 2 fs to be used. Long range electrostatic interactions were accounted for using the Particle Mesh Ewald 55 method with a cut-off of 10 Å. Systems were equilibrated for 700 ns prior to production runs of 700 ns. DNA conformational analysis was performed using the software from the w3DNA web server. 44

Rotational autocorrelation functions
The autocorrelation function of each vector, l, associated with either the magnetic axes of the nitroxide head group of each spin label or the principal axes of DNA, can be calculated from an MD trajectory using the following expression: 34 where P 2 (x) is the second order Legendre polynomial: and the bracket in (1) denotes the average taken over time.
The local motional component of the probes and the global tumbling of the DNA were separated as follows. The global DNA motion was approximated by the dynamics of the principal axes of the tensor of inertia of each DNA fragment. Because of the effectively nearly cylindrical shape of DNA the X/Y principal components of the rotational tensor become poorly resolved for different moments of time. Thus they were estimated using the vector between the C1 0 (glycosidic link) atoms of two complementary base pairs. The local motions of the probes (motions in the DNA fixed frame) were extracted using a mass-weighted RMS structural fit with the Ptraj module of AmberTools. 56 In the M-F framework the motional and order parameters were obtained from the fitting of autocorrelation functions for the local and global motions using eqn (3) and (4), respectively, derived using the M-F formalism of Lipari and Szabo. 34,57 According to these authors, for Markovian type motions the correlation function for internal dynamics can be generally expressed as a series of exponentials. 57 Indeed, for many molecular systems the autocorrelation function for local dynamics can be well approximated by three different motional contributions and one time independent term S 2 which is the square of the generalised local order parameter (eqn (3)). In eqn (3) index i corresponds to x, y, or z magnetic axes and w i are the weighting factors associated with the motional contributions. The effective local correlation time is calculated as t i = w 1 t i1 + w 2 t 2i + w 3 t 3i . Eqn (4) corresponds to the model of the free axial rotational diffusion of an DNA fragment where the components of the rotational diffusion tensor D > and D 8 are related to correlation times according to t i = 1/(6D i ). 34,57 Direct prediction of EPR spectra from MD trajectories A previously reported trajectory-based method that employs the numerical solution of the Stochastic Liouville Equation (SLE) in the Langevin form for the spin density matrix has been used for the simulation of the EPR line shapes. 25,28 A program developed and described previously by one of us 26 has been employed. The relatively long MD trajectories generated in this work allowed the simulation of CW EPR spectra directly by propagation of the spin density matrix along the entire sampling time without further approximations. In the program, single MD trajectories are used to calculate the variation in time of the averaged transverse magnetisation and, eventually, the EPR line shapes. 26 At each time increment the propagation of the density matrix was carried out in Hilbert space using both eigenvalues and eigenvectors of the Spin-Hamiltonian as reported previously. 29 Statistical averaging was achieved by the ''sliding time window technique'', allowing the use of single MD trajectories for predicting EPR line shapes. 25,28 At the X-band the spectrum is dominated by the anisotropic hyperfine coupling A tensor and, under the condition where the relevant elements of the g L and A L tensors in the laboratory frame are determined from the principal values for electron g and hyperfine coupling A tensors 28 in the frame of the nitroxide using the following Cartesian transformations R(O(t)) In (5) o 0 , b, h and B are the resonance frequency, Bohr magneton, Planck's constant and magnetic field respectively and m = AE1, 0. Note that there is no hyperfine contribution to the central line (m = 0). The orientational history of the magnetic axes in the fixed frame of the simulation box is calculated and processed. The EPR spectral line shapes of nitroxide spin labels are determined entirely by the variation with time of two angles that define the orientation of the applied magnetic field to the principal axis of the nitroxide group. The z axis of the nitroxide ring (coincident with the direction of the p z -orbital of N) is calculated from the crossproduct of the unit vectors of two N-C bonds of the nitroxide ring (see Fig. 1). 28,33 The x axis is calculated as a projection vector of the N-O bond on the nitroxide plane (defined by the C-N-C atoms) and the y axis is taken as a cross-product of the z and x vectors.

Prediction of EPR spectra using a combination of MD simulations and the Model-Free approach
Most recently we have reported the simulation of EPR line shapes using the M-F approach with the motional parameters extracted from the MD trajectories of lyotropic liquid crystals doped with a paramagnetic spin probe in a range of different aggregation states, namely, micro-aggregates, micelles, rods and lamellar states. 34 EPR line shapes were simulated by solving the SLE in the Fokker-Planck (F-P) form. 58,59 Thus a combined MD-EPR methodology allowed us to test directly the validity of the application of the M-F approach coupled with the partial averaging of magnetic tensors due to fast local motions 34,60,61 to systems with complex multi-component molecular dynamics by comparing the resulting EPR spectra to those simulated directly from MD trajectories and also to the experimental ones.
The advantage of using all-atom MD simulations is that both the motional and order parameters employed in the M-F approach can be readily calculated from the MD outputs. In most of lyotropic aggregate states reported in ref. 34 the internal dynamics of the probe was sufficiently fast for partial averaging of the magnetic tensors. In addition, because of the bulk structures the orientations of the principal components of the partially averaged magnetic tensors of the spin probe were well defined in each of the aggregate states. In the spin labelled DNA fragments the orientation of the partially averaged principal tensor components is not a priori obvious making them ideal systems to test the application of the M-F approach combined with the partial averaging of A and g by fast local motions in a general case. Both the values and the orientations of partially averaged principal components of the magnetic tensors were calculated from MD using the following procedure. Firstly, both tensors were averaged according to the following equations: where the averages are performed on the products of the projection cosines l ij of the three magnetic axes. This was followed by the diagonalization of both hAi and hgi by performing the Cartesian transformation % R whose columns are the eigenvalues (projection cosines) that define the orientations of the averaged principal tensor components in the DNA fixed frame.
For the g tensor equations similar to (6) and (7) are used. Finally, the rotational angles that transform the principal orientations into the DNA fixed frame are calculated in accordance with: Transformation from the DNA fixed frame into the frame defined by the directions of the principal components of the partially averaged magnetic tensor is shown in Fig. S4 of ESI. † Our results showed that % A and % g were nearly collinear and as such only the rotational angles for % A were used in the simulations of EPR. They are reported in Table 4 for all systems. Together with the partially averaged values of % A and % g and the global rotational diffusional coefficients (correlation times), obtained from the fitting of the autocorrelation function using eqn (4), they were used for the simulation of the EPR spectra using M-F approach.

Calculation of magnetic parameters
Magnetic parameters, namely the principal values of both A and g tensors, were calculated using the B3LYP functional with the N07D 62,63 basis set as implemented in the Gaussian 09 software package. 46 Structures were first optimised in the gas phase followed by optimisation in water using the polarizable continuum solvent scheme. Hyperfine coupling tensor A is calculated by combining contributions from the Fermi contact and anisotropic spin dipole coupling. 64,65 Importantly, the calculation of magnetic parameters from first principles makes the entire MD-EPR simulation approach fully predictive.

Results and discussion
The principal values of the g and A tensors calculated by DFT methods indicate that the magnetic parameters of both C* and Q probes are close to each other (Table 2). Excellent agreement is found between the predicted g and A values and those obtained from EPR measurements on the immobilised Q spin label (g xx = 2.0076, g yy = 2.0061, g zz = 2.0028; A xx = 6.68 G, A yy = 5.41 G, A zz = 33.90 G) available from the literature, 10 confirming the accuracy of the B3LYP/N07D model for calculation of these parameters for second row elements and nitroxide radicals. 62,63,66 The calculated values were, therefore, used without scaling for the direct prediction of EPR spectra.

Effect of attached spin label on duplex DNA structure
Within the duplex DNA both spin labels Q and C* were found to remain base-stacked for the entire MD trajectories. The calculated RMSD of the DNA backbone with terminal base pairs removed (C* 1.72 Å, Q 2.16 Å) gives values in agreement with those typically reported for B-DNA (1.6-2.2 Å). 40 The small fluctuations observed for this parameter in labelled DNA fragments confirm the structural stability of both duplexes over the course of the production trajectory (see Fig. S5 of ESI †).
For the C* label, which is a modified cytosine base, comparison with the corresponding unlabelled sequence has been performed by calculating several geometric parameters of the helical base pair step using the analysis tools of the 3DNA program. 44 The results are presented in Table S6 of ESI. † Good agreement between the twist, roll, slide and rise is observed between the average structures for the labelled and unlabelled sequences denoted by the red and blue lines, respectively, in Fig. S6 of ESI. † The shift and tilt of the average structures display slightly lower agreement, however in each case the geometry of selected random frames shows this to be within the expected deviation of the DNA structure caused by motions. Sugar pucker angles (phase), which are closely related to the backbone conformation, are within one degree difference between the unlabelled and labelled sequences (Table S6 of ESI †). Additionally, the average amplitude values are within the range of 251-451, as observed for standard B-DNA structures. 67 Both the calculated structural parameters and the RMSD results unambiguously confirm that the C* label has a negligible perturbing effect on the geometry and flexibility of the host sequence, in agreement with conclusions based on experimental studies. 11 Since the Q probe and its complementary base pair P both have no natural equivalents to draw a comparison with, calculation of their geometric parameters have not been attempted.
Contributions to the dynamics of spin labels in duplex DNA fragments The re-orientational autocorrelation functions of the magnetic axes of Q label are presented in Fig. 2a. The autocorrelation functions of the magnetic axes of Q obtained by excluding the global DNA tumbling are presented in Fig. 2b. The motional and order parameters obtained from the fitting of autocorrelation functions are given in Table 3 (for all DNA models reported in this paper the fitted curves are presented in the ESI †). Due to the rigidity of Q, in duplex DNA the local re-orientational motions of all three magnetic axes demonstrate a high degree of order (see Table 3). Among the three the y axis exhibits the greatest local order with the local motions arising mainly from tilting about this axis. A less pronounced twisting motion about the x axis causes rotation of y and z. The z axis exhibits the lowest order. The results of the fitting of the autocorrelation functions of both the local motions of the probe and the global DNA tumbling (Fig. 2c) are presented in Table 3. The effective correlation times for local motions are calculated from the three motional contributions obtained from the fitting of the relevant autocorrelation functions using eqn (3). The results of the fitting are given in Table S7 of ESI. † Fitting of the autocorrelation function for the principal axes of the labelled DNA duplex with  an axial model (eqn (4)) yields the correlation times for the global diffusion, namely, t > = 6.21 ns and t 8 = 2.99 ns. Since the local motion of Q is highly restrained (S i 4 0.9 for all three magnetic vectors), the autocorrelation functions of the magnetic axes ( Fig. 2a) closely resemble those of the principal DNA axes (Fig. 2c), confirming that the dynamics of Q in duplex DNA adequately represents the global motion of the entire DNA fragment. 10 Notably, t > E 2t 8 , as would be expected for a 14-mer fragment where the long cylinder axis is roughly twice the length of the short axis. 11 Both components of the global motion are approximately 1.5 times faster than those reported by Okonogi et al.
(t > = 9.89 ns, t 8 = 4.43 ns) 8 using the hydrodynamic theory of Tirado and de la Torre for a duplex DNA of this size. 18,19 The principal difference is that the hydrodynamic theory treats DNA as a rigid rod, 18 whereas in our all-atom MD simulations DNA flexibility is naturally present and has an impact on the global motional rates. The value of t > is close to that reported by Miller et al. using an isotropic rotational diffusion model for the fitting of the EPR spectrum (t Iso = 7.6 ns). 10 The EPR spectrum predicted directly by propagation of the MD trajectory ( Fig. 3 red line) is found to be in good agreement with the experimental one (black line) reported by Miller et al. 10 Such a validation against a highly sensitive spectroscopic technique confirms the reasonable accuracy of the MD model employed in this work. In addition, we have performed a simulation of the EPR spectrum using the M-F approach with the motional parameters extracted from MD (Tables 3 and 4), assuming an axially symmetric rotational diffusion model of the rigid rod of the DNA fragment and partial averaging of magnetic parameters of the probe. The resulting spectrum, simulated using Easyspin software, 68 is presented as the blue line in Fig. 3 and found to be in good agreement with the line shape predicted directly and completely from the MD trajectory. Firstly, this suggests that the internal motions of the Q probe are highly restrained and thus have a negligible contribution to the overall motion and that Q indeed serves as an adequate reporter of the DNA motions. Secondly, this also confirms the validity of the application of the M-F approach to duplex DNA with the Q label. The apparent absence of the split in the low field region of the lineshape by both prediction methods is attributed to a slight overestimation of the global rotation diffusion of the DNA duplex by the SPC/E water model resulting in a somewhat reduced value of global t > . This is in agreement with the results from the fitting of the EPR spectrum reported previously. 10 The autocorrelation functions calculated from MD for the total and local motions of the C* spin probe and the dynamics of the principal rotation axes of the labelled DNA duplex are presented in subpanels (a), (b) and (c), respectively, of Fig. 4. Similar to the Q label, C* demonstrates the complex multicomponent dynamics with the very fast decay on the 1-100 ps time scale attributed to the local motion of the probe. The situation is however principally different from Q in several aspects. Firstly, autocorrelation functions with the global motion excluded (Fig. 4b) confirm that the rotation and bending of the alkyne tether leads to C* experiencing a considerably lower local order than Q for all three magnetic vectors. Secondly, rotation along the tether is associated with the most prominent re-orientational motion leading to the highest order parameter calculated for the x magnetic axis. However, measurement of the dihedral angle between the plane of the cytosine base and the nitroxide ring of the C* label (Fig. S7 of ESI †) indicates that, as has been inferred in the study of the related T* label, rotation about the alkyne linkage is not free.
This confirmed the previous assumption that due to a short tether the nitroxide group is effectively trapped in the major groove 11 (Fig. 4). Rotation along the single bond of the tether leads to the averaging of the A zz and A yy magnetic principal values of the nitroxide moiety resulting in a narrower predicted EPR line shape at 293 K (Fig. 5 top red curve) compared to the one corresponding to the Q label at the same temperature, confirming the higher level of rotational flexibility of the C* label in duplex DNA.
The EPR spectrum simulated using the M-F approach with the motional parameters extracted from MD (Tables 3 and 4) is presented as the top blue line in Fig. 5. It has reasonably good agreement with the one simulated directly and completely from  the relevant MD trajectory confirming the validity of the M-F approach for this system. Indeed, according to the calculated effective correlation times shown in Table 3 the partially restrained local dynamics of C* remains in the fast motional regime (o1 ns) at 293 K. Previously it has been shown that the EPR spectra of the closely related thymine derived spin label T* could be fitted with a simplified model which assumes that the fast axial local motion (t L o 1 ns 22 ) partially averages the magnetic parameters. 43 Robinson and co-workers using both T* and the double alkyne-bridged analogue T** demonstrated that the location of the nitroxide group in the major groove hinders local motions for alkyne-bridged labels in duplex DNA. 11 In order to compare the predicted EPR spectrum with the experimental one available from the literature we have performed an additional MD simulation at 273 K. 20 The results are compared as bottom lines in Fig. 5. A reasonable agreement between the predicted directly from MD (red line) and experimental (solid black line) EPR line shapes is observed. As in the case with the Q label (Table S7 of ESI †), fitting of the autocorrelation function of the local dynamics of the magnetic axes for the C* label required a minimum of three motional components on the 10 ps, 100 ps and 1-5 ns timescales (Table S9 of ESI †). The value of the S x local order parameter for C* (0.80 at 273 K) was found to be very close to that reported by Fischhaber (0.77 at 273 K) 20 and decreases slightly with increasing temperature (Table 3). Simulation using the M-F approach with partially averaged local magnetic tensors is presented as a blue line and is in good agreement with both the spectra predicted using the direct propagation method and the experimental one. Some discrepancy (narrowing down of some spectral features) is attributed to the presence of the slow motional mode in the local dynamics of C* (third component in the fitting of autocorrelation functions, see Table S9, ESI †) which is outside the fast motional regime, the condition that is required for the use a Since the principal axes of partially averaged % A and % g tensors are nearly collinear only the rotational angles corresponding to % A are shown. Angles are given in degrees.  of the partial averaging of the magnetic tensors. As a result, the slow motional contribution becomes underestimated in the simulated spectrum. Apparently, such effects are less pronounced when the local order parameter is high as evident from the agreement between the simulated EPR spectra by two methods of the DNA duplex with the Q label. Also, as expected, there was little difference between the motional parameters extracted for the global dynamics of both DNA duplexes labelled with Q and C* (see Fig. 2c, 4c and Table 3).
Contributions to the dynamics of spin labels in single-strand DNA fragments The autocorrelation functions calculated from MD for the total and local motions of both probes and the dynamics of the principal rotation axes of the single-strand DNA structures are presented in subpanels (a), (b) and (c), respectively, of Fig. 6. The subpanels on the left and on the right correspond to the strands labelled with Q and C* probes, respectively. In the single-strand sequences folding of the outer base pairs and temporary disruption of stacking of individual bases both enable a greater degree of local motional flexibility for both Q and C* labels than in duplex DNA. This is confirmed by the autocorrelation functions calculated from MD for local motions of the probes (Fig. 6b) where the local order of the magnetic axes of both Q and C* is lower than that in the corresponding duplex forms. As in the case of the duplex DNA the autocorrelation function of the local motion is well represented by three motional modes (Table 3). For both C* and Q-labelled singlestrand fragments the global diffusion remains highly axial with t > 4 3.8t 8 but with both correlation times smaller compared to the duplex form, as expected, (Table 3).
Interestingly, the effective correlation times for the local motions of the probes were somewhat slower than those observed in the duplex forms. At the same time the order parameters for both spin labels attached to the single-stand DNA fragments are noticeably reduced compared to the cases of labelled duplexes. Both differences can be explained by the emergence of additional modes of motion in both labels with larger amplitudes but slower correlation times when not fully stacked.
As with the duplex DNA, the autocorrelation functions of the magnetic axes of Q in the single strand DNA in the laboratory frame (Fig. 6a left) are very close to the ones corresponding to the principal axes of the single-strand DNA (Fig. 6c left). This indicates that Q label in the single-strand structure bears a significant motional contribution from the latter thus serving as an adequate reporter of the DNA motions. In contrast, in the case of C*-labelled singe-strand DNA the nitroxide group of the spin probe is no longer trapped within the major groove with the rotation around the alkyne tether becoming relatively unrestricted. This is confirmed by the calculated time evolution of the dihedral angle between the plane of the cytosine base and the nitroxide ring of C* showing an increased frequency of flips between the 01 and 1801 angles in the case of a single strand (Fig. S7 of ESI †). The greater mobility of C* label results in a more noticeable difference observed between the correlation functions of the magnetic axes of C* (Fig. 6a right) and the principal axes of the single-strand DNA fragment (Fig. 6c right). EPR spectra predicted directly and completely from relevant MD trajectories of single strand DNA are shown as red lines in (a), (b) and (c) of Fig. 7 for the Q labelled fragment at 293 K, the C* labelled fragment at 293 K and the C* labelled fragment at 273 K, respectively. In order to compare the predicted EPR spectrum with the experimental one available from the literature 20 we have performed an additional MD simulation at 273 K. The result is presented in Fig. 7c demonstrating very good agreement with experiment thus confirming the accuracy of the force field employed in this study.
As one would expect, the line shapes for both spin labels are much narrower compared to their counterparts in the duplex form. Two factors contribute to the narrowing of the line shapes in both cases, namely, the decreased correlation times of the global DNA tumbling by a factor of B1.5 and the reduction of the order parameters of the attached probes.
Simulations of EPR spectra using the M-F approach assuming axially symmetric rotational diffusion of the DNA fragments combined with partial averaging of magnetic axes of the probes are shown by a blue line for each of the cases (a), (b) and (c). For single-strand DNA labelled with Q the simulated spectrum is in perfect agreement with the one predicted directly from MD (Fig. 7a) thus confirming the validity of the F-M approach in this case. The same conclusion can be drawn for the case of single-stranded DNA labelled with C* at 273 K (Fig. 7c). The situation is noticeably different in the case of the single-strand C* labelled DNA fragment at 293 K where the simulation by the M-F approach appears to result in much narrower EPR features compared to the ones by the direct method (Fig. 7b). This can be explained as follows. Although the local dynamics of the probes in all three cases (a), (b) and (c) does not strictly satisfy the fast motional regime due to the presence of a slow motional mode (t 3 B 4 ns) (see Table S13 of ESI †) in the case (a) because of the high motional restriction of Q (S = 0.84) the impact of such a slow mode on the spectrum becomes negligible. In the case (c) the global tumbling of the DNA fragment is much slower (t B 9 ns) compared to the label's local motion thus minimising the impact of its slowest local mode on the EPR line shape. The situation is different for the case (b) where the correlation time of the local slow mode (t 3 B 4 ns), as well as the total local effective correlation time (t B 1 ns), become comparable to the correlation time of global motion (t B 3.75 ns) making the impacts from local and global motional both equally significant for the EPR line shape. As a result, the partial averaging of A and g magnetic tensors, used in the EPR simulation procedure, clearly underestimates the effect of local motions of the C* label on the EPR spectrum.
Finally, it is instructive to inspect the sensitivity of the predicted from MD EPR spectra to the choice of force fields employed in DNA modelling. For that reason we have performed an additional simulation of the DNA duplex with the Q spin label at 293 K with the parm99 force field. 69 Parm99, combined with the TIP3P water model, has been used previously to study DNA conformations as well as DNA interactions with proteins. 70,71 The motion and order parameters extracted from the MD simulation are presented in Tables S15 and S16 (ESI †). Fig. S20 of ESI † shows that the predicted from MD EPR spectrum is significantly narrower compared to both the experimental spectrum and the one simulated using the force field parameters employed in this work. Narrowing of the hyperfine coupling lines in the EPR spectrum is attributed mainly to the use of the TIP3P water model which noticeably underestimates the viscosity of water molecules 72 and consequently overestimates the rotational diffusion of DNA. In contrast, the SPC/E water model is known for a reasonably adequate representation of water diffusion. 72 This was recently confirmed in the MD-EPR combined study of different micellar aggregates in water, 34 demonstrating the high feasibility of using our MD-EPR simulation methodology as a test bed for MD force field models. We have therefore also performed an MD simulation using parm99 combined with the SPC/E water model on single-strand DNA labelled with C* at 273 K. Single strand DNA structures have wider amplitudes of local motions that are expected to be more sensitive to DNA force field parameters describing base-base and base-probe interactions. Motional and order parameters extracted from MD run are presented in Tables S17 and S18 (ESI †). Fig. S21 of the ESI, † compares the predictions from the MD EPR spectrum to both the experimental one and the one simulated using the parmbsc1 force field. The results confirm that parmbsc1 provides a better agreement with the EPR experiment compared to parm99.

Conclusions
This study reports the first simulation of the CW EPR spectra of spin labelled DNA fragments in both duplex and single-strand forms from fully atomistic MD simulations. Force field models were developed for two structurally different spin probes, namely, Q and C*, that were the first ones introduced to the studies of DNA structures by EPR spectroscopy. Firstly, EPR spectra predicted directly and completely from the resulting MD trajectories using the direct propagation method demonstrate good agreement with the experimental spectra thus confirming the accuracy of the recently improved version of the parmbsc1 force field for DNA MD modelling. Secondly, structural analysis concludes that the effect of the modified cytosine spin label C* on the helical geometry of the duplex DNA is minimal with the nitroxide group partially restrained within the major groove of duplex DNA. As expected, the data obtained at the fully atomistic level confirm the higher mobility of C* compared to the quinolonyl derived probe Q with the latter shown to be an accurate reporter on the global DNA tumbling. Thirdly, our combined MD-EPR methodology allowed us to test the validity of the application of the M-F approach combined with the partial averaging of the local motions of the probe in the simulation and interpretation of EPR spectra. The 14-mer spin labelled DNA fragments with both Fig. 7 (a and b) Comparison between EPR spectra at 293 K simulated directly and completely from MD (red) and indirectly using the M-F approach with the parameters extracted from MD (blue) of single strand DNA fragments labelled with Q and C*, respectively; (c) comparison among EPR spectra of C*-labelled single strand DNA at 273 K simulated directly and completely from MD (red), using the M-F approach with the parameters extracted from MD (blue) and the experimental one (solid black). The experimental EPR spectrum is reproduced from ref. 20 with permission from Taylor and Francis.
local and global motional contributions having a prominent effect on the EPR line shapes serve as an ideal test bed for the M-F approach. Our results conclude that the M-F approach coupled with partial averaging of magnetic tensors provides an adequate simulation of the EPR spectra when the local motions fall within the fast motional regime and/or are highly restrained. The methods, however, become inadequate when the correlation times of the local and global motional contributions become comparable. It is important to note that parameters employed in the simulation of EPR line shapes by the M-F approach can be readily generated from MD trajectories. The calculation of the relevant autocorrelation functions usually requires much shorter trajectory lengths (o100 ns for DNA fragments in this work) compared to the ones required in the direct propagation method. Thus, in many cases the M-F based EPR simulation methodology using relatively short MD trajectories would be advantageous for the prediction and analysis of EPR spectra compared to direct propagation techniques. This would be crucial in cases of large higher order DNA structures where long scale MD simulations would be challenging or impractical at the all-atom level. Recent advances in click chemistry have led to the design of flexible, base-independent methods for the spin labelling of nucleic acids, 12,73-75 making the EPR studies of DNA increasingly attractive. Several nitroxide based probes representing an improvement over Q and C* in terms of DNA labelling have been recently reported including base independent labelling for both DNA and RNA. [73][74][75] The MD-EPR simulation methodology reported in this work is transferable to the novel DNA labels thus broadening the potential of EPR applications to study the assembly and conformational changes of higher order DNA structures such as the four way Holliday junction. 76 For instance, as have been recently discovered, the assembly of such structures can be induced by small molecules with potential for medical and nanotechnology applications. 77,78 Such MD-EPR applications are currently in progress.

Conflicts of interest
There are no conflicts to declare.