Chimera states in multi-strain epidemic models with temporary immunity

TL;DR

This paper explores how time delays and cross-immunity influence the complex dynamics of multi-strain epidemic models, revealing phenomena like chimera states and various synchronization patterns.

Contribution

It introduces a comprehensive analysis of delay-induced bifurcations and the emergence of chimera states in multi-strain epidemic models with different coupling topologies.

Findings

Delay causes Hopf and saddle-node bifurcations affecting disease dynamics.

Cross-immunity induces complex behaviors like traveling waves and chimera states.

Chimera states occur with narrower cross-immunity kernels, impacting disease spread understanding.

Abstract

We investigate a time-delayed epidemic model for multi-strain diseases with temporary immunity. In the absence of cross-immunity between strains, dynamics of each individual strain exhibits emergence and anni- hilation of limit cycles due to a Hopf bifurcation of the endemic equilibrium, and a saddle-node bifurcation of limit cycles depending on the time delay associated with duration of temporary immunity. Effects of all-to-all and non-local coupling topologies are systematically investigated by means of numerical simulations, and they suggest that cross-immunity is able to induce a diverse range of complex dynamical behaviors and synchro- nization patterns, including discrete traveling waves, solitary states, and amplitude chimeras. Interestingly, chimera states are observed for narrower cross-immunity kernels, which can have profound implications for understanding the dynamics of multi-strain diseases.

Equations

$$\begin{array}[]{l}\dot{S}_{i}(t)=\eta-\eta S_{i}(t)-S_{i}(t)\sum\limits_{j=1}^{N}\beta_{j}\sigma_{ij}I_{j}(t)+\gamma_{i}I_{i}(t-\tau)e^{-\eta\tau},\\ \dot{I}_{i}(t)=\beta_{i}S_{i}(t)I_{i}(t)-(\gamma_{i}+\eta)I_{i}(t),\\ \\ \dot{R}_{i}(t)=\gamma_{i}I_{i}(t)-\gamma_{i}I_{i}(t-\tau)e^{-\eta\tau}-\eta R_{i}(t)\\ \\ \hskip 36.98866pt+S_{i}(t)\sum\limits_{j=1,j\neq i}^{N}\beta_{j}\sigma_{ij}I_{j}(t),\\ \end{array}$$

$${\dot{S}}_i(t)+{\dot{I}}_i(t)+{\dot{R}}_i(t)$$

$$\iff {\dot{N}}_i(t)$$

$$\begin{array}[]{l}\dot{S}_{i}(t)=\rho[1-S_{i}(t)]-rS_{i}(t)I_{i}(t)+(1-\rho)I_{i}(t-\tau)e^{-\rho\tau}\\ \\ \hskip 36.98866pt-rS_{i}(t)\sum\limits_{j=1,j\neq i}^{N}\sigma_{ij}I_{j}(t),\\ \\ \dot{I}_{i}(t)=rS_{i}(t)I_{i}(t)-I_{i}(t),\end{array}$$

$$\begin{array}[]{l}\dot{S}(t)=\rho[1-S(t)]-rS(t)I(t)+(1-\rho)I(t-\tau)e^{-\rho\tau},\\ \\ \dot{I}(t)=rS(t)I(t)-I(t).\end{array}$$

$$E_{\tau }^{\ast }=(S_{\tau }^{\ast },I_{\tau }^{\ast })=\left(\frac{1}{r},\rho \frac{r-1}{r}\frac{1}{1-(1-\rho )e^{-\rho \tau }}\right).$$

$$E_{\infty }^{\ast }=\left(S_{\infty }^{\ast },I_{\infty }^{\ast }\right)=\left(\frac{1}{r},\rho \frac{r-1}{r}\right).$$

$${\lambda }^2+\lambda (\rho +r{I}_{\tau }^{\ast })-[r{I}_{\tau }^{\ast }-\rho (r-1)e^{-\lambda \tau }]=0,$$

$$-{\omega }^2+r{I}_{\tau }^{\ast }=[r{I}_{\tau }^{\ast }-\rho (r-1)\cos (\omega \tau )],$$

$$\omega (\rho +r{I}_{\tau }^{\ast })=[r{I}_{\tau }^{\ast }+\rho (r-1)\sin (\omega \tau )].$$

