Dynamics of epidemic diseases without guaranteed immunity

The pandemic of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) suggests a novel type of disease spread dynamics. We here study the case where infected agents recover and only develop immunity if they are continuously infected for some time τ. For large τ, the disease model is described by a statistical field theory. Hence, the phases of the underlying field theory characterise the disease dynamics: (i) a pandemic phase and (ii) a response regime. The statistical field theory provides an upper bound of the peak rate of infected agents. An effective control strategy needs to aim to keep the disease in the response regime (no ‘second’ wave). The model is tested at the quantitative level using an idealised disease network. The model excellently describes the epidemic spread of the SARS-CoV-2 outbreak in the city of Wuhan, China. We find that only 30% of the recovered agents have developed immunity.


Introduction
The rapid spread of a disease across a particular region or regions (epidemic) or the global outbreak of a disease (pandemic), see Porta [17], can have a detrimental effect on health systems, on local and global economies including the financial markets and the socioeconomic interactions, ranging from the city to the international level. Measures to reduce the pandemic spread include curtailing interactions between infected and uninfected parts of the population, reducing infectiousness or the susceptibility of members of the public, see e.g., Ferguson et al. [5]. The two major strategies governments use to handle an outbreak are to slow down an outbreak (mitigation) or to interrupt the disease spread (suppression). Since each of those interventions bears itself significant risks for the societal and economic well-being, it is crucial to understand the effectiveness of these strategies (or any hybrid of them).
Mathematical methods provide essential input for governmental decision making that aims at controlling the outbreak. Among those are statistical methods, Unkel et al. [22], Becker and Britton [2], deterministic state-space models, Brauer et al. [3] with its prototype developed by Kermack et al. [12], and a variety of complex network models, e.g., Hwang et al. [10], Shirley and Rushton [19]. The different mathematical approaches have different objectives: A significant application of the statistical methods frequently aims at the early detection of disease outbreaks as described by Unkel et al. [22], while modelling either tries to develop a model as realistic as possible for a given outbreak or to design a simplistic model, which, however, reveals some universal truth about the outbreak dynamics.
In the simplest version, the so-called compartmental models (see Kermack et al. [12], Hethcote [9]) consider the fraction of the population which is either susceptible (S), infected (I) or removed (R) from the disease network. Coupled differential equations capture the dynamics of the disease that determine the time dependence of S, I and R. Extensions add more compartments to the Susceptible Infected Removed (SIR) model such as (E) exposed. For example, such an Susceptible Exposed Infected Removed (SEIR) model was used by Lekone and Finkenstädt [15] for a description of the Ebola outbreak in the Democratic Republic of Congo in 1995. Compartmental models have been applied to describe the recent SARS-CoV-2 outbreaks. Selected publications are: Giordano et al. [7], Krishna and Prakash [13], Tagliazucchi et al. [21], Lin et al. [16], Anastassopoulou et al.  [4] in order to capture stochastically unknown influences, such as changing behaviours. Such models were recently used to analyse the COVID-19 outbreak in Wuhan by Kucharski et al. [14]. A novel extended epidemic SEIR model, taking into account by a socio-political classification of different interventions, was proposed by Proverbio et al. [18] for assessing the value of several suppression approaches.
Compartmental models address global quantities such as the fraction of susceptible individuals and assume that heuristic rate equations can describe the disease dynamics. In cases of a strongly inhomogeneous (social) network, e.g. taking into account different population densities, the above assumption seems not always be justified. In these cases, spatial disease spread patterns can be described by a stochastic network model with Monte-Carlo simulations a common choice for the simulation.
In this paper, we consider a disease dynamics for which the duration (severity) of the illness depends on the amount of exposure. Using an elementary (social) network, we are looking for universal mechanisms describing a pandemic spread. We will reveal a connection to statistical field theory, enabling us to characterise an outbreak with the tools of critical phenomena. We will discuss the impact of the findings on policies to curb an outbreak and will draw conclusions from the COVID-19 outbreak in Wuhan, Hubei province, China.

