Predicting antibiotic resistance in complex protein targets using alchemical free energy methods

Abstract Drug resistant Mycobacterium tuberculosis, which mostly results from single nucleotide polymorphisms in antibiotic target genes, poses a major threat to tuberculosis treatment outcomes. Relative binding free energy (RBFE) calculations can rapidly predict the effects of mutations, but this approach has not been tested on large, complex proteins. We use RBFE calculations to predict the effects of M. tuberculosis RNA polymerase and DNA gyrase mutations on rifampicin and moxifloxacin susceptibility respectively. These mutations encompass a range of amino acid substitutions with known effects and include large steric perturbations and charged moieties. We find that moderate numbers (n = 3–15) of short RBFE calculations can predict resistance in cases where the mutation results in a large change in the binding free energy. We show that the method lacks discrimination in cases with either a small change in energy or that involve charged amino acids, and we investigate how these calculation errors may be decreased.


| INTRODUCTION
Tuberculosis is a difficult disease to treat; the standard regimen is four antibiotics, rifampicin, isoniazid, pyrazinamide, and ethambutol, for 6 months. An infection that is resistant to both rifampicin and isoniazid is called multi-drug resistant tuberculosis (MDR-TB) and the treatment regimen recommended by the World Health Organization (WHO) is complex but always includes levofloxacin or moxifloxacin, which are fluoroquinolones. 1 Rifampicin (RIF) acts by binding to the β-subunit of the RNA polymerase (RNAP, encoded by the rpoB gene), preventing the extension of the RNA ( Figure 1A). The most common resistance-conferring mutation is rpoB S450L, however a wide range of mutations have been observed clinically. [2][3][4][5] The majority of these are found in amino acids 428 to 452 which pack against the drug (usually known as the "rifampicin resistance determining region" or RRDR), enabling the development of nucleic acid amplification tests, such as the Cepheid GeneXpert MTB/RIF system which is endorsed by the WHO for diagnosis of MDR-TB. 6,7 Not all non-synonymous mutations in the RRDR, however, confer resistance, for example rpoB L443F. 5 Nor does resistance arise purely within the RRDR: rpoB I491F and V170F are proximal to S450L and the former was suspected to be behind an outbreak of MDR-TB in Eswatini since it is not detected by GeneXpert. 8 The fluoroquinolones target the DNA gyrase (DNAG), a tetrameric enzyme which unwinds DNA by forming and re-ligating double stranded DNA breaks prior to transcription and replication ( Figure 1B). Specifically, two fluoroquinolone molecules intercalate into DNA breaks and bind specific gyrA residues via a coordinated Mg 2+ ion. This stabilizes DNA-DNA gyrase covalent linkages and prevents re-ligation of DNA double stranded breaks.
The most common DNA gyrase mutations found in MDR-TB samples are gyrA D94G and gyrA A90V and these mutations are strongly associated with fluoroquinolone resistance. 2,3,[9][10][11][12] These residues are part of the gyrA "quinolone resistance determining region" (QRDR), defined as gyrA codons 74-113. 13 However, again, not all mutations in this region confer resistance, leading to false positive resistance results in genotypic assays. 14 Rarely seen DNA gyrase mutations in gyrB are also associated with fluoroquinolone resistance, and a gyrB QRDR from residues 461 to 501 has also been proposed. 15 The residues of the two QRDR regions make up the fluoroquinolone binding pocket, and gyrB A642P is the only mutation significantly associated with an increase in minimum inhibitory concentration (MIC) to fluoroquinolones that was found outside this region. 9 We assume that mutations cause resistance by reducing the affinity of an antibiotic ligand for its target. Since we are only interested in whether a mutation increases or decreases the antibiotic's affinity for the target, the difference in binding free energy (ΔΔG) between the wild type and mutant systems is calculated. This can be achieved by employing relative binding free energy (RBFE) methods, whereby a wild type amino acid is transmuted into the mutant along a nonphysical pathway defined by a progress coordinate, 0 ≤ λ ≤ 1. RBFE approaches have so far been successfully applied in small molecule drug design where the effect of perturbations to a lead compound on the binding affinity to its target can be predicted to within 1 kcal mol À1 error. [16][17][18] The regular nature of amino-acids compared to small molecules means that forcefields are well parameterized for amino acids 19 but potentially less so for small drug-like molecules 20 and therefore applying RBFE to predict the effects of amino acid mutations may be more accurate. Indeed, RBFE methods have been shown to yield accurate resistance predictions for genetic mutations associated with disease. [21][22][23][24] However, these studies have focused on small, monomeric protein targets, for example, we previously used RBFE methods to successfully predict trimethoprim resistance associated with mutations in the Staphylococcus aureus dihydrofolate reductase protein, which comprises 157 residues. 21,25 In this paper, to assess how well the method can be applied to much larger systems, we shall apply the same approach to two large protein complexes, the RNA polymerase (4671 residues) and the DNA gyrase cleavage complex (1473 residues), to assess how well we can predict the effect of seven and five mutations on the action of rifampicin and moxifloxacin, respectively. We emphasize that we define success as the ability of the method to rapidly predict whether each mutation confers resistance or not to the relevant drug. This qualitative approach is different to most RBFE studies which instead assess F I G U R E 1 Structures of Mycobacterium tuberculosis (A) RNA polymerase (RNAP) 28 and (B) DNA gyrase (DNAG) 27 cleavage complex, showing the selected clinical mutations associated with antibiotic resistance and susceptibility relative to the antibiotic binding sites. For clarity, RNAP subunits (excluding rpoB) are shown in surface view and nucleic acids are hidden in close-up visualizations. Resistance-conferring mutations are drawn in red, those associated with susceptibility blue and those residues where different mutations confer different resistance phenotypes purple. An asterisk (*) indicates a gyrB mutation.
how well the method can calculate the quantitative change in binding free energy (which are unknown for these proteins in any case). To retain the ability for the method to rapidly return results in any future implementation, we have used standard parameters (e.g., for ions), even if this potentially reduces our accuracy or precision.