$${\omega }^4+{\omega }^2[{\rho }^2-2r{I}_{\tau }^{\ast }(1-\rho )+r^2(I_{\tau }^{\ast }{)}^2]$$

$$-\rho (r-1)[\rho (r-1)-2r{I}_{\tau }^{\ast }]=0,$$

$${\omega }_{\pm }^2=\frac{1}{2}\left[-{\rho }^2+2r{I}_{\tau }^{\ast }(1-\rho )-r^2(I_{\tau }^{\ast }{)}^2\right]\pm \sqrt{{\left[{\rho }^2-2r{I}_{\tau }^{\ast }(1-\rho )+r^2(I_{\tau }^{\ast }{)}^2\right]}^2+4\rho (r-1)[\rho (r-1)-2r{I}_{\tau }^{\ast }]}.$$

$${\tau }_c=\frac{1}{\omega }\left[\arctan \left(\frac{\omega (\rho +r{I}_{{\tau }_c}^{\ast })}{{\omega }^2-r{I}_{{\tau }_c}^{\ast }}\right)+n\pi \right],n\in \mathbb{N}.$$

$$\begin{array}[]{l}\dot{S}_{i}(t)=\rho[1-S_{i}(t)]-rS_{i}(t)I_{i}(t)+(1-\rho)I_{i}(t-\tau)e^{-\rho\tau},\\ \\ \dot{I}_{i}(t)=rS_{i}(t)I_{i}(t)-I_{i}(t),\\ \\ \dot{S}_{j}(t)=\rho[1-S_{j}(t)]-r\sigma_{ji}S_{j}(t)I_{i}(t),\hskip 8.5359ptj\neq i.\\ \end{array}$$

$$\begin{array}[]{l}\displaystyle{S_{i}^{*}=\frac{1}{r},\hskip 8.5359ptI_{i}^{*}=\rho\frac{r-1}{r}\frac{1}{1-(1-\rho)e^{-\rho\tau}},}\\ \\ \displaystyle{S_{j}^{*}=\frac{1-(1-\rho)e^{-\rho\tau}}{1-(1-\rho)e^{-\rho\tau}+\sigma_{ij}(r-1)},\hskip 8.5359ptI_{j}^{*}=0,\hskip 8.5359ptj\neq i.}\end{array}$$

$$S_1^{\ast }=...=S_N^{\ast }=S_{\mathrm{end}}^{\ast },I_1^{\ast }=...=I_N^{\ast }=I_{\mathrm{end}}^{\ast },$$

$$S_{\mathrm{end}}^{\ast }=\frac{1}{r},I_{\mathrm{end}}^{\ast }=\rho \frac{r-1}{r}\frac{1}{1-(1-\rho )e^{-\rho \tau }+{\sigma }_c},$$

$${\sigma }_{ij}=\exp \left(-\frac{{\left[\frac{N}{2}-\min (|j-i|,N-|j-i|)\right]}^2}{d^2}\right),$$

$${\sigma }_{ij}=\frac{\sigma }{2}\left\{1-\cos \left[2\pi d_p\left(\min \left[|z_j-z_i|,z_N-|z_j-z_i|\right]\right)\right]\right\}$$

$$d_p(z)=z+pz\left(z-\frac{1}{2}\right)(z-1).$$

Full text

multi-strain epidemic model, cross-immunity, coupled oscillators, non-local coupling, chimera states

One of the most fascinating phenomena that has intrigued researchers in the area of nonlinear dynamics for the last fifteen years is a very peculiar pattern of behavior known as chimera states, which is characterized by the simultaneous coexistence of regions of coherent and incoherent dynamics. This pattern was found when identical oscillators were connected with a non-local coupling of high symmetry. In the following years chimera states have attracted a lot of interest and have been studied theoretically and experimentally in a variety of different contexts. This paper investigates how chimera states can appear in epidemic models, and it also explores wider dynamics of multi-strain diseases with time delay and non-local coupling.

I Introduction

