Vaccination shapes evolutionary trajectories of SARS-CoV-2
Matthijs Meijers${}^{\,\rm a}$, Denis Ruchnewitz${}^{\,\rm a}$, Marta Łukszab,
Michael Lässig${}^{\,\rm a,\ast}$
${}^{\,\rm a}$ Institute for Biological Physics, University of Cologne,
Zülpicherstr. 77, 50937 Köln, Germany
${}^{\,\rm b}$ Tisch Cancer Institute, Departments of Oncological Sciences and
Genetics and Genomic Sciences, Icahn School of Medicine at Mount Sinai, New
York, NY, USA
∗ To whom correspondence should be addressed. Email<EMAIL_ADDRESS>
## Abstract

The large-scale evolution of the SARS-CoV-2 virus has been marked by rapid
turnover of genetic clades. New variants show intrinsic changes, notably
increased transmissibility, as well as antigenic changes that reduce the
cross-immunity induced by previous infections or vaccinations[1, 2, 3, 4]. How
this functional variation shapes the global evolutionary dynamics has remained
unclear. Here we show that selection induced by vaccination impacts on the
recent antigenic evolution of SARS-CoV-2; other relevant forces include
intrinsic selection and antigenic selection induced by previous infections. We
obtain these results from a fitness model with intrinsic and antigenic fitness
components. To infer model parameters, we combine time-resolved sequence
data[5], epidemiological records[6, 7], and cross-neutralisation assays[9, 8,
10]. This model accurately captures the large-scale evolutionary dynamics of
SARS-CoV-2 in multiple geographical regions. In particular, it quantifies how
recent vaccinations and infections affect the speed of frequency shifts
between viral variants. Our results show that timely neutralisation data can
be harvested to identify hotspots of antigenic selection and to predict the
impact of vaccination on viral evolution.

## Introduction

Two classes of molecular adaptation have been observed in the evolution of
SARS-CoV-2 to date. Multiple mutations carry intrinsic changes of viral
functions, such as increasing the binding affinity to human receptors[1], the
efficiency of cell entry[2, 3], or the stability of viral proteins[11, 12].
Other mutations, referred to as antigenic changes, decrease the neutralizing
activity of human antibodies[4, 8, 9, 10], thereby reducing the immune
protection against secondary infections[13, 14]. The strains that inherit a
given mutation define a clade of the evolving viral population. Several of
these molecular changes had drastic evolutionary and epidemiological impact,
inducing global turnover of viral clades and concurrent waves of the pandemic.
Over the last two years, three genetic variants and their associated clades
successively gained global prevalence: Alpha ($\alpha$) from March to June in
2021, Delta ($\delta$) from June to December in 2021 and Omicron ($o$) in
2022. These were named Variants of Concern (VOCs) by the World Health
Organization[15]; other VOCs gained temporary regional prevalence. Several
studies reported fitness advantages of VOCs inferred from epidemiological
trajectories and comparative functional studies[16, 17, 18, 3, 19].
Importantly, however, the evolutionary impact of antigenic changes is time-
dependent, because it depends on previously acquired population immunity: a
larger amount of previous infections or vaccinations increases the global
fitness advantage of an antigenic escape mutation. Specifically, multi-strain
epidemiological models and simulations suggest that vaccinations can favour
the emergence of escape variants[20, 21, 22, 23] and influence the turnover of
circulating clades[24, 25]; effects of this kind have been reported for some
clades of human influenza[26]. In the case of SARS-CoV-2, pandemic infection
and massive vaccination programs, with a global count of 4.5 billion
vaccinations and $>$200 million confirmed cases in 2021[6], have built up
partial population immunity, but its feedback on viral evolution has not been
quantified. This leads to the central question of this paper: what is the
impact of vaccination and infection rates on the turnover of SARS-CoV-2
clades? To address this question, we infer a data-driven fitness model for
SARS-CoV-2 variants with distinct components of intrinsic fitness and
antigenic fitness by vaccination and infection.

## Results

### Trajectories and speed of clade turnover

As a first step, we map the evolutionary trajectories of the three global
clade shifts in the last two years. To track circulating clades, we analyse a
set of $>$5M quality-controlled SARS-CoV-2 sequences obtained from the GISAID
database[5]. We assign these sequences to genetic clades using a standard set
of amino acid changes[27]; then we infer time-dependent clade frequencies from
strain counts smoothened over a period of $\sim$30 days (Methods). To obtain
accurate, time-resolved data, we record frequency trajectories at the level of
regions (countries and US states). Including all regions satisfying uniform
criteria of data availability (Methods), we obtain frequency trajectories of
the ${\sl 1}-\alpha$, the $\alpha-\delta$, and the $\delta-o$ shift for 11,
16, and 14 regions, respectively. Here, 1 denotes the set of clades
circulating prior to $\alpha$, including the wild type (wt) and the early 614G
mutation in the spike protein. Fig. 1a shows trajectories of the ancestral and
the invading clade for the $\alpha-\delta$ and $\delta-o$ shifts in Italy;
trajectories for all regions of this study are reported in Fig. S1 and S2.

Assuming that large-scale frequency shifts of viral clades are adaptive
processes, we can infer the underlying selective force from frequency
trajectories. Specifically, the fitness difference (selection coefficient)
between invading and ancestral clades takes the form

$\hat{s}(t)=\frac{\rm d}{{\rm d}t}\log\frac{x_{\rm inv}(t)}{x_{\rm anc}(t)},$
(1)

where $x_{\rm inv}(t)$ and $x_{\rm anc}(t)$ are the corresponding frequencies
(here and below, empirical selection coefficients inferred from frequency
trajectories are marked by a hat). We note that this relation is independent
of other co-circulating clades (Methods). In Fig. 1b, we show time-resolved,
regional selection coefficients of the invading and ancestral clade for the
$\alpha-\delta$ and $\delta-o$ shifts. These data reveal two opposing trends:
During the $\alpha-\delta$ shift, selection increases with time in 16 of 16
regions. Conversely, selection driving the $\delta-o$ shift decreases with
time in 12 of 14 regions. Compared to a reference of time-independent
selection, the $\alpha-\delta$ shift runs at an accelerating speed, the
$\delta-o$ shift at a decelerating speed. The time dependence of selection is
statistically significant ($P<10^{-15}$ for $\alpha-\delta$, $P<10^{-5}$ for
$\delta-o$; two-sided Wald test). In contrast, the earlier ${\sl 1}-\alpha$
shift does not show a significant signal of time-dependent selection
($P>0.01$, Fig. S3). In what follows, we will relate this pattern to feedback
of vaccination on viral evolution.

Fig. 1: Evolutionary, epidemiological, and immune tracking of SARS-CoV-2.
Time-dependent trajectories are shown for the clade shift from $\alpha$ to
$\delta$ (left column) and from $\delta$ to $o$ (right column). (a) Observed
frequency trajectories of relevant clades, $x_{i}(t)$, for the clade shifts in
Italy. (b) Empirical selection coefficient (fitness difference) between
invading and ancestral clade, $\hat{s}(t)$, for all regions. Selection
trajectories are derived from the frequency trajectories of (a) and plotted
against time counted from the midpoint. Summary statistics: cross-region
linear regression (black solid line), cross-region average (black diamond),
and rms cross-region variation of selection (black dashed line). (c)
Population immunity functions of the ancestral and invading variant, $C_{\rm
anc}^{k}(t)$ and $C_{\rm inv}^{k}(t)$, in relevant immune channels $k$ for
Italy (coloured lines). Cross-immunity differences in a given channel, $C_{\rm
inv}^{k}(t)-C_{\rm anc}^{k}(t)$, are highlighted by shading (colours indicate
which variant receives a fitness advantage). See Figs. S1–S3 for tracking of
all shifts in all regions and reporting of rms statistical errors.

### Cross-immunity trajectories

Cross-immunity induced by a primary infection against subsequent infections by
related pathogens is routinely tested by neutralisation assays, which measure
the minimum antiserum concentration required to neutralise the second antigen.
Relative, inverse concentrations are reported as serum dilution titers; here
we use logarithmic titer values, $T$ (with base 2). For SARS-CoV-2, recent
work[9, 8, 10, 28, 29] has established a matrix of titers, $T_{i}^{k}$,
measuring neutralisation of variant $i$ in immune channel $k$ (Fig. 2a, Table
S1). Here and below, immune channels label primary challenges inducing
specific antisera, including infections by different variants
($k=\alpha,\delta,o,\dots$), as well as primary and booster vaccinations
($k={\rm vac},{\rm bst}$; titers shown here are for mRNA vaccines). Together,
these data provide a first, coarse-grained cross-immunity landscape of SARS-
CoV-2. Infection-induced cross-immunity titers are maximal when primary
infection and secondary challenge are by the same variant (Fig. 2a).
Similarly, titers induced by primary vaccination, $T_{i}^{\rm vac}$, are
maximal against strains from the ancestral clade, which contains the strain
used for vaccination[30]. Differences of neutralisation titers, $\Delta
T_{ij}^{k}=T_{i}^{k}-T_{j}^{k}$, measure differences in functional antibody
binding between strains of different variants; evolved titer reductions are
also referred to as antigenic advance. Notably, each of the global clade
shifts observed to date, ${\sl 1}-\alpha$, $\alpha-\delta$, and $\delta-o$,
has decreased neutralisation by vaccination, i.e., generated antigenic
advance, $\Delta T_{{\sl 1}\alpha}^{\rm vac},\Delta T_{\alpha\delta}^{\rm
vac},\Delta T_{\delta o}^{\rm vac}>0$ (Fig. 2b). Moreover, the in-vivo
concentration of neutralising antibodies decays exponentially with time after
immunisation[31, 32]. This translates into a linear titer reduction,
$T_{i}^{k}(\Delta t)=T_{i}^{k}-\Delta t/\tau$, with an estimated half life
$\tau\sim 65$ days (Methods).

Fig. 2: Cross-neutralisation and antigenic fitness. (a) Neutralisation titers,
$T_{i}^{k}$, of human antisera induced by different primary challenges
(columns: infection by strains from clade $\alpha$, $\delta$, $o$, mRNA
primary and booster vaccination) assayed against different test strains (rows:
strains from clade $\alpha$, $\delta$, $o$); see Table S1. (b) Strain tree of
SARS-CoV-2 with lineages $i$ colored by vaccine neutralisation titers,
$T_{i}^{\rm vac}$. WHO Variants of Concern are marked by bars. The ancestral
clade (${\sl 1}$) has the highest neutralisation (yellow); the successive
clade shifts ${\sl 1}-\alpha$, $\alpha-\delta$, and $\delta-o$ decrease
neutralisation, i.e., induce antigenic advance (see text). (c) The cross-
immunity, $c_{i}^{k}$, induced by a primary immunisation in channel $k$ (top:
infection by $\alpha$, bottom: infection by $\delta$) against a secondary
infection of clade $i$ (blue dot: $\alpha$, red dot: $\delta$) is a Hill
function of the neutralising titer, $T_{i}^{k}$ (ref. [13, 14], Methods).
Cross-immunity decreases with increasing antigenic advance $\Delta T_{ik}^{k}$
(bars). (d) Cross-immunity decays with time after primary immunisation,
$\Delta t$ (in units of the characteristic decay time $\tau$; ref. [31, 32]).
According to the fitness model, cross-immunity induces a proportional
antigenic fitness cost; the resulting time-dependent selection coefficient
(fitness difference) between clades is marked by shading.