| Selecting the mutations studied
To test the ability of RBFE to predict antibiotic resistance we selected a small number of mutations in the RNAP and DNAG that confer resistance; to act as negative controls we added several more F I G U R E 2 Free energy cycles for (A) rifampicin binding RNAP and (B) moxifloxacin binding DNAG gyrase cleavage complex. The subscripts qoff , vdW , and qon describe the process of first removing the electrical charge from atoms being perturbed, followed by transforming their van der Waals parameter, before finally recharging the atoms being perturbed. Double headed arrows represent the restraint used to prevent rifampicin from leaving the binding pocket. In all cases we are making use of the fact that free energy is a state function and therefore we can write the difference binding free energy (ΔΔG binding ) as a sum of so-called alchemical free energies (e.g., ΔG 4 -ΔG 3 )  Figure 1A). rpoB I491F is one of the socalled "disputed" mutations which either have variable or borderline rifampicin minimum inhibitory concentrations (MICs). 12,26 For moxifloxacin we also tested E501D in gyrB 12 which is close to the antibiotic binding site but not in the gyrA QRDR (Table 1, Figure 1B).
When choosing negative controls, we prioritized mutations that were observed multiple times in clinical samples, are close to the drug binding site and do not involve a charge change or a proline residue.
For the RNAP, L443F was selected since it lies within the RRDR and is close to the rifampicin binding site yet does not confer resistance 11 and therefore is a good negative control (Table 1, Figure 1A). We also selected S388L and T585A which are further from the binding site and are seen in clinical samples. Finally, we chose an amino acid (Ser428) at which non-synonymous mutations are expected to confer resistance, since it lies in the RRDR, but for which no firm statistical association has yet been made, and chose a mutation (S428C) which minimally chemically perturbs the sidechain. We expect this to not confer resistance, since it has not been observed clinically and the sidechain points away from the drug and S428C is therefore a good, if somewhat artificial, negative control.
For the DNA gyrase negative controls, we chose gyrA S95T ( Figure 1B) since it is very common-it is found in almost all samples except the H37Rv reference genome-and is within the QRDR. Testing different mutations at the same position which have different effects is a particularly stringent test of the ability of RBFE methods to predict antibiotic resistance. We therefore also tested the gyrA A90S mutation ( Figure 1B)this is not seen clinically but a serine is present at the equivalent position in the DNA gyrase of other bacterial species and is suggested to help stabilize the gyrasefluoroquinolone complex via participation in water-ion bridging interactions with the drug coordinated Mg 2+ . Mycobacterium tuberculosis has some innate immunity to fluoroquinolones which has been suggested is due to the alanine at this position. 27 The gyrA A90S mutation is therefore expected to strengthen the binding of moxifloxacin, thereby conferring hyper-susceptibility. Tyr as found in the forcefield. Due to the "loss" of the hydrogen atom from the tyrosine the system was left with a non-integer charge, so a solvent chloride ion was modified to provide a balancing charge. This is clearly not optimal, however we felt it preferable to manually redistributing charge. Mutations were then introduced into each of these structures using pmx. 33 To reduce the likelihood of clashes between the "new" sidechain and the remainder of the protein (i.e., in simulations with λ $ 1) we then applied a short Alchembed procedure 34 to each structurethis involved a 1000 step simulation where λ was increased from 0 to 1 using a soft-core van der Waals potential. For the RNAP simulations, we found it necessary to then apply a very short 0.5 ps simulation during which the temperature was increased to 310 K, presumably because this system had not been subjected to same degree of equilibration as DNAG. This created a pool of presumed independent mutated structures that could be used to seed alchemical thermodynamic integration simulations.
Following best practice, 35  were then also calculated. In hindsight, rifampicin would likely have remained bound without a restraintsince the restraint was not necessary to keep the drug bound the free energies calculated for removal of the restraint were minimal (supplementary file "rpob_summary.csv"), however, to ensure consistency between repeats we kept the restraint throughout.
Files and scripts necessary to reproduce the above steps, starting with the alchembed step, for the gyrA A90V DNAG and rpoB S450L RNAP mutations can be found here: https://github.com/fowler-lab/ tb-rbfe-setup.

