Molecular Simulations Matching Denaturation Experiments for N6-Methyladenosine

Post-transcriptional modifications are crucial for RNA function and can affect its structure and dynamics. Force-field-based classical molecular dynamics simulations are a fundamental tool to characterize biomolecular dynamics, and their application to RNA is flourishing. Here, we show that the set of force-field parameters for N6-methyladenosine (m6A) developed for the commonly used AMBER force field does not reproduce duplex denaturation experiments and, specifically, cannot be used to describe both paired and unpaired states. Then, we use reweighting techniques to derive new parameters matching available experimental data. The resulting force field can be used to properly describe paired and unpaired m6A in both syn and anti conformation, which thus opens the way to the use of molecular simulations to investigate the effects of N6 methylations on RNA structural dynamics.


Introduction
Post-transcriptionally modified nucleotides are crucial for RNA function [1,2].Methylation of adenine in the N6 position (m 6 A) is the most prevalent chemical modification in messenger RNAs and has been observed both in coding and noncoding RNAs [1,2,3,4].m 6 A can fine regulate the interaction of RNA with specific proteins known as m 6 A readers.In addition, similarly to other chemical modifications [5], it can directly affect RNA stability and structural dynamics (see, e.g., [6,7,8]).Specifically, m 6 A has been suggested to weaken Watson-Crick pairings due to the incompatibility of its most stable conformation (syn) with duplex formation [9,10,11] (see also Fig. 1).Interestingly, recent nuclear magnetic resonance (NMR) experiments have identified syn/anti dynamics in both paired and unpaired m 6 A, recapitulating the effect of N6 methylation on RNA conformational kinetics [11].Molecular dynamics (MD) simulations give access to structural dynamics at the atomistic resolution [12] and are thus an ideal tool to complement NMR studies.The capability of classical force fields to predict the dynamics of difficult structural motifs is steadily increasing [12,13,14].However, the number of applications of MD simulations to N6-methylated RNAs reported to date is still limited [11,15,16,17,18,19,20].Force fields parameters for m 6 A developed by Aduri et al. [21] were determined in a bottomup fashion and are compatible with the AMBER force field, which is widely adopted for RNA simulations [12].These parameters are the default choice for a structure refinement tool [22].However, they have been quantitatively validated against a limited set of experimental results [18].The recent publication of denaturation experiments on a number m 6 A-containing duplexes [9,23] calls for a more extensive validation of force-field parameters and, ideally, for fitting force-field parameters directly on experiments [24].
In this paper, we validate and improve the parameters introduced in Ref. [21] by using alchemical free-energy calculations (AFECs) [25].To this end, an unmodified adenine is converted to a modified one by switching on/off nonbonded interactions of specifically chosen atoms, in both single-stranded and double-stranded RNAs [18] (see also Fig. 2).We then develop a reweighting technique that can be used to predict results corresponding to a different set of charges without the need to perform new MD simulations.Additionally, we extend a recently-introduced force-field fitting strategy [26] to be usable in the context of alchemical simulations.The introduced approach allows training six charges and a dihedral potential so as to quantitatively reproduce methylation effects in denaturation experiments.The resulting force field can be used to properly describe paired and unpaired m 6 A in both syn and anti conformation.