Importantly, recent work for SARS-CoV-2 has also shown that neutralisation
titers predict the cross-immunity $c_{i}^{k}$, defined as the relative drop of
secondary infections in human cohorts. Specifically, $c_{i}^{k}=H(T_{i}^{k})$
is a Hill function[13, 14] (Fig. 2c, details are given in Methods), consistent
with the underlying biophysics of antibody-antigen binding and with results
for other viral pathogens[33, 34, 35, 36]. The post-immunisation decay of
antibody concentration induces a decay of cross-immunity, $c_{i}^{k}(\Delta
t)=H(T_{i}^{k}-\Delta t/\tau)$. Together, cross-immunity depends in a
predictable, nonlinear way on neutralisation titer and on time since primary
immunisation. Fig. 2 shows two examples of this pattern. Primary infection by
an $\alpha$ strain induces a high cross-immunity against other $\alpha$
strains and a reduced cross-immunity against $\delta$ strains
($c_{\alpha}^{\alpha}>c_{\delta}^{\alpha}$) (Fig. 2c, top). Both factors
decrease by antibody decay; their difference has a maximum at an intermediate
time since primary infection (Fig. 2d, top). Infection by a $\delta$ strain
induces cross-immunity factors of opposite ranking
($c_{\delta}^{\alpha}<c_{\delta}^{\delta}$) and similar decay (Fig. 2cd,
bottom).

To track population immunity over time, we combine these cross-immunity
factors with infection and vaccination data. In each region, we record
cumulative fractions of immunised individuals, $y_{k}(t)$, in each channel $k$
(clade-specific infections, primary and booster vaccinations; see Figs. S1 and
S2). Their derivatives $\dot{y}_{k}(t)$ are the rates of new immunisations in
channel $k$. Clade-specific infection data are obtained by multiplying the
total rate of new infections reported in each region with the simultaneous
viral clade frequencies $x_{k}(t)$ (Fig. 1a). In the regions included in our
analysis, vaccination has been predominantly by mRNA vaccines (Methods). By
weighting with the time-dependent cross-immunity factors $c_{i}^{k}(\Delta
t)$, we infer the population cross-immunity against clade $i$ by immunisation
in channel $k$,

${C}_{i}^{k}(t)=\int^{t}c_{i}^{k}(t-t^{\prime})\,\dot{y}_{k}(t^{\prime})\,dt^{\prime}.$
(2)

In Fig. 1c, we plot the cross-immunity trajectories relevant for the
$\alpha-\delta$ and $\delta-o$ shifts in Italy; trajectories for all regions
are reported in Figs. S1 and S2. The $\alpha-\delta$ shift shows sizeable and
increasing immunity induced by primary vaccination, while infection-induced
immunity remains small. During the $\delta-o$ shift, immunity by primary
vaccination declines, while booster- and infection-induced immunity components
increase. During the earlier ${\sl 1}-\alpha$ shift, population cross-immunity
is still small in most regions. We conclude that the joint dynamics of new
immunisations and antibody decay can produce complex and opposing cross-
immunity patterns.

### Inference of intrinsic and antigenic selection

To quantify the feedback of cross-immunity on viral evolution, we use a
minimal, computable fitness model,

$f_{i}(t)=f_{i}^{0}-\sum_{k}\gamma_{k}C_{i}^{k}(t),$ (3)

where $f_{i}(t)$ is the absolute fitness, or epidemic growth rate, of a viral
strain. Fitness is proportional to the log of the effective reproductive
number, $f_{i}(t)=\tau_{0}^{-1}\log R_{i}(t)$, where $\tau_{0}$ denotes the
infectious period (Methods). Here, we write fitness as the sum of a time-
independent intrinsic component, $f_{i}^{0}$, and of time-dependent antigenic
components, $f_{i}^{k}(t)=-\gamma_{k}C_{i}^{k}(t)$ (Methods). Each component
is proportional to the corresponding cross-immunity factor $C_{i}^{k}(t)$ with
a weight factor $\gamma_{k}$ for each immune channel $k$. Hence, selection is
generated by cross-immunity differences between competing strains (shading in
Fig. 1c and Fig. 2c). This type of fitness model has been established for
predictive evolutionary analysis of human influenza[24, 37, 38] and is
grounded in multi-strain epidemiological models[39]. The minimal fitness model
does not account for differences in cross-immunity between human hosts (for
example, through differences in immunodominance[40]) and for correlations
between multiple prior infections (antigenic sin[41]).

For SARS-CoV-2, we compute fitness at the level of variants, neglecting
fitness differences between strains within a clade. Similarly, we evaluate
cross-immunity at the level of variant-specific prior infection and of primary
and booster vaccination, using the trajectories $C_{i}^{k}(t)$ calculated
above (Fig. 1c). To compare model and data, we compute the fitness difference
between invading and ancestral strain for each regional trajectory:
$s(t)=f_{\rm inv}(t)-f_{\rm anc}(t)=s_{0}+s_{\rm ag}(t)$, where intrinsic
selection, $s_{0}$, and antigenic selection, $s_{\rm ag}(t)=\sum_{k}s_{k}(t)$,
are given by equation (3). Then we decompose the model-based selection
trajectories into mean and change, $s(t)=\langle s\rangle+\Delta s(t)$
(brackets denote time averages over the trajectory for a given region). The
empirical trajectories $\hat{s}(t)$ are decomposed in the same way (Fig. 1b).
Cross-region selection differences, measured by the rms deviation of
$\langle\hat{s}\rangle$, reflect inhomogeneous conditions of contact
limitations, surveillance, geography, and population structure that are not
included in the minimal model. In a given region, however, variants compete
under more homogeneous conditions. Therefore, we infer antigenic selection
from the regional selection change, $\Delta\hat{s}(t)$. We use a minimal model
with just 3 antigenic parameters: a uniform $\gamma_{\rm vac}$ for vaccination
and boosting (downweighted by a factor $a$ in the $\delta-o$ shift to account
for double infections[42]) and a uniform $\gamma_{k}=b\gamma_{\rm vac}$ for
all infection channels $k$ (upweighted by a factor $b$ to correct for relative
underreporting; Methods). We infer maximum-likelihood (ML) values of these
parameters by calibrating computed and empirical trajectories, $\Delta s(t)$
and $\Delta\hat{s}(t)$, for the $\alpha-\delta$ and $\delta-o$ shifts. The
intrinsic selection coefficients, $s_{0}$, are then obtained as the time-
independent part of selection. Details of the inference procedure are given in
Methods; ML model parameters and selection coefficients for all clade shifts
are reported in Table S2 and S3.

Fig. 3: Antigenic and intrinsic selection drive SARS-CoV-2 evolution. We
compare empirical selection trajectories and model predictions for the clade
shift from $\alpha$ to $\delta$ (left column) and from $\delta$ to $o$ (right
column). (a) Empirical selection change, $\Delta\hat{s}$, obtained from the
selection trajectories of Fig. 1a are plotted against model predictions,
$\Delta s$, for all regions. Rms statistical errors are reported in Figs.
S1–S2. (b) Breakdown of fitness model components. Intrinsic selection
coefficients (black) and antigenic selection coefficients in marked immune
channels (coloured), as inferred from the ML fitness model (bars: region- and
time-averaged value for each crossover; arrows: region-averaged rms temporal
change, $\langle(\Delta s)^{2}\rangle^{1/2}$, with marked direction;
confidence intervals are given in Table S3).

In Fig. 3a, we plot $\Delta s$ from the ML fitness model against the
corresponding empirical selection change, $\Delta\hat{s}$, obtained from the
trajectories of Fig. 1b. We obtain a remarkable data compression: for most
regions, the antigenic fitness computed from equation (3) reproduces the
empirical fitness changes (intrinsic selection drops out of this comparison).
This can be further quantified: the covariance between data and ML model,
$\langle\Delta s\Delta\hat{s}\rangle$, explains $\sim 50$% of the empirical
variance of selection, $\langle(\Delta\hat{s})^{2}\rangle$; this level of
covariance is found on average and in most individual regions. A detailed
comparison of data and model trajectories, $\Delta\hat{s}(t)$ and $\Delta
s(t)$, for all regions is shown in Figs. S1 and S2. As a control, the model
predicts only small selection change for the ${\sl 1}-\alpha$ shift,
consistent with the weak time dependence of the empirical selection
trajectories (Fig. S3). We conclude that time-dependent cross-immunity
explains the time-dependence of selection governing SARS-CoV-2 variant shifts.

### Impact of vaccination and infection on evolution

From the ML fitness model, we obtain a breakdown of intrinsic and antigenic
selection components relevant for each clade shift. Intrinsic selection is
strong and positive in all three major clade shifts, with average selection
coefficients $s_{0}=0.05-0.08$, consistent with strong functional differences
observed between the $\alpha$, $\delta$, and $o$ variants[3, 43] (Fig. 3b,
Table S3). Antigenic selection becomes equally strong in the $\alpha-\delta$
and $\delta-o$ shifts. Its two main components, vaccination- and infection-
induced selection, are statistically significant parts of the fitness model,
partial models with only one component have a strongly reduced posterior
likelihood (differences in model complexity are accounted for by a Bayesian
information criterion; see Methods and Table S2).

Vaccination induces cross-immunity differences between variants, resulting in
positive antigenic selection of average strength $s_{\rm vac}=0.04$ in the
$\alpha-\delta$ shift and $s_{\rm vac}=0.06$ in the $\delta-o$ shift (Fig. 3b,
Table S3). These selection coefficients quantify the evolutionary impact of
primary SARS-CoV-2 vaccination: they measure the relative increase in
effective reproduction number of the invading variant by partial escape from
vaccination-induced immunity ($\tau_{0}s_{\rm vac}=R_{\rm inv}/R_{\rm
anc}-1$). Vaccination-induced antigenicity also explains the observed time-
dependence of selection (Fig. 1b, Figs. S1 and S2): $s_{\rm vac}$ increases
during the $\alpha-\delta$ shift because of increasing vaccination levels, but
decreases during the $\delta-o$ shift because vaccination-induced immunity
fades. In both shifts, primary vaccination generates the dominant components
of antigenic selection (Fig. 3b). Booster vaccinations have increased breadth;
they induce higher neutralisation $T_{\delta}^{\rm bst},T_{o}^{\rm bst}$ and
reduced antigenic advance $\Delta T_{\alpha\delta}^{\rm bst},\Delta T_{\delta
o}^{\rm bst}$ compared to primary vaccinations[44, 10, 9, 45] (Fig. 2a, Table
S1). Hence, booster vaccinations generate higher cross-immunity but weaker
selection for antigenic escape (Fig. 1c, Fig. 4ab). The net effect of boosters
in the $\delta-o$ shift is opposite to that of primary vaccinations: we infer
a negative selection coefficient $s_{\rm bst}=-0.01$. This is because boosters
remove cross-immunity differences and antigenic selection generated by the
preceding primary vaccination (Fig. 3b).

