Rational Design of Antibiotic Treatment Plans: A Treatment Strategy for Managing Evolution and Reversing Resistance

The development of reliable methods for restoring susceptibility after antibiotic resistance arises has proven elusive. A greater understanding of the relationship between antibiotic administration and the evolution of resistance is key to overcoming this challenge. Here we present a data-driven mathematical approach for developing antibiotic treatment plans that can reverse the evolution of antibiotic resistance determinants. We have generated adaptive landscapes for 16 genotypes of the TEM β-lactamase that vary from the wild type genotype “TEM-1” through all combinations of four amino acid substitutions. We determined the growth rate of each genotype when treated with each of 15 β-lactam antibiotics. By using growth rates as a measure of fitness, we computed the probability of each amino acid substitution in each β-lactam treatment using two different models named the Correlated Probability Model (CPM) and the Equal Probability Model (EPM). We then performed an exhaustive search through the 15 treatments for substitution paths leading from each of the 16 genotypes back to the wild type TEM-1. We identified optimized treatment paths that returned the highest probabilities of selecting for reversions of amino acid substitutions and returning TEM to the wild type state. For the CPM model, the optimized probabilities ranged between 0.6 and 1.0. For the EPM model, the optimized probabilities ranged between 0.38 and 1.0. For cyclical CPM treatment plans in which the starting and ending genotype was the wild type, the probabilities were between 0.62 and 0.7. Overall this study shows that there is promise for reversing the evolution of resistance through antibiotic treatment plans.


Introduction
Antibiotic resistance is an inevitable outcome whenever antibiotics are used. There are many reasons for this: 1) As humans (also as eukaryotes), we are vastly outnumbered by bacteria in nearly all measures, including total population size, biomass, genetic diversity, emigration, and immigration [1]; 2) bacteria can use horizontal gene transfer to share resistance genes across distantly related species of bacteria, including non-pathogens [2]; 3) compared to humans, bacteria have relatively few vulnerable target sites [3]; 4) microbes are the sources of nearly all antibiotics that are used by humans [4]. Given the overwhelming numbers of bacteria, the limited number of target sites, the numerous ways that they can infect humans, and that they have been exposed to naturally occurring antibiotics for billions of years, resistance to antibiotics used by human populations is unavoidable.
Once resistance is present in a bacterial population, it is exceedingly difficult to remove for several reasons. If any amount of antibiotic is present in the environment, antibiotic resistance genes will confer a large fitness advantage [5], and even when antibiotics are not present in an environment, the fitness costs for carrying and expressing resistance genes are small to non-existent [6]. In addition to it being difficult to remove antibiotics from the environment [7], even if humans were to completely abandon the use of antibiotics, resistance would persist for years [8].
Efforts to remove resistance genes from clinical environments by either discontinuing or reducing the use of specific antibiotics for some period of time, either through general reduction of antibiotic consumption or periodic rotations of antibiotics (cycling) have not worked in any reliable or reproducible manner [9]; indeed it would have been surprising if they had worked [10], [11].
Since antibiotic resistance is unavoidable, it only makes sense to accept its inevitability and develop methods for mitigating the consequences. One reasonable approach is to rotate the use of antibiotics. This has been implemented in many ways and there are recent studies to model the optimal duration, mixing versus cycling, and how relaxed antibiotic cycles may be and still function as planned [12,13]. However, none of those models have focused on developing a method for designing an optimal succession of antibiotics.
In a previous publication [14], we proposed that susceptibility to antibiotics could be restored by rotating consumption of multiple antibiotics that are a) structurally similar, b) inhibit/kill bacteria through the same target site, and c) result in pleiotropic fitness costs that reduce the overall resistance of bacteria to each other. We presented a proof-of-principle example [14] of how this might work with a series of β-lactam antibiotics in which some of them would select for new amino acid substitutions in the TEM β-lactamase and others that would select reversions in TEM ultimately leading back to the wild type (un-mutated) state. We have focused particularly on β-lactamases because there is often no fitness cost associated with their expression, and they are particularly difficult to remove from clinical microbial populations.
Our current work seeks to identify β-lactam treatment plans that have the highest probability of returning a population expressing a small number of variant TEM genotypes to the wild type state. The wild type TEM-1, and a handful of its descendants, confers resistance to penicillins alone. However, most of the descendants confer resistance to either cephalosporins or penicillins combined with β-lactamase inhibitors (inhibitor resistance), and a few confer resistance to both. Of the 194 clinically identified TEM genotypes that encode unique amino acid sequences [15], 174 (89.7%) differ from the wild type TEM-1 by at most four amino acid substitutions (see Table 1). Our choice of a system that includes four amino acid substitutions is based upon an apparent threshold for amino acid substitutions among functional TEM genotypes. The rarity of the co-existence of cephalosporin resistance and inhibitor resistance and the fact that no single substitution confers both phenotypes suggested that sign epistasis (i.e. reversals of substitutions from beneficial to detrimental) exists as the substitutions that contribute to this dual phenotype are combined. We have assumed that substitutions arise according to the strong selection weak mutation model (SSWM) [16] in which single substitutions reach fixation before the next substitution occurs. Recent work [17] in addition to past phylogenetic analysis [18] and competition experiments [19] suggest that this is a valid model for TEM evolution.
The ability to apply selective pressures that favor reversions of substitutions within an evolved TEM genotype would increase the number of antibiotics that could be used. To embark upon our effort of determining the best way to do this, we decided to create a model system based upon the TEM-50 genotype, which differs from TEM-1 by four amino acid substitutions. All four substitutions by themselves confer clearly defined resistance advantages in the presence of certain antibiotics. Additionally, TEM-50 is one of the few genotypes that simultaneously confers resistance to cephalosporins and inhibitor combined therapies.

