Variable Mutation Rates as an Adaptive Strategy in Replicator Populations

For evolving populations of replicators, there is much evidence that the effect of mutations on fitness depends on the degree of adaptation to the selective pressures at play. In optimized populations, most mutations have deleterious effects, such that low mutation rates are favoured. In contrast to this, in populations thriving in changing environments a larger fraction of mutations have beneficial effects, providing the diversity necessary to adapt to new conditions. What is more, non-adapted populations occasionally benefit from an increase in the mutation rate. Therefore, there is no optimal universal value of the mutation rate and species attempt to adjust it to their momentary adaptive needs. In this work we have used stationary populations of RNA molecules evolving in silico to investigate the relationship between the degree of adaptation of an optimized population and the value of the mutation rate promoting maximal adaptation in a short time to a new selective pressure. Our results show that this value can significantly differ from the optimal value at mutation-selection equilibrium, being strongly influenced by the structure of the population when the adaptive process begins. In the short-term, highly optimized populations containing little variability respond better to environmental changes upon an increase of the mutation rate, whereas populations with a lower degree of optimization but higher variability benefit from reducing the mutation rate to adapt rapidly. These findings show a good agreement with the behaviour exhibited by actual organisms that replicate their genomes under broadly different mutation rates.


Introduction
Mutations constitute the main source of genetic diversity in asexual populations. Although most of them have deleterious effects on fitness [1][2][3], natural selection increases the representation of those having beneficial effects, which can become fixed in a population. This combined action of mutation and selection promotes the adaptation to environmental changes, and in the long term leads to the evolution of populations.
In the framework of the fitness landscape described by Wright [4], populations placed near the top of a fitness peak will experience less beneficial mutations than populations placed far from the adaptive optimum [2,[5][6][7][8][9]. This dependence of the mutation effects on the degree of adaptation of populations led to the theoretical prediction that mutation rates would be reduced in constant environments, in which the population has had enough time to adapt. Once the optimum has been attained, a homogeneous population of individuals with the optimal phenotype is the best adaptive solution, so the generation of further diversity is not necessary. Nevertheless, the actual situation is that environments never remain static. They continuously undergo changes that alter the fitness landscapes, displacing populations towards suboptimal fitness regions, where the amount of mutations with positive effects increases. These poorly-adapted populations could benefit from having higher than standard mutation rates.
The variation of mutation effects with fitness, together with the fact that error rates can be easily modified as a consequence of mutations producing genotypes with variable capacity to cause errors, suggest that mutation rates are a character subjected to the action of natural selection [10,11]. Stable environments would favour low mutation rates (anti-mutator genotypes), constrained only by the costs of error-repair mechanisms [12,13]. In contrast to this, environments subjected to frequent changes would select for increased mutation rates (mutator genotypes) that permit faster adaptation to the new conditions [14][15][16]. However, the optimization of the mutation rate is not only determined by its impact on adaptation but also by the consequences that the variation of this character has on fitness. High mutation rates can increase the number of deleterious mutations, whereas low mutation rates can have metabolic costs associated. The existence of these opposing forces causes that natural selection often fails to fully optimize this character [17]. The study of the evolution of mutation rates has been addressed theoretically [13,[18][19][20][21], and using digital organisms [17]. There are also many reported examples of natural and experimental bacterial populations with higher than standard mutation rates [22][23][24][25][26], showing that there are multiple situations in Nature in which being a mutator confers a selective advantage. Although mutator variants have also been isolated in the DNA phage T4, they have been rarely observed in the case of RNA viruses [27,28]. A possible explanation is that RNA viruses replicate their genomes at the maximum error rate compatible with the preservation of genetic information, and additional increases would lead to fitness losses that could cause the extinction of the population [29,30]. Despite this almost absence of mutators in RNA viruses, it has been observed that low fitness clones of an RNA bacteriophage increase their replicative ability when infections take place in the presence of a mutagen, a clear example of the adaptive advantages arising from an increase in the mutation rate in low fitness populations [31].
Anti-mutator mutants, with lower than normal mutation rates, have been observed in bacteria [32], in the phage T4 [33], and in RNA viruses evolving in the presence of mutagens [34]. In the latter case, the anti-mutator phenotype can be produced by single changes in the viral polymerase, without requiring the expression of corrector activities. These observations suggest that RNA viruses could easily evolve to lower mutation rates. If they do not, it could be due to the major adaptive advantages provided by high mutation rates. The finding that high-fidelity genotypes of an RNA virus have lost some of their adaptive properties in mice constitutes a strong support of this hypothesis [35,36]. Other studies, however, point to the existence of a trade-off between rapid replication and fidelity to explain the high mutation rates of RNA viruses [37].
Asexual populations of replicators, such as RNA molecules evolving in silico, with selection acting on their folded conformation, constitute a simple system to study how the variation of the mutation rates influences adaptation. After a sufficiently long time, these virtual populations reach a stationary state characterized by mutation-selection equilibrium and a quasispecies structure [38]. Populations of RNA molecules have been very successfully used as a computational model for the study of evolutionary processes [39][40][41]. The influence of the mutation rate on the degree of adaptation attained at the stationary state, and on the genotypic and phenotypic diversity of the population are questions that have been addressed previously with this model [9,42]. In this work we focus on the adaptability of populations of RNA molecules that reached the stationary state at different error rates, and that are affected by a sudden environmental change. To this end, we determine those mutation rates promoting maximal adaptation after a short number of generations. In practice, our population evolves under selection for folding into a given secondary structure until mutation-selection equilibrium is reached. At that point, it is confronted with a new selective pressure, represented by a new target structure, towards which it evolves under a second mutation rate. Our results show that, before reaching the new equilibrium, and especially at the early stages of the adaptive process, there is no simple relationship between the momentary degree of adaptation and the new error rate at which a population evolves. There is also a strong influence of the mutation rate at which populations evolved towards the previous stationary state in their ability to adapt to new selective pressures. Our results are of relevance to understand the adaptive process in changing environments when variations of the mutation rate are allowed.
Some actual examples of this situation are the in vitro evolution of structural or catalytic RNA molecules and proteins -where the experimenter can manipulate the extension of the genetic diversity generated-the selection of mutator variants of pathogenic bacteria in response to antibiotics, hampering the treatment of many diseases [43,44], and also RNA viruses in which even very mild mutator or antimutator phenotypes can have important consequences in shaping not only virus evolution, but also pathogenesis, transmission, and emergence [28,35,36]. A deeper knowledge of how mutation rates can affect fitness and adaptive ability can be of great importance to evaluate the effectiveness and long term consequences of therapies, especially those based on the increase of the mutation rate through the use of mutagens [45].

Evolutionary algorithm
The system used in our simulations consists of a population of N = 1000 RNA sequences, each of length l = 50 nucleotides. At the beginning of the simulation, each molecule of the population is initialized with a random sequence. Every time that a sequence replicates, each of its nucleotides has a probability (defined by the mutation rate m) to be replaced by another nucleotide, randomly chosen among the four possibilities. We define a target secondary structure which is endowed with the highest replication rate. After each replication event, the molecules are folded into secondary structures with help of the Vienna package, version 1.5 [46] and the base pair distance d i between each molecule i in the population and the target is calculated. The base pair distance is defined as the number of base pairs that have to be opened and closed to transform a given structure into the target structure. The probability p(d i ) that a molecule i in the population replicates depends on the distance d i according to the following equation: where Z is the overall normalization factor Z = The parameter b denotes the selection pressure and takes the value b = 1 for all the simulations carried out in this study; l is the length of the molecule. Generations are non-overlapping and the offspring generation is selected following a Wright-Fisher sampling at each time step. At any time point, the population can be characterized by two main quantities: the fraction of molecules correctly folded (those with d i = 0) and the average distance of the population to the target structure. Because of the dependence of the probability of replication of each molecule on its d i value, we chose the average distance to the target as an estimator of the degree of adaptation of a population (the higher the distance the lower the degree of adaptation).

Stationary, non-adapted, and adapting populations
A general scheme showing the main evolutionary characteristics of the different populations used in our simulations is represented in Fig. 1.
The initial random populations R replicate with different mutation rates (m 1 = 0.001, 0.002, 0.005, 0.01, 0.02, and 0.05) during 8000 generations (g), a number large enough to reach a statistically stationary state, which is determined by the selective pressure represented by the target structure S 1 (the hairpin structure shown in Fig 1). We performed r = 50 independent realizations of this process for each value of the mutation rate. The resulting equilibrated populations (populations E) experience a sudden change in the selective pressure implemented as a change of the target structure from S 1 to S 2 (the hammerhead represented in Fig. 1b). At this point, populations become non-adapted (populations NA), as they have not been optimized in the new environment defined by S 2 . We randomly choose one of the 50 non-adapted populations as the starting point of a new adaptive process towards S 2 . The population evolves during 200 generations under 6 different values of the mutation rate (m 2 = 0.001, 0.002, 0.005, 0.01, 0.02, and 0.05) giving rise to the adapting populations (populations A) shown in Fig. 1. We performed r = 1000 independent realizations for each value of the mutation rate. The average distance was determined for each adapting population at every generation, and averaged over the 1000 independent runs.

Nomenclature
We focus our analyses in the three populations E, NA, and A. In the following, we will call d i,j E (S 1 ) the set of base pair distances between each molecule i in populations E and the target structure S 1 for each realization j. Likewise, we call d i,j NA (S 2 ) the set of base pair distances between each of the molecules in populations NA and S 2 for each of the realizations. It is important to remark that populations E and NA have the same composition. The only difference between them is the structure used to calculate the d i,j values, which in populations E has been the target of a completed optimization process whereas in populations NA has not. Analogously, we define d i,j A (S 2 ) as the base pair distance between molecules in the adapting populations A and the target structure S 2 for each realization j. Note that the latter set of distances is a timedependent quantity, while d i,j E (S 1 ) correspond to populations at the mutation-selection equilibrium, and are thus independent of time. Average values are calculated over the N = 1000 molecules in the population and over r independent realizations of the process. We first define the average distance to target for each realization j as, vd E j (S 1 )w. As indicated, this is an important quantity used as a measure of the degree of adaptation of populations. Similar quantities are defined for nonadapted NA and adapting populations A. The number of realizations is r = 50 for E and NA populations and r = 1000 for A populations. For adapting populations A we calculate the standard deviation s A (S 2 ) of the average distances vd A j (S 2 )w with respect to d A (S 2 ) in order to measure the variability among realizations. The statistical significance of the differences between average distance values was determined with the Student's t-test.
To get a deeper characterization of populations E and NA, we also determined several additional statistical parameters. We call SD E j (S 1 ) the standard deviation of the distances d E i,j (S 1 ) with respect to vd E j (S 1 )w. The average over standard deviations for the r realizations is SD E (S 1 )~r {1 P r j~1 SD E j (S 1 ). The minimum and maximum distances in each set d E i,j (S 1 ) are called d E Min,j (S 1 ) and d E Max,j (S 1 ), respectively. The averages over realizations are The skewness of the set d E i,j (S 1 ) is called Sk E j (S 1 ), and the average over realizations is Sk Sk E j (S 1 ). Analogous quantities are defined for NA populations. All statistical calculations were performed with the program Mathematica 5.0 (Wolfram Research).

Distributions of distance values in equilibrated and nonadapted populations
The structure of the populations used in this study can be described through the distribution of distance values (see Methods), which depends on the value of the mutation rate, the selective pressure dominant, and the time of evolution. The statistical parameters characterizing the distribution of distances to targets S 1 and S 2 in equilibrated and non-adapted populations are reported in Tables 1 and 2, respectively. Given the population structure at the stationary state, d E i,j (S 1 ), and the new selective pressure, S 2 , the set d NA i,j (S 2 ) is completely determined, illustrating how the composition of populations at equilibrium influences its subsequent evolution.
Equilibrated populations are optimized at a mutation rate m 1 towards the target structure S 1 (Fig. 1). At the stationary state, the average distance of these populations to the target structure S 1 depends on m 1 (see values for d E (S 1 ) in Table 1). The higher m 1 , the larger is the average distance to the target at the mutationselection equilibrium. This result illustrates the dependence of the degree of adaptation on the mutation rate, as previously reported [42]. When the selective pressure is changed by choosing a new target structure (S 2 , see Fig. 1), populations optimized to reach target structure S 1 are non-adapted when confronted to S 2 . As a rule, d d NA (S 2 ) values in non-adapted populations (see Table 2) are   1).  Table 1), although the two populations differ much more at low than at high error rates. Actually, if the target structures S 1 and S 2 are at a large base pair distance of each other, optimized populations are typically farther from S 2 than a population of randomly chosen sequences, especially for small m 1 (results not shown).
In addition to the the average distance, we have evaluated the average standard deviations [SD E (S 1 ) and SD NA (S 2 )] as a measure of the relative phenotypic diversity, the average minima

Adaptive dynamics of stationary-state populations
We have analyzed the adaptive dynamics of the stationary-state populations described in Table 1 when the selective pressure is changed by choosing a new target structure (S 2 , see Fig. 1), as described. Populations optimized to reach target structure S 1 are non-adapted when confronted to S 2 , and constitute the initial condition (g = 0) of a new adaptive process (see Table 2). With the aim of determining the optimal mutation rates that promote adaptation to S 2 , each population was allowed to replicate under a range of values of the new error rate (m 2 ) between 0.001 and 0.05. The variation of the value of d d A (S 2 ) was evaluated through 200 generations. The results obtained for three representative populations differing in the value of m 1 at which they reached the previous stationary state are shown in Fig. 2. In the three cases considered we observe large differences in the adaptive dynamics depending on the value of m 2 used to adapt to the new selective pressure. For the largest value of m 2 considered (0.05), d d A (S 2 ) decreases only slightly as the number of generations increases, showing that too high error rates strongly hinder adaptation. For m 2 values lower than 0.05, there is a noticeable decrease in the average distance as time elapses. In general, the variation of d d A (S 2 ) is faster at the beginning of the adaptive process, slowing down later. The initial decay occurs more rapidly for populations that reached the previous stationary state at moderate to high mutation rates than for those that were previously optimized at low mutation rates.
An interesting result observed in the evolution of the three populations represented in Fig. 2 is that, during the transient, the value of m 2 that promotes the maximum degree of adaptation to S 2 (or the minimum value of d d A (S 2 ) at a given time) depends on the number of generations elapsed. There are also differences among the three populations depending on the value of m 1 at which they were optimized towards S 1 . We can draw the general conclusion that the optimal mutation rate promoting maximal adaptation to a new selective pressure before reaching the new stationary state depends on the number of generations elapsed under the new conditions and on the previous mutation rate at which the population had evolved. In other words, the evolutionary history of populations has important consequences in their posterior adaptive capacity.

Previous state of the population and optimal mutation rates at early stages of adaptation
To explore the relationships between the value of m 1 at which populations adapted to S 1 and the value of m 2 that promotes maximal adaptation to S 2 after a short number of generations, we evaluated d d A (S 2 ) after 40 generations in populations that differed in m 1 and evolved at increasing values of m 2 (Fig. 3). We have chosen g = 40 to observe populations that are still far enough from equilibrium but have undergone a significant degree of adaptation. Other values of g meeting these requirements would be appropriate as well and do not qualitatively change our results.
In all populations the variation of d d A (S 2 ) as a function of m 2 shows a non-monotonic behaviour, with important differences depending on the value of m 1 at which they reached the previous stationary state.  Another observation emerging from the results of Fig. 3 is that the influence of m 1 is maximal at low values of m 2 . As this parameter increases, populations behave more similarly, and at m 2 = 0.1 all of them converge to approximately the same value of d d A (S 2 ) (results not shown). That is, the lower the value of the mutation rate m 2 , the higher the influence of the previous state of the population.
The lowest values of d d A (S 2 ) reached in all the cases after 40 generations, and the value of m 2 at which it is obtained, are shown in Table 3. We observe that the maximal degree of adaptation appears in two quite different situations. The first one corresponds to a population that reached the previous stationary state at  Table 1). That population, when confronted with the new selective pressure S 2 , displays the almost lowest value of d d A (S 2 ) obtained in our simulations at g = 40 for all m 2 ,0.01. Interestingly, the population that performs worst for any value of m 2 is that with the lowest value of m 1 (Table 1), which permitted the highest degree of adaptation to S 1 . As indicated above, a small increase in the value of m 1 (from 0.001 to 0.002) yields a considerable increase in the adaptive capacity of this population (p%0.001 for the difference of d d A (S 2 ) between populations evolved at the same m 2 and optimized at m 1 = 0.001 or 0.002), highlighting once more the significative nonproportional effects that small changes in the mutation rate might produce.