Chimera is a hybrid state with coherent and incoherent dynamics, which was first described by Kuramoto and Battogtokh in a system of coupled identical oscillators Kuramoto and Battogtokh (2002). This unusual dynamical pattern was called a chimera state by Abrams and Strogatz Abrams and Strogatz (2004) in light of analogy with a mythological creature with three heads of three different animals. Chimera states have been subsequently discovered in various contexts: SQUID materials Lazarides, Neofotistos, and Tsironis (2015), quantum systems Bastidas et al. (2015), electronic oscillators Gambuzza et al. (2014), and many more. It is currently debated that the dynamics observed, for instance, in uni-hemispheric sleep in mammals and birds Rattenborg, Amlaner, and Lima (2000), and blackouts in power-grids Filatrella, Nielsen, and Pedersen (2008); Martens et al. (2013); Motter et al. (2013) can be interpreted as chimera states. Whilst chimera states have been observed in a number of natural phenomena, they are quite complicated to implement experimentally for several reasons. Firstly, only small networks can be realized in laboratory conditions, and identical oscillators with identical intrinsic frequencies are required Omelchenko et al. (2015). Secondly, chimera states can be very sensitive to initial conditions and often occur only in a small region of the parameter space, and, thus, an experimental setup has to be very precisely controlled in terms of all parameters. Recent studies on two coupled populations of phase oscillators have also demonstrated the possibility of extended basins of attraction Martens, Panaggio, and Abrams (2016), and the existence of chimeras even for small numbers of elements Panaggio et al. (2016). Thirdly, chimera states are a transient state that collapses after a finite period of time into a state of full synchrony Wolfrum and Omel’chenko (2011). Although the lifetime of chimeras has been reported to increase exponentially Wolfrum and Omel’chenko (2011) or as a power-law Olmi (2015); Olmi et al. (2015) in dependence on the number of oscillators, it can be very short for small networks. Despite these challenges, chimera states have been robustly produced in several experiments, including chemical oscillators Nkomo, Tinsley, and Showalter (2013), optical systems Hagerstrom et al. (2012), time-delayed laser networks Larger, Penkovsky, and Maistrenko (2015), electrochemical oscillators Larger, Penkovsky, and Maistrenko (2013), and mechanical oscillators Martens et al. (2013). For a recent review, see Panaggio and Abrams Panaggio and Abrams (2015).

The formation and properties of chimera states have been studied in a number of theoretical models represented as networks of FitzHugh-Nagumo Omelchenko et al. (2013), Kuramoto Olmi (2015); Wang and Li (2011), Ginzburg-Landau Sethia, Sen, and Johnston (2013), van der Pol Ulonska et al. (2016), leaky integrate-and-fire Tsigkri-DeSmedt et al. (2016), Stuart-Landau Zakharova, Kapeller, and Schöll (2016), Hindmarsh-Rose Hizanidis et al. (2014), Hodgkin-Huxley Glaze, Lewis, and Bahar (2016), and SNIPER Vüllings, Schöll, and Lindner (2014) oscillators, as well as many other models. Whilst originally chimera states were discovered in the case of non-local coupling Kuramoto and Battogtokh (2002), subsequently a number of other topologies have been identified that can result in chimera states, including global Böhm et al. (2015) and local Laing (2015); Bera and Ghosh (2016) coupling. There is a large variety in manifestations of chimera states and how they can appear in different systems. If one considers amplitude-phase representation of individual node dynamics, it is possible to distinguish between phase chimeras and amplitude chimeras. The phase chimera is defined as the coexistence of coherent and incoherent regions in the space of phases of different oscillators Bogomolov et al. (2016). In this case, the average phase velocity of different oscillators exhibits a characteristic arc-shape profile, with a pronounced increase or decrease in the average frequency for the incoherent region associated with the chimera state. In contrast, an amplitude chimera appears as a sudden increase or decrease in the average amplitude of oscillations Sethia, Sen, and Johnston (2013); Zakharova, Kapeller, and Schöll (2014); Bogomolov et al. (2016).

