Plateaus, rebounds and the effects of individual behaviours in epidemics

Plateaus and rebounds of various epidemiological indicators are widely reported in Covid-19 pandemics studies but have not been explained so far. Here, we address this problem and explain the appearance of these patterns. We start with an empirical study of an original dataset obtained from highly precise measurements of SARS-CoV-2 concentration in wastewater over nine months in several treatment plants around the Thau lagoon in France. Among various features, we observe that the concentration displays plateaus at different dates in various locations but at the same level. In order to understand these facts, we introduce a new mathematical model that takes into account the heterogeneity and the natural variability of individual behaviours. Our model shows that the distribution of risky behaviours appears as the key ingredient for understanding the observed temporal patterns of epidemics.


Supplementary information
The classical SIR model with time dependent transmission rate The classical SIR compartmental model 1 uses a dynamics governed by mass-action laws. It assumes that individuals are homogeneously mixed and that every individual is equally likely to interact with every other individual. It reads: (1)

Single location growth and behavioural changes.
Changes of behaviour driven by phenomena such as awareness of fatalities or fatigue with respect to mobility restrictions can be modelled by a time modulation of the infection rate β .
To start with, we analyse this problem on synthetic data. Namely, we consider an exact dynamics given by the SIR model with parameter changes on two given dates. We now take part of the resulting points as given observations. The goal is to carry an optimisation procedure where we try to identify the dates and magnitudes of these changes i.e. the different values of β that come into play. Later on we will consider the Thau lagoon data. Figure S1 represents a piece-wise constant modulation in time with reduction of social interactions down to 70% starting Day t = 25 followed by an increase to 50% on Day t = 50. The effect on such modulation leads to the interruption of the growth phase followed by a lower decrease for later times. We developed a toolbox based upon PYGMO 2 particle swarm optimisation algorithms in order to solve the type of inverse problem we encounter here. We consider a time sampled version of infection rates in the total population. We allow a given number of transition times associated with behaviour changes, in between which the coefficient β is constant. The problem then is to determine in an optimal way the SIR model parameters, initial conditions, the transition times and the values of the β 's in between. Using a least square objective function, we obtained a satisfactory convergence towards the dynamics of the initial model and associated modulation function. We developed a similar approach with noise added to sampled synthetic data, with satisfactory convergence both in terms of dynamics of infection rate ( Figure S5) and modulation functions (levels and time changes in Figure S6). We applied the preceding algorithmic approach on real data of the four WWTP measurement campaigns of the Thau lagoon. In order to keep the model parsimonious, we used the same time modulation function of β for all the four cities together. The transition dates and levels have to be determined. Figure S7. Change of behaviour in logarithmic scale in the Thau lagoon (initial day t = 0 being associated with May 12 th , 2020). The red curve represents the model output (rate of infected population) and the blue dots correspond to measurements. Zones I 1 , I 2 , I 3 and I 4 respectively correspond to Sète, Mèze, Frontignan and Pradel-Marseillan.
The above results show that the optimal capture of the transition from exponential growth to plateau type dynamics occurs 4/15 on August 10 th , 2020. The resulting piecewise constant modulation function in time is represented in Figure S8 below. Even though the plateau dynamics seems to be correctly represented in Figure S7, the main shortcoming of this is that it does not really explain the observations of the Thau lagoon data. Indeed, the model does not explain why the four zones reach the same infection rate level associated with the plateau regime and why there is a two weeks delay. Furthermore, it does not yield a shoulder effect, as seen in the data.
One may use spatial diffusion to describe the delay between the zones (Frontignan, Mèze) on the one hand and (Sète, Pradel-Marseillan) on the other hand. But it also fails to explain why there are two parallel curves. We plan to come back to spatial diffusion in future work.

