Antigenic escape selects for the evolution of higher pathogen transmission and virulence

Oligomorphic dynamics (OMD) of antigenic escape

We considered a model of the antigenic escape of a pathogen from host herd immunity on a one-dimensional antigenicity space (x). We tracked the changes in the density (S(t,x)) of hosts that are susceptible to antigenicity variant x of pathogen at time t, and the density (I(t,x)) of hosts that are currently infected and infectious with antigenicity variant x of pathogen at time t:

where β, α and γ are the transmission rate, virulence (additional mortality due to infection) and recovery rate of pathogens, which are independent of antigenicity. The function σ(x−y) denotes the degree of cross immunity: a host infected by pathogen variant y acquires perfect cross immunity with probability σ(x−y), but fails to acquire any cross immunity with probability 1−σ(x−y) (this is called polarized cross immunity by Gog and Grenfell²⁵). The degree σ(x−y) of cross immunity is assumed to be a decreasing function of the distance |x−y| between variants x and y. When a new variant with antigenicity x = 0 is introduced at time t = 0, the initial host population is assumed to be susceptible to any antigenicity variant of pathogen: S(0,x) = 1. In equation (6), (D = mu sigma _mathrm{m}^2/2) is one half of the mutation variance for the change in antigenicity, representing random mutation in the continuous antigenic space.

Susceptibility profile moulded by the primary outbreak

We first analysed the dynamics of the primary outbreak of a pathogen and derived the resulting susceptibility profile, which can be viewed as the fitness landscape subsequently experienced by the pathogen. For simplicity, we assumed that mutation can be ignored during the first epidemic initiated with antigenicity strain x = 0. The density ({{{S}}}_0left( {{{t}}} right) = {{{S}}}({{{t}}},0)) of hosts that are susceptible to the currently prevailing antigenicity variant x = 0, as well as the density ({{{I}}}_0left( {{{t}}} right) = {{{I}}}({{{t}}},0)) of hosts that are currently infected by the focal variant change with time as

with (S_0left( 0 right) = 1), (I_0left( 0 right) approx 0) and (R_0left( 0 right) = 0). The final size of the primary outbreak,

is determined as the unique positive root of

where (rho _0 = beta /left( {gamma + alpha } right) > 1) is the basic reproductive number⁶. Associated with this epidemiological change, the susceptibility profile (S_xleft( t right) = S(t,x)) against antigenicity x ((x ne 0)) other than the currently circulating variant (x = 0) changes by cross immunity as

Integrating both sides of equation (11) from t = 0 to (t = infty), we see that the susceptibility profile (sleft( x right) = S_x(infty )) after the primary outbreak at x = 0 is

where the last equality follows from equation (10). The susceptibility can be effectively reduced by cross immunity when the primary variant has a large impact (that is, when the fraction of hosts remaining uninfected, 1−ψ₀, is small) and when the degree of cross immunity is strong (that is, when σ(x) is close to 1). With a variant antigenically very close to the primary variant (x ≈ 0), the cross immunity is very strong ((sigma left( x right) approx 1)) so that the susceptibility against variant x is nearly maximally reduced: (s(x) approx 1 – psi _0). With a variant antigenically distant from the primary variant, σ(x) becomes substantially smaller than 1, making the host more susceptible to the variant. For example, if the cross immunity is halved ((sigma left( x right) = 0.5)) from its maximum value 1, then the susceptibility to that variant is as large as (left( {1 – psi _0} right)^{0.5}). If a variant is antigenically very distant from the primary variant, then (sigma left( x right) approx 0), and the host is nearly fully susceptible to the variant ((sleft( x right) approx 1)).

Threshold antigenic distance for escaping immunity raised by primary outbreak

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

or

where we used equation (12). With a specific choice of cross-immunity profile,

the threshold antigenicity beyond which the virus can increase in the susceptibility profile s(x) after the primary outbreak is obtained, by substituting equation (14) into equation (13)

and taking the logarithm of both sides twice:

OMD