In this paper, we consider the emergence and behavior of chimeras in the specific context of epidemic models of multi-strain diseases. A number of effective mathematical frameworks have been developed over the years for the analysis of various aspects of strain interactions Andrews, Halpern, and Purves (1997); Gupta, Ferguson, and Anderson (1998); Gog and Grenfell (2002); Gog and Swinton (2002); Gomes, Medley, and Nokes (2002); Koelle et al. (2006); Koelle, Kamradt, and Pascual (2009), with particular attention being paid to cross-immunity and its effects Calvez, Korobeinikov, and Maini (2005); Adams and Sasaki (2007); Minayev and Ferguson (2008); Cobey and Pascual (2011). Multi-strain epidemic models have been shown to exhibit a wide range of behaviors, including (partially) synchronized dynamics, anti-phase oscillations, as well as chaotic dynamics Minayev and Ferguson (2008); Recker et al. (2009); Kucharski, Andreasen, and Gog (2016). Group-theoretical analysis of multi-strain models has yielded significant inroads to systematic classification of steady states and periodic solutions in terms of their symmetry Blyuss and Kyrychko (2012); Blyuss (2013, 2014); Charles and Baca (2013); Chapman and Mesbahi (2013). Motivated by the recent work on chimeras in locally coupled, delayed oscillators Bera and Ghosh (2016), we explore the dynamics of a multi-strain network, in which coupling between strains quantifies the degree of their cross-immunity, while the dynamics of each individual strain is represented by a compartmental model, with the time delay representing a period of temporary immunity upon recovery from infection.

The remainder of this paper is organized as follows. In the next Section we introduce the model and discuss its basic properties. Section  III contains analytical and numerical bifurcation studies of single-strain dynamics for completely antigenically distinct strains. In Section IV different types of dynamics are investigated in the presence of all-to-all and non-local cross-immunity coupling kernels. The paper concludes in Section V with the discussion of results.

II Model

We consider a multi-strain disease, in which recovery from an infection with any single strain results in a certain period of temporary immunity against subsequent infections with that strain. To analyze the dynamics of such a disease, one can combine an SIRS-type model of temporary immunity proposed by Kyrychko and Blyuss Kyrychko and Blyuss (2005); Blyuss and Kyrychko (2010) with the status-based approach of Gog and Grenfell Gog and Grenfell (2002) for multi-strain diseases, which gives the following model

[TABLE]

where $S_{i}$, $I_{i}$ and $R_{i}$ represent the number of people in the population that are susceptible, infected or recovered from strain $i=1,2,\dots,N$, with $N$ being the total number of disease strains in circulation, $\eta>0$ is a constant birth rate and death rate assumed to be the same for all strains, $\beta_{i}>0$ and $\gamma_{i}>0$ are the transmission rate and the recovery rate of strain $i$, respectively. This model assumes that after recovery, individuals remain the class of recovered from strain $i$ for a period of temporary immunity $\tau>0$, upon which they return to the class of susceptible. For simplicity, we assume that the transmission and recovery rates for all strains are the same, namely, $\beta_{i}=\beta$ and $\gamma_{i}=\gamma$. The factor $0\leq\sigma_{ij}\leq 1$ denotes the reduction in the susceptibility to strain $i$ due to immune response to a previous infection with strain $j$ Gomes, Medley, and Nokes (2002), with zero denoting the complete cross-immunity, that is, the same immunological response between two strains $i$ and $j$, and unity denoting the complete absence of cross-immunity, i.e., absolutely distinct immunological responses against the two strains $i$ and $j$. In this paper, we will consider all-to-all coupling, i.e., $\sigma_{ij}\equiv 1$, as well as two types of non-local coupling kernels that represent more realistic immunological relations between disease strains.

Summation of the left- and right-hand sides of Eqs. (1) yields

[TABLE]

where $N_{i}(t)=S_{i}(t)+I_{i}(t)+R_{i}(t)$ denotes the total population of strain $i$. Since the birth and death rates are equal, the total population for each strain is asymptotically constant Gog and Grenfell (2002); Kyrychko and Blyuss (2005), that is, all $N_{i}$ tend to unity. The observations that $R_{i}(t)=1-S_{i}(t)-I_{i}(t)$ and that $R_{i}(t)$ does not feature in equations for $S_{i}$ and $I_{i}$, suggest that it is sufficient to focus on the dynamics of variables $S_{i}$ and $I_{i}$ only. To reduce the number of free parameters, we rescale time with $(\eta+\gamma)^{-1}$, and introduce a basic reproduction number $r=\beta/(\eta+\gamma)$ and a rescaled mortality rate $\rho=\eta/(\eta+\gamma)$. This gives the following rescaled model

[TABLE]

where the self-coupling term is written out explicitly with $\sigma_{ii}=1$.

III Antigenically distinct strains

