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.

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 [8-10], conclusive empirical evidence has been very recently found that measles can wipe out acquired immunity to other infectious diseases through a mechanism called IA [11-13]. 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,11-13]. 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 [19-22]; (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.

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—[23-27] 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 [23-26], 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, β, β and γ, γ denote the infectivity and recovery rates for measles and disease X, respectively. Note that IA is explicitly implemented within the term +γρIR in equation (3.1), which drags X-recovered individuals back into the pool of X-susceptible ones, in the case they had been infected with measles. The parameters, v and v 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, v = 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 [29–31] 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.

	symbol	base-line value
measles infectivity rate	βM	2.1 d−1
measles recovery rate	γM	1/8 d−1
COVID-19 infectivity rate	βX	0.25 d−1
COVID-19 recovery rate	γX	1/14 d−1
immigration/emigration rate	μ	1/365 d−1
COVID-19 vaccination coverage	vX	0.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.

Results

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 ρ is the population fraction in state ij.

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 v and v is to determine the initial fraction of vaccinated population of M and X diseases, respectively. To simplify the forthcoming mathematical analyses, let us define i = ρIS + ρIR and s = ρSS + ρSR as the total fraction of infectious and susceptible individuals of measles, and their counterparts i = ρSI + ρRI and s = ρ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 v vaccinated against X is sufficiently high to initially avoid spreading the disease (i.e. i ≈ 0 prior to the onset of a measles outbreak). Note that, under this assumption, from equations (3.1), (3.1), (3.1) and (3.1) one readily derives the standard SIR-model mean-field equations for measles dynamics and 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 di/dt > 0, and decays to 0 if this is negative. Therefore, since s(0) = (N (1 − v) − ε)/N ≈ 1 − v 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 —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, 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 and , 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.1), (3.1), (3.1) and (3.1) an epidemic-threshold condition analogous to equation (3.4) which results into and 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-infectious individuals, we let an outbreak of measles spread through the population (for which we set ); once it fades out, we add at some initial time T a second seed consisting of ε X-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 v is only slightly above the herd-immunity threshold (i.e. ) 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, v: as naively expected, (i) the total outbreak size grows and (ii) the time elapsed between measles and X outbreak peaks becomes shorter as v is reduced. Observe also that X epidemics can only break out (i.e. become supercritical) if v is set below a minimum, critical vaccination threshold , 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 , with ɛ = 0.001. Conversely, was set to guarantee the spreading of measles outbreaks. At time T = 100 days, a seed 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 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 , below which the disease is always endemic, independently of IA effects (black-dashed line). Note that a triple point appears for and (marked by the white arrow), at which the three phases collide.

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.1)–(3.1) with μ > 0 (see table 1) and X-vaccination coverage , 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 v drops below a certain measles herd-immunity level . 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 v drops below a certain critical threshold value . Such an X-disease endemic state is purely induced by IA effects, i.e. it is a IA-induced endemic phase.

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 v (vanishing values in the disease-free steady states are not shown). Blue and green dashed lines mark the analytically found critical points and , 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.1)–(3.1). (b) Minimum fraction of M-vaccinated individuals, , necessary to prevent an endemic state of disease X as a function of the reproduction number . The black curve represents the theoretical approximation for 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.

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 , below which the disease is always endemic, independently of IA effects (black-dashed line) and (marked by the pink arrow), above which the epidemics does not propagate despite the IA contribution. Note that a triple point appears for and (marked by the white arrow), at which the three phases collide.

We also analysed the dependence of on by solving for the stationary states of equations (3.1)–(3.1); 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 . Observe in particular, that as , approaches that of the critical point for measles outbreak , 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 (). Similarly, figure 5 shows the fraction of infected individuals in the stationary state (colour coded) as a function of both, v and v. 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 , but in both cases—even in the limit of stringent vaccination policies for the X disease, v ≈ 1—an IA-induced endemic state can appear if M-vaccination drops below a critical threshold level . The exact value of this threshold (which marks the boundary of the stable endemic region) will depend on v. 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 v value, as given by equation (3.10); population size N = 105. (a) Results for a mildly infectious disease, with 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. Each point in the grid was obtained numerically by solving for the fixed points of equations (3.1)–(3.1).

In summary, the value of the transition point has been numerically shown to depend on both the infectivity and the vaccination coverage v 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.

SIR-IA model with demography: theoretical results

Let us first analyse the possible stable fixed points of equations (3.1)–(3.1) as a function of v. Starting from equations (3.1), (3.1), (3.1) and (3.1) one can write and Assuming, as above, initial herd immunity for the X disease, i.e. , one can approximate i ≈ 0 and readily find the two steady-state solutions of equations (3.6) and (3.7): a disease-free one, (i*, s*) = (0, 1 − v) and an endemic one . A standard linear stability analysis allows us to recover the existence of a transcritical bifurcation at 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 i and s and

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 i of measles infectious population. However, assuming that a steady state for measles has been reached, it is possible to show that ; using this and searching for steady-state solutions of equations (3.8) and (3.9), one can compute the minimum value of v preventing the existence of a stable X-endemic state, thus obtaining an expression for (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, μ ≪ γ, 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 rather that —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 plane.

Observe that, since we assumed , then trivially must hold, with both critical points coinciding for . Similarly, as , one readily sees that the existence of a non-trivial second threshold is limited by the constraint (these two limits are marked with arrows in figure 6). Therefore, an immune-amnesia-induced endemic phase for .

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.1)–(3.1), 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 , 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 and v (figure 5).

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,33-39]. 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,40-42] (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 [43-46]). (ii) As in the previous analyses, we impose that 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 β = β0〈k〉, with and . Within this convention, the values presented in table 1 are recovered for an average network degree 〈k〉 ≈ 20. These same values of and 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 v. It can be readily noticed that the vaccination threshold 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 v = 0.9 to v = 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, , 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 p ∼ k−, numerical simulations reveal that reducing α leads to larger values of , 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 [47-49].

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 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 . Defining the auxiliary function ζ(q) = q v(1 + 〈k〉q) e−〈 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 . 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 —where u is now the minimum degree of a node in the network—and , 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 v = v = 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,51-53]: 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 (see figure 7d, where for v = 0.9 a minimum M-vaccination level preventing an X outbreak is defined only in networks with ). 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 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 [54-56].

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.

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 [17-22]; 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.

Material and methods

The steady-state solutions reported in figures 4b and 5 were obtained by solving equations (3.1)–(3.1) 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.1)–(3.1) 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 − v) (1 − v) − (ɛ + ɛ), ρSI = ɛ, ρSR = 0, ρIS = ɛ, ρIR = 0, ρSR = (1 − v)v, ρRS = v(1 − v), ρRR = vv}. To obtain the deterministic solutions (figures 4b and 5), we set ɛ = ɛ = 50/N, while for the stochastic simulations we chose ɛ = 1/N and ɛ = 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].