Simulated systems
We simulated the isolated m 6 A nucleoside, 9 m 6 A-methylated duplexes for which denaturation experiments are available in literature [9,23] (see Table 1), and the corresponding single stranded RNAs.For the isolated m 6 A nucleoside, we computed the ∆G syn/anti by taking the difference in the ∆Gs obtained with AFEC by methylating the adenosine in syn or anti conformations.We then chose a number of systems from Ref. [9].For systems A4 and A5, where m 6 A is present as dangling end and thus unpaired, we only performed AFEC corresponding to the syn conformation.For the other systems, we performed AFEC in the expected anti conformation.For the A2 and A3 systems we additionally performed AFEC in the unexpected syn conformation as a validation (population reported in Ref. [11] is ≈ 1%).In addition, we chose 5 more systems from Ref. [23], with the following criterion: they have a single methylation per-strand and the methylation occurs in an internal position of the duplexes.For all these systems, we performed AFEC in the expected anti conformation.
Starting structures for MD simulations were built using the proto-Nucleic Acid Builder [27].Single strands were generated by deleting one of the chains from duplex structures.All the MD simulations were performed using a modified version of GROMACS 2020.3 [28] which also implements the stochastic cell Figure 1: N 6 -methyladenosine (m 6 A) nucleobase in anti (less stable) and syn (more stable) conformations (panel a) [9,11].Atom names in red correspond to charges reparametrized in this work.Example of a duplex containing m 6 A in anti conformation, which is the expected conformation for the nucleotide when Watson-Crick paired (panel b).A 6 is used to denote m 6 A in secondary structures for compactness.Example of a duplex with m 6 A as a dangling end in syn conformation (panel c).The m 6 A methyl group is highlighted in yellow.
Figure 2: Thermodynamic cycle.Alchemical free-energy calculations (AFEC) allow computing ∆G by integrating along an alchemical path λ describing the transformation of non methylated adenosine into m 6 A, by switching on/off nonbonded interaction of specifically chosen atoms.The relative free-energy change due to the modification can be estimated as the ∆∆G between AFECs performed on a duplex and on the corresponding single strand.This quantity can be directly compared to the difference in thermodynamic stability of duplexes with or without the modification, that can be measured experimentally through denaturation experiments.
Table 1: List of systems involved in the fitting and relative experimental ∆∆G.The first system is the single nucleotide and the experimental value corresponds to ∆G syn/anti .In the A2-A5 and B1-B5, the "6" in the double strand sequences is used to identify m 6 A for compactness.∆∆Gs for systems A2-A5 were measured by Roost et al. [9].In A4 and A5, the m 6 A is positioned as a dangling end and has a stabilizing effect on the duplex.Experiments for systems B1-B5 were performed by Kierzek et al. [23].In B2-B5 systems, the methylation occurs in both strands, however the ∆∆Gs reported are to be intended "per-methylation" .
rescaling barostat [29].The AMBER force field was used for RNA [30,31,32], TIP3 model for water [33], and Joung and Cheatham parameters for ions [34].As a starting parametrization for m 6 A, we used AMBER adenosine parameters combined with modrna08 [21] charges for the nucleobase, adjusted so as to preserve the total charge of the nucleoside.Details on the implementation of these parameters and on the initial tests are reported in Section S1.Charges are given in Table S2.We refer to this parametrization as Aduri force field.

Alchemical Simulations
For AFECs, we included a hybrid adenosine with double topology in the forcefield definition: the first topology corresponding to standard adenosine, and the second one corresponding to m 6 A. We used 16 replicas in which Lennard-Jones parameters and partial charges were simultaneously interpolated.In order to avoid singularities due to electrostatic interaction when the repulsive LJ potential is switched off [25], we used soft-core potentials as implemented in GRO-MACS (sc-alpha=0.5,sc-sigma=0.3, and sc-power=1) [35].Simulation boxes consist in rhombic dodecahedrons containing RNA, water, Na + and Cl − ions with an excess salt concentration of 0.1 M. For a subset of the systems, further simulations were performed for a salt concentration of 1 M.The systems were energy minimized and subjected a multi-step equilibration procedure for each replica: 100 ps of thermalization to 300 K in the NVT ensemble was conducted through the stochastic dynamics integrator (i.e., Langevin dynamics) [36], and other 100 ps were run in the NPT ensemble simulations using the Parrinello-Rahman barostat [37].In production runs, the stochastic dynamics integrator was used in combination with the stochastic cell rescaling barostat [29] to keep pressure at 1 bar.Equations of motion were integrated with a time-step of 2 fs.Long-range electrostatic interactions were handled by particle-mesh Ewald [38].During production, Hamiltonian replica exchange was used proposing exchanges every 200 fs.The set of λ values defining replica's Hamiltonians were chosen in such a way to guarantee transition probabilities above 20% and as homogeneous as possible (see Section S2), ensuring a sufficient phase space overlap between replicas.Each replica was simulated for 10 ns, for a total of 16 × 10 ns = 160 ns for each system.To decrease numerical errors in energy recalculations, trajectories were saved in uncompressed format.At the end of the production phase, the 16 independent "demuxed" (i.e., continuous) trajectories were processed to recompute energies for each of the 16 Hamiltonian functions so as to compute ∆G via binless weighted histogram analysis method (WHAM) [39,40,41].Specifically, for each trajectory a weight w was found for each snapshot x that allows computing statistics for the unmodified adenine as a weighted average over the set of concatenated replicas.These weights were then used to compute the ∆G associated to the methylation as: where ∆E(x) = E λ=1 (x) − E λ=0 (x) is the difference between the total energy computed with the Hamiltonian energy functions associated to m 6 A and adenosine, respectively.We used a bootstrapping procedure to compute the statistical error on ∆G estimates by resampling the 16 continuous trajectories 200 times with replacement [42].As a control, we computed ∆Gs using the standard Bennett-acceptance-rate estimate implemented in Gromacs [43,44].∆∆Gs were obtained taking the difference between ∆Gs obtained by methylating the adenosine in anti or syn conformation on the duplex or dangling end, respectively, and the ∆G obtained methylating in syn conformation on the relative single strand.Transitions between syn and anti states were never detected during the alchemical simulations.In this way, the contribution to the free energy given by the syn (anti ) conformation in the duplex (single strand or dangling end) were ignored.Indeed, we expect these contributions to be negligible based on the experimental evidences [9,11], which show a syn/anti isomer preference when paired (≈1:100) versus unpaired (≈10:1).This was additionally verified with supplementary simulations performed on the A2 and A3 systems.Moreover, we computed ∆G syn/anti by performing the alchemical transformations on the isolated nucleoside in solution for the two isomers and computing their difference.

