{'fraction_non_alphanumeric': 0.022854200499069197, 'fraction_numerical': 0.05489761159741752, 'mean_word_length': 5.057581573896353, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A computational method for predicting the most likely 1 evolutionary trajectories in the stepwise accumulation of 2 resistance mutations 3

31 Pathogen evolution of drug resistance often occurs in a stepwise manner via accumulation of 32 multiple mutations that in combination have a non-additive impact on fitness, a phenomenon 33 known as epistasis. Epistasis complicates the sequence-structure-function relationship and 34 undermines our ability to predict evolution. Here we test how well Rosetta Flex ddG is able to 35 estimate epistasis between resistance mutations and observe good agreement with experimental 36 data. By using these predictions to parameterize an existing mechanistic model of molecular 37 evolution, we also investigate how well such a model is able to capture evolution to resistance 38 where the order of fixation is influenced by epistasis. We apply this method to predict evolutionary 39 trajectories to known quadruple mutations associated with resistant phenotypes for malaria. 40 Resistance to the combination drug sulfadoxine-pyrimethamine (SP) in malaria-causing species 41 Plasmodium falciparum (Pf) and Plasmodium vivax (Pv) has arisen via the accumulation of 42 multiple point mutations in the DHFR and DHPS genes. The evolution of resistance in the DHFR 43 genes of Pf and Pv has been studied extensively. Four known Pf DHFR pyrimethamine resistance 44 mutations are highly prevalent in field-isolates and multiple studies have shown epistatic 45 interactions between these mutations determine the accessible evolutionary trajectories to the 46 highly resistant quadruple mutation N51I,C59R,S108N,I164L. We simulated the possible 47 evolutionary trajectories to this quadruple Pf DHFR mutation as well as the homologous Pv DHFR 48 mutations. In both cases, the most probable pathways agreed well with those determined 49 experimentally. This method can be applied to any drug-target system where the drug acts by 50 binding to the target. resulting of antimicrobial drugs. Our ability to monitor emerging resistance and direct the most appropriate treatment strategy (stewardship) would be strengthened by the development of methods to accurately predict the mutational pathways to resistance. Here, we present a computational method to predict the most likely pathways to resistance that result from combinations of mutations that alter drug binding affinity. We utilized the Rosetta Flex ddG protocol to parameterize a thermodynamic evolutionary model to predict the order of fixation of each mutation in the pathway. This approach accurately captures known resistance trajectories and can be applied to any system where resistance arises via drug binding affinity changes

combinations of these four mutations [6,12]. This means that mutations which on their own are 140 not associated with a resistance phenotype, can be when in combination with other mutations.

141
Epistasis between these mutations has been shown to determine the evolutionary trajectories to 142 the quadruple mutation N51I,C59R,S108N,I164L, which is strongly associated with 143 pyrimethamine resistance [6].

167
The development of computational methods to predict resistance trajectories would enable fast 168 and efficient predictions and would be more widely accessible than lab-based methods.

169
Computational tools could help narrow down the pool of mutations to be studied experimentally 170 and would also be applicable to difficult to study targets. Some target-specific computational tools 171 to predict individual resistance mutations have been developed [36,37]. However, such tools are 172 target specific and so not generalizable. Furthermore, they only consider independent mutations 173 on a single structure and so ignore epistasis between resistance mutations. Therefore, they are 174 not suitable for predicting evolutionary trajectories to resistance.

176
To determine a generalizable computational method to predict evolutionary trajectories to 177 resistance, we need to consider the main determinants of resistance. Rodrigues et al. (2016) [38] 178 investigated three mutations in Escherichia coli DHFR associated with trimethoprim resistance 179 and considered activity, binding affinity, fold stability, and intracellular abundance. They found 180 that whilst resistance is a trade-off between these factors, binding affinity is the single most 181 predictive trait of resistance, especially at later points in evolution. Therefore, we decided to 182 investigate if predictions of binding affinity change can be used to predict the order of fixation of 183 resistance mutations involved in evolutionary trajectories to resistance.

185
Rosetta Flex ddG [39] is the current state-of-the art method for predicting changes in protein-186 protein and protein-ligand binding free energy. Rosetta is a software suite for macromolecular 187 modelling and design that uses all-atom mixed physics-and knowledge-based potentials, and 188 provides a diverse set of protocols to perform specific tasks, such as structure prediction, 189 molecular docking and homology modelling [40]. The Flex ddG protocol has been found to 190 perform better than machine learning methods and comparably to molecular dynamics methods 191 when tested on a large dataset of ligand binding free energy changes upon protein mutation [41, 192 42]. However, its ability to capture epistasis has not yet been tested. Therefore, we investigated 193 how well Flex ddG can capture epistasis between resistance mutations in PfDHFR and observed 194 a good agreement with experimental data.

196
Next, we used the Flex ddG predictions to parameterize a fitness function and applied this fitness 197 function in an existing model of molecular evolution. We used this method to predict evolutionary 198 trajectories to known resistant quadruple mutants in both PfDHFR and PvDHFR, where the 199 evolutionary trajectories have been studied experimentally [6,7]. Good agreement was observed 200 between the most likely trajectories to the quadruple mutations predicted by our model and those 201 predicted experimentally.

203
The main advantage of this approach is that it does not require access to an experimental 'wet' 204 lab and can be carried out by anyone with access to a high-performance computer (to run the 205 Flex ddG predictions). It is generalizable to any antimicrobial drug that acts by binding to its target 206 and can be easily applied to any drug-target complex for which there is an available structure.

207
Therefore, it can be used to study complexes and systems that might be problematic 208 experimentally. It enables accurate assessment of the predictability of the evolutionary 209 landscape, in other words it can predict whether we would expect to see constrained evolutionary 210 trajectories on a fitness landscape as a result of epistatic interactions in drug binding free energy.

212
In addition, we analyzed if evolutionary pathways can be inferred from the frequency of mutations 213 found in isolate data. We determined the frequency of mutations in PfDHFR and PvDHFR, and 214 inferred the most likely evolutionary pathways under the assumption that the most likely mutation 215 at each step corresponds to the most frequent mutation. We carried out this analysis first upon a 216 combined set of global isolates from all regions and then upon isolates from each region 217 individually. In general, the most likely pathways inferred from the global isolate data agreed 218 remarkably well with both the experimentally determined pathways and the pathways predicted 219 by our computational method. When considering geographical regions separately, the inferred 220 pathways from several regions agreed well with the experimental pathways and our predicted 221 pathways, however the most likely pathways inferred in some regions differed from the main 222 pathways, highlighting the importance of considering the evolution in different regions separately.

225
Rosetta Flex-ddG accurately predicts binding free energy changes and captures epistasis

226
We investigated if Flex ddG predictions agree with experimentally measured binding free energy 227 and if these predictions can be used to calculate the non-additivity in binding free energy 228 (interaction energy), which for a double mutant defines the epistasis between the two mutations 229 and, for a triple mutant or higher, captures the level of epistatic interactions. We calculated the 230 interaction energy by finding the difference between the predicted change in binding free energy 231 of a multiple mutation and the sum of the predictions of their independent binding free energy 232 changes. A positive value of the interaction energy indicates the sum of the independent impacts 233 is more destabilizing than the impact of the multiple mutation and a negative value indicates the 234 sum is less destabilizing than the combined impact.

236
The change in binding free energy was predicted using Rosetta Flex ddG for the combinatorically 237 complete set of the four PfDHFR pyrimethamine resistance mutations N51I, C59R, S108N and 238 I164L. (Note on notation: lists of single mutations are written X, Y, Z, multiple mutations are 239 written X,Y,Z (i.e. no space between the commas and the mutations) and pathways are written 240 X/Y/Z to denote the order of fixation).

242
We compared the predictions to the data from [12] in which they present the binding free energy 243 changes for a subset of the possible combinations of mutations, the sum of the independent 244 mutations (calculated for multiple mutants to compare to the experimentally determined binding 245 free energy changes of multiple mutants) and the interaction energy of the multiple mutants 246 (

254
The authors of the Flex ddG protocol suggest conducting a minimum of 35 runs and to take the 255 average of the distribution as the prediction for that mutation [39]. To determine the optimum 256 number of runs per mutation, we ran predictions for 35 ( Figure S1, Table S1) and 150 ( Figure S2, 257 Table 1) runs and a single prediction at 300 runs ( Figure S3), and concluded 150 runs provides 258 the best trade-off between accuracy and efficiency (see Supplementary text for detailed 259 discussion). Therefore, we will be discussing the predictions for n=150 going forward

261
Examining the distributions of the predicted change in binding free energy for each of the 262 mutations considered in [12] ( Figure S2) we can see that the distributions are not Gaussian but 263 have a defined peak, and so the average of the distribution may not be the best summary metric 264 to use to compare to the data. We therefore chose the peak of the distribution as the summary 265 metric of each prediction. We compared the predictions for 150 runs and the data from [12] 266 (Table 1) and observed a correlation of 0.635 between for the binding free energy data, 0.435 for 267 the sum of the independent predictions for multiple mutants and 0.765 for the interaction energy.

268
We found 8/9 binding free energy predictions were correctly classified, 4/5 of the sum of the 269 independent predictions were correctly classified and 4/5 of the interaction energies were 270 correctly classified.

272
Mutation S108N was the only single mutation to destabilize pyrimethamine binding in both the 273 experimental data and the Flex ddG predictions. However, in the experimental data the double 274 mutation N51I,S108N is more destabilizing to binding than single mutation S108N, but the peak 275 of the Flex ddG distribution was stabilizing. The triple mutation C59R,S108N,I164L was found 276 experimentally to be the most destabilizing of the triple mutations, however the peak of the Flex 277 ddG distribution predicted it to be only mildly destabilizing and the least destabilizing of the triple 278 mutations. Furthermore, the quadruple mutation was found experimentally to have the most 279 destabilizing impact out of all combinations of single and multiple mutations, however, Flex ddG 280 predicted it to be less destabilizing than the double mutation C59R,S108N.

282
Considering the interaction energy, the incorrectly classified mutation was again N51I,S108N,

283
where the interaction energy was predicted to be positive, but found experimentally to be 284 negative, because the sum of the individual predictions was destabilizing but the double mutation

290
The triple mutation C59R,S108N,I164L was found to have a large positive interaction energy in 291 both the experimental data and our predictions, where the triple mutation was found to be 292 destabilizing whilst the sum of the individual impacts was stabilizing to pyrimethamine binding.

293
We also observed large negative interaction energy between S108N and C59R, where C59R is 294 stabilizing in the wildtype background but destabilizing in the background of S108N, an example 295 of sign epistasis and in agreement with the observations of both [12] and [6]. However, whilst the 296 interaction energy of the triple mutation N51I,C59R,S108N was positive for both the experimental 297 data and predictions, in our predictions its magnitude was much smaller compared to the data.

298
Both single mutations N51I and C59R were predicted to be only marginally stabilizing -almost 299 neutral -to pyrimethamine binding, whilst in the experimental data both mutations have a large 300 stabilizing impact. Furthermore, the triple mutation was predicted to be only marginally more 301 destabilizing than single mutation S108N, resulting in the small negative interaction energy.

303
We use these predictions in the following sections to parameterize an evolutionary model to

327
To determine how well our simulations reflect the evolutionary process to the PfDHFR quadruple 328 mutation N51I,C59R,S108N,I164L, we compared our results to experimentally determined 329 evolutionary trajectories as presented by [6]. In our simulations, the quadruple mutation was 330 reached in approximately 12% of runs (see file 'PfDHFR_endpoint_numbers.csv' in the 331 Supplementary data). The majority of the runs (75%) terminated at a double mutation, with 332 S108N/C59R the most likely trajectory over all. The algorithm was often unable to proceed 333 passed S108N/C59R in many runs because the Flex ddG distribution for C59R,S108N is likely trajectory to the quadruple mutation in [6] was predicted to be S108N/N51I/C59R/I164L, 350 however this pathway was predicted to be very unlikely in our simulations, due to the fact that the 351 distribution of Flex ddG predictions for double mutation N51I,S108N was mostly stabilizing to 352 pyrimethamine binding ( Figure S2e), whereas all of the S108N distribution was destabilizing to 353 pyrimethamine binding, so this step was unlikely to be chosen by the evolutionary algorithm.

355
Considering the frequency at which the single mutations were chosen as the first step in all 356 simulated pathways (see file 'PfDHFR_total_pathway_probabilities.csv'), S108N was the most 357 likely single mutation and C59R was the second most likely single mutation, in agreement with 358 the two most likely first steps in the pathways predicted in [6]. In our simulations, N51I was the 359 third most likely single mutation and I164L has zero probability of being chosen as the first 360 mutation, due to the fact that all 150 predicted binding affinity changes were stabilizing to 361 pyrimethamine binding. The most likely pathway to a double mutation realized in all trajectories in 362 both our simulations and the simulations in [6] is S108N/C59R. Similarly, the most likely pathway 363 to a triple mutation realized in all our simulations and in [6] was S108N/C59R/N51I.

365
To simulate the evolutionary pathways for PvDHFR, we also carried out predictions of binding 366 free energy changes for the homologous set of four mutations in PvDHFR, (N50I, S58R, S117N 367 and I173L

380
However, the triple mutation S58R,S117N,I173L was predicted to have a smaller destabilizing 381 impact than S58R,S117N, making pathway S58R/I173L/N50I unlikely, whilst N50I,S58R,I173L 382 was predicted to be stabilizing to pyrimethamine in all Flex ddG runs and therefore pathway 383 S58R/I173L/N50I had a zero probability of occurring in the simulations. This resulted in many 384 runs terminating at step S58R/I173L.

452
Considering the total frequency of each mutation in the set of four PfDHFR mutations (N51I,

491
Single mutation S58R was the second most frequent single mutation in our isolate data (20/847).

492
This supports the predicted first evolutionary steps in [7] which were predicted to be S58R for the 493 highest pyrimethamine concentration and S117N for three lower concentrations of 494 pyrimethamine, which suggests both single mutations are possible, but S117N is more likely for a 495 lower pyrimethamine concentrations. An alternative pathway to the triple mutation could therefore 496 be S58R/S117N/I173L, which may be more likely under higher pyrimethamine concentrations.

497
This corresponds to the most likely pathway to a triple mutation in our simulations and is part of

524
As in the previous section, we inferred the most likely pathway by assuming the most frequent 525 mutation at each step corresponds to the most likely evolutionary trajectory. Considering the 526 frequency of mutations found in PfDHFR isolates from the different geographical regions ( Figure   527 3), the most likely evolutionary pathway to the quadruple mutation N51I,C59R,S108N,I164L was 528 inferred to agree with the main pathway, S108N/C59R/N51I/I164L, predicted by our evolutionary model and the data presented in [6], as well as the most likely pathway inferred by considering 530 the frequency of all worldwide PfDHFR isolate data in West Africa (S108N: 31/2594; 531 C59R,S108N: 211/2594; N51I,C59R,S108N: 1739/2594; N51I,C59R,S108N,I164L: 1/2594),

543
However, in East Africa S108N/N51I/C59R/I164L was the most likely pathway to the quadruple 544 mutation (S108N,N51I: 67/904), which corresponds to the fifth most likely pathway to the 545 quadruple in our simulations and the third most likely pathway presented in [6] and inferred from 546 the total frequency data. This was also an alternative pathway in West Africa (S108N,N51I: 547 32/2594) .

549
Furthermore, the quadruple mutation was not observed in the isolate data from Central Africa,

558
Melanesia, evolution up to the double mutation C59R,S108N was observed in the isolate data 559 and was inferred using the frequency data to have followed the pathway S108N/C59R (S108N: 560 1/119; C59R,S108N: 23/119), which corresponds to the most likely pathway of all our 561 evolutionary runs.

563
Next, we considered the frequency of combinations of the four PvDHFR mutations N50I, S58R, 564 S117N and I173L in different geographical regions and used these frequencies to infer 565 evolutionary trajectories ( Figure 4). As mentioned previously, the quadruple mutation is not 566 observed in the isolate data and so we will infer trajectories up to triple mutant combinations of 567 the four mutations where possible, and compare to the most likely pathways to triple mutations in 568 our simulations and in the experimental data. The triple mutation S58R,S117N,I173L was the only 569 triple mutant combination observed in the isolate data and was only found in South America.

570
Analyzing the frequency of the constituent mutations (S117N: 34/257; S58R,S117N: 74/257; 571 S59R,S117N,I173L: 2/257), the most likely inferred pathway to this mutation is 572 S117N/S58R/I173L. This was the third most likely triple mutation pathway in our simulations, 573 however it was not an observed pathway in [7]. There are two other possible pathways inferred 574 from the frequency data to this triple mutation: S117N/I173L/S58R (S117N,I173L: 1/257), and 575 S58R/S117N/I173L (S58R: 9/257). Pathway S117N/I173L/S58R was the fourth most likely 576 pathway to a triple mutation in our simulations when considering all possible pathways but was 577 not observed in the simulations in [7] and S58R/S117N/I173L the most likely pathway to a triple 578 mutation in our simulations and was observed as part of the second most likely pathway to the 579 quadruple mutation at the highest pyrimethamine concentration considered in [7].

581
Evolution only up to double mutation S58R,S117N was found in Southeast Asia, East Africa and

582
Melanesia and S117N/S58R was the most fixation order in East Africa, whilst S58R/S117N was 583 the most likely fixation order in Southeast Asia and Melanesia. Pathways S58R/S117N and 584 S117N/S58R were the second and third most frequently observed double mutation pathways 585 observed in our simulations when considering all possible pathways observed. In South Asia 586 evolution up to double mutations S58R,S117N and N50I,S117N was observed, following 587 pathways S117N/S58R and S117N/N50I, respectively. In East Asia, evolution only up to single 588 mutation S58R was observed and in Central America no steps in the evolutionary pathway 589 including combinations of these four mutations were found.

613
Indeed, Rodrigues et al. (2016) [38] observed a clear tradeoff between catalytic activity and binding affinity for increased drug resistance and found that whilst many molecular features affect 615 drug resistance, drug binding affinity was the key determinant of drug resistance in later stages of 616 evolution. This may be why our predictions agree well their predictions for high pyrimethamine

643
Mutations occurring at a drug-binding site may also reduce the protein's thermodynamic stability

644
[50] and therefore may not be selected for, even if they improve the resistance phenotype.
645 Therefore, our model may also be improved by including selection for mutations that do not 646 reduce thermodynamic stability relative to the wild-type enzyme. However, it must also be noted

716
This highlights the need to distinguish between geographical regions and account for existing 717 resistance alleles within that region and trace their lineages when attempting to predict the next 718 step in evolutionary trajectories to highly resistant multiple mutants. Given the current dominant 719 resistance allele from a specific region, our method could be used to predict the most likely next 720 steps from a subset of likely mutations.

722
We have presented a computational method for predicting the most likely evolutionary trajectories 723 that has demonstrated good agreement with trajectories predicted experimentally and has the 724 advantage of being much quicker and more cost-effective than laboratory-based methods. This 725 method can be applied to any system in which a drug binds to a target molecule, provided a

742
To produce a complete structural model of PvDHFR, PDB entry 2BLB was used as a template,

743
which provides the X-ray crystal structure of wild-type P. vivax DHFR in complex with 744 pyrimethamine. This structure was only missing a loop section between residues 87-105 and so

745
Modeller was used to build this missing loop.  or destabilising (Barlow et al. 2018 [39]). Therefore, we assume Flex ddG can provide 796 approximate predictions of binding free energy changes comparable to experimental changes in 797 kcal/mol and can therefore be used to parameterize a thermodynamic model.

799
To predict the most likely evolutionary trajectories to reach a quadruple mutant we used a model 800 based in thermodynamics and statistical mechanics where the fitness of a protein is determined 801 by the probability it would not be bound to a ligand, . We consider a two where c0 is 807 a reference ligand concentration (set here arbitrarily to 1M), ∆ is the protein-ligand binding free 808 energy, k is the Boltzmann constant and T is the temperature in Kelvin.
Starting from the wild-type protein, with binding free energy ∆ and fitness , we extract 811 one sample i from the 150 values of the predicted binding affinity changes for the single 812 mutations to determine the binding free energy after mutation X, ∆ = ∆ + ∆∆ , and 813 calculate the fitness of each single mutant protein ( ) . We can calculate the probability the 814 mutation will fix in the population using the Kimura fixation probability for a haploid organism where Ne is the effective population size (set to 10 as previous models in [43,44]) and s is the probabilities of all mutations sampled at a step are effectively zero, no mutation is chosen at that 827 step and the algorithm begins again by choosing a single mutation. Therefore, not all of the runs 828 produce a complete trajectory and some will terminate before reaching the quadruple mutation.

829
The algorithm was written in R.

831
We calculate the probabilities, , of each realized pathway (even those which don't reach the 832 quadruple mutation) by dividing the total number of times that specific pathway occurs by the 833 number of runs. We calculate the probability of a particular step by dividing the number of times 834 that step occurs in all realized pathways by the total number of runs.  London) for introducing her to the evolutionary algorithm used in this paper.

865
Data Availability