Jacopo Marchi

Michael Lässig

Aleksandra M Walczak
email: aleksandra.walczak@phys.ens.fr

Thierry Mora
email: thierry.mora@phys.ens.fr

Antigenic waves of virus-immune co-evolution

The evolution of many microbes and pathogens, including circulating viruses such as seasonal influenza, is driven by immune pressure from the host population. In turn, the immune systems of infected populations get updated, chasing viruses even further away. Quantitatively understanding how these dynamics result in observed patterns of rapid pathogen and immune adaptation is instrumental to epidemiological and evolutionary forecasting. Here we present a mathematical theory of co-evolution between immune systems and viruses in a finite-dimensional antigenic space, which describes the cross-reactivity of viral strains and immune systems primed by previous infections. We show the emergence of an antigenic wave that is pushed forward and canalized by cross-reactivity. We obtain analytical results for shape, speed, and angular diffusion of the wave. In particular, we show that viral-immune co-evolution generates a new emergent timescale, the persistence time of the wave's direction in antigenic space, which can be much longer than the coalescence time of the viral population. We compare these dynamics to the observed antigenic turnover of influenza strains, and we discuss how the dimensionality of antigenic space impacts on the predictability of the evolutionary dynamics. Our results provide a concrete and tractable framework to describe pathogen-host co-evolution.

Introduction

The evolution of viral pathogens under the selective pressure of its hosts' immunity is an example of rapid co-evolution.

Viruses adapt in the usual Darwinian sense by evading immunity through antigenic mutations, while immune repertoires adapt by creating memory against previously encountered strains. Some mechanisms of in-host immune evolution, such as the affinity maturation process, are important for the rational design of vaccines. Examples are the seasonal human influenza virus, where vaccine strain selection can be informed by predicting viral evolution in response to collective immunity [START_REF] Dh Morris | Predictive Modeling of Influenza Shows the Promise of Applied Evolutionary 659[END_REF], as well as chronic infections such as HIV [START_REF] Wang | Manipulating the Selection Forces during Affinity Maturation to Generate Cross-Reactive HIV Antibodies[END_REF][START_REF] Jp Barton | Relative rate and location of intra-host HIV evolution to evade cellular immunity are predictable[END_REF][START_REF] Nourmohammad | Host-Pathogen Coevolution and the Emergence of Broadly Neutralizing Antibodies in Chronic Infections[END_REF][START_REF] Nourmohammad | Optimal evolutionary control for artificial selection on molecular phenotypes[END_REF], where coevolution occurs within each host. Because of the relatively short time scales of selection and strain turnover, these dynamics also provide a laboratory for studying evolution and its link to ecology [START_REF] Gandon | Forecasting Epidemiological and Evolutionary Dynamics of Infectious Diseases[END_REF].

It is useful to think of both viral strains and immune protections as living in a common antigenic space [START_REF] Gandon | Forecasting Epidemiological and Evolutionary Dynamics of Infectious Diseases[END_REF], corresponding to an idealized "shape space" of binding motifs between antibodies and their cognate epitopes [START_REF] Segel | Shape space: an approach to the evaluation of cross-reactivity effects, stability and controllability in the immune system[END_REF]. While the space of molecular recognition is high-dimensional, projections onto a low-dimensional effective shape space have provided useful descriptions of the antigenic evolution. In the example of influenza, neutralization data from hemagglutination-inhibition assays can be projected onto a two-dimensional antigenic space [START_REF] Smith | Mapping the Antigenic and Genetic[END_REF][START_REF] Bedford | Integrating influenza antigenic dynamics with molecular evolution[END_REF][START_REF] Jm Fonville | Antibody landscapes after influenza virus infection or vaccination[END_REF]. Mapping historical antigenic evolution in this space suggests a co-evolutionary dynamics pushing the virus away from its past positions, where collective immunity has developed. Importantly, the evolution of influenza involves competitive interactions of antigenically distinct clades in the viral population, generating a "Red Queen" dynamics of pathogen 31 evolution [START_REF] Van Valen | A new evolutionary law[END_REF][START_REF] Yan | Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens[END_REF]. Genomic analysis of influenza data has 32 revealed evolution by clonal interference [START_REF] Strelkowa | Clonal interference in the evolution of influenza[END_REF]; this mode of 33 evolution is well-known from laboratory microbial populations 34 [START_REF] Pj Gerrish | The fate of competing beneficial mutations in an asexual population[END_REF]. In addition, the viral population may split into subtypes. [START_REF] Keiner | Using NFFT 3-a software library for various nonequispaced fast Fourier transforms[END_REF] Such splitting or "speciation" events, which are marked by 36 a decoupling of the corresponding immune interactions, hap-37 pened in the evolution of influenza B [START_REF] Pa Rota | Lineages of Influenza Type B Virus since 1983[END_REF] and of noroviruses 38 [START_REF] White | Evolution of norovirus[END_REF].

39