Fitting Procedure
We employ a fitting strategy based on reweighting [26] where a subset of the partial charges and a dihedral potential are adjusted to match experimental data.Specifically, we decided to fit the atoms that are closer to the methyl group (N6, C6, H61, C10, H101/2/3, and N1, see Fig. 1).The total charge was maintained, leading to 5 free parameters associated to the partial charges.A single cosine was added to the η 6 torsional angle identified by atoms N1-C6-N6-C10: . This angle controls the syn/anti relative populations, leading to a total of 6 parameters, and the shift is chosen so that a positive value of V η favors syn configurations over anti.
To optimize the calculation of the total energy of the system at every iteration of our fitting procedure, where up to 6 charges were possibly modified, we notice that the total energy of the system is a quadratic function of the charge perturbations ∆Q i .Without loss of generality, one can write the energy change associated to charges and torsion perturbation as In total, for every analyzed snapshot (x), 20 coefficients (K i and K ij ) can be precomputed that allows obtaining the energy change for arbitrary choices of ∆Q with simple linear algebra operations, without the need to recompute electrostatic interactions explicitly.The coefficients were obtained by using GRO-MACS in rerun mode for 20 sets of test charge perturbation that were extracted from a Gaussian with zero average and variance set to 1 e.The perturbations were constructed to maintain the total charge constant.Importantly, this approach correctly takes into account the effect of charge perturbations on 1-4 interactions, where electrostatics is scaled with a force-field-dependent fudge factor, as well as on 1-2 and 1-3 interactions, for which it is discarded, and interaction with all the periodic images.The second order expansion above is exact if one neglects roundoff errors.The magnitude of charge perturbations was chosen to minimize such errors.Eq. 1 should then be suitably modified replacing ∆E with ∆E + ∆U .Its derivatives with respect to the free parameters (charge and dihedral potential coefficient) can be computed as well (see Section S6).
Our fitting is based on the minimization of a L2-regularized cost function defined as follows: Here the χ 2 measures the discrepancy between computations and experiments, whereas the regularization terms on the charges and on the torsional η 6 are governed by the hyperparameters α and β and are needed in order to avoid overfitting on the training set.This function is minimized using the L-BFGS-B method [45] as implemented in SciPy [46].
The result crucially depends on the choice of the hyperparameters α and β.Lower values for the hyperparameters imply that larger corrections are allowed, with the risk of overfitting, and thus lower transferability to new experiments.Higher values for the hyperparameters imply that lower corrections are allowed, with the risk of underfitting, and thus lower accuracy in reproducing experimental data.The sweet point could be in principle found with a cross validation (CV) procedure and a scan over possible values for α and β [26,24].For the smallest dataset (set A in Table 1), we used a leave-one-out CV strategy, i.e., we trained the parameters on all systems except one.For the largest dataset (set AB in Table 1), we used a leave-3-out strategy, iteratively training the parameters on 7 randomly chosen experiments and validating on the 3 left out experiments.In both cases, we then assessed the transferability of the model by evaluating its average χ 2 on the system (or the subset of systems) that was left out.