Integrating both sides of equation (6) over the whole space, we obtained the dynamics for the total density of infected hosts, ({{{bar{ I}}}}left( {{{t}}} right) = {int}_{ – infty }^infty {{{{I}}}left( {{{{t}}},{{{x}}}} right){{{dx}}}}):

where

is the relative frequency of antigenicity variant x in the pathogen population circulating at time t, and

is the mean susceptibility experienced by currently circulating pathogens. The dynamics for the relative frequency (phi left( {t,x} right)) of pathogen antigenicity is

As in Sasaki and Dieckmann²⁷, we decomposed the frequency distribution to the sum of several morph distributions (oligomorphic decomposition) as

where p_i(t) is the frequency of morph i and (phi _i(t,x)) is the within-morph distribution of antigenicity. By definition, and ({int}_{ – infty }^infty {phi _i} (t,x)dx = 1). Let

be the mean antigenicity of a morph and

where O is order be the within-morph variance of each morph, which is assumed to be small, of the order of ({it{epsilon }}^2). We denoted the mean susceptibility of host population for viral morph (i) by (bar S_i = {int}_{ – infty }^infty {Sleft( {t,x} right)phi _i(t,x)dx}). As shown in Sasaki and Dieckmann²⁷, the dynamics for viral morph frequency is expressed as

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

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

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

are derived, where (xi _i = x – bar x_i) and (Eleft[ {xi _i^4} right] = {int}_{ – infty }^infty {left( {x – bar x_i} right)^4phi _ileft( {t,x} right)dx}) are the fourth central moments of antigenicity around the morph mean. Assuming that the within-morph distribution is normal (Gaussian closure), (Eleft[ {xi _i^4} right] = 3V_i^2), and hence equation (25) becomes

Second outbreak predicted by OMD

Equations (22), (24) and (26) are general, but they rely on a full knowledge of the dynamics of the susceptibility profile S(t,x). To make further progress, we used an additional approximation by substituting equation (13), the susceptibility profile, over viral antigenicity space after the primary outbreak at x = 0 and before the onset of the second outbreak at a distant position. We kept track of two morphs at positions x₀(t) and x₁(t), where the first morph is that caused by the primary outbreak at x = 0, and the second morph is that emerged in the range x > x_c beyond the threshold antigenicity x_c defined in equation (13) (and equation (15) for a specific form of σ(x)) as the source of the next outbreak.

As (sleft( x right) = left( {1 – psi _0} right)^{sigma (x)} = exp [sigma left( x right)log (1 – psi _0)]), we have

and

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

The predicted change in the mean antigenicity was plotted by integrating equation (27). As initial condition, we chose the time when a seed of second peak in the range x > x_c first appeared, and then computed the mean trait as

In the case of Fig. 2, where β = 2, γ + α = 0.6, D = 0.001 and ω = 2, the final size of epidemic for the primary outbreak, defined as equation (7), was ψ = 0.959, and the critical antigenic distance for the increase of pathogen variant obtained from equation (26) was x_c = 2.795. The initial conditions for the oligomorphic dynamics (equation 27) for the second morph were then (p_1left( {t_0} right) = 1.6 times 10^{ – 8}), (bar x_1left( {t_0} right) = 3.239), (V_1left( {t_0} right) = 0.2675) at t₀ = 41. In Fig. 2, the predicted trajectory for the mean antigenicity (equation 28) is plotted as a red curve, together with the mean antigenicity change observed in simulation (blue curve).

Accuracy of predicting the antigenicity with OMD and the timing of the second outbreak

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

We divided the antigenicity space into two at x = x_c, above which the pathogen can increase under the given susceptibility profile after the primary outbreak, but below which the pathogen cannot increase. We then took the relative frequencies of pathogens above x_c and below x_c, and the conditional mean and variance in these separated regions to set the initial frequencies, means and variances of the morphs at time t₀ when we started integrating the oligomorphic dynamics to predict the second outbreak:

We then compared the trajectory for mean antigenicity change observed in simulation (blue curve in Fig. 2) and the predicted trajectory (red curve in Fig. 2) for mean antigenicity (equation 28) by integrating oligomorphic dynamics (equation 27) with the initial condition (equation 29) at time t = t₀. Extended Data Fig. 2 shows how the accuracy of prediction, measured by the Kullback–Leibler divergence between these two trajectories, depends on the timing t₀ chosen for the prediction. The second outbreak occurs around t = 54.6, where mean antigenicity jumps from around 0 to around 5. The prediction with OMD is accurate if it is made for t₀ > 40. Figure 2 is drawn for t₀ = 41 where the second peak is about to emerge (see Extended Data Fig. 2). Even for the latest prediction for t₀ = 51 in Extended Data Fig. 2, the morph frequency of the emerging second morph was only 0.3% off, so the prediction is still worthwhile to make.

Extended Data Fig. 2 shows that the prediction power is roughly constant (albeit with a wiggle) for (5 < t_0 < 30) (the predicted timings are 10–15% longer than actual timing for (5 < t_0 < 30)), and steadily improved for t₀ > 30. When the prediction was made very early (t₀ < 5), the deviations were larger.

OMD for the joint evolution of antigenicity and virulence

Let s(x) be the susceptibility of the host population against antigenicity x. A specific susceptibility profile is given by equation (12), with cross-immunity function σ(x) and the final size ψ₀ of epidemic of the primary outbreak. Note that, as above, the susceptibility profile is, in general, a function of time. The density (I(x,alpha )) of hosts infected by a pathogen of antigenicity x and virulence α changes with time, when rare, as

The change in the frequency (phi left( {x,alpha } right) = Ileft( {x,alpha } right)/{int!!!!!int} I left( {x,alpha } right)dxdalpha) of a pathogen with antigenicity x and virulence α follows

where

is the fitness of a pathogen with antigenicity x and virulence α and (bar w = {int!!!!!int} w left( {x,alpha } right)dxdalpha) is the mean fitness.

We decomposed the joint frequency distribution ϕ(x, α) of the viral quasi-species as (oligomorphic decomposition):

where ϕ_i(x, α) is the joint frequency distribution of antigenicity x and virulence α in morph i (({int!!!!!int} {phi _idxdalpha = 1})) and p_i is the relative frequency of morph i ((mathop {sum}nolimits_i {p_i = 1})). The frequency of morph i then changes as

where (bar w_i = {int!!!!!int} w left( {x,alpha } right)phi _ileft( {x,alpha } right)dxdalpha) is the mean fitness of morph i.

Assuming that the traits are distributed narrowly around the morph means (bar x_i = {int!!!!!int} x phi _ileft( {x,alpha } right)dxdalpha) and (bar alpha _i = {int!!!!!int} alpha phi _i(x,alpha )dxdalpha), so that (xi _i = x – bar x_i = O({it{epsilon }})) and (zeta _i = alpha – bar alpha _i = O({it{epsilon }})) where ({it{epsilon }}) is a small constant, we expanded the fitness w(x, α) around the means (bar x_i) and (bar alpha _i) of morph i,

Substituting this and

into equation (34), we obtained

where (w_i = wleft( {bar x_i,bar alpha _i} right)), (left( {frac{{partial w}}{{partial x}}} right)_i = frac{{partial w}}{{partial x}}left( {bar x_i,bar alpha _i} right)), (left( {frac{{partial w}}{{partial alpha }}} right)_i = frac{{partial w}}{{partial alpha }}left( {bar x_i,bar alpha _i} right)), (left( {frac{{partial ^2w}}{{partial x^2}}} right)_i = frac{{partial ^2w}}{{partial x^2}}left( {bar x_i,bar alpha _i} right)), (left( {frac{{partial ^2w}}{{partial xpartial alpha }}} right)_i = frac{{partial ^2w}}{{partial xpartial alpha }}left( {bar x_i,bar alpha _i} right)) and (left( {frac{{partial ^2w}}{{partial alpha ^2}}} right)_i = frac{{partial ^2w}}{{partial alpha ^2}}left( {bar x_i,bar alpha _i} right)) are fitness and its first and second derivatives evaluated at the mean traits of morph i, and

