How immune dynamics shape multi-season epidemics: a continuous-discrete model in one dimensional antigenic space

We extend a previously published model for the dynamics of a single strain of an influenza-like infection. The model incorporates a waning acquired immunity to infection and punctuated antigenic drift of the virus, employing a set of coupled integral equations within a season and a discrete map between seasons. The long term behaviour of the model is demonstrated by examples where immunity to infection depends on the time since a host was last infected, and where immunity depends on the number of times that a host has been infected. The first scenario leads to complicated dynamics in some regions of parameter space, and to regions of parameter space with more than one attractor. The second scenario leads to a stable fixed point, corresponding to an identical epidemic each season. We also examine the model with both paradigms in combination, almost always but not exclusively observing a stable fixed point or periodic solution. Adding stochastic perturbations to the between season map fails to destroy the model’s qualitative dynamics. Our results suggest that if the level of host immunity depends on the elapsed time since the last infection then the epidemiological dynamics may be unpredictable.


Introduction
Respiratory viral illnesses typically cause regular seasonal epidemics in temperate climates.The influenza viruses in particular display strong seasonal behaviour (Webster et al. 2013), but other viruses such as RSV, coronaviruses and rhinoviruses are also etiological causes of seasonal acute respiratory illness (Howard et al. 2013;Mandell 2005).As SARS-CoV-2, the virus that causes Covid-19, transitions out of the pandemic phase-characterised by transient dynamics exhibiting large scale oscillations driven by a range of biological, behavioural and environmental factors-it too may settle into a seasonal cyclic pattern.
Viral infection stimulates the human immune system.Taking influenza as an example, infection elicits a strain-specific antibody immune response, that contributes to the resolution of infection and provides long lasting strain-specific protection such that subsequent exposure is unlikely to result in a productive infection.A similar immune response occurs for other pathogens, including SARS-CoV-2 (Dan et al. 2021;Guo et al. 2022), albeit with evidence that strain-specific immunity may also wane through time (Khoury et al. 2021).Strong strain-specific immunity, combined with the inherent erroneous replication of RNA viruses (Duffy 2018) drives viral evolution (entitled 'antigenic drift' for influenza), selecting for immune-escape variants (Guarnaccia et al. 2013).In consequence, for many viral infections, while strain-specific immunity is long (perhaps even life-long), protection against the circulating viral strains may be short lived (Fonville et al. 2014;Grenfell et al. 2004;Smith et al. 2004).We are thus infected with influenza or the common cold on multiple occasions during our life.SARS-CoV-2 shows clear signs of also establishing such transmission patterns (Callaway 2022), which are unsurprising on theoretical grounds (Lavine et al. 2021).
From a dynamical systems perspective, multi-season epidemic dynamics arise from the gain of immunity as a particular strain spreads, and the loss of protection due to the induced selection pressure on the virus leading to the emergence of new immune variants, and/or waning of the host's strain specific immune response.Efforts to understand this characteristic gain and loss of immunity, and the characteristic dynamics of these systems have been the focus of a number of mathematical models.Kucharski et al. (2015Kucharski et al. ( , 2018) ) investigated antibody dynamics, and how previous exposure influences a host's response to a new strain, but did not consider the consequential epidemic dynamics.Andreasen (2003) modelled multi-season epidemic dynamics, and in par-ticular a scenario in which variable infectiousness was determined by the time since exposure.In earlier work, through construction of a hybrid continuous-time-discretemap model, we have considered how immunity obtained in a previous season, and possibly lost in the inter-epidemic period, influences multi-season epidemic dynamics (Roberts et al. 2019).For this model the within season epidemic was represented by a system of ordinary differential equations and the between season changes by a discrete map.Hence a fixed point of the map was equivalent to a periodic solution of the model.Our analysis demonstrated that epidemic dynamics were highly sensitive to parameter choices, with steady-state (i.e.period 1), periodic and complicated dynamics displayed.In another study Kucharski et al. (2016) introduced a general framework for multi-strain epidemic dynamics, in which cross-immunity between strains was modelled based on past exposure history.Noting the combinatorial complexity of considering all possible exposure histories (a challenge that is playing out with SARS-CoV-2 right now (Vattiatio et al. 2022)), they proposed two alternative approaches to reducing the dimensionality of the general system, known as 'exposure' and 'history' based models (Kucharski et al. 2016).
In the present paper we provide a general model for the within-and between-season dynamics of a respiratory virus.We maintain the structure used by Andreasen (2003) and in our previous work (Roberts et al. 2019): the within season dynamics are represented by a continuous model on t ∈ (0, ∞), and the between season dynamics by a discrete map.First we present a within-season Kermack-McKendrick type model, with the host's susceptibility to infection taking values from a continuous density.We present a stochastic model for the between-season dynamics, reflecting population turnover and loss of immunity due to waning protection and antigenic drift.
We then specialise the model, with the host susceptibility taking discrete values and deterministic between-season dynamics.In the same vein as Kucharski et al. (2016), the long term behaviour of this model is demonstrated by examples where immunity to infection depends on the time since the host was last infected, and where immunity depends on the number of times that a host has been infected.