Before investigating the collective behavior in the full multi-strain system, it is instructive to consider what happens in the absence of cross-immunity, i.e., when each strain is genetically distinct, so as to cause a completely distinct immunological response to infection, which is represented by $\sigma_{ij}=\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta. In this case, the system (4) decouples into $N$ independent copies, and the dynamics of each individual strain is described by the following system of equations

[TABLE]

This system always has the disease-free steady state $E_{0}=(S_{0},I_{0})=(1,0)$, and it can also possess an endemic steady state

[TABLE]

The endemic equilibrium $E^{*}_{\tau}$ is only biologically feasible if $r>1$, which, in terms of original parameters, corresponds to the transmission rate $\beta$ being larger than the sum of the natural death rate $\eta$ and the recovery rate $\gamma$. In the case of very long immunity period, i.e. for $\tau\rightarrow\infty$, the SIRS model (5) transforms into a standard SIR model with vital dynamics and permanent immunity, and the endemic steady state then reduces to

[TABLE]

Linearization of the system (5) near the disease-free steady state $E_{0}$ gives the characteristic eigenvalues as $\lambda_{1}=-\rho$ and $\lambda_{2}=r-1$, thus implying that the disease-free steady state is stable, provided $r<1$. For the endemic steady state $E^{*}_{\tau}$, the characteristic equation has the form

[TABLE]

which, for a vanishing delay $\tau=0$, always gives stable eigenvalues due to $r>1$. One root of this equation is $\lambda=-\rho<0$, which is stable independently of the time delay. For non-zero immunity period, the endemic steady state can lose its stability in a Hopf bifurcation, giving rise to periodic solutions.

Since for $\tau=0$ the eigenvalues $\lambda$ of the characteristic equation (8) are stable, and $\lambda=0$ is never a solution of this equation, the only possibility how the stability of the endemic steady state can change is if a pair of complex conjugate eigenvalues crosses the imaginary axis for some value of $\tau$. To find this critical time delay, we substitute $\lambda=i\omega$ into Eq. (8) and separate real and imaginary parts, which yields

[TABLE]

Squaring and adding these two equations gives an implicit equation for the Hopf frequency

[TABLE]

which can be readily solved to give

[TABLE]

Alternatively, by dividing the equations (III), we find the critical value of the time delay at which the Hopf bifurcation occurs

[TABLE]

Unfortunately, due to the fact that the steady-state value of the infected fraction $I^{*}_{\tau}$ itself explicitly depends on the time delay $\tau$ as shown in Eq. (6), it does not prove possible to find a closed form expression for the Hopf frequency or the critical time delay.

To get a better understanding of the bifurcations of the endemic fixed point, we perform numerical bifurcation continuation using DDE-Biftool Engelborghs, Luzyanina, and Samaey (2001), choosing $\tau$ as the continuation parameter. Figure 1 illustrates regions of stability and instability of this steady states, together with multiple branches of characteristic eigenvalues. For the chosen parameter values, this figure shows that a single branch escapes the stable region from $\tau_{1}\approx 8.88$ to $\tau_{2}\approx 38.49$, and in this interval of time delays, the endemic steady state is unstable.

Having identified the points at which the endemic equilibrium loses/gains its stability, we now focus on the limit cycle that emerges at these bifurcation points. Figure 2 shows the period and amplitude of the limit cycle depending on the time delay. This figure indicates that at $\tau=\tau_{1}\approx 8.88$, the endemic steady state undergoes a supercritical Hopf bifurcation, giving rise to a stable limit cycle, whereas at $\tau=\tau_{2}\approx 38.49$ it undergoes a subcritical Hopf bifurcation, at which the endemic equilibrium regains its stability, and an unstable limit cycle is born. These two limit cycles coexist for $\tau>\tau_{2}$ until they merge at a point $\tau_{3}\approx 43.97$ and annihilate in a saddle-node bifurcation of limit cycles.

IV Multi-strain dynamics

As a next step, we consider the network of $N$ coupled strains (4), where in the absence of coupling the dynamics of each strain is described by a delayed SIR model (5). Before proceeding with numerical simulations, it is worth noting that for any form of the coupling $\sigma_{ij}$, the system (4) admits a one-strain solution with $I_{i}(t)\neq 0$ and $I_{j}(t)=0$ for $j\neq i$ that defines an invariant manifold (cf. Blyuss & Gupta Blyuss and Gupta (2009) for a similar type of behavior in a $\mathbb{D}_{4}$-symmetric model of antigenic variation), and whose behavior is described by the following system