Statistical Significance
When recomputing energies through a reweighting procedure, particular attention must be taken towards the statistical significance that may be lost during the computation, by reducing the effective sample size of the data set.This is usually monitored computing the Kish effective sample size [47,48].In our case, the most affected ensemble is the one corresponding to m 6 A (λ = 1).We thus monitor the Kish size computed using weights corresponding to the λ = 1 ensemble, defined as We then compare it with the Kish size obtained with the original force field, defined as To quantify how much statistical efficiency is lost due to the reweighting to a modified set of parameters we use the Kish size ratio (KSR), that we define as 3 Results In this work, we fit point charges and a single torsional potential correction for a m 6 A RNA residue using alchemical MD simulations and a set of experimental data.In all the fittings, charges and torsional potential were subject to L2 regularization with hyperparameters α and β, respectively.We initially employed only the first 5 experimental datapoints of Table 1, namely (A1) ∆G syn/anti for a nucleobase and (A2-A4) ∆∆G in melting experiments [9].Thus, we first report the results obtained with such set of charges, including a validation done on a more recent set of melting experiments (B1-B5) [23].We then report results obtained with charges that were fitted on the entire dataset (A1-A5 and B1-B5).As a reference, results obtained with the Aduri et al. [21] modifications (modrna08) for the commonly used AMBER force field are also reported, either as is or complemented with a custom torsional correction that result in a ∆G syn/anti matching experiment A1.All the calculated ∆Gs are reported in Table S3.A complete list of the performed alchemical simulations is reported in Section S9.

Fitting on the smaller dataset
For this first fitting, we only employed data set A1-A5 from Table 1.χ 2 errors were computed assuming the experimental error of each data point to be equal to each other and to 1 kJ/mol.Figure 3a reports the results of a cross-validation test performed with a leave-one-out procedure.Namely, we fit the whole experimental dataset leaving out one experimental data point at a time and report the average error on the left-out experiment.In this leave-one-out procedure we decided not to iterate on the ∆G syn/anti experiment (A1), since this is expected to be crucial to correctly reproduce the conformation of non-Watson-Crick-paired residues (mostly syn).From this map, we can hardly appreciate any variation of the χ 2 along the vertical axis corresponding to the β hyperparameter.This suggests that β could be set to zero, thus simplifying all subsequent hyperparameter scans.Conversely, the χ 2 grows significantly for low α values.This implies that regularization on charges is required to avoid overfitting.In general, one should expect a minimum to be observed in this type of hyperparameter scans [26,24].This is not the case here for the α scan (see also Fig. S3, showing projection on α for β = 0), implying that the performance of the parameters on a given system are not improved when excluding that system from the training set.This is likely due to the small data set employed.
Figure 3b shows the optimized parameters (charge and torsional corrections) as a function of the regularization hyperparameter α while fixing β = 0.A transition can be seen at α ≈ 10.Namely, when α > 10, parameters have a smooth dependence on α, whereas when α < 10, both the charges and the torsional potential change suddenly.In the limit α → ∞, it can be seen that charge corrections tend to zero with an inverse law dependence, which is expected for L2 regularization, and the torsional correction tends to V η ≈ 1.5 kJ/mol, which corresponds to the amplitude of the torsional potential that optimizes the χ 2 without modifying the charges of the reference Aduri et al. model.We notice that ∆G syn/anti obtained when using the Aduri et al. force field is ≈ 1.7 kJ/mol, and thus this correction results in ∆G syn/anti ≈ 1.7 + 2 × 1.5 = 4.7 kJ/mol, which is still smaller than the experimental reference ≈ 6 kJ/mol.The obtained parameters indeed strike a balance between favoring the syn state in the isolated nucleoside and not favoring it too much in the single-stranded calculations used to predict the ∆∆G from melting experiments, which would lead to too large destabilizations associated to the methylation.When α is decreased, the optimal torsional correction changes, since all the parameters are coupled.This confirms that charges and torsional parameters should be fitted simultaneously.
Figure 3c shows the individual χ 2 associated to the same hyperparameter scan.The average χ 2 error is, by construction, monotonically increasing with α, and most of the individual errors follow the same trend.Figure 3c also shows the statistical efficiency of the analysis, quantified by the relative reduction of the Kish effective sample size associated to reweighting.A low number here indicates that the tested charges are so different from those employed in the simulation to make the result statistically not significant.The Kish size displays a significant drop for α < 10, indicating that results in this regime might be not significant.This is a likely explanation for the discontinuous behavior observed in Fig. 3b.
We then tested the charges obtained with this reduced training set on the newer data set B1-B5, see Table 1, which was not included in the training phase.This set of data involves 5 recently published melting experiments [23], 4 of which have m 6 A occurring in both chains of the duplex.We notice that double methylations are expected to lead to an even lower statistical efficiency of the reweighting procedure.We thus performed this analysis by reweight-ing simulations that were generated using an alternative set of parameters that was derived with a similar fitting procedure shown in this paper, but with an incorrect regularization on the charges as described in Section S4.Since this parametrization is closer to the right solution of the fitting with respect to Aduri, it obtains higher Kish size values in the relevant α range (see Fig. 3d).The χ 2 computed on the second data set shows that an optimal result can be obtained setting α ≈ 10.We also compared with results obtained using the original Aduri charges and optionally including a torsional correction to fix the syn/anti balance.These results are obtained with direct simulation, that is without reweighting.It can be seen that the results with the parameters trained on systems A1-A5 largely outperform those obtained with Aduri parameters on systems B1-B5, thus confirming the transferability of the parameters.Aduri+tors parametrization corresponds to setting V η =2.35 kJ/mol in such a way to perfectly fit experiment A1 (single nucleoside) without modifying charges.χ 2 computed for Aduri+tors demostrates that acting exclusively on the torsional is not sufficient to reproduce both ∆G syn/anti and melting experiments.It is also important to note that the improvement in reproducing experiments is obtained by changes in the partial charges that are small when compared to differences between charges derived with the standard restrained electrostatic potential protocol [49] in different conformations (see Section S5).