Discussion
The observation that mutation rates per nucleotide vary by orders of magnitude across species suggests that this character has not an optimal universal value [10,47,48]. Each species  evolves under a mutation rate arising from many factors that are not universal. Among the most relevant we find the variability of the environment, the effect that mutations have on fitness, the metabolic costs of having more faithful replication machinery, the population size, and the replication rate [49]. Variations in any of these factors can modify the optimal value of the mutation rate, with the result that natural selection has to carry out a continuous fine tuning of this parameter and under the action of non-compatible trends may fail to find an optimal solution [11,17].
In this paper we have focused on the study of the values of the mutation rates that promote maximal adaptation in a short number of generations when populations previously optimized at different error rates experience a single environmental shift. The model system we have used is constituted by ensembles of RNA molecules that evolve through mutation and selection towards a defined target structure. This system permits to establish direct correspondences between the genotype (the sequence of the RNA molecule), the phenotype (the structure into which it folds), the replicative ability (inversely related to the distance to the target), and the degree of adaptation of the whole population. Although in our simulations, mutation rates are imposed by the researcher, the relative amount of beneficial and deleterious mutations is not a fixed parameter, as it varies through the evolutionary process as a consequence of the variation in the degree of adaptation [9]. Consequently, the optimal mutation rate at the stationary state does not necessarely coincide with the optimal mutation rate before mutation-selection equilibrium is reached.
The populations considered in this paper have first evolved with different error rates (between 0.001 and 0.05) for a number of generations, large enough to reach mutation-selection equilibrium. Their degree of adaptation at equilibrium decreases with the mutation rate. However, when confronted with a new selective pressure, these populations can experience adaptive advantages if they vary their mutation rate. This means that the optimal mutation rate for an adapting population can be quite different from the optimal mutation rate under conditions that remain constant for a long time. In particular, for populations optimized at low mutation rates an increase of this parameter may be favourable, while for populations replicating under high mutation rates a decrease would be advantageous. These results can be partially explained as a consequence of the different dynamics of the adaptive process at different mutation rates, and by the influence of the composition of the population in its subsequent ability to adapt to new selective pressures. Adaptation is a complex phenomenon in which, in addition to the diversity generated de novo, the nature and distribution of existing mutants plays an important role. In our simulations, populations able to attain a high degree of adaptation to a new selective pressure in a short time were those previously optimized at moderate (0.002) to high (0.05) values of m 1 (Fig. 3 and Table 3). After changing the target structure, populations optimized at high mutation rates give rise to highly diverse, non-adapted populations, which contain in their mutant distributions molecules closer to S 2 at g = 0 than populations optimized at lower error rates. These populations respond better when the mutation rate is decreased, thus enhancing the presence of structures close to S 2 . On the other hand, populations optimized at m 1 = 0.002 benefit from an increase in the mutation rate to rapidly adapt to the new target structure. At g = 0, these populations display lower diversity, and the molecules closest to S 2 have typical distances larger than those in populations optimized at high error rates. Increases of the error rate promote the appearance of better adapted structures, the substrate of further optimization. These results provide good examples of the importance of pre-existent and newly generated diversity in the adaptive capacity of populations, and illustrate how the relative amount of each of them determines whether a population needs to increase or reduce the mutation rate to get rapid adaptation.
We have used an exponential fitness function, such that replication of molecules with any value of d i is possible. This is one reason why we do not observe extinction in any of the cases studied. The situation would be different should we work with a truncated landscape, for instance, where molecules folding too far from the target structure would not have the minimal functionality required for replication. More restrictive landscapes of this kind, together with a population that could vary its size, would yield extinction in cases of low m 1 , especially. Still, the conditions for extinction in truncated landscapes would be alleviated in large enough populations, with a broader range of d i values and hence of diversity.
Our results can be discussed in the context of actual organisms that replicate using error rates that differ by orders of magnitude (from 10 28 -10 29 in DNA organisms to 10 24 -10 26 in RNA viruses). Although the error rates for the RNA molecules used in this study are apparently much higher than the error rates of actual organisms, they become much more similar if the values are expressed per genome, or replicating molecule. Taking genome size into account, our RNA molecules replicate with error rates between 0.05 and 2.5 errors per molecule and generation, a value quite similar to the error rates found in Nature, which range between roughly 0.003 errors per genome and generation in DNA based microbes to 1-5 in RNA viruses [10]. These values can be further increased by mutagenic agents. Populations optimized at low mutation rates present a high degree of adaptation and a low phenotypic diversity (see Table 1). After a limited number of generations, under a new selective pressure, the highest degree of adaptation reached by these populations takes place when they increase the mutation rate to values of 0.01 or 0.02 ( Fig. 3 and Table 3). These populations behave similarly to DNA organisms that maintain a certain degree of constancy in the intracellular medium, which allows a reduced mutation rate for replication. However, under conditions of environmental stress, a more convenient evolutionary strategy for these organisms would be to increase the mutation rate to get rapid adaptation [14][15][16]. This strategy is particularly important when the new selective pressure is strong enough to extinguish the population in case it is not able to adapt in a short time span. A well known example is the selection of hypermutator strains when bacterial populations are infected with phages or exposed to antibiotics [43,44,50]. There is however a limit for the beneficial effects derived from the increase of the error rate. In our system, when the mutation rates are increased above m 2 ,0.01, the degree of adaptation reached diminishes (see Fig. 3), a result that suggests that there must be also a limit for the increase of the mutation rate that can be attained by a hypermutator strain. In our case, this limit is set by the mutation rate above which fixation of the optimal phenotype becomes impossible.
There are important differences in the behaviour of populations evolving at low error rates. In our system, the population that performs worst under a new selective pressure at any of the values of m 2 assayed is that one previously optimized at the lowest mutation rate considered (m 1 = 0.001; see Fig. 3). This population reached the highest degree of adaptation at the stationary state (see Table 1) and resembles specialist organisms that perform optimally under a narrow range of very well established conditions, but have difficulties to adapt when these conditions are modified [51]. Populations optimized at a slightly higher error rate (m 1 = 0.002) experience a substantial increase in their adaptive ability ( Fig. 3; Table 3). This population could represent less specialized organisms, not so well adapted to a concrete selective pressure, but able to perform optimally in a wider diversity of environments [52] due to their higher degree of genotypic and phenotypic diversity. These results also suggest that the transition from a rigid population, with little ability to respond to environmental changes, to a more flexible population, able to adapt rapidly, may occur through small increases of the mutation rate that produce concomitant increases in the pre-existent diversity.
Populations optimized at high mutation rates (m 1 = 0.01) can adapt rapidly to a new selective pressure keeping this parameter constant. These cases resemble in some aspects RNA virus populations, which, together with viroids, replicate at the highest error rates found in Nature and that are able to adapt rapidly to changing environments without significantly altering the error rate. The disadvantages experienced by these populations when the error rate is further increased [30,53] agree with the fact that only mild mutator mutants have been isolated in the case of RNA viruses [27,28]. In normal conditions, RNA viruses also do not select mutants with lower than standard mutation rate. It is likely that the associated reduction in the genotypic and phenotypic diversity strongly challenges their ability to undergo adaptation [35,36].
As could have been expected, populations optimized at values of m 1 $0.02 (above the fixation threshold) can dramatically increase their adaptive potential if they reduce the mutation rate. These populations sustain a high degree of diversity, with molecules already close to other possible new target structures. If the mutation rate is kept high, molecules with low distance values are lost due to the strength of mutation. However, when the mutation rate is reduced, these molecules can become fixed, permitting in this way the optimization of the whole population. The equivalent in Nature of these populations could be RNA viruses replicating at higher than standard mutation rates, a condition that has been explored experimentally by exposing RNA virus populations to mutagens. RNA viruses can escape the negative consequences of the increase of the mutation rate by selecting anti-mutator mutants [34]. The consequences of the selection of an anti-mutator mutant in a population previously mutagenized have not been explored, but they could be quite negative from the viewpoint of the host if, as predicted by our model, they are associated to a short term increase in the adaptive ability. The treatment of RNA virus infections with mutagens is being investigated as a new therapeutic approach known as lethal mutagenesis [53]. One of the main criticisms to this alternative therapy is that the increase of diversity caused by the mutagen could induce a parallel increase in the adaptive ability of the virus [54]. Our results suggest that if the mutagen is withdrawn before infection clearance, or if antimutator mutants emerge, the resulting populations could adapt more easily to new selective pressures. A great care should be taken when manipulating the error rate of pathogenic organisms. The associated variations in fitness and adaptive capacity could result in the generation of strains better suited to resist new treatments or the action of the immune system of the host. Therefore, a continuous research in this field combining both experimental and computational approaches is highly promising.