The joint dynamics of viral strains and the immune systems 40 of the host population can be modeled using agent-based sim-41 ulations [START_REF] Nm Ferguson | Ecological and immunological determinants of influenza evolution[END_REF][START_REF] Bedford | Canalization of the evolutionary trajectory of the human influenza virus[END_REF] that track individual hosts and strains. Such 42 approaches have been used to study the effect of competition 43 on viral genetic diversity [START_REF] Zinder | The Roles of Competition and Mutation in Shaping Antigenic and Genetic Diversity in Influenza[END_REF], to study geographical effects 44 [START_REF] Wen | Explaining the geographical origins of seasonal influenza A (H3N2)[END_REF], and the effect of vaccination [START_REF] Ft Wen | The beneficial effects of vaccination on the evolution of seasonal influenza[END_REF]. Alternatively, systems 45 of coupled differential equations known as Susceptible-Infected-46 Recovered (SIR) models may be adapted to incorporate evolu-47 tionary mechanisms of antigenic adaptation [START_REF] Gandon | Forecasting Epidemiological and Evolutionary Dynamics of Infectious Diseases[END_REF][START_REF] Jr Gog | Dynamics and selection of many-strain pathogens[END_REF][START_REF] Koelle | Understanding the dynamics of rapidly evolving pathogens through modeling the tempo of antigenic change : In fl uenza as a case study[END_REF]. Agent-48 based simulations in 2 dimensions were used to recapitulate 49 the ballistic evolution characteristic of influenza A [START_REF] Bedford | Canalization of the evolutionary trajectory of the human influenza virus[END_REF], and 50 to predict the occurence of splitting and extinction events [START_REF] Marchi | Multi-lineage evolution in viral populations driven by host immune systems[END_REF]. 51 In parallel, theory was developed to study the Red Queen effect 52 [START_REF] Yan | Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens[END_REF][START_REF] Im Rouzine | Antigenic evolution of viruses in host populations[END_REF], based on the well established theory of the traveling 53 fitness wave [START_REF] Im Rouzine | The solitary wave of asexual evolution[END_REF][START_REF] Cohen | Front propagation up a reaction rate gradient[END_REF][START_REF] Desai | Beneficial mutation-selection balance and the effect of linkage on positive selection[END_REF]. While effectively set in one dimension, 54 this class of models can nonetheless predict extinction and 55 splitting events assuming an infinite antigenic genome [START_REF] Yan | Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens[END_REF].

56

In this work, we propose a co-evolutionary theory in an 57 antigenic interaction space of arbitrary dimension d, which is 58 described by joint non-linear stochastic differential equations 59 coupling the population densities of viruses and of protected 60 hosts. We show that these equations admit a d-dimensional 61 antigenic wave solution, and we study its motion, shape, and

62

Significance Statement

Viruses, such as influenza, evolve under the selection of host immune systems. Previously infected individuals become immune, forcing the virus to find susceptible hosts or mutate, chasing it away in antigenic space. We formulate this viral escape process in terms of a low-dimensional wave moving in antigenic space. The dimensionality of the antigenic space impacts both the persistence, as well as stability of viral evolution. We uncover a new characteristic timescale for the persistence of the viral strain, which is an order of magnitude longer than individual host immunity, and emerges collectively from the pressure of the chasing immune systems. These results offer intuition about the antigenic turnover of viruses, and highlight the importance of the effective dimensionality of co-evolution.

D R A F T

stability, using simulations and analytical approximations.

Based on these results, we discuss how canalization and predictability of antigenic evolution depend on the dimensionality d.

Results

A. Coarse-grained model of viral-immune co-evolution. Our model describes the joint temporal evolution of populations of viruses and immune protections in some effective antigenic space of dimension d. Both viral strains and immune protections are labeled by their position x = (x1, . . . , x d ) (or "phenotype") in that common antigenic space (Fig. 1A). In that space, viruses randomly move as a result of antigenic mutations and proliferate through infections of new hosts. Immune memories are added at the past positions of viruses.

Immune memories distributed across the host population provide protection that reduces the effective fitness of the virus.

We coarse-grain that description by summarizing the viral population by a density n(x, t) of hosts infected by a particular viral strain x, and immunity by a density h(x, t) of immune memories specific to strain x in the host population.

At each infection cycle, each host may infect R0 unprotected hosts, where R0 is called the basic reproduction number.

However, a randomly picked host is susceptible to strain x with probability (1 -c(x, t)) M , where c(x, t) is the coverage of strain x by immune memories of the population, and the number M of immune memories carried by each host. Because of cross-reactivity, which allows immune memories to confer protection against closeby strains, immune coverage is given as a function of the density of immune memories:

c(x, t) = 1 M dx h(x , t)H(x -x ), [1]
where H(x -x ) = exp(-|x -x |/r) is a cross-reactivity kernel describing how well memory x protects against strain x, and r is the range of the coverage provided by cross-reactivity. In summary, the effective growth rate, or "fitness", of the virus is given by f

(x, t) ≡ ln[R0(1 -c(x, t)) M ].
The coupled dynamics of viruses and immune memories is then described by the stochastic differential equations (with time in units of infection cycles throughout):

∂tn(x, t) = f (x, t)n(x, t) + D∂ 2 x n + n(x, t)η(x, t) [2] ∂th(x, t) = 1 N h n(x, t) -N (t) h(x, t) M . [3]
Here η is a Gaussian white noise in time and space, η(x, t)η(x , t ) = δ(x -x )δ(t -t ), accounting for demographic noise [START_REF] Hallatschek | The noisy edge of traveling waves[END_REF]. This stochastic term is crucial, as it will drive the evolution of the wave. The diffusion constant D describes the effect of infinitesimal mutations on the phenotype,