Fitting on the full data set
Next, we perform a fitting using the full data set reported in Table 1.Since the variability of error in this data set is larger, we here computed χ 2 using the experimental errors reported in Table 1.For the ∆G syn/anti experiment, for which an experimental error is not reported, we used a nominal σ = 0.5 kJ/mol so as to assign to this experiment a larger weight when compared to the other data points corresponding to melting experiments.
Figure 4a reports the results of a cross-validation test performed with a leave-three-out procedure.Namely, we randomly select seven systems to be used in training and we report the average χ 2 error obtained for the remaining three systems.This time also system A1 was allowed to be left out from the training set.Results are qualitatively consistent with those obtained with the smaller data set (see Fig. 3).It is difficult to appreciate any variation of the χ 2 along the vertical axis corresponding to the β hyperparameter, suggesting that we can safely set β = 0. We also do not find any clear minimum when scanning over α (see also Fig. S4, showing projection on α for β = 0).Figure 4b shows the parameters as a function of the regularization hyperparameter α while fixing β.A clear transition can be seen at α ≈ 20.The average χ 2 error is monotonically increasing with α, but some of the systems have a non-trivial behavior (Fig. 4c).The Kish size shows a significant drop for α < 50, showing that results in this regime might be not statistically reliable.We thus select the parameters obtained with α = 50 as the optimal ones trained on the entire data set.
We then compare the performance of a number of different sets of param-  Table 2: Charge modifications (∆Qs) and torsional potential (V η ) for the fitting performed on the smaller dataset (fit A, α = 10) and for the fitting performed on the larger dataset (fit AB, α = 50).
eters in reproducing all the available experimental data points.Namely we compare (a) the original Aduri parameters (Aduri), (b) the Aduri parameters augmented with a torsional correction to enforce the correct syn/anti balance in a nucleobase (Aduri+tors), (c) the parameters obtained fitting on the initial dataset (A1-A5), with hyperparameter α = 10 (fit A), and (d) the parameters obtained fitting on the full dataset (A1-A5 and B1-B5), with hyperparameter α = 50 (fit AB).Results are reported in Fig. 4d.The quality of the fit is also summarized in the reported χ 2 values.The addition of a simple torsional correction to the Aduri parameters results in an overall worsening of the agreement with experiment.On the other hand, the two sets of parameters obtained in this work (fit A and fit AB) display a significantly better agreement with experimental data.Points are computed directly from energies generated during the MD, except for fit A case in which energies are computed through a reweighting procedure that gives a KSR about 0.65.Note that even if fit A has the best performance, fit AB allows to get very similar χ 2 with a stronger regularization (higher α), keeping the parametrization closer to the reference one, as shown in Figs.S6-S7.Finally, as observed in the previous subsection, we note that the improvement in reproducing experiments is obtained by relatively small changes in the partial charges (see also Section S5).The fitted parameters are summarized in Table 2.