Infection-induced antigenic selection increased in net strength from 0.01 in
the $\alpha-\delta$ shift to 0.03 in the $\delta-o$ shift. Notably, it always
contains components of opposite sign: primary infections by the ancestral
clade generate positive selection, while infections by the invading clade
generate negative selection. This frequency-dependent negative feedback acts
to prolong the coexistence of ancestral and invading clade. Together,
antigenic selection can produce complex but computable patterns of time
dependence.

These results require careful interpretation. They show that vaccination and
previous infections induced sizeable antigenic selection on circulating SARS-
CoV-2 variants and modulated the speed of successive clade shifts. However,
antigenic selection did not cause or prevent any of these shifts, because
intrinsic functional changes generated sizeable fitness advantages of the
invading variants independently of population immunity. The breakdown of
selection given in Fig. 3b applies to the set of regions accessible to our
analysis; the relative weights of vaccination- and infection-induced selection
components are expected to be different in other regions. The availability of
comparable data precludes a fully global model-based analysis. An additional,
model-free inference of selection in regions with low vaccination coverage is
given in Methods. Most importantly, the fitness model and our data analysis do
not predict any simple relation between vaccination coverage and speed of
evolution. This is because cross-immunity channels are correlated: fewer
vaccinations lead to more infections, generating buildup of cross-immunity in
other channels and complex long-term effects.

Fig. 4: Antigenic landscapes and immune weights generate selection hotspots.
(a) Antigenic landscapes. Cross-immunity factors of major and recent clades
(colored dots) are plotted against neutralisation titers in different immune
channels: (a) primary vaccination, (b) booster vaccination, (c) $o$-induced
infection. (d) Long-term immune weight trajectories of different channels,
$Q_{k}(t)$. (e, f) Long-term fitness trajectories of major and recent clades.
Periods of clade shifts are highlighted by shading. (e) Antigenic selection
between successive variants, $s_{\rm ag}(t)$. Model-based trajectories for
each variant pair are shown up to the end of the corresponding clade shift.
(f) Time-dependent fitness gap, $\delta f(t)$ (Methods). Model-based
trajectories for each variant $i$ (lines) are shown from in the time interval
where $0.01\leq x_{i}(t)\leq 0.99$; empirical selection is marked by dots.
Selection hotspots: when sizeable cross-immunity drop on the flank of an
antigenic landscape (arrows in a-c) coincides with large immune weights
(arrows in d), the fitness model predicts time windows of strong selection for
antigenic escape (arrows in e, f mark clade shifts starting in a selection
hotspot). Immune weight and fitness trajectories are averaged over regions
(see Fig. S4 for regional trajectories).

### Fitness trajectories and selection hotspots

The ML fitness model can be applied to the long-term turnover of viral clades
up to date, including recent frequency changes between the variants BA.1,
BA.2, and BA.4/5 within the $o$ clade (we use shorthands $o1$, $o2$, and
$o45$). First, we look at two building blocks of antigenic fitness: antigenic
landscapes and immune weights. We define antigenic landscapes for each immune
channel $k$ by plotting all cross-immunity factors $c_{i}^{k}$ against their
corresponding titers $T_{i}^{k}$ (using a fixed time delay $\Delta t=\tau$ to
account for antibody decay in an approximate way; Methods). These landscapes
visualise antigenic drift, that is, the partial escape from population
immunity by gradual evolutionary steps[46, 47] (Fig. 4a-c; arrows mark
sizeable steps between successive variants). In this picture, the time-
dependence of cross-immunity is captured by immune weight functions
$Q_{k}(t)$, which measure recent infections or vaccinations in channel $k$,
again over a time window of order $\tau$ (Fig. 4d, Methods). Next, we
juxtapose these immune trajectories to long-term trajectories showing the
antigenic selection between successive variants, $s_{\rm ag}(t)$ (Fig. 4e),
and the fitness gap of each variant, $\delta f_{i}(t)=f_{i}(t)-\bar{f}(t)$
(Fig. 4f); these trajectories are computed from equation (3). Fitness gaps are
shifted by the mean population fitness $\bar{f}(t)=\sum_{i}x_{i}(t)f_{i}(t)$
and include the intrinsic component (Methods). As expected from the analysis
above, the ML model is in quantitative agreement with empirical selection
(dots in Fig. 4f). The trajectories of Fig. 4ef are averaged over 14 regions
(for regional trajectories, see Fig. S4).

Together, the trajectories of Fig. 4 show a pattern of selection hotspots. The
fitness model predicts time windows of strong antigenic selection when
antigenic advance on the flank of a Hill landscape generates sizeable cross-
immunity loss and coincides with high immune weight. We now trace this pattern
through successive clade shifts. At early stages of evolution, until spring
2021, all immune weights were small (Fig. 4d). Hence, intrinsic selection
governed the ${\sl 1}-\alpha$ shift; consistently, the antigenic advance
(neutralisation titer drop) $\Delta T_{{\sl 1}\alpha}^{{\rm vac}}$ was small
(Fig. 4a). Between spring 2021 and spring 2022, primary and booster
vaccination generated the dominant immune weights and induced directional
antigenic selection for escape from the vaccine strain, while infection-
induced immunity remained relatively small (Fig. 4d). The clade shifts
$\alpha-\delta$ and $\delta-o$ carried increasing antigenic advance $\Delta
T_{\alpha\delta}^{{\rm vac}}$ and $\Delta T_{\delta o}^{{\rm vac}}$, and
smaller advance with respect to boosting (Fig. 4ab). These changes mark the
onset of antigenic drift. The fitness model identifies a first clear hotspot
of antigenic selection driving the $\delta-o$ shift. At the start of this
shift, large cross-immunity change and large immune weight coincided in the
primary vaccination channel (orange arrows in Fig. 4ade). This is consistent
with the observed dynamics: $\delta-o$ was faster than the previous shifts and
under exceptionally high initial selection, $\hat{s}=0.15$ (orange arrow in
Fig. 4f, Fig. 1b). The following shift, $o1-o2$, involved subclades with
similar neutralisation by primary and booster vaccination. This shift is
inferred to be governed predominantly by intrinsic selection (Table S1 and
S3); consistently, selection is only weakly time-dependent (Fig. S3).

The most recent viral-immune co-evolution shows two important novelties. In
spring 2022, infection immunity increased, while vaccination-induced immunity
decreased; both components are reaching comparable weights (Fig. 4d). Whether
vaccination remains at sizeable immune weight will depend on availability and
acceptance of vaccines in the future. These immune weight changes also mark
the onset of immune drift, that is, the response of population immunity to
antigenic drift of the viral population. Recently, population immunity has
shifted its center of mass from wt towards $o$; components cognate to each of
these clades have reached comparable weight (red arrows in Fig. 4d).

At this point, our fitness model predicts the next selection hotspot for novel
variants carrying antigenic advance away from vaccination _and_ from $o$
infections. The recent antigenic evolution within the $o$ clade follows this
scenario: the emerging variants $o4$ and $o5$ (BA.4 and BA.5) combine
antigenic advance in three channels[48, 45] (red arrows in Fig. 4bce).
Consistently, these variants show fast initial growth with empirical selection
$\hat{s}=0.12$ (red arrow in Fig. 4f). Moreover, near-future mutations
carrying antigenic advance in the same channels are predicted to be in the
same hotspot (dashed arrows in Fig. 4bc).

If these emerging variants develop into major clade shifts, they will further
increase the $o$ immune weight factors. On the other hand, population immunity
against earlier variants could be maintained by backboosting[49] of $o$
infections or by bivalent vaccines with a wt component. The future
evolutionary trajectories of new variants will also depend on their mutual
antigenic relations, which have not yet been assayed comprehensively to date.
Together, the recent evolutionary dynamics signals the unfolding of antigenic
complexity towards coexistence of multiple antigenic variants and immune
classes.

## Discussion

Here we have established a data-driven, multi-component fitness model for the
evolution of SARS-CoV-2. By applying this model to recent evolutionary
trajectories in multiple regions, we have quantified intrinsic and antigenic
selection driving the genetic and functional evolution of the virus. In
particular, primary vaccination impacted on the speed of global clade shifts
in 2021. Booster vaccination generated higher cross-protection, but weaker
selection for antigenic escape in the same period (Fig. 3). These results
underscore that vaccine breadth is important for constraining antigenic escape
evolution. More broadly, they highlight the need to integrate evolutionary
feedback into vaccine design.

In the recent evolution of SARS-CoV-2, two general trends are revealed by our
analysis. Antigenic selection has increased in strength and has broadened its
target: primary infection by distinct viral variants has generated an
increasing number of antigenic selection components (Fig. 3b, Fig. 4). These
trends mark the transition from initial, post-zoonotic adaptation of the virus
to evolution to an endemic state, where antigenic evolution continues to be
fuelled by the buildup of population immunity to circulating viral variants. A
plausible end point of this transition becomes clear by comparison with
influenza, a long-term endemic virus in humans. In influenza, the viral escape
from population immunity follows a specific mode of antigenic drift, where
multiple variants with different cross-immunity profiles compete for
prevalence[50]. This mode is marked by continuous, adaptive clade turnover
with characteristic time scales of several months, which is substantially
slower than the recent prevalence shifts of SARS-CoV-2 variants (Fig. 1a). In
contrast, non-antigenic mutations in influenza proteins are under broad
negative selection; observed changes often compensate the deleterious
collateral effects of antigenic evolution on conserved molecular traits
(including protein stability and receptor binding)[51, 50, 52]. If SARS-CoV-2
reaches a similar endemic state, antigenic evolution is expected to slow down
and most intrinsic changes, e.g., in binding affinity to human receptors, will
become compensatory. Recent findings of compensatory evolution leading to
Omicron support this scanario[53].

The expected transition of SARS-CoV-2 to gradual, multi-faceted antigenic
evolution will open the possibility to predict the future evolution of the
viral population by data-driven fitness models[24, 37, 38] and to inform
preemptive vaccination strategies[54]. Previous work has established an
important prerequisite of predictions: neutralisation assays of human antisera
against viral strains quantify the immune protection of human cohorts against
secondary infections[14, 13]. Here, we have shown that this data can be
harvested at the population scale, to compute immune drift and inform
antigenic fitness models. As a first step of short-term predictions, we have
identified emerging variants in antigenic selection hotspots, in quantitative
agreement with their observed clade growth (Fig. 4). This and future
predictions of SARS-CoV-2 evolution require integrated analysis of genome
sequences, epidemiological records, and increasingly complex antigenic data.
While sequence and epidemiological data are already collected in large
amounts, our analysis calls for world-wide, real-time tracking of antigenic
evolution by cross-neutralisation assays. This will be critical for our
ability to predict antigenic escape evolution and to integrate such
predictions into vaccine design.