The within season model
We assume a population of fixed size.Let the proportion of hosts who have never been infected with the virus at time t be S ∅ (t).We assume that those that have been infected have a lesser susceptibility to infection with the virus.We index that susceptibility by a variable θ , so that when experiencing the same force of infection, the probability that an individual becomes infected is k(θ ) times the probability that an individual in S ∅ would be infected if exposed to the same degree.We assume k(0) ≤ 1 and k (θ ) ≤ 0, where the prime denotes the derivative.We restrict θ to the range [0, 1], so that the proportion of the population with θ in the range [a, b] at time t is b a S(t, θ) dθ and S ∅ (0)+ 1 0 S(0, θ) dθ = 1.The model does not have a removed compartment as such, although we can model removal by setting k(θ ) = 0 in some interval θ ∈ (1 − , 1] for 0 < 1. Let an epidemic take place within a time period represented by t ∈ [0, ∞), in the sense that the first transmission of infection occurs at or after t = 0, and all transmission has ceased as t → ∞.During an epidemic the proportion of hosts who are fully susceptible and the density of hosts with susceptibility index θ decrease according to where ı ∅ (t) is the incidence of infection in fully susceptible hosts, and ı(t, θ) is the incidence of infection in hosts with susceptibility θ at time t.If the epidemic is precipitated by imported cases with incidence j ∅ (t) or j (t, θ), according to their prior status, and the probability mass function of infectiousness (PMF) with time since infection is p(τ ), then where R 0 is the basic reproduction number (Diekmann et al. 2013) and the force of infection is For small t we approximate The susceptibility profile of the population at the end of an epidemic may be found from S ∅ (∞) = S ∅ (0)e −R 0 P and S(∞, θ) = S(0, θ)e −R 0 k(θ)P , where P is the final size of the epidemic (proportion of the population infected throughout the epidemic) and P solves The derivation of Eq. 3 may be found in Appendix A.