D = µ δx 2 1 /2
, where µ is the mean number of mutations per cycle, and δx 2 1 the mean squared effect of each mutation along each antigenic dimension (assuming that mutations do not have a systematic bias, δx1 = 0). The continuous-diffusion assumption implied by Eq. 2 is only valid when there are many small mutation effects, µ 1 and δx r, in constrast with regimes where mutations are rare but have a substantial fitness effect drawn from a distribution [START_REF] Im Rouzine | Antigenic evolution of viruses in host populations[END_REF][START_REF] Bh Good | Distribution of fixed beneficial mutations and the rate of adaptation in asexual populations[END_REF]. Our choice is simpler in that it describes the mutation process through a single parameter D. Along with the choice of the cross-112 reactivity kernel H, it also naturally preserves the isotropy of 113 the antigenic space.

114

The total viral population size, or number of infected hosts, 115 N (t) = dx n(x, t) is subject to fluctuations. At the same 116 time, the host population size N h , remains constant because 117 newly added memories (first term of right-hand side of Eq. 3) 118 overwrite existing ones picked uniformly at random (second 119 term of r.h.s. of Eq. 3). Since each host carries M immune 120 receptors, we have dxh(x, t) = M .

121

If we assume that the system reaches an evolutionary steady 122 state, with stable viral population size N (t) = N , then Eq. 3 123 can be integrated explicitly:

124 h(x, t) = M N t -∞
dt τ e -t-t τ n(x, t ), [START_REF] Nourmohammad | Host-Pathogen Coevolution and the Emergence of Broadly Neutralizing Antibodies in Chronic Infections[END_REF] 125

with τ = M N h /N . This equation shows how the density of 126 protections reflects the past evolution of the viral population. 127 B. Antigenic waves. We simulated Eq. (2)-Eq. (3) on a square 128 lattice (Methods) and found a stable wave solution (Fig. 1B-129 D). The wave has a stable population size N , and moves 130 approximately ballistically through antigenic space, pushed 131 from behind by the immune memories left in the trail of 132 past viral strains (Fig. 1B). These memories exert an immune 133 pressure on the viruses, forming a fitness gradient across the 134 width of the wave (Fig. 1C), favoring the few strains that are 135 furthest from immune memories, at the edge of the wave.

D R A F T

We assume that the solution of the coupled evolution equa-137 tions Eq. ( 2)-Eq. ( 3) takes the form of a moving quasispecies 138 in a d-dimensional antigenic space,

139 n(x, t) = N √ 2πσ 2 exp - (x1 -vt) 2 σ 2 ρ(x2, . . . , x d ). [5]
140

Here, we have written the solution in a co-moving frame, in There are several models of fitness waves that differ in the assumptions on the statistics of mutational effects. Our assumption of diffusive motion makes our projected dynamics equivalent to that studied in ref. [START_REF] Neher | Genealogies of rapidly adapting populations[END_REF], which itself builds on earlier work [START_REF] Cohen | Front propagation up a reaction rate gradient[END_REF]. This equivalence results from the two key assumptions of the mutation model in antigenic space: mutations have a small effect, and their distribution is isotropic, meaning that there are as many deleterious as beneficial mutations. In the limit where the wave is small compared to the adaptation time scale, vτ σ, the wave may be replaced by a Dirac delta function at x = (vt, 0, . . . , 0) in Eq. 4. One can then calculate explicitly the immune density (upstream of the wave) and coverage (downstream of the wave, using Eq. 1):

h(x, t) ≈ M vτ e -vt-x 1 vτ Θ(vt -x1)δ(x2) • • • δ(x d ), [6] c(x, t) ≈ e -(x 1 -vt)/r 1 + vτ /r , x1 ≥ vt, xi>1 r [7]
where Θ(x) = 1 for x ≥ 0 and 0 otherwise. This idealized exponential trail of immune protections h(x, t) corresponds to 175 the blue trace of Fig. 1B, and the coverage or fitness gradient to the isolines of Fig. 1C.

In the moving frame of the wave, (u, x2, . . . , x d ), with u = x1 -vt, the local immune protection and viral fitness can be expanded locally for u, xi vτ (see [START_REF] Im Rouzine | Antigenic evolution of viruses in host populations[END_REF] for a similar treatment in a one-dimensional antigenic space): [START_REF] Smith | Mapping the Antigenic and Genetic[END_REF] where f0 = ln R0 -M ln[1 + r/(vτ )] is the average population fitness, and

f ((u, xi>1); t) ≈ ln R0 1 - e -u/r 1 + vτ /r M ≈ f0 + su,
s = |∂x 1 f | = M r R 1/M 0 -1 [9]
is the fitness gradient. Rescaling the antigenic variable x1 as sx1, this process is equivalent to the evolution of a population where mutation effects are described by diffusion in fitness space with coefficient Ds 2 . This is precisely the model from which the fitness wave solution of Ref. [START_REF] Cohen | Front propagation up a reaction rate gradient[END_REF][START_REF] Neher | Genealogies of rapidly adapting populations[END_REF] was described (see SI Appendix). In the following we will use results from these works to describe the antigenic wave. However, we note that in the usual fitness wave theory, population is kept constant by construction, which implies that fitness is only relevant when compared to the mean of the population. By contrast, in our model population size is itself a dynamical variable, and fitness is defined as an absolute growth rate.

In this version of the model, the fitness of the whole viral population undergoes continuous negative drift due to the constant adaptation of immune systems, encoded in the -svt term in Eq. 8. This negative fitness drift has an analogous effect to subtracting the mean fitness in models with constant population size, making the equivalence possible.

