Evolution into chaos - Implications of the trade-off between transmissibility and immune evasion

Abstract

Predicting viral evolution presents a significant challenge and is a critical public health priority. In response to this challenge, we develop a novel model for viral evolution that considers a trade-off between immunity evasion and transmissibility. The model selects for a new strain with the highest invasion fitness, taking into account this trade-off. When the dominant strain of the pathogen is highly transmissible, evolution tends to favor immune evasion, whereas for less contagious strains the direction of evolution leads toward increasing transmissibility. Assuming a linear functional form of this trade-off, we can express the long-term evolutionary patterns following the emergence of subsequent strains by a non-linear difference equation. We provide sufficient criteria for when evolution converges, and successive strains exhibit similar transmissibility. We also identify scenarios characterized by a two-periodic pattern in upcoming strains, indicating a situation where a highly transmissible but not immune-evasive strain is replaced by a less transmissible but highly immune-evasive strain, and vice versa, creating a cyclic pattern. Finally, we show that under certain conditions, viral evolution becomes chaotic and thus future transmissibilites become unpredictable in the long run. Visualization via bifurcation diagrams elucidates our analytical findings, revealing complex dynamic behaviors that include the presence of multiple periodic solutions and extend to chaotic regimes. Our analysis provides valuable insights into the complexities of viral evolution in the light of the trade-off between immune evasion and transmissibility.

Keywords: Chaos; Global stability; Immune evasion; Periodic solutions; Trade-off; Transmissiblity; Viral evolution.

Fig. 1

The model outline flowchart: Beginning with the SIR model with a transmission rate β, while the system is in an endemic steady state (red). New strains with immune-evasive property p, where p ∈ [0, 1], emerge. Assuming a trade-off between transmissibility and immune evasion, denoted by f(p), the transmission rates of these new strains are β_v = β + f(p) (blue). As the evolution tends to favor the strain with maximal invasion fitness (maximizing the invasion reproduction number), we consider the newly emerging strain which is the fittest with a transmission rate β_v = β + f(p^max) (green). The system transitions to a new endemic steady state characterized by a novel transmission rate determined by the trade-off function. The process returns to the initial step and the evolutionary process is continued.

Fig. 2

Panel (a) shows the invasion reproduction number with respect to p for three distinct values of the transmission rate parameter β, 0.26, 0.296, and 0.35 (curves red, green, and blue respectively). Each value corresponds to distinct behaviors exhibited by the invasion reproduction function when its maximum value occurs at p^max = 0, p^max = 1, and p^max ∈ (0, 1), respectively. Panel (b) shows the function p^max(β) with respect to transmission rates of the resident strain. When the transmission rate β surpasses ${\tilde{\beta }}_2$, the maximum invasion reproduction number exhibits an increasing trend with maximum at p = 1(p^max = 1), and for β close enough to γ + μ, p^max = 0. In both panels, μ = 3.5 ⋅ 10⁻⁵, a = 0.02, b = 0.1 and γ = 0.2.

Fig. 3

Graph of g(β) = β + a − bp^max plotted with the enveloping function ϕ(β) = 2β∗ − β (red line) and identity line β_v = β (dashed line). In panel (a), the parameter b = a + 0.05 is chosen such that $g({\tilde{\beta }}_1)>g({\tilde{\beta }}_2)$, and in panel (b), $g({\tilde{\beta }}_1)<g({\tilde{\beta }}_2)$ for b = 0.55. In this example, a = 1.6(γ + μ) > μ + γ, thus global asymptotic stability follows from Theorem 3.2. Other parameters are γ = 0.2, μ = 3.5 ⋅ 10⁻⁵.

Fig. 4

Bifurcation diagram for equations (3.2) and (3.1). The system was iterated for n = 1000 steps from ten initial values for each b between 0.02 and 0.21 with step-size 0.0005, and the last 50 iterations are displayed in the plot. Here, μ = 3.5 ⋅ 10⁻⁵, γ = 0.2, and a = 0.02.

Fig. 5

The cobweb plots correspond to three distinct behaviors of (3.2). μ = 3.5 ⋅ 10⁻⁵, γ = 0.2, and a = 0.02.

Fig. 6

The figure delineates four distinct regions, each corresponding to a unique dynamical behavior of the system as defined by the difference equation (3.2). The regions colored in blue and yellow represent domains where the fixed point exhibits global and local asymptotic stability, respectively. The green regions are zones of instability, as given in Theorem 3.2. The area depicted in gray encapsulates chaotic dynamics of (3.2). Here γ = 0.2 and μ = 3.5 ⋅ 10⁻⁵.

Fig. 7

Prevalence dynamics of infection under varying transmission rate scenarios. The initial value of S, I, and R are 0.999, 0.001, and 0 respectively; μ = 3.5 ⋅ 10⁻⁵, γ = 0.2.