Relative stability of syn and anti conformations
One piece of the experimental information that we implicitly used in our fitting procedure is the relative stability of syn and anti conformations in a nucleoside (data point A1, in Table 1).We indeed assumed a predominant population of syn conformation for the unpaired nucleotides used in the reference single-stranded systems.We also assumed that m 6 A adopts exclusively its anti conformation when paired, in agreement with experiments [9,11].In particular, Ref. [11] reports that, for the most common G6C sequence, m 6 A forms a Watson-Crick base pair with uridine that transiently exchanges on the millisecond time-scale between a main substate (anti ) and a low populated (1%), singly hydrogen-bonded and mismatch-like conformation through isomerization of the methylamino group to the syn conformation.This population corresponds to a ∆G duplex syn/anti ≈ −11 kJ/mol.We a posteriori validated this population by performing alchemical transformations on the duplex systems enforcing the syn  3: Free-energy differences between syn and anti isomer states in systems A1-A3.The last column corresponds to experimental estimates, whereas other columns corresponds to computed ∆∆G for different parametrization.Energies are given in kJ/mol units.
conformation.The predicted ∆G syn/anti for a nucleotide and for two of the tested duplexes are reported in Table 3, where the corresponding experimental values are also included.For the A1 experiment, as expected, the proposed sets of parameters closely match the experimental value, that was used during training.The Aduri et al. force field underestimates the ∆G syn/anti , resulting in a relatively high population of the unexpected anti conformation in a nucleoside.This difference can be directly corrected with a torsional potential applied on the η torsion (Aduri+tors).However, when analyzing duplexes A2 and A3 with the Aduri+tors parameters, we found that the predicted ∆G syn/anti would be close to zero, in fact making the assumption of neglecting the syn conformation in duplexes in our alchemical calculations difficult to justify, and in disagreement with experimental findings.In other words, the original Aduri charges allow to reproduce the relative stability of syn and anti conformations either in the paired state (Aduri parameters) or the unpaired state (with torsional correction), but not in both simultaneously.Remarkably, the sets of parameters proposed here, that also contain a torsional term penalizing the anti conformation, result in a significantly higher value for ∆G dup syn/anti , much closer to a qualitative agreement with experiment.This suggests that the proposed parameters better describe the interactions of the m 6 A nucleobase with the surrounding environment and are thus more transferable.We notice that the relative stability of syn and anti conformations is predicted to be sequence dependent, being different for system A3 (sequence U6G).
To gain insight about how the m 6 A-U pairings occurs in the duplexes, we analyzed snapshots of system A2, both for m 6 A in syn and anti, together with histograms of distances between atoms belonging to the two nucleobases (Fig. 5).The reported histograms are unimodal and with an increased average associated to the distortion of the A-U Watson-Crick pairings due to the steric clash induced by the methylation.However, the hydrogen bond between A-N1 and U-H3 is present, in agreement with what has been suggested previously [11].

Interpretation of parameters and dependence on ionic strength
To provide an interpretation for the obtained parameters, we performed a few additional fittings.In particular we investigated which charges have a major Figure 5: Interfacing atom distances for m 6 A-U pairing in system A2 in Table 1, for anti conformation (left) and syn (right).Histograms show unimodal distributions, and the averaged value are indicated in the box.Distances are sampled from the alchemical trajectories considering only the λ = 1 replica.
In the syn conformation, m 6 A-H10 correspond to the hydrogen of the methyl group closest to the uracil oxygen O4.
impact on enforcing agreement with experiments.To this purpose we performed fitting on exclusively two charges at a time plus the torsional term, considering couples of charges that in fit A and fit AB have systematically positive and negative ∆Qs.Results are discussed in Section S7 and Fig. S9.Overall, the results suggest that the main contribution of the fitted correction is to increase the stability of Watson-Crick hydrogen bonds by making N1 and H61 more polar and at the same time using the η torsional potential to control the syn/anti relative population.As a further test, we simulated a subset of the systems using a salt concentration of 1 M, which is consistent with that used in experiments.As shown in Section S8, results are equivalent to those obtained at 0.1 M salt concentration.We interpret this result with the fact that the methylation is not altering sufficiently the electrostatic environment to be sensitive to changes in ion concentration.This implies that a training performed using simulations performed at a different ion concentration would result in an equivalent set of parameters and further confirms the robustness of our results.