Model basics
Each individual interacts with four 'neighbours' of the social network. The disease spread is described as a stochastic process. At each time step (say 'day'), the probability that an individual gets infected (or recovers) depends on the status of the neighbours in the social network. Here, we only study the simple case of a homogeneous network with four neighbours for each site. We also consider periodic boundary conditions to minimise edge effects.

Immunity
We study two closely related scenarios.
(i) There is no immunity. Every individual can be reinfected and can recover only to be susceptible again (Susceptible Infected Susceptible (SIS) model). (ii) Individuals can be reinfected and recover. Only if individuals stay infected for τ consecutive days, they are considered immune. In case (ii), the sites of immune individuals are removed from the disease network.

Disease dynamics
If x is a site of the disease network, at every time step the state u x ∈ {0, 1} is randomly chosen with probability where xy is an elementary link on the lattice joining sites x and y and, hence, n is the number of infected neighbours, and N x = 1 + exp{4βn x + 2h} is the normalisation. The parameter h is linked to the probability to contract the disease from outside the network.
In fact, if no-one of the network is infected (n x = 0, ∀x), the probability p that any individual contracts the diseases, is connected to h by The parameter β describes the contagiousness of the disease. The probability that any individual gets infected (u x = 1) monotonically increases with 4βn x + 2h. The parameter β hence describes how sensitive this probability depends on exposure, i.e., the number n x of infected neighbours. If the lattice contains N individuals (i.e., sites), one time step is said to be completed if we have considered N randomly chosen sites for the update.

Peak infection rate
Scenario (ii) shows the typical time evolution of an epidemic with the infections rate approaching zero for large times due to agents recovering and an increasing number being immune. By contrast, scenario (i) has an asymptotic state independent from the initial state and described by statistical field theory. After the change of variable z x = 2u x -1, the asymptotic state is described by the partition function of the Ising model, Ising [11], Friedli and Velenik [6]: For a vanishing external field H, the model shows a critical behaviour with a phase transition at β = β c = ln(1 + √ 2)/2 ≈ 0.44. In the ordered phase for β > β c , a small seed probability p > 0 triggers an infection rate close to 100 % of the population. The model is in the 'pandemic' phase. For β < β c , the model is in the 'response' phase, i.e., the infection rate is in response to the seed probability p, but no outbreak occurs. The asymptotic infection rate can be calculated using Markov Chain Monte-Carlo (MCMC) methods. Starting, e.g., with no infected agents (u x = 0 or z x = -1), each time step (see Sect. 2) creates one sample of the disease spread. Each disease pattern only depends on that of the previous day, and the sequence of the days form a Markov chain. After some time, which is called 'thermalisation time' in statistical physics, the daily infection rate starts to fluctuate around the average, i.e, the asymptotic rate. The asymptotic rate is independent of the simulation details if certain MCMC conditions are satisfied. Among those, ergodicity is easily violated in the pandemic phase for the so-called local update algorithms, most prominently the Metroplis-Hastings one, Hastings [8]. Here, we used the state-of-the-art Swendsen-Wang cluster algorithm, which performs well across both phases, see Swendsen and Wang [20]. Our numerical findings are summarised in Fig. 1, left panel. Curve (2) clearly separates both phases-the pandemic phase and the response regime.

Immunity
Let us now study scenario (ii), where individuals can develop immunity if they are infected for τ consecutive days. For τ > t th , the peak infection rate is that of the asymptotic state of the corresponding model (i) and, hence, inherits the classification 'pandemic' or 'response' phase. This is illustrated in Fig. 1, right panel, for the pandemic phase for several values of τ . Figure 2 illustrates the vastly different behaviour of the disease spread in the pandemic phase (β = 0.41, p = 5%) and in the response regime (β = 0.38, p = 4%). Results are for a N = 100 × 100 network and τ = 11. Note that the curve for 'infected + immune' ('triangle' symbol) in the pandemic phase is not monotonically increasing with time since the infected individuals can return to 'susceptible' state, i.e., not every infected individual becomes immune. Note that in the response regime ('circle' and 'square' symbol), the 'pan- demic' peak is absent altogether. However, on the downside, the so-called 'herd immunity' slowly develops over time.

