Svoboda | Graniru | BBC Russia | Golosameriki | Facebook
Journal of The Royal Society Interface
Open AccessResearch articles

Immune amnesia induced by measles and its effects on concurrent epidemics


    It has been recently discovered that the measles virus can damage pre-existing immunological memory, destroying B lymphocytes and reducing the diversity of non-specific B cells of the infected host. In particular, this implies that previously acquired immunization from vaccination or direct exposition to other pathogens could be partially erased in a phenomenon named ‘immune amnesia’, whose effects can become particularly worrisome given the actual rise of anti-vaccination movements. Here, we present the first attempt to incorporate immune amnesia into standard models of epidemic spreading by proposing a simple model for the spreading of two concurrent pathogens causing measles and another generic disease. Different analyses confirm that immune amnesia can have important consequences for epidemic spreading, significantly altering the vaccination coverage required to reach herd immunity. We also uncover the existence of novel propagating and endemic phases induced by immune amnesia. Finally, we discuss the meaning and consequences of our results and their relation with, e.g. immunization strategies, together with the possibility that explosive types of transitions may emerge, making immune-amnesia effects particularly dramatic. This work opens the door to further developments and analyses of immune-amnesia effects, contributing also to the theory of interacting epidemics on complex networks.

    1. Introduction

    The measles virus is among the most contagious human pathogens; it can cause severe symptoms and death, mostly during childhood and, as such, it represents a serious problem for global public health, targeted by the World Health Organization (WHO) [1,2]. In spite of the 73% global drop in measles deaths achieved thanks to improved vaccination policies in the period 2000–2018, measles is still common in many developing countries; indeed, over 500 000 cases were reported worldwide in 2019, more than a half of which in Africa [3]. Also in the USA as well as in Europe, where measles is considered endemic in at least 10 countries, outbreaks are becoming ubiquitous in recent years—with a 30% overall increase from 2017 to 2018—mostly as a consequence of anti-vaccination movements [4,5]. Moreover, the WHO has recently raised the alarm over the increasing chance of measles outbreaks due to poor vaccination coverage as the COVID-19 pandemic progresses, with millions of children at risk of missing out on measles vaccines [6].

    Given the magnitude of the problem, it should come as no surprise that, as of 2016, there were already over 100 mathematical models proposed in the literature to specifically reproduce and predict the evolution of measles outbreaks [7]. In spite of this wealth of modelling approaches, there is a crucial aspect of measles that is still systematically neglected and that has the potentiality to be more harmful than the outbreaks themselves: immune amnesia (IA).

    Building over a previous series of works that linked childhood mortality and severe immunosuppression with preceding measles-virus infection [810], conclusive empirical evidence has been very recently found that measles can wipe out acquired immunity to other infectious diseases through a mechanism called IA [1113]. More specifically, measles infection has been shown to destroy B lymphocytes (specific to whichever other pathogens) and to reduce the diversity of non-specific B cells, thus limiting severely the acquired defenses in the adaptive immune system (regardless of whether these have been achieved by means of vaccination or direct contact with a pathogen) [4,1113]. In fact, previous studies with rhesus macaques [8], as well as with unvaccinated children [11,13], had measured a depletion of up to 70% of the existing antibody repertoire across individuals after measles infection, even if there is a large subject-to-subject variability. As a matter of fact, it has been long reported that the majority of measles-related deaths are not due to the measles virus itself, but to secondary infections caused by the associated immunosuppression [14,15], hence the importance of taking into account IA into the broader field of epidemic mathematical modelling [16].

    In this work, we give a first step towards bridging this gap by incorporating the possibility of measles-induced IA into standard models of epidemic spreading for an arbitrary infectious disease co-occurring with measles outbreaks. Starting from an initial situation where vaccination coverage is assumed to grant herd immunity for a certain infectious disease—i.e. a sufficient number of vaccinated individuals so that the disease can hardly spread across the population—could measles outbreaks wipe out such immunity to the point where sizeable secondary epidemics are unleashed? If that was the case, the aforementioned recent increase in measles outbreaks worldwide could be a greater threat than previously thought. Even worse, a potential herd-immunity strategy relying on vaccination for COVID-19 could be hindered by the effects of measles outbreaks, all the more in countries where measles vaccination coverage during this health crisis is at its lowest.

    To analyse these issues, here we develop a relatively simple mathematical model that sheds light on the effects of IA over the dynamics of a second epidemic disease coexisting with outbreaks of measles. In particular, we perform mathematical analyses and extensive computer simulations of a modified susceptible–infectious–recovered (SIR) model that accounts for two coupled diseases, vaccination coverage, and possible demographic effects, as well as, crucially, the possibility of IA.

    We start considering homogeneously mixed populations, i.e. fully connected networks, and perform standard mean-field calculations that allow us to derive, e.g. analytical estimates for the minimum measles-vaccination coverage needed to maintain herd immunity for the second epidemics. Then, we extensively analyse, both theoretically and computationally, the impact that the structure of the underlying network of contacts can have on the results. In all cases, we elucidate the possible emergence of IA-induced phases, where the X disease becomes propagating/endemic just as a consequence of IA effects.

    Before closing, let us emphasize that for the sake of mathematical tractability our model considers a number of simplifying assumptions. Among others, let us underline that: (i) it assumes that infection with the measles virus erases the totality of the pre-existing memory cells in all cases, while this percentage has been shown to vary across individuals [10,11]; (ii) it assumes that there is no spontaneous waning of immunity neither for measles nor for the secondary infection, even though waning immunity is a well-documented fact [17,18] that has already been analysed in simple mathematical modelling approaches [1922]; (iii) it does not consider any explicit age-structure or population heterogeneity. We remark that all these ingredients can be straightforwardly implemented as additional layers of complexity in our model, and in particular we analyse the effects of considering an ‘imperfect’ IA in one of the appendices.

    2. The SIR-IA model

    In order to mimic the effects of IA on a given population, we extend the SIR model—either with or without demographic dynamics—[2327] to account for two co-occurring diseases: measles (M) and a second generic infectious disease to which we refer by X hereon. For the sake of illustration, we consider X to be COVID-19 as a guiding example, and use its associated epidemic parameters. We name this SIR-like model with IA, SIR-IA model.

    As in the standard SIR dynamics [2326], in the SIR-IA model each of the N individuals within the focus population can be either susceptible to be infected (S), infected (I), or resistant (R), for each of the two diseases. Thus, there are a total of nine possible states: i, j ∈ {SS, SI, SR, IS, II, IR, RS, RI, RR}, denoting the state of individuals that are simultaneously in state i ∈ {S, I, R} for measles and j ∈ {S, I, R} for disease X. It is important to remark that:

    State II representing individuals infected simultaneously of both diseases will be dismissed in first approximation as highly unlikely, given the short recovery periods.

    Resistant populations include not only recovered individuals but also those who achieved immunity through vaccination.

    The model dynamics is defined by a master equation including the set of possible transition rates between these states. This can be numerically integrated in an exact way by using the Gillespie algorithm [28] (see below). The set of possible transitions between states, together with the corresponding rates at which they occur, are schematically depicted in figure 1 (see also table 1 for a definition of model parameters together with their base-line values). In particular, following the standard notation, βM, βX and γM, γX denote the infectivity and recovery rates for measles and disease X, respectively. Note that IA is explicitly implemented within the term +γMρIR in equation (3.1e), which drags X-recovered individuals back into the pool of X-susceptible ones, in the case they had been infected with measles. The parameters, vX and vM represent the vaccination coverage for each disease, i.e. the probability that a new individual added into the system is vaccinated against measles and disease X, respectively. For COVID-19, it has been estimated that around a 65% of the population should be resistant (either through vaccination or naturally acquired immunity) in order to reach herd immunity [32]. Hence, to study to what extent IA effects can affect a potentially achieved COVID-19 herd immunity, we consider a hypothetical large vaccination coverage value, vX = 0.9.

    Figure 1. Sketch of the transitions between the eight allowed system states. Green (orange) cells correspond to measles (X) disease states. For the sake of clarity, demographic processes are not included.

    Table 1. Epidemic parameters used to model measles (M) and COVID-19 (X) epidemics. Infectivity and recovery periods for each disease were taken from [2931] to match the reported R0 values (e.g. for COVID-19 we took R0 = 3.4; although recent studies suggest a shorter infectivity period for COVID-19, this does not significantly affect the forthcoming conclusions). COVID-19 vaccination coverage was set to a high value to ensure initial herd immunity in the homogeneous-mixing approach.

    symbolbase-line value
    measles infectivity rateβM2.1 d−1
    measles recovery rateγM1/8 d−1
    COVID-19 infectivity rateβX0.25 d−1
    COVID-19 recovery rateγX1/14 d−1
    immigration/emigration rateμ1/365 d−1
    COVID-19 vaccination coveragevX0.9

    As for ‘demographic’ parameters, the death and birth rates for all individuals—regardless of their possible disease state—have been set to a common value μ. These rates can also be interpreted—looking at the problem from a meta-population perspective—as describing emigration and immigration processes. In particular, this latter interpretation justifies the use of relatively large rate values (table 1). In individual-based stochastic simulations of the model, any removed individual is instantaneously replaced by a new-arrived one, thus keeping a fixed population size.

    For the sake of simplicity, we begin by studying the case of homogeneously mixed populations and then analyse more structured populations with a non-trivial underlying network of contacts. We study versions of the model with either no explicit demography (i.e. μ = 0) or explicit demographic effects μ ≠ 0. In the first case, much as in the standard SIR model, there cannot possibly be any non-trivial stationary endemic state, while in the second such states can possibly exist [16]. We investigate in parallel all these possible scenarios to illustrate the generality of the conclusions from complementary perspectives.

    3. Results

    3.1. Homogeneously mixed populations

    To gain insight into the model key features, we employ a standard mean-field approximation which, as usual, is exact in the limit of infinitely large homogeneously mixed populations. This, leads rather straightforwardly to the following set of eight differential equations (sometimes called ‘rate equations’) [23]:

    where ρij is the population fraction in state ij.

    3.1.1. SIR-IA model without demography

    Let us begin by analysing the case with no demography, i.e. with μ = 0. In this scenario, the only role of vM and vX is to determine the initial fraction of vaccinated population of M and X diseases, respectively. To simplify the forthcoming mathematical analyses, let us define iM = ρIS + ρIR and sM = ρSS + ρSR as the total fraction of infectious and susceptible individuals of measles, and their counterparts iX = ρSI + ρRI and sX = ρSS + ρRS for the disease X. Let us remark that states SI and IS are not counted as susceptible ones (neither for measles nor for X disease) since states II have been neglected, as previously said.

    To analyse the situation in which herd immunity for the disease X is potentially erased by the effect of IA, we assume that the population fraction vX vaccinated against X is sufficiently high to initially avoid spreading the disease (i.e. iX ≈ 0 prior to the onset of a measles outbreak). Note that, under this assumption, from equations (3.1a), (3.1c), (3.1e) and (3.1f) one readily derives the standard SIR-model mean-field equations for measles dynamics

    From equation (3.2), it is clear that a small seed (of size εN) of M-infectious individuals introduced in a population of initially susceptible and resistant individuals grows exponentially if diM/dt > 0, and decays to 0 if this is negative. Therefore, since sM(0) = (N (1 − vM) − ε)/N ≈ 1 − vM is the initial fraction of M-susceptible to measles, the epidemic threshold separating the above two regimes is specified by the condition
    where, as usual, the basic reproduction number R0M=βM/(μ+γM)—defined as the average number of secondary infections generated by a primary case in a completely susceptible population [16]—has been introduced for measles. Therefore, vM is the herd-immunity threshold and represents the minimum fraction of M-vaccinated population needed to prevent measles spreading. For the considered parameter values (table 1) it follows that R0M17 and vM0.95, numbers that emphasize the well-known high infective power of the measles virus.

    In what respects the X disease and in the absence of IA, one can easily derive from equations (3.1a), (3.1b), (3.1d) and (3.1g) an epidemic-threshold condition analogous to equation (3.4)

    which results into R0X3.5 and vX0.65 for the parameters in table 1. Let us underline that, in this IA-free case, both diseases are uncoupled and, hence, their respective thresholds are independent of each other.

    To start scrutinizing the full problem, including immune-amnesia—which turns out to be much more intricate from an analytical viewpoint—we start by performing computational analyses of the stochastic model (see Methods for technical details). In particular, we run simulations of the SIR-IA stochastic model by implementing a Gillespie algorithm as follows: beginning with an initial seed of εM M-infectious individuals, we let an outbreak of measles spread through the population (for which we set vM<vM); once it fades out, we add at some initial time T a second seed consisting of εXX-infectious individuals who will potentially spread disease X through the system. Figure 2 illustrates the resulting time courses of epidemics as obtained from stochastic simulations averaged over many realizations. In particular, for the case in which the vaccination coverage vX is only slightly above the herd-immunity threshold vX (i.e. vXvX) the figure clearly shows that—on average—much larger outbreaks occur under the influence of IA. It also reveals that, not surprisingly, the duration of the outbreaks is shorter when IA effects are considered, as the disease takes over the susceptible population in a much faster way. The inset of figure 2 illustrates how these results change quantitatively with the M vaccination coverage, vM: as naively expected, (i) the total outbreak size grows and (ii) the time elapsed between measles and X outbreak peaks becomes shorter as vM is reduced. Observe also that X epidemics can only break out (i.e. become supercritical) if vM is set below a minimum, critical vaccination thresholdvM(μ=0)<vM, which represents the minimum population fraction that needs to be vaccinated against measles in order to preserve herd immunity for disease X even in the presence of IA effects.

    Figure 2. Time evolution of the fraction of: (i) M-infectious individuals (green curve), (ii) X-infectious individuals in the absence of IA (purple), and (iii) X-infectious individuals under the presence of IA effects (orange). X-vaccination coverage was set slightly over its herd-immunity threshold vX=vX(1+ε), with ɛ = 0.001. Conversely, vM=0.8<vM was set to guarantee the spreading of measles outbreaks. At time T = 100 days, a seed ε=104N of X-infectious individuals was inserted into the system. In the inset, we plot the maximum size of the X outbreak (measured as the total fraction of infectious individuals integrated in time) and the time elapsed between the green and orange peaks (in years) against M-vaccination coverage. The blue-dashed line marks the theoretical approximation to the vaccination threshold (see equations (A 30) and (A 31) in appendix A). Simulations were performed within the mean-field, homogeneous-mixing approximations, using the Gillespie algorithm with a total population size N = 105. Shaded areas around the curves indicate one standard deviation, as determined from over 50 independent runs. Parameters were chosen according to table 1.

    Last but not least, let us stress that, in spite of the fact that one needs to deal with a set of eight differential equations, we have been able to find an analytical solution for vM(μ=0) as a function of other parameters. The detailed calculations, which are obtained as a limiting case of the more general problem of a general network architecture, can be found in appendix B.1. Figure 2 confirms that the analytically derived threshold (blue dotted line in the inset), is in excellent agreement with computational results, a conclusion that remains true for other choices of parameter values. A summary of the analytical and computational results is provided by figure 3 which explicitly illustrates the existence, for a broad range of parameter values, of a propagating phase which emerges as a mere consequence of IA effects.

    Figure 3. Analytically determined phases for the X disease in the SIR-IA model without demography. The white-dashed line depicting the vaccination threshold under IA effects was determined by solving equations (A 30) and (A 31) in appendix A. A lower bound for the IA-induced endemic phase is given by vX=vX, below which the disease is always endemic, independently of IA effects (black-dashed line). Note that a triple point appears for vX=vX and vM=vM (marked by the white arrow), at which the three phases collide.

    3.1.2. SIR-IA model with demography

    In general, introducing demographic dynamics such as birth/death and/or emigration/immigration processes into a SIR-like model, opens up the possibility for stable endemic states, with a non-zero fraction of infectious individuals, to appear. For this, it is necessary that such processes occur at a fast-enough pace so that a flux of new susceptible individuals is constantly generated to ‘feed’ the contagion process, otherwise the epidemics necessarily vanish [16].

    As a first test to analyse the demographic version of the model, with μ ≠ 0, we verified the existence of endemic states by numerically solving the mean-field equations (3.1a)–(3.1h) with μ > 0 (see table 1) and X-vaccination coverage vX>vX, for which the X-disease-free state is stable in the absence of IA effects. In particular, the green curve in figure 4a shows that a measles endemic state is found as soon as the fraction vM drops below a certain measles herd-immunity level vM0.95. This result has been also verified by means of direct Gillespie simulations of the full stochastic model as well as proved analytically, as shown below. Figure 4a also reveals that—even if it has been obtained for a large X-vaccination coverage value—an X-disease endemic state exists (orange curve) provided vM drops below a certain critical threshold value vM(μ0). Such an X-disease endemic state is purely induced by IA effects, i.e. it is a IA-induced endemic phase.1

    Figure 4. Analysis of the stationary states and epidemic thresholds in the SIR-IA with demography. (a) Stable stationary state for the X-infectious (orange line) and M-infectious (purple line) population fractions as a function of the measles vaccination coverage vM (vanishing values in the disease-free steady states are not shown). Blue and green dashed lines mark the analytically found critical points vM=0.72 and vM=0.94, respectively, according to equation (3.10). In the insets, we illustrate the time evolution (in days) of the X-infectious (orange) and M-infectious (green) population fractions for the three different regimes found in figure 6. Curves were obtained by direct integration of equations (3.1a)–(3.1h). (b) Minimum fraction of M-vaccinated individuals, vM, necessary to prevent an endemic state of disease X as a function of the reproduction number R0X. The black curve represents the theoretical approximation for vM according to equation (3.10), while blue markers correspond to the numerical solutions obtained by studying the fixed points of the system and their stability (see Methods). All parameters were chosen according to table 1.

    We also analysed the dependence of vM on R0X by solving for the stationary states of equations (3.1a)–(3.1h); the endemic states are represented by blue diamonds in figure 4b. Clearly, the larger the virulence of the X disease the larger the value of vM. Observe in particular, that as R0X, vM approaches that of the critical point for measles outbreak vMvM, implying that as soon as an outbreak of measles takes place, an epidemic of X may emerge taking over the newly generated pool of IA-induced susceptible individuals, in spite of the fact that the population was massively vaccinated for the X disease (vXvX). Similarly, figure 5 shows the fraction of infected individuals in the stationary state (colour coded) as a function of both, vM and vX. Results for two different X diseases are shown: a mildly infectious one with (a) R0 = 1.3 and (b) COVID-19, with an estimated R0 = 3.4. It can be observed that, not surprisingly, the region in the vaccination parameter space where X-disease endemic states appear becomes larger with increasing R0X, but in both cases—even in the limit of stringent vaccination policies for the X disease, vX ≈ 1—an IA-induced endemic state can appear if M-vaccination drops below a critical threshold level vM. The exact value of this threshold (which marks the boundary of the stable endemic region) will depend on vX. Note also that the shift between endemic and disease-free states is a gradual one, as it corresponds to a continuous phase transition in the present case of homogeneously mixed populations.

    Figure 5. Fraction of infectious individuals for disease X (colour coded) in the stationary state, against vaccination coverage for both diseases. White-dashed lines indicate the theoretical approximation to the critical point for each vX value, as given by equation (3.10); population size N = 105. (a) Results for a mildly infectious disease, with R0X=1.3 slightly over SIR critical threshold. (b) Results for a more virulent disease, with parameters akin to those of COVID-19, as recorded in table 1, i.e. R0X3.4 Each point in the grid was obtained numerically by solving for the fixed points of equations (3.1a)–(3.1h).

    In summary, the value of the transition point vM(μ0) has been numerically shown to depend on both the infectivity R0X and the vaccination coverage vX for the disease X (as well as on other parameters such as μ). In what follows, in order to get a deeper understanding of the phenomenon as well as a better quantitative description of the phase diagram we derive analytical expressions for such dependencies.

    3.1.3. SIR-IA model with demography: theoretical results

    Let us first analyse the possible stable fixed points of equations (3.1a)–(3.1h) as a function of vM. Starting from equations (3.1a), (3.1c), (3.1e) and (3.1f) one can write

    Assuming, as above, initial herd immunity for the X disease, i.e. vX>vX, one can approximate iX ≈ 0 and readily find the two steady-state solutions of equations (3.6) and (3.7): a disease-free one, (iM*, sM*) = (0, 1 − vM) and an endemic one (iM,sM¯)=(R0M((1vM)1)μ/βM,1/R0M). A standard linear stability analysis allows us to recover the existence of a transcritical bifurcation at vM=11/R0M where these two fixed points exchange their stability; thus the system shifts in a continuous or smooth way from a non-propagating to an endemic state. Observe that this transition point is a natural extension of the threshold for μ = 0, even when no endemic state existed in that case, but just a separation between propagating and quiescent phases. Similarly, for the X disease, using the joint variables iX and sX

    Note that it is not possible to solve exactly the above equations for their fixed points, as they do not form a closed set: equation (3.8) depends on the fraction ρSS of individuals susceptible to both, measles and X disease, as well as on the fraction iM of measles infectious population. However, assuming that a steady state for measles has been reached, it is possible to show that ρSS=(1vX)/R0M; using this and searching for steady-state solutions of equations (3.8) and (3.9), one can compute the minimum value of vM preventing the existence of a stable X-endemic state, thus obtaining an expression for vM (see appendix A for a detailed derivation).

    This can be further simplified in the typical case under consideration where the immigration rate is much smaller than the measles recovery rate, μγM, resulting in:

    which depends not only on the R0’s of both diseases but also on the vaccination coverage for disease X, as computationally observed above. Note also that the rightmost expression in equation (3.10)—written in terms of vX rather that R0X—underlines the relationship between the two vaccination thresholds. This result is illustrated in figure 6, which shows the three resulting phases for a generic X disease: (i) disease-free, (ii) IA-induced endemic, and (iii) endemic, as well as their phase boundaries in the (vM,vX/vX) plane.

    Figure 6. Analytically determined phases for disease X in the SIR-IA model with demography. The white-dashed line depicting the vaccination threshold under IA effects was determined through equation (3.10). The IA-induced endemic phase is bounded by vX=vX, below which the disease is always endemic, independently of IA effects (black-dashed line) and vX=vXR0M (marked by the pink arrow), above which the epidemics does not propagate despite the IA contribution. Note that a triple point appears for vX=vX and vM=vM (marked by the white arrow), at which the three phases collide.

    Observe that, since we assumed vXvX, then trivially vM11/R0M=vM must hold, with both critical points coinciding for vX=vX. Similarly, as 0vM1, one readily sees that the existence of a non-trivial second threshold vM is limited by the constraint vX/vXR0M (these two limits are marked with arrows in figure 6). Therefore, an immune-amnesia-induced endemic phase for X disease exists under broad conditions for realistic parameter values.

    All the above analytical results, derived from linear stability analyses with some additional approximations, have been confirmed by numerically determining the fixed points for the full system of mean-field equations (3.1a)–(3.1h), without invoking any approximation beyond numerical accuracy (see Material and methods). Moreover, one can also cross-check the consistency between these analytical approaches and the previous results from computational simulations, for example, by looking at figures 4a and 5, which reveal that the analytical predictions for vM, as given by equation (3.10), explain well both the onset of the X outbreak in the deterministic calculation (figure 4a) and its dependence on R0X and vX (figure 5).

    3.2. SIR-IA model on structured networks

    After two decades of frantic activity on the development of the theory of complex networks, by now it is broadly recognized that the structure of the underlying network of contacts plays a crucial role in spreading phenomena such as epidemics [26,3339]. Thus, to have a broader view on IA effects, here we scrutinize the behaviour of the SIR-IA model beyond the homogeneous-mixing approach, considering more structured topologies such as Erdős–Rényi (ER) random networks and power-law degree-distributed networks [34,4042] (see Methods).

    In order to make further progress, we make some simplifying assumptions: (i) vaccination for both measles and X is considered to be performed in a random way across the network (i.e. there is no ‘targeted-immunization’ programme selecting preferentially specific nodes for vaccination according to, e.g. their network centrality or connectivity [4346]). (ii) As in the previous analyses, we impose that vx>vx to analyse how herd immunity can be potentially lost by IA effects. (iii) For simplicity, we limit ourselves here to the analytically more-tractable non-demographic version of the SIR-IA model (i.e. μ = 0).

    We first report on computational findings for stochastic (Gillespie) simulations of the SIR-IA model on structured networks. Let us remark that, in order to compare ER networks with different average connectivity (or ‘degree’) 〈k〉, we defined an infectivity per contactβ0 so that β = β0k〉, with β0X=0.017 and β0M=0.14. Within this convention, the values presented in table 1 are recovered for an average network degree 〈k〉 ≈ 20. These same values of β0X and β0M were also used in power-law degree distributed networks, but in this case the average connectivity 〈k〉 = u (1 − α)/(2 − α) was kept fixed by changing the minimum node degree u according to the chosen exponent α.

    Figure 7a shows results of the mean epidemic size for X-disease outbreaks occurring on ER networks with different average connectivity, as a function of the measles vaccination coverage vM. It can be readily noticed that the vaccination threshold vM grows and the transition becomes sharper as 〈k〉 is increased. Remarkably, the transition becomes very abrupt for large mean degrees, implying that a small variation in the measles vaccination level can induce a dramatic effect on the typical size of the subsequent X-disease outbreaks. For instance, for 〈k〉 ≥ 30, lowering the vaccination coverage from vM = 0.9 to vM = 0.82 (a reduction on the number of vaccinated individuals of just an 8% of the population size) entails an increase of two orders of magnitude in the subsequent X-outbreaks.

    Figure 7. Analysis of epidemic size and thresholds in ER and power-law degree distributed networks. Outbreak sizes (as measured by the difference between X-resistant individuals before and after the introduction of one single X-infectious individual in the system) are presented for ER networks with different average degrees (a) and power-law networks with different degree exponents α (c). The dependence of the vaccination threshold, vM, with the average network degree for ER networks (b) and α exponent for power-law networks (d), respectively, is compared with the corresponding analytical predictions given by the heterogeneous mean-field approximation. Stochastic simulations were performed in networks of size N = 105 and N = 5 × 105 for ER and power-law networks, respectively. Error bars are computed as the standard deviation over 200 realizations. The rest of parameters were chosen according to table 1.

    On the other hand, for power-law degree distributions pkkα, numerical simulations reveal that reducing α leads to larger values of vM, as illustrated in figure 7c. The figure also shows that, in this case, the transition becomes more abrupt as larger values of α are considered. In light of these results, one could wonder whether the transition could eventually become discontinuous or explosive for larger values of 〈k〉 or α, as it has been shown to occur in other models with cooperative contagion [4749].

    Remarkably, the dependence of the vaccination threshold on the average network degree and exponent α can be described using the heterogeneous mean-field approach, a generalization of the mean-field theory that groups all nodes with the same degree into a common class [39], which has been shown to be accurate (up to finite size corrections) when applied to epidemic models that do not posses non-trivial steady states [50]. For the SIR-IA model, the calculations turn out to be a mathematical ‘tour the force’, so the details are deferred to appendix B, and here we just present the final results for both ER and power-law degree-distributed networks.

    In particular, for ER networks, calling q=s1M(T)/(1vM) to the fraction of degree-1 M-susceptible individuals at the time T of inserting an X-infectious seed with respect to the initially M-susceptible population fraction, we found the following equation:

    where RMeff=β0M/γM. Defining the auxiliary function ζ(q) = q vX(1 + 〈kq) e−〈k〉(1−q) we were able to derive a closed equation for the vaccination threshold as a function of parameter values
    This expression, when inserted into equation (3.11), can be numerically solved, leading to an analytical determination of vM. Figure 7b shows such a theoretical solution (black-dashed line) as a function of the network average degree, together with the computational results obtained by averaging over many runs of the stochastic simulations of the model running on ER networks with different averaged connectivities.

    On the other hand, for power-law degree-distributed networks with α > 3 we defined ξ=log((1vM)/suM)—where u is now the minimum degree of a node in the network—and vX=1(k/RXeffk2), which is the corresponding threshold in a standard SIR model with vaccination [46]. In terms of these quantities, we obtained (see appendix B.3),

    with the constraint
    needed to determine ξ. Conversely, for scale-free networks (i.e. 2 < α ≤ 3) one finds that outbreaks do not propagate only if vM = vX = 1, i.e. the whole population needs to be vaccinated for both diseases to prevent the epidemic. This result is in consonance with the well-known phenomenon of vanishing epidemic threshold in BA networks [48,5153]: all single nodes need to be vaccinated to prevent epidemic propagation, a result that stems from the existence of super-spreaders, i.e. network hubs. However, our analytical results predict also that, for the SIR-IA model, the M-vaccination threshold for the propagation of an X-disease outbreak depends also on the X-vaccination coverage and, hence, the vaccination threshold can saturate even in scale-rich networks, with α > 3 (see figure 7d, where for vX = 0.9 a minimum M-vaccination level preventing an X outbreak is defined only in networks with α3.4). Nevertheless, it is important to remark that these result are strictly valid only in the infinite size limit (N → ∞), in which the second-moment of the degree distribution truly diverges in scale-free networks. In fact, for 〈k〉 = 15 and α = 3, we found a vaccination threshold clearly below 1 in stochastic simulations (figure 7d). This phenomenon that should come as no surprise, since finite size effects in the SIR model have been shown to be responsible for the appearance of non-trivial thresholds in scale-free networks at sufficiently low transmission rates [52].

    As shown in figure 7b,d, the analytical results match quite well the values of vM found in stochastic simulations. In both ER and power-law networks, small disagreements with computational results are most likely rooted in finite-size effects and the associated possibility of stochastic fade-out that occurs in finite networks and not in infinitely large ones; however, performing a proper finite-size analysis to study such effects is beyond the scope of the present work [5456].

    Finally, to the question of whether discontinuous transitions could be found in power-law networks with sufficiently high values of α, application of the heterogeneous-mean-field theory to our model predicts that—provided a non-trivial vaccination threshold exist in the large-N limit—transitions between a quiescent and a propagating phase are always continuous independently of the α exponent (see appendix B.3) and a similar conclusion holds for ER networks even in the limit of large average connectivities. Still, it remains to be carefully investigated how other structural aspects such as clustering, geography and network modularity, to name but a few, could affect such a conclusion.

    4. Conclusion and discussion

    We have seen that, when measles vaccination policies are relaxed, the expected herd immunity for any secondary infectious disease X can be lost owing to the proliferation of individuals affected by IA. In particular, under IA effects, the epidemic threshold is shifted so that severe outbreaks can take place even under extensive X vaccination. We have studied the conditions under which measles vaccination can prevent such outbreaks. To support the generality of our findings, we considered two different variants of the SIR-IA model: one in which all the state transitions are given by infectious/recovery processes, and a second version in which demographic effects are also taken into account. For both models, we were able to derive—under homogeneous-mixing assumptions—analytical expressions for the epidemic threshold in terms of the fraction of people vaccinated against measles, which were able to reproduce with significant accuracy the results obtained through simulations. Remarkably, the analytical results also allowed us to construct a full phase diagram in both models, where three distinct phases were found: quiescent, endemic/propagating and, more remarkably, IA-induced endemic/propagating phases.

    We have also studied the persistence of the generic X disease under IA when more realistic network architectures—beyond the homogeneously mixed population paradigm, were considered. In particular, when the SIR-IA model is implemented in random ER networks, it was shown that larger fractions of M-vaccination coverage are necessary to prevent outbreaks as the average network connectivity increases. We remark that this dependence is highly nonlinear (as also shown by our analytical results), with the epidemic threshold decaying very fast at low connectivity values. Thus, when fighting an outbreak in ER networks under IA effects, it is likely that measures taken to lower the average connectivity, even without imposing full confinement—e.g. limiting the allowed number of individuals in gatherings—can indeed have a profound impact on spreading of the disease.

    At the other end of the spectrum, for scale-free networks the vaccination threshold is equal to unity, implying that the whole population needs to be vaccinated against both diseases to prevent epidemic propagation; this is the counterpart of the well-known phenomenon of vanishing epidemic threshold, and is predicted by our analytical calculations in the limit of very large networks. The fact that outbreaks persist even when a large percentage of the population is vaccinated manifests not only the key role of hubs (i.e. super-spreaders) in the spreading process but also the scarce effectiveness of vaccine uptake measures when these are randomly administrated (as opposed to, e.g. targeting the most connected nodes [44,46]). Aiming for a more profound understanding of this effect, we considered general power-law degree-distributed networks which drift from the scale-free to the random regime by considering values of α > 3. As shown by stochastic computer simulations, the vaccination threshold becomes non-trivial once the network presents a non-divergent second moment, scaling with the exponent α as predicted by the analytically derived expression in the limit of N → ∞.

    Let us also discuss the nature of the transition between the propagating and quiescent regimes. It was shown in [57] and further investigated in [47,48,58,59], that cooperative epidemics can show hybrid-type phase transitions for large enough α in the limit of large system sizes. Although performing a detailed analysis of this important issue is beyond our scope here (and is left for a further work), the implications of a possible discontinuous transition are enormous from the dynamical perspective, opening the door to catastrophic regime shifts [60,61]: under such conditions—with just a slight reduction in the fraction of the vaccinated population—the system could suddenly undergo a transition from a quiescent state with overall herd immunity to a state in which anomalously large pandemics could surge. As discussed for the homogeneous-mixing case, at the epidemic threshold we find a second-order phase transition, which holds at first sight when more realistic networks are considered. Nevertheless, we see that the transition becomes more abrupt as we consider larger mean degrees in ER networks, or greater exponents α in power-law degree distributed networks. For the latter, however, application of the heterogeneous mean-field theory to our model predicts a continuous phase transition independently of the value of α. More realistic network architectures including clustering, modularity or embedding into a metric space, as well as possibly temporally changing networks, still need to be analysed to have a full view of the potential effects of IA on real populations. Furthermore, the possibility of finding actual first-order transitions when further cooperative interactions are implemented is left for future work. In particular, let us note that during lockdown periods to fight COVID-19, vaccination coverage for measles decreases, thus creating a feedback loop that could enormously enhance the impact of a future M-outbreaks and thus of IA.

    Our results not only pave the way to the study of cooperative contagions from the perspective of immunization (in comparison with other works, in which cooperativity is modelled as changes in infectivity or recovery rates [57,59,62]), but also open a branch of new exciting questions. For example, it is known that vaccination strategies targeting the hubs in scale-free networks can effectively reduce to finite values the epidemic threshold [63]. It would be therefore very relevant to study the effects that using different immunization strategies for each disease (e.g. random and targeted vaccination) have on the existence or absence of an epidemic threshold in scale-free networks. Let us also emphasize that we have developed an advanced version of the heterogeneous mean-field approach which can be applied to derive analytical results in similar problems of cooperative contagion or, more in general, when different epidemics coexist, influencing each other.

    The present work can be extended in a number of ways to include further biological and epidemiological details. For instance, in appendix C, we consider a version of the model in which IA is only partial, as infection with the measles virus does not usually remove all existing memory cells, but just a fraction of them; as naively expected, and explicitly demonstrated in appendix C, this effect diminishes the impact of IA reducing the size of the IA-induced endemic phase. Similarly, one could also consider waning immunity for both measles and the secondary disease [1722]; such effects can be easily implemented in our model by allowing for spontaneous transitions from recovered to susceptible states. In this case, these effects would add up to IA, further enlarging the endemic phase and reducing the quiescent one. It remains a challenging goal for further studies to quantitatively analyse these diverse types of immunity loss when they occur in concomitance. Finally, additional aspects such as age-structured populations, spatial distributions, seasonality, etc., could be implemented in our model. Careful analyses of all these possibilities will certainly contribute to build a more clear picture of the quantitative epidemiological aspects of immune amnesia, but are beyond the scope of the present work.

    It is our hope that this work makes it clear the importance of keeping measles vaccination (and vaccination in general) at levels that are as high as possible, to prevent IA effects to have a strong negative impact at a global level. We also hope that this work fosters further investigations along these lines as well as novel developments in other directions taking advantage of the techniques we have set up.

    5. Material and methods

    The steady-state solutions reported in figures 4b and 5 were obtained by solving equations (3.1a)–(3.1h) for their fixed points. The eigenvalues of the associated Jacobian matrix were then analysed to determine their possible stability [64]. Only the resulting stable fixed points are plotted in such figures. Likewise, the deterministic trajectories in the standard mean-field approximation shown in figure 4b (insets) were obtained by integrating equations (3.1a)–(3.1h) with Matlab ode23 function, which implements an explicit Runge–Kutta (2,3) pair algorithm [65]. On the other hand, simulations of the stochastic dynamics were performed through the standard Gillespie algorithm in the homogeneous mean-field approach [28] and through a network-adapted Gillespie for any other network structure, as described in [66]. Unless otherwise specified, initial conditions for the simulations were set in agreement with the vaccination rates following ρ0 = {ρSS = (1 − vM) (1 − vX) − (ɛX + ɛM), ρSI = ɛX, ρSR = 0, ρIS = ɛM, ρIR = 0, ρSR = (1 − vM)vX, ρRS = vM(1 − vX), ρRR = vMvX}. To obtain the deterministic solutions (figures 4b and 5), we set ɛM = ɛX = 50/N, while for the stochastic simulations we chose ɛM = 1/N and ɛX = 0, introducing one X-infectious individual only after the measles outbreak has fade out. In all cases, the total population size was fixed to N = 105.

    In what respect structured networks, ER graphs were constructed following the standard algorithm as described in [40], while for power-law degree-distributed networks we used the configuration model [67] to generate a first graph from which we then removed all multiple and self-connections. Although this last step may introduce correlations within the networks, they are negligible for all purposes when α > 3 and large system sizes are considered [68].

    Data accessibility

    This article has no additional data.

    Authors' contributions

    Both authors conceived and designed the study, and they all drafted, read and approved the manuscript. G.B.M. carried out the computational simulations and mathematical analyses.

    Competing interests

    We declare we have no competing interests.


    We acknowledge the Spanish Ministry and Agencia Estatal de investigación (AEI) through grant no. FIS2017-84256-P (European Regional Development Fund), as well as the Consejeria de Conocimiento, Investigación Universidad, Junta de Andalucía and European Regional Development Fund, Ref. A-FQM-175-UGR18 and SOMM17/6105/UGR for financial support.


    We thank V. Buendía and M. Sireci, for very helpful discussions.


    1 Let us remark that vM depends on μ and, in particular, it does not coincide for the cases with and without demography.

    Appendix A. SIR-IA model with demography in homogeneous networks

    We begin the analytical study by assuming that, initially, the number of X-infectious individuals is very low due to a high-vaccination coverage vX1. Under this approximation, equations (3.2) and (3.3) of the main text describing the evolution of the total fraction of susceptible and infectious individuals of measles, can be written as

    diMdt=βMiMsMγMiMμiManddsMdt=βMiMsMμsM+μ(1vM).}A 1
    In this way, the dimensionality of the problem is largely reduced, allowing us to compute the quiescent (Q) and endemic (E) stationary states in the now independent set of variables {iM, sM}
    ξQ=(iMQ,sMQ)=(0,1vM)andξE=(iME,sME)=(μ(R0M(1vM)1)βM,1R0M).}A 2
    To analyse the stability of the above fixed points, we compute the Jacobian associated with the linearization of equations (3.2) and (3.3) around a point in which there is no infectious individuals of X (i.e. iX = 0), as we assumed to be above the herd-immunity threshold for this disease. Therefore, linearizing around any point (iM*, sM*, iX = 0), we have
    JM=(βMsM(γM+μ)βMiMβMsMβMiMμ).A 3
    If vM is chosen as the control parameter, one can perform a stability analysis by evaluating the determinant (Δ) and the trace (τ) of JM at the each of the fixed points [64]. Thus, the disease-free fixed point ξQ is a stable focus of the linearized, X-disease-free system if vM>11/R0M , and a saddle-node otherwise. On the other hand, the endemic stationary state ξE only exists provided vM<11/R0M, and it is always a stable focus. Therefore, the epidemic threshold for the propagation of measles is given by the critical value vM=11/R0M.

    Now, to study the epidemic threshold for X using a similar approach requires of some approximations. In particular, we assume that a stationary state for the measles variables sM and iM is reached (regardless if endemic or quiescent) before the evolution of disease X can significantly progress. This a reasonable assumption because, on one hand, the dynamics of measles occurs at a faster timescale (shorter recovery period and faster infection rate). Moreover, we start with an initial state for X that is over herd immunity, and thus activity (in terms of the fraction of X infectious individuals) is very low until IA effects can convert a large fraction of resistant individuals into susceptibles. Within this approximation, it is possible to solve equations (3.8) and (3.9) in the main text for the stationary states, obtaining

    ψQ=(iXQ,sXQ)=(0,iMμ(γMβMρSS)+(1vX))andψE=(iXE,sXE)=(iM(γMβMρSS)+μ(1vX(1/R0X))γX+μ,1R0X),}A 4
    where iM* and ρSS* denote the steady-state fraction of M-infectious and fully susceptible individuals, respectively.

    Beginning with a scenario where the system is in a disease-free stationary state ξQ for measles, we find

    ψQ,1=(iX,1Q,sX,1Q)=(0,1vX)andψE,1=(iX,1E,sX,1E)=(μβX(R0X(1vX)1),1R0X).}A 5
    Note that these stationary states are the X-counterparts of ξQ and ξE. This symmetry reflects the fact that, without a preceding outbreak of measles and its consequent IA effects, the evolution of the X infection should follow that of an independent SIR model with vaccination and demographic dynamics, where the endemic state ψE is stable provided vX<vX=11/R0X.

    Let us move now to the more interesting case in which there is a endemic measles stationary state, given by ξE. Since now equations (3.8) and (3.9) do not form a closed set, we first proceed to estimate ρSS using sM = ρSS + ρSR and equation (1c). Assuming stationarity in the number of M-infectious individuals and low initial fraction of X-infectious (iX ≈ 0), one can write

    dρSRdtβMρSRiM2μρSR+μvX(1vM)=0,A 6
    from where ρSR=vX/R0M and ρSS=(1vX)/R0M. It is important to remark that these are not truly stationary states, but serve us as an estimation of the expected distribution of M-susceptibles near the measles endemic state ξE before the onset of an X outbreak. Let us nevertheless carry on with the analysis and impose fixed values ρSS* and iM=iME in equation (A 4). One then obtains the following pair of stationary states:
    ψQ,2=(iX,2Q,sX,2Q)=(0,γMβM(R0M(1vM)1)+(1vX)(vM+1R0M))andψE,2=(iX,2E,sX,2E)=(μμ+γX(sX,2Q1R0X),1R0X).}A 7

    As before, one can linearize the above system around each of these solutions. Studying the determinant and trace of the resulting Jacobian matrix in each point

    JX=(μβXiXβXsXβXiXβXsX(γX+μ)),A 8
    allows us to conclude that the X-disease-free steady state ψQ,2 is a stable focus if sX,2Q<1/R0X, and a saddle-node otherwise; whereas ψE,2 is a stable focus only when sX,2Q>1/R0X and does not exist otherwise. Therefore, writing the explicit form found for sX,2Q and taking vM as the control parameter, one finally obtains the epidemic threshold for the X disease in terms of the fraction of M-vaccinated population
    vM=(1/R0X)(1vX)(1/R0M)+(γM/βM)(1R0M)1vX(R0MγM/βM).A 9
    In the typical case in which μγM, one can approximate γM/βM1/R0M, and the vaccination threshold presented in equation (3.10) of the main text is recovered
    vM1vX(11R0X)1R0M=vXvX1R0M.A 10

    A.1. SIR-IA model without demography in heterogeneous networks

    The heterogeneous mean-field (HMF) approach [39] is a degree-block approximation carried out by assuming that nodes with the same degree k are statistically equivalent, i.e. belong to the same connectivity class. Within this framework, the total fraction of X-infectious individuals at each degree-class, ikX, obeys

    dikX(t)dt=β0XkskXθXγXikX,A 11
    where θX=(1/k)kkpkikX is the density function, defined as the fraction of infected nodes in the neighbourhood of a susceptible node with degree k, and can be shown to be independent of k if the network has no degree-correlations [51]. Multiplying both sides of the above equation by the excess degree, kpk/〈k〉, and summing over all possible values of k, we get an equation for the evolution of the density function for X-infectious
    dθXdt=1kβ0Xk2sXθXγXθX.A 12
    Placing a single X-infectious individual at time T, the condition for an outbreak to spread in the mean-field approximation is trivially given by dθX(T)/dt > 0. Hence, the epidemic threshold is determined by the condition
    k2sX(T)=kRXeffA 13
    where RXeff=β0X/γX is an effective R0 per contact. To estimate the critical value 〈k2sX(T)〉 for which the epidemic breaks out, we need an expression for the minimum fraction of X-susceptible individuals in each degree-class at the time T of inserting the X-seed
    skX(T)=ρkSS(T)+ρkRS(T).A 14
    Using the same approximation as in the case with μ ≠ 0, we consider that an outbreak of measles propagates and dies out before a seed of X-infectious individuals is placed in the system at time T (i.e. ikM(t)=0,tT and ikX(t)=0,t<T). If this is the case, then ρRS is a ‘dead-end’ state for tT, and three sources contribute to its value at t = T:

    Initial fraction of M-resistants that are not resistant for X: ρkRS(0)=vM(1vX).

    Fully susceptible individuals that became resistant for measles after undergoing an infection (ρSSρRS): ρkSS(0)ρkSS(T)=(1vM)(1vX)ρkSS(T).

    X-resistants who became infected with measles, and lost their immunity due to IA effects (ρSRρRS): ρkSR(0)ρkSR(T)=(1vM)vXρkSR(T).

    Hence, summing up all contributions one obtains

    skX(T)=1vMvXρkSR(T).A 15
    The epidemic threshold for X depends therefore on the fraction ρkSR(T) of M-susceptible individuals who are X-resistant at the end of the measles outbreak. During the latter, ikX=0k, and thus we can write
    dρkSS(t)dt=β0MkρkSSθManddρkSR(t)dt=β0MkρkSRθM.}A 16
    Dividing the above equations and integrating
    ρkSR(T)=(ρkSR(0)ρkSS(0))ρkSS(T).A 17
    Now, since vaccination is randomly performed with the same probability across all nodes, we have ρkSR(0)=(1vM)vX and ρkSS(0)=(1vM)(1vX)k, whence
    ρkSR(T)=vX1vXρkSS(T)=vX1vX(skM(T)ρkSR(T))=vXskM(T).A 18
    This relation between ρkSR(T) and skM(T) allows us to finally write equation (A 15) as an expression linking the fraction of X and M-susceptibles at time T
    skX(T)=1vMvXvXskM(T).A 19
    It follows then that the value vM at which the epidemic threshold takes place, according to equation (A 13), is given by
    k2(1vMvX)vXk2sM(T)=kRXeff.A 20
    All that remains now is to determine the fraction of M-susceptibles at the end of the measles outbreak. To do so, we follow a similar approach to the one presented in [69] for a standard SIR model. Let us define the minimum degree of any node in a network by u = min(k) and write the equation for the rate of change of M-susceptible individuals in the class of degree k
    dskMdt=β0MkskMθM.A 21
    Noticing that dskM/dsuM=kskM/usuM, one can then write
    skM(t)=C0(suM(t))k/u.A 22
    with C0(k)=skM(0)/(suM(0))k/u=(1vM)1(k/u). Therefore, the fraction of susceptible individuals for each degree class at the end of the measles outbreak can be expressed as
    skM(T)=(1vM)(suM(T)1vM)k.A 23
    Note how, with the above relation, it is possible to determine the fraction of M-susceptible individuals at any time and degree-class k just by studying the minimum-degree class. We introduce now the variable z(t)=log(suM(t)), so that
    z˙=s˙uM(t)suM(t)=βM0uθM.A 24
    In terms of this new quantity and using equation (A 23), one can write the density function for M-infectious nodes as
    θ˙M=1kkkpkdikMdt=βM0kkpkkC0(k)(suM)k/uθMγMθM=θM(βM0kkpkk2C0(k)e(kz/u)γM).A 25
    Dividing by equation (A 24)
    dθdz=kpkkC0(k)e(kz/u)ukγMβM0u,A 26
    and noticing that iM(0) ∼ 0 implies θM(0) ∼ 0, with z(0) = −log (1 − vM), one can integrate the above equation to obtain
    θM=(1vM)(1kpkkk((1vM)ez)(k/u))γMβM0u(z+log(1vM)).A 27
    At the end of the measles outbreak θM(z(T)) = 0, and one finally obtains an equation for suM(T) that will obviously depend on the degree distribution pk of the network under consideration
    0=(1vM)(1kpkkk(1vMsuM(T))(k/u))γMβM0ulog(1vMsuM(T)).A 28
    Armed with equations (A 20), (A 23) and (A 28), we are now in a position to confidently explore the effects of different networks topologies on the epidemic threshold.

    A.2. Homogeneous networks

    In a homogeneous network, all nodes have the same degree 〈k〉 (i.e. pk = δ(k − 〈k〉)), so all connectivity classes are equal. Thus, we have skM=sM, with u = 〈k〉, and equation (A 28) is reduced to

    0=(1vM)sM(T)1R0Mlog(1vMsM(T)).A 29
    Solving for sM(T), we find the following self-consistency equation
    sM(T)=(1vM)eR0M((1vM)sM(T)).A 30
    Similarly, if all nodes have the same degree 〈k〉, equation (A 20) for the epidemic threshold can be easily reduced to
    vM=vXvXsM(T).A 31
    Note that, in the case in which μ ≠ 0, the stationary state for the fraction of M-susceptibles is given by sM(T)=1/R0X, and we recover the epidemic threshold found in equation (A 10). Equations (A 30) and (A 31) allow us to reproduce an equivalent—and in fact, very similar—phase diagram to the one found for the demographic model (see figure 3 in the main text).

    A.3. Erdős–Rényi

    For an ER network and if 〈k〉 ≪ N, the binomial distribution typical of random networks can be approximated by a Poisson distribution, pk = (〈kk/k!) e−〈k. Moreover, in the limit of very large networks, one can assume u = 1 to be the minimum degree in the network. In terms of this distribution, we can write

    kpkkk((1vM)ez)k=ekkkk1(k1)!((1vM)ez)k=e(z+k)1vMexp(kez1vM).A 32
    Therefore, equation (A 28) can be expressed as
    0=s1M(T)ek(s1M(T)/(1vM)1)1RMefflog(s1M(T)1vM)(1vM),A 33
    which has a trivial solution s1M(T)=skM(T)=1vM, representing the case in which no measles-epidemic has taken place. On the other hand, making use of equation (A 23) we have
    where 1F0(a, b; z) is the generalized hypergeometric function [70]. Thus, using that 〈k2〉 = 〈k〉(〈k〉 + 1) for ER networks, the epidemic threshold condition is written as
    1RXeff=(k+1)(1vMvX)vX[s1M(T)(ks1M(T)1vM+1)ek((s1M(T)/(1vM))1)],A 34
    with s1M(T) given by
    s1M(T)ek(s1M(T)/(1vM)1)1RMefflog(s1M(T)1vM)(1vM)=0.A 35
    In order to simplify these expressions further, one can observe that s1M(T) must necessarily be a fraction q of the initial number of susceptibles (i.e. s1M(T)=qs1M(0)=q(1vM)) and then rewrite the above equations as
    1RXeff=(k+1)(1vMvX)vXek[q(1vM)(kq+1)ekq]A 36
    q(1vM)ek(q1)1RMefflog(q)(1vM)=0.A 37
    Regrouping terms in the epidemic threshold condition and defining ζ(q) = qvX(1 + 〈kq) e−〈k〉(1−q), we finally obtain
    vM(q)=1/RXeff(1+k)+ζ(q)ζ(q)vX(1+k).A 38
    which can be inserted in equation (A 37) to obtain a mathematical expression depending solely on q. We can then solve this equation to finally study the dependence of vM with the network average degree, as done in the upper-right panel of figure 6.

    A.4. Power-law degree-distributed networks

    These networks are characterized by a probability distribution of the form pk = C0kα, where C0 is determined through the normalization condition k=ukmaxpk=1. Furthermore, for networks with a large number of nodes, one can resort to the continuum formalism and replace the sums by integrals over k. Working in the limit of N → ∞, we then have pk=(α1)uα1kαΘ(ku). In terms of this probability distribution function, one can write

    kpkkk((1vM)ez)ku(α1)ukξα2ξexx(2α)1dx=(α2)ξα2Γ(2α,ξ),A 39
    where we have used x = (k/u)ξ, with ξ=log((1vM)/suM); 〈k〉 = ((1 − α)/(2 − α))u is the average degree of the network, and Γ(a,z) is the upper incomplete gamma function [70]. Inserting the above result into equation (A 28) one obtains
    ξuRMeff=(1vM)(1(α2)ξα2Γ(2α,ξ)).A 40
    On the other hand, equation (A 23) can be used again to perform an equivalent analysis of 〈k2sM(T)〉 to the one for ER networks
    k2sM(T)=(1vM)(α1)uα1k=uk2α(suM(T)1vM)k/u(1vM)(α1)uα1uk2αeξ(k/u)dk=(1vM)(α1)u2ξα3ξexx(3α)1dx=u2(1vM)(α1)ξα3Γ(3α,ξ).A 41
    Finally, the vaccination threshold condition reduces to
    kRXeff=k2(1vMvX)vXu2(1vM)(α1)ξα3Γ(3α,ξ),A 42
    which, using k2=(1α3α)u2 and noticing that vX=1(k/RXeffk2) is just the immunization threshold in the standard SIR model (derived in the thermodynamic limit) for power-law degree distributed networks [46], can be written as
    vM=(α3)ξα3Γ(3α,ξ)vX/vX(α3)ξα3Γ(3α,ξ)1.A 43
    with a condition over suM(T) given in terms of ξ by equation (A 40). Although one could in principle solve the system given by equations (A 40) and (A 42), we show in what follows that it is possible to derive the epidemic threshold from a slightly different approach that we believe is more informative. Beginning with the equations for the total fraction of susceptibles and resistants for X, and assuming as usual that measles outbreaks precede X outbreaks, we can write
    dsXkdt=βX0ksXkθXanddrXkdt=γXiXk.}A 44
    Solving each equation between the time T at which the X seed is introduced and t
    sXk(t)=sXk(T)eβX0kTtθX(t)dtandrXk(t)=rXk(T)+γXTtiXk(t)dt.}A 45
    Defining ϕX=TtθX(t)dt one can then write
    ϕX=kpkkkTtiXk(t)dt=1kkpkk(rXk(t)rXk(T)γX).A 46
    Using iXk(t)=1sXk(t)rXk(t) and sXk(T)=1rXk(T)
    One can now make use of equation (A 19) and 〈k2〉 = ((1 − α)/(3 − α))u2, replacing again the infinite sums by integrals over k to obtain
    ϕ˙X=γXϕX+1vMvX+vX(1vM)(α2)ξα2Γ(2α,ξ)(1vMvX)(α2)xα2Γ(2α,x)+vX(1vM)(α2)(x+ξ)α2Γ(2α,x+ξ).A 47
    where for ease of notation we have defined x=βX0uϕX. On the other hand, provided s ≠ 0, − 1, − 2, …; it is possible to expand the upper incomplete gamma function as [71]
    xα2Γ(s,x)=xα2Γ(s)n=0(1)nxnn!(s+n),A 48
    where Γ(a) can be naturally extended to negative arguments a < 0 by means of the recursive relation Γ(a+1)=aΓ(a). Applying the former identity, and after several manipulations, one can show that
    (x+y)α2Γ(s,x+y)=xα2Γ(s,x)+yα2Γ(s,y)+1s+Γ(s)[(x+y)α2xα2yα2]+xn=1(1)nynn!(s+1+n)12x2n=1(1)nynn!(s+2+n)+O(x3),A 49
    which allow us to re-express equation (A 47) as
    ϕ˙X=γXϕX+1vX(1vX)(α2)xα2Γ(2α,x)+vX(1vM)(α2)xn=1(1)nξnn!(3α+n)12vX(1vM)(α2)x2n=1(1)nξnn!(4α+n)+vX(1vM)(α2)Γ(2α)[(x+ξ)α2xα2ξα2]+O(ϕX3).A 50
    Expanding Γ(2α,βX0uϕX) and regrouping everything we are finally left with
    ϕ˙X=γXϕX(1vX)(α2)[x(3α)x22(4α)]+vX(1vM)(α2)x[n=1(1)nξnn!(3α+n)]+12vX(1vM)(α2)x2[n=1(1)nξnn!(4α+n)](1vXvM)(α2)Γ(2α)xα2+vX(1vM)(α2)Γ(2α)[(x+ξ)α2ξα2]+O(ϕX3).A 51
    The goal now is to rewrite the infinite series back in terms of gamma functions using equation (A 48). To do so, one can add and subtract the missing n = 0 terms that complete the series and use the identity Υ(s,x)=sΓ(s)Γ(s,x) [70] to obtain
    ϕ˙X=γXϕX(1vXvM)[x(α2)(3α)x2(α2)2(4α)]+vX(1vM)(α2)xξα3Υ(3α,ξ)12vX(1vM)(α2)x2ξα4Υ(4α,ξ)(1vXvM)(α2)Γ(2α)xα2+vX(1vM)(α2)Γ(2α)[(x+ξ)α2ξα2]+O(ϕX3).A 52
    We remark that, although the lower incomplete gamma function Υ(s,x)=0xts1etdt is in principle defined for Re(s) > 0, it is possible to extend its domain to any non-integer s < 0 by using Υ(s+1,x)=sΥ(s,x)xsex [71]. Finally, although we have expanded ϕ˙ up to second order, the above expression can be easily generalized to any desired cut-off order imax, writing it in a more compact way
    ϕ˙X=γXϕX+(α2)i=2imax+1(1)ixi1(i1)![(1vXvM)(αi1)+vX(1vM)ξαi1Υ(i+1α,ξ)]+(1vMvX)xα2Γ(3α)vX(1vM)Γ(3α)[(x+ξ)α2ξα2]+O(ϕXimax).A 53
    Now, to obtain the epidemic threshold we analyse the stability of the quiescent fixed point ϕ0 = 0, knowing that if ϕ˙X=f(ϕX), then λ=fϕ|ϕ0=0 gives the bifurcation point. Before, we proceed to study the different cases, note that in the vicinity of ϕ0 one can write the binomial sum (for any real α) as (x+ξ)α2=ξα2+k=1C(α2,k)xkξα2k, where C(s, k) = (s(s − 1)(s − 2) · · · (s + 1 − k))/k!. This can simplify a bitequation (A 53), which can be re-arranged—after grouping together terms of the same order—as
    ϕ˙X=γXϕX+(α2)uβX0[(1vXvM)(α3)vX(1vM)tξα3Γ(3α,ξ)]ϕX(α2)i=3imax+1(1)i(uβX0)i1(i1)![(1vXvM)(αi1)vX(1vM)ξαi1Γ(i+1α,ξ)]ϕXi1+(1vMvX)(uβX0ϕX)α2Γ(3α)+O(ϕXimax).A 54
    Two different cases need to be distinguished depending on whether the networks are scale-free or not

    If 2 < α < 3 (i.e. the network is scale-free) then the term ϕXα2 dominates with

    λ(1vMvX)(α2)xα3Γ(3α),A 55
    and since ϕ0α3, we will find λ = 0 only if, trivially, vX = vM = 1. This is the counterpart of saying that the epidemic threshold vanishes, in agreement with the well-known result in scale-free networks [26,34].


    If otherwise α > 3, the leading order is given by the linear term and we can write

    λ(α2)uβX0[(1vXvM)(α3)vX(1vM)ξα3Γ(3α,ξ)]γX.A 56
    One could, for instance, fix vM and write the epidemic threshold in terms of vX
    vX=vXvM+(1vM)(α3)ξα3Γ(3α,ξ),A 57
    where vX=1(k/RXeffk2)=1((α3)/((uRXeff)(α2))) is the expected epidemic threshold for a SIR model with random vaccination in a power-law degree-distributed network [46]. As a sanity check, note that one can retrieve the expected epidemic threshold for an independent disease in two different ways: by setting vM = 1, or by letting ξ ≈ 0 and taking the asymptotic behaviour of the upper-incomplete gamma function, ξα3Γ(3α,ξ)(α3)1. It is also not difficult to show that (as expected) vXvXSIR. Using ξα3Γ(3α,ξ)eξα3 and sMu(T)1vM one can write
    vM+(1vM)(α3)ξα3Γ(3α,ξ)vM+(1vM)eξ=vM+sMu(T)1,A 58
    which proves that the denominator in equation (A 57) is always smaller than or equal to one. On the other hand, it is also possible to derive the exact value of the epidemic threshold in terms of the fraction of M-vaccinated individuals, recovering equation (A 43). Remarkably, both vX and vM can be computed numerically using the relation between vM and ξ given by equation (A 40), although in the main text the analysis is limited to the latter quantity.

    Finally, equation (A 54) can help us to discern the nature of the bifurcation within our HMF approach. Developing the sum up to second-order terms in ϕX

    ϕ˙X=γXϕX+(α2)uβX0[(1vXvM)(α3)vX(1vM)ξα3Γ(3α,ξ)]ϕX12(α2)(uβX0)2[(1vXvM)(α4)vX(1vM)ξα4Γ(4α,ξ)]ϕX2+(1vMvX)(uβX0ϕX)α2Γ(3α)+O(ϕX3).A 59

    If 3 < α < 4: ϕ˙aϕ|c|ϕα2 with b ≠ 0, which is the form of a transcritical bifurcation. ϕ˙=0 admits therefore a solution with arbitrary small ϕ ≠ 0 when we are slightly above the epidemic threshold, and the transition is continuous.


    If 4 < α: ϕ˙aϕ+bϕ2. Using that ξα4Γ(3α,ξ)eξ/α4 [71] and sMu(T)1vM one can write

    (1vXvM)(α4)vX(1vM)ξα4Γ(4α,ξ)1vXvMvXsMu(T)(α4)1vX(α4)0.A 60
    It follows that b < 0 and the transition is therefore continuous.

    Appendix B. The effect of partial immune amnesia

    In the following section, we show that a more general model in which immunity to disease X is lost with some probability p (rather than with certainty) is qualitatively similar to the case p = 1 that has been discussed throughout the paper. In particular, we prove mathematically that an extended version of the SIR-IA model (in the case with demography and homogeneous underlying networks) presents a vaccination threshold v~M that is analogous in form to the one derived in the main text. Figure 8a shows a sketch of the extended version of the model: an individual in the IR state now moves to the RS state (i.e. loses immunity to X disease due to IA) with a rate M or becomes fully immune (RR) with rate (1 − p)γM.

    Figure 8. Analysis of the extended version of the model in which immune amnesia effects occur only with some probability. (a) Sketch of transitions between all possible states. As in the main text, green (orange) cells correspond to measles (X) disease states, and demographic processes are not included. (b) Analytically determined phases for the X disease in the extended SIR-IA model for the case with demography on homogeneous networks. The white-dashed lines which would separate in each case the quiescent and IA-induced endemic states were determined through equation (B 6) for different values of the probability p.

    Thus, it is straightforward to see that only equations (3.1e) and (3.1h) from the main set of dynamical equations in the main text need to be modified

    ρ~˙RS=βXρRS(ρSI+ρRI)+γM(ρIS+pρIR)μρRS+μvM(1vX)andρ~˙RR=γXρRIμρRR+μvMvX+γM(1p)ρIR}B 1

    Since the dynamics of the measles outbreak is left unchanged under this modification, the stationary states ξQ = (0, 1 − vM) and ξE=((μ(R0M(1vM)1))/βM,1/R0M) remain as calculated in appendix A. For disease X on the other hand, one can write

    ψ~Q=(i~XQ,s~XQ)=(0,sXQγM(1p)μρIR)andψ~E=(i~XE,s~XE)=(iXEγM(1p)γX+μρIR,1R0X),}B 2
    The stationary states for the above equations are
    ψ˜Q=(i˜XQ,s˜XQ)=(0,sXQγM(1p)μρIR)andψ˜E=(i˜XE,s˜XE)=(iXEγM(1p)γX+μρIR,1R0X),} B 3
    where sXQ and iXE denote the corresponding values in the same stationary states when p = 1, and ρIR* is the fraction of measles-infectious individuals that are resistant for disease X at the end of the measles outbreak. To estimate this last fraction, let us note that
    dρIRdt=0ρIR=βMρSRγM+μiM=vXiM,B 4
    where ρSR=vX/R0M, as derived in appendix A. Taking the more interesting situation in which there is an endemic measles stationary state given by ξE, we can finally obtain the epidemic threshold condition by setting s~XQ=sXQ(γM(1p)/μ)vXiM=1/RoX with sXQ=sX,2Q and iM=iME (for a detailed derivation of all these previously calculated quantities see appendix A)
    1RoX=γMβM(R0M(1v~M)1)+(1vX)(v~M+1R0M)(γM(1p)μ)(vXμ(R0M(1v~M)1)βM).B 5
    Finally, solving for v~M
    v~M=(1/R0X)(1vX)(1/R0M)+(γ~M/βM)(1R0M)1vX(R0Mγ~M/βM),B 6
    where we have defined γ~MγM(1vX(1p)) to recover an expression for the vaccination threshold that has the same form as equation (A 9) in the main text. Figure 8b shows the different dynamical regimes that can be found for different values of p, using the above equation to determine the boundary between the purely endemic and IA-induced endemic states. Observe how, as p becomes smaller (i.e. losing immunity through IA effects becomes less likely), the region in which one can find an IA-induced endemic state (delimited by the white dashed line in each case) shrinks. Nevertheless, from a qualitative point of view the system behaves in the same way as before and one can still find an IA-induced endemic phase in which herd immunity can be lost despite vaccination efforts.

    Finally, let us remark that one could also consider partial IA in such a way that affected individuals have a reduced (but non-zero) re-infection rate to other diseases. Similarly, it is possible to account for spontaneous waning immunity (at the risk of losing mathematical tractability) just by introducing new rates for transitions from recovered to susceptible states. These extensions are nevertheless beyond the scope of the present work.

    Published by the Royal Society under the terms of the Creative Commons Attribution License, which permits unrestricted use, provided the original author and source are credited.