[TABLE]

Effectively, the system decouples into the single-strain dynamics (5) for strain $i$, which then drives the evolution of $S_{j}$ variables, while all $I_{j}$ remain zero. The equivalent one-strain endemic steady state is given by

[TABLE]

In the case of all-to-all coupling with $\sigma_{ij}=1$, the system (4) possesses a $\mathbb{Z}_{N}$ symmetry, hence it has $N$ identical one-strain steady states given by Eq. (14) for any $i=1,...,N$. Furthermore, for such coupling the system (13) reduces to just strain $i$ with the dynamics given by Eq. (5), and all other strains, whose dynamics is exactly the same and is fully driven by the strain $i$. Techniques of equivariant bifurcation theory can be used to systematically characterize various steady states and periodic solutions in terms of their symmetry Blyuss and Kyrychko (2012); Blyuss (2013); Blyuss and Gupta (2009); Blyuss (2014).

Besides one-strain steady states, the system (4) also has a fully symmetric endemic steady state

[TABLE]

where

[TABLE]

with $\sigma_{c}=\sum\limits_{i\neq j}\sigma_{ij}$.

For each type of coupling, we have used the dde23 solverShampine and Thompson (2001) to numerically integrate the system (4) with the initial conditions taken as follows: $S_{i}$ are uniformly distributed random numbers between [math] and $1$ independent for each strain, and random $I_{i}\in[0,1-S_{i}]$ being constant in $t\in[\tau,0)$. We investigate possible dynamical behavior for three different types of coupling between strains: the all-to-all coupling, a Gaussian kernel based on the model of Gog and Grenfell Gog and Grenfell (2002), and a functional cosine kernel suggested by Gomes et al. Gomes, Medley, and Nokes (2002) Since the last two kernels are non-local, in principle, one can expect to observe chimera states in such multi-strain systems Kuramoto and Battogtokh (2002); Abrams and Strogatz (2004), and below we investigate the appearance of such states and transitions between them and other dynamical regimes.

IV.1 All-to-all coupling

In the case of global all-to-all coupling $\sigma_{ij}=1$, the same amount of cross-immunity is present between all interacting strains, which biologically means that every strain is related to all other strains in the same way.

Figure 3 shows the dynamics of system (4) with $N=60$ strains for an all-to-all coupling and time delay $\tau=25$, for which a stable limit cycle is observed in the single-strain dynamics. The time series, as well as the snapshot and the space-time plot, indicate that in this case all nodes become synchronized, except for one strain ($i=40$ here), as shown in Fig. 3. The latter strain exhibits large-amplitude oscillations in both $S$ and $I$ variables, which then drive smaller amplitude oscillations in the $S$ variable for all other strains. As discussed earlier, the dynamics of such a solitary state can be effectively described by a reduced two-strain model: one delayed model (5) for the solitary strain, and one for all other synchronized strains, as given in Eq. (13).

One should note that due to the above-mentioned $\mathbb{Z}_{N}$ symmetry of the system, the fact that the system has settled on the strain $i=40$ being the main driving strain is completely random and is purely determined by the initial conditions, as for the same parameter values, any of the other solitary states is equally possible. The other observation is that since the system starts with random and independent initial conditions for all strains, the fact that eventually it settles on a solitary state suggests that a one-strain invariant manifold described by Eq. (13) is stable. Moreover, since this corresponds to a situation where in the absence of coupling all individual strains have the dynamics of a stable limit cycle, effectively the coupling appears to suppress these oscillations in a manner similar to symmetry-breaking oscillations death that has been recently studied in time-delayed systems Zakharova et al. (2013).

IV.2 Non-local kernels

By analogy with non-local coupling kernels for which chimera states have been observed in various systems of coupled oscillators Kuramoto and Battogtokh (2002); Abrams and Strogatz (2004); Omel’chenko, Wolfrum, and Maistrenko (2010), we focus our attention on two kernels that represent the biologically realistic scenario where the more related strains are, the higher is the level of cross-immunity between them Gog and Grenfell (2002); Gomes, Medley, and Nokes (2002). The first example is a slightly modified Gaussian kernel introduced in Gog & Grenfell Gog and Grenfell (2002)

