Epidemic spreading under pathogen evolution
Xiyun Zhang, Zhongyuan Ruan, Muhua Zheng, Jie Zhou, Stefano Boccaletti, Baruch Barzel
EEpidemic spreading under pathogen evolution
Xiyun Zhang, , ∗ Zhongyuan Ruan, Muhua Zheng, Jie Zhou, Stefano Boccaletti , , , , † & Baruch Barzel , , †
1. Department of Physics, Jinan University, Guangzhou, Guangdong 510632, China2. Institute of Cyberspace Security, Zhejiang University of Technology, Hangzhou, Zhejiang 310023, China3. School of Physics and Electronic Engineering, Jiangsu University, Zhenjiang, Jiangsu, 212013, China4. School of Physics and Electronic Science, East China Normal University Shanghai 200241, China5. CNR - Institute of Complex Systems, Via Madonna del Piano 10, I-50019 Sesto Fiorentino, Italy6. Unmanned Systems Research Institute, Northwestern Polytechnical University, Xi’an 710072, China7. Moscow Institute of Physics and Technology (National Research University), 9 Institutskiy per., Dolgo-prudny, Moscow Region, 141701, Russian Federation8. Universidad Rey Juan Carlos, Calle Tulip´an s/n, 28933 M´ostoles, Madrid, Spain9. Department of Mathematics, Bar-Ilan University, Ramat-Gan, 5290002, Israel10. Gonda Multidisciplinary Brain Research Center, Bar-Ilan University, Ramat-Gan, 5290002, Israel * Correspondence : [email protected] † These Authors equally contributed to the manuscript
Battling a widespread pandemic is an arms race between our mitigation efforts, e.g. , social distancing or vaccination, and the pathogen’s evolving persistence. Thisis being observed firsthand during the current COVID-19 crisis, as novel mutationschallenge our global vaccination race. To address this, we introduce here a generalframework to model epidemic spreading under pathogen evolution, finding thatmutations can fundamentally alter the projection of the spread. Specifically, wedetect a new pandemic phase - the mutated phase - in which, despite the fact thatthe pathogen is initially non-pandemic ( R < ), it may still spread due to theemergence of a critical mutation. The boundaries of this phase portray a balancebetween the epidemic and the evolutionary time-scales. If the mutation rate is toolow, the pathogen prevalence decays prior to the appearance of a critical mutation.On the other hand, if mutations are too rapid, the pathogen evolution becomesvolatile and, once again, it fails to spread. Between these two extremes, however,we observe a broad range of conditions in which an initially sub-pandemic pathogenwill eventually gain prevalence. This is especially relevant during vaccination, whichcreates, as it progresses, increasing selection pressure towards vaccine-resistance.To overcome this, we show that vaccination campaigns must be accompanied byfierce mitigation efforts, to suppress the potential rise of a resistant mutant strain. Evolutionary time-scales are often considered to be vast, occurring gradually over the courseof millions of years. However, if prevalent enough, a species may undergo even rare mutationsat relatively short time-scales. This is especially relevant during the course of a widespreadand prolonged pandemic. The global spread ensures a sufficiently large pool of pathogens formutations to occur, and on top of that, the long duration of the pandemic affords the pathogenssufficient time to evolve.Such troubling scenario is currently unfolding in the case of COVID-19, where novel mutationsof the SARS-CoV-2 virus continue to challenge our mitigation efforts, but are equally relevantin other infections, such as influenza A, where they are, in fact, part of our suppression routine,1 a r X i v : . [ q - b i o . P E ] F e b aving to distribute a dedicated vaccine in each yearly cycle. Another notable example isnorovirus, whose enhanced transmission, likely due to mutation, led to an observable spike ingastric flu patients in England and Wales from 1991 to 2006, and finally, beyond viruses,artemisinin-resistance, a parasite mutation, rendered void the common treatment of malaria inAfrica. The common approach for tracking the spread of evolving pathogens is to introduce severalcompeting strains and extract their interacting contagion dynamics.
This captures the pat-terns of spread of already evolved pathogens, overlooking the dynamics, and most importantly,the time-scales, of the evolution itself. Indeed, in an ongoing pandemic, mutations represent agradual random process, in which an originally unfit pathogen mutates step-by-step via a seriesof small changes, until reaching a critical mutation that allows it to efficiently spread. Suchprocess may take a significant amount of time, and, in some cases, the disease may taper offbefore such critical mutation has the opportunity to take over.Another crucial aspect, absent in the present treatment based on pre-mutated strains, is the factthat pathogen evolution is responsive . As we tighten our mitigation, either through prophylacticmeasures or via pharmaceutical interventions, we induce a selective pressure for mitigationresistance. For example, if one enforces social distancing to push the reproduction rate R belowthe pandemic threshold, the pathogen becomes naturally pressured towards higher transmissi-bility. Similarly, if one employs therapeutic treatment to expedite recovery, natural selectionwill push the pathogen to higher drug-persistence.To address this, we introduce here an evolving pathogen model, designed to capture the delicateinterplay between the pathogen’s spread and its developing fitness. The evolution, a randomwalk in fitness space , is driven by the pathogen’s mutation rate. At the same time the naturalselection , in which the fitter strains proliferate, is pushed by the epidemiological parameters,characterizing how fast a mutated strain propagates. Together, we identify a rather broad setof conditions - the mutated phase - in which a non-pandemic pathogen will eventually reach anevolved pandemic state.We find that besides the classic epidemiological parameters, i.e . infection/recovery rates, twoadditional components factor in - the mutation rate governing the evolutionary time-scales, andthe number of infected individuals, which determines the likelihood of a critical mutation tooccur within the relevant time-frame. Therefore, the current prevalence ρ ( t ) of the pathogenhas direct impact on its projected spread. This has significant implications pertaining to our twomain mitigation strategies • Social distancing suppresses the reproduction number R to belowthe pandemic threshold. However, if many have already been infected, i.e . ρ ( t ) is large, then astricter suppression may be required to avoid the emergence of a critical mutation. This indicatesthat the projection of the spread, and hence its mitigation, depends on the current prevalence ρ ( t ) - a hysteresis phenomenon, unobserved in the classic modeling frameworks • Vaccinationcampaigns create strong evolutionary pressure towards a vaccine resistant mutation, whose risk,once again, is directly related to the current pathogen prevalence. Hence, to succeed, we showthat vaccine roll-out, no matter how efficient, must be coupled with fierce suppression via socialdistancing.
Evolving pathogen model
Consider a social network of N individuals, linked through the adjacency matrix A ij , a randomnetwork with average degree k . At t = 0 the network experiences an outbreak, which then2preads via the susceptible-infected-susceptible (SIS) dynamics. In the classic SIS formulation,the projected spread is driven by two time-independent parameters: the recovery rate µ and theinfection rate β , whose ratio R = kβ/µ , the reproduction number , determines the state of thesystem - pandemic ( R ≥
1) or healthy ( R < µ i ( t ) = 1 F i ( t ) µ, (1)where the fitness F i ( t ) stands for the level of mutation of the pathogen carried by individual i at time t , hence the unmutated pathogen has F i (0) = 1. The above equation models the factthat (i) each individual i may carry a distinct version of the virus; (ii) this version my graduallychange in time t due to mutations. The smaller is µ i ( t ), the higher is the transmissibility of thepathogen, as described by the evolving reproduction number R i ( t ) = kβµ i ( t ) = R F i ( t ) . (2)Indeed, a low rate of recovery µ i ( t ) extends the duration of the infectious state, providingindividual i with more opportunities to infect their peers. Hence, as the r.h.s. of (2) indicates,pathogens with increased F i ( t ) exhibit higher reproduction, and therefore they spread moreefficiently than their lower fitness competitors.While mutation may also impact the transmissibility of the pathogen directly (by altering thevalue of β ), in (1) we only evolve µ i ( t ), as indeed only the ratio R i ( t ) affects the spread inthe SIS framework, and not the specific values of µ or β . Therefore, for simplicity, here weonly track the pathogen evolution through µ i ( t ), and its subsequent R i ( t ), and set β stationary.Complementary to this, in Supplementary Section XXX we examine the case of β -mutations.The spread is driven by the following set of transitions: the process of infection between a pairof individuals i and j is modeled by S i + I j A ij β −−−→ I i + I j (3) F i ( t ) = F j ( t ) , (4)in which a susceptible ( S ) individual i interacts with their infected ( I ) neighbor j ( A ij = 1)at rate β . This leads to both individuals becoming infected. The newly infected individual i inherits j ’s pathogen, and hence in (4) we set i ’s fitness at the time of infection equal to that of j . Both fitness parameters, F i ( t ) and F j ( t ) may later change via mutation. Next, we considerthe process of recovery I i µ i ( t ) −−−→ S i , (5)in which an infected individual I i transitions to S i at the evolved recovery rate µ i ( t ) of (1).Finally, the process of mutation follows 3 F i (0) = 1 F i ( t + 1) = max (cid:16) F i ( t ) + δ i ( t ) , (cid:17) , (6)capturing a random walk with variable step size δ i ( t ), i.e . a sequence of random shifts in fitness,caused by small discrepancies in the pathogen’s reproduction. Note that F i ( t ) is prohibited frombecoming negative, as, indeed, a below zero fitness in (2) is meaningless. The case where F i ( t )does approach zero corresponds to µ i ( t ) → ∞ in (1), a limit in which recovery is instantaneous,and hence the pathogen is unfit for reproduction. Such strains will be rapidly eliminated fromthe pathogen pool.The magnitude of each mutation step is extracted from a zero-mean normal distribution, namely δ i ( t ) ∼ N (0 , σ ). Consequently, in the limit where σ = 0, we have δ i ( t ) = 0 at all times,mutations are suppressed, and Eqs. (3) - (6) converge to the classic SIS model, with R i ( t ) = R as the constant reproduction number. In contrast, as σ is increased, significant mutations becomemore frequent and the pathogens rapidly evolve. We therefore vary σ to control the mutationrate of the pathogens.Taken together, our modeling framework accounts for the dynamics of infection and recovery(SIS) under the effect of pathogen mutation. As the spread progresses, pathogens evolve via Eq.(6), blindly altering their epidemiological parameters at random. Natural selection, however,will favor the positive mutations, in which δ i ( t ) >
0. Indeed, such mutations lead to higherfitness, reducing the recovery rate µ i ( t ), and consequently increasing R i ( t ). Such pathogens,with increased R i ( t ), will proliferate more rapidly, and will eventually dominate the population. Critical mutation . Consider an outbreak of a pathogen with R < i.e . below the epidemicthreshold. This can be either due to the pathogen’s initial sub-pandemic parameters, or aresult of mitigation, e.g ., social distancing to reduce β . In the classic SIS formulation, suchpathogen with fail to penetrate the network. However, in the presence of mutations ( σ >
0) thepathogen may potentially undergo selection, reach R i ( t ) >
1, and from that time onward beginto proliferate. This represents a critical mutation , which, using (2), translates to F c = 1 R , (7)the critical fitness that, once crossed, may lead an initially non-pandemic pathogen to becomepandemic. The smaller is R the higher is F c , as, indeed, weakly transmissible pathogens requirea long evolutionary path to reach pandemic spread.Next, we analyze the spreading patterns of our evolving pathogens, seeking the conditions forthe appearance of the critical mutation. Phase-diagram of evolving pathogens
To examine the behavior of (1) - (6) we constructed an Erd˝os-R´enyi (ER) network with N =5 ,
000 nodes and k = 15, as a testing ground for several epidemic scenarios (Fig. 1). Each scenariois characterized by a different selection of our model’s three epidemiological parameters: µ and β , which determine the pathogen’s unmutated reproduction R , and σ , which controls the rateof mutation. We then follow the spread by measuring the prevalence ρ ( t ), which monitorsthe fraction of infected individuals vs. time. We also track the pathogen’s evolution via the4opulation averaged fitness F ( t ) = (1 /N ) (cid:80) Ni =1 F i ( t ). Pandemic (Fig. 1a, red). In our first scenario we set µ = 0 . β = 8 × − and σ = 10 − . Thiscaptures a pandemic pathogen, which, using k = 15, has R = 1 . >
1, namely it can spreadeven without mutation. Indeed ρ ( t ) rapidly climbs to gain macroscopic coverage, congruent withthe prediction of the classic SIS model, but this time constantly growing, due to the gradual,but continuous, increase in fitness F ( t ). Mutated (Fig. 1c, green). Next we reduce the infection rate to β = 1 . × − , an initialreproduction of R = 0 . <
1. This describes a pathogen whose transmissibility is significantlybelow the epidemic threshold, and therefore, following the initial outbreak we observe a declinein ρ ( t ), which by t ∼
50 almost approaches zero, as the disease seems to be tapering off. In thisscenario, however, we set a faster mutation rate σ = 1. As a result, despite the initial remission,at around t ∼
15, the pathogen undergoes a critical mutation as F ( t ) crosses the critical F c =1 /R = 4 (grey dashed line) and transitions into the pandemic regime. Consequently, ρ ( t )changes course, the disease reemerges and the mutated pathogens successfully spreads. Lethargic (Fig. 1b, blue). We now remain in the sub-pandemic regime, with R = 0 .
25, butwith a much slower mutation rate, set again to σ = 10 − . As above, ρ ( t ) declines, however thepathogen evolution is now too slow - it is lethargic , and cannot reach critical fitness on time.Therefore, the disease fails to penetrate the network, lacking the opportunity for the criticalmutation to occur.Taken together, the dynamics of the spread are driven by three parameters: the initial epidemi-ological characteristics of the pathogen, µ and β , which determine R , and the mutation rate σ , which governs the time-scale for the appearance of the critical mutation. Therefore, to de-termine the conditions for a mutation-driven contagion we investigate the balance between thedecay in ρ ( t ) vs. the gradual increase in F ( t ). The mutated phase
To understand the dynamics of the evolving pathogen model, we show in Supplementary SectionXXX that at the initial stages of the spread, the prevalence ρ ( t ) follows ρ ( t ) = ρ (0) e ξ ( t ) . (8)The time-dependent exponential rate ξ ( t ) is determined by the epidemiological/mutation ratesvia ξ ( t ) = − µ (1 − R ) t + 12 σ µ R t , (9)whose two terms characterize the pre-mutated vs. post-mutated spread of the pathogen. Thefirst term, linear in t , captures the initial patterns of spread, which are determined by theoriginal pathogen parameters, µ, R . For R < a l`a theSIS dynamics in the sub-pandemic regime. At later times, however, as t becomes large, thesecond term begins to dominate, and the exponential decay is replaced by a rapid proliferation,this time driven by the mutation rate σ . The transition between these two behaviors - decay vs.proliferation - occurs at τ c = (cid:112) − R ) / µσ R , which quantifies the anticipated time-scalefor the appearance of the critical mutation F c in (7).This analysis portrays the mutated phase as a balance between two competing time-scales: on5he one hand the exponential decay of the sub-pandemic pathogen, and on the other hand theevolutionary time-scale τ c for the appearance of the critical mutation. For the evolution to win this race the pathogen must not vanish before t = τ c . This imposes the condition (Fig. 2a-c) ρ ( τ c ) ≥ N , (10)ensuring that at τ c there are still one or more individuals hosting the pathogen. Indeed, ρ ( τ c ) < /N indicates that on average , at t = τ c less than a single individual is left in the infected pool.Under this condition, the critical mutation is too late, the spread has already tapered off, andthe exponential spread driven by the positive term in (9) is averted.Taking ρ ( τ c ) from (8), we can now use (10) to express the boundary of the mutated phase,predicting the critical mutation rate as σ c ∼ (cid:32) (cid:112) µ (1 − R ) R (cid:33) I ) , (11)where I = N ρ ( t = 0) is the number of individuals infected at t = 0. Equation (11) describesthe minimal mutation rate required for the pathogen to evolve a pandemic strain. For R = 1it predicts σ c = 0, as such pathogen can indeed spread even without mutation. However, as R is decreased, for example under mitigation, the pathogen prevalence rapidly declines, and henceit must evolve at an accelerated rate to reach critical fitness. This is expressed in (11) by anincreased σ c , which approaches infinity as R → ,
050 realizations of Eqs. (1) -(6), representing different epidemiological scenarios. We varied R from 0 to 1 . i.e . from non-transmissible to highly contagious, and scanned a spectrum of mutation rates from σ = 10 − to σ = 10, spanning four orders of magnitude. Simulating each scenario 50 times we observethe probability P for the disease to spread. This is done by tracking the pathogen’s long-term prevalence ρ = ρ ( t → ∞ ) and counting the realizations in which ρ → ρ >
0. As predicted, we find that the pandemic state, classically observed only at R ≥ R in the presence of sufficiently rapid mutations. This gives rise to the mutated phase (green), in which an initially decaying contagion suddenly turns pandemic. Thetransition between the lethargic and the mutated states (grey zone) is well-approximated by ourtheoretical prediction of Eq. (11), depicted by the black solid line.Equation (11) shows that σ c depends not only on the epidemiological characteristics of thepathogen ( µ, R ), but also on the initial condition, here captured by the number of infectedindividuals I = ρ ( t = 0) N . If I is large the critical rate σ c becomes lower, in effect expandingthe bounds of the mutated phase. To understand this consider the evolutionary paths followedby the pathogens as they reproduce. These paths represent random trajectories in fitness space ,each starting from F i (0) = 1, and with a small probability crossing the critical fitness F c . Themore such attempts are made, the higher the chances that at least one of these paths will besuccessful. Therefore, a higher initial prevalence I of the pathogen increases the probabilityfor the appearance of a critical mutation, enabling a mutated phase under a lower σ . Indeed,in Fig. 2d we find that the phase boundary shifts towards lower σ c as the initial prevalence isincreased (grey shaded lines). Hysteresis . This dependence on I indicates that the transition of Eq. (11) behaves differently6f we approach it from the pandemic state or from the healthy state. To observe this let us fixthe mutation rate at σ = 0 . R , seeking the critical point where thesystem shifts to the mutated phase. This is mapped to a vertical trajectory in the σ, R plane(Fig. 2d, yellow dashed line). At each value of R we instigate an outbreak with ρ (0) = 0 .
2, andobserve its long-term prevalence ρ . For small R this outbreak decays and the system revertsto the healthy state ρ = 0. However, as we transition into the mutated phase, here predictedat R = R High = 0 .
6, the pathogen turns pandemic and its prevalence abruptly changes to ρ ≈ . R slightly below this critical point,for instance, by practicing social distancing to reduce transmission. The challenge is that now,moving in the opposite direction - from large to small R - our initial condition is pandemic,with prevalence of order unity ( ∼ I ∼ N . Under these conditions, Eq. (11)predicts that, for a fixed σ , the critical R is now lower, at R Low = 0 .
35. This results in ahysteresis phenomenon, in which criticality occurs at different points depending on the statefrom which we approach the transition (Fig. 2e).We find, therefore, that pathogen evolution fundamentally changes the phase space of epidemicspreading. First it predicts a broad range of conditions - the mutated phase - in which a sub-pandemic pathogen can gain prevalence. On top of that, it also predicts that this phase exhibitsa discontinuous transition, characterized by hysteresis, a phenomenon unobserved in the classicSIS dynamics, yet congruent with with other models that incorporate feedback betweena pathogen’s prevalence ( ρ ( t )) and its potency ( R i ( t )). These two observations have directimplications on mitigation: • Soft mitigation is risky . Most mitigation strategies seek a minimal approach, aimingto drive R just below unity. This is understandable as (i) major restrictions on socialinteractions are costly and difficult to sustain for extended periods; (ii) having R < ρ ( t ) todecay exponentially towards zero. Our analysis, however, shows that this is insufficient.For R (cid:46) σ c →
0, indicating that even a relatively stable pathogen, with a lowmutation rate, may eventually break through. Our analysis accounts for this effect andthrough Eq. (11), provide us, given σ , the level of tolerable R that sufficiently mitigatesthe mutated phase risk. • The sooner the better . Another common assumption, driven by the classic epidemicphase-diagram, is that the projected state ρ ( t → ∞ ) depends only on R , i.e . the epi-demiological parameters. The current state of the spread ρ (0) plays no role. The observedhysteresis, however, shows that successful mitigation strongly depends on the prevalenceat the time of instigation. If the pathogen has already gained sufficient ground, we willneed to suppress the reproduction number below R Low , namely the lower phase-boundaryin Fig. 2e. It is, therefore, crucial to respond early, and instigate our mitigation when ρ ( t )is still small, eradicating the pandemic before mutations risk its reemergence. Bounded fitness . Our mutation process in Eq. (6) allows the pathogen an unbounded randomwalk in fitness space. In reality, however, there are practical restrictions on fitness, as R i ( t )cannot grow ad infinitum . Therefore, we now consider our evolving pathogen model, substitutingthe mutation in (6) with F i ( t + 1) = min (cid:16) F max , max (cid:0) F i ( t ) + δ i ( t ) , (cid:1)(cid:17) , (12)7n which the pathogen fitness is bounded from above by F max and from below by zero. Setting F max = 20 we now revisit our phase-diagram (Fig. 3a). For small σ , mutations are slow, andthe evolution path is unaffected by the upper bound on F i ( t ). Therefore, we continue to observethe same transition as in the unbounded model of Fig. Fig1d. As we increase σ , however, wewitness a second transition, this time back to the healthy state, indicating that now, mutationsare too rapid. This captures the final phase of our evolving pathogen model - the volatile phase: Volatile (Fig. 3, blue). When the mutation rate is too high the pathogen fitness becomesunstable. On the one hand it can rapidly reach critical fitness, yet, on the other hand, dueto the random nature of its frequent mutations, it fails to sustain this fitness - resulting in anirregular F ( t ), that fluctuates above and bellow the critical F c (Fig. 3c).To gain deeper insight into the volatile phase, consider the natural selection process, here drivenby the reproduction benefit of the fitter strains. This process is not instantaneous, and requiresseveral reproduction instances, i.e . generations, to gain a sufficient spreading advantage. With σ to high, natural selection is confounded, the pathogen shown no consistent gain in fitness and,as Fig. 3c indicates, ρ ( t ) decays exponentially to zero. In Supplementary Section XXX we use atime-scale analysis, similar to the one leading to Eq. (11), to show that the volatile phase occurswhen σ exceeds σ c ∼ F max R − R . (13)This prediction is, indeed, confirmed by our simulated phase diagram (Fig. 3a, black solid line).Together, our phase-diagram illustrates the different forces governing the spread of pathogensin the presence of mutations. While classically, spread is prohibited for R <
1, here weobserve a new, previously undocumented pandemic phase, in which the disease can successfullypermeate despite having an initially low reproduction rate. The conditions for this phaserequire a balance between three separate time-scales: (i) The time for the initial outbreak ρ (0) to reach near zero prevalence τ r ; (ii) The time for the pathogen to evolve beyond criticalfitness τ c ; (iii) The time for the natural selection to lock-in the fitter mutations τ s . Pathogenswith small R , we find, can still spread provided that τ r > τ c > τ s . (14)The l.h.s. of (14) ensures that the pathogen can reach critical fitness before reaching zeroprevalence. This gives rise to the first transition in Eq. (11), between the lethargic and the mutated phases. The r.h.s. of (14) is responsible for the second transition, from mutated to volatile . It ensures that fitter pathogens do not undergo additional mutation before they havetime to proliferate via natural selection. Therefore, we observe a Goldilocks zone , in which themutation rate σ is just right: on the one hand, enabling unfit pathogens to cross the Rubicontowards pandemicity, but on the other hand, avoiding aimless capricious mutations. Mutation risk in vaccine distribution
Vaccination during an ongoing pandemic is, by nature, a competition between the rate of thevaccine roll-out and the spread of the pathogen.
Therefore, na¨ıvely, to win this race all onehas to do is disseminate the vaccine as efficiently as possible, aiming to reach the majority of the8opulation before the pathogen does. This, however, ignores the role of mutations, which maygravely impact even the most efficient vaccination campaign. Such mutations may, generally, beless fit epidemiologically, i.e . have a lower F and consequently a lower R i ( t ). Therefore, absent avaccine, they will be rapidly overcome by the faster spreading pathogen strains. However, oncethe vaccine becomes widespread, resistance, even if less contagious, becomes a highly desirabletrait, and a resistant mutation, if occurs, will take over the population.To examine this we consider a pathogen spreading via the susceptible-infected-recovered (SIR)model, and mutating via (6). However, this time we add two additional components: • Vaccination . The population is vaccinated at a rate ν , quantifying the percentage of the(susceptible) population that receives the vaccine per unit time (day). • Resistance . At each time-step, in addition to the mutations considered in Eq. (6), thepathogen may undergo a vaccine-resistant mutation with probability p . This mutation hasno bearing on its epidemiological parameters µ, β , hence providing no additional spreadingadvantage, other than being resistant to the vaccine. We set p = 5 × − , a rare mutation.In Fig. 4 we track a two vaccination scenarios, capturing slow ( ν = 2 × − ) and rapid ( ν =2 × − ) roll-outs. In the rapid scenario we set a quite optimistic ν , assuming a coverage of ∼ .
2% of the population per day. Still, in both cases we find that vaccination cannot guaranteethe eradication of the pathogen. For example, in our rapid scenario, a resistant pathogen emergesin 44% of realizations, capturing a significant risk for vaccination failure (Fig. 4d). Hence, wefind that rapid vaccination, in and of itself, may be insufficient.The point is that as the vaccine prevalence rises, so does the selective pressure for resistance.Indeed, during the spread, there may be several instances of resistant mutations. But, in apredominantly unvaccinated population, such mutations rarely take over, since they are notnecessarily fitter than all other coexisting strains. However, as the vaccination campaign pro-gresses, and a significant portion of the population is immune to the non-resistant strain, thesesame mutations provide a fitness boost, and hence, when they occur, they will eventually gainprevalence. This is clearly observed in Fig. 4c,f, as the serendipitous resistant mutation (blacknodes) rapidly penetrates the network, rendering the vaccine moot.To mitigate this risk we must strongly suppress the probability of the resistance mutation tooccur. As we have no control over the microscopic biological processes driving the pathogenmutation, our only means to reduce this risk is by decreasing the prevalence of the pathogen ρ ( t ). Indeed, as our analysis has shown, critical mutations are more likely to appear if thepathogen is widespread. Therefore, we find that during a vaccination campaign it is crucialto continue the suppression of the pandemic. This, we emphasize is not just to overcome thespread - that can indeed be accomplished by vaccination alone - but mainly to mitigate the riskof vaccine-resistant mutations.To observe this, in Fig. 4g-l we revisit our vaccination scenarios, but this time reinforced themwith standard mitigation efforts, such as social distancing to reduce R below unity. In case ofslow vaccination we continue to observe a large portion of failed campaigns (Fig. 4g, 52%). Thisis because the slow build-up of immunity afforded the pathogen sufficient time for evolving acritical mutation (red nodes). It then gained prevalence, increasing the chances for the resistancemutation (black nodes). Indeed, in Fig. 4i we observe first the appearance of red nodes, andonly then black nodes, indicating that initially a critical mutation a l`a Eq. (7) boosts the spread,after which resistance is evolved. 9onsequently, our only successful vaccination scenario combines both strategies - mitigationtogether with rapid roll-out. The latter ensures the pathogen has no sufficient time for spread-ing via critical mutation, while the former suppresses ρ ( t ), alleviating the risk of a resistancemutation. As a result, the few sporadic instances of resistance mutations fail to penetrate thenetwork thanks to the mitigation ( R < Discussion and outlook
The phase diagram of epidemic spreading is a crucial tool for forecasting and mitigating pan-demic risks. First, it identifies the relevant control parameters, such as µ, β and k in our SISframework, or additional parameters in more complex contagion processes, whose value deter-mines R and hence the expected the patterns of spread. The phase boundaries, then, help usassess the state of the system - healthy or pandemic - and provide guidelines for our response.For example, social distancing to reduce k , therapeutic treatment to increase µ or mask wearingto suppress β , all aimed to navigate the system’s location along the pandemic phase-diagramtowards the desired healthy state.The common thread binding all of these strategies is the assumption that the epidemiologicalcontrol parameters themselves are constant in time, and hence our intervention must just pushthem beyond the static phase-boundary, form which point on the epidemic will decay towards ρ → µ, β is slow comparedto the epidemic spreading dynamics - as observed in the case of our lethargic phase. However,once the epidemiological parameters can change at sufficient rate, it fundamentally changes the rules of the game . This is because now, not only are the parameters dynamic , but, thanks tonatural selection, they also become responsive . If, for instance, we develop drug-based treatmentto increase the recovery rate µ , we inevitably also generate selection pressure towards drugpersistence. Similarly, if we vaccinate or practice distancing to reduce k, β , we initiate anevolutionary race towards higher transmissibility or vaccine resistance.The result is a complex interplay between the spreading dynamics ( R ), the instantaneousprevalence of the pathogen ( ρ ( t )), and the dynamic evolution of its parameters ( σ ), whichreshapes the pandemic phase diagram. It not only expands the pandemic risk to a range of R <
1, but also predicts an explosive transition pattern, i.e . the hysteresis of Fig. 2e, that isnot observed in standard epidemiological transitions. This altered phase diagram, and its abruptfirst-order like transition, we have shown, has crucial implications pertaining to our mitigationstrategies. Yet, more broadly, as a physical phenomenon, it offers an interesting mechanism forexplosive transitions. Most often, such abrupt phase-shifts are caused by internal suppressionrules, that hold back the transition until it breaks through in an explosive fashion.
Incontrast, here what holds back the transition is the waiting time for the critical mutation. Untilits appearance the system behaves in one way ( R < R > R , I , σ ). This local event thenchanges fundamentally the state of the system - capturing a feedback between the system’s phaseand its intrinsic control parameters. We believe this describes a unique mechanism, inherent tothe basic ingredients of our biological system, reproduction, mutation and selection.10 cknowledgements X.Z. thanks Dr. Xiaobo Chen, Dr. Tingting Shi, Dr. Xing Lu and Prof. Weirong Zhong for usefuldiscussions and supports in numerical calculations. This research was supported by the IsraelScience Foundation (grant No. 499/19) and the Bar-Ilan University Data Science Institute grantfor COVID-19 related research.
Author contribution
X.Z. developed the concept. X.Z., B.B. and S.B. designed the framework. X.Z. and Z.R.performed the numerical simulations. All authors jointly analyze the results and developed theanalytic framework. X.Z., B.B. and S.B. wrote the paper.
Code availability
All code to study and reproduce the results shown here will be made freely available online uponpublication. 11 eferencesReferences [1] H. Yao, X. Lu, Q. Chen, K. Xu, Y. Chen, L. Cheng, F. Liu, Z. Wu, H. Wu, C. Jin,M. Zheng, N. Wu, C. Jiang and L. Li, Patient-derived mutations impact pathogenicity ofSARS-CoV-2, medRxiv 2020.04.14.20060160 (2020).[2] Y. Konno, I. Kimura, K. Uriu, M. Fukushi, T. Irie, Y. Koyanagi, S. Nakagawa and K. Sato,SARS-CoV-2 ORF3b is a potent interferon antagonist whose activity is further increasedby a naturally occurring elongation variant, bioRxiv 2020.05.11.088179 (2020).[3] J. Hu, C. He, Q. Gao, G. Zhang, X. Cao, Q. Long, H. Deng, L. Huang, J. Chen, K. Wang,N. Tang and A. Huang, The D614G mutation of SARS-CoV-2 spike protein enhances viralinfectivity and decreases neutralization sensitivity to individual convalescent sera, bioRxiv2020.06.20.161323 (2020).[4] B. Korber, W. Fischer, S. Gnanakaran, H. Yoon, J. Theiler, W. Abfalterer, N. Hengart-ner, E. Giorgi, T. Bhattacharya, B. Foley, K. Hastie, M. Parker, D. Partridge, C. Evans,T. Freeman, T. de Silva, C. McDanal, L. Perez, H. Tang, A. Moon-Walker, S. Whelan,C. LaBranche, E. Saphire, D. Montefiori, on behalf of the Sheffield COVID-19 GenomicsGroup, Tracking changes in SARS-CoV-2 Spike: evidence that D614G increases infectivityof the COVID-19 virus, Cell (2020).[5] N.D. Grubaugh, W.P. Hanage, A.L. Rasmussen, Making sense of mutation: what D614Gmeans for the COVID-19 pandemic remains unclear, Cell (2020).[6] S. Ojosnegros and N. Beerenwinkel, Models of RNA virus evolution and their roles in vaccinedesign, Immunome Research 6(Suppl 2):S5 (2010).[7] A. Nougair`ede, R. N. Charrel and D. Raoult, Models cannot predict future outbreaks:A/H1N1 virus,the paradigm, Eur J Epidemiol 26:183, 186 (2011).[8] T. Bedford, M. A. Suchard, P. Lemey, G. Dudas, V. Gregory, A. J. Hay, J. W. McCauley,C. A. Russell, Derek. J. Smith and A. Rambaut, Integrating influenza antigenic dynamicswith molecular evolution, eLife 3:e01914 (2014).[9] M. I. Nelson and E. C. Holmes, The evolution of epidemic influenza, Nature Reviews Ge-netics 8, pages196,205 (2007).[10] E. Ghedin, N. A. Sengamalay, M. Shumway, J. Zaborsky, T. Feldblyum, V. Subbu, D. J.Spiro, J. Sitz, H. Koo, P. Bolotov, D. Dernovoy, T. Tatusova, Y. Bao, K. St George, J.Taylor, D. J. Lipman, C. M. Fraser, J. K. Taubenberger and S. L. Salzberg, Large-scalesequencing of human influenza reveals the dynamic nature of viral genome evolution, Nature437, pages 1162, 1166 (2005).[11] D. J.D. Earn, J. Dushoff and S. A. Levin, Ecology and evolution of the flu, TRENDS inEcology & Evolution 17, Pages 334-340 (2002).[12] B. Lopman, M. Zambon and D. W. Brown, The Evolution of Norovirus, the Gastric Flu,PLoS Medicine 5(2): e42. (2008). 1213] A. Uwimana et al., Emergence and clonal expansion of in vitro artemisinin-resistant Plas-modium falciparum kelch13 R561H mutant parasites in Rwanda, Nature Medicine (2020).[14] J. Birnbaum et al., A Kelch13-defined endocytosis pathway mediates artemisinin resistancein malaria parasites, Science 367, 51-59 (2020).[15] R. Pastor-Satorras and A. Vespignani, Epidemic spreading in scale-free networks, Phys.Rev. Lett. 86, 3200 (2001).[16] Z. Liu and B. Hu, Epidemic spreading in community networks, Europhys. Lett. 72, 315(2005).[17] Z. Ruan, M. Tang, and Z. Liu, Epidemic spreading with information-driven vaccination,Phys. Rev. E 86, 036117 (2012).[18] R. Pastor-Satorras, C. Castellano, P. V. Mieghem, and A. Vespignani, Epidemic processesin complex networks, Rev. Mod. Phys. 87, 925 (2015).[19] H. Alexander, T. Day, Risk factors for the evolutionary emergence of pathogens. J. R. Soc.Interface 7, 1455?1474 (2010).[20] F. Uekermann and K. Sneppen, Spreading of multiple epidemics with cross immunization.Phys. Rev. E, 86, 036108 (2012).[21] L. H´ebert-Dufresne, O. Patterson-Lomba, G. M. Goerg, and B. M. Althouse, PathogenMutation Modeled by Competition Between Site and Bond Percolation, Phys. Rev. Lett.110, 108103 (2013).[22] C. Alexandrou, V. Harmandaris, A. Irakleous, G. Koutsou and N. Savva, Modeling theevolution of COVID-19, arXiv:2008.03165 (2020).[23] R. Eletreby et al., The effects of evolutionary adaptations on spreading processes in complexnetworks, Proc. Natl. Acad. Sci. U. S. A. 117,5664 (2020).[24] R.M. Anderson, H. Heesterbeek, S. Klinkenberg and T.D. Hollingsworth, How will country-based mitigation measures influence the course of the COVID-19 epidemic? The Lancet,395,10228 (2020).[25] C. Shen and Y. Bar-Yam. COVID-19: How to win (2020).[26] D. Meidan, N. Schulmann, R. Cohen, et al. Alternating quarantine for sustainable epidemicmitigation. Nature Communications 12,220 (2021).[27] D. Pride and C. Ghose, Meet the trillions of viruses that make up your vi-rome, The Conversation, https://theconversation.com/meet-the-trillions-of-viruses-that-make-up-your-virome-104105[28] J. L. Mokili, F. Rohwer and B. E. Dutilh, Metagenomics and future perspectives in virusdiscovery, Current Opinion in Virology 2(1), 63-77 (2012).[29] J. G´omez-Garde˝nes, L. Lotero, S. N. Taraskin, and F. J. P´erez-Reche, Explosive contagionin networks, Sci. Rep. 6, 19767 (2016). 1330] Q. Liu, W. Wang, M. Tang, T. Zhou, and Y. Lai, Explosive spreading on complex networks:The role of synergy, Phys. Rev. E 95, 042320 (2017).[31] J. Wu, M. Zheng, K. Xu, and C. Gu, Effects of two channels on explosive informationspreading, Nonlinear Dynamics, 99, 2387-2397, (2020).[32] L. H´ebert-Dufresne, S. V. Scarpino, and J. Young, Macroscopic patterns of interactingcontagions are indistinguishable from social reinforcement, Nature Physics 16, 426 (2020).[33] T. Gross, C. J. D. Lima, and B. Blasius, Epidemic dynamics on an adaptive network, Phys.Rev. Lett. 96, 208701 (2006).[34] L. B¨ottcher, O. Woolley-Meza, E. Goles, D. Helbing, and H. J. Herrmann, Connectivitydisruption sparks explosive epidemic spreading, Phys. Rev. E 93, 042315 (2016).[35] X. Zhang, Z. Ruan, M. Zheng, B. Barzel, S. Boccaletti, Epidemic spreading under infection-reduced-recovery, Chaos, Solitons and Fractals 140, 110130 (2020).[36] J.M. Epstein, Modelling to contain pandemics. Nature 460,687 (2009).[37] U. Harush, and B. Barzel, Dynamic patterns of information flow in complex net- works,Nat. Comm. 8, 2181 (2017).[38] C. Hens, U. Harush, R. Cohen, and B. Barzel, Spatiotemporal signal propagation in complexnetworks, Nat. Phys. 15, 403 (2019).[39] A. Hacohen, R. Cohen, S. Efroni, B. Barzel, and I. Bachelet, Digitizable therapeutics fordecentralized mitigation of global pandemics, Scientific Reports 9, 14345 (2019).[40] S. Boccaletti, J. A. Almendral, S. Guan, I. Leyva, Z. Liu, I. Sendi˝na-Nadal, Z. Wang andY. Zou, Explosive transitions in complex networks? structure and dynamics: Percolationand synchronization, Physics Reports 660, 1-94 (2016).[41] I. Leyva et al., Explosive first-order transition to synchrony in networked chaotic oscillators,Phys. Rev. Lett. 108 (16), 168702 (2012).[42] X. Zhang, X. Hu, J. Kurths, and Z. Liu, Explosive synchronization in a general complexnetwork, Phys. Rev. E 88, 010802(R) (2013).[43] X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Explosive synchronization in adaptive andmultilayer networks, Phys. Rev. Lett. 114, 038701 (2015).14 usceptibleInfected (
𝐹 < 𝐹 𝑐 )Infected ( 𝐹 > 𝐹 𝑐 ) 𝑅 = 𝑅 σ PandemicMutatedLethargic (d)
Time 𝜌 ( 𝑡 ) tt ത 𝐹 ( 𝑡 ) Pandemic 𝑅 = 𝜎 = (a) Time 𝜌 ( 𝑡 ) tt ത 𝐹 ( 𝑡 ) Lethargic 𝐹 𝑐 𝑅 = 𝜎 = (b) Time 𝜌 ( 𝑡 ) tt ത 𝐹 ( 𝑡 ) 𝑅 = 𝜎 = 𝐹 𝑐 (c) Mutated 𝑃 Figure 1:
The phases of a pandemic under pathogen mutation . (a)
Pandemic . For R > ρ ( t ) vs. t (top) grows continuously as the fitness F ( t ) (bottom) increasesdue to mutation and natural selection. (b) Lethargic . For R < ρ ( t ) exponentially decaying to zero.The mutation rate σ = 0 .
01 is too slow F ( t ) remains almost constant (bottom), and the pathogen fails to reachcritical fitness F c (grey dashed line) on time. (c) Mutated . We now remain in the sub-pandemic regime R < σ = 1. For small t we observe ρ ( t ) rapidly decaying (top). However, thanks tothe rapid mutations F ( t ) reaches critical fitness (grey dashed line) within a short time. Following this point thedisease reemerges and ρ ( t ) changes course, turning pandemic. This is observed in the network example throughthe appearance of sporadic instances of high fitness pathogens (middle, dark red nodes), which then spread toinfect the majority of the population. (d) σ, R phase diagram . To systematically observe the different phaseswe varied R ∈ (0 , .
5) and σ ∈ (10 − , ,
050 epidemiological scenarios, with differentinitial µ, β and σ . For each scenario we ran 50 stochastic realizations and measured the probability P to have ρ ( t → ∞ ) > i.e . pandemic. We observe three phases with sharp boundaries between them. First, the pandemicphase (red) for R >
1, independent of σ , as predicted by the classic SIS model. In addition to that the sub-pandemic regime R < σ , P tends to zero (blue) and the pathogen fails tospread, giving rise to the lethargic phase. For large σ , the spreading probability becomes almost certain, as p ∼ P → P →
1, a dramatic shift occurring within a narrow range of R , σ values. This greyrange is well-approximated by our theoretical prediction (solid black line) as appears in Eq. (11). All simulations,here and throughout, were done on a random network of N = 5 ,
000 nodes and k = 15. The disease parameterswere set to µ = 0 . β = µR /k , to obtain the different values of R .The mutation rate σ is specified in each figure. In each scenario we set the initial condition to ρ ( t ) = 0 . = (d) 𝜌 𝑅 (e) σ = 𝑅 H i g h 𝑅 L o w 𝑅 High 𝑅 Low 𝑅 σ PandemicMutatedLethargic 𝑡 𝑡 𝑐 (c) Mutated 𝜌 ( 𝑡 ) 𝑡 𝑡 𝑐 (b) Critical 𝜌 ( 𝑡 ) 𝑡 𝑡 𝑐 (a) Lethargic 𝜌 ( 𝑡 ) 𝑅 = 𝜎 = 𝑅 = 𝜎 = 𝑅 = 𝜎 = Figure 2:
The transition to the mutated phase . For the mutated phase a critical mutation must arise beforethe pathogen is eliminated, namely before ρ ( t ) crosses 1 /N (grey dashed lines), capturing the unit line in whichthere is a single infected individual among the N node population. (a) ρ ( t ) vs. t (grey solid line) as obtained fromEq. (8) in the lethargic phase ( R = 0 . , σ = 0 . t c ), whichis below the unit line. Therefore the epidemic decays prior to the appearance of the critical mutation. Indeed,the stochastic simulation (blue solid line) approaches zero prevalence, never reaching the positive branch of ρ ( t ).(b) Setting σ = 0 .
16 the system is at criticality. ρ ( t c ) is adjacent to the unit line, and hence we observe criticalbehavior: some realizations decay (blue), whereas others successfully mutate. (c) Under σ = 0 .
5, the system is inthe mutated phase, ρ ( t c ) is sufficiently above the unit line and the critical mutation is reached with probability P →
1. (d) The lethargic-mutated phase boundary in Eq. (11) depends on the initial size of the infected population I . Here we show this boundary for I = 10 , . . . , . (e) The long term prevalence ρ = ρ ( t → ∞ ) vs. R under σ = 0 . R we begin with an initial small infection of I = 10 and observe a transition to the mutates phase at R = R High . In the opposite direction, however, as we begin withlarge R we approach the transition from an already pandemic state with I ∼ . Now the phase boundarytraverses through R = R Low . We, therefore arrive at a hysteresis phenomenon, in which the critical transitionpoint depends on the current state of the spread. Consequently, preemptive mitigation, done when the spread isstill at its embryonic stage ( I small), is more effective than reactive mitigation, applied when I is already large. 𝑅 = 𝑅 (a) Pandemic
MutatedLethargic 𝜌 𝑅 (b) σ = Time
Susceptible Infected (
𝐹 < 𝐹 𝑐 ) Infected ( 𝐹 > 𝐹 𝑐 ) t ത 𝐹 ( 𝑡 ) 𝐹 𝑐 𝜌 ( 𝑡 ) t 𝑅 = 𝜎 = (c) Volatile
Figure 3:
The volatile phase . (a) The σ, R phase diagram under the bounded fitness of Eq. (12). We nowobserve a volatile phase, in which ρ → σ is too large. Hence, the mutated phase (green) now onlyappears in the Goldilocks zone in which the mutation rate in not too high nor too low. The theoretical predictionof (13) is also shown (black solid line on right). (b) ρ vs. R under σ = 3 (yellow path in panel (a)). As opposedto the lethargic-mutated phase change, in the shift from volatile to mutated we observe a continuous secondorder transition. (c) ρ ( t ) vs. t in the volatile phase decaying, as predicted, to the healthy state ρ = 0. (d) F ( t )changes rapidly thanks to the large σ , and crosses the critical F c (grey dashed line) early on. However the rapidmutations prevent the slower processes of the natural selection from securing a steady increase in F ( t ). Hence,the achieved fitness cannot be stably sustained for the pathogen to continually spread. (e) Indeed, we observemultiple instances of critical fitness (dark red) that fail to reproduce and dominate the pathogen population. ime t 𝜌 ( 𝑡 ) t t R e s i s t a n t (a) (b)(c) Time t 𝜌 ( 𝑡 ) t t R e s i s t a n t (d) (e)(f) Time t 𝜌 ( 𝑡 ) t t R e s i s t a n t (j) (k)(l) Time t 𝜌 ( 𝑡 ) t t R e s i s t a n t (g) (h)(i) ResistantCritical
Slow vaccination 𝑣 = × -4 Fast vaccination 𝑣 = × -4 N o m i t i ga t i on 𝑅 = . M i t i ga t i on 𝑅 = . Figure 4:
Vaccination under the threat of mutation . We consider a pandemic pathogen ( R = 1 . ν = 2 × − ) and fast ( ν = 2 × − ) vaccination roll-out. (a) The prevalence ρ ( t ) vs. t asobtained from 50 stochastic realizations. In 72% of the cases vaccination fails to overcome the spread. (b Thereason is that all these cases a resistant mutant strain has overtaken the pathogen population. (c) As vaccinesare disseminated the selective pressure for resistance increases, until a single resistant strain (black node) spreadsto the entire network. (d)-(f) Rapid vaccination is insufficient, still witnessing 44% failure. (g)-(i) Combiningvaccines with mitigation ( R = 0 .
8) suppresses the spread, and delays the appearance of the resistant mutation.Still, due to the slow roll-out, the pathogen is afforded time to evolve higher fitness, reaching a critical mutation(red nodes), then gaining sufficient prevalence to also develop the rare vaccine resistance. Indeed, we find in panel(i) that first a critical mutation occurs, then only once it spreads, we observe instances of resistance mutations.(k)-(l) Therefore, the only strategy that ensures ∼ R below unity. Together allowing us to eradicate the epidemic before a strain with both R >1 (critical) andvaccine resistance ever evolves.