| Calculation of errors
In previous studies of S. aureus DHFR 21,25 all alchemical free energy calculations were repeated the same number of times which, since n values of the final difference in binding free energy (ΔΔG) were then obtained, simplified the calculation of errors. Both simulation unit cells studied here were over an order of magnitude larger and we therefore instead calculated the SEM at the level of each individual alchemical free energy (e.g., ΔG vdW ), with the final error in ΔΔG estimated by adding these in quadrature. Throughout a 95% confidence limit was estimated by multiplying the SE by the appropriate t-statistic. We arbitrarily decided that at least three independent values of each alchemical free energy would be calculated, and then additional repeats would be run with the aim of reducing the magnitude of the overall 95% confidence limit to less than 1 kcal mol À1 . Achieving the latter was not always possible even when large numbers of repeats were run (n ≥ 10, see Supplementary Information S1). To avoid equilibration effects, the first 0.25 ns of each λ simulation was discarded before measuring the ΔG using thermodynamic integration calculated with the trapezoidal rule. Specifically, we used the numpy.

| Simulations run
trapz tool from NumPy for the integration calculation. Since the decision to discard 250 ps from the start of each simulation was arbitrary, we investigated whether discarding different amounts of data altered the results for one of the simple susceptible mutations, gyrA S95T, and the most complex mutation, gyrA D94G. Discarding either 125 or 375 ps, compared to 250 ps, led to no significant difference in the final calculated ΔΔG measurement for gyrA S95T or gyrA D94G ( Figure S1).

