In the long term steady state, the solutions of this dynamical system are obtained by imposing the conditions on pone.0052742.e067pone.0052742.e068Eqs. (1), solving the ensuing algebraic equations, and checking for the stability of the solutions, by looking at the eigenvalues of the Jacobian matrix, evaluated at the respective solution. The solutions obtained this way in the general case turn out to be quite complex, so in order to simplify the resulting expressions, we choose particular values of the fitnesses, namely , and , respecting their natural ordering. Solutions for other values can be obtained using the same steps.

. Densities of (blue), (green) and (gray) genomes as a function of the environmental mutation rate for fixed , and as a function of for fixed . The left panel corresponds to , and the right panel to . Dashed orange lines represent the average fitness of the population.

The analytical solutions given by Eqs. (3) and (4) are quite complex, and it is difficult to extract direct interpretations from them. However, the behavior of the solutions can be understood in the biologically relevant region of small . Fig. 2 shows the densities as a function of or , at fixed and respectively, for two values of , along with the value of the total fitness . In the upper left corner plot for each value of , we consider the case , representative of the biologically relevant scenario of small mutation rate. When is very small, most of the population stays aligned with the environmental feature and . As increases, maladapted genomes appear and eventually overcome the specialized genes. At a definite value of , however, a discontinuity takes place and neutral individuals suddenly appear and become the majority of the population, while the density of both maladapted and specialists decreases. The decrease in specialists is larger for larger values of . Thus, for sufficiently large and small , trying to catch up with the rapidly evolving environmental feature is not a viable strategy, since the risk of producing a maladapted offspring becomes destructive. Interestingly the strategy adopted by the majority of the individuals guarantees the maximum average fitness in any given region of the plane, for every value of . For large values of (lower left plots), the situation is qualitatively different. For small , maladapted and neutral individuals are almost equally numerous. When increases beyond a threshold, neutral individuals again appear suddenly, but they are unable to overcome maladapted genomes. Only for large values of are neutrals capable to prevail over the specialists. The right plots for each value of in Fig. 2 show the evolution of the species' densities as a function of for fixed . In this case, for small neutrals are absent from the system, and there is a simple competition between specialists and maladapted, the former being predominant for small genetic mutation rates, but going extinct for large . On the other hand, when the rate of environmental change is sufficiently large, we enter a new scenario in which neutrals are predominant for small . Beyond a mutation rate threshold, however, neutrals suddenly become extinct, and their population is replaced by maladapted genomes, while specialists decrease their density for large . Interestingly, in this region of large , specialists can survive even for very large mutation rates, close to , due to the effect of a nonzero that prevents their complete elimination. Fig. 3 shows the complete picture of the relative species' abundances as as a function of and for the previously considered values of .

Solution describes a density that is negative in the parameter region , i.e., it is an “unphysical” solution that does not describe any realistic scenario. Solution has nonzero densities in the whole parameter space, while solution is physical only in the region (2)In order to find the relative stability of the physical solutions and , we consider the Jacobian matrix of the equation system pone.0052742.e067 pone.0052742.e068 Eqs. (1) , taking the form We then compute the eigenvalues of matrix , evaluated for the different solutions and . In any given region of parameter space, the stable solution (in the stationary limit) is the one possessing a negative largest eigenvalue. Examination of these eigenvalues, operation performed again with the help of standard computational software packages, leads to the solution:

Case

A more precise mathematical characterization of the phenomenology discussed above, and in particular of position of the transitions taking place for different values of and , can be obtained in the particular case . Here, qualitative arguments allow us to solve the model in a much simpler way, for general values of , and . This analysis, moreover, reveals the role of in the dynamics of our model.

The relevant equations in the case read (5a) (5b) (5c)To find their solution, we argue as follows: If is very close to , all terms in square brakets in pone.0052742.e162pone.0052742.e163Eqs. (5) will be negative. Therefore, the only stable solution will be , , for any value of . Under this conditions, the quantities within square brackets in Eqs. (5a) and (5b) take the form and , respectively. Decreasing the value of , the first solution with will take place for the first of these values that become zero. This occurs when is smaller than either or , respectively. Since , we have , and this transition will always be physical. However, for , we have and it is not physical. In this case, when , if , decays exponentially, and in the long time limit ; the existence of a non-zero solution imposes , or , from where the restricted normalization condition leads to the solution , and . In the case , which density or becomes first non-zero depends on which threshold, or is larger. Thus, if , then . Therefore, when decreasing , the first density to take a non-zero value is . decays again exponentially, so the solution is the same as in the case . Finally, for , is the first density to become non-zero when . The steady state solution comes from imposing , leading to . In this region, the factor in square brackets in Eq. (5b) becomes negative, indicating an exponential decay and a corresponding steady state value . We are therefore led to the solution, using the normalization condition, , .