[TABLE]

where $d$ is the characteristic length associated with cross-immunity, and the distance between strains $i$ and $j$ is measured as the smallest difference on the interval $[1,N]$ with periodic boundary conditions. Strains that are genetically close to each other have a higher value of cross-immunity $1-\sigma_{ij}$, leading to a decrease in the inflow of the infected population for the strain at hand. This effect is a combination of the reduced susceptibility and reduced infectivity due to various immunological interactions between strains Gog and Grenfell (2002); Ballesteros, Vergu, and Cazelles (2009). Figure 4 illustrates the shape of the kernel $\sigma_{ij}$ for different characteristic lengths $d$.

A similar reasoning, but with a different biological rationale, is used in the model of Gomes et al. Gomes, Medley, and Nokes (2002) who considered strains as being distributed on the unit circle with positions $z_{i}=(2i-1)/2N$ along the circle, with the kernel being given by

[TABLE]

with

[TABLE]

The profile of $\sigma_{ij}$ depending on the distance between strains is illustrated in Fig. 5 for different values of parameter $p$.

In the coupling kernel (17) there are two different parameters that characterize the strain space. Firstly, there is $\sigma$ ($0\leq\sigma\leq 1$), which plays the role of the bound on the range of the strain diversity. Secondly, there is $p$ which represents antigenic differences between strains for the given genetic range. Gomes et al. Gomes, Medley, and Nokes (2002) focused on the specific values of $p=-2$, [math] and $2$, but one can prove that parameter $p$ must lie in the range $p\in[-2,4]$ to ensure $\sigma(z)$ has a single maximum at $z=0.5$ and two minima at $z=0$ and $z=1$, which biologically means that the strain most genetically different from the current strain experiences the smallest amount of cross-immunity.

IV.3 Emergent dynamical scenarios

Below we present and discuss different patterns observed in the case of a non-local Gaussian coupling kernel (16). Figures 6-11 illustrate a modulated-amplitude profile, a solitary state, a traveling wave, (multi-headed) amplitude chimeras, and a transition state, respectively. To get a better insight into the dynamics, in each case the actual time series is plotted for all $N$ strains, accompanied by a snapshot at a fixed moment in time, a space-time plot, as well as plots of the average amplitude of oscillations for both dynamical variables. The amplitude is computed as the difference between maximum and minimum values of the respective variable for each strain.

Figure 6 shows a regime where all strains oscillate with the same frequency and without phase shift, but with different amplitudes, as is clear from the space-time plots and the plots of the amplitude. Since many of the $I$ variables stay equal to zero in a manner similar to all-to-all coupling, while the frequency of oscillations is the same for $S$ variable for all oscillators, for $I$ variables it gets adjusted to the frequency of $S$ variables for those strains that do exhibit oscillations. The highest amplitude of oscillations occurs in the middle of modulated profile, suggesting the potential for amplitude rather than phase chimeras, but since the snapshot of the modulate profile is smooth, this state cannot be interpreted as a proper chimera state Omelchenko et al. (2011, 2012). From epidemiological perspective, this is an interesting state in that all non-zero strains follow synchronous oscillations, namely, they appear and disappear at the same time. On the other hand, infected fractions have substantially different magnitudes, which means that immunological interactions between strains results in some of them always being more dominant (i.e. having a significantly larger amplitude), whereas other strains are more suppressed, and this relation between different strains is repeated with every oscillation.

An exemplary case, where only a few strains exhibit oscillations of considerable amplitude, is shown in Fig. 7 for a larger value of time delay $\tau=42$ and a smaller characteristic length $d=4$.

For small coupling strength $\sigma$ and narrow, that is, local, coupling kernels, we find a traveling wave pattern shown in Fig. 8. This observation is important from a biological point of view, as it illustrates a regime of sequential strain dominance, which is often observed in epidemiological data Recker et al. (2009); Blyuss (2014).

In Fig. 9 we illustrate the dynamical regime of an amplitude chimera. The different strains still oscillate with the same frequency, but in contrast to the previous pattern, the snapshots do not exhibit a smooth profile anymore but rather are represented by two different regions: coherent and incoherent. It is worth noting that whilst the coupling is still non-local, the amplitude chimera is observed even when the characteristic length of the coupling is quite small $(d=2)$. For the same parameter values but a smaller time delay, the system can also exhibit a multi-headed amplitude chimera, characterized by several coherent and incoherent regions with almost no variation in terms of frequency, but showing the amplitude profile typical for chimera states. An example of such state is shown in Fig. 10.