Discussion
In this work, we proposed a protocol to parametrize charges in modified nucleobases using available melting experiments.The approach is applied to m 6 A and leads to a set of charges that can reproduce a set of 10 independent experimental values.The approach is based on the force-field fitting strategies introduced in earlier works [50,51,26], which are here extended with a number of technical improvements.
A first methodological contribution is a formalism that allows alchemical calculations to be used as a reference.Previous works were only using observables computed with a single set of force-field parameters [50,51,26,14].The method introduced here allows free-energy differences between different sets of parameters to be evaluated and compared with experiment.This opens the way to the optimization of parameters based on experimentally measured ∆∆Gs.We based our analysis on optical melting experiments, which are commonly employed in the nucleic acids community [52], but other types of experiments might be considered.In our specific application, only the parameters of one of the two end states were refined, but one could similarly fit parameters for both adenosine and m 6 A, at the price of increasing the number of parameters and thus the risk of overfitting.A second improvement is that we developed a way to efficiently recompute the total energy of the system using test charges.This is achieved by precomputing the total electrostatic energy of the system with a set of randomly perturbed charges.Given the high cost of electrostatic calculations, this makes the cost of each of the iterations performed during force-field fitting significantly faster, and implicitly takes into account combination rules, non-bonded exclusions, and periodicity.These two improvements can be readily integrated in other MD-based force-field optimization strategies.
A limitation of optimizing charges with the introduced procedure is that the statistical efficiency of reweighting is significantly decreased even by small charge perturbations.This implies that simultaneously parametrizing many copies of the same nucleotide, or parametrizing a larger number of charges for the same nucleotide, would be more difficult.In our case, we had to include at most two m 6 A residues in the same simulation.If more copies of the same reparametrized nucleotide are needed in the same system, one might have to design strategies where only a few copies at a time are reparametrized, or follow an iterative procedure where modifications are included in consecutive steps [26].In this application, this was not necessary.
Overfitting was avoided by using a standard L2 regularization term on the charge increments.This penalty does not depend the charge location.Importantly, the regularization hyperparameters tune the relative weight of the experimental data and of the reference charges, here taken from Ref. [21], thus allowing to achieve a meaningful set of parameters also in regimes where the number of data points is very limited.More effective regularization strategies might be designed based on the molecular dipole, as done in Ref. [53], so as to minimize the perturbation of the electrostatic potential at a large distance from the molecule.Alternatively, one might directly use as a regularization term the deviation from the quantum-mechanical electrostatic potential at short distance.In the limit of a large regularization hyperparameter, this would lead to ESP charges [54].Finally, other regularization criteria might be used [14].When comparing our procedure with standard ESP charge fitting, it is important to realize that we are aiming to reproduce experimentally observed ∆∆G, which are non-linear functions of the energy of each configuration, which in turn is a quadratic function of the charges.These non-linearities make it possible for multiple local minima of the cost function to exist, and could thus make the minimization not reproducible.However, when sufficiently regularized, the fitting procedure results in reproducible charges that depend smoothly on the control parameters.In standard ESP fitting, instead, the electrostatic potential is fitted, thus resulting in a linear fit with a unique solution.
We notice that the parameters of the unmethylated force field were not modified.This was based on the assumption that the employed set of force-field parameters is already capable to reproduce ∆∆G experiments associated to mutations between non-modified nucleobases [55].The m 6 A charge optimization could be easily repeated using another set of initial parameters, and the parameters of non-modified nucleobases might be adjusted as well, although with the caveat discussed above.
Another possible limitation of the employed alchemical simulations is the sufficient sampling of the end states.The duplex is expected to be stable and well structured, so that sampling multiple structures should not be necessary.For selected cases, we also explored the possibility to include the unlikely syn paired state, which, as expected, gives a negligible contribution to the stability of the duplex.For single strands, instead, we only sampled the syn state.More importantly, our simulations were short enough to avoid any significant reconformation of the single strand.Sampling the conformations of flexible, single stranded RNAs is notoriously difficult [12].In addition, the generated ensemble might contain artificially stabilized intercalated structures, whose population is known to be overestimated by the RNA and water force fields adopted here [56,57].This would make the correct sampling of the single stranded state unfeasible.We also notice that the experimental results that we aimed to reproduce were performed on systems designed to have the isolated strands unstructured, so as to capture the effect of methylation on hybridization.Putting everything together, we conclude that the approximation of a single strand ensemble that do not depart too much from the initial A-form helix is a sensible choice for this specific application.
An important finding of this work is that the parameters of Aduri et al. cannot reproduce the syn/anti balance expected for m 6 A residues.This balance is extremely important, and is related to the mechanism by which m 6 A modifications modulate duplex stability [9].This could not be rectified with a straightforward correction on the single involved torsion.The optimized charges, instead, allow the correct syn/anti balance to be recovered both in paired and unpaired nucleobases, as well as an heterogenous set of optical melting experiments to be reproduced.Interestingly, the Aduri et al. parameters were tested in a recent work [18], with results for system A2 in Table 1 consistent with ours and with experiments.However, systems A1 and A3 were not tested, and thus the problem that we observed here could not be identified.Another interesting finding is that the ∆∆G associated to N6 methylation are here predicted to be independent of ion concentration.We are not aware of any experimental validation of this finding, which could be obtained by comparing melting experiments at different ion concentration.Finally, our results suggest that the relative population of the syn excited state in duplexes [11] might significantly depend on the identity of the neighboring nucleotides.The precise hybridization kinetics could thus be quantitatively different for RNAs with different sequences.
A convenient property of our approach is that it does not require changing the functional form of the interaction potential, so that new parameters can be readily incorporated in existing MD software.This is not the case if ad hoc corrections are employed [58,14].In addition, it is worth noting that the charge modifications obtained are very small, and in particular they are smaller than the typical difference between sets of charges derived with slightly different procedures or using different reference conformations.In spite of this small difference, the effect on experimental observables is significant.These observations imply that there is still significant space to improve the performance of current force fields without necessarily modifying the functional form, if experimental information is used during training [24].
Using our approach it is possible to dissect the individual contribution of the modified force-field parameters.The main factors playing a role on the change of duplex stability induced by m 6 A methylation are (a) the penalty for switching to the unfavored anti isomer [9], (b) the stabilization induced by hydrophobic shielding of the methyl group against surrounding bases (see also Fig. S8) [59,60], (c) the impact of partial charges on stacking interactions [59], and (d) the impact on the strength of Watson-Crick hydrogen bonds.Since, on average, experimental ∆∆G for denaturation experiments performed on duplexes are smaller than the anti isomer penalty, we could expect that the sum of the other factors has a stabilizing effect on the majority of the considered duplexes.We notice that Aduri charges for N1 and H61, which are involved in Watson-Crick pairings with the complementary uridine, have partial charge absolute values significatively lower compared to the standard adenine parameters (0.28948 vs 0.41150 for H61, −0.675968 vs. −0.76150for N1).This may lead to a weakening of hydrogen bonds which may cause an overestimation of destabilization induced on duplexes, as we observed in Aduri+tors cases (see Fig. 4d).The results of our fitting systematically increase the absolute value of H61 and N1 partial charges, hence resulting in a stronger Watson-Crick pairing.At the same time, the torsional term allows to reproduce the correct anti isomer penalty.Parameters are coupled, so that it is necessary to fit them simultaneously so as to avoid double counting effects.
To the best of our knowledge, this is the first attempt to tune partial charges of a biomolecular force field based on experiments performed on macromolecular complexes.We expect that this approach could be used in the future to improve the capability of biomolecular force fields to match experimental observations exploiting a part of the functional form that has been traditionally derived in a bottom-up fashion.Remarkably, the parameters derived here for m 6 A allow to properly describe paired and unpaired m 6 A in both syn and anti conformation, and thus open the way to the use of molecular simulations to quantitatively investigate the effects of N6 methylations on RNA structural dynamics.