are within-morph variances and covariance of the traits of morph i. Here (E_ileft[ {fleft( {x,alpha } right)} right] = {int!!!!!int} f left( {x,alpha } right)phi _ileft( {x,alpha } right)dxdalpha) denotes taking expectation of a function f with respect to the joint trait distribution (phi _i(x,alpha )) of morph i.

Substituting equation (36) into the change in the mean antigenicity of morph i

we obtained

Similarly, the change in the mean virulence of morph i was expressed as

Equations (38) and (39) from the mean trait change was summarized in a matrix form as

where

is the variance-covariance matrix of the morph i.

Substituting equation (36) into the right-hand side of the change in variance of antigenicity of morph i,

and those in the change in the other variance and covariance, we obtained

If we assume that antigenicity and virulence within a morph follow a 2D Gaussian distribution for given means, variances and covariance, we should have (E_i(xi _i^4) = 3left( {V_i^x} right)^2,E_i(xi _i^3zeta _i) = 3V_i^xC_i), (E_i(xi _i^2zeta _i^2) = V_i^xV_i^alpha + 2C_i^2), (E_i(xi _izeta _i^3) = 3V_i^alpha C_i) and (E_i(zeta _i^4) = 3left( {V_i^alpha } right)^2), and hence

Equations (43) and (44) were rewritten in a matrix form as

where

is the Hessian of the fitness function of morph i.

In our equation (30) of the joint evolution of antigenicity and virulence of a pathogen after its primary outbreak, the fitness function is given by (w(x,alpha ) = beta (alpha )s(x) – alpha ,) and hence (w_i = beta left( {bar alpha _i} right)sleft( {bar x_i} right) – bar alpha _i), (left( {frac{{partial w}}{{partial x}}} right)_i = beta left( {bar alpha _i} right)sprime left( {bar x_i} right)), (left( {frac{{partial w}}{{partial alpha }}} right)_i = beta prime left( {bar alpha _i} right)sleft( {bar x_i} right) – 1), (left( {frac{{partial ^2w}}{{partial x^2}}} right)_i = beta left( {bar alpha _i} right)sprimeprime left( {bar x_i} right)), (left( {frac{{partial ^2w}}{{partial xpartial alpha }}} right)_i = beta prime left( {bar alpha _i} right)sprime left( {bar x_i} right)), (left( {frac{{partial ^2w}}{{partial xpartial alpha }}} right)_i = beta prime left( {bar alpha _i} right)sprime left( {bar x_i} right)) and (left( {frac{{partial ^2w}}{{partial alpha ^2}}} right)_i = beta primeprime left( {bar alpha _i} right)sleft( {bar x_i} right)), where a prime on β(α) and s(x) denotes differentiation by α and x, respectively. Substituting these into the dynamics for morph frequencies (equation 35), for morph means (equations 38 and 39), and for within-morph variance and covariance (equations 43–45), we obtained

Equations (48)–(53) describe the oligomorphic dynamics of the joint evolution of antigenicity and virulence of a pathogen for a given host susceptibility profile s(x) over pathogen antigenicity.

Of particular interest is whether antigenicity or virulence evolve faster when they jointly evolve than when they evolve alone. After the primary outbreak at a given antigenicity, for example x = 0, the susceptibility s(x) of the host population increases due to cross immunity as the distance x > 0 from the antigenicity at the primary outbreak increases. Hence, (sprime left( {bar x_i} right) > 0.) Combining this with the positive trade-off between transmission rate and virulence, we see that (left( {partial ^2w/partial xpartial alpha } right)_i = beta prime (bar alpha _i)sprime (bar x_i) > 0), and then from equation (52), we see that the within-morph covariance between antigenicity and virulence becomes positive starting from a zero initial value:

If all second moments are initially sufficiently small for an emerging morph, a quick look at the linearization of equations (51)–(53) around ((V_i^x,C_i,V_i^alpha ) = (0,0,0)) indicates that both (V_i^x) and (V_i^alpha) become positive due to the random generation of variance by mutation, D_x > 0 and D_α > 0, while the covariance stays close to zero. Then, equation (54) guarantees that the first move of the covariance is from zero to positive, which then guarantees that C_i > 0 for all t. Therefore, the second term in equation (38) is positive until the mean virulence reaches its optimum ((beta prime (alpha )s(x) = 1)). This means that joint evolution with virulence accelerates the evolution of antigenicity. The same is true for virulence evolution: the first term in equation (39) (which denotes the associated change in virulence due to the selection in antigenicity through genetic covariance between them) is positive, indicating that joint evolution with antigenicity accelerates virulence evolution.