| Predictions
The simplest approach is to assume that a positive value of the change in binding free energy of the antibiotic (ΔΔG > 0) indicates that the antibiotic binds less well to the target following the mutation and therefore would be predicted to confer resistance to that drug. Clinically, however, a sample is categorized as "resistant" if its minimum inhibitory concentration (MIC) is greater than a critical con- For rifampicin, only one of the four negative controls (S388L) was correctly predicted to have no effect on the action of rifampicin ( Figure 3A); since the confidence limits of S428C, L443F, and T585A all bracket the ECOFF threshold no definite prediction could be made for these mutations. Clinically the method as implemented would therefore return an "Unknown" phenotype for these mutations. All three rifampicin-resistance conferring mutations, including the disputed mutation I491F, not only have positive values of ΔΔG but also lie above the clinical threshold derived from the ECOFF/ECV. These mutations are therefore correctly predicted to confer resistance to rifampicin.
Both moxifloxacin negative controls (gyrA S95T and A90S) were correctly predicted to not affect the binding of moxifloxacin to the DNA gyrase ( Figure 3B). Although hyper-susceptibility is expected for A90S, the magnitude of the confidence limits prevents us drawing any conclusions. No definite prediction could be made for any of the three mutations associated with moxifloxacin resistance since the confidence limits of all three mutations straddled the clinical threshold.
Unlike the RNA polymerase, two of the mutations to the DNA gyrase involved charged residues ( gyrB E501D and gyrA D94G) and not surprisingly these had the largest estimated errors.
To see how our ΔΔG values compared with clinical resistance measurements, we calculated an estimated "expected ΔΔG" from the geometric mean of MICs associated with each of the resistance conferring mutations, using previously described methods. 21 However, the errors in both the "expected ΔΔG" and the ΔΔG values calculated by RBFE were too large to enable us to draw any conclusions about how well the values compare with one another ( Figure S2). We also note that the expected free energy for the resistance conferring rpoB I491F mutation does not cross the clinical threshold we applied for resistance. Although the mutation is accepted as conferring resistance, 37 we had few clinical examples of this mutation, several of which had low MICs leading to a low expected ΔΔG value.  for the susceptible mutation rpoB L443F, which is also in the RRDR ( Figure 1A) and, whilst this also involves the introduction of a larger sidechain, in the crystal structure this is directed away from rifampicin.

| Investigation of sources of error
Differences in ΔG vdW between apo-and complexed DNA gyrase also appear mainly responsible for the positive value of ΔΔG for gyrA A90V, however the net effect is reduced.
Hence the variation in ΔΔG arises mainly from the apo and complexed values of ΔG vdWthe notable exceptions being gyrB E501D and gyrA D94G. This is despite our efforts to minimize the overall error by running up to 4Â the number of repeats for those transitions (Table 2) to reduce their individual estimated errors. For gyrB E501D and gyrA D94G all three transitions contribute significant error, which since they add in quadrature, leads to a large overall error in ΔΔG.
This is not surprising since both mutations involve turning off (and on) electrical charge and D94G involves a net charge change that must be compensated for elsewhere in the system. To investigate how far we might reduce the errors, let us now consider the individual values of ΔG qoff , ΔG vdW and ΔG qon ( Figure 5).
By starting each simulation from a different structural seed and discarding the first half of the alchemical free energy simulations and then applying statistics to the resulting values of ΔG we are assuming that they are independent. If true, then one would also expect the values to be normally distributed which would appear to be the case for most sets of ΔG values ( Figure 5). Applying the Shapiro-Wilks test of normality to the rpoB data confirms that, despite the small numbers  Figure S4b). In general, the individual λ values are much larger in magnitude for the gyrA D94G mutation than S95T, which is not unexpected due to the magnitude of change, but may explain why there is larger error associated with the gyrA D94G prediction ( Figure 3B).
The accuracy and precision of the free energy calculations will also be dependent on the convergence of the individual calculations.
Convergence of the calculations can be judged by the amount of overlap between energy distributions from neighboring λ windows, with overlap of >1% being considered reliable. 38 For the complex D94G mutation we found instances of zero overlap between neighboring λ windows for both the ΔG vdW and ΔG qon calculations and therefore these calculations were not converged (Table S1). This is not unsurprising given the use of very short calculations but will affect the accuracy of quantitative results and could affect the accuracy and confidence intervals for the overall qualitative prediction.
Conversely, for the less complex DNAG mutations there was no issues with the overlap (overlap data for S95T and the resistance conferring mutation A90V are shown in Table S1, data for A90S and E501D not shown).