The fitness wave theory allows us to make analytical prediction about the properties of the antigenic wave. Let us start with its population size N , which is regulated by how fast the immune system catches up with the wave. The immune turnover time τ in Eq. 4 is inversely proportional to N : the larger the population size, the faster immune memories are updated, increasing the immune pressure on current viral strains (lower f0), and thus decreasing N . As the moving wave reaches a stable moving state, its size N becomes stable over time, giving the condition (1/N )dN/dt = f0 = 0, which in turn constraints the ratio between the wave's size and speed:

N v = M N h r R 1/M 0 -1 = N h s. [10]
But the fitness wave theory predicts that the speed of the wave itself depends on the population size. The larger N , the more outliers at the nose of the fitness wave, and the further out they may jump in antigenic space, establishing fitter ancestors of the future population. This results in a fitness wave whose speed depends only weakly on population size and mutation rate (see [START_REF] Neher | Genealogies of rapidly adapting populations[END_REF] and SI Appendix),

vF ≈ D 2/3 F 24 ln(N D 1/3 F ) 1/3 , [ 11
]
where DF = s 2 D and vF = sv are the diffusivity and wave speed in fitness space, which are related to their counterparts in antigenic space through the scaling factor s. Replacing this scaling into Eq. 11 yields a relation between antigenic speed and population size, To explain the width σ ⊥ of the wave in the other phenotypic dimensions than that of motion (xi>1), we note that in these directions evolution is neutral. Two strains taken at random in the bulk are expected to have drifted, or 'diffused' in physical language, by an average squared displacement ∆x 2 i = 2DT2 from their common ancestor, so that their mean squared distance is 4D T2 = 2ασ 2 along xi. If one assumes an approximately Gaussian wave of width σ ⊥ , the mean square distance between two random strains along xi should be equal to 2σ 2 ⊥ . Equating the two estimates yields σ 2 ⊥ = ασ 2 . Fig. 2D checks the validity of this prediction against simulations.

v ≈ D 2/3 s 1/3 24 ln(N (Ds 2 ) 1/3 ) 1/3 , [12] D R A F T
Both longitudinal and transversal fluctuations in antigenic space are instances of quantitative traits under interference selection generated by multiple small-effect mutations. The width of these traits is governed by the common relation ∆x 2 i = 2D T2 ∼ σ 2 , which expresses the effective neutrality of the underlying genetic mutations [START_REF] Held | Survival of the simplest in microbial evolution[END_REF]. This relation says that antigenic variations in all dimensions scale in the same way with the model parameters, and the wave should have an approximately spherical shape. Consistently, here we find a wave with a fixed ratio α ≈ 1.66 between transverse and longitudinal variations. This implies a slightly asymmetric shape (which may be non-universal and depend on the microscopic assumptions of our mutation model).

In what parameter regime is our theory valid? The fitness wave theory we built upon is meant to be valid in the large population size, N 1. In addition, we assumed that the fitness landscape was locally linear across the wave. This approximation should be valid all the way up to the tip of wave, given by uc, since this is where the selection of future founder strains happen. This condition translates into uc r, implying D r 2 / ln(N ) 2 (using uc = v 2 /(4Ds) and Eqs. 9,12), where D is in antigenic unit squared per infection cycle. This result means that one infection cycle will not produce enough mutations for the virus to leave the cross-reactivity range. In that limit, another assumption is automatically fulfilled, namely that the width of the wave be small compared to the span of immune memory: σ vτ . Our simulations, which run in the regime of very slow effective diffusion (D/r 2 10 -6 ) and have relatively large population sizes (N 104 ), satisfy these conditions. This explains the good agreement between analytics and numerics.

E. Equations of motion of the wave's position.

The wave solution allows for a simplified picture. The wave travels in the direction of the fitness gradient (or equivalent the gradient of immune coverage) with speed v (Fig. 3A). Occasionally the population splits into two separate waves that then travel away from each other and from their common ancestor (Fig. 3B). The tip of the wave's nose, which contains the high-fitness individual that will seed the future population, determines its future position in antigenic space. In the directions perpendicular to the fitness gradient, this position diffuses neutrally with coefficient D. This motivates us to write effective equations D R A F T D Fig. 3. Stochastic behaviour of the wave: diffusive motion, splits, and extinctions. A. The wave moves forward in antigenic space but is driven by its nose tip, which undergoes antigenic drift (diffusion) in directions perpendicular to its direction of motion. These fluctuations deviate that direction, resulting in effective angular diffusion. B. When antigenic drift is large, the wave may randomly split into subpopulations, creating independent waves going in different directions. Each wave can also go extinct as size fluctuations bring it to 0. C. Cartoon illustrating the wave's angular diffusion. Selection and drift combine to create a inertial random walk of persistence time tpersist. D. Analytical prediction (Eq. 17) for the persistence time, versus estimates from simulations. Symbols and colors are the same as in Fig. 2.

of motion for the mean position of the wave:

dx dt = -v + 2D ξ (t) ∂xc |∂xc| + √ 2Dξ ⊥ (t), [13] c(x, t) = t -∞ dt τ e -t-t τ - |x-x(t )| r , [ 14
]
where ξ and ξ ⊥ are Gaussian white noises in the directions (as illustrated in Fig. 3C). Eqs Eq. ( 13)-Eq. ( 14) simplify to 325 (see SI Appendix):