Numerical example

Figure 5 shows the oligomorphic dynamics prediction of the emergence of the next variant in antigenicity–virulence coevolution. To make progress numerically, we assumed s(x) to be constant in the following analysis because we are interested in the process between the end of the primary outbreak and the emergence of the next antigenicity–virulence morph. The partial differential equations for the density of host S(t,x) susceptible to the antigenicity variant x at time t, and the density of hosts infected by the pathogen variant with antigenicity x and virulence α are

with the boundary conditions (left( {partial S/partial x} right)left( {t,0} right) = left( {partial S/partial x} right)left( {t,x_{{{{mathrm{max}}}}}} right) = 0), (left( {partial I/partial x} right)left( {t,0,alpha } right) = left( {partial I/partial x} right)left( {t,x_{{{{mathrm{max}}}}},0} right) = 0), (left( {partial I/partial x} right)left( {t,x,alpha _{{{{mathrm{min}}}}}} right) = left( {partial I/partial x} right)left( {t,x,alpha _{{{{mathrm{max}}}}}} right) = 0), and the initial conditions (Sleft( {0,x} right) = 1) and (Ileft( {0,x,alpha } right) = {it{epsilon }}delta left( x right)delta left( alpha right)), where (delta ( cdot )) is the delta function and ({it{epsilon }} = 0.01). The trait space is restricted in a rectangular region: (0 < x < x_{{{{mathrm{max}}}}} = 300) and (alpha _{{{{mathrm{min}}}}} = 0.025 < alpha < 10 = alpha _{{{{mathrm{max}}}}}). Oligomorphic dynamics prediction for the joint evolution of antigenicity and virulence was applied for the next outbreak after the outbreak with the mean antigenicity around x = 108 at time t = 102. The susceptibility of the host to antigenicity variant x at t₀ = 104.8 after the previous outbreak peaked around time t = 102 came to an end is

This susceptibility profile remained unchanged until the next outbreak started, and hence the fitness of a pathogen variant with antigenicity x and virulence α is given by

We bundled the pathogen variants into two morphs at time t₀ at the threshold antigenicity x_c, above which the net growth rate of the pathogen variant under the given susceptibility profile s(x) and the mean antigenicity become positive:

The initial frequency and the moments of the two morphs, the variant 0 with (x < x_c) and the variant 1 with (x > x_c) were then calculated respectively from the joint distribution (I(t_0,x,alpha )) in the restricted region (left{ {left( {x,alpha } right);0 < x < x_c,alpha _{{{{mathrm{min}}}}} < alpha < alpha _{{{{mathrm{max}}}}}} right}) and that in the restricted region (left{ {left( {x,alpha } right);x_c < x < x_{{{{mathrm{max}}}}},alpha _{{{{mathrm{min}}}}} < alpha < alpha _{{{{mathrm{max}}}}}} right}). The frequency p₁ of morph 1 (the frequency of morph 0 is given by (p_0 = 1 – p_1)), the mean antigenicity (bar x_i) and mean virulence (bar alpha _i) of morph i, and the variances and covariance, (V_i^x), and (V_i^alpha) C_i of morph i (i = 0,1) follow equations (48)–(53), where the dynamics for the morph frequency (equation 48) is simplified in this two-morph situation as

with (p_0left( t right) = 1 – p_1(t)). This is iterated from (t = t_0 = 104.8) to (t_e = 107). The frequency p₁ of the new morph, the population mean antigenicity (bar x = p_0bar x_0 + p_1bar x_1), virulence (bar alpha = p_0bar alpha _0 + p_1bar alpha _1), variance in antigenicity (V_x = p_0V_0^x + p_1V_1^x), covariance between antigenicity and virulence (C = p_0C_0 + p_1C_1), and variance in virulence (V_alpha = p_0V_0^alpha + p_1V_1^alpha) are overlayed by red thick curves on the trajectories of moments observed in the full dynamics (equation 55).

In Fig. 5a, the dashed vertical line represents the threshold antigenicity x_c, above which (R_0 = beta s(x)/(gamma + bar alpha ) > 1) at (t = t_s = 104.8), where oligomorphic dynamics prediction was attempted. Two morphs were then defined according to whether or not the antigenicity exceeded a threshold x = x_c: the resident morph (morph 1) is represented as the dense cloud to the left of x = x_c and the second morph (morph 2) consisting of all the genotypes to the right of x = x_c with their R₀ greater than one. The within-morph means and variances were then calculated in each region. The relative total densities of infected hosts in the left and right regions defined the initial frequency of two morphs in OMD. A 2D Gaussian distribution was assumed for within-morph trait distributions to have the closed moment equations as previously explained. Using these as the initial means, variances, covariances of the two morphs at t = t_s, the oligomorphic dynamics for 11 variables (relative frequency of morph 1, mean antigenicity, mean virulence, variances in antigenicity and virulence and their covariance in morphs 0 and 1) was integrated up to t = t_e. The results are shown as red curves in Fig. 5c–h, which are compared with the simulation results (blue curves).

Fig. 5c–e respectively show the change in total infected density, mean antigenicity and mean virulence. Red curves show the predictions by oligomorphic dynamics from the initial moments of each morph at t = t_s to the susceptibility distribution (s(x) = S(t_s,x)), which are compared with the simulation results (blue curves). The OMD-predicted mean antigenicity, for example, is defined as

where p₁(t) is the frequency of morph 1, (bar x_0) and (bar x_1) are the mean antigenicities of morphs 0 and 1.

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

where (C_0(t)) and (C_1(t)) are the antigenicity–virulence covariances in morphs 0 and 1, and (bar alpha _0(t)) and (bar alpha _1(t)) are the mean virulence of morphs 0 and 1.

Selection for maximum growth rate

We next show that a pathogen that has the strategy of maximizing growth rate in a fully susceptible population is evolutionarily stable in the presence of antigenic escape.

At stationarity, the travelling wave profiles of (hat I(z)) and (hat S(z)) along the moving coordinate, (z = x – vt), that drifts constantly to the right with the speed v are defined as

with (hat Ileft( { – infty } right) = hat Ileft( infty right) = 0), (hat Sleft( infty right) = 1).

Let j(t,x) be the density of a mutant pathogen variant, with virulence α′ and transmission rate β′, that is introduced in the host population where the resident variant is already established (equation 50). For the initial transient phase in which the density of mutants is sufficiently small, we have an equation for the change in (Jleft( {t,z} right) = j(t,x)):

with the initial condition (Jleft( {0,z} right) = {it{epsilon }}delta (z)), where ({it{epsilon }}) is a small constant and (delta ( cdot )) is Dirac’s function.

Consider a system

with (wleft( {0,z} right) = Jleft( {0,z} right) = {it{epsilon }}delta (z)). Noting that (hat Sleft( z right) < 1), we have (Jleft( {t,z} right) le w(t,z)) for any (t > 0) and (z in {Bbb R}) from the comparison theorem. The solution to equation (52) is

where (rprime = beta prime – left( {gamma + alpha prime } right)). This follows by noting that (wleft( {t,x} right)e^{ – rprime t}) follows a simple diffusion equation (partial w/partial t = Dpartial ^2w/partial x^2). By rearranging the exponents of equation (53),

where

Here (vprime = 2sqrt {rprime D}) is the asymptotic wave speed if the mutant variant monopolizes the host population. Therefore, if (vprime < v), then (a < 0), and hence (w(t,z)) for a fixed z converges to zero as t goes to infinity; this, in turn, implies that (J(t,z)) converges to zero because (Jleft( {t,z} right) le wleft( {t,z} right)) for all t and z. Therefore, we conclude that any mutant that has a slower wave speed than the resident can never invade the population, implying that a variant that has the maximum wave speed (v = 2sqrt {rD}) is locally evolutionarily stable.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com