## Methods

### Sequence data and primary sequence analysis.

The study is based on sequence data from the GISAID EpiCov database[5]
available until 06-22-2022. For quality control, we truncate the 3’ and 5’
regions of sequences and remove sequences that contain more than 5% ambigous
sites or have an incomplete collection date. We align all sequences against a
reference isolate from GenBank[55] (MN908947), using MAFFT v7.490[56]. Then we
map sequences to Variants of Concern/Interest (VOCs/VOIs), using the set of
identifier amino acid changes given in Outbreak.info[27]. As a cross-check, we
independently infer a maximum-likelihood (ML) strain tree from quality-
controlled sequences under the nucleotide substitution model GTR+G of
IQTree[57], using the reference isolate hCoV-19/Wuhan/Hu-1/2019 (GISAID-
Accession: EPI ISL 402125) as root. For assessment of the tree topology, we
use the ultrafast bootstrap function[58] with 1000 replicates. Internal nodes
are timed by TreeTime[59] with a fixed clock rate of $8\times 10^{-4}$ under a
skyline coalescent tree prior[60]. Consistently, variants are mapped to unique
genetic clades (subtrees) of the ML tree (Fig. 2b).

### Frequency trajectories of variants.

For a given variant $i$, we define the smoothened count
$n_{i}(t)=Z^{-1}\sum_{\nu\in i}\exp[-(t-t_{\nu})^{4}/\delta^{4}]$, where the
sum runs over all sequences $\nu$ mapped to variant $i$, $t_{\nu}$ is the
collection date of sequence $\nu$, and $Z$ is a normalisation constant. We use
a smoothening period $\delta=33$d. The corresponding variant frequency is then
defined by normalisation over all co-existing variants,
$x_{i}(t)=n_{i}(t)/\sum_{i}n_{i}(t)$. These frequency trajectories, evaluated
separately for each region of this study, are shown in Fig. 1 and Figs. S1,
S2.

### Inference of empirical selection.

In a population of different variants, the absolute fitness of each variant is
defined as the growth rate of its population,

$f_{i}(t)=\frac{\dot{N}_{i}(t)}{N_{i}(t)}$ (4)

($i=1,\dots,n$). The absolute fitness is related to its reproductive number,
defined as the mean number of new infections generated by an individual during
its infectious period $\tau_{0}$,

$R_{i}(t)=\exp[\tau_{0}f_{i}(t)].$ (5)