326 ∂tx ⊥ (t) = +∞ 0 dt T x ⊥ (t) -x ⊥ (t -t ) t e -t /T + √ 2Dξ ⊥ (t), [ 15
]
where T = (v/r+1/τ ) -1 = (r/v)R

-1/M 0 is an effective memory timescale combining the host's actual immune memory, and the cross-reactivity with strains encountered in the past.

Eq. 15 may be solved in Fourier space. Defining x⊥ (ω) = +∞ -∞ dte iωt x ⊥ (t), it becomes:

-iωx ⊥ (ω) 1 + ln(1 -iωT ) iωT = √ 2D ξ⊥ (ω). [16]
To understand the behaviour at long times T , we expand at small ω: -ω 2 x ⊥ (ω) ≈ √ 8D ξ⊥ (ω)/T or equivalently in the temporal domain ∂ 2 t x ⊥ ≈ √ 8Dξ ⊥ (t)/T . This implies that the direction of motion, ê ∼ ∂xf /|∂xf | ∼ ∂tx/|∂tx|, undergoes effective angular diffusion in the long run: ∂tê = √ 8Dξ ⊥ (t)/(vT ). The persistence time of that inertial motion,

tpersist = v 2 T 2 4D = r 2 4D R -2/M 0 , [ 17
]
does not depend explicitly on speed, population size, or the dimension of antigenic space. However, a larger diffusivity implies larger N and v while reducing the persistence time.

Likewise, a larger reproduction number R0 or smaller memory capacity M speeds up the wave and increases its size, but also reduces its persistence time. This implies that, for a fixed number of hosts N h , larger epidemic waves not only move faster across antigenic space, but also change course faster.

This persistence time scales as the time it would take a single virus drifting neutrally to escape the cross-reactivity range, r 2 /D. For comparison, the much shorter timescale for a population of viruses to escape from the cross-reactivity range r,

tescape = r v = T R 1/M 0 = N h M N (R 1/M 0
-1), [START_REF] Bedford | Canalization of the evolutionary trajectory of the human influenza virus[END_REF] scales with the inverse incidence rate N h /N . This is consistent with the whole population having been infected at least one every ∼ N h /N infection cycles. This separation of time scales is consistent with the observation that evolution in the transverse directions is driven by neutral drift, which is much slower than adaptive evolution in the longitudinal direction.

Both tpersist and tescape are longer than the coalescence time of the viral population, uc/v ∼ σ 2 /4D, since they reflect long-term memory from the immune system. However, while tescape ∼ N h /N is related to the re-infection period and is thus bounded by the hosts' immune memory (itself bounded by their lifetime, which we do not consider), tpersist can be longer than that. This is possible thanks to inertial effects, which are allowed by the high-order dynamics of Eq. 15 generated by the immune system. This very much like when, in mechanics, a massive object set in motion in a given direction will keep that direction without the need for an external force to maintain it.

The high-frequency behaviour of Eq. ( 16) has a logarithmic divergence, meaning that the total power of ê is infinite unless we impose a (ultraviolet) cutoff. Such a regularization emerges from the fine structure of the wave. While the motion of the wave is driven by its nose tip, the immune pressure only extends back to the recent past of the bulk of the distribution, which stands at a distance uc away from the nose. In other words, there is a lag (and thus an gap uc in antigenic space) between the most innovative variants that drive viral evolution, and the majority of currently circulating variants which drive host immunity. Mathematically, this implies that the domain

D R A F T

of integration of the first term in the right-hand side of Eq. ( 15) should start at tc = uc/v, which regularizes the divergence.

A more careful analysis provided in the SI Appendix shows that this regularization does not affect the long-term diffusive behaviour of the wave. G. Canalization, speciations, and predictability of antigenic evolution. We now examine how deflections of the wave in the transverse direction determines the predictability and stability of the viral quasi-species. Assuming t T , angular diffusion causes motion to be deflected as (see SI Appendix)

x 2 ⊥ = 8(d-1)D 3T 2 t 3 .
Crucially, this deflection depends on the dimension of the antigenic space, because the displacement acts additively in each of the transversal coordinates. Higher dimension means more deviation from the predictable course of the wave, and thus less predictability. We can define a predictability time scale

t predict ∼ [8(d -1)/3] -1/3 T 2/3 (r 2 /D) 1/3 , [19]
which is the time it takes for prediction errors to become of the order of the cross-reactivity range. In low dimensions, this time scales as a weighted geometric mean between tescape ∼ T and tpersist ∼ r 2 /D. However, at high dimensions t predict may be significantly reduced, causing loss of predictability even below tescape. The prediction timescale is distinct from the previously discussed persistence time: t predict involves the integrated displacement in the transversal direction, while tpersist quantifies the diffusion of the tangent velocity vector. Thus, t predict may be interpreted as quantifying the predictability of the actual location of the next viral population in antigenic space, while tpersist gives the predictability of the general direction of evolution, which changes more slowly. Therefore, the persistence time is both harder to extract from data and less relevant for actionable predictions.