Derivation of the Model
This section is concerned with the formal derivation of the additional diffusion term in the dynamics of susceptible individuals. As previously emphasised, this diffusion term stems from the combination of the heterogeneity of behaviours and their variability in time. A natural case to consider here is to assume that the shuffling of behaviours happens according to Brownian motion. Namely, individual risk traits move according to the process where σ > 0 is the possibly a-dependent diffusion, and W t is a Brownian motion in (0, 1) with reflection conditions at the end-points of the interval. By the Fokker-Planck equation, in the absence of epidemic, an initial distribution of population S(0, a) gives rise to S(t, a) = E[S(0, a t )] governed by To consider the dynamics of the epidemic, we are thus led to the following more general version: The first equation involves drift and other additional terms. The numerical and theoretical analysis of the case where σ depends on a, and the corresponding extensions of the Fokker-Planck equation with drift terms, seem relevant from a modelling viewpoint.
In the simple case where σ 2 (a) = d where d is a constant, one ends up with the first equation of our system (see System (4) below).
One could also envision a multi-dimensional variable a = (a 1 , . . . , a p ) to encompass various behavioural characteristics that could possibly be related to sociological observations. In this case, one would be led to a higher dimensional partial differential equation. where the term ∂ 2 aa S is replaced by the Laplace operator ∆ a S as we did in the general model presented in Methods section. Such developments would be useful to improve the model and its applicability.

Simple properties of Model (4)
For the convenience of the reader, we recall the model introduced in the main paper The case d = 0 (heterogeneous but without variability) writes: We list below some elementary mathematical properties of the model. First, we note that this model contains the traditional SIR model when the initial profile of susceptible a → S 0 (a) is uniform in a. Indeed, it is straightforward to see that then, S(t, a) also does not depend on a.
The dynamics of total infectious individuals is governed by so that the equivalent of the so called "effective R" coefficient in traditional SIR model, namely β S(t)/(γN), is replaced by Unlike traditional SIR model, in which the evolution of infectious t → I(t) begins by an exponentially growing phase, has a unique maximum and tends to zero for large time, Model (4) may exhibit more sophisticated behaviours such as rebounds in I with multiple local maxima as well as plateaus. This property can be intuitively understood by analysing the evolution of the growth factor in (6): In the absence of diffusion, d = 0, this growth factor is non increasing in time. For heterogeneous profiles and non zero diffusion, the first term of the right hand side is always negative whereas the second term may be positive if for instance a → β (a) is increasing and a → S(t, a) is a decreasing function. Below we show that it is always the case under some conditions. Thus, depending on the amplitude of the above two terms, the growth factor may have increasing and decreasing phases. It naturally leads one to consider initial conditions characterised by the property S 0 (a)β (a) < 0 such as It means that individuals are sorted by increasing infection rate, according to the variable a, and that the vast majority of individuals have low infectiousness whereas few individual close to a = 1 have extremely risky behaviour. To some extent, social contact surveys 3-5 justify an initial distribution of susceptible individuals like (10). However, we are not aware of a rigorous justification of the expression (9) of the sharply increasing profile in β .
Let us now show that the evolution preserves the property of being decreasing in a for the profiles of the susceptible population.
Proposition 1. If S 0 is a decreasing function of a, then S(t, ·) is also decreasing in a for all time t. To prove this property, we simply note that the derivative v(t, a) := ∂ S(t, a)/∂ a of S with respect to a satisfies the following equation Since the right hand side of the first equation of (11) is negative and the initial condition v(0, a) ≤ 0, the parabolic maximum principle then shows that v(t, a) ≤ 0 for all further time.

Further mathematical properties of Model (4)
Let us first notice that we can describe the dynamics of susceptible individuals in (4) in terms of its probability density where we denoteψ(t) = 1 0 ψ(t, a) f (t, a) da. This equation shows the effect of infections on the distribution of susceptible population as a function of a. Indeed, it shows that, for instance in the absence of diffusion (d = 0), the proportion of the population of high a goes down as a result of infection affecting it more in relative terms than the remaining of the population. The presence of diffusion (d > 0) mitigates this effect by fuelling as it were the epidemic with individuals who had initial lower risk trait.
The equivalent of (8) for the average transmission rateβ is In other words, in the absence of behavioural variability (d = 0), the average of the transmission rate decreases under the effect of its variance, the latter being directly linked to behavioural heterogeneity. Thus a higher heterogeneity promotes a faster decay of the average transmission rate. This is a remarkable property of the system in absence of diffusion or with small diffusion.
In the presence of diffusion d > 0, this effect is in competition with the second term of the right hand side of (13) which may be positive if for instance the profile a → β (a) is increasing while the distribution f decreases along variable a.
Our next result concerns the effect of heterogeneity on herd immunity. We analyze it in the absence of diffusion.
) be the solution of (5) (that is when d = 0) with initial conditions S(0, a) =S 0 f 0 (a), I(0) =Ĩ 0 > 0 and R(0) = 0 (for some positive constantsS 0 andĨ 0 ). Then, t → S(t, ·) tends as t → +∞ to a limit profile a → S ∞ (a) in such a way that whereS ∞ > 0 is the limit value as t goes to +∞ of susceptible individualsS(t) associated to the homogeneous SIR model (1) with initial conditions (S 0 ,Ĩ 0 , 0) and transmission rateβ 0 = 1 0 β (a) f 0 (a) da. Thus, this theorem asserts that heterogeneity lowers the herd immunity rate needed to stop an epidemic. Note that Several works study the asymptotic states of similar models without diffusion, in a discrete class framework. We mention in particular Magal et al. 6 , Dolbeault and Turinici 7 and Almeida et al. 8 .