The final solution in this case can thus be summarized as follows:

For If (6) If (7)

For If (8) If (9)



Fig. 4 sketches the phase diagram, as a function of and , resulting from the previous equations. The different scenarios for small and large values of are now explicit. For small , specialist individuals (in the class) are able to survive, and even dominate the population, as long as the mutation rate is small. In fact, for , the density of specialists is larger than the density of maladapted individuals. For larger , the density of specialized genomes decreases, until it reaches the -dependent threshold , leading to a continuous, second order, phase transition (akin to the error catastrophe in quasispecies models [28]) beyond which the whole population becomes maladapted and thus prone to eventual extinction. In all of this region of small , neutral individuals are irrelevant. For small , specialists perform much better, while for large only maladapted individuals survive.

When increases, a different picture emerges. For fixed, small , the explicit behavior of as a function of , reads we can obtain from Eq. (6) (10)for , and zero otherwise. Thus, when crossing , the density of specialized individuals experiences a first order transition to extinction, with a jump of magnitude (11)The sudden extinction of specialized individuals coincides with the abrupt emergence of neutrals, in a related first order transition for with an associated jump (12)In this large region, neutrals are able to cope with environmental change if the mutation rate is sufficiently small, again up to a maximum mutation rate , after which continuously becomes zero and only maladapted individuals can survive. These discontinuous transitions as a function of at fixed are reminiscent of the phenomenology observed in quasispecies models with higher order replication mechanisms [35]. We note, however, that transition as a function of at fixed are all continuous.

Fig. 5 shows the proportions of the three genomes along with the average fitness as function of for fixed , and as a function of at fixed (left panel), and the general scenario as a function of both and (right panel). As it is clear, while an abrupt transition occurs at the level of the genome frequencies at ( in Fig. 5), the average fitness exhibits a continuous behavior, decreasing monotonously from the maximum for to for large values of (when all the other genomes simply mutate into ). When is fixed (left column, left panel), increasing causes an increase of genomes and a simultaneous decrease in and . As , however genomes disappear as the neutral genomes abruptly appear. The latter guarantees a constant value of . genomes are constantly created due to the genetic mutation rate, but their fitness is lower so they do not reproduce frequently. This scenario is stable, and any further increase in does not produce any effect. The role of is better understood at fixed values of . When (top panel) increasing deteriorates the fitness of the population since genomes are substituted by ones, which eventually become fixed ( as in figure). A similar behavior is observed for (bottom panel), but here the genomes take the place of the genomes, till the latter disappear at for the values of the simulations.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 5. Species densities at the steady state for small . Left panel: Genome densities as a function of the environmental mutation rate for fixed , and as a function of for fixed . Fitness values are chosen as , and . Hence, , and , implying (see main text). Right panel: Densities of (blue), (green) and (gray) genomes as a function of the environmental mutation rate and the genetic mutation rate . https://doi.org/10.1371/journal.pone.0052742.g005

The crucial difference between the cases and , as can be observed from the comparison of Figs. 2 and 5 is the effect of a positive density of specialists for large and . In the case , the density goes to zero after the corresponding transition, specialist being unable to cope with extreme genetic and/or environmental rates of change. In the presence of a non-zero , signaling the possibility of collateral beneficial effects of an environmental change to previously maladapted individuals, susceptible individuals are still able to thrive in an situation combining both fast genetic and environmental change (see lower right plots in Fig. 2). This effect is due to the feedback mechanism induced by the parameter , that allows the replenishment of the individuals from previously maladapted individuals. Their prevalence is however relatively small, and comparatively negligible with respect to the predominant species, either or , especially in the case of small populations. Finally, it is worth noting that, while the prevalence of genomes is stable in our model, it can be interpreted as a metastable state leading to extinction in a multi-species scenario.