key: cord-0262435-wudh9n9d
authors: Sasaki, Akira; Lion, Sébastien; Boots, Mike
title: The impact of antigenic escape on the evolution of virulence
date: 2021-01-19
journal: bioRxiv
DOI: 10.1101/2021.01.19.427227
sha: c00e2dd2f63c3b690fb3511446eb513f1d122b0f
doc_id: 262435
cord_uid: wudh9n9d

Understanding the evolutionary drivers determining the transmission rate and virulence of pathogens remains an important challenge for evolutionary theory with clear implications to the control of human, agricultural and wildlife infectious disease. Although disease is often very dynamic, classical theory examines the long-term outcome of evolution at equilibrium and, in simple models, typically predicts that R0 is maximized. For example, immune escape may lead to complex disease dynamics including repeated epidemics, fluctuating selection and diversification. Here we model the impact of antigenic drift and escape on the evolution of virulence and show analytically that these non-equilibrium dynamics select for more acute pathogens with higher virulence. Specifically, under antigenic drift and when partial cross immunity leads to antigenic escape, our analysis predicts the long-term maximization of the intrinsic growth rate of the parasite resulting in more acute and virulent pathogens than those predicted by classic R0 maximization. Furthermore, it follows that these pathogens will have a lower R0 leading to implications for epidemic, endemic behavior and control. Our analysis predicts both the timings and outcomes of antigenic shifts leading to repeated epidemics and predicts the increase in variation in both antigenicity and virulence before antigenic escape. There is considerable variation in the degree of antigenic escape that occurs across pathogens and our results may help to explain the difference in virulence between related pathogens most clearly seen in the human A, B and C influenzas. More generally our results show the importance of examining the evolutionary consequences of non-equilibrium dynamics.

Recent epidemics emphasize that infectious disease remains a major problem for human health and agriculture [1] [2] [3] [4] and they are increasingly recognized as important in ecosystems and conservation [5, 6] . This importance has led to the development of extensive theoretical literature on the epidemiology, ecology and evolution of host-pathogen interactions [7] [8] [9] [10] . Understanding the drivers of the evolution of virulence, typically defined in the evolutionary literature as the increased death rate of individuals due to infection, is a key motivator of this theoretical work [8, [10] [11] [12] [13] [14] . Generally, models assume that a higher transmission rate trade-offs against the intrinsic cost of reducing the infectious period due to higher death rates (virulence) of rapidly dividing strains and classically predict that evolution maximizes the parasite epidemiological R 0 leading to an optimal virulence [8, [10] [11] [12] [13] [14] . In fact, this result only holds in models where ecological feedbacks take a constrained form, such that even relatively simple processes such as densitydependent mortality, multiple infections and spatial structure may lead to diversification or different optima [10, 12, 13, 15] . Moreover, this classic evolutionary theory examines the long-term equilibrium evolutionary outcome in the context of stable endemic diseases, but in nature infectious diseases often exhibit complex dynamics, with potentially important impacts on pathogen fitness [16] [17] [18] [19] .

Antibody-mediated immunity is a critical factor driving the dynamics of important infectious diseases such as seasonal influenza, leading to selection for novel variants that can escape immunity to the current predominant strain [20] [21] [22] . Such antigenic escape typically causes the optimal strain of the parasite to change through time as it moves through antigenic space. Moreover, partial cross-immunity between the different parasite strains may lead to recurrent epidemics, fluctuations in parasite strains and potentially strain coexistence [23] [24] [25] [26] . Current theory has shown that the evolution of immune escape can lead to dramatic disease outbreaks [24] [25] [26] , but the implications of these epidemiological dynamics for the evolution of disease virulence are unknown. This question is challenging in part because much of the theoretical framework used to study virulence evolution typically considers diseases that are at an endemic equilibrium [8, [10] [11] [12] [13] [14] . As such we currently lack a broad theoretical understanding of the evolution of virulence in the presence of antigenic escape, despite its importance as an epidemic process and the likely implications of its inherently dynamical epidemic nature.

Here we examine the implications of antigenic escape for the evolution of infectious disease in the context of the well-studied transmission/virulence trade-off [10, 27] . For simplicity we first examine analytically the case without cross-immunity and then apply a new 'oligomorphic' analysis that combines quantitative genetic and game theoretical approaches [28] to examine the impact of cross-immunity that leads to antigenic jumps and epidemic outbreaks. With this analysis we are able to study the evolutionary outcome across a range of ecological and evolutionary time scales allowing us to examine evolutionary outcomes under non-equilibrium conditions. Our key result is that antigenic escape selects for higher transmission and virulence due to the repeated epidemics caused by immune escape, leading to the long-term persistence of acute pathogens. Indeed, antigenic escape has the potential to select for infectious diseases with substantially higher virulence than that predicted by the maximization of R0 in classic disease models leading to the evolution of much more acute diseases.

We consider a population of pathogens structured by a one-dimensional antigenic trait , so that ( , ) is the density of hosts infected with antigenicity strain at time . Following Gog and Grenfell [29] , we assume that an individual is either perfectly susceptible or perfectly immune to a strain. A strain of pathogen can infect any host, but will be infectious only when the host is susceptible to that strain. When a strain of pathogen infects a host that is susceptible to a strain , the host may become (perfectly) immune to the strain with probability ( − ). This is the partial cross immunity function between strain and , that takes a value between 0 and 1 and is a decreasing function of antigenic distance | − | between strain and . The density of hosts susceptible to antigenicity strain at time is noted ( , ).

Assuming that all pathogen strains have the same transmission rate and virulence , we can describe the dynamics with the following structured Susceptible-Infected-Recovered model:

where is the recovery rate, and = $ # /2 is one half of the mutation variance for the change in antigenicity, representing random mutation in the continuous antigenic space. The dynamics for the density of recovered hosts is omitted from (1) as it does not affect the dynamics (1) for the densities of susceptible and infected hosts.

Invasion of a single pathogen: in our first scenario, we start with a population where all hosts are susceptible to any strain ( (0, ) = 1) and a pathogen strain with antigenicity trait = 0 is initially introduced, so that (0, ) = ( ), where is a small positive constant and (⋅) is Dirac's delta function. The system then exhibits travelling wave dynamics in antigenicity space. At the front of the travelling wave, ( , ) is sufficiently small and ( , ) is sufficiently close to 1. Eq. (1a) can then be linearized as

where = − ( + ) is the rate of increase of an antigenicity strain before it spreads in the population and causes the build-up of herd immunity. The system (1) asymptotically approaches travelling waves of both pathogen antigenicity quasispecies ( , ), which have an isolated peak around the current antigenicity, and host susceptibility profile ( , ), that smoothly steps down towards a low level after pathogen antigenicity quasispecies passes through, with a common constant wave speed [30] (Fig 1A) = 2√ = 2FG − ( + )H . Evolution of antigenic escape with cross-immunity: To predict how cross-immunity affects the evolution of antigenic escape, we use an oligomorphic dynamics analysis [28] . We consider a population composed of different antigenicity strains (or morphs), that can be viewed as quasispecies. The analysis in Appendix A allows us to track the dynamics of morph frequencies, % , and mean trait values, ̅ % , as:

where ( ) is the susceptibility profile of the population, which depends on the cross-immunity function , ̅ % is the mean susceptibility perceived by viral morph , and ̅ the mean susceptibility averaged over the different viral morphs. Note that, in general, ( ), ̅ % and ̅ will be functions of time, as the susceptibility profile is molded by the epidemiological dynamics of ( , ) and ( , ).

Equation (4a) reveals that, as expected, morph will increase in frequency if the susceptibility of the host population to this strain is higher on average. Equation (4b) shows that the increase in the mean antigenicity trait of morph depends on (i) the variance of the morph distribution, % , (ii) the transmission rate, and (iii) the slope of the susceptibility profile close to the morph mean ̅ % . Together with an equation for the dynamics of variance under mutation and selection (SI: Appendix S1), equations (4a) and (4b) allow us to quantitatively predict the change in antigenicity after a primary outbreak, as shown in Figure 2 .

For instance, after a primary outbreak caused by a strain with antigenicity ̅ ' = 0 at = 0, the susceptibility profile is approximately constant and given by ( ) = (1 − ' ) ((*) , where ' is the final size of the epidemic of the primary outbreak at antigenicity = 0 (SI: Appendix S1). Thus, for a decreasing cross-immunity function, ( ), the slope of the susceptibility profile is positive, which selects for increased values of the mean antigenicity trait ̅ , of a second emerging morph (SI: Appendix S1). As the process repeats itself, this leads to successive jumps in antigenic space. In addition, a more peaked cross-immunity function, , yields larger slopes to the susceptibility profile and thus selects for higher values of the antigenicity trait. High cross immunity (ω=0.6) A Low cross immunity (ω=0.2)

We now extend our analysis to account for mutations affecting pathogen life-history traits such as transmission and virulence. To simplify, we use the classical assumption of a transmission-virulence trade-off [8,10-14] and consider that a pathogen morph, , has frequency, % , mean antigenicity trait, ̅ % , and mean virulence Y % . In Appendix S2, we show that the morph's mean traits change as

where G % is the genetic (co)variance matrix, and the vector on the right-hand side is the selection gradient. Note that, while the selection gradient on antigenicity depends on the slope of the antigenicity profile at the morph mean, the selection gradient on virulence depends on the slope of the transmission-virulence trade-off at the morph mean, weighted by the susceptibility profile at the morph mean.

Assuming we can neglect the build-up of correlations between antigenicity and virulence due to mutation and selection, the genetic (co)variance matrix is diagonal with elements % * and % -. Then, as shown previously, antigenicity increases if the slope of the susceptibility profile is locally positive, while mean virulence increases as long as ′( Y % ) > 1/ ( ̅ % ). However, when antigenicity can evolve, selection will also lead to the build-up of a positive covariance between antigenicity and virulence, resulting in a synergistic effect (SI: Appendix S2; Figure S3 ). As the antigenicity trait increases, the evolutionary trajectory of virulence converges to the solution of & ( ) = 1 which corresponds to maximizing the rate of increase of pathogen ( ) = ( ) − ( + ) in a fully susceptible population. This is equivalent to maximizing the wave speed ( ) = 2`( ) , as shown in Appendix S3. Figure 3a shows that, in the absence of cross-immunity, the ES virulence is well predicted by r maximization. With cross-immunity (Figure 3b ), virulence evolution is characterized by jumps that reflect the sudden shifts in antigenicity due to cross-immunity.

As such antigenic escape selects for higher transmission and virulence and more acute infectious diseases. This has parallels with the results that show that there is a transient increase in virulence at the start of an epidemic with r rather than R0 being maximized [16,17,19,32], but it is important to note that here we predict the long-term persistence of highly transmissible and virulent disease strains due to the existence of antigenic escape.

Although we have so far assumed a never-ending antigenic escape process, it is easy to extend our analysis to consider that antigenic escape is constrained by pleiotropic effects. Then, once the antigenicity trait has stabilized, the ES virulence would satisfy

where = / -measures the correlation between antigenicity and virulence. Thus, the slope to the transmission-virulence trade-off at the ESS now takes an intermediate value between /( + ) and 1, as shown in figure 4 .

Although our analysis allows us to understand the long-term evolution of pathogen traits, it can also be used to accurately predict the short-term dynamics of antigenicity and virulence. We now consider that a primary outbreak has molded a susceptibility profile ( ) that we assume constant. Although this assumption will cause deviations from the true susceptibility profile, it allows us to decouple our evolutionary oligomorphic dynamics from the epidemiological dynamics. Figure 5 shows that the approximation accurately predicts the jump in antigenic space and joint increase in virulence during the secondary outbreak. The accuracy of the prediction depends on the time at which we seed the oligomorphic dynamical system, as detailed in Appendix S1, but remains high for a broad range of values of this initial time. Hence, our analysis can be used to successfully predict the trait dynamics after the emergence of a new antigenic strain. The oligomorphic dynamics describing the changes in the frequency ' ( ) = 1 − , ( ) of the currently prevailing morph at time 1 and the frequency , ( ) of upcoming morph, the mean antigenicity ̅ % ( ) and mean virulence Y % ( ) of the two morphs ( = 0,1), and the within-morph variances % * and % in antigenicity and virulence, as well as the within-morph covariance % (t) between antigenicity and virulence in each morph ( = 0,1) are defined as (A44)-(A49) in Appendix S2.

We have shown how antigenic escape selects for more acute, virulent infectious diseases with higher transmission rates that cause increased mortality in infected hosts. This result is important given the number of important infectious diseases such as seasonal influenza that have epidemiology driven by antigenic escape. Until recently the evolution of virulence literature has mostly focused on equilibrium solutions that in simple models lead to the classic idea that pathogens evolve to maximize their basic reproductive number R 0 [8,10-14]. Our results show that the process of antigenic escape leading to the continual replacement of strains [23-26], creates a dynamical invasion process that in and of itself selects for more acute, fast transmitting, highly virulent strains that do not maximize R 0 . This has parallels with the finding that more acute strains are selected temporally at the start epidemics [24-26], but critically, in our case the result is not a short-term transient outcome. Rather, the eco-evolutionary process leads to the long-term persistence of more acute strains with higher transmission rates causing higher virulence. As such, antigenic escape may be an important driver of high virulence in infectious disease.

In the simpler case where there is no cross immunity, there is a travelling wave of new strains invading due to antigenic escape. In this case we can use established methods to gain analytical results that not only predict the speed of change of the strains, but also the evolutionarily stable virulence. With our model's assumptions, without antigenic escape we would get the classic result of the maximization of the reproductive number R 0 [8,10-14], but once there is antigenic escape we show analytically that the intrinsic growth rate of the infectious disease r is maximized. Maximizing the intrinsic growth rate leads to selection for higher transmission leading to a much higher virulence. Effectively this is the equivalent of an infectious disease "live fast, die young" strategy. The outcome is due to the dynamical replacement of strains, with new strains invading the population continually leading to a continual selection for the strains that invade better [24-26]. As such we predict that any degree of antigenic escape will in general select for more acute faster transmitting strains with higher virulence in the presence of a transmission-virulence tradeoff.

Partial cross-immunity leads to a series of jumps in antigenic space that are characteristic of the epidemiology of a number of diseases and in particular, in humans, the well-known dynamics of influenza A & B [23-26]. Here a cloud of strains remains in antigenic space until there is a jump that, on average, overcomes the cross immunity and leads to the invasion of a new set of strains that are distant enough to escape the immunity of the resident strains [23] [24] [25] [26] . In order to examine the evolutionary outcome in this scenario we applied a novel oligomorphic analysis [28] . Again, we find that antigenic escape selects for higher virulence with more acute faster transmitting strains being favored with again selection leading toward the maximization of the intrinsic growth rate r. Both our analysis and simulations show that in the long term the virulence increases until it reaches a new optima potentially of an order of magnitude higher than would be expected by the classic prediction of maximizing R 0 . We show that virulence increases between each antigenic jump, falling slightly at the next jump before increasing again until it reaches this new equilibrium. It is also important to note that diversity in both antigenicity and virulence increases as we move towards the next epidemic, maximizing just at the point when the jump occurs. This increase in diversity could in principle be used as a predictor of the next jump in antigenic space.

Clearly the virulence of any particular infectious disease depends on multiple factors, including both host and parasite traits, and critically the relationship between transmission and virulence. This makes comparisons of the virulence across different infectious diseases problematic since the specific trade-off relationship between transmission and virulence is often unknown. However, our model shows that antigenic escape will, all things being equal, be a driver of higher virulence favoring more acute strains. It is also important to note that since antigenic escape is a very general mechanism that selects for higher virulence it follows that we may see high virulence in parasites even when the costs in terms of reduced infectious period are substantial. Amongst the influenzas, influenza C does not show obvious antigenic escape [33] and is typically much less virulent than the other influenzas which are the classic examples of infectious disease with antigenic escape. Furthermore, Influenza A shows much more antigenic escape than Influenza B and again in line with our predictions typically influenza A is the more virulent [34, 35] . Clearly these differences can be ascribed to multiple factors and indeed the higher virulence of influenza A is often posited to be due to a more recent zoonotic emergence [36] . However, our models suggest that the differences in the degree to which they show antigenic escape per se may also contribute to these differences. Clearly there are also highly virulent pathogens that do not show antigenic escape and a formal comparative analysis is confounded by multiple factors, but the evidence from influenza is consistent with our predictions.

An important implication of our work is that it follows that antigenic escape selects for strains with a higher virulence than the value that maximizes R 0 and therefore leads to the evolution of infectious diseases with lower R 0 . One simplistic implication of this is that, from this point of view, diseases with antigenic escape may be easier to eliminate and control with vaccination. Of course, in practice the opposite is often true since producing an effective vaccine is much more problematic when there is antigenic escape [33, 37] . In addition, it follows that epidemics will tend to be less explosive than they otherwise would have been, having a lower peak but lasting longer, with evolution here effectively 'flattening the curve'. Infectious diseases that show jumping antigenic escape are characterized by repeated epidemics, but our work suggests that due to the selection for a lower R 0 the eco-evolutionary feedback will have significantly impacted the pattern of these epidemics. This is an important example of a dynamical ecological/epidemiology feedback into evolutionary outcomes that in turn then feedbacks into the epidemiology characteristics of the disease.

We have used oligomorphic dynamics [28] to make predictions on the waiting times and outcomes of the antigenic jumps in our model with cross immunity. This approach tracks both mean trait values and the changes in trait variances in models with explicit ecological dynamics. As such it combines aspects of eco-evolutionary theory [38,39] and quantitative genetics approaches [40,41] to provide a more complete understanding of the evolution of quantitative traits. Our approach can assume a wide range of different ecological and evolutionary time scales and therefore allows us to address fundamental questions on ecoevolutionary feedbacks and on the separation between ecological and evolutionary time scales. This is important since it allows us to test the implications of the different assumptions of classical evolutionary theory and to better understand the role of eco-evolutionary feedbacks on evolutionary outcomes. Furthermore, the approach can be applied widely to model transient dynamics, and to predict the waiting times and extent of diversification that occurs in a range of contexts [28, 42] . Here it allowed us to predict evolutionary outcomes analytically in the dynamical context of antigenic escape.

Broadly, our results emphasize that epidemiological dynamics may have important implications to the evolution of infectious disease. Here repeated invasions driven by antigenic escape mutants were shown to lead to the long term evolution of more acute virulent pathogens. It has previously been shown that newly invading parasites will temporarily not be at the classic R 0 optima [19,32], but our model shows that this more acute parasite state is a long term optima when there is antigenic escape. Furthermore, this long term optima leads to a more acute severe disease which is clearly of particular importance in human health. Human coronaviruses can evolve antigenically to escape antibody immunity [43] and moreover in principle epidemics of new variants as a novel disease adapts to the host will lead to equivalent dynamics and evolutionary outcomes. It also follows that interventions that impact epidemiological dynamics may also have impacts on the evolution of pathogen traits such as virulence or transmission and it is therefore critical that we have a broader theory of the evolution of virulence that predicts the implications of these epidemiological effects. Furthermore, these dynamical feedbacks are likely to be important in a range of contexts beyond infectious disease and the analytical approaches we use here are likely to be useful in understanding dynamical evolutionary outcomes. 

We consider a model of antigenic escape of a pathogen from host herd immunity on a onedimensional antigenicity space ( ), which describes the changes in the density ( , ) of hosts that are susceptible to antigenicity strain of pathogen at time , and the density ( , ) of hosts that are currently infected and infectious with antigenicity strain of pathogen at time :

where , , and are the transmission rate, virulence (additional mortality due to infection), and recovery rate of pathogens which are independent of antigenicity. ( − ) denotes the degree of cross immunity: a host infected by strain of pathogen acquires perfect cross immunity with probability ( − ) but fails to acquire any cross immunity with probability 1 − ( − ) (this is called polarized cross immunity by Gog and Grenfell 2002) . The cross-immunity kernel ( − ) is assumed to be a decreasing function of the distance | − | between strains and . When a new strain with antigenicity = 0 is introduced at time = 0, the initial host population is assumed to be susceptible to any antigenicity strain of pathogen: (0, ) = 1. = $ # /2 is one half of the mutation variance for the change in antigenicity, representing random mutation in the continuous antigenic space.

We first analyze the dynamics of the primary outbreak of a pathogen and derive the resulting susceptibility profile, which can be viewed as the fitness landscape subsequently experienced by the pathogen. For simplicity we assume that mutation can be ignored during the first epidemic started with the antigenicity stain = 0. The density ' ( ) = ( , 0) of hosts that are susceptible to the currently prevailing antigenicity strain = 0, as well as the density ' ( ) = ( , 0) of hosts that are currently infected by the focal strain change with time as

with ' (0) = 1, ' (0) ≈ 0, and ' (0) = 0. The final size of the primary outbreak,

is determined as the unique positive root of

where ' = /( + ) > 1 is the basic reproductive number [Anderson and May, 1991] . Associated with this epidemiological change, the susceptibility profile * ( ) = ( , ) against antigenicity ( ≠ 0 ) other than the currently circulating strain ( = 0 ) changes by cross immunity as

Integrating both sides of (A7) from = 0 to = ∞, we see that the susceptibility profile ( ) = * (∞) after the primary outbreak at = 0 is

where the last equality follows from (A6). The susceptibility can be effectively reduced by crossimmunity when the primary strain has a large impact (the fraction of hosts remained uninfected, 1 − ' , is small) and when the degree of cross immunity is strong ( ( ) is close to 1). With a strain antigenically very close to the primary strain ( ≈ 0), the cross immunity is very strong ( ( ) ≈ 1) so that the susceptibility against strain is nearly maximally reduced ( ) ≈ 1 − ' . With a strain antigenically distant from the primary strain, ( ) becomes substantially smaller than 1, making the host more susceptible to the strain. For example, if the cross immunity is halved ( ( ) = 0.5) from its maximum value 1, then the susceptibility to that strain is as large as (1 − ' ) '. 6 . If a strain is antigenically very distant from the primary strain, ( ) ≈ 0, and the host is nearly fully susceptibility to the strain ( ( ) ≈ 1). 

Of particular interest is the threshold antigenicity distance 7 that allows the antigenic escape, i.e. any antigenicity strain more distant than this threshold from the primary strain ( > 7 ) can increase when introduced after the primary outbreak. Such a threshold is determined from

With a specific choice of cross-immunity profile,

B High cross immunity (ω=0.6) A Low cross immunity (ω=0.2) the threshold antigenicity beyond which the virus can increase in the susceptibility profile ( ) after the primary outbreak is, from

, and taking logarithm of both sides twice,

Integrating both sides of (A2) over the whole space, we obtain the dynamics for the total density of infected hosts, ̅ ( ) = ∫ ( , ) . As shown in Sasaki and Dieckmann (2011) , the dynamics for viral morph frequency is expressed as

while the dynamics for the within-morph distribution of antigenicity is

From this, the dynamics for the mean antigenicity of a morph,

and the dynamics for the within-morph variance of a morph 

Equations (A18), (A20) and (A22) are general, but they rely on a full knowledge of the dynamics of the susceptibility profile ( , ). In order to make further progress, we use an additional approximation by substituting Eq. (A9), the susceptibility profile over viral antigenicity space after the primary outbreak at = 0 and before the onset of the second outbreak at a distant position. We keep track of two morphs at positions ' ( ) and , ( ), where the first morph is that caused the primary outbreak at = 0, and the second morph is that emerged in the range > 7 beyond the threshold antigenicity 7 defined in Eq. (A9) (and (A11) for a specific form of ( )) as the source of the next outbreak.

Therefore, the frequency, mean antigenicity, and variance of antigenicity of an emerging morph ( = 1) change respectively as

The predicted change in the mean antigenicity is plotted by integrating Eq. (A23). As initial condition, we choose the time when a seed of second peak in the range > 7 first appears, and then compute the mean trait as

In the case of Figure 1 , where = 2, + = 0.6, = 0.001, and = 2, the final size of epidemic for the primary outbreak, defined as (7) was = 0.959, and the critical antigenic distance for the increase of pathogen strain obtained from (A22) was 7 = 2.795. The initial condition for the oligomorphic dynamics (A23) for the second morph was then , ( ' ) = 1.6 × 10 "< , ̅ , ( ' ) = 3.239, , ( ' ) = 0.2675 at ' = 41. In Figure 1 , the predicted trajectory for the mean antigenicity (A24) is plotted as red curve, together with the mean antigenicity change observed in simulation (blue curve).

Here we describe how the initial conditions for oligomorphic dynamics, i.e., the frequency, the mean antigenicity and the variance in antigenicity of the morph that caused the primary outbreak and the morph which may cause the second outbreak are determined. We then show how the accuracy in prediction of the second outbreak depends on the timing for prediction.

We divide the antigenicity space into two at = 7 above which the pathogen can increase under the given susceptibility profile after the primary outbreak, but below which the pathogen cannot increase. We then take relative frequencies of pathogens above 7 and below 7 , and the conditional mean and variance in these separated regions to set the initial frequencies, means, and variances of the morphs at the time ' when we predict the second outbreak:

We then compare the trajectory for mean antigenicity change observed in simulation (blue curve in Figure 1 ) and the predicted trajectory (red curve in Figure 1 ) for mean antigenicity (A24) by integrating oligomorphic dynamics (A23) with initial condition (A25) at time = ' . Figure S1 shows how the accuracy of prediction, measured by the Kullback-Leibler divergence or L2 norm between these two trajectories depends on the timing chosen for the prediction, ' . The second outbreak occurs around t = 54.6 where mean antigenicity jumps from around 0 to around 5. The prediction with OMD is accurate if it is made for ' > 40. Figure 1 is drawn for ' = 41 where the second peak is about to emerge (see Figure S2 ). Even for the latest prediction for ' = 51 in Figure S2 , the morph frequency of the emerging second quasispecies was only 0.3%, so the prediction is still worthwhile to make. Figure S2 shows that the prediction power is roughly constant (albeit with a wiggle) for 5 < ' < 30 (the predicted timings are 10-15% longer than actual timing for 5 < ' < 30), and steadily improved for ' > 30. When the prediction is made very early <( ' < 5) the deviations are larger. 

Let ( ) be the susceptibility of the host population against antigenicity a specific susceptibility profile is given by (A8) in Appendix S1 with cross-immunity function ( ) and the final size ' of epidemic of the primary outbreak. Note that as in appendix S1, the susceptibility profile is in general a function of time. The density ( , ) of hosts infected by a pathogen of antigenicity and virulence changes with time, when rare, as is the fitness of a pathogen with antigenicity and virulence and Š = ∫ ∫ ( , )

Numerical example Figure 5 in the main body shows the oligomorphic dynamics prediction of the emergence of next strain in antigenicity-virulence coevolution. In order to make progress numerically we assume ( ) to be constant in the following analysis because we are interested in the process between the end of the primary outbreak and the emergence of the next antigenicity virulence morph. A The partial differential equations for the density of host ( , ) susceptible to the antigenicity strain at time , and the density of hosts infected by pathogen strain with antigenicity and virulence are The trait space is restricted in a rectangular region: 0 < < max = 300 and min = 0.025 < < 10 = max . Oligomorphic dynamics prediction for the joint evolution of antigenicity and virulence is applied for the next outbreak after the outbreak with the mean antigenicity around = 108 at time = 102 . The susceptibility of host to antigenicity strain at ' = 104.8 after the previous outbreak peaked around time = 102 came to an end is ( ) = ( ' , ). This susceptibility profile remains unchanged until the next outbreak starts, and hence the fitness of a pathogen strain with antigenicity and virulence is given by ( , ) = ( ) ( ) − ( + ) We divide the pathogen quasispecies into two morphs at time ' at the threshold antigenicity 7 above which the net growth rate of pathogen strain under the given susceptibility profile ( ) and the mean antigenicity becomes positive: The initial frequency and the moments of two morphs, the strain 0 with < 7 and the strain 1 with > 7 are then calculated respectively from the joint distribution ( ' , , ) in the restricted region {( , ); 0 < < 7 , min < < max } and that in the restricted region {( , ); 7 < < max , min < < max }. The frequency , of the morph 1 (the frequency of the morph 0 is given by ' = 1 − , ), the mean antigenicity ̅ % and mean virulence Y % of the morph , and variances and covariance, % * , % , % of the morph ( = 0,1) follow (A44)-(A49), where the dynamics for the morph frequency (A44) is simplified in this two morph situation as

with ' ( ) = 1 − , ( ). This is iterated from = ' = 104.8 to 2 = 107. The frequency , of the new morph, the population mean antigenicity ̅ = ' ̅ ' + , ̅ , , virulence Y = ' Y ' + , Y , , variance in antigenicity * = ' ' * + , , * , covariance between antigenicity and virulence = ' ' + , , , and variance in virulence -= ' ' -+ , , are overlayed by red thick curves on the trajectories of moments observed in full dynamics (A51).

In the panel (a) of Figure 5 , the dashed vertical line represents the threshold antigenicity 7 above which ' = ( )/( + Y) > 1 at = 1 = 104.8 where oligomorphic dynamics (OMD) prediction is attempted. Two morphs are then defined according to whether or not the antigenicity exceeds a threshold = 7 : the resident morph (morph 1) is represented as the dense cloud to the left of = 7 and the second morph (morph 2) consisting of all the genotypes to the right of = 7 with their ' greater than one. The within-morph means and variances are then calculated in each region. The relative total densities of infected hosts in the left and right regions defines the initial frequency of two morphs in OMD. A 2D Gaussian distribution is assumed for within-morph trait distributions to have the closed moment equations as explained before. Using these as the initial means, variances, covariances of the two morphs at = 1 . The oligomorphic dynamics for 11 variables (relative frequency of morph 1, mean antigenicity, mean virulence, variances in antigenicity and virulence and their covariance in morph 0 and 1) is integrated up to = 2 . The results are shown in red curves in the panels of second and third rows, which are compared with the simulation results (blue curves).

The panels (c)-(e) in Figure 5 show the change in total infected density (c), mean antigenicity (d), and mean virulence (e). Red curves show the prediction by oligomorphic dynamics from the initial moments of each morph at = 1 to the susceptibility distribution ( ) = ( 1 , ). Red curves in the (d) and (e) show the OMD prediction, which is compared with the simulation results (blue curves). The OMD predicted mean antigenicity, for example, is defined as

where , ( ) is the frequency of morph 1, ̅ ' and ̅ , are the mean antigenicity of morph 0 and 1.

The red curves in the third row of Figure 5 show the OMD-predicted changes in the variance in antigenicity, variance in virulence, and covariance between antigenicity and virulence, which are compared with the simulation results (blue curves). The OMD predicted covariance, for example, is defined as 

Zika virus in the Americas-yet another arbovirus threat

Coronavirus disease 2019 (COVID-19): situation report

2020 COVID-19 and Italy: what next? The Lancet

Dynamics of the 2001 UK foot and mouth epidemic: Stochastic dispersal in a heterogeneous landscape

Prevention of population cycles by parasite removal

DISEASE AND CONSERVATION

Infectious Diseases of Humans

Dynamics and selection of many-strain pathogens

Oligomorphic dynamics for analyzing the quantitative genetics of adaptive speciation

This study was supported by grants R01 GM122061-03 and NSF-DEB-2011109 to MB, ANR JCJC grant ANR-16-CE35-0012-01 to SL and the ESB Cooperation Program, The Graduate University for Advanced Studies, SOKENDAI.

is the mean fitness.Let us decompose the joint frequency distribution ( , ) of the viral quasispecies as (oligomorphic decomposition):where % ( , ) is the joint frequency distribution of antigenicity and virulence in morph (∫ ∫ % = 1) and % is the relative frequency of morph (∑ % % = 1). The frequency of morph then changes aswhere Š % = ∫ ∫ ( , ) % ( , ) is the mean fitness of morph .Assuming that the traits are distributed narrowly around the morph means ̅ % = ∫ ∫ % ( , ) and Y % = ∫ ∫ % ( , ) , so that % = − ̅ % = ( ) and % = − Y % = ( ) where is a small constant, we expand the fitness ( , ) around the means ̅ % and Y % of morph ,

18A-% ( ̅ % , Y % ) are fitness and its first and second derivatives evaluated at the mean traits of morph , and Substituting (A32) into the change in the mean antigenicity of morphwe haveSimilarly, the change in the mean virulence of morph is expressed asEquations (A34)-(A35) from the mean trait change is summarized in a matrix form asis the variance-covariance matrix of the morph .Substituting (A32) into the right-hand side of the change in variance of antigenicity of morph ,and those in the change in the other variance and covariance, we haveIf we assume that antigenicity and virulence within a morph follow two-dimensional Gaussian distribution for given means, variances and covariance, we should have % ( %Eqs. (A39)-(A41) are rewritten in a matrix form aswhereis the Hessian of the fitness function of the morph .In our case (A26) of the joint evolution of antigenicity and virulence of a pathogen after its primary outbreak, the fitness function is given by ( , ) = ( ) ( ) − , and hence % = (, where a prime on ( ) and ( ) denotes differentiation by and , respectively. Substituting these into the dynamics for morph frequencies (A31), for morph means (A34)-(A35), and for within-morph variance and covariance (A39)-(A41), we haveEquations (A44)-(A49) describe the oligomorphic dynamics of the joint evolution of antigenicity and virulence of a pathogen for a given host susceptibility profile ( ) over pathogen antigenicity.Of particular interest is whether the antigenicity evolution or virulence evolution accelerates each other by allowing them evolving simultaneously than when they evolve alone. After the primary outbreak at a given antigenicity, say = 0, the susceptibility ( ) of host population increases due to cross-immunity as the distance > 0 from the antigenicity at the primary outbreak.Combining this with the positive tradeoff between transmission rate and virulence, we see that ( # / ) % = ′( Y % ) ′( ̅ % ) > 0, and then from (A48) we see that the within-morph covariance between antigenicity and virulence stays positive if its initial value is nonnegative:If all second moments are sufficiently small initially for an emerging morph, a quick look at the linearization of (A47)-(A49) around ( % * , % , % -) = (0,0,0) indicates that both % * and % become positive due to the random generation of variance by mutation, * > 0 and -> 0, while the covariance stays close to zero. Then (A50) guarantees that first move of covariance is from zero to positive, which then guarantees that % > 0 for all . Therefore, the second term in (A34) is positive until the mean virulence reaches its optimum ( ′( ) ( ) = 1). This means that joint evolution with virulence accelerates the evolution of antigenicity. The same is true for virulence evolution: the first term in (A35) (which denotes the associated change in virulence due to the selection in antigenicity through genetic covariance between them) is positive, indicating that joint evolution with antigenicity accelerates the virulence evolution.

In this appendix, we show that a pathogen that has the strategy maximizing the growth rate in a fully susceptible population is evolutionarily stable in the presence of antigenic escape.At stationarity, the travelling wave profiles of ž ( ) and ž ( ) along the moving coordinate, = − , that drifts constantly to right with the speed are defined aswith ž (−∞) = ž (∞) = 0, ž (∞) = 1.Let ( , ) be the density of a mutant pathogen strain with virulence ′ and transmission rate ′ that is introduced in the host population where the resident strain is already established (A46). For the initial transient phase in which he density of mutant is sufficiently small, we have an equation for the change in ( , ) = ( , ):with the initial condition (0, ) = ( ), where is a small constant and (⋅) is Dirac function. Here & = 2√ & is the asymptotic wave speed if the mutant strain monopolizes the host population. Therefore, if & < , then < 0, and hence ( , ) for a fixed converges to zero as goes to infinity; which in turn implies that ( , ) converges to zero because ( , ) ≤ ( , ) for all and . Therefore, we conclude that any mutant that has a slower wave speed than the resident can never invade the population, implying that a strain that has the maximum wave speed = 2√ is locally evolutionarily stable.