Comparison with data
We stress that the model assumption of a homogeneous (social) network with 'four neighbours' is unrealistic. A study of an heterogeneous disease netork is work in porgress. The knowledge of the underlying disease network is essential to make quantitative predictions for e.g. the critical value β c of the contagiousness. Here, we adopt a different approach: we assume that qualitative time evolution of bulk quantities such as the fraction of infected individuals is within the grasp of model scenario (ii) and use those as fit functions to determine the model parameters such as β, p and τ by comparison with actual data.
For this study, we used data from the COVID-19 outbreak in 2020 in the city of Wuhan, Hubei province in China, Yu [24] (accessed April 16, 2020). The data of the number of infected individuals shows a jump at day 73 (on the arbitrary time scale) by 40%, which is due to a change in reporting. We assume that the same 'under-reporting' has occurred in the days before. Guided by the fact that the probability distribution (the infected rate) is a continuous function, we have corrected the data by multiplying the number of infected (and infected + recovered) by a factor 1.4 for times t ≤ 73. Let D(t, τ , β, p) be the fraction of the population of infected individuals as a function of time t and depending on the parameters τ (time to develop immunity), β (contagiousness) and p (seed probability) to get infected. We have calculated D(t, τ , β, p) using a N = 250 × 250 lattice. We checked that the result is independent of the lattice size in the percentage range for the parameters relevant in this study. If D wuhan (t) quantifies the measured values for the number of infected in the Wuhan outbreak, we want to approximate these data, i.e., with a suitable choice of the parameter N pop , t s , β and p. Since the offset of the time axis in the Wuhan data is arbitrary, we have chosen the shift t s such that the peaks of simulated data and measures data coincide. All other parameters are treated as fit parameters. Those parameters have been obtained by fittng the model to the infected data only. Altogether, Figure 3 Comparison of scenario (ii) results with actual data from the COVID-19 outbreak in Wuhan we find a good agreement with the data for (see Fig. 3): The results for 'recovered + infected' and 'immune' are then model predictions. The former can be compared with actual data to gauge the viability of the model. The model data overshoot the data in the early days of the epidemic spread, which could be related to underreporting due to limited testing capabilities. It is interesting to observe that the curve of the infection rate is asymmetric: the slope of the raise at the beginning is larger than the slope of the decline after the maximum. Also, the number of infected seem to level off at a non-zero value. In the present model, this is explained as follows: with more agents being immune, it is harder for susceptible agents to be continuously infected for time greater or equal τ and, thus, to develop immunity. We also find that only about 30% of the infected (and recovered) develop an immunity.

Conclusions and interpretations
A new type of stochastic disease model is proposed: agents can recover from an infection and are susceptible again. They only develop immunity if their infection lasts longer than a characteristic time τ . For τ → ∞, the infection rate is described by statistical field theory. For finite τ , the infection rate of the field theory provides an upper bound of the infection rate of the dynamical model. This opens up the possibility to characterise the disease dynamics in the light of critical phenomena of the underlying field theory: a pandemic spread corresponds to the ordered phase of the field theory, and the critical value for the contagiousness is that for the phase transition. The disease is in controllable response mode if the corresponding field theory is in the disordered phase.
Quantitative results, reported here, are derived with an unrealistic homogeneous disease network for which each agent interacts with four neighbours. Nevertheless, we find that the COVID-19 data of the Wuhan outbreak are well represented. For this case, we find that only 30% of infected develop an immunity.
The heavy tail of the decline of the number of infected, which levels off at non-zero values, is an inherent feature of the model and can be traced back to the fact that agents can be reinfected. In a network with a sizeable portion of immune agents, it is increasingly challenging to develop immunity. If these model assumptions were underpinned by medical investigations, achieving 'herd immunity' would be difficult. This should influence the decision to what extent efforts focus on developing a cure or a vaccine.