7/15
Sketch of proof: The existence of the limit profile a → S ∞ (a) follows from the monotonicity of t → S(t, a) for fixed a ∈ (0, 1). Since t → I(t) can be easily proved to vanish for large time, the recovered compartment R(t) also tends to a limit R ∞ . Using R as a time variable leads to: Equation (14) expressesĨ 0 as an increasing function F β of R ∞ , which is therefore uniquely defined. Application of Jensen's convexity inequality then leads to We deduce that F β (x) ≥ Fβ 0 (x) for all x > 0, where Fβ 0 is associated with the homogeneous SIR model (1) with initial conditions (S 0 ,Ĩ 0 , 0) and homogeneous parameterβ 0 , which allows to conclude.
An interesting property of System (5) is the conservation of the following entropy-like quantity The property V (t) = V (0) can be deduced from (5). Next, we consider more specifically the role of diffusion d > 0. First, we note that in the presence of diffusion, the above entropy V d based on solutions (S d (t, a), I d (t), R d (t)) of (4) is non-decreasing in a similar manner as in the framework of viscous perturbations of scalar conservation laws: The following result emphasises the differences between the case with diffusion (d > 0) and the case without (d = 0).
Proposition 2. Assume d > 0 and consider solutions (S d (t, a), I d (t), R d (t)) of (4). The following properties hold: (i) t → S d (t, a) is not necessarily decreasing in time for a ∈ (0, 1) given.
(ii) If R 0 = 1 0 S d 0 (a)β (a) da/(Nγ) > 1, then there exists t 0 > 0 such that t → I d (t) is increasing for t ∈ (0,t 0 ). We will give the proof of this result in a forthcoming work.
We now compare the combined impact of heterogeneity and variability compared with the case of homogeneous populations.
Theorem 2. For given d > 0, let (S d (t, a), I d (t), R d (t)) be solution of (4) with initial conditions (S d 0 (a) = S 0 (a), I d (0) = I 0 , R d (0) = 0). The following results holds (i) S d (t, ·) converges as t goes to +∞ to a positive constant S d ∞ independent of a.

8/15
We already know from Theorem 1 above that heterogeneity is beneficial in terms of herd immunity. From this point of view, the effect of variability is a priori not so clear because the constant shuffling of population keeps fuelling as it were the epidemic. However, part (ii) in the previous result shows that heterogeneity and variability still have a positive effect on herd immunity (that is S d ∞ >S ∞ ). The convergence of solutions to states that do no depend on a is natural because of the role of diffusion. Indeed, the diffusion term takes over once the I(t) term becomes small and thus the first equation describes a diffusion. But this only happens in the long term and the relevance here is rather at intermediate times.
Sketch of proof: Assuming that the initial profile S 0 and β are regular enough in the a variable, integration by parts of the S equation of (4) multiplied by ∂ t S leads to: Hence, so that we have It follows that ∂ a S d (t, ·) converges to 0 in L 2 (0, 1): if it does not hold, there exists (t n ) n∈N → +∞ and α > 0 such that ∂ a S(t n , ·) 2 L 2 (0,1) ≥ α. On the other hand, there exists T such that From (21) generalised between two positive times, for all t > T , there exists n such that t n > t and which contradicts the integrability of ∂ a S in L 2 (R + × (0, 1)). Next, from (21), we deduce that S d (t, ·) is bounded in H 1 (0, 1) uniformly in time, so that it is compact in C 0,1/2 (0, 1): a sequence (t n ) n∈N → +∞ exists such that S d (t n , ·) converges in 9/15 C 0,1/2 (0, 1) to some S d ∞ ∈ H 1 (0, 1) as n goes to +∞. Since ∂ a S d ∞ = 0, S d ∞ is a constant. In order to prove the uniqueness of S d ∞ , we observe that for all (t,t ) ∈ R + 2 and use the fact that ∂ a S ∈ L 2 (R + × (0, 1)), S/N ≤ 1, and I ∈ L 1 (R + ).
In order to prove (ii), using the fact that dR d (t)/dt = γI d (t) ≥ and R d (t) ≤ N, we deduce that R d (t) converges to some R d ∞ > 0. The conservation of total population 1 0 S d (t, a) da + I d (t) + R d (t) = 1 then yields the convergence of I d (t) to 0 as t → +∞.
Integrating (17) between 0 and t, and letting t go to +∞ leads to which yields the estimate S d ∞ >S ∞ . The detailed proofs and further developments are postponed to a forthcoming work.