Figure 3 :
Figure 3: Results obtained with parameters fitted on the initial dataset, A1-A5 in Table1.Cross validation error obtained with a leave-one-out-procedure, shown as a function of the two regularization hyperparameters α, for charges, and β, for the torsional potential (panel a).Darker green colors correspond to lower values of the average χ 2 computed on the systems left out iteratively from the fitting.Parameters (∆Q and V η ) obtained from the entire initial dataset as a function of α, with β = 0 (panel b).χ 2 errors for individual experiments and Kish size ratio (KSR, see text for definition) obtained using parameters fitted on the entire initial dataset as a function of α, with β = 0 (panel c).Validation on the second dataset (B1-B5 in Table1) of the parameters obtained on the first dataset (panel d).Results using Aduri parameters are shown as horizontal lines, either as reported in the original paper (green) or including a single torsional correction to obtain the correct syn/anti population (data point A1).

Figure 4 :
Figure 4: Results obtained with parameters fitted on the full dataset, A1-A5 and B1-B5 in Table 1.Cross validation error obtained with a leave-three-outprocedure, shown as a function of the two regularization hyperparameters α, for charges, and β, for the torsional potential (panel a).Darker green colors correspond to lower values of the average χ 2 computed on the systems left out iteratively from the fitting.Parameters (∆Q and V η ) obtained from the entire dataset as a function of α, with β = 0 (panel b).χ 2 errors for individual experiments and Kish size ratio (KSR, see text for definition) using parameters fitted on the entire initial dataset as a function of α, with β = 0 (panel c).∆∆G computed for each of the ten analyzed systems with 4 different sets of parameters (panel d).fit A are parameters obtained fitting on the first data set (A1-A5) with regularization α = 10.fit AB are derived fitting on the entire data set (A1-A5 and B1-B5) for α = 50.χ 2 obtained for each force field set of parameters are shown in the table inside panel d.

Table 1
) of the parameters obtained on the first dataset (panel d).Results using Aduri parameters are shown as horizontal lines, either as reported in the original paper (green) or including a single torsional correction to obtain the correct syn/anti population (data point A1).