The stochastic between season model
Let the PMF of θ post-infection for those that were infected in a season and had previously not been infected be f (θ ).Let the PMF of θ post-infection for those that were infected in a season and had previous level of immunity ξ be g(θ, ξ ).Finally, let the PMF at the end of the season for those that were not infected and had previous level of immunity ξ be h(θ, ξ ).By definition , and 1 0 h(θ, ξ )dθ = 1.As all those who are infected within a season have increased (strictly not decreased) their level of immunity to infection, we have f (0) = 0 and g(θ, ξ ) = 0 if θ < ξ.Those who are not infected within a season do not increase their level of immunity to infection, hence h(θ, ξ ) = 0 if θ > ξ.In the absence of population turnover (i.e.demographic processes of births and deaths) and antigenic drift the PMF of susceptibility at the beginning of the next season would be where Now use the subscript n to signify this season (prior to the epidemic if it occurs), and n + 1 for next season.Allow a proportion 1 − c of the population to be replaced with fully susceptible hosts between seasons, to allow for population turnover and antigenic drift.The initial conditions for next season are The lower limit of the integral above reflects that any individual for whom θ = 0 is assumed to have the same susceptibility as those in S ∅ , and the conservation rule allocates them to that compartment.That situation can only arise in the model through the action of the integral with kernel h(θ, ξ ).The process outlined above determines the profile of population susceptibility at the start of year n + 1.The relative susceptibility function k(θ ) then determines the dynamics during that year.The conditions on k are k( 0 Some examples of the functions f , g and h are discussed below.We first suggest suitable functions that could apply if the time since the most recent infection determines a host's susceptibility.We then suggest functions that could apply if the number of past infections determines a host's susceptibility.Finally, we discuss the general case when both of these mechanisms apply.

The time since most recent infection determines susceptibility
For this example all those who are infected during the year increase their level of immunity to θ = 1, the maximum level, and those not infected reduce their level of immunity.An example function for f could be where 0 < 1, a and b are positive constants, and B is a beta function with the added constraint that < 1 − ξ to ensure that θ > ξ.An example function for h could be This specifies a function in the (θ, h) plane centred at θ = yξ with base 2 ξ that integrates to one.The positive constants a and b do not necessarily take the same values as in the definitions of f (θ ) and g(θ, ξ ) above.In general, y could be a function of ξ .To ensure that θ > 0 we require y > , and to ensure that θ is reduced in the absence of infection we require y + < 1.In the limit → 0 we obtain h(θ, ξ ) = δ(θ − yξ).
The actions of f , g and h are illustrated in Fig. 1a.We have not discussed how the value of θ may be determined by the PMFs f , g and h.For a purely deterministic model we would take the limits → 0, realising delta functions.In this situation, and where y is a constant, the only levels of θ realised are θ = 1 and θ = y n , where the last infection occurred n + 1 seasons ago.For this limiting case, if θ does not vary continuously at some initial time, it can only take discrete values at future times.

The number of infections determines susceptibility
For this example an individual's immunity level is increased at each infection until a maximum value θ = 1 is reached.Those who were infected during the year increase their level of immunity to θ > ξ, and those not infected retain their level of immunity.An example function for f could be for some x < 1/ (1 + ).We have lim →0 f (θ ) = δ(x), and require < 1 to ensure θ > 0. An example function for g could be where z is a function with 0 < z(ξ ) < 1/ (1 + ), and z (ξ ) > 0 to ensure that θ > ξ.
One suitable function would be z(ξ ) = x + (1 − x) ξ , and more generally x could be a function of ξ .The function h(θ, ξ ) = ξ .The actions of f , g and h are illustrated in Fig. 1b.

The general case
In the general case an individual's immunity level is increased at each infection until a maximum value θ = 1 is reached, those who were infected during the year increase their level of immunity to θ > ξ, and those not infected reduce their level of immunity.The functions f and g could be as in the second example (Eqs.5 and 6), and the function h could be as in the first example (Eq.4).Their actions are illustrated in Fig. 1c.

The deterministic continuous-discrete model of multi-season epidemics
For the situations considered in Sects. 2 and 3, if θ does not vary continuously at some initial time, θ can only take discrete values.Hence the model can be recast with m discrete compartments S (t) instead of S(t, θ) depending on a continuous variable θ .
The function k(θ ) is then replaced with a number of factors k ≤ 1, with k +1 ≤ k .The functions f (θ ), g(θ, ξ ), and h(θ, ξ ) that map points to PMFs are replaced with functions F , G and H that map points to points, as per the illustration in Fig. 1.

Within-season
The within-season model becomes The epidemic takes off if R > 1 where The total proportion of the population infected in the epidemic is where S ∅ (∞) = S ∅ (0)e −R 0 P , S (∞) = S (0)e −R 0 k P , and the final size of the epidemic solves This is a generalisation of the model analysed in Roberts et al. (2019).Differentiating Eq. 9 by S j (0

Between-seasons
Let those who were infected for the first time this season enter the compartment S u , those who were previously in compartment S who were infected this season enter the compartment S v( ) with v( ) ≥ , and those who were previously in compartment S who were not infected this season enter the compartment S w( ) with w( ) ≤ .In the absence of population turnover and antigenic drift, the proportion of the population in the S σ compartment at the beginning of the next season would be F σ + G σ + H σ where otherwise and To account for population turnover and antigenic drift, we allow a proportion 1 − c of the population to be replaced with fully susceptible hosts between seasons.As in Sect. 3 we use the subscript n to signify this season, and n + 1 for next season.The initial conditions for season n + 1 are We can write this mapping in vector form.Define s n to be the vector whose th component is S n (0) for = 1 . . .m, and E(s n ) = e −R 0 P n where P n is the final size of the within season epidemic with initial conditions s n .There is no need to explicitly model the initial condition for the fully susceptible host compartment, as S ∅ n (0) = 1 − s n 1 .The mapping becomes where C is an m × m matrix, and q is an m dimensional vector valued function.The fixed point of the map solves where E = exp(−R 0 P ) and Stability of the map depends on the eigenvalues of the Jacobian matrix, evaluated at s = s .Differentiating Eq. 10 we find that J has components where s is the th component of the vector s n .The fixed point is stable if the eigenvalues of J(s , E ) have absolute value less than one.

Numerical simulations of the deterministic continuous-discrete model
We now explore four scenarios.For each we provide a general specification, with explicit expressions and numerical results for a host population that is either fully susceptible or has one of four levels of immunity (m = 4).For the numerical results we present orbit diagrams showing the annual final size P (see Eq. 8) and effective reproduction number R (see Eq. 7) as the parameter k 1 varies from zero to one, with k 2 = k 2 1 , k 3 = k 3 1 and k 4 = 0; with c = 0.9 or c = 0.7; and with R 0 = 2.0 or R 0 = 4.0.Recall that the parameter k specifies the relative susceptibility of the compartment S and c accounts for between season population turnover and antigenic drift.Hence c = 0.9 implies that 10% of the population is effectively replaced with fully susceptible hosts between each season, and k 4 = 0 means that the S 4 class is immune to infection.In most cases the initial conditions are S ∅ 1 (0) = 1, S 1 (0) = 0 for = 1 . . .m, but we also investigate alternative conditions to check for the existence of multiple attractors.

The time since most recent infection determines susceptibility
For this example we choose u = m, v( ) = m for all and w( ) = − 1 for = 2 . . .m.This leads to where we have used S ∅ n (0) + m =1 S n (0) = 1 and k m = 0.The between season map s n+1 = C(s n )s n + q(s n ) is defined as follows.The matrix C(E) has components otherwise and the vector q(E) has all components zero apart from q m = c (1 − E).The fixed point has components and The map is illustrated for m = 4 below (recall E k 4 = 1).
The fixed point has components Numerical results are presented in Fig. 2. The orbit diagram Fig. 2a shows the proportion of the population infected each year P as k 1 varies from zero to one.The diagram reveals that solutions tend to a fixed point (same size epidemic each year) for 0.39 < k 1 < 0.63, but have complicated (high-period, non-recurrent or chaotic) dynamics for 0.06 ≤ k 1 ≤ 0.39.For k 1 > 0.63 if an epidemic occurs it has the same size as in previous years, but some years there is no epidemic.Figure 2b shows the variability of the effective reproduction number R as k 1 varies from zero to one.The figure reveals that for k 1 > 0.63 there is an epidemic in alternate years, but for k 1 < 0.06 there is a five year cycle, with four years without an epidemic, followed by an epidemic in the fifth year.Figure 2c, d show the dynamics with the same parameter values but with a different set of initial conditions.These diagrams reveal the existence of an alternative attractor for small values of k 1 .Figure 2e-g show the dynamics with alternative values of R 0 or c.For these figures only the dynamics of the effective reproduction number R are shown.No sensitivity on initial conditions was observed with these parameter values.A region of complicated dynamics was observed for the example with c = 0.9 and R 0 = 4.0 (Fig. 2f), but no such region was found for the examples with c = 0.7 (Fig. 2e, g).

The number of infections determines susceptibility
We define a finite number of compartments S , with 1 ≤ ≤ m.If each infection results in a host increasing its immunity to the next level, and not being infected results in immunity staying at the same level, then we have u = 1, v( ) = + 1 for = 1 . . .m − 1 and w( ) = for = 1 . . .m.This leads to where we have used S ∅ n (0) + m =1 S n (0) = 1 and k m = 0.The between season map s n+1 = C(s n )s n + q(s n ) is defined as follows.The matrix C(E) has components : otherwise and the vector q(E) has all components zero apart from q 1 = c (1 − E).The map is illustrated for m = 4 below.
A fixed point of this map is where s 1 can be found by substituting in Numerical results are presented in Fig. 3.The figure is structured to enable direct comparison with the results shown in Fig. 2. In particular, for all of the parameter values chosen the attractor was a fixed point (period one) and no dependency on initial conditions was observed.

The number of infections and the time since last infection determine susceptibility
We define a finite number of compartments S , with 1 ≤ ≤ m.Then define u > 0, v( ) ≥ with v(m) = m, and w( ) < ; u, v, w all integer.First we consider a scenario where infection increases immunity more than noninfection results in a decrease.For example, suppose that having been infected increases the immunity level by two steps, but without infection immunity drops one level between seasons.We would then have Numerical results are presented in Fig. 4a-d.We only present orbit diagrams for the changes in effective reproduction number R, as these are sufficient to provide a description of the dynamics.It can be seen that for the examples with c = 0.9 (Fig. 4a, c) there is a period three solution for low values of k 1 , comprising two seasons with no epidemic (R < 1) followed by one season with an epidemic.For the example with R 0 = 4.0 and c = 0.9 there is also a region of k 1 with complicated (high period) dynamics (Fig. 4c).For all other examples and values of k 1 the observed attractor is a period one solution (fixed point of the discrete map).Now consider a scenario where infection increases immunity less than non-infection results in a decrease.For example having been infected increases the immunity level by one step, but without infection immunity drops two levels between seasons.We would then have u = 1, v( ) = + 1 for = 1 . . .m − 1, and w( ) = − 2 for > 2. For example, with m = 4.
Numerical results are presented in Fig. 4e-h.The orbit diagrams show epidemics in alternate years for some values of k 1 , otherwise the fixed point is stable with an epidemic every year.None of the orbit diagrams show complicated high period dynamics.

Adding a stochastic component
In order to approximate solutions of the stochastic model, we simulated a hybrid model by calculating the final size of the deterministic within-season model (P, Eq. 9) and then adding a small perturbation P → P + δP.The magnitude of the perturbation was taken at random from a uniform distribution δP ∈ [−δ, δ], truncated if necessary so that P ∈ [0, 1].Realisations of the orbit diagrams for the model with immunity determined by the time since last infection, as presented in Fig. 2a and b, but now with the perturbations δP added at each season, are presented in Fig. 5. Comparing Fig. 2a, b with Fig. 5 we see that when δ = 0.05 the stochastic perturbation obscures the period five solution for k 1 < 0.06 and the period one solution for 0.39 < k 1 < 0.63 (Fig. 5a,  b).Orbits in these regions become virtually indistinguishable from those in the region of k 1 that leads to complicated dynamics.The period five and period one solutions are discernible when δ = 0.02 (Fig. 5c, d) and are clearly visible when δ = 0.01 (Fig. 5e,  f).In all of the orbit diagrams shown in Fig. 5 the period two solutions for k 1 > 0.63 are preserved.

Discussion
We previously investigated a model for a seasonal infectious disease where having been infected in one season an individual host was fully protected for the next season, and partially protected for the season after, before becoming fully susceptible again if not infected (Roberts et al. 2019).That model exhibited complicated dynamics for a range of parameter values, as have a number of other models for seasonal infections with different structures (Andreasen 2003;Bacaer and Ouifki 2007;Earn et al. 2000;He and Earn 2007;Mathews et al. 2007;Viboud et al. 2006).The focus of our present study was to investigate if different dynamics would be observed if the level of immunity was determined by the elapsed time since a host was last infected, or by the number of infections a host had experienced.When immunity was determined by the time since infection we observed complicated dynamics over a range of parameter values, especially if acquired immunity resulted in stronger protection (lower values of k 1 in Fig. 2).In contrast, when immunity was determined by the number of infections experienced the model exhibited a steady state, corresponding to an epidemic of the same size every year (Fig. 3).When the two mechanisms were combined the steady state attractor was predominant, although periodic or complicated dynamics were observed for lower values of k 1 in some instances (Fig. 4).These findings have potential implications for predictive modelling studies.For example, our results suggest that if the level of host immunity depends on the elapsed time since the last infection then the epidemiological dynamics may be unpredictable.
The model as presented includes the assumption that epidemiological parameters, and hence R 0 are constant from season to season.However, virus evolution may result in variants presenting with different transmissibility in subsequent seasons.Similarly, the factor allowing for population turnover and antigenic drift has been taken as constant, but for influenza it has been shown that drift can be punctuated by some years of increased change (Smith et al. 2004).Should Covid-19 become seasonal then it is likely that a similar phenomenon would be observed, as different variants of the virus have already exhibited different transmissibility (Manathunga et al. 2023).
Where annual epidemics involve multiple strains of a virus, the timing of the introduction of each strain influences the within season dynamics (Roberts 2012).This would be another source of unpredictability when modelling the long term dynamics of seasonal infections such as influenza, for which some degree of cross-protection between strains has been determined (Ferguson et al. 2003).Modelling epidemics with multiple interacting viruses introduces complications (Gog and Grenfell 2002) and is the subject of ongoing research.Now let s → 0. We have p( 0

Fig. 1
Fig. 1 Illustration of the modes of action used to model the immune response.a the time since last infection determines immunity; b the number of infections determines immunity; c the general case where the two mechanisms combine.The function f (θ ) models the response of previously uninfected hosts to infection, the function g(θ, ξ ) models the response of previously infected hosts to infection, and the function h(θ, ξ ) models the loss of immunity in the absence of infection

Fig. 5
Fig. 5 Orbit diagrams for the model where the time since last infection determines immunity and stochastic perturbations have been added.A,C,E: the annual proportion of hosts infected (P) as a function of the level of partial immunity (k 1 ), the broken line shows an unstable fixed point of the corresponding deterministic model; B,D,F: the effective reproduction number (R) as a function of the level of partial immunity (k 1 ), the horizontal line is at R = 1.A,B: δ = 0.05; C,D: δ = 0.02; E,F: δ = 0.01.Initial conditions s 0 = (0, 0, 0, 0) t) dt = S ∅ (0) − S ∅ (∞) θ) dt = S(0, θ) − S(∞, θ)