Some simulations with rebounds and plateaus
Some simulations show that the model adequately captures dynamical features such as epidemiological rebounds and plateaus. The initial conditions considered here are (10) where κ = 0.44, S 0 = 90 and S 1 = 119, 000 (total initial susceptible individuals are around 46, 398), initial infected people I(0) = 1, and β is given by (9) for β 0 = 0.008 and η is taken such that β (1) = 26. Diffusion of susceptible individuals seem to be the most sensitive parameter in order to obtain either traditional SIR type behaviour or rebounds and plateau-type dynamics. Figure S9 illustrates the changes of dynamics depending on diffusion parameter d. The associated dynamics of susceptible individuals are represented in Figure S10 when variable a ∈ (0, 1) is discretized into n g = 50 groups through a finite difference scheme. Coloured curves represent each discretized group and the black curve the average of the groups (total susceptible individuals). Figure S10. Model behaviour depending on diffusion parameter values: 50 groups of susceptible individuals, in logarithmic scale. The black curve corresponds to total susceptible individuals and the bottom light blue curve is associated with the most risky behaviour (a close to 1) and may not be monotonically decreasing in time. From left to right and then top to bottom, graphs are associated with d = 10 −5 , d = 5 · 10 −5 , d = 10 −3 and d = 5 · 10 −3 (in day −1 unit) One can observe the effect of diffusion: while individuals with the riskiest behaviour are rapidly consumed, their group is renewed through diffusion from the majority of less infectious individuals. Depending on the parameter range of the diffusion coefficient d in (4), it leads to multiple epidemic waves, plateaus, or classical single wave SIR-like dynamics.

Stability of plateaus: the key role of diffusion
Let us illustrate in a heuristic manner the impact of diffusion on the existence and stability of plateaus in the dynamics of infectious individuals in Model (4). First, the second derivative of x = log I satisfies In the case d = 0, we deduce from (24) that d 2 x/dt 2 (t) < 0 for all t ∈ R + , so that x = log I is a strictly concave function of time. The consequence is that only slow decrease at the beginning of the decay phase can be obtained. Plateaus or shoulder-like patterns cannot arise in the case d = 0, nor rebounds since x(t) has a single maximum after which it keeps decreasing.
In the case d = 0, Equation (24) can be reformulated as

11/15
where Assuming that β and S 0 are respectively increasing and decreasing functions of a, we deduce from Proposition 1 that η(t) and I p (t) are always positive. Introducing x p (t) = log I p (t), Equation (25) can be rewritten as It means that d 2 x/dt 2 (t) may be positive (resp. negative) depending on whether x(t) is lower (resp. greater) than x p (t). In particular, in the neighborhood of times t such that x(t) is close to x p (t), so that oscillations may arise, leading to patterns such as the ones seen on Figure S9. These oscillations also explain rebounds, or shoulders-like patterns before plateaus. Of course, this is not a proof but an indication to this effect and it warrants further mathematical investigation.

More general systems
First, it should be noted that the coefficient β (a) can be time dependent: β = β (t, a). We can thus incorporate public health policies changes such as social distancing, lockdowns etc. For instance, we can consider a transition between two profiles β 1 (a) and β 2 (a) at a given date. Then, we observe that System (4) is a particular case of a more general class of models. Indeed, we may want to also consider heterogeneity in the compartment of infected. This for instance reflects the division between symptomatic and asymptomatic individuals, or other differences. In particular, in terms of behaviour one can then include effects such as "super-spreaders". We then get a model where the population of infected also depends on the same parameter a, that is I = I(t, a). We assume that the probability of interaction between a susceptible of class a and an infected of class b is given by a coefficient of transmission B(a, b). These considerations lead us to the following general system.
We observe that our model above (4) is derived from this more general class by assuming that B(a, b) is independent of b. In fact, the same is true for a finite number of states. This type of modelling of heterogeneity can be found in Magal et al. 6 and Arino et al. 9 in the discrete setting (a has a finite number of values) and without diffusion (d = 0). These works do not seem to yield shoulder, plateau and rebound-type patterns. For simplicity we assume that the decay rate γ is uniform. Note that the second equation states that new infected individuals of class a result from an infection of a susceptible individual of the same class. Other choices are possible, but then, in order to be consistent, one would need to introduce new compartments such as symptomatic/asymptomatic, hospitalised etc. Likewise, we could consider other variants of the SIR model, such as the SEIR model. It would then be natural to assume that exposed individuals also move around. We are then led to a system of the following kind.

12/15
We can also consider a more general diffusion process on the behavioural variables than Brownian diffusion. We already mentioned above (3) the version involving a diffusion ∂ 2 aa [σ 2 (a)S(t, a)] rather than d ∂ 2 aa S(t, a). With a somewhat different point of view, we can assume that there is a probability kernel K(a, b) that represents the probability distribution that an individual with trait b jumps to behaviour a. Then, the variation of S involves an integral operator rather than a heat operator. We are then led to the following system.
System (30) is supplemented with the same type of initial conditions as above (4). Note that there are no boundary conditions in this case.
It is classical that one can recover System (4) as a certain limit of the more general class of systems (30) (compare e.g. the article by the first author et al. 10 ). However, here it is interesting to look at more behavioural assumptions on the kernel K. For instance, it is natural to assume that K(a, b) is singular when b is close to a, that K(a, b) is rather large when b is close to 1, whereas it is very small when b is small. We leave the study of such kernels for future study.

Digital PCR in the framework of wastewater based epidemiology
Wastewater-based surveillance has been used over the past two decades to assess in real time the exposure of people to chemicals (in particular pollutants, licit and illicit or misused therapeutic or other drugs) and pathogens (bacteria or viruses) at the community level 11 . It provides a global picture of population health at the scale of an urban area connected to wastewater treatment plants (WWTP). Fully complying with the privacy of individuals, this approach rests on an explicit, objective observation. It also has the advantage to serve as an early warning system. Nonetheless, there are also uncertainties in the quantitative interpretation of these measurements. In particular, dilution due to storm-water, infiltration of parasite inflow water, varying population, bio-marker degradation, varying time spent in the flow..., all generate a variety of biases.
In spite of these uncertainties, there is a renewed interest in wastewater based epidemiology (WBE), in the context of the Covid-19 outbreak. There are numerous initiatives over the world, in Australia 12 , Japan 13 , Netherlands 14 , Detroit 15 , Paris 16 , Spain 17 and Canada 18 . The recent periodic sampling campaigns in WWTPs rely on classical quantitative real time PCR (qRT-PCR) to estimate the concentration of SARS-Cov-2. At this stage, these programs serve as early warning platforms based on the observation that they detect SARS-Cov-2 in sewers at least seven days ahead of individual testings 15,17,19 .
Only a few of them also used digital PCR on sub-samples to assess the measurement uncertainties of qRT-PCR measurements. Even if the benefit of dPCR in particular for low viral loads is mentioned 18 , it is not currently used routinely for systematic measurements, mainly for economic reasons.

Plateau-type patterns at country level
In most European countries, Covid-19 data exhibit similar plateau-type features, as illustrated in the so called effective R coefficient computed from the daily number of deaths 20 . This effective R coefficient exhibits phases where it stabilises around 1, which corresponds to an epidemic plateau.