Immune amnesia induced by measles and its effects on concurrent epidemics

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 [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.

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—[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:

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, v_X and v_M 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_X = 0.9.

• Download figure
• Open in new tab
• Download PowerPoint

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]:

$$\begin{array}{rl}{\dot{\rho }}_{\mathrm{SS}}=& -{\beta }_X{\rho }_{\mathrm{SS}}({\rho }_{\mathrm{SI}}+{\rho }_{\mathrm{RI}})-{\beta }_M{\rho }_{\mathrm{SS}}({\rho }_{\mathrm{IS}}+{\rho }_{\mathrm{IR}})-\mu {\rho }_{\mathrm{SS}}\\ & +\mu (1-v_X(1-v_M)-v_M(1-v_X)-v_M{v}_X),\end{array}$$3.1a

$${\dot{\rho }}_{\mathrm{SI}}={\beta }_X{\rho }_{\mathrm{SS}}({\rho }_{\mathrm{SI}}+{\rho }_{\mathrm{RI}})-{\gamma }_X{\rho }_{\mathrm{SI}}-\mu {\rho }_{\mathrm{SI}},$$3.1b

$$\begin{array}{rl}{\dot{\rho }}_{\mathrm{SR}}=& -{\beta }_M{\rho }_{\mathrm{SR}}({\rho }_{\mathrm{IS}}+{\rho }_{\mathrm{IR}})+{\gamma }_X{\rho }_{\mathrm{SI}}-\mu {\rho }_{\mathrm{SR}}\\ & +\mu v_X(1-v_M),\end{array}$$3.1c

$${\dot{\rho }}_{\mathrm{IS}}={\beta }_M{\rho }_{\mathrm{SS}}({\rho }_{\mathrm{IS}}+{\rho }_{\mathrm{IR}})-{\gamma }_M{\rho }_{\mathrm{IS}}-\mu {\rho }_{\mathrm{IS}},$$3.1d

$$\begin{array}{rl}{\dot{\rho }}_{\mathrm{RS}}=& -{\beta }_X{\rho }_{\mathrm{RS}}({\rho }_{\mathrm{SI}}+{\rho }_{\mathrm{RI}})+{\gamma }_M({\rho }_{\mathrm{IS}}+{\rho }_{\mathrm{IR}})\\ & -\mu {\rho }_{\mathrm{RS}}+\mu v_M(1-v_X),\end{array}$$3.1e

$${\dot{\rho }}_{\mathrm{IR}}={\beta }_M{\rho }_{\mathrm{SR}}({\rho }_{\mathrm{IS}}+{\rho }_{\mathrm{IR}})-{\gamma }_M{\rho }_{\mathrm{IR}}-\mu {\rho }_{\mathrm{IR}},$$3.1f

$${\dot{\rho }}_{\mathrm{RI}}={\beta }_X{\rho }_{\mathrm{RS}}({\rho }_{\mathrm{SI}}+{\rho }_{\mathrm{RI}})-{\gamma }_X{\rho }_{\mathrm{RI}}-\mu {\rho }_{\mathrm{RI}}$$3.1g

$$\mathrm{and}\hspace{0.17em}{\dot{\rho }}_{\mathrm{RR}}={\gamma }_X{\rho }_{\mathrm{RI}}-\mu {\rho }_{\mathrm{RR}}+\mu v_M{v}_X$$3.1h

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 v_M and v_X 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_M = ρ_IS + ρ_IR and s_M = ρ_SS + ρ_SR as the total fraction of infectious and susceptible individuals of measles, and their counterparts i_X = ρ_SI + ρ_RI and s_X = ρ_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_X vaccinated against X is sufficiently high to initially avoid spreading the disease (i.e. i_X ≈ 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

$$\frac{\mathrm{d}{i}_M}{\mathrm{d}t}={\beta }_M{i}_M{s}_M-{\gamma }_M{i}_M$$3.2

and

$$\frac{\mathrm{d}{s}_M}{\mathrm{d}t}=-{\beta }_M{i}_M{s}_M.$$3.3

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_M/dt > 0, and decays to 0 if this is negative. Therefore, since s_M(0) = (N (1 − v_M) − ε)/N ≈ 1 − v_M is the initial fraction of M-susceptible to measles, the epidemic threshold separating the above two regimes is specified by the condition

$$v_M^{†}=1-\frac{1}{R_0^M},$$3.4

where, as usual, the basic reproduction number $R_0^M={\beta }_M/(\mu +{\gamma }_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, $v_M^{†}$ 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 $R_0^M\approx 17$ and $v_M^{†}\approx 0.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)

$$v_X^{†}=1-\frac{1}{R_0^X},$$3.5

which results into $R_0^X\approx 3.5$ and $v_X^{†}\approx 0.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 $v_M<v_M^{†}$); once it fades out, we add at some initial time T a second seed consisting of ε_X 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_X is only slightly above the herd-immunity threshold $v_X^{†}$ (i.e. $v_X\gtrsim v_X^{†}$) 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_M: as naively expected, (i) the total outbreak size grows and (ii) the time elapsed between measles and X outbreak peaks becomes shorter as v_M is reduced. Observe also that X epidemics can only break out (i.e. become supercritical) if v_M is set below a minimum, critical vaccination threshold $v_M^{††}(\mu =0)<v_M^{†}$, 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.

• Download figure
• Open in new tab
• Download PowerPoint

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 $v_M^{††}(\mu =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.

• Download figure
• Open in new tab
• Download PowerPoint

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 $v_X>v_X^{†}$, 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_M drops below a certain measles herd-immunity level $v_M^{†}\approx 0.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 v_M drops below a certain critical threshold value $v_M^{††}(\mu \ne 0)$. Such an X-disease endemic state is purely induced by IA effects, i.e. it is a IA-induced endemic phase.¹

• Download figure
• Open in new tab
• Download PowerPoint

We also analysed the dependence of $v_M^{††}$ on $R_0^X$ 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 $v_M^{††}$. Observe in particular, that as $R_0^X\to \mathrm{\infty }$, $v_M^{††}$ approaches that of the critical point for measles outbreak $v_M^{††}\approx v_M^{†}$, 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 ($v_X\gg v_X^{†}$). Similarly, figure 5 shows the fraction of infected individuals in the stationary state (colour coded) as a function of both, v_M and v_X. Results for two different X diseases are shown: a mildly infectious one with (a) R₀ = 1.3 and (b) COVID-19, with an estimated R₀ = 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 $R_0^X$, but in both cases—even in the limit of stringent vaccination policies for the X disease, v_X ≈ 1—an IA-induced endemic state can appear if M-vaccination drops below a critical threshold level $v_M^{††}$. The exact value of this threshold (which marks the boundary of the stable endemic region) will depend on v_X. 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.

• Download figure
• Open in new tab
• Download PowerPoint

In summary, the value of the transition point $v_M^{††}(\mu \ne 0)$ has been numerically shown to depend on both the infectivity $R_0^X$ and the vaccination coverage v_X 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 v_M. Starting from equations (3.1a), (3.1c), (3.1e) and (3.1f) one can write

$$\frac{\mathrm{d}{i}_M}{\mathrm{d}t}={\beta }_M{i}_M{s}_M-{\gamma }_M{i}_M-\mu i_M$$3.6

and

$$\frac{\mathrm{d}{s}_M}{\mathrm{d}t}=-{\beta }_M{i}_M{s}_M-{\beta }_X{i}_X{\rho }_{\mathrm{SS}}-\mu s_M+\mu (1-v_M).$$3.7

Assuming, as above, initial herd immunity for the X disease, i.e. $v_X>v_X^{†}$, one can approximate i_X ≈ 0 and readily find the two steady-state solutions of equations (3.6) and (3.7): a disease-free one, (i_M^*, s_M^*) = (0, 1 − v_M) and an endemic one $(\overline{i_M^{\ast },s_M^{\ast }})=(R_0^M((1-v_M)-1)\mu /{\beta }_M,1/R_0^M)$. A standard linear stability analysis allows us to recover the existence of a transcritical bifurcation at $v_M^{†}=1-1/R_0^M$ 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_X and s_X

$$\begin{array}{rl}\frac{\mathrm{d}{s}_X}{\mathrm{d}t}& =-{\beta }_X{s}_X{i}_X-{\beta }_M{\rho }_{\mathrm{SS}}{i}_M+{\gamma }_M{i}_M-\mu s_X+\mu (1-v_X)\end{array}$$3.8

and

$$\frac{\mathrm{d}{i}_X}{\mathrm{d}t}={\beta }_X{s}_X{i}_X-{\gamma }_X{i}_X-\mu i_X.$$3.9

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_M of measles infectious population. However, assuming that a steady state for measles has been reached, it is possible to show that ${\rho }_{\mathrm{SS}}^{\ast }=(1-v_X)/R_0^M$; using this and searching for steady-state solutions of equations (3.8) and (3.9), one can compute the minimum value of v_M preventing the existence of a stable X-endemic state, thus obtaining an expression for $v_M^{††}$ (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:

$$v_M^{††}\approx \frac{1}{v_X}(1-\frac{1}{R_0^X})-\frac{1}{R_0^M}=\frac{v_X^{†}}{v_X}-\frac{1}{R_0^M}$$3.10

which depends not only on the R₀’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 $v_X^{†}$ rather that $R_0^X$—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 $(v_M,v_X/v_X^{†})$ plane.

• Download figure
• Open in new tab
• Download PowerPoint

Observe that, since we assumed $v_X\ge v_X^{†}$, then trivially $v_M^{††}\le 1-1/R_0^M=v_M^{†}$ must hold, with both critical points coinciding for $v_X=v_X^{†}$. Similarly, as $0\le v_M^{††}\le 1$, one readily sees that the existence of a non-trivial second threshold $v_M^{††}$ is limited by the constraint $v_X/v_X^{†}\le R_0^M$ (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 $v_M^{††}$, 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 $R_0^X$ and v_X (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,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 $v_x>v_x^{†}$ 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 β₀ so that β = β₀〈k〉, with ${\beta }_0^X=0.017$ and ${\beta }_0^M=0.14$. Within this convention, the values presented in table 1 are recovered for an average network degree 〈k〉 ≈ 20. These same values of ${\beta }_0^X$ and ${\beta }_0^M$ 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_M. It can be readily noticed that the vaccination threshold $v_M^{††}$ 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_M = 0.9 to v_M = 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.

• Download figure
• Open in new tab
• Download PowerPoint

On the other hand, for power-law degree distributions p_k ∼ k^−α, numerical simulations reveal that reducing α leads to larger values of $v_M^{††}$, 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 $q=s_1^M(T)/(1-v_M)$ 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:

$$q(1-v_M)\hspace{0.17em}{\mathrm{e}}^{⟨k⟩(q-1)}-\log (q)/R_M^{\mathrm{eff}}-(1-v_M)=0,$$3.11

where $R_M^{\mathrm{eff}}={\beta }_0^M/{\gamma }_M$. Defining the auxiliary function ζ(q) = q v_X(1 + 〈k〉q) e^{−〈k〉(1−q)} we were able to derive a closed equation for the vaccination threshold as a function of parameter values

$$v_M^{††}(q)=\frac{1/R_X^{\mathrm{eff}}-(1+⟨k⟩)+\zeta (q)}{\zeta (q)-v_X(1+⟨k⟩)}.$$3.12

This expression, when inserted into equation (3.11), can be numerically solved, leading to an analytical determination of $v_M^{††}$. 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 $\xi =\log ((1-v_M)/s_u^M)$—where u is now the minimum degree of a node in the network—and $v_X^{†}=1-(⟨k⟩/R_X^{\mathrm{eff}}⟨k^2⟩)$, which is the corresponding threshold in a standard SIR model with vaccination [46]. In terms of these quantities, we obtained (see appendix B.3),

$$v_M^{††}=\frac{(\alpha -3){\xi }^{\alpha -3}\mathit{\Gamma }(3-\alpha ,\xi )-v_X^{†}/v_X}{((\alpha -3){\xi }^{\alpha -3}\mathit{\Gamma }(3-\alpha ,\xi )-1)}.$$3.13

with the constraint

$$\begin{array}{rl}(1-v_M)({\mathrm{e}}^{-\xi }& -(\alpha -2){\xi }^{\alpha -2}\mathit{\Gamma }(2-\alpha ,\xi )-1)\\ & +\frac{{\gamma }_M}{u{\beta }_M^0}\xi =0\end{array}$$3.14

needed to determine ξ. Conversely, for scale-free networks (i.e. 2 < α ≤ 3) one finds that outbreaks do not propagate only if v_M = v_X = 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 α > 3 (see figure 7d, where for v_X = 0.9 a minimum M-vaccination level preventing an X outbreak is defined only in networks with $\alpha \gtrsim 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 $v_M^{††}$ 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.

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 [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.

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 ρ₀ = {ρ_SS = (1 − v_M) (1 − v_X) − (ɛ_X + ɛ_M), ρ_SI = ɛ_X, ρ_SR = 0, ρ_IS = ɛ_M, ρ_IR = 0, ρ_SR = (1 − v_M)v_X, ρ_RS = v_M(1 − v_X), ρ_RR = v_Mv_X}. 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 = 10⁵.

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.

Funding

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.