| DISCUSSION
We have shown how relative binding free energy (RBFE) techniques can be applied to large protein complexes to predict, with some success, the effect of individual protein mutations on the binding of an antibiotic, and thence whether resistance is conferred. When the size of the signal is large and the mutations do not involve significant changes in the electrical charge, as is the case for the rpoB mutations, one can successfully predict whether a mutation confers resistance to the antibiotic (in this case rifampicin). If the fold increase in minimum inhibitory concentration is small and/or there are significant charge changes, as is the case for most of the resistance conferring mutations in the DNA gyrase, then the estimated error of ΔΔG will likely be so large that no definite prediction can be made. In addition, the observed non-normality and lack of overlap between λ windows for constituent free energies for gyrA D94G indicates that these values are also not independent and not converged: to solve this either the λ simulations would have to be extended or the equilibration simulations would need to be more numerous as well as longer.
Despite the focus on resistance, it is more useful to be able to accurately and reproducibly predict susceptibility since clinically that will lead to immediate action, that is, starting the patient on the appropriate treatment regimen. A prediction of resistance will likely result in the sample being sent for further testing, at which point any incorrect predictions (false positives) will be detected. Due to the magnitudes of the errors and the fact that we are constrained to start from the wild-type structure, it may be more difficult, unfortunately, to predict susceptibility than resistance using RBFE. For rpoB only one susceptible mutation could be confidently predicted. If we assume that most susceptible mutations will not affect the binding affinity for the antibiotic, then they would have a ΔΔG of zero. The predicted ΔΔG of such mutations would hence require the estimated error to be at least less than the value of the ECOFF (for rpoB and DNA gyrase 1.2 and 0.9 kcal mol À1 , respectively) to make a confident susceptible prediction. The magnitude of error for the mutations in this study was greater than the relevant ECOFF in all but one case (rpoB S388L, Finally, we assume that the conformations used to seed each calculation are independent of one another and/or that the λ simulations are long enough to allow the initial state to be "forgotten". Given we chose to use very short λ simulations the latter is almost certainly not true and whilst the majority of our calculated ΔG values appear to be normally distributed, some are not which is concerning. One would expect to have to run 4Â the number of simulations to reduce the estimated error to half its original value if the simulations are independent. This makes simulations of these size prohibitively computationally intensive using the software and compute that we employed. Significant speed up could however be achieved by use of updated GROMACS software, where yearly updates have been shown to increase ns/day performance on a range of computational resource. 39,40 It is not in doubt that how the structure and dynamics of a protein change upon mutation contains valuable information that can, in theory, be used to predict whether individual mutations confer antibiotic resistance. Although RBFE may not (yet) be an appropriate tool for resistance prediction in the two complex systems we studied, several alternative routes exist. For instance, machine learning algorithms or energy scoring software such as Rosetta can predict binding free energies and these methods are less computationally intensive. 22,41,42 Ultimately a combination of RBFE/MD and other approaches may not only complement one another but also form part of a larger toolkit that helps us to tackle antimicrobial resistance by improving diagnosis.

DATA AVAILABILITY STATEMENT
The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials.