key: cord-0782554-ppim2p2z authors: Ghanbarnejad, Fakhteh; Seegers, Kai; Cardillo, Alessio; Hovel, Philipp title: Emergence of synergistic and competitive pathogens in a co-evolutionary spreading mode date: 2021-08-10 journal: Phys Rev E DOI: 10.1103/physreve.105.034308 sha: 8ba25710006e1102af76c595927fec0a062c733c doc_id: 782554 cord_uid: ppim2p2z Cooperation and competition between pathogens can alter the amount of individuals affected by a co-infection. Nonetheless, the evolution of the pathogens' behavior has been overlooked. Here, we consider a co-evolutionary model where the simultaneous spreading is described by a two-pathogen susceptible-infected-recovered model in an either synergistic or competitive manner. At the end of each epidemic season, the pathogens species reproduce according to their fitness that, in turn, depends on the payoff accumulated during the spreading season in a hawk-and-dove game. This co-evolutionary model displays a rich set of features. Specifically, the evolution of the pathogens' strategy induces abrupt transitions in the epidemic prevalence. Furthermore, we observe that the long-term dynamics results in a single, surviving pathogen species, and that the cooperative behavior of pathogens can emerge even under unfavorable conditions. Understanding the diffusion of pathogenic agents is important as its aftermath reverberates on many aspects of our lives from health policies to economy, from politics to the transportation of people and goods [1] [2] [3] [4] [5] . Lately, the scientific community has devoted tremendous efforts in elucidating the dynamics behind these phenomena [6] [7] [8] [9] . Despite the many achievements of computational epidemiology, the spreading of multiple pathogens has received less attention than expected. This comes as a surprise as phenomena like comorbidity and cross-immunity constitute the norm rather than the exception. The former indicates the simultaneous presence of multiple diseases within the same host. The latter denotes the acquisition of immunity towards a certain disease as a result of infection by another one. In the epidemiology jargon, pathogens supporting comorbidity are indicated as cooperators, whereas the others as competitors. Cooperative pathogens may mutually promote their contagion like the Spanish flu and pneumonia [10] , or HIV co-infections [11] . Some types of influenza are, instead, examples of competitors [12] [13] [14] [15] [16] [17] . Moreover, comorbidity and cross-immunity are not observed exclusively among distinct pathogenic strains. For example, high levels of genetic diversity can provide a substrate for selection and rapid adaption, which are crucial to escape immune system recognition, * fakhteh.ghanbarnejad@gmail.com developing resistance to drugs, and adapt to new host types (spillover) [18] . Curiously, cases of cooperation have been reported also between distinct variants (quasi-species) of human H3N2 influenza [19, 20] . Consequently, the amount of literature about the spreading of cooperative [21] [22] [23] [24] [25] [26] [27] and competitive [28] [29] [30] [31] [32] [33] pathogens has grown over the years. Despite these evidences, the models developed hitherto assume that the strategy of a pathogen to cooperate -or not -is costless. In the evolutionary game framework, however, cooperation has a rather different meaning and implies always the payment of some costs [34] [35] [36] . For instance, the cooperative behavior of the Spanish flu exhibited towards secondary infections led to high death rates of the host subjects, making this synergistic epidemic -or syndemic -prone to relatively "quick" disappearance [10] . Notably, the evolution of spreading pathogens has been rarely considered [37] [38] [39] [40] . In this paper, we envision an evolutionary scenario where interacting pathogens have two different strains with cooperative and defective strategies, and evolve while maximizing their own benefit. We extend the two-strain susceptible-infected-recovered (SIR) model of Chen et al. [23] by intertwining it with an evolutionary infection game, and show under which conditions synergy and competition emerge. We observe that the evolution of the pathogens' strategies reverberates on the epidemic prevalence. Moreover, the dynamics foster more the survival of single strains, i.e., strategy, rather than the co-existence of multiple strains. We combine two distinct processes: (i) the simultaneous spreading of two different pathogens and (ii) the evolution of their strategies. The former takes place over a short, continuous, time-scale (within the season) indicated by variable t; whereas the latter occurs on a longer, discrete, one (between two seasons) denoted by variable T . The two processes are intertwined by using the outcome of the spreading as the input of the evolution and vice-versa (see Fig. 4 ). The two phases of spreading and evolution of strategies take place sequentially in a cyclic way until a global stationary state (GSS) is reached. We explain first the spreading and evolutionary game dynamics separately, and then describe how these processes (named phases 1 and 2) are combined together. In each season, the model describes the simultaneous spreading of two pathogens A and B that obey an extended SIR compartmental dynamics [23] . Each pathogen, σ X , has two strains corresponding to its cooperator (C) and defector (D) strategies. Cooperators (σ X C ): A strain of pathogen σ X that cooperates with other strains. Defectors (σ X D ): A strain of pathogen σ X that competes with other strains. Following the above definitions, the pathogenic population is composed by four different species, namely: A C , A D , B C , and B D . These species follow an SIR dynamics and interact with each other yielding a 25 states model whose transitions' diagram is shown in Fig. 1 . In agreement with the case without strategies, and without loss of generality, the recovery rates are set equal to one (i.e., r = r = 1). Moreover, we assume that the rates of infection depend only on the strategy of the pathogen occupying the host, which is equivalent to imposing a symmetry in the parameters' space. According to such a symmetry, we can write the infection rates of empty hosts as α C = α D = α. Accordingly, for a host already occupied by one pathogen, we pass from four rates (i.e., β CC , β CD , β DC , and β DD ) denoting all the possible combinations of pairs of strategies to just two of them, namely: β CC = β DC = β C and β DD = β CD = β D . After imposing the aforementioned symmetries, we can write the system of ODEs describing the spreading dynamics using an approach similar to that introduced in [23] . Hence, we can define the groups of individuals actively infected by a species as: (1) Transition scheme of the two strategies double SIR spreading dynamics. We display all the possible transitions among compartments in the multi disease propagation of pathogens A and B. The introduction of strategies for pathogens expands the original transition scheme introduced in [23] , which now translates into a 25 states diagram. Capital letters denote infected states, whereas small letters denote recovered ones. The infection transition is denoted by a single arrow for simple pathogen contagion, and by double arrows for multiple pathogens infection. Capital letters denote infected states, whereas small letters denote recovered ones. The notation [·] indicates the density of individuals in a given state, while the sum of all 25 states is normalized to 1. Considering the definitions of infected groups displayed in Eqs. (1), we can write the complete set of equations modeling the spreading dynamics [Eqs. (2) ]. Despite the high number of equations, the system in Eq. (2) depends only on three parameters: α, β C , and β D . The next section illustrates how these parameters are related with each other. [ , As mentioned above, the dynamics of the double SIR with cooperative and defective pathogens (i.e., the quadruple SIR) depends on three parameters: the empty host infection rate α, and the infection rates for hosts occupied either by a cooperator (β C ) or by a defector (β D ) pathogen. In the following, we show how these quantities are related to one another. According to the definition of cooperative and defective strategies made in the previous section, it is natural to impose that β C > β D . Moreover, as discussed in [22] , the hierarchy between the empty host and occupied host infection rates determines whether two pathogens act in symbiosis (i.e., their strategy is cooperation) or in antagonism (i.e., their strategy is defection). The former case occurs when α < β, whereas the latter when α > β. These hierarchies between the infection rates denote the existence of a relationship between them. Amidst the plethora of possible choices, we decided to adopt the following one: with c ∈ ]0, ∞[ being a parameter. Imposing a relationship between the infection rates reduces the number of free parameters in Eqs. (2) from three to two: α and c. Figure 2 displays the values of β C and β D as a function of c for four different values of α: two sub-critical (i.e., α ≤ 1) and two super-critical (i.e., α > 1). The visual inspection of Fig. 2 reveals that, in agreement with Eq. (3), for c = 1 one gets β C = β D = α. Such a case corresponds to the neutral spreading scenario in which the rate of infection is independent on the occupation state of the host, as well as on the strategy of the pathogen occupying it. Hence, the only factor influencing the outcome of the dynamics is the value of α which controls the epidemic prevalence 1−S ∞ . The value c = 1 corresponds also to the point where the hierarchy between β C and β D changes. In Fig. 3 we plot the hierarchy between β C and β D as a function of α and c. In particular, the region c < 1 corresponds to the case β D > β C corresponding to the scenario in which it is easier to infect a host occupied by a pathogen acting as a defector than a cooperator. Although such a scenario is biologically not meaningful, we nevertheless explore it for the sake of completeness. In summary, assuming that it is easier to infect a host already infected by a cooperator pathogen than by a defector one, we set β CC = β DC = β C = α c and The quadruple SIR model encoded by Eqs. (2) and the hierarchies between the infection rates analyzed in Sec. II B describe the spreading of pathogens acting in symbiosis or antagonism with the other pathogens. However, the sole consequence of adopting a certain strategy is how easily a pathogen can infect a host occupied by another pathogen. Nevertheless, from the biological point of view, facilitating (or not) the infection from another pathogen could be thought -among other things, -as a proxy for the will of the occupying pathogen to share (or not) the host's resources with the invader. Following this hypothesis, one could model the act of infecting a host as a game and, consequently, assign a payoff to it [35] . The payoff accumulated by all the pathogens of species X and strategy Y , Π σ X Y , constitutes the fitness of that pathogen's type which, in turn, determines its ability to reproduce. Hence, we can use the fitness of each pathogen's type computed at the end of the T -th season of the spreading process to compute its abundance in the initial seed of the (T + 1)-th season. Such an approach allows to describe features like comorbidity and cross-immunity as the byproduct of natural selection. In light of the above reasoning, the act of infecting a host splits into two main scenarios: one in which the host is empty, and another in which it is already occupied by another pathogen. Assuming that the total amount of resources available in the host is equal to one, let us discuss the two scenarios separately. If the host is empty, then the pathogen infecting the host will have access to all of its resources regardless of its strategy. Therefore, the payoff, π, associated to the event of infecting an empty host (i.e., single infection) is equal to π = 1. If the host is occupied instead (i.e., in the secondary infection event), the value of the payoff depends on the combination of the strategies of the pathogen infecting the host, and of the pathogen already present within the host. Since we assume that pathogens can either cooperate or defect, and that cooperators are more keen to share the host's resources than defectors, we can adopt the payoff matrices of the Hawk and Dove (HD) game [34, 41] , given by: with γ ∈ 0, 1 2 . The payoff, π X,Y , of a pathogen with strategy X infecting a host infected by another pathogen with strategy Y corresponds to the element a X,Y (X, Y ∈ {C, D}) of A. Analogously, the payoff of a pathogen with strategy Y occupying a host that gets infected by another pathogen with strategy X corresponds to the element a X,Y of A . The fitness of species σ, Π σ , depends on the history (i.e., the sequence) of its contagion record (events). Finally, the concentration, ρ i , of species i in the spreading seed of season T + 1 is regulated by the so-called replicator equation [35, 36] , given by: , and t ∞ denotes the time t at which the spreading dynamic reached its stationary state. According to Eq. (5), species with fitness higher than the average will proliferate, whereas those with fitness lower than the average will become extinct [36] . to Eqs. (2) (whose reaction kinetics is summarized in the top left panel of Fig. 4 and details are shown in Fig. 1 ) and accumulate payoff via infection events according to the infection game introduced in Sec. II C (Fig. 4 top right panel). Specifically, a pathogen X Λ infecting a susceptible host S receives a payoff π XΛ = 1 regardless of its strategy. If, instead, the host is already infected by another pathogen Y Γ -i.e., in a secondary infection,the payoff of the infecting pathogen, π XΛ , is given by the (Λ, Γ) element of matrix A of Eq. (4) according to the four possible combinations of the (infecting, infected) pathogens' strategies: (C, C), (C, D), (D, C), and (D, D) (similar reasoning applies to payoff π YΓ using matrix A ). The choice of the payoff scheme (i.e., the game) associated with the secondary infection event is dictated by the idea of how the pathogens share the host's resources. Moreover, the idea that cooperators will share the host's resources while defectors will try to fight to seize them all, reverberates also on the properties of the secondary infection rates β C and β D which: i) depend exclusively on the strategy of the pathogen already occupying the host, and ii) controls the hierarchy existing among them (i.e., β C > β D ). For these reasons, a payoff scheme like that of the HD game fits quite well with our model. Considering other payoff schemes like the Stag-Hunt one which are typical of "coordination dynamics" (i.e., where the best strategy is to coordinate with the other player and adopt the same strategy) clashes with the idea of having finite host's resources. A similar reasoning applies to the case of the Prisoner's Dilemma (PD) [36] . Once the spreading process reaches its stationary state, i.e., t = t ∞ = ∞, the values of ρ, for the next season T +1 are given by Eq. (5) (Fig. 4, phase 2 ). The mixture of species in the seed evolves season after season until ρ T i reaches a stationary state, i.e., ρ T +1 i = ρ T i ∀i. Then, the co-evolutionary dynamics is considered at the GSS which is labeled as evolutionary season T = ∞. We start our analysis by looking at the state of the system at the end of the spreading dynamics, t ∞ , and compare it at the first (T = 1) and the last (T = ∞) season. In Fig. 5 we portray the outcome of the co-evolutionary dynamics from both the pathogens (left panels) and the host (right panels) perspectives. If the composition of the spreading seed, , is richer in cooperation, i.e., ∆ CD | t0 > 0, then cooperator strains will easily take over. Nevertheless, there is a region of the (α, c) space where cooperation still has the chance to thrive even though the seed is dominated by defectors (i.e., ∆ CD | t0 < 0). Such regions correspond to the brighter areas displayed in Figs. 5(a) and (b) , where the colors encode the density of the pure cooperator and defector strains at the end of the spreading season, ables refer to the recovered compartments. We find that the co-evolutionary dynamics amplifies the differences of ∆ CD | t∞ observed for the first season T = 1 [panel (a)], and a region where cooperation thrives emerges, i.e., the bright area in panel (b) corresponding to ∆ CD | t∞ > 0. In addition, we observe another region where defection prevails (∆ CD | t∞ < 0). Figures 5(c) and (d) capture the effects of the evolution from the hosts' perspective through the lens of disease incidence (or transitivity), 1 − S ∞ . As we vary α [panel (c)], we observe a transition from lower incidences for α ≤ 1, towards complete infection, i.e., 1 − S ∞ 1. In the absence of evolution, that is, for a single season only (T = 1), the transition between these regimes is smooth, whereas the catalytic effect of the evolutionary dynamics, i.e., over many seasons (T = ∞), triggers the appearance of a gap in the transition. Even though at first glance the gap resembles the discontinuous transitions observed in the cooperative SIR co-infection models [23, 42] , they are not the same. In fact, in non-evolutionary SIR coinfection models the gap occurs exclusively for cooperator only pathogens, whereas purely defector pathogens do not make any gap (see the competitive regime in Ref. [23] ). Instead, in our case both cooperators and defectors spread at the same time, thus highlighting the inability of the latter to suppress the gap. Also, the transition occurs at smaller values of c and greater values of α, compared to cooperative only SIR dynamics, (see Eq. (7) in Ref. [42] ). The effect of varying c on the disease incidence [panel (d)] is less strong. Notwithstanding, we still observe a gap in the transition for T = ∞. In the neutral spreading setup β C = β D = α (i.e., c = 1), we do not observe any difference between the incidences measured in the absence and presence of seed evolution. Such a similarity suggests that the sole accumulation of payoff is not enough to induce differences in the phenomenology. Finally, we stress that the region corresponding to c < 1 is biologically not meaningful as β C < β D . A way to study the evolution of the dynamics is to look at the evolution across the seasons of the spreading seed. More specifically, the sum of the densities of pathogens' species at the beginning of each spreading season -i.e., the size of the spreading seed, -is fixed and equal to ω ∈ ]0, 1[. As the pathogens' densities are all positive definite, and their sum is constant, the equation: certain conditions it is possible to project the density space onto a bidimensional surface. To this aim, let us consider two variables 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1 fulfilling the following conditions: = ω x y , The x coordinate accounts for the pathogen's species concentration with x = 0 (x = 1) denoting a population fully made of pathogens of species B (A). The y coordinate, instead, accounts for the pathogen's strategy, with y = 0 (y = 1) denoting full defection (cooperation). Not all the possible combinations of the seed components (only three of them are truly independent) can be mapped using the present approach. Still, the bidimensional representation allows one to scrutinize the dynamics between the most relevant points. Figure 6 provides a visual summary of the bidimensional projection. The diagram's corners correspond to pure seeds (i.e., made of only one species with one strategy), with coordinates ( 1 2 , 1 2 ). Therefore, the path connecting point P ≡ (x P , y P ) with point Q ≡ (x Q , y Q ) denotes a variation in the composition of the spreading seed induced by the evolutionary component of our model. We continue our analysis by studying the evolution of the epidemic seeds between seasons, to investigate how the interplay between the two phases of our model impacts the overall dynamics. To this aim, Fig. 7 displays how the concentrations of species ρ A C , ρ B C , ρ A D , and ρ B D in the seed change season after season until reaching the GSS. The main observations are summarized as follows: • Contrary to the expected behavior of the HD game, the pure states correspond to stable fixed points of the dynamics, whereas the balanced states are saddle points, i.e., they are stable in one direction and unstable in the other. Hence, no matter what the initial composition of the seed is, only one species with one strategy survives at the GSS, while all other species become extinct. • Two unstable manifolds split the phase portrait according to each of the pathogens' features, i.e., species and strategy. • Figures 7(a) and (b) show that, when the simple infection rate is sub-critical (i.e., α < 1), the structure of the phase portrait and the epidemic prevalence are not affected by variations of c. • For supercritical values of the simple infection rate [α > 1, i.e., Figs. 7(c) and (d)], increasing c translates into an expansion of cooperation's basin of attraction. Hence, seeds initially composed of more defectors than cooperators can evolve towards fully cooperative strains. • Even though changing the value of γ displaces the intermediate stable fixed point in the "classical" HD game dynamics; we have verified that both the position of the unstable manifolds and the epidemic prevalence are robust to variations of γ. In this paper, we have studied the emergence of synergistic and competitive traits in a co-evolutionary model intertwining epidemic spreading and evolutionary game theory. To this aim, we have extended the SIR model of Chen et al. [23] by assuming that pathogens have two strategies (cooperate or defect), and allowing that their concentration inside the epidemic seed (whose total size is kept fixed throughout the whole dynamics) evolves according to the replicator equation [34] . For the latter, the pathogens accumulate payoff according to an infection game whose payoff scheme for secondary infections events maps onto the hawk-and-dove game. The resulting co-evolutionary model displays features i ρ T i |t 0 . The colors encode the transitivity at the GSS. The red, dashed line is an approximation to guide the eye finding the manifold separating cooperation from defection. The arrows indicate the displacement of the seed between epidemic seasons. We set i ρ T i |t 0 = 0.05 and γ = 0.25, and use four sets of (α, c) pairs. See Fig. 6 of for the details of the bi-dimensional mapping. that are not present in either spreading or evolutionary game dynamics taken singularly. In particular, the game reverberates on the outcome of the spreading dynamics. It alters the disease prevalence in the host and can trigger the emergence of a gap in the transition. This phenomenon has also been observed for non-evolutionary SIR models where two pathogens only cooperate and sharp transitions occur for c > 2 and α < 1 [23, 42] . Moreover, the interplay between the spreading and game dynamics alters the stability of the state space's fixed points of the game. Pure seeds are stable solutions of the dynamics, contrary to the expected stable mixed scenario typical of the hawk-and-dove game [36] . We have also observed the emergence of a region of the phase space where cooperation thrives even under unfavorable conditions. Such a behavior is in agreement with the phenomenology reported for stochastic games where cooperation emerges even in conditions where usually it should not [43] [44] [45] . Summing up, the overall phenomenology exhibited by our model cannot be traced back to one of the two processes separately but emerges from their interplay. The framework presented here constitutes a way to understand under which conditions synergy, i.e., comorbidity, and competition, i.e., cross-immunity, can emerge through an evolutionary framework where the hosts and the pathogens mutually interact. The present work is a first attempt towards understanding multiple pathogen spreading and opens new roads for theoretical modeling of multistrain epidemic processes from the perspective of evolutionary theory. Despite the richness of the dynamical scenarios observed, we made simplifying assumptions. For instance, we have considered a symmetric scenario for the primary infection rates and have chosen a specific set up for the secondary infection rates. Moreover, the values of the parameters do not evolve across time, and we have considered a well mixed population. Relaxing these assumptions can lead to even more realistic dynamics. Proceedings of the National Academy of Sciences of the United States of America Proceedings of the National Academy of Sciences of the United States of America International Health Regulations (IHR) Evolution and the Theory of Games Evolutionary games and population dynamics Evolutionary Dynamics (Harvard University Press Modeling Infectious Diseases in Humans and Animals Parasite community ecology and biodiversity The authors thank the members of the GOTHAM lab (and in particular D. Soriano Paños) for helpful comments and discussions. AC acknowledges the financial support of SNSF through the project CRSII2 147609 and of the Spanish Ministerio de Ciencia e Innovación (MICINN) through Grant IJCI-2017-34300. F.Gh. acknowledges the partial support by Deutsche Forschungsgemeinschaft (DFG) under the grant project: 345463468 (idonate). Graphics have been prepared using the matplotlib python package [46] .