Evolution of joint cooperation under phenotypic variations

Effects of phenotypic variation on the species-environment systems and the evolution of cooperation under prescribed phenotypic diversity have been well addressed respectively. Interspecies interactions in the context of evolvable phenotypic diversity remain largely unconsidered. We address the evolutionary dynamics by considering evolvable phenotypic variations under group interactions. Each individual carries a capacitor of phenotypes and pays a cost proportional to its volume. A random phenotype from the capacitor is expressed and the population is thus divided into subpopulations. Group interactions happen in each of these subpopulations, respectively. Competition is global. Results show that phenotypic diversity coevolves with cooperation under a wide range of conditions and that tradeoff between expanding capacitor and rising cost leads to an optimal level of phenotypic diversity best promoting cooperation. We also find that evolved high levels of phenotypic diversity can occasionally collapse due to the invasion of defector mutants, suggesting that cooperation and phenotypic diversity can mutually reinforce each other.

To understand the emergence and persistence of cooperative behavior among selfish individuals, a variety of mechanisms have been proposed. One intensively studied is the tag based cooperation. Riolo et al. have originally constructed a model with continuous tags 1 , but does not count the full defection as a possible strategy 2 . In ref. 3 , Traulsen and Schuster have simplified and presented the basic aspect of, the original model 1 . Later they extended this simplified model 3 to the spatial settings 4 . Jansen and van Baalen 5 have explored the chromodynamics on lattice-structured populations. One gene encodes for tag and another for strategic behavior. It has been shown that tight coupling of strategy and tag inhibits tag diversity, leading to the low level of altruism. Loosening the coupling to some degree but not too much induces the co-existence of many tags and thus upgrading altruism to high levels. In another work, Traulsen and Nowak have considered the chromodynamics in well-mixed populations of finite size 6 .
These studies have involved two key aspects, tags (or phenotypes) and contingent cooperation [1][2][3][4][5][6][7][8][9][10] . Phenotypes serving as signals are observable and thus individuals are able to base their strategic behaviors on observed phenotypes of opponents. In these models 1-10 , a set of phenotypes are available. Once the set is pre-assigned, it does not evolve over time and thus the phenotypic diversity is fixed. Individuals have an equal access to each of these phenotypes. Natural questions arise which levels of phenotypic diversity would be selected if the phenotypic diversity itself evolves, and if cooperation can get stabilized under these selected levels?
On the other hand, recent experimental studies [11][12][13][14] have shown that genetically identical cells are capable of exhibiting different phenotypes, and a single cell can switch between distinct phenotypes [15][16][17][18][19] . By phenotype switching, the population can either survive fluctuating environments 16,18 , or maximize fitness 19 , or preserve some properties 15,17 . These observations have at least three implications: individuals carry redundant phenotypes; individuals have the ability to switch phenotypes; some phenotypes are more fit than others in some environment while it turns around when environment changes. Moreover, these studies 11,12,[15][16][17][18][19] have well captured the species-environment interactions but seldom dealt with direct inter-species interactions.
In fact, when cooperators cooperate probabilistically, the cooperative act itself seen as phenotype exhibits variations 20,21 . In ref. 20 phenotypic variations differentiate the cooperation gene carriers into those who have expressed (i.e., having cooperated) the gene and those who have not. The later can benefit from their genetically identical species's contributions as defectors do. When the population is structured by an infinite number of demes, this phenotype-mediated interactions suffice to evolve cooperation. This claim is borne out experimentally by regarding TTSS-1 expression as phenotypic trait. Another experimental research 21 has also confirmed that cooperative virulence can be stabilized by the expression of an avirulent phenotype. Even cheaters cheat facultatively in the social amoeba Dictyostelium discoideum 22 . We would like to emphasize that in these studies phenotype describes strategic behavior but does not serve as signal. Interestingly, the authors 23 have pointed out that under greenbeard mechanism, an allele can play dual roles of recognition and thus preferentially directing the benefits to copies of itself in others.
We shall construct a model incorporating the evolvable phenotypic variations and strategic behaviors to address the resulting coevolutionary dynamics. We shall consider a well-mixed population of finite size (=N). Each individual is well defined by three properties: a capacitor of potentially expressible phenotypes, phenotype expressed, and strategic behavior. Capacitor prescribes which phenotypes can be expressed. Each individual randomly expresses one phenotype from her capacitor. This phenotype expressed serves as signal and is observable. Strategic behavior defines whether or not to cooperate. Each individual needs to bear a cost to retain capacitor and express phenotype. According to the phenotype expressed, the population is divided into a number of subpopulations, each consisting of individuals expressing the same phenotype. Cooperators therefore can tunnel their help towards those of the same phenotype. This is perfectly captured by limiting public goods interactions to each subpopulations (For more details see Model description and Fig. 1). Unlike previous models, our model does not preassign the level of phenotypic diversity. Rather we let it be an evolvable trait. This setting allows us to probe and answer questions previously asked. In the limit of small mutation rates, the population dynamics can be approximated by an embedded Markov chain. Our results show that cooperation and phenotypic diversity can coevolve within a wide range of conditions, and moreover, natural selection favors an optimum level of phenotypic diversity.