To get a sense of numbers, we can compare our results with epidemiological data, taking the evolution of influenza as an Large deflections may also cause speciations, or splits, which occur when two substrains co-exist long enough to become independent from the immune standpoint. This happens when two sub-lineages see the difference of their transverse positions ∆x ⊥ become larger than ∆x0 ∼ r, within some limited period given by the coalescence time. We estimated the rate of such splitting events using a saddle-point approximation (see SI Appendix): [START_REF] Wen | Explaining the geographical origins of seasonal influenza A (H3N2)[END_REF] with α some numerical factor. Simulations confirmed the validity of this scaling (Fig. 4a). The splitting rate grows with the dimension (Fig. 4b), 441 consistent with the intuition that departure from canalized 442 evolution is easier when more directions of escape are available. 443 Splitting events are expected to strongly affect our ability to 444 predict the future course of the wave. However, the rarity 445 of such events (exponential scaling of k split ) means that they 446 will have a lower impact on predictability than deflections. 447 These results provide a theoretical and quantitative basis from 448 which to assess the effect of dimension on predictability, and 449 possibly estimate d from antigenic time course data of real 450 viral populations.

k split ≈ 3 8 v 2 4D e -L , L = α s 3 R -2/M 0 D 2 r 4 (d -1)v 5 1/4

Discussion

452

In this work, we have developed an analytical theory for 453 studying antigenic waves of viral evolution in response to 454 immune pressure. We showed that predictabilty is limited 455 by two features of antigenic evolution, transversal diffusion 456 and lineage speciations of the antigenic wave, both of which 457 explicitly depend on the dimensionality of antigenic space.

458

To derive these results, we explicitly embedded the anti-459 genic phenotype in a d-dimensional Euclidean space. This 460 description is different from previous work that considered one-461 [START_REF] Im Rouzine | Antigenic evolution of viruses in host populations[END_REF] or infinite-dimensional antigenic spaces [START_REF] Yan | Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens[END_REF]. It allows 462 for the possibility of compensatory mutations, and makes it 463 easier to compare results with empirical studies of viral evo-464 lution projected onto low-dimensional spaces [START_REF] Smith | Mapping the Antigenic and Genetic[END_REF][START_REF] Bedford | Integrating influenza antigenic dynamics with molecular evolution[END_REF]. Unlike 465 these studies, however, our work does not address the question 466 how an effective dimension of antigenic space arises from the 467 molecular architecture of immune interactions. Rather, we 468 focused on the implications of the dimensionality of antigenic 469 space for phenotypic evolution and its predictability.

470

Our results suggest a hierarchy of time scales for viral 471 evolution. The shortest is the coalescence time T2 , which 472 determines population turnover. Then comes tescape, which 473 is the time it takes the viral population to escape immunity 474 elicited at a previous time point. The longest timescale is the 475 persistence time tpersist, which governs the angular diffusion 476 of the wave's direction, but has no bearing on the prediction 477 of the actual position of the dominant strain in antigenic 478 D R A F T space. That time scale is due to inertial effects. It does not rely directly on the hosts' immune memories, and may thus exceed their individual lifetimes. Finally, the prediction timescale t predict , beyond which prediction accuracy falls below the resolution of cross-reactivity, scales between tescape and tpersist at low dimensions. This time scales measures the predictability of transversal fluctuations, and is thus the most relevant for actual predictions of future dominant strains in antigenic space. Importantly, it decreases with the dimension of the antigenic space, and may become arbitrarily low at very high dimensions. The fact that the evolution of influenza strains are hard to predict beyond a year suggests that the effective dimension may indeed be large.

Our solution builds on the fitness wave solution for a diffusion model of mutation effects [START_REF] Cohen | Front propagation up a reaction rate gradient[END_REF][START_REF] Neher | Genealogies of rapidly adapting populations[END_REF]. It implies a particular dependence of the wave's speed on the population size, Eq. 12.

General distribution of non-infinitesimal mutational effects, such as considered in [START_REF] Bh Good | Distribution of fixed beneficial mutations and the rate of adaptation in asexual populations[END_REF], would yield different expressions for the speed. However, we expect most of our other results to hold-in particular, all expression that do not carry an explicit logarithmic dependence on N , as well the effective equations of motion for the wave. Our results strongly rely on the assumption of a homogeneous, isotropic antigenic space.

We expect our results to be affected by anisotropies (e.g. in the mutational or the cross-reactivity Kernels), or by structure in the intrinsic fitness landscape (i.e. not linked to immunity). Such structure may funnel the wave in preferred directions, hinder it, or favor its splitting. Generally, the local geometry and metric of the space is expected to determine the evolutionary behaviour. For instance, Yan et al. ( 12) assumed a Hamming distance metric in an effectively infinite antigenic space, meaning that any mutation is both an escape mutation and a candidate for a lineage split. By contrast, in our geometry, escape happens only in the direction of the wave, while splits originate from mutations perpendicular to that direction, due to the choice of a Euclidean metric. While our results emphasize the role of the effective dimension d, studying other geometrical effects is an interesting topic for future work.

Despite these caveats, it is interesting to ask whether the effective antigenic dimension d can be extracted from data. A possible scheme for doing so starts by inferring the effective model parameters. R0 may be estimated from exponential epidemic growth in a susceptible population. Dependence of key quantities on M such as s is weak. M may be assumed to be of the order of the number of antigenically distinct infections encountered during a host's lifetime, ∼ 4-6 (every 15 years). D/r 2 may be inferred from v/r, which can be estimated from cross-immunity assays or from the incidence rate N/N h . Alternatively, since D/r 2 is the inverse time it takes for mutations to neutrally evade immunity, it could be estimated directly from genomic data by computing the time for unselected mutations (whose rate is inferred from synonymous mutations) to affect antigenic sites. Interestingly, if v/r and D/r 2 can be inferred independently, predictions about the wave's shape, width, angular diffusion and splitting do not depend on the particular choice of fitness wave theory.

Assuming that all these parameters are known, the splitting rate, which depends sensitively on d (Fig. 4b), could be used to infer an effective dimension. Since splitting is rare and may not be observed in practice, one could define instead partial splits, where a sublineage diverges an antigenic distance ∆x0 < r from the main lineage, for which the same scaling as Eq. 20 540 holds (see SI Appendix). Alternatively, our results could 541 be used to check the consistency of dimensionality-reduction 542 schemes based on serological assays [START_REF] Smith | Mapping the Antigenic and Genetic[END_REF][START_REF] Bedford | Integrating influenza antigenic dynamics with molecular evolution[END_REF][START_REF] Jm Fonville | Antibody landscapes after influenza virus infection or vaccination[END_REF], by testing our 543 predicted relations between the speed of the wave, its width 544 and length, and angular diffusion properties, and ask what 545 choice of dimension best agrees with our theory.

546

Our framework should be applicable to general host-547 pathogens systems. For instance, co-evolution between viral 548 phages and bacteria protected by the CRISPR-Cas system 549 [START_REF] Er Westra | The ecology and evolution of microbial CRISPR-Cas adaptive immune systems[END_REF] is governed by the same principles of escape and adap-550 tation as vertebrate immunity. Even more generally, our 551 theory (Eqs. 2,3) may be relevant to the coupled dynamics 552 of predators and preys interacting in space (geographical or 553 phenotypic), opening potential avenues for experimental tests 554 of these theories in synthetic microbial systems. Given the 555 current context of the global SARS-CoV-2 pandemic, it is 556 natural to ask whether our results could be applicable to pre-557 dict its evolution. While our theory describes the long-term 558 co-evolution of viral strains with the hosts' immune systems, 559 in which most hosts have been exposed to at least on strain of 560 the virus, SARS-CoV-2 is still in a phase of growth, and has 561 not exhausted the reservoir of susceptible hosts. As the situa-562 tion develops, it will be interesting to see whether its future 563 evolution follows a Red Queen type of evolution like influenza, 564 goes extinct, or splits into many antigenically independent sub-565 lineages. While our model may shed light on these questions, 566 fine microscopic details such as geographical and population 567 structure impose additional challenges for predictions.

Methods

569

We simulated discrete population dynamics of infected hosts 570 n(x, t)) and immune protections n h (x, t) ≡ N h h(x, t) (all inte-571 gers) on a 2D square lattice with lattice size ∆x ranging from 572 10 -5 r to 0.1r. Each time step corresponds to a single infection 573 cycle, ∆t = 1. At each time step: (1) viral fitness f is com-574 puted at each occupied lattice site from the immune coverage 575 Eq. 1; (2) viruses at each occupied lattice site are grown accord-576 ing to their fitness, n(x, t + 1) ∼ Poisson[(1 + f ∆t)n(x, t)]; (3) 577 viruses are mutated by jumping to nearby sites on the lattice; 578 (4) the immune system is updated according to a discrete ver-579 sion of Eq. 3, by implementing n h (x, t + 1) = n h (x, t) + n(x, t) 580 and then removing N (t) protections at random (so that N h 581 remains constant).

582

To implement (1), we used a combination of exact compu-583 tation of Eq. 1 and approximate methods, including one based 584 on non-homogeneous fast Fourier transforms [START_REF] Keiner | Using NFFT 3-a software library for various nonequispaced fast Fourier transforms[END_REF][START_REF] Potts | Fast convolution with radial kernels at nonequispaced knots[END_REF]. Details 585 are given in the SI Appendix.

586

To implement (3), we drew the number of mutants at each 587 occupied site from a binomial distribution Binomial(n(x, t), 1-588 e -µ∆t ). The number of new mutations m affecting each 589 of these mutants is drawn from a Poisson distribution of 590 mean µ∆t conditioned on having at least 1 mutation. The 591 new location of each mutant is drawn as x + δx, with 592 δx = round( m i=1 i) (rounding is applied to each dimen-593 sion), where i is a vector of random orientation and modulus 594 drawn from a Gamma distribution of mean δ ∼ 2∆x and shape 595 parameter 20. This distribution was chosen so as to maximize 596 the number of non-zero jumps while maintaining isotropy. We 597 then define D = µ δx 2 1 /2.

D R A F T

was initialized as a Gaussian distribution centered at (0, 0)

Fig. 1 .

 1 Fig. 1. A simple model of viral-host co-evolution predicts the emergence of an antigenic wave. A. Schematic of the co-evolution model. Viruses proliferate while effectively diffusing in antigenic space (here in 2 dimensions) through mutations, with coefficient D. Past virus positions are replaced by immune protections (light blue). Immune protections create a fitness gradient for the viruses (green gradient) favoring strains at the front. Both populations of viruses and immune populations are coarse-grained into densities in antigenic space. B. Snapshot of a numerical simulation of Eq. 2-3 showing the existence of a wave solution. The blue colormap represents the density of immune protections h(x, t) left behind by past viral strains. The current virus density n(x) is shown in red. C. Close-up onto the viral population, showing fitness isolines. The wave moves in the direction of the fitness gradient (arrow) through the enhanced growth of stains at the edge of the wave (black dots). D. Distribution of fitness across the viral population (corresponding to the projection of B. along the fitness gradient). Parameters for B-D: D/r 2 = 3•10 -9 , N h = 10 8 , ln R0 = 3, M = 1.

141 which a motion

  with constant speed v takes place in the di-142 rection of the coordinate x1, and fluctuations in the other 143 dimensions, ρ(x2, . . . , x d , t), centered around xi = 0 for i > 1, 144 are assumed to be independent. In the next sections, we will 145 analyse solutions of this form. First, we will project the d-146 dimensional antigenic wave onto the one-dimensional fitness 147 space; this projection produces a travelling fitness wave (26-148 28, 31, 32) that determines the antigenic speed v and the mean 149 pair coalescence time T2 of the viral genealogy. Second, we 150 will study the shape of the d-dimensional quasispecies and 151 determine the fluctuations in the transverse directions. These 152 fluctuations produce a key result of this paper: immune in-153 teractions canalize the evolution of the antigenic wave; this 154 constraint can be quantified by characteristic time scales gov-155 erning the transverse antigenic fluctuations. Canalization is 156 most pronounced in spaces of low dimensionality d and, as we 157 discuss below, affects the predictability of antigenic evolution. 158 C. Speed of antigenic evolution. Projected onto the fitness 159 axis f = f (x, t), the solution is approximately Gaussian 160 (Fig. 1D). This representation suggests a strong similarity to 161 the fitness wave solution found in models of rapidly adapting 162 populations with an infinite reservoir of beneficial mutations 163 (26-28, 31, 32). To make the analogy rigorous, we must assume 164 that the fitness gradient in antigenic space is approximately 165 constant, meaning that fitness isolines are straight and equidis-166 tant. Mutations along the gradient direction have a fitness 167 effect that is linear in the displacement, while mutations along 168 perpendicular directions are neutral and can be treated inde-169 pendently. Note that while we will use this projection onto fitness to compute the speed of the antigenic wave, the under-171 lying antigenic wave remains in d dimensions; we will come 172 back to transverse fluctuations in the next sections.

DFig. 2 .

 2 Fig.2. Analytical prediction of wave properties. Shown are the numerical versus analytical predictions for the wave's population size N (A), speed v (B), width σ along the wave's direction of motion (C), and width σ ⊥ in the direction perpendicular to motion (D), with d = 2 dimensions. Length are in units of the crossreactivity range (so that r = 1, with no loss of generality). Parameters: N h = 10 8 (squares), 10 10 (circles), or 10 12 (triangles); ln R0 = 1 (filled symbols) or 3 (empty symbols); M = 1 (small symbols) or 5 (large symbols).

F. Angular diffusion and persistence of the antigenic wave. 314

 314 305 along, and perpendicular to, the fitness gradient ∂xf /|∂xf | = 306 -∂xc/|∂xc|. D is an effective diffusivity in the direction of 307 motion resulting from the fluctuations at the nose tip. These 308 fluctuations are different than suggested by D, as they involve 309 feedback mechanisms between the wave's speed v, size N , and 310 advancement of the fitness nose uc. In the following, we do 311 not consider these fluctuations, and focus on perpendicular 312 fluctuations instead. 313 In the description of Eqs. 13-14, the viral wave is pushed by 315 immune protections left in its trail. The fitness gradient, and 316 thus the direction of motion, points in the direction that is set 317 by the wave's own path. This creates an inertial effect that 318 stabilizes forward motion. On the other hand, fluctuations in 319 perpendicular directions are expected to deviate the course of 320 that motion, contributing to effective angular diffusion. To 321 study this behaviour, we assume that motion is approximately 322 straight in direction x1 = vt, and study small fluctuations in 323 the perpendicular directions, x ⊥ = (x2, . . . , x d ), with |x ⊥ | r 324

  example, with an infection cycle time of 3 days. It is assumed that individuals lose immunity to the circulating strain of the flu within ∼ 5 years ∼ 500 cycles, meaning that the wave would travel a distance r in t = 500, i.e. v/r ∼ 2 • 10 -3 . For instance, with N h = 10 9 -10 10 , R0 = 2, and M = 1, we may choose D/r 2 = 3 • 10 -6 to get a speed of the same order, v/r ∼ 1.3 • 10 -3 , and tpersist ∼ 2 • 10 4 ∼ 200 years. By contrast, the predictability timescale t predict is much shorter and depends on dimension, albeit slowly, ranging from ∼ 20 years for d = 2 to about 2 years for d = 1000. We stress that these numbers are obtained by scaling laws, and should not be taken as precise quantitative predictions.

Fig. 4 .

 4 Fig. 4. Rate of speciation. A. Rescaled rate of splitting events, defined as the emergence of two substrains at distance ∆x0 = 0.1r from each other in antigenic space, meaning that they are becoming antigenically independent. The predicted scaling, k split ∼ (v 2 /D)e -L , as well as the definition of the collective variable L as a function of the model parameters, are given by Eq. 20. The line shows a linear fit of the logarithm of the ordinate. B. Predicted rate of splitting as a function of the dimension d, for R0 = 2, M = 1, N h = 10 9 , and D/r 2 = 3 • 10 -6 , with ∆x0 = r.

  451

  568

| Marchi et al.

with size N and width σ in all dimensions, to which 0.1% 601 additional viruses are randomly added within the interval 602 (0; uc) along x1 (N , σ, and uc being all given by the theory 603 prediction). Immune protections are placed according to Eq. 6.