A pattern of transition between an amplitude chimera and a modulated profile is demonstrated in Fig. 11. Whilst there is an incoherent region in the middle of strain domain, the edges of the chimera have a smooth profile similar to that of the modulated profile, indicating that being a transition, this regime features the characteristics of both the chimera and the modulated profile. Similar to the amplitude chimera, the largest amplitude of oscillations for the transition state occurs in the incoherent regime. It should be noted that transition states can be found for a whole range of parameter values between modulated profile and amplitude chimeras, making them closer in terms of dynamics to either of those states.

Figure 12 provides a summary of different dynamical states that can be observed in the system (4) depending on the time delay $\tau$ and the cross-immunity length scale $d$. Larger values of the cross-immunity length scale, i.e. broader coupling kernels, are associated with modulated amplitude profiles, while, surprisingly, chimera states (single- and multi-headed) are found for narrower, i.e. more local, coupling kernels. Solitary states in which infections with only a single strain are present, can occur for any lengths of cross-immunity $d$, provided the time delay $\tau$ is sufficiently large.

We have also performed extensive simulations for the case of cosine kernel (17), and a summary of results is shown in Fig. 13. Unlike the Gaussian kernel, in this case only modulated profiles, solitary states, and transition states are observed, while traveling waves and amplitude chimeras were never found. The most likely explanation for this lies in the fact that amplitude chimeras are associated with quite narrow Gaussian kernel (as described by small values of $d$), whereas for the biologically feasible values of parameter $p\in[-2,4]$, the distribution (17) is quite broad. In fact, Figure 5 suggests that the narrowest width of the cosine distribution corresponds to $p=-2$, which, for a system of $N=60$ strains is equivalent to $d=12$, and for large values of $p$, the coupling is very broad, making it more similar to the situation described by an all-to-all coupling. As a result, the dynamics is dominated by modulated amplitude profiles for smaller durations of temporary immunity, and by solitary states with single-strain dynamics for larger values of the time delay.

V Discussion

In this paper we have studied an important question about the range of dynamical behaviors that can be exhibited by multi-strain epidemic models with temporary immunity and various types of cross-immunity. Whilst the time delay associated with temporary immunity provides a simple mechanism supporting stable oscillations in the susceptible and infected populations for individual disease strains, this dynamics undergoes major changes under the influence of long-range coupling. Under assumption of all-to-all coupling, the system settles on the dynamical regime of solitary states, or single-strain oscillations, where infections with only one strain are present, while all other strains remain equal to zero. Interestingly, the system approaches such a state for random initial conditions, suggesting that this state is, in fact, a stable invariant manifold of the model, which dynamically represents a symmetry-breaking suppression of oscillations. The complete symmetry between strains means that the surviving strain is determined purely by the initial conditions, and for the same parameter values, all other strains are equally possible.

For the case of Gaussian cross-immunity kernel, the model exhibits a wide range of dynamical scenarios that include solitary states, traveling waves, and, most interestingly, single- and multi-headed amplitude chimeras, characterized by some groups of strains oscillating coherently, while others are performing incoherent oscillations. Whilst the cosine kernel is also non-local, by virtue of being very broad, the range of different behaviors for this kernel is smaller and is more reminiscent of the case of all-to-all coupling. The fact that chimera states were observed only for sufficiently narrow cross-immunity kernels suggests that in epidemiological data these types of solutions would only be observed in the cases where individual strains or serotypes elicit cross-reactive immune responses against very genetically similar strains. For multi-strain diseases with a wide antigenic repertoire, chimera states could be interpreted as dynamical regimes where a number of closely immunologically related strains appear to have similar dynamics and show up concurrently, while other strains have irregular and unsynchronized oscillations. Understanding parameter regimes that result in chimera states can provide useful insights for design and deployment of multi-valent vaccines.

Acknowledgements.

PH and JB acknowledge support by Deutsche Forschungsgemeinschaft in the framework of the Collaborative Research Center SFB 910. YK and KB acknowledge the hospitality of the Institute of Theoretical Physics, TU Berlin, where part of this work was completed.