Results
Let us start with the pairwise invasion dynamics, which shall be profitable for understanding the full population dynamics. When the mutant defectors attempt to invade the resident defectors, they incur no cost and generate no public goods in terms of the game interaction, but still burden the cost of retaining capacitor and expressing phenotypes. Larger volume of capacitor leads to a higher cost. Thus those defectors who carry more potentially expressible phenotypes put themselves in a more unfavorable position. This explains that the plot peaks at the bottom right corner and cascades downward towards the top left corner (Fig. 2a).
In fact, this property can be rigorously verified. Suppose that the invading defectors carry capacitor of volume K Y and the invaded defectors carry capacitor of volume K X . Irrespective of the population composition and their actual phenotype expressions, the fitness is βθ − e K Y and βθ − e K X , for an invader and a resident, respectively. It can be easily obtained that the fixation probabilities read as ρ ) . Using Equation (1) as presented in Methods Section, we can get the transition rate from state X to state Y as r X Y K K ( , ; , ) for both K X > K Y and K X ≤ K Y . This transition rate decreases with K Y and increases with K X , respectively. Similar trends are observed when the mutant cooperators compete with the resident defectors. Obviously, before cooperators increase to the minimal number as required for the public goods game to happen, these cooperators fare like defectors. Due to neutral drift, once cooperators grow to the extent that they can play games, they Figure 1. Phenotypic diversity and local interactions. As phenotypic variation is considered individuals may vary in capacity of carrying phenotypes. Each individual expresses one phenotype at random from those he has carried. Only phenotypes expressed are visible. According to the phenotype expressed, the population is divided into a certain number of subpopulations, each consisting of individuals expressing the same phenotype. In this schematic illustration, there are 11 individuals, say G1, G2, …, G10 and G11. They are located in three subpopulations based on phenotypes expressed: Subpopulations Red, Blue and Yellow. For instance, individual G1 just carries one phenotype Red and expresses this phenotype. Individual G11 carries three potentially expressible phenotypes, say Red, Blue,and Orange and happens to express Red. So G1 and G11 belongs to the same subpolulation Red and have chance to play the public goods games. Following the same logic, one can grasp other individuals in term of phenotype carrying and expressing. Obviously, individuals carrying more phenotypes are more likely to enter each of these subpopulations. But as each individual needs to pay a cost proportional to the number of phenotypes carried, those individuals who carry too many phenotypes would be selected against. When the size of the subpopulation is less than the group size as required for the public goods game to happen, no interaction takes place in this subpopulation. At this time, fitness is uniquely determined by the cost of retaining capacitor and expressing phenotype.
can accrue payoffs. These payoffs compensate the cost of retaining the capacitor and expressing phenotypes and thus give them an upper hand in competing with defectors. As a result, cooperators share a higher chance when attempting to invade defectors than defectors do (Fig. 2b).
Interesting scenarios appear when the mutant defectors compete with the resident cooperators. As is well established, natural selection favors defection over cooperation in well-mixed populations. Inherent phenotypic diversity offers cooperators with potential opportunities to escape the chase and exploitation of defectors. To do this, cooperators need to expand their capacitors, which undoubtedly raises the cost of carrying increasing phenotypes. So how many phenotypes should cooperators carry? Fig. 2c shows that there exists an optimal interval in terms of capacitor's volume. The interval is optimal in two senses that mutual breed between cooperators are able to counteract additional cost due to carrying more phenotypes, and that cooperators can effectively ward off direct competition with defectors. On the whole, cooperators possessing such capacitors enjoy the strongest resistance against defectors' invasion, represented by the dependence of the transition rate on K C in a U-shaped way. Meanwhile, the transition rate increases on either side of this interval. Reasons differ. At the left-hand side, cooperators carry fewer phenotypes, which can be quickly 'deciphered' and thus invaded. At the right-hand side, rising cost of expanding capacitor destroys the survival of these cooperators themselves.
The invasion dynamics shows similar sensitivity on K C when the mutant cooperators attempt to invade the resident cooperators. Before the invading cooperators are likely to play the public goods games, they compete with the resident cooperators as if they were defectors. Once they expand to the group size, the evolutionary race goes another way. From then on, these invading cooperators can also benefit from reciprocity. Especially when the resident cooperators possess large-volume capacitors, the competitive advantage of the invaded cooperators over the invading cooperators is not significant or even non-existential, thus offering the invading cooperators much time of space to expand. Therefore at this time as invaders, cooperators can invade the resident cooperators more easily than defectors do (Fig. 2d).
With the elaborate elucidation of the pairwise invasion dynamics, we are now able to analyze the full population dynamics. Consider the competition of an arbitrary number of strains (=2M) in the well-mixed population of finite size. In the limit of rare mutations, the transition rates between any two strains have been analytically derived (see Equation (1) in Methods Section), defining an embedded Markov chain. By the chain we can compute the stationary distribution of the population, corresponding to the fraction of the time that the population spends in each of these 2M homogeneous states in the long run. We have corroborated by numerical simulations that our main results show no qualitative change for a wide range of mutation rate. In particular, the overall cooperation level is 0.5024, 0.6992, 0.8107 corresponding to the mutation rate 0.5, 0.05, and 0.01 respectively, as other parameters are kept the same as in Fig. 3c.
Results are presented in Fig. 3, where we plot the equilibrium level of these 2M strains as the parameter θ varies. Except for θ = 0, the distribution reveals general features. For defective strains, the larger capacitor they possess, the lower their fraction in the long. In other words, fractions monotonically decrease with expanding capacitor, a feature invariable with the change of θ. When it comes to cooperative strains, there exists an optimal solution in terms of phenotypic variations, associated with the capacitor's volume K C whose corresponding strain accounts for the highest fraction of all cooperative strains. Cooperative strains next to K C also enjoy reletively high fractions. Moving away from K C towards either side, the fraction of cooperative strain ramps down. This optimal volume lies in between one and M, the capacitor's upper boundary. It is subject to θ. Growth in θ reduces the fraction of each type of cooperative strain and this optimal volume. Of interest, as θ expands to as high as 0.3 (Fig. 3d), defective strain with capacitor just covering one phenotype, especially outperforming the cooperative strain with capacitor covering K C phenotypes, attains the highest fraction of all strains. Even so, the overall cooperation level can be still higher than the overall defection level. Therefore, the cooperation can be preserved and selected by diversifying inherent phenotypes.
Intuitions responsible for the emergence of an optimal level of phenotypic diversity are interesting and fundamental. As mutation is uniform and unbiased, it is the evolutionary force that uniquely determines the eventual fate of strains. Resident defective strains possessing capacitors of large volume can be easily invaded by defective strains, or by cooperative strains, when both just possess capacitor covering a small number of phenotypes in the evolutionary process. Occasionally it may happen that defective strain attempts to upgrade its capacity in expressing phenotypes, but as soon as it does so it will be pulled back. As a consequence, defective strains most of the time are stuck in low levels with respect to phenotypic variations. This undoubtedly limits the ability of defective strains in chasing and exploiting cooperative strains.
This, however, does not mean that the larger capacitors cooperative strains possess, the more they will be favored by selection. In fact, when carrying large numbers of phenotypes is prohibitively expensive, though cooperative strains are able to escape stalking of defective strains, mutual breed between themselves does not suffice to offset the cost of carrying so many phenotypes. These cooperative strains will be in a disadvantage place. When cooperative strains possess capacitors covering a small number of phenotypes, cost of carrying phenotypes does decrease. But the phenotypes they can express are so few that they can be readily 'deciphered' by defective strains who also possess volume-similar capacitors, and thus quickly be invaded. In between comes the optimal phenotypic variations for cooperative strains to carry. When cooperative strains' capacitors with volume around K C they can effectively dodge defective strains' chasing and exploitation. Meanwhile, the cost of carrying such phenotypes can still be eclipsed by their mutual breed. So once these cooperative strains dominate the population, the evolutionary dynamics are quite stable in the sense that the dominance persists longer than other strains do.
These explanations can be better understood in a quantitative way with the help of typical evolutionary paths, one of which is illustrated in Fig. 4. For conveniently describing this evolutionary path, we denote by Low, Middle, and High level when the average phenotypic diversity is located in the interval (0, 5], (5,40], and (40,50], respectively. Denote by C L the state of cooperators with Low average phenotypic diversity, and the like. With respect to this typical evolutionary path, five states emerge: D M , C M , C M + D L , C M + D M , and C M + D H , and their corresponding fractions are respectively 0.659%, 99.099%, 0.0139%, 0.1899%, and 0.0389%. When the population is occupied by C M , the dynamics are extremely stable and stay put with a probability as high as 99.997%. It can but very rarely enter into either of the states C M + D L , C M + D M , C M + D H with probabilities 0.0002%, 0.0023%, and 0.0006%, respectively. Whenever residing in either of these three states, the population transits to the state C M with the probabilities 1.4347%, 1.2113%, and 1.4918%, correspondingly. Switching between these four states constitutes the core component of the population dynamics. As the rate of flowing into the state C M is trivial compared with that of flowing away from it, the population spends extremely high fraction of time in this state, and thus explains the macroscopic observations. We would like to point out that the population starting with the state D M can be stabilized at it with a very high likelihood of 99.997%, at least in this typical path, but transitions to this state are so meagrely few that this state can produce no essential effect in the long run. So far, our results have shown that cooperation coevolves with phenotypic diversity under a wide range of conditions and that there exists an optimal level of phenotypic diversity best promoting cooperation. It pays for cooperators to expand their capacitors (see Fig. 3a), thereby improving their opportunity of establishing cooperation via novel phenotypes, as these new phenotypes serve as secret handshakes that are difficult for defector to discover and chase after. Linearly increasing cost of expanding capacitor prevents capacitor's volume rising to infinity. Tradeoff between these two forces leads to symbiosis of optimal phenotypic diversity and highest fraction of such cooperative strain. Our results also show that evolved high levels of phenotypic diversity can occasionally collapse due to the invasion of defector mutants, illustrating that cooperation and phenotypic diversity can mutually reinforce each other.
We have also considered the effects of different time scale of interaction to selection on the stationary distributions of strains. Figure 5 shows that our main results still hold true as long as public goods interactions happen fast enough. One can also find that whenever interactions happen slowly, effects of game interactions are weakened. Strains carrying many phenotypes cannot offset the cost of carrying so many phenotypes and shall be eliminated soon. Evolutionary competition mainly unfolds between strains carrying a small number of phenotypes. Even so, cooperative strains are likely to dominate the evolutionary race (see Fig. 5a). In contrast, whenever rate of interaction is comparable (see Fig. 5b), or even faster than that of selection (see Fig. 5c and d), the evolutionary dynamics exhibit similar properties. Fraction of defective strain monotonically decreases as phenotypes this kind of strain carries increases. For cooperative strain, there exists an optimal level of phenotype to carry and corresponding cooperative strain accounts for the highest fraction. Thus, our results provide new insights into better understanding the coevolution of cooperation and phenotypic diversity.

Discussion
How to understand the rise and maintenance of altruistic behavior is a key problem in evolutionary biology 22,24-33 . Our study provides a possible path for cooperation to get established. Our study falls into the scope of chromodynamics of cooperation 5,6 , but is decisively different from previous studies concerning this topic. Generally speaking, the essence of chromodynamics is through decelerating the cooperator-defector interaction to promote cooperation. Ways achieving this purpose are various 3-6,10 , such as loosing the coupling of tag and strategy 5 and prescribing high levels of phenotypic diversity 6 . The mechanism we proposed that inherently carrying diverse phenotypes while expressing one plays a similar role as theirs 5 . In this sense, our mechanism is paralleled to the ones proposed in refs 5,6 and thus adds to the 'other mechanisms that can accomplish the same stabilizing effect' as the authors of ref. 5

have suggested.
Furthermore, our model is verily like the one proposed in ref. 34 in that interactions are local and competition is global. Group (or subpopulation) size changes in significantly different ways. In their model 34 , group size arriving at a threshold divides with certain probability and thus the group selection is involved. Our model does not involve group selection. Subpopulation expands or shrinks, entirely relying on the frequency-dependent selection. Though the significant difference, group splitting and inherent phenotypic variations both serve as avenues orienting the evolution of cooperation. Whether natural selection acts on the population just on the individual level as suggested in our model, or also by invoking group selection as assumed in the model of Traulsen et al., should be carefully treated.
Under natural selection, phenotypic diversity may be maintained by various mechanisms including mutant games 35 , sexual selection 36 , coevolving host-parasite population 37 , occasional recombinations of tags and of strategies 5 , and phenotype noise 11,[15][16][17][18][19]38 . In our model, the population equilibrates at high levels of phenotypic diversity when driven by the frequency-dependent evolutionary dynamics. Our model does not prescribe any given phenotypic diversity but lets individuals carry a capacitor of potentially expressible phenotypes. Capacitor Blue denotes cooperative strain, and red defective strain. The abscissa value represents the volume of capacitor. The evolutionary process is fully characterized in the main text. In panel a, whenever interactions happen slowly, accumulated payoffs resulting from game interactions are low. High cost of carrying many phenotypes cannot be counteracted. Therefore, strains carrying too many phenotypes shall be easily eliminated. Evolutionary competition mainly progresses between strains carrying a small number of phenotypes. Even so, cooperative strains are likely to dominate the evolutionary race. In panels b, c and d, whenever rate of interaction is comparable, or even faster than that of selection, the evolutionary dynamics exhibit similar properties. Fraction of defective strain monotonically decreases as phenotypes this kind of strain carries increases. For cooperative strain, there exists an optimal level of phenotype to carry and corresponding cooperative strain accounts for the highest fraction. Parameters: N = 40, r = 3.2, β = 0.1, and θ = 0.12. In a,b,c, and d, χ(x) = 0.01x, 0.5x, 2x, 5x, correspondingly, and the overall cooperation levels are 0.675, 0.805, 0.784, and 0.787, respectively.
SCiEnTifiC REPORTS | (2018) 8:4137 | DOI:10.1038/s41598-018-22477-5 volume is subject to evolutionary force. Phenotypic diversity is therefore an evolvable trait. On this front, our results provide an adaptive explanation for the phenotypic diversity.
Our model marries two recurrent themes in evolutionary biology: the evolution of cooperation and inherent phenotypic variations. This combination may well capture the rationale underlying many observations in biological circles. For example, Salmonella Typhimurium can be either cooperative by expressing virulence, or non-cooperative by not expressing 20,21 . This expression does determine strategy. But most likely it also sends signal to others who have not expressed the virulence. In fact, by a mathematical model Ackermann et al. have pointed out that phenotype noise itself is insufficient to evolve cooperation 20 . The side-blotched male lizards exhibiting diverse colors in throats may be another pertinent example. Our model philosophy may be conceived in the process of evolving this property 39,40 . These conjectures have implications for and awaits the confirmation of field studies.

Methods
Model description. Consider a well mixed population of finite size (=N). These N individuals compete to survive by reproducing asexually. Each individual i is characterized by a triplet (H i , S i , K i ), where H i is the current phenotype expressed, S i is the behavioral strategy, and K i is the volume of individual i's capacitor, which contains a number of phenotypes individual i can potentially express. Here S i is the probability that i cooperates when the public goods game happens. We limit our attention to two strategies, cooperation (S i = 1) and defection (S i = 0). It is straightforward to extend strategy to the full space as done in ref. 41 .
Each individual randomly expresses one phenotype from her capacitor while the rest are phenotypically silent. Only phenotypes expressed are visible and thus the population is divided into a number of subpopulations. Each subpopulation consists of individuals expressing the same phenotype (see Fig. 1).
Individuals in the same subpopulation interact by playing the public goods games. In this sense, phenotypes play the role of assortment. We assume that the number of group interactions happening in each subpopulation is proportional to the size of this subpopluation. When a subpopulation size is less than the group size as required for the public goods game to happen, no interaction takes place. Individuals need to bear the cost of retaining capacitor and expressing phenotypes. The cost is assumed to be proportional to the volume of her capacitor, that is, κ i (K i ) = θK i . Here we choose the simplest possible phenotype cost function. The net payoff π i determines the reproductive success (fitness) f i of an individual i. Here the fitness is an exponential function of payoff f e i i = βπ , where β is the intensity of selection.
The evolutionary updating occurs according to a frequency-dependent Moran process. At each time step, an individual is chosen with probability proportional to its fitness to reproduce an offspring. Following birth, a randomly selected individual in the population dies. The population size is thus constant throughout the evolution. Reproduction is however subject to mutation. As a first attempt, we would like to address whether cooperation can emerge and persist when individuals can carry redundant expressible phenotypes. Therefore, we assume that mutation is completely random so that any two different kinds of strains are mutaully accessible by mutation, and thus effects of mutation mode are unbiased. When mutation happens (with probability μ), the offspring randomly adopts one of the two behavioral strategies and also acquires a capacitor of volume K i ′ at a cost K i θ ′ . Here ′ K i is picked up from the set {1, 2, ..., M} equiprobably and M is the largest allowed capacitor. The mutant expresses one phenotype at random from this new capacitor, which we call random phenotype switching.
When it comes to the public goods game, g individuals decide whether or not to contribute to a common pool. A cooperator contributes c to the pool, while a defector contributes nothing. The sum contributed is multiplied by the factor r and then equally distributed to these g individuals, no matter whether they have contributed. Thus the payoff for a cooperator is c rjc g − and that for a defector is rjc g when there are j cooperators in the group. Without loss of generality and for simplicity 42 , we fixate the cost c of cooperation as 1. Under the classic constraint 1 < r < g, group of cooperators outperforms group of defectors whereas defectors are better off than cooperators in any given mixed group. Social dilemma arises as the best strategy for a player conflicts with the best strategy for the group.
We would like to point out that our present study can recreate qualitatively similar results presented in our previous work 38 by setting the group size as g = 2. However, group interactions cannot be simply regarded as the sum of pairwise interactions. Group interaction cannot proceed in this way that all possible groups would play the public goods game 43,44 . Otherwise payoffs from games are likely to be unrealistically high. More importantly, occurrence of group interaction requires higher complexity 10,45-50 . Fixation probability. We here present the general procedure of calculating the fixation probability in brevity. For sufficiently small mutation rate, the population simultaneously admits at most of two different types of individuals, say A and B. Here 'different' means that at least an element of the triplet (H A , S A , K A ) differs from that of the triplet (H B , S B , K B ). We should bear in mind that K A denotes the volume of A's capacitor, and S A her strategic behavior. K B is B's volume of capacitor and S B her strategic behavior. S A takes one if A is a cooperator, and zero otherwise. It is the same case with S B . Denote by i the number of A s in the population. The Moran process describing the evolutionary race has two absorbing states: i = 0 and i = N. When the population arrives at either of the two states, it will stay there forever. Denote by φ i the fixation probability that the population eventually ends up at the state i = N when starting with the state i.
When A and B have expressed the same phenotype, they are located in the subpopulation (also the whole population) and would play the public goods game. At this time, we can easily get their expected payoffs respect ively as located in different subpopulations. Public goods interactions just happen in each of the two subpopulations. In this situation, the payoff for an A and a B can be respectively written as P A = χ(i) · (r − 1)S A I (i>=g) − θK A and P B = χ(N − i) · (r − 1) S B I (N−i>=g) − θK B . Indication function I (i>=g) means that only when the subpopulation size arrives at least the group size can interactions take place. Otherwise, no interaction happens. The function χ(i) describes how the number of interactions depends on the subpopulation size. We choose linear function χ(x) = 0.1x. We shall consider other time scales 51 = β , respectively. The intensity of selection β measures how much payoff contributes to fitness. In an updating event, the population can add one, reduce one, remain the same in terms of the number of any two competing strains and thus the transition matrix. The normalized left eigenvector associated with the eigenvalue 1 of the transition matrix A provides the stationary distribution of these 2M full states. The overall cooperative level can be obtained by summing all the elements with even indices in normalized eigenvector 54,55 .
For slightly larger mutation rates, we can arrive at the population dynamics by the agent-based simulation.