The fitness difference (selection coefficient) between a given pair of
variants, $s_{ij}(t)=f_{i}(t)-f_{j}(t)$, is given by $s_{ij}(t)=({\rm d}/{\rm
dt})(\log N_{i}(t)-\log N_{j}(t))=({\rm d}/{\rm dt})\log((N_{i}(t)/N_{j}(t))$,
independently of the other co-circulating variants. This relation can also be
written in terms of the population frequencies
$x_{i}(t)=N_{i}(t)/\sum_{k}N_{k}(t)$, leading to equation (1) of the main
text.

For each clade shift and each region included in the study, we infer a
trajectory of empirical selection,

$\hat{\bf s}=(\hat{s}(t_{1}),\hat{s}(t_{2}),\dots,\hat{s}(t_{n})),$ (6)

which records the time-dependent fitness difference between invading and
ancestral strain, $\hat{s}(t_{i})=f_{\rm inv}(t_{i})-f_{\rm anc}(t_{i})$
($i=1,\dots,n$). A hat distinguishes these empirical selection coefficients
from their model-based counterparts introduced below. At each point of the
trajectory, we evaluate the selection gradient of equation (1),

$\hat{s}(t_{i})=\frac{1}{\Delta t}\left[\log\left(\frac{x_{\rm
inv}(t_{i}+\Delta t/2)}{x_{\rm anc}(t_{i}+\Delta
t/2)}\right)-\log\left(\frac{x_{\rm inv}(t_{i}-\Delta t/2)}{x_{\rm
anc}(t_{i}-\Delta t/2)}\right)\right]\qquad(i=1,\dots,n),$ (7)

using a time window $\Delta t=30$d for the ${\sl 1}-\alpha$ and $\delta-o$
shifts and $\Delta t=40$d for the $\alpha-\delta$ shift (which extends over a
longer period). Increasing $\Delta t$ reduces the statistical error of
$\hat{s}(t_{i})$ but reduces the time span covered by a trajectory
$\mathbf{\hat{s}}$. We evaluate equation (6) for the maximal time interval
such that $x_{\rm anc}(t_{i}\pm\Delta t/2)>0.01$ and $x_{\rm
inv}(t_{i}\pm\Delta t/2)>0.01$ along the entire trajectory. The start point
$t_{1}$ is the first day when $x_{\rm inv}(t-\Delta t/2)>0.01$. From this
point, selection is recorded weekly, $t_{i}-t_{i-1}=7$d ($i=2,\dots,n$), and
$t_{n}$ is the last point of this sequence where $x_{\rm anc}(t+\Delta
t/2)>0.01$. Single measurements $\hat{s}(t_{i})$ are excluded when at least
one of the sequence counts $n_{\rm anc}(t_{i}\pm\Delta t/2)$ or $n_{\rm
inv}(t_{i}\pm\Delta t/2)$ is $<10$. Statistical errors for selection
trajectories are evaluated by binomial sampling of counts $n_{\rm anc}(t)$ and
$n_{\rm inv}(t)$ with a pseudocount of 1. Empirical selection trajectories are
reported in Fig. 1b and Figs. S1-S3.

For the subsequent analysis, we grade the complete clade shifts ${\sl
1}-\alpha$, $\alpha-\delta$, $\delta-o1$, $o1-o2$ by the time dependence of
their empirical selection trajectories (Fig. S3). We evaluate two summary
statistics: (i) the amount of systematic time-dependent variation of
selection, defined as $\mbox{Var}(s_{\rm lin})$, averaged over regions, where
$s_{\rm lin}(t)$ is a linear regression to the ensemble of trajectories; (ii)
the statistical significance of the linear regression, $P$ (two-sided Wald
test). This identifies two shifts with substantial, statistically significant
time-dependent variation of selection, $\alpha-\delta$ and $\delta-o$.

### Infection and vaccination trajectories.

Daily vaccination and infection rates for individual regions have been
obtained from Ourworldindata.org[6] and from CDC COVID Data Tracker[7] for US
states (download date: 06-22-2022). Clade-specific infection rates $y_{k}(t)$
are computed by multiplying the total daily infection rates reported in each
region with the simultaneous viral clade frequencies $x_{k}(t)$. The resulting
cumulative population fractions of infected individuals, $y_{k}(t)$, together
with cumulative population fractions of primary and booster vaccinations,
$y_{\rm vac}(t)$ and $y_{\rm bst}(t)$, are reported in Figs. S1-S2.

### Data integration for regional analysis.

This study is based on sequence data and epidemiological data from multiple
regions (countries and US states) for parallel analysis. Sequence data is used
to infer empirical selection trajectories for individual clade shifts, as
defined in equations (6) and (7). Epidemiological records provide input to the
antigenic fitness model, equations (2) and (3). Evaluation of the fitness
model, which is detailed below, integrates data of both categories and
requires stringent criteria of data availability and comparability.

To enable this analysis, we choose the set of countries to be included in
model inference based on uniform criteria. Additionally, we include 3 US
states (New York, Texas, California), each representative of a different
geographic region, that satisfy the same criteria. For each clade shift ${\rm
anc}-{\rm inv}$, we require the following: (i) anc and inv are majority
variants at times $t$ and $t^{\prime}>t$ of the clade shift, respectively;
i.e., $x_{\rm anc}(t)>0.5$ and $x_{\rm inv}(t)>0.5$. This criterion excludes
regions where other variants are prevalent during the shift ${\rm anc}-{\rm
inv}$ (e.g., Brazil and South Africa have $x_{\alpha}<0.5$ throughout the
${\sl 1}-\alpha$ shift). (ii) anc and inv have a combined, smoothened sequence
count $n_{\rm anc}(t)+n_{\rm inv}(t)>n_{0}$ throughout the clade shift. This
criterion ensures that the empirical frequencies $x_{\rm anc}(t)$ and $x_{\rm
inv}(t)$, especially minority frequencies, can be estimated with reasonable
statistical errors. We use threshold values $n_{0}=500$ for ${\sl 1}-\alpha$
and $\alpha-\delta$ and $n_{0}=750$ for $\delta-o$ (reflecting the increased
sequence availability). (iii) The empirical selection trajectory
$\hat{\mathbf{s}}$ contains at least 4 (${\sl 1}-\alpha$, $\delta-o$) or 6
($\alpha-\delta$) measured points $\hat{s}(t_{i})$; the threshold values
reflect the relative duration of shifts. This criterion ensures a sufficient
signal-to-noise ratio for inference of temporal variation along the
trajectory. (iv) In the $\delta-o$ shift, the cumulative fraction of $o$
infections exceeds a threshold value, $y_{o}>0.01$. The $o$ variant, which is
characterised by many less severe cases, is likely to be particularly affected
by underreporting. This criterion excludes regions with very low $o$ count
($y_{o}$ is less than $\sim$ 20% of the remaining regions) and ensures that
cross-immunity trajectories, as given by equation (2), can be evaluated across
regions with sufficient consistency. (v) Vaccinations have been predominantly
by mRNA vaccines and epidemiological records in the database[6] are complete.
This criterion ensures that antigenic data for mRNA vaccines can be used
uniformly (Table S1). It excludes regions with substantial use of viral vector
vaccines (e.g., the UK) and with partial records (e.g., for booster
vaccinations in Sweden and Croatia).

Based on these criteria, our analysis includes (i) 11 regions for the ${\sl
1}-\alpha$ shift, (ii) 16 regions for the $\alpha-\delta$ shift, and (iii) 14
regions for the $\delta-o$ shift (Fig. 1, Fig. 3, Figs. S1-S3). Regions
analysed for both $\alpha-\delta$ and $\delta-o$ are used for the long-term
trajectories (Fig. 4, Fig. S4). Scope and limitations of this set of regions
for the inference of selection are described below.

### Antigenic data.

Neutralisation assays for SARS-CoV-2 test the potency of antisera induced by a
given primary immunisation to neutralise viruses of different variants. Log
dilution titers measure the minimum antiserum concentration required for
neutralisation,

$T_{i}^{k}=\log_{2}\frac{K_{0}}{K_{i}^{k}},$ (8)

relative to a reference concentration $K_{0}$. Hence, $\log_{2}$ titer
differences, or neutralisation fold changes, $\Delta T_{ij}^{k}\equiv\Delta
T_{i}^{k}-T_{j}^{k}$, measure differences in antigenicity between variants,
$\Delta T_{ij}^{k}=\log_{2}(K_{j}^{k}/K_{i}^{k})$. We note that these
differences are specific to each primary challenge (immune channel) $k$. For
example, the inequality $T_{\alpha\delta}^{\rm bst}<T_{\alpha\delta}^{\rm
vac}$ reflects the increased breadth of booster vaccinations compared to
primary vaccinations. In contrast, uni-valued antigenic distances between
variants, $d_{ij}$, can be computed from the titer matrix $(T_{i}^{k})$ by
multi-dimensional scaling methods[61, 62]. Such distance measures average over
inhomogeneities between immune channels.

Here we define a matrix of titer drops $\Delta T_{i}^{k}$,

$\Delta T_{i}^{k}=T_{*}^{k}-T_{i}^{k}\qquad\mbox{($i=\alpha$, $\delta$, $o$
($o1$), $o2$, $o45$; \; $k=\alpha$, $\delta$, $o$ ($o1$), $o2$, $o45$, vac,
bst)},$ (9)

with respect to a reference for each immune channel, $T_{*}^{k}=T_{k}^{k}$
($k=\alpha,\delta,o(o1),o2,o45$) and $T_{*}^{k}=T_{\sl 1}^{k}$ ($k={\rm
vac},{\rm bst}$). This procedure eliminates technical differences between
assays in absolute antibody concentration. We assemble this matrix in Table
S1, using primary data from different sources[63, 72, 65, 66, 64, 67, 68, 3,
69, 8, 75, 70, 71, 44, 9, 28, 29, 74, 76, 77, 78, 45, 73, 48]. We proceed as
follows: (i) For matrix elements with available data, $\Delta T_{i}^{k}$ is
the average of the corresponding primary measurements. This procedure
eliminates technical differences between assays in absolute antibody
concentration. As appropriate for the analysis in our set of regions, all
vaccination titers refer to mRNA vaccines. (ii) If no data are available for
$\Delta T_{i}^{k}$ but the conjugate titer $\Delta T_{k}^{i}$ has been
measured, we use the approximate substitution $\Delta T_{i}^{k}\approx\Delta
T_{k}^{i}$, as discussed in ref. [79]. (iii) If no data are available for
$\Delta T_{i}^{k}$ but the titer $\Delta T_{j}^{k}$ of a closely related clade
has been measured, we use the approximate substitution $\Delta
T_{i}^{k}\approx\Delta T_{j}^{k}$, which should be understood as a lower bound
(this applies to the recent variants $o2$ and $o45$).

The matrix of absolute neutralisation titers, $T_{i}^{k}$, is then computed by
equation (9), combining the titer drops $\Delta T_{i}^{k}$ of Table S1 and the
reference titers $T_{*}^{k}=6.5$, ($k=\alpha,\delta,o(o1),o2,o45$),
$T_{*}^{\rm vac}=7.8$, $T_{*}^{\rm bst}=9.8$ reported in ref. [30]. A titer
difference between vaccination and booster, $T_{*}^{\rm bst}-T_{*}^{\rm
vac}\approx 2.0$, has been observed in several studies[9, 44, 10]. The titers
$T_{i}^{k}$ enter the cross-immunity functions $c_{i}^{k}$, $C_{i}^{k}(t)$,
and $\bar{c}_{i}^{k}$ defined below and are shown in Fig. 2ab.

The decay of antibody concentration after primary immunisation has been
characterised in recent work[31, 32]. Here we describe this effect by a linear
titer reduction with time after primary challenge,

$T_{i}^{k}(\Delta t)=T_{i}^{k}-\frac{\Delta t}{\tau},$ (10)

corresponding to an exponential decay of antibody concentration, with a
uniform decay time $90$d (i.e., half life $\tau=65$d). This is broadly
consistent with experimental data; we infer decay times in the range
$[60,170]$d from several studies[31, 32, 8, 9, 44]. In addition, we check that
varying $\tau$ in this range does not affect our results (in particular, the
rank order of variants with respect to antigenic fitness remains unchanged).

### Cross-immunity trajectories.

The cross-immunity factor $c_{i}^{k}$ is defined as the relative reduction in
infections by variant $i$ induced by (recent) immunisation in channel $k$. As
shown in recent work[13, 14], absolute titers of SARS-CoV-2 neutralisation
assays can predict cross-immunity, $c_{i}^{k}=H(T_{i}^{k})$ with

$H(T)=\frac{1}{1+\exp[-\lambda(T-T_{50})]}.$ (11)

This relation has been established in ref. [13] with constants $T_{50}=4.2$
and $\lambda=0.9$. The resulting cross-immunity factors $c_{i}^{k}(\Delta t)$
include antibody decay, as given by equation (10). Hence, they depend on the
time since primary immunisation,

$c_{i}^{k}(\Delta t)=H\big{(}\mbox{$T_{i}^{k}-\frac{\Delta t}{\tau}$}\big{)}.$
(12)

These factors enter the population cross-immunity functions $C_{i}^{k}(t)$,
equation (2), which become

$C_{i}^{k}(t)=\int^{t}H\big{(}T_{i}^{k}-\mbox{$\frac{t-t^{\prime}}{\tau}$}\big{)}\,\dot{y}(t)\,dt^{\prime}.$
(13)

These functions enter all evaluations of the fitness model, equation (3) (Fig.
1c, Fig. 3, Fig. 4ef, Figs. S1, S2, and S4).

To display the emergence of selection hotspots, we approximate the cross-
immunity functions, equation (13), by time-independent effective factors and
immune weights. (i) The effective cross-immunity factors $\bar{c}_{i}^{k}$ are
obtained from equation (12) at a fixed time delay $\Delta t=\tau$ after
primary immunisation, $\bar{c}_{i}^{k}=H(T_{i}^{k}-1)$, which accounts for the
decay of immune response in an approximate way. These factors define the
antigenic landscapes shown in Fig. 4a-c. (ii) The immune weight functions,

$Q_{k}(t)=\int_{-\infty}^{t}H\big{(}T_{0}-\mbox{$\frac{t-t^{\prime}}{\tau}$}\big{)}\,\dot{y}(t^{\prime})\,dt^{\prime},$
(14)

measure the effective population fractions of immune individuals in channel
$k$. They account for immune decay from a fixed reference titer $T_{0}=6.5$
and follow the time-dependence of cross-immunity functions $C_{i}^{k}(t)$ and
selection coefficients $s_{k}(t)$ in an approximate way (Fig. 4de, Fig. S4).
Selection hotspots emerge if large steps on an antigenic landscape,
$\bar{c}_{\rm inv}^{k}-\bar{c}_{\rm anc}^{k}$, coincide with sizeable immune
weights $Q_{k}(t)$ in one or more immune channels (Fig. 4).

### Fitness model.

Equation (3) expresses the fitness of viral variants as a sum of intrinsic and
antigenic fitness components, $f_{i}(t)=f^{0}_{i}+\sum_{k}f_{i}^{k}(t)$.
Intrinsic fitness, $f_{0}$, integrates contributions from several molecular
phenotypes, including protein stability, host receptor binding, and traits
related to intra-cellular viral replication. The antigenic components
$f_{i}^{k}(t)$ describe the impact of antibody binding on viral growth, summed
over the immune repertoire components of different channels of primary
infection or vaccination. Importantly, the input of this fitness model can be
learned by integration of sequence data, epidemiological records, and
antigenic assays.

The additive form the fitness model neglects epistasis between fitness
components. The additivity assumption is justified between intrinsic and
antigenic fitness, because these components are associated to different stages
of the viral replication cycle. The additivity of antigenic fitness components
rests on the approximation of a well-mixed host population and short infection
times. In this approximation, each viral lineage is subject to a dense
sequence of random encounters with hosts of different immune channels $k$,
leading to averaging of antigenic fitness effects. Multiple infections in an
individual can generate additional immune channels; however, these effects are
relatively small over the short periods of SARS-CoV-2 evolution studied in
this paper.

Our analysis of the fitness model focuses on selection coefficients between
co-existing variants,

$s_{ij}(t)\equiv
f_{i}(t)-f_{j}(t)=s_{ij}^{0}-\sum_{k}\gamma_{k}\big{[}C_{i}^{k}(t)-C_{j}^{k}(t)\big{]},$
(15)

because these can directly be compared with their empirical counterparts
$\hat{s}_{ij}(t)$. Of equal importance, selection coefficients within a region
decouple from the changes in viral ecology within that region. Specifically,
seasonality and contact limitations can generate strongly time-dependent
reproductive numbers. However, any modulation of the form
$R(t)\to\alpha(t)R(t)$ leaves the selection coefficients $s_{ij}$ invariant,
as can be seen from equation (5). Our inference of empirical selection, as
described above, is also independent of the underlying infectious period
$\tau_{0}$, which may itself be under evolutionary pressure and change with
time[80, 81]. To keep this independence, we report all selection coefficients
in fixed units [1/d].

### Inference of fitness model parameters.

The free parameters $\gamma_{k}$ ($k=1,\dots,n$) measure the fitness effect of
each cross-immunity component. These parameters calibrate the model to data of
complex real populations differing, for example, in population structure
(including incidence structure), infection histories, and monitoring of
infections. To avoid overfitting, we use a minimal model with just 3 global
antigenic parameters: (i) A basic rate $\gamma_{\rm vac}=\gamma_{\rm bst}$
translates cross-immunity generated by vaccination into units of selection.
(ii) This rate is downweighted to a value $\gamma_{\rm
vac}^{\prime}=a\gamma_{\rm vac}$ for the shift $\delta-o$ and later shifts.
This can be seen as a heuristic to account for the effect of double
infections[42], which increase cross-immunity and decrease cross-immunity
differences between variants. (iii) Cross-immunity in all infection channels
is uniformly upweighted, $\gamma_{k}=b\gamma_{\rm vac}$, to account for
underreporting of infections relative to vaccinations.

Our inference proceeds in two steps. First, we train the antigenic fitness
model using data from the clade shifts $\alpha-\delta$ and $\delta-o$. These
shifts are suitable because they carry sizeable antigenic advance $\Delta
T_{\alpha\delta}^{\rm vac}$ and $\Delta T_{\delta o}^{\rm vac}$ (Fig. 2a),
selection shows a substantial and statistically significant time dependence
(Fig. 1b), and population immunity has started to pick up (Fig. 4d). We infer
the ML likelihood model by aggregation of log likelihood scores over the sets
of regional selection trajectories for the clade shifts $\alpha-\delta$ and
$\delta-o$. We use the score function

$L(\hat{\bf s},{\bf s})=-\sum_{i=1}^{n}\frac{(\Delta\hat{s}(t_{i})-\Delta
s(t_{i}))^{2}}{2\sigma^{2}(t_{i})}$ (16)

for a single empirical selection trajectory $\hat{\bf s}$, equation (6), and
its model-based counterpart ${\bf s}$. This score evaluates selection change,
$\Delta s(t)=s(t)-\langle s\rangle$, where brackets denote averaging over
time. Hence, the fitness model is trained on the time-dependence of selection
in each region, in order to avoid the confounding factor of heterogeneity
across regions. The expected square deviation is
$\sigma^{2}(t_{i})=\sigma^{2}_{s}(t_{i})+\sigma_{0}^{2}$; the first term
describes the sampling error of sequence counts, which enters frequency and
empirical selection estimates, the second term summarises fluctuations
unrelated to sequence counts. The total log likelihood score is the sum
$L=\sum L({\bf s},\hat{\bf s})$, which runs over both shifts and all included
regions. Table S2 lists the ML parameters $\gamma_{\rm vac},a,b$ and the ML
score $L$ relative to a null model of time-independent selection (see below).
The 95% confidence intervals of the inferred parameters are computed by
resampling the empirical selection data with fluctuations $\sigma^{2}$. We
note that the ML values $a<1,b>1$ are consistent with the interpretation of
these parameters as weighting factors accounting for double-infections and
underreporting (see above). Second, we infer the intrinsic selection for each
shift as the difference between empirical selection and ML antigenic
selection, $s_{0}=\langle\\!\langle\hat{s}-s_{\rm ag}\rangle\\!\rangle$, where
the double brackets denote averaging over time and regions. The ML antigenic
selection coefficients, $\langle\\!\langle s_{k}\rangle\\!\rangle$, and the
intrinsic selection coefficient $s_{0}$ between invading and ancestral variant
are listed in Table S3; see also Fig. 3b. Confidence intervals are computed by
resampling model parameters with their confidence intervals. Consistently, we
infer weak antigenic selection for the shifts ${\sl 1}-\alpha$ and $o1-o2$,
which also show only weak time dependence of selection (Fig. S3).

### Significance analysis of the fitness model.

To assess the statistical significance of our inference, we compare four
fitness models of the form (3): the full model used in the main text (VI:
antigenic selection by vaccination and infection, intrinsic selection), two
partial models (V: antigenic selection only by vaccination, intrinsic
selection; I: antigenic selection only by infection, intrinsic selection), and
a null model (0: intrinsic selection only). We infer conditional ML parameters
for each model and we rank models by their ML score difference to the null
model, $\Delta L=L-L_{0}$ (Table S2). An alternative ranking by BIC score[82],
which contains a score penalty for the number of model parameter, leads to the
same result. We observe the following: (i) All antigenic fitness models have
significantly higher scores than the null model, which shows that the
empirical selection data are incompatible with time-independent selection.
(ii) The full model has a significantly higher score than any of the other
models; both vaccination and infection are significant components of antigenic
selection. (iii) Vaccination explains a larger part of the time-dependent data
than infection ($\Delta L_{V}>\Delta L_{I}$), which is consistent with the
ranking of selection coefficients inferred from the full model (Fig. 3b, Table
S3). (iv) The score gain of the full model is less than the sum of its parts
($\Delta L_{IV}<\Delta L_{V}+\Delta L_{I}$). This can be associated with
statistical correlations in the input data for both antigenic model
components. For example, the fraction of vaccinated individuals $y_{\rm vac}$
is weakly anti-correlated with the fraction of $\delta$ infections,
$y_{\delta}$.

### Fitness trajectories.

Long-term fitness trajectories display clade turnover of multiple successive
shifts (Fig. 4f, Fig. S4). For each of the variants
$\alpha,\delta,o(o1),o2,o45$, we plot the time-dependent fitness gap, $\delta
f_{i}(t)=f_{i}(t)-\bar{f}(t)$, where $\bar{f}(t)=\sum_{j}x_{j}(t)f_{j}(t)$ is
the mean population fitness. Like selection coefficients, fitness gaps
decouple from ecological factors affecting absolute growth (see the discussion
above). Assuming that the fitness difference between ancestral and invading
variant, $s(t)$, is dominant during each crossover, we obtain the fitness gap
trajectories $\delta f_{\rm anc}(t)=-s(t)x_{\rm inv}(t)$ and $\delta f_{\rm
inv}(t)=-s(t)[1-x_{\rm inv}(t)]$, as well as their empirical counterparts
$\delta\hat{f}_{\rm anc}(t)$ and $\delta\hat{f}_{\rm inv}(t)$. For each
variant, we patch trajectories from origination to near-fixation (here in the
time interval $(t_{0,i},t_{f,i})$ given by $x_{i}(t_{0,i})=0.01$ and
$x_{i}(t_{f,i})=0.99$). The long-term trajectories display selection hotspots
and confirm the quantitative agreement between empirical and model-based
fitness.

### Model-based inference of selection across regions.

As shown by the preceding analysis, we can infer a statistically significant
fitness model with few, global parameters from sequence and epidemiological
data aggregated over a set of regions and combined with antigenic data. The
model describes common time-dependent patterns of selection in these regions
and serves two main purposes: to provide a breakdown of selection in intrinsic
and antigenic components (Fig. 3) and to display selection hotspots in long-
term trajectories (Fig. 4). Our inference procedure rests on stringent
criteria for the joint availability of sequence and epidemiological data in
each of these regions (as listed above). A number of points support this
procedure: (i) The results are robust under variation of the inclusion
criteria for regions. In particular, the signal of antigenic selection in data
and model is broadly distributed over regions (Figs. S1-S3). Hence, the
selection averages reported in Fig. 3b and Table S3 are reproducible in
subsampled sets of regions. (ii) Within the set of regions included, the model
is applicable beyond the $\alpha-\delta$ and $\delta-o$ shifts used for
training. The early ${\sl 1}-\alpha$ shift and the recent $o1-o2$ shifts serve
as controls. In both cases, we infer weak antigenic selection, consistent with
weak time dependence of empirical selection (Fig. S4). For the emerging
$o2-o45$ shift, strong antigenic selection is consistent with fast initial
growth of the new variants.

Our model-based inference of selection excludes a number of regions that do
not fulfil the criteria of joint data availability. (i) For the
$\alpha-\delta$ shift, several regions are excluded because VOCs other than
$\alpha$ were majority variants prior to the shift to $\delta$ (for example,
Beta in South Africa and Gamma in Brazil). Unlike $\alpha-\delta$, these
shifts do not involve antigenic advance in the vaccination channel; i.e.,
$\Delta T^{\rm vac}_{{\rm anc}\,\delta}<0$. However, we lack comprehensive
antigenic data on other VOCs as input for the fitness model. (ii) For the
$\delta-o$ shift, several regions are excluded because of low reported
incidence (e.g., Brazil, California, India, Mexico, Poland, South Korea,
Turkey, Texas). Most of these regions show a signal of time-dependent
selection consistent with the regions included; however, much lower reported
incidence counts prevent reliable immune tracking of infection channels during
the $\delta-o$ shift. Variation in reported incidence can, in principle, be
incorporated into the fitness model by region-dependent $\gamma_{\rm vac}$
factors, but this would likely lead to overfitting. We conclude that at
current levels of data availability, a comprehensive cross-regional analysis
is not feasible.

### Model-free inference of selection in regions with low vaccination
coverage.

Countries with low vaccination coverage during the $\alpha-\delta$ shift
($y_{\rm vac}<0.1$) disqualify for the model-based analysis because $\alpha$
was not a majority variant (India, Malaysia, Russia, Philippines, Indonesia,
South Africa, South Korea) or sequence counts are too low for the inference of
selection trajectories (Australia). For this set of countries, we can still
infer a selection coefficient $s$ by fitting a sigmoid function to the
frequency trajectory $x_{\delta}(t)$ (to be interpreted as the growth
difference between $\delta$ and the average of all other coexisting variants).
We find lower region-averaged selection in the set of low-vaccination
countries compared to other countries
($\langle\\!\langle\hat{s}\rangle\\!\rangle=0.08$ vs.
$\langle\\!\langle\hat{s}\rangle\\!\rangle=0.12$). This is qualitatively
consistent with our model-based inference of vaccination-induced selection;
however, the genetic heterogeneity of this clade shift prevents a systematic
breakdown of antigenic selection into immune channels.

### Data availability.

The datasets analysed in this study are available in published work.

### Code availability.

The code used in this study is available at
https://github.com/m-meijers/vaccine_effect

### Acknowledgements.

We thank Florian Klein and Kanika Vanshylla for discussions. This work has
been partially funded by Deutsche Forschungsgemeinschaft grant CRC 1310
Predictability in Evolution.

##

  * 1\.  Ozono, S. _et al._ SARS-CoV-2 D614G spike mutation increases entry efficiency with enhanced ACE2-binding affinity.  _Nature Communications_ 12 848 (2021).
  * 2\.  Meng, B. _et al._ Recurrent emergence of SARS-CoV-2 spike deletion H69/V70 and its role in the Alpha variant B.1.1.7.  _Cell Reports_ 35 13 (2021).
  * 3\.  Mlcochova, P. _et al._ SARS-CoV-2 B.1.617.2 Delta variant replication and immune evasion.  _Nature_ 599 114-119 (2021).
  * 4\.  Harvey, W. T. _et al._ SARS-CoV-2 variants, spike mutations and immune escape  _Nature reviews microbiology_ 19 409-424 (2021).
  * 5\.  Shu, Y. & McCauley, J.  GISAID: Global initiative on sharing all influenza data – from vision to reality  _EuroSurveillance_ 22 13 (2017).
  * 6\.  Ritchie, H. _et al._ Coronavirus pandemic (covid-19).  _Our World in Data_ (2020).  Https://ourworldindata.org/coronavirus.
  * 7\.  Centers for Disease Control and Prevention. COVID Data Tracker. Atlanta, GA: US Department of Health and Human Services, CDC.  https://covid.cdc.gov/covid-data-tracker.  Accessed: 22-06-2022.
  * 8\.  Planas, D. _et al._ Reduced sensitivity of SARS-CoV-2 variant Delta to antibody neutralization.  _Nature_ 596, 276–280 (2021).
  * 9\.  Planas, D. _et al._ Considerable escape of SARS-CoV-2 Omicron to antibody neutralization.  _Nature_ 602 671–675 (2021).
  * 10\.  Garcia-Beltran, W. F. _et al._ mRNA-based COVID-19 vaccine boosters induce neutralizing immunity against SARS-CoV-2 Omicron variant.  _Cell_ 185, 457–466 (2022).
  * 11\.  Wrobel, A. G. _et al._ SARS-CoV-2 and bat RaTG13 spike glycoprotein structures inform on virus evolution and furin-cleavage effects.  _Nature Structural and Molecular Biology_ 27, 763–767 (2020).
  * 12\.  Zeng, C. _et al._ Neutralization and stability of SARS-CoV-2 Omicron variant.  _bioRxiv:_ (2021).  URL https://doi.org/10.1101/2021.12.16.472934.
  * 13\.  Khoury, D. S. _et al._ Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection.  _Nature Medicine_ 27, 1205–1211 (2021).
  * 14\.  Feng, S. _et al._ Correlates of protection against symptomatic and asymptomatic SARS-CoV-2 infection.  _Nature Medicine_ 27, 2032–2040 (2021).
  * 15\.  Tracking SARS-CoV-2 variants.  https://www.who.int/en/activities/tracking-SARS-CoV-2-variants.  Accessed: 22-06-2022.
  * 16\.  Davies, N. G. _et al._ Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England.  _Science_ 372 (2021).
  * 17\.  Dhar, M. S. _et al._ Genomic characterization and epidemiology of an emerging SARS-CoV-2 variant in Delhi, India.  _Science_ 374, 995–999 (2021).
  * 18\.  Kepler, L., Hamins-Puertolas, M. & Rasmussen, D. A.  Decomposing the sources of SARS-CoV-2 fitness variation in the united states.  _Virus Evolution_ 7 (2021).
  * 19\.  Ulrich, L. _et al._ Enhanced fitness of SARS-CoV-2 variant of concern Alpha but not Beta.  _Nature_ 602, 307–313 (2022).
  * 20\.  Grenfell, B. T. _et al._ Unifying the epidemiological and evolutionary dynamics of pathogens.  _Science_ 303, 327–332 (2004).
  * 21\.  Rella, S. A., Kulikova, Y. A., Dermitzakis, E. T. & Kondrashov, F. A.  Rates of SARS-CoV-2 transmission and vaccination impact the fate of vaccine-resistant strains.  _Scientific Reports_ 11 (2021).
  * 22\.  Saad-Roy, C. M. _et al._ Epidemiological and evolutionary considerations of SARS-CoV-2 vaccine dosing regimes.  _Science_ 372, 363–370 (2021).
  * 23\.  Lobinska, G. _et al._ Evolution of resistance to COVID-19 vaccination with dynamical social distancing.  _Nature_ 507, 57–61 (2022).
  * 24\.  Łuksza, M. & Lässig, M.  A predictive fitness model for influenza.  _Nature Human Behaviour_ 6, 57–61 (2014).
  * 25\.  Wen, F. _et al._ The potential beneficial effects of vaccination on antigenically evolving pathogens.  _American Naturalist_ 199 2 193-206 (2022).
  * 26\.  Wen, F. T., Bell, S. M., Bedford, T. & Cobey, S.  Estimating vaccine-driven selection in seasonal influenza.  _Viruses_ 10 (2018).
  * 27\.  Mullen, J. L. _et al._ Outbreak.info (2020).  URL https://outbreak.info.
  * 28\.  van der Straten, K. _et al._ Mapping the antigenic diversification of SARS-CoV-2.  _bioRxiv_ (2022).  URL https://doi.org/10.1101/2022.01.03.21268582.
  * 29\.  Wilks, S. H. _et al._ Mapping SARS-CoV-2 antigenic relationships and serological responses.  _bioRxiv_ (2022).  URL https://doi.org/10.1101/2022.01.28.477987.
  * 30\.  Polack, F. P. _et al._ Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine.  _New England Journal of Medicine_ 383, 2603–2615 (2020).
  * 31\.  Iyer, A. S. _et al._ Persistence and decay of human antibody responses to the receptor binding domain of SARS-CoV-2 spike protein in COVID-19 patients.  _Science Immunology_ 5 (2020).
  * 32\.  Israel, A. _et al._ Large‐scale study of antibody titer decay following BNT162b2 mRNA vaccine or SARS‐CoV‐2 infection.  _Vaccines_ 10 (2022).
  * 33\.  Coudeville, L. _et al._ Relationship between haemagglutination-inhibiting antibody titres and clinical protection against influenza: development and application of a bayesian random-effects model.  _BMC Medical Research Meth._ 10 (2010).
  * 34\.  Dunning, A. J. _et al._ Correlates of protection against influenza in the elderly: Results from an influenza vaccine efficacy trial.  _Clinical and Vaccine Immunology_ 23, 228–235 (2016).
  * 35\.  Rotem, A. _et al._ Evolution on the biophysical fitness landscape of an rna virus.  _Molecular Biology and Evolution_ 35, 2390–2400 (2018).
  * 36\.  Meijers, M., Vanshylla, K., Gruell, H., Klein, F. & Laessig, M.  Predicting in vivo escape dynamics of HIV-1 from a broadly neutralizing antibody.  _PNAS_ 118 (2021).
  * 37\.  Morris, D. H. _et al._ Predictive modeling of influenza shows the promise of applied evolutionary biology (2018).
  * 38\.  Huddleston, J. _et al._ Integrating genotypes and phenotypes improves long-term forecasts of seasonal influenza A/H3N2 evolution.  _eLife_ 9, 1–48 (2020).
  * 39\.  Gog, J. R. & Grenfell, B. T.  Dynamics and selection of many-strain pathogens.  _PNAS_ 99, 17209–17214 (2002).
  * 40\.  Lipsitch, M. _et al._ Mapping person-to-person variation in viral mutations that escape polyclonal serum targeting influenza hemagglutinin (2019).
  * 41\.  Cobey, S. & Hensley, S. E.  Immune history and influenza virus susceptibility.  _Current Opinion in Virology_ 22, 105–111 (2017).
  * 42\.  Bates, T. A. _et al._ Vaccination before or after SARS-CoV-2 infection leads to robust humoral response and antibodies that effectively neutralize variants _Science Immunology_ 7 68 (2022).
  * 43\.  Yuan, S. _et al._ Pathogenicity, transmissibility, and fitness of SARS-CoV-2 Omicron in Syrian hamsters.  _Science_ (2022).  URL https://www.science.org/doi/10.1126/science.abn8939.
  * 44\.  Gruell, H. _et al._ mRNA booster immunization elicits potent neutralizing serum activity against the SARS-CoV-2 Omicron variant.  _Nature Medicine_ 28, 477–480 (2022).
  * 45\.  Hachmann, N. P. _et al._ Neutralization escape by SARS-CoV-2 Omicron subvariants BA.2.12.1, BA.4, and BA.5.  _New England Journal of Medicine_ (2022).
  * 46\.  Earn, D. J. D., Dushoff, J. & Levin, S. A.  Ecology and evolution of the flu.  _Trends in Ecology and Evolution_ 17, 334–340 (2002).
  * 47\.  Boni, M. F., Gog, J. R., Andreasen, V. & Feldman, M. W.  Epidemic dynamics and antigenic evolution in a single season of influenza a.  _Proceedings of the Royal Society B: Biological Sciences_ 273, 1307–1316 (2006).
  * 48\.  Khan, K. _et al._ Omicron sub-lineages ba.4/ba.5 escape ba.1 infection elicited neutralizing immunity  _bioRxiv_ (2022).  URL https://doi.org/10.1101/2022.04.29.22274477.
  * 49\.  Fonville, J. M. _et al._ Antibody landscapes after influenza virus infection or vaccination.  _Science_ 346, 996–1000 (2014).
  * 50\.  Strelkowa, N. & Lässig, M.  Clonal interference in the evolution of influenza.  _Genetics_ 192, 671–682 (2012).
  * 51\.  Gong, L. I., Suchard, M. A. & Bloom, J. D.  Stability-mediated epistasis constrains the evolution of an influenza protein.  _eLife_ 2013 (2013).
  * 52\.  Lässig, M., Mustonen, V. & Walczak, A. M.  Predicting evolution.  _Nature Ecology and Evolution_ 1 (2017).
  * 53\.  Moulana, A. _et al._ Compensatory epistasis maintains ace2 affinity in sars-cov-2 omicron ba.1  _bioRxiv_ (2022).  URL https://doi.org/10.1101/2022.06.17.496635.
  * 54\.  Mustonen, V. & Lässig, M.  Eco-evolutionary control of pathogens.  _PNAS_ 117, 19694–19704 (2020).
  * 55\.  Benson, D. A. _et al._ GenBank.  _Nucleic Acids Research_ 43, D30–D35 (2015).
  * 56\.  Katoh, K. & Standley, D. M.  MAFFT multiple sequence alignment software version 7: Improvements in performance and usability.  _Molecular Biology and Evolution_ 30, 772–780 (2013).
  * 57\.  Minh, B. Q. _et al._ IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era.  _Molecular Biology and Evolution_ 37, 1530–1534 (2020).
  * 58\.  Hoang, D. T. _et al._ UFBoot2: Improving the ultrafast bootstrap approximation.  _Mol. Biol. Evol_ 35, 518–522 (2017).
  * 59\.  Sagulenko, P., Puller, V. & Neher, R. A.  TreeTime: Maximum-likelihood phylodynamic analysis.  _Virus Evolution_ 4 (2018).
  * 60\.  Kingman, J. F. C.  The coalescent.  _Stochastic Processes and their Applications_ 13, 235–248 (1982).
  * 61\.  Smith, D. J., Forrest, S., Hightower, R. R. & Perelson, A. S.  Deriving shape space parameters from immunological data.  _J. theor. Biol_ 189, 141–150 (1997).
  * 62\.  Smith, D. J., Lapedes, A. S. & Jong, J. C. D.  Mapping the antigenic and genetic.  _Science_ 305, 371–377 (2004).
  * 63\.  Bates, T. A. _et al._ Neutralization of SARS-CoV-2 variants by convalescent and BNT162b2 vaccinated serum.  _Nature Communications_ 12 (2021).
  * 64\.  Zhou, D. _et al._ Evidence of escape of SARS-CoV-2 variant B.1.351 from natural and vaccine-induced sera.  _Cell_ 184, 2348–2361.e6 (2021).
  * 65\.  Garcia-Beltran, W. F. _et al._ Multiple SARS-CoV-2 variants escape neutralization by vaccine-induced humoral immunity.  _Cell_ 184, 2372–2383.e9 (2021).
  * 66\.  Liu, C. _et al._ Reduced neutralization of SARS-CoV-2 B.1.617 by vaccine and convalescent serum.  _Cell_ 184, 4220–4236.e13 (2021).
  * 67\.  Wang, P. _et al._ Antibody resistance of SARS-CoV-2 variants B.1.351 and B.1.1.7.  _Nature_ 593, 130–135 (2021).
  * 68\.  Liu, Y. _et al._ BNT162b2-elicited neutralization against new SARS-CoV-2 spike variants.  _New England Journal of Medicine_ 385, 472–474 (2021).
  * 69\.  Planas, D. _et al._ Sensitivity of infectious SARS-CoV-2 B.1.1.7 and B.1.351 variants to neutralizing antibodies.  _Nature Medicine_ 27, 917–924 (2021).
  * 70\.  Uriu, K. _et al._ Neutralization of the sars-cov-2 mu variant by convalescent and vaccine serum.  _New England Journal of Medicine_ 385, 2395–2397 (2021).
  * 71\.  Rössler, A., Riepler, L., Bante, D., von Laer, D. & Kimpel, J.  SARS-CoV-2 Omicron variant neutralization in serum from vaccinated and convalescent persons.  _New England Journal of Medicine_ 386, 698–700 (2022).
  * 72\.  Cameroni, E. _et al._ Broadly neutralizing antibodies overcome sars-cov-2 omicron antigenic shift.  _Nature_ 602, 664-670 (2022).
  * 73\.  Wang, Q. _et al._ Antibody evasion by sars-cov-2 omicron subvariants ba.2.12.1, ba.4; ba.5.  _Nature_ (2022).
  * 74\.  Iketani, S. _et al._ Antibody evasion properties of sars-cov-2 omicron sublineages.  _Nature_ 604, 553–556 (2022).
  * 75\.  Muik, A. _et al._ Neutralization of SARS-CoV-2 lineage B.1.1.7 pseudovirus by BNT162b2 vaccine-elicited human sera  _bioRxiv_ (2021).  URL https://doi.org/10.1101/2020.12.30.20249034.
  * 76\.  Bowen, J. E. _et al._ Omicron ba.1 and ba.2 neutralizing activity elicited by a comprehensive panel of human vaccines.  _bioRxiv_ (2022).  URL http://www.ncbi.nlm.nih.gov/pubmed/35313570.
  * 77\.  Cao, Y. _et al._ Ba.2.12.1, ba.4 and ba.5 escape antibodies elicited by omicron infection division of hiv/aids and sex-transmitted virus vaccines, institute for biological product  _bioRxiv_ (2022).  URL https://doi.org/10.1101/2022.04.30.489997.
  * 78\.  Mykytyn, A. Z. _et al._ Omicron ba.1 and ba.2 are antigenically distinct sars-cov-2 variants  _bioRxiv_ (2022).  URL https://doi.org/10.1101/2022.02.23.481644.
  * 79\.  Neher, R. A., Bedford, T., Daniels, R. S., Russell, C. A. & Shraiman, B. I.  Prediction, dynamics, and visualization of antigenic phenotypes of seasonal influenza viruses.  _Proceedings of the National Academy of Sciences_ 113, E1701–E1709 (2016).
  * 80\.  Hart, W. S. _et al._ Generation time of the alpha and delta SARS-CoV-2 variants: an epidemiological analysis.  _The Lancet Infectious Diseases_ 22, 603–610 (2022).
  * 81\.  Hart, W. S. _et al._ Inference of the SARS-CoV-2 generation time using UK household data.  _eLife_ 11 (2022).
  * 82\.  Schwarz, G.  Estimating the dimension of a model.  _The Annals of Statistics_ 6, 461–464 (1978).

## Supplementary Tables and Figures

Table S1: Antigenic data.

| $\alpha$ | $\delta$ |  $o$ ($o1$) | $o2$ | $o45$ | vac | bst
---|---|---|---|---|---|---|---
$\alpha$ | 0 | 1.8 | _5.0_ | _5.0_ | _5.0_ | 0.8 | _0.8_
$\delta$ | 1.5 | 0 | _4.8_ | _4.8_ | _4.8_ | 1.7 | 1.5
$o$ ($o1$) | 5.0 | 4.8 | 0 | 2.1 | _2.2_ | 5.6 | 2.7
$o2$ | $<\\!\emph{5.0}$ |  $<\\!$ _4.8_ | 1.4 | 0 | _1.2_ |  $<\\!$ _5.6_ | 2.5
$o45$ |  $<\\!$ _5.0_ |  $<\\!$ _4.8_ | 2.2 | 1.2 | 0 |  $<\\!$ _5.6_ | 3.9

We list log titer drops, or neutralisation fold changes, $\Delta
T_{i}^{k}=T_{*}^{k}-T_{i}^{k}$, of strains from variant $i$ assayed against
human antisera induced by primary immunisation (infection or vaccination) with
strains of variant $k$ (columns). Numbers are average values of primary data
from ref. [63, 72, 65, 66, 64, 67, 68, 3, 69, 8, 75, 70, 71, 44, 9, 28, 29,
74, 76, 77, 78, 45, 73, 48]. All vaccination titers refer to mRNA vaccines.
Where no primary data is available, titer drops are inferred by symmetry or
(as lower bounds) by genetic similarity (numbers in italics, Methods).
Absolute titers $T_{i}^{k}$ are shifted by the reference titers
$T_{*}^{k}=6.5$, ($k=\alpha,\delta,o(o1),o2,o45$), $T_{*}^{\rm vac}=7.8$,
$T_{*}^{\rm bst}=9.8$ obtained from ref. [30, 9, 44, 10]; see Methods and Fig.
2a.

Table S2: Ranking of fitness models.

model | antigenic parameters | posterior scores
---|---|---
| $\hskip 14.22636pt\gamma_{\rm vac}\hskip 14.22636pt$ |  $a$ |  $b$ |  $\Delta L\hskip 14.22636pt$ |  $\Delta H$
VI | $1.22\pm 0.03$ | $0.24\pm 0.03$ | $2.0\pm 0.5$ | 947 | -1883
V | 1.22 | 0.34 | - | 882 | -1758
I | 4.9 | 0.21 | - | 378 | -750
0 | - | - | - | 0 | 0

We compare the full fitness model used in the main text (VI: vaccination +
infection + intrinsic selection) with partial models (V: vaccination +
intrinsic selection, I: infection + intrinsic selection) and a null model (0:
intrinsic selection only). Columns from left to right: model parameters,
$\gamma_{\rm vac}$, $a$, $b$, ML values and 95% confidence intervals
(definitions are given in Methods); log likelihood difference to the null
model, $\Delta L$; BIC score difference to the null model, $\Delta H$.

Table S3: Intrinsic and antigenic selection components.

clade shift | selection coefficients
---|---
| ${\alpha}$ | ${\delta}$ | ${o\,(o1)}$ | ${o2}$ | ${\rm vac}$ | ${\rm bst}$ | 0 | $s$
${\sl 1}-\alpha$ | $<$ | - | - | - | $<$ | - | $.08\pm.001$ | $.08\pm.01$
$\alpha-\delta$ | $.01\pm.001$ | $<$ | - | - | $.04\pm.002$ | - | $.05\pm.002$ | $.09\pm.02$
$\delta-o$ | $<$ | $.02\pm.006$ | $-.01\pm.002$ | - | $.06\pm.01$ | $-.01\pm.002$ | $.06\pm.01$ | $.14\pm.03$
$o1-o2$ | $<$ | $<$ | $.01\pm.004$ | $<$ | $<$ | $<$ | $.08\pm.002$ | $.08\pm.01$
$o2-o45$ | $<$ | $<$ | $.01\pm.002$ | $.01\pm.002$ | $<$ | $.04\pm.01$ | $.06\pm.008$ | $.12\pm.02$

Selection coefficients between the invading and the ancestral variant,
$s=f_{\rm inv}-f_{\rm anc}$, and their decomposition into antigenic and
intrinsic components are inferred for the full fitness model; all values are
time averages for each clade shift. Rows from top to bottom: major clade
shifts, ${\sl 1}-\alpha$, $\alpha-\delta$, $\delta-o$; recent clade shifts,
$o1-o2$, $o2-o45$ (shift incomplete, entries refer to initial period). Columns
from left to right: average antigenic selection in immune channels $k=\alpha$,
$\delta$, $o$ ($o1$), $o2$, ${\rm vac}$, ${\rm bst}$; intrinsic selection (0);
total selection ($s$). Selection coefficients are given in units [1/day]; the
symbol “$<$” marks values $s<0.01$. We list ML values with 95% confidence
intervals (for selection components) or with rms cross-region variation of
selection (for $s$; cf. Fig. 1b).

(continued on next page)

Fig. S1: Empirical and model-based trajectories of the $\alpha-\delta$ shift.
Evolutionary, epidemiological, and cross-immune trajectories are shown for all
regions of this study. (a) Observed frequency trajectories of relevant clades,
$x_{i}(t)$; rms sampling error is indicated by shading. (b) Cumulative
coverage of primary vaccination, $y_{\rm vac}(t)$ (light gray), and of booster
vaccination, $y_{\rm bst}(t)$ (dark gray); cumulative population fraction of
$\alpha$ infections, $y_{\alpha}(t)$ (purple), and of $\delta$ infections,
$y_{\delta}(t)$ (green). (c) Population immunity functions, $C_{i}^{k}(t)$ (as
in Fig. 1c). (d) Empirical selection change, $\Delta\hat{s}(t)$ (dots, with
rms statistical errors indicated by bars), together with ML model prediction,
$\Delta s(t)$ (dashed line). Criteria for inclusion of regions are given in
Methods.

Fig. S1: (continued)

(continued on next page)

Fig. S2: Empirical and model-based trajectories of the $\delta-o$ shift.
Evolutionary, epidemiological, and cross-immune trajectories are shown for all
regions of this study. (a) Observed frequency trajectories of relevant clades,
$x_{i}(t)$; rms sampling error is indicated by shading. (b) Cumulative
coverage of primary vaccination, $y_{\rm vac}(t)$ (light gray), and of booster
vaccination, $y_{\rm bst}(t)$ (dark gray); cumulative population fraction of
$\delta$ infections, $y_{\delta}(t)$ (green), and of $o$ infections,
$y_{o}(t)$ (orange). (c) Population immunity functions, $C_{i}^{k}(t)$ (as in
Fig. 1c). (d) Empirical selection change, $\Delta\hat{s}(t)$ (dots, with rms
statistical errors indicated by bars), together with ML model prediction,
$\Delta s(t)$ (dashed line). Criteria for inclusion of regions are given in
Methods.

Fig. S2: (continued)

Fig. S3: Selection tracking in multiple regions and clade shifts. Empirical
selection change between invading and ancestral clade,
$\Delta\hat{s}(t)=\hat{s}(t)-\langle s\rangle$, for all complete clade shifts
and all regions of this study (brackets denote time averages for each
trajectory). Selection trajectories are derived from regional frequency
trajectories and plotted against time counted from the midpoint (colored
lines); rms statistical error is indicated by shading. Summary statistics:
cross-region linear regression, $s_{\rm lin}(t)$ (black solid line, length
gives r.m.s. time span of trajectories). (a) ${\sl 1}-\alpha$ shift: small,
statistically insignificant time dependence, ${\rm Var}(s_{\rm lin})=3.6\times
10^{-4}$, $P>0.01$; (b) $\alpha-\delta$ shift: substantial, statistically
significant time dependence, ${\rm Var}(s_{\rm lin})=2.\times 10^{-3}$,
$P<10^{-16}$; (c) $\delta-o$ shift: substantial, statistically significant
time dependence, ${\rm Var}(s_{\rm lin})=1.4\times 10^{-3}$, $P<10^{-5}$; (a)
$o1-o2$ shift: small, but statistically significant time dependence, ${\rm
Var}(s_{\rm lin})=2.9\times 10^{-4}$, $P<10^{-4}$. All $P$ values are computed
using a two-sided Wald test. The statistical grading of shifts is described
and criteria for inclusion of regions are given in Methods.

Fig. S4: Regional long-term trajectories of immune weight and fitness. (a)
Time-dependent weight factors of different immune classes, $Q_{k}(t)$. (b)
Time-dependent fitness gap, $\delta f_{i}(t)$. Criteria for inclusion of
regions are given in Methods; see Fig. 4 for averaged trajectories.