From experimental data to mathematical models
We created all 16 variant genotypes of the four amino acid substitutions found in TEM-50 using site directed mutagenesis (Table 2) and measured the growth rates of 12 replicates of E. coli DH5α-E expressing each genotype in the presence of one of fifteen β-lactam antibiotics (Table 3). Each genotype was grown in each antibiotic in 12 replicates. We computed the mean growth rate of those replicates (Table 4) and the variance of each sample, as well as the significance between adjacent genotypes that differed by one amino acid substitution. This was done using one-way ANOVA analysis. The results are summarized in Figs 1-15, where the arrows in the fitness graphs connect pairs of adjacent genotypes. For each comparison of adjacent genotypes, we indicate the one whose expression resulted in the faster growth by directing the arrowhead towards that genotype, and implying that evolution would proceed in that direction if the two genotypes occurred simultaneously in a population [20,21]. In other words, the node indicated by the arrowhead would increase in frequency and reach fixation in the population, while the other would be lost. Red arrows indicate significance, and black arrows indicate differences that were not statistically significant by ANOVA, but that may still exist if a more sensitive assay was used.
Figs 1-15. These figures present fitness graphs, which are a visual summary of the adaptive landscape 2x2x2x2 tensors in which each resistance phenotype conferred by each TEM genotype is enumerated. Arrows pointing upward represent selection for the addition of a substitution. Arrows pointing downward represent selection for reversions. Red arrows indicate significance between adjacent growth rates as determined by one way ANOVA. Genotypes that confer the most resistance to each antibiotic are shown in red.
We rank ordered the genotypes (Table 5) in each landscape diagram with a score from 1 to 16, with the genotype promoting the fastest growth receiving a score of "1" and the genotype with the slowest growth a score of "16". This analysis shows that all genotypes have a score of 5 or better and a score of 13 or worse, in at least one landscape, indicating that there is abundant pleiotropy as antibiotic selective pressures change. That pleiotropy provides a basis for effectively alternating antibiotic to restore the wild type.
We considered the 15 antibiotics previously mentioned in Table 3: AMP, AM, CEC, CTX, ZOX, CXM, CRO, AMC, CAZ, CTT, SAM, CPR, CPD, TZP, and FEP and their interactions with a bi-allelic 4-locus TEM system {0,1} 4 where four functionally important amino acid residues involved in the evolution of TEM-50 are considered. The number "1" denotes an amino acid substitution, whereas "0" denotes no substitution at the site. We experimentally determined growth rates for all genotypes in our TEM system at a selected concentration of each antibiotic. Those growth rates depend upon the states of the four amino acid residues. The growth rates for all genotypes in one antibiotic can be represented by a real 2×2×2×2 tensor f = (f ijkl ), where f(a r ) is the fitness landscape for the antibiotic r. We can identify f(a r ) with a vector whose coordinates are indexed by {0,1} 4 . The resulting 15 vectors, one for each antibiotic, are the rows in Table 4. Our substitution model M(f) is a function M : R 16 ! R 16Â16 that assigns a transition matrix to each fitness landscape. The rows and columns of M(f) are labeled by the genotypes {0,1} 4 according to the degree lexicographic order. The entries in M(f(a r )) u,v represent the probability that that genotype u is replaced by genotype v under the presence of antibiotic a r . For that reason, the rows of our transition matrices have nonnegative entries and their rows sum to 1.
We require that our transition matrices respect the adjacency structure of the 4-cube, that is, M(f) u,v = 0 unless u and v are vectors in {0,1} 4 that differ in at most one coordinate. In other words, we reasoned that resistant strains are most likely to be in competition with those that express resistance genotypes that are immediately adjacent (vary by a single amino acid substitution).  Based on the strong patterns of pleiotropy we observed, we reasoned that the choice and the succession of antibiotics were at least as important as other cycling considerations. We formalized our approach to identifying optimal antibiotic treatment paths as follows. doi:10.1371/journal.pone.0122283.t005

Rational Design of Antibiotic Treatment Plans
The combined effect of a sequence a 1 ,. . .,a k of k antibiotics is described by a new transition matrix Mðf a 1 Þ Á Mðf a 2 Þ Á ::: Á Mðf a k Þ obtained as the product of the transition matrices for each drug.
For any genotype u other than 0000, our goal is to find a sequence of antibiotics which maximizes the probability of returning to the wild type. In other words, if we restrict to sequences of length k our goal is to find a sequence of antibiotics a 1 ,. . .,a k which maximizes the matrix entry M(f(a 1 ))ÁM(f(a 2 ))Á. . .ÁM(f(a k )) u,0000 . For each u this requires searching over all 15 k antibiotic sequences of length k.

Finding optimal sequences of antibiotics
We used two substitution models to determine the optimal (highest probability) sequences of β-lactams for returning TEM genotypes back to their wild type state. Briefly, the Correlated Probability Model (CPM) allows probabilities to be based upon the actual growth rates. It is given by applying Eq (3) to the growth rates in Table 4. The Equal Probability Model (EPM) assumes that beneficial substitutions are equally likely and that only the direction of the arrows in Figs 1-15 is important. This means that the matrix entry M(f) u,v is 1/N if genotype u has N outgoing arrows and there is an arrow from u to v.
A visual summary of the highest probabilities according to the CPM is provided in Fig 16. The CPM provides good estimates if fitness differences between genotypes are small [14,[22][23][24]. The EPM has been used in settings where only rank order (as in Table 5) is available [25].
From the graph, it is possible to find candidate treatment plans. For example, when starting at genotype 1010 the graph shows that the probability for ending at 0000 is 0.71for the sequence ZOX-TZP (0.71 is the product of the arrow labels). Similarly, when starting at 1111 the probability for ending at 0000 is 0.62 for the sequence CEC-CAZ-TZP-AM. When starting at 0001 the graphs shows that a single drug gives probability at most 0.29, whereas the probability for ending at 0000 for the sequence AMC-CRO-AM (one arrow up, two arrows down) is at least 0.96Á0.62 = 0.6.
This graph can also be used to generate treatment paths that start and end at the same genotype, making possible the development of a fixed treatment plan. For example, from a starting point 0000, the probability for ending at 0000 is 0.62 for the sequence: CEC-CTX-ZOX-CPD-CPR-CAZ-TZP-AM For all sequences of antibiotics of a fixed length (two, three, four, five, and six), we examined the probability that a given genotype is returned to the wild type state. It is worth noting that within these paths, a single genotype can be repeated consecutively with different antibiotics, thus making it possible to have an odd number of steps in the treatment paths when an even number of subtitutions are introduced. For every starting genotype, we found we were able to return to the wild type genotype with a probability between 0.6 and 1.0 when using the CPM model and a probability of 0.375 and 1.0 when using the EPM model. These results are summarized in Tables 6-9 and Fig 17. These results show the number of paths and their probabilities (Tables 6 and 7) and the substitutions selected through the optimal treatment plans (Tables 8-11) for returning to the wild type state from various starting points.
Once returned to the wild type state, we identified cycles that would allow for alternation of antibiotics, and allow for some variation through amino acid substitution, but then rapidly return bacteria to the wild type state (Table 12 and Fig 18). Such cycles were possible for path length of two, four, and six and the probabilities of those paths were respectively 0.704, 0.617, 0.617. We found that in the most probable cases, the genotype varied by only one amino acid Rational Design of Antibiotic Treatment Plans substitution before reverting back to the wild type state. However, when treatment plans with lower probabilities are considered, we find that more amino acid substitutions in the genotype are allowed.

Discussion
In this study, we have developed an experimental approach for measuring pleiotropy and a computational mathematical approach for optimizing antibiotic treatment paths. The experimental approach we developed is rapid and high throughput, consistent with previous work [26], and should be applicable to many species of resistant bacteria. The mathematical model  we created expresses the problem of antibiotic resistance in general terms, and can therefore be applied to other resistance phenotypes where pleiotropy occurs to identify the antibiotic treatment plans that have the highest probability of reversing the evolution of resistance. The purpose of this study was to determine whether it is possible to use selective pressures to return TEM-genotypes to the wild type state, as observed in 1963 when TEM-1 was originally isolated. The methods may also be used to select for any particular genotype within our data set. As such, we may select with reasonably high probabilities, resistance genotypes that existed at some prior point in time. To highlight this feature, we have named our software package "Time Machine".
Once given growth rates of adjacent genotypes, Time Machine returned treatment plans that restored the wild type state as observed in 1963 with probabilities greater than 0.6 when using the CPM model and greater than 3/8 (>0.375) when using EPM. These results suggest that when possible it is desirable to use actual growth rates rather than rough ranking data.

Substitutions
Drugs associated with substitutions in optimal paths (probability) While these treatment methods may have clinical value, we have yet to determine the ideal duration of each therapy. Additionally the antibiotics included in our study may have different applications in the clinic. A further issue is that if new genotypes arise, the treatment plan may fail. The inclusion of more resistance genes in this type of approach may aid in the creation of robust treatment plans that are effective even when unexpected genotypes arise.
The discrete optimization problem motivated by our goal to reverse resistance, or the challenge to build a better time machine, is of independent mathematical interest. Tables 6 and 7 suggest that the maximum probabilities in each row stagnate after a limited number of steps. This is not always the case. We have created an example (see supplemental information) of two substitution matrices on a 3-locus system where the maximum probabilities can be increased indefinitely (S1 Fig). These results show that great potential exists for remediation of antibiotic resistance through antibiotic treatment plans when pleiotropic fitness costs are known for an appropriate Table 9.

CPM Reversions of Substitutions And Associated β-lactam Antibiotics From Optimal Six
Step Treatment Plans (*Maximum Probability for Path).

Reversions
Drugs associated with substitutions in optimal paths (probability)

Experimental methods
Strains and Cultures. We expressed 16 mutant constructs of the bla TEM gene in plasmid pBR322 from strain DH5-αE. The 16 genotypes differ at all combinations of four amino acid residues and have been previously described [14]. We grew them overnight (16 hours) in standing cultures and diluted them to a concentration of 1.9X10 5 as described elsewhere [14].
We transferred 80 μl of each culture to a 384-well plate with one genotype present in each of the 16 rows. The first 12 wells of each row were antibiotic free (controls) and the last 12 wells contained a single antibiotic at an inhibitory, sublethal concentration. We tested many concentrations and used those that maximized our ability to make comparisons between alleles.
After plating, a membrane is placed over the plate and simultaneously incubated/measured in the Eon Microplate Spectrophotometer at a temperature of 25.1°C for 22 hours. This relatively cool (<37°) temperature is used because degradation of the antibiotics is much slower, while the growth rate of the bacteria is still sufficient to capture the complete exponential period of growth over the duration of the experiment. Overall, we have found that a temperaturẽ 25°C yields more reliable and consistent measurement of growth rates in the presence of antibiotics.
Measurements of cell density (light scattering) at a wavelength of 600 nanometers were automatically collected every 20 minutes after brief agitation to homogenize and oxygenate the culture. Step CPM and EPM Treatment Paths beginning at genotype 1111 and ending at genotype 0000. An arrow indicates that the substitution is included in a path that starts at 1111 and ends at 0000, where the pathway has non-zero probability. Black arrows show substitutions present in six step paths computed using both the CPM and the EPM. Red arrows signify substitutions found only in optimum paths computed using the CPM whereas blue signify substitutions only found using the EPM. Growth Rates. The data obtained from the microplate spectrophotometer is exported to the GrowthRates program to derive the growth rates. In essence, by measuring the optical density at frequent intervals the GrowthRates program can estimate the growth rate,α, through a linear regression algorithm fitting the data from the exponential growth phase. Details can be found in [27] in the section entitled "The Growth Curve" located on pages 233-4. There is not a direct or simple correlation between this method and other methods such as minimum inhibitory (MIC) or disk diffusion testing. The output of this program for the data we collected was a list f(a 1 ),f(a 2 ),. . .,f(a k ) of 15 tensors, each of format 2×2×2×2. These are the rows in Table 4. So if u ∊ {0,1} 4 is a genotype, then f(a i ) u is the fitness of genotype u in the presence of antibiotic a i . This fitness is a growth rate, so we are here using the letter f for a quantity often denoted by α.  One-Way Analysis of Variance (ANOVA) was then used to compare the means of the growth rates obtained, and to determine if there were significant differences between the growth rates of adjacent genotypes.
Correlated Probability Model (CPM): Once the growth rates have been determined under various experimental conditions, the next step is to use them to compute fixation probabilities.
If the (multiplicative) absolute fitnesses W u and W v of two neighboring genotypes u and v, differ by a small quantity then the (additive) relative fitness ln W u W v can be approximated by   Step Antibiotic Cycles. In this figure, cycles are distinguished from paths in that TEM-1 (0000) is the first and last genotype visited, thus creating circular paths. An arrow indicates a substitution included in a mutational pathway which starts and ends at 0000, where the mutational pathway has a non-zero probability for the optimal treatment cycle. The substitutions that are included in optimal two steps cycles are shown in red. Substitutions that are included in optimal four and six step cycles are shown in blue. Four and six step cycles differ only in the number of substitutions and reversions that occur within each cycle. Their probabilities are identical. where T is the generation time. Using a Taylor series approximation, is the probability for v to substitute u, where uj are the neighbors of u with higher fitness than u [23]. Equal Probability Model (EPM): According to the EPM model, the probabilities are equal for all beneficial substitutions, so that one needs the fitness graphs only for computing the probabilities. The matrix entry M(f) u,v is 1/N if genotype u has N outgoing arrows and there is an arrow from u to v.
CPM is accurate if fitness differences between genotypes are small, while EPM may provide better estimates if fitness differences are substantial. Indeed, if the fitness effects of all available beneficial mutants exceed some threshold, then fixation probabilities are independent of fitness values [28]. We applied both CPM and EPM, since no complete theory for substitution probabilities exists. Additionally, comparison of two models is useful in learning how sensitive our results are for variation in substitution probabilities.
Time Machine Programs. Optimal antibiotic sequences and pathways of genotypes: Let M (f(a)) denote the 16×16 transition matrix we derived for the antibiotic labeled a (S1 File EPM Prepare and S2 File CPM Prepare). For any sequence a 1 ,a 2 ,. . .a k of k antibiotics, we consider the matrix product M(f(a 1 ))M(f(a 2 ))M(f(a 3 )). This product is also a 16×16 transition matrix. Its entry in row u and column v is the fixation probability of genotype u mutating to genotype v under the antibiotic sequence a 1 ,a 2 ,. . .,a k . That probability is a sum of products of entries in the individual matrices M(f(a i )), with one sum for each possible pathway of genotypes from u to v. The Time Machine enumerates all 15 k antibiotic sequences of length k, and it selects all sequences that maximize the entry in row u and column v of the matrix product (S3 File EPM Run and S4 File CPM Run). In a subsequent step we then analyze these optimal antibiotic sequences, and for each such sequence, we extract the full list of genotype pathways that contribute (S5 File EPM Out and S6 File CPM Out).
We implemented this algorithm in the computer algebra software Maple, and we ran it for k = 2,3,4,5,6. The running time of the program is slow because of the exponential growth in the number of sequences. At present we do not know whether an efficient algorithm exists for solving our optimization problem for larger values of k.
Cycles of antibiotics: We also used this method to compute cyclical treatment paths in which the starting and ending genotypes were the wild type 0000 (S7 File EPM CyclingOut and S8 File CPM CyclingOut). The problem we solved was somewhat different from the previous one, in that we focused on obtaining the maximal probabilities of a cycle that includes some substitutions and then returns to the wild type without halting. Halting means that adjacent genotypes in a mutational pathway coincide, which is undesirable.
Supporting Information S1 Fig. Locus Model. For any biallelic system and set of drugs, the maximum probabilities for returning to the wild-type depend on how many steps one allows in the treatment plan. The following example demonstrates that the maximum probabilities may increase by the number of steps indefinitely. Consider a three-loci system where the genotypes are ordered as 000; 100; 010; 001; 110; 101; 011; 111: Assume that the starting point is the genotype 100 and that Drugs A and B (see the next page) are available. For the sequence A, the probability for ending at 000 is 0.9, for A-B-A 0.99, for A-B-A-B-A 0.999, and so forth. (TIFF) S1 File. EPM Prepare. File used to convert growth rate averages into data matrices.