The effects of heterogeneity and stochastic variability of behaviours on the intrinsic dynamics of epidemics

We report here on a campaign of weekly measurements of concentration of SARS-Cov-2 in wastewater in several treatment plants around the Thau lagoon in the Southwest of France over a nine month period of time. The use of Digital PCR yielded very precise measurements. The observations thus generated exhibit a rough stabilization on plateaus of the epidemic and other remarkable features. Such plateaus are widely reported in the setting of the Covid-19 pandemics. In this paper we raise the question of why such plateaus and other features of epidemics dynamics arise. Indeed, the classical SIR model and its extensions hardly provide an explanation for such behavior. To address this question we introduce here a new model, which takes into account heterogeneity and natural variability of behaviors in populations. Owing to this model, we show that features such as plateaus, rebounds, and shoulders are part of the intrinsic dynamics of an epidemic. In particular, in the context of the Thau lagoon, we argue that they are not generated by public health policy measures or psychological reactions of the population. We then show that this model fits very well the measures obtained around the Thau lagoon.


Introduction
The onset of plateaus for various indicators of the current outbreak of Covid-19 appears to be a rather general feature of its dynamics, along with periods of exponential growth or decay, rebounds etc. Public health and governments officials have kept mentioning them all along the pandemics. Among a large number we just quote a handful of such declarations in the SI Appendix. Nonetheless, there are few theoretical explanations offered to understand this phenomenon. Yet, such plateaus hardly agree with the classical SIR paradigm of epidemics. We show here that plateaus arise intrinsically in the unfolding of an epidemic. For this, we take into account two elements: an underlying heterogeneity and a random variability of behavior in the population. These features are of course more realistic than assuming that the population is perfectly homogeneous with an unwavering behavior. To this end, we propose in this paper a new mathematical model that extends the classical and simple SIR model of epi-demics. It takes the form of a system of reaction-diffusion equations, where one variable represents the behavior of individuals. More precisely, we assume that the population is structured into classes according to some continuous trait variable a. This implicit variable describes the more or less risky behavior individuals adopt. It can represent for instance the number of social contacts per day of an individual, taking into account their length, whether social distancing is observed, if it takes place in a indoor or outdoor, more or less crowded environment, or if individuals wear a mask etc. Thus, it can be seen as a lumped variable that represents a global score of risk individuals are taking. In this context, the impact of political restrictions can easily be represented by a modification of the distribution law of the behaviors. Furthermore, it is natural to consider that individuals may change their behavior from one day to the next one. We assume here that individuals' behaviors move randomly according to Brownian motion among these classes. We now have a model that includes heterogeneity and variability of behaviors. We show here that such a system exhibits a richness of dynamics and in particular gives rise to intrinsic formations of plateaus. Up to now, there are two main alternative explanations for the onset of plateaus. The first one is political (1). By managing the epidemic and keeping the exponential growth at bay, without destroying the economy, a plateau appears as some kind of optimal compromise under constraint. Another approach appears in a very recent work of Weitz et al. (2). It develops a SEIR-type compartmental model with a global transmission rate coefficient that varies in time. The authors argue that plateaus are caused by the evolution of behaviors due to two competing main psychological reactions to the epidemic. These play out awareness of fatalities and fatigue of the public facing regulatory mobility restrictions, which impact the transmission rate coefficient of the disease. In recent works, Arthur et al. (3) and Radicchi et al. (4) proposed models in the same spirit. We refer to Section 6 and SI Appendix for other references and more detailed discussion of related literature. It is to be expected that the heterogeneity we include in the model we present here, together with stochastic variability, will yield an even richer dynamics. We describe and dis-cuss the model in Section 3. To discuss the validity of our approach, we rely on observations akin to a natural experiment that we refer to here as the "Thau lagoon experiment". We show the outcomes in Section 2. The key to obtain such precise results is the use of the "digital PCR" method of analysis. This series, over an extended period, strikingly reveals the formations of plateaus, in some cases after a "shoulder" pattern. The appearance of such plateaus and other features need not be explained neither by psychological reactions like in Weitz et al. (2) nor by public health policy effects. Indeed, the policy regulations were roughly constant during the measurement campaign and awareness or fatigue effects do not seem to have altered the dynamics over this long period of time. Witness to this is the fact that two groups of towns saw the same evolution, but two weeks apart one from the other.
In Section 4, we report on the calibration of our model on the data of the Thau lagoon. It yields a remarkable fit. Thus, we show here that the model involving heterogeneity of the population behavior with respect to risk and the variability of behavior exhibits an intrinsic dynamics giving rise to plateaus. We further discuss our model in more detail in the light of the measurements in Section 5 below.

The Thau lagoon experiment
The measurement campaign concerned four wastewater treatment plants in the Thau lagoon area, serving the cities of Sète, Pradel-Marseillan, Frontignan and Mèze. The measurements were obtained by using digital PCR (5) (dPCR) to estimate the concentration of SARS-Cov-2 virus in samples taken weekly from 2020-05-12 to 2021-01-12. We provide further details about the measurement method in the Methods section.  Observations of the covid-19 outbreak unfolding in the Thau lagoon. Figures 1 and 2 show the outcomes in a logarithmic scale. The above measurements over a nine months period show some striking features that we summarize now.
1. An exponential phase starts simultaneously in Mèze and Frontignan WWTPs in early June.
2. The genomic units concentration curves in these two places reach, again simultaneously, a plateau. It has stayed essentially stable or slightly decreasing since then.
3. The evolution at Sète and Pradel-Marseillan remarkably followed the previous two zones in a parallel way, with a two weeks lag. The measurements at Sète and Pradel-Marseillan continued to grow linearly (recall that this is in log scale, thus exponentially in linear scale), while the Mèze and Frontignan figures had stabilized ; then, after two weeks, they too stabilized at a plateau with roughly the same value as for the other two towns.
4. The measurements seem to show a tendency to increase over the very last period.

An epidemiology model with heterogeneity and natural variability of population behavior
Numerous papers (6-9) describe the number of daily social contacts as a key variable in the spread of infectious diseases like Covid-19 insofar as it is closely related to the transmission rate. Daily social contacts are usually described in terms of age, gender, income, type of job, household size (10). . . Parameters that are particularly relevant in the context of viral outbreaks are also studied (11) such as cumulative duration of such contacts, social distance, indoor/outdoor environment. . . . The very fine grain microscopic models (12) aim at identifying such parameters in the most precise way possible. A very recent study by Di Domenico et al. (13) takes up the data of hospitalizations in France for the past six months. Again, these exhibit striking epidemic plateaus since the beginning of 2021. The authors of this paper provide a microscopic insight of the propagation, emphasizing the role of two different strains of the virus and the role of public health measures such as the curfew, school closing etc. These approaches are different from ours, as we rather adopt here a mesoscopic point of view. Given the complexity and multiplicity of behavioral factors favoring the spread of the epidemic, we assume that the transmission rate involves a normalized variable a ∈ (0, 1) that defines an aggregated indicator of risky behavior within the susceptible population. Thus, we represent the heterogeneity of individual behaviors with the variable a. We think of it as encompassing not only the number of daily social contacts, but also their duration, distance, environment etc. Thus, we take a as an implicit parameter that we do not seek to calculate. The classical SIR model is macroscopic and the type of model we discuss here can be viewed as intermediate between macroscopic and microscopic and this is why we call it mesoscopic. The initial distribution of susceptible individuals S 0 (a) in the framework of a SIR-type compartmental description of the epidemic can be reasonably taken as a decreasing function of a. We take the infection transmission rate a → β(a) to be an increasing function of a. In the SI Appendix, the interest reader will find a more general version of this model involving a probability kernel of transition from one state to another. The model here can be derived as a limiting case of that more general version. Likewise, the behavior of individuals usually changes from one day to another (14). Many factors are at work in this variability: social imitation, public health campaigns, opportunities, outings, the normal variations of activity (e.g. work from home certain days and use of public transportation and work in office on others) etc. Therefore, the second key feature of our model is variability of such behaviors: variations of the population for a given a do not only come from individuals becoming infected and leaving that compartment but also results from individuals moving from one stratus a to another (14). In the simplest version of the model, variability is introduced as a diffusion term in the dynamics of susceptible individuals.
The model. We denote by S(t, a) the density of individuals at time t associated with risk parameter a, by I(t) the total number of infected, and by R(t) the number of removed individuals. We are then led to the following system: where γ −1 denotes the inverse of typical duration (in days) of the disease and d a positive diffusion coefficient. System [1] is supplemented with initial conditions and homogeneous Neumann (zero flux assumption) boundary conditions Note that the total population N is conserved by the dynamics: In the SI Appendix, we describe a more general version where the term in the first equation of [1], ∂ 2 aa S(t, a) is replaced by a heterogeneous term ∂ 2 aa [σ 2 (a)S(t, a)], with a rate of change of behavior that depends on a. We observe that in this model, we adopt the point of view that the behavior distribution affects the risk takers in the susceptible population rather than the behavior of infected individuals. Indeed, it appears more natural to consider that the susceptible face an ambient distribution I(t). Recall that, for instance, the choice to go to a crowded bar where there might be a super-spreader, is reflected in the variable a. Nonetheless, we can also write a similar model where the risky behavior structures the infected population. More generally, in the SI Appendix, we describe a model where both the susceptible and the infected populations are structured by risk. Several earlier works have considered SIR type systems with heterogeneity. In particular, Arino et al. (15), and, more recently, Dolbeault and Turinici (16,17), Magal et al. (18) have studied models with a finite number of different coefficients β. These systems have a discrete set of classes and do not involve variability. Almeida et al. (19) considered the case of continuous classes associated with a multidimensional trait x to study the influence of variability of infectious individuals on the final size of an epidemic. DiMarco et al. (20) have recently proposed a model close in spirit to ours but more complex. Based on kinetic theory, it describes the heterogeneity of individuals in terms of a variable x ≥ 0 corresponding to the number of daily social contacts. Their model consists of a system of three SIR equations coupled with Boltzmann or Fokker-Planck type equations. The authors emphasize a collision type term in the resulting Boltzmann equations. This term represents changes of behaviors (that is of the values of x) when two individuals meet. Consequently, it ends up as quite different from the present work. As a closure of their model, DiMarco et al. (20) formally derive a so called S-SIR model. There, the variability of behaviors involves an explicit dependence of the transmission rates upon the current number of infectious individuals (21). Because of this feature, as a matter of fact, when applied to real data, it is eventually very similar to the approach in Weitz et al. (2) 4 Application to the Thau lagoon measurements The above model [1] describes the dynamics of infected fraction of population t → I(t)/N . In order to represent the dynamics of concentration of SARS-Cov-2 as measured in WWTPs, the effective proportionality coefficient λ minimizing the least square distance between measured concentrations and λ multiplied by the rate of infected individuals is computed in the calibration process over the time interval from 9 June 2020 to 12 January 2021. Calibration of model [1] also requires fitting the values of γ, the profiles a → β(a) and the initial distribution of susceptible individuals in terms of a. Satisfactory capture of the exponential phase and the subsequent plateau leads to coefficient λ = 111, 230, 001 (in genomic unit per liter). The underlying dynamics of the rate of susceptible individuals is given in Figure 4 below for n g = 20 groups. The lower curve illustrates the competition phenomenon between diffusion and sink term due to new infections, depending on the level of risk a of each stratus. Note that Model [1] assumes that the modeled area is a closed system, with constant population (no outgoing or incoming population). The application on Sète area (it may be equivalently done on the other three WWTP measurements) is therefore questionable, since the measurement period coincides with the summer holiday season in France. Still, the above calibration seems satisfactory to the extent that it captures both the dynamics of the exponential growing phase as well as the subsequent shoulder-like and plateau phase. The resulting effective coefficient λ somehow embeds the combination of complex phenomena involved in Equation [4], including virus degradation in sewer environment.

Discussion
We claim that the formation of plateaus and rebounds in the dynamics of the outbreak is explained by Model [1] through the heterogeneity and variability of population behavior with respect to epidemiological risk. Figure 4 shows that susceptible individuals with the riskiest behavior, characterized by the highest β transmission coefficient, are rapidly transferred to the infected compartment. Variability of behaviors modeled by diffusion with respect to a parameter then re-feeds the fringe of riskiest susceptible individuals. Parameter regimes where those two phenomena have the same order of magnitude generate patterns such as plateaus, shoulders and rebounds. Figures S9, S10 and S11 in the SI Appendix illustrate these various dynamics.
In Figures S9 and S10 in SI Appendix, we explore the types of patterns arising when we vary the diffusion coefficient d. We note that plateaus occur for an intermediate range of values of d > 0, which represents the amplitude of variability. Namely, large values of d lead to standard SIR-type dynamics, since diffusion quickly homogenizes the behaviors. When d decreases, plateaus starting with shoulder-like patterns appear. However, for even smaller values of d, oscillations arise, which can be interpreted as epidemic rebounds. From numerical experiments, the amplitude of rebounds always seem to be of smaller amplitude than the first epidemic peak. However, higher secondary peaks arise when we significantly modulate in time the transmission rate. This may represent a progressive exit from lockdown or the effect of new and more contagious variants (see Figure S11 in SI Appendix).
The dynamics of the Covid-19 outbreak in the cities of Thau lagoon appears almost insensitive to public health regulations. In particular the second lockdown in France from 28 October to 14 December 2020 had hardly any effect. Likewise, the Christmas Holiday season also seemed to have had little influence on the observed plateaus.
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 (22). This effective R coefficient exhibits phases where it stabilizes around 1, which corresponds to an epidemic plateau. Effective R coefficient as a function of time (23) https://shiny.biosp.inrae.fr/app_direct/mapCovid19/.
Campo et al. (24) study more specifically the dynamics of Covid-19 in Italy and put emphasis on plateaus and multiwave patterns, for which they discuss meta-stability properties. The number of hospitalized individuals in France over the last quarter of 2020 is represented in Figure 6. There, the dynamics shows a growing phase followed by a shoulder and a plateau, very similar to the pattern observed in Thau lagoon. Several noteworthy observations came out from the Thau lagoon experiment. First, we can see two distinct exponential growths in two separate subsets of two towns. The two graphs are parallel with a constant delay of two weeks. The first group reaches a plateau at a certain level of infected and stays thereafter essentially flat, while the second group continues to grow until it reaches about the same value of infected and then becomes essentially flat too. These remarkable observations call for interpretations. Indeed, first, they cannot simply result from public policy measures as these would have affected all these neighboring towns in a similar fashion. What we observed, however, is that the second group continued to grow for two weeks until it reached about the same level of infected as the first group, which had stayed flat in the meantime. Second, there likely is a spatial diffusion effect that triggered the growth in the second group coming from the first one. However, this diffusion would not explain the fact that the two groups reached the same level of plateaus. It is worth noting that these observations are inconsistent with the classical SIR model: in this model, once an epidemic reaches a peak, it then decreases steadily by another exponential factor. In other words, a plateau requires an effective R t number (23, 25) of approximately one: R t ∼ 1. However this would be hard to sustain over such a long period of time as observed because the susceptible population is depleted. It also appears difficult to explain that this occurs exactly at the same plateau level for two distinct populations.

Conclusions
In this paper we reported on precise concentration measurements of SARS-CoV-2 in wastewater upstream of four WWTPs around the Thau lagoon in South of France during the nine months period May 2020 to January 2021. These were made possible owing to dPCR technique, used here in a systematic fashion for intensive epidemiological analyses. Such virus concentration time series, essentially proportional to the fraction of people infected by Covid-19, provides an accurate quantitative method to monitor the epidemic. Furthermore, as mentioned in (26-28), the appearance of the virus in wastewater precedes the observations of the disease and therefore yields a remarkable early warning system, ahead of hospital counts. Moreover, it reflects all infectious people regardless of whether they are symptomatic or asymptomatic. It transpires from the Thau lagoon observations that plateaus and other characteristics result from the intrinsic dynamics of epidemics. To wit, there are many reports of such plateaus being observed in the unfolding of epidemics, for instance at country levels. We provide here an explanation of this phenomenon by considering that the population's behavior is heterogeneous in that individuals differ by the level of risk in their behavior. Next, individuals (or part of them) change randomly their behaviors from day to day. We show here that such heterogeneity and variability favors the possibility of observing a stabilization at a near 1 effective R t number. The reason is that the combination of heterogeneity and variability leads to a constant replenishment of the population of susceptible individuals that yields plateaus in the dynamics of the epidemics. To prove this claim, we propose a mathematical model for epidemics that explicitly involves heterogeneity and variability. The model takes the form of a diffusion SIR system with diffusion of behaviors. Here, to simplify as much as possible, we only assume the diffusion of risks among susceptible, disregarding the similar effect among infectious. The numerical simulations of this system indeed bring to light exactly this type of shape that we observed with shoulders, long plateaus and also rebounds. We were able to calibrate this model on the Thau lagoon observations and producing a curve that fits remarkably well the data.
Perspectives and extensions of the model. Besides the mathematical extensions we discuss in SI Appendix, there are many directions where one would want to go further in the development of this model. First, one might envision a more precise approach to the parameter a. The notion of unit of contact requires further discussion. One can take into account e.g., duration and circumstances of contacts. It appears natural, to use "risky sociability" rather than mere sociability in computing the transmission rate. Understanding quan-titatively the change of behavior as a function of this variable a appears to be an interdisciplinary avenue of research involving sociological and biological studies to get a handle on possibly more precise quantitative forms. In order to achieve a greater precision, one could also envision a multidimensional variable a = (a 1 , . . . , a p ) to encompass various behavioral 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 ∆ a S. Such developments would be very useful to improve the model and its applicability. Likewise, it is quite natural to assume that the variability of behavior also depends on a rather than being uniformly distributed. We discuss such an extension in the SI Appendix where we show that it leads to an equation with drift terms. We leave to further work developments along this direction. Another important aspect that transpires in the Thau lagoon data, is the spatial spreading that takes place in the epidemic. In a recent paper, the first author of this study and collaborators (22) have proposed a model at the country level with a quantitative approach to spatial diffusion in France. The study of diffusion at a smaller scale that we might call mesoscopic and its inclusion in the framework we propose here are promising perspectives.

SARS-Cov-2 concentration measurement from
wastewater with digital PCR. The teams of the local authorities of the Thau lagoon, Syndicat Mixte du Bassin de Thau (or SMBT) collected samples every Tuesday from each of the four WWTPs. Each sample consisted of a compound of 24 hourly samples. I.A.G.E. 1 developed a diagnostic method to detect very low concentrations of SARS-CoV-2 in such wastewater samples. This method combines an optimized extraction process with a DNA quantification based on a digital PCR (dPCR) targeting region of the RdRp (IP2/IP4) (5) 2 . The measures produced identical results within three significant digits between the two targets. In contrast to classical quantitative real-time PCR (qRT-PCR), dPCR allows the absolute quantification of low concentration levels of target sequences of nucleic acid molecules from DNA or RNA samples. dPCR outperforms qRT-PCR with respect to accuracy (29) and repeatability of measurements; it is also has a much lower detection threshold (about 20 times lower) (5). Among recently achieved wastewater measurement campaigns in sewers (26-28, 30-34), it seems that none of them exploit these measurements in a dynamical model. This is probably due to the uncertainties associated with qRT-PCR measurements and to the difficulty to translate genomic unit concentrations into numbers of infected individuals. As out- The approach developed here is different from these earlier works and attempts to take advantage of an original modeling approach of population heterogeneity in order to estimate the effective parameter λ relating the rate of infected people to the measured concentration of SARS-Cov-2. A satisfactory capture of the dynamics of concentrations by the model, including the exponential growth followed by the formation of shoulder-like and plateau pattern leads to such effective parameter estimation, as described in Section 4.
An epidemiological model with heterogeneity and variability. The model we develop here extends the classical SIR compartmental approach (compare with Equation [1] in SI Appendix) by taking into account heterogeneity and variability of behaviors on the susceptible population. First, we describe the heterogeneous behaviors by a variable a that we discuss below. The total population is assumed to be constant equal to N . That is, we do not take into account incoming or outgoing populations, nor demographic changes. This is a standard assumption in the Covid-19 studies. Actually, one may want to dispense with it if one is to consider a significant amount of travelers, especially during vacation periods or because there is a lockdown that brings many people to leave large cities. But here, we keep this constant population assumption.
As usual, the population is divided into the three compartments of susceptible, infected and removed individuals. At time t, the population is thus structured by a and t and is described by its density S(t, a). We assume that the total number of infected is given by I(t) and that of the removed by R(t). We do not distinguish from which population strata the infected individuals come from. In the Conclusion section among the perspectives, we discuss some extensions of this framework with a heterogeneity of behaviors of the infected.
One can think of a as a trait parameter roughly describing the level of risk a given susceptible individual is taking with respect to Covid-19. It is normalized so that 0 ≤ a ≤ 1, a = 0 being associated with very cautious people and a = 1 to people with very risky behavior. This abstract variable a is an aggregation of multiple risk factors for a disease such as Covid-19: number of contact per day, duration and distance of daily social contacts, indoor/outdoor environment, respect of barrier gestures and mask wearing. . . Thus, taking into account heterogeneity leads to the following system: Note that the total population is constant: for all times. The previous model assumes a "static" behavior of individuals. Once it is in a given class a, an individual remains as it were in this class forever or until it becomes infected. But this not realistic and we now take into account variability. Individuals often vary their behavior, because for instance of fatigue effects for people who have heeded too strongly social distancing calls or, on the contrary, for people who have been reckless and see other people fall sick and consequently become somewhat more cautious. But also, its variability has more mundane causes which make this variability rather usual. For example, certain days, a person may be secluded while on other days, she may have professional meetings. Other persons may want to go shopping certain days or may want to go to the hairdresser or physiotherapy some days even if trying to minimize these occasions. What we are saying is that the potential reservoir of individuals for a given stratus level a is not static and there is more important than it would appear at first glance. So, here comes the second ingredient of our model, the variability. That is, we take into account that the variations of S(t, a) do not only come from individuals becoming infected and leaving that compartment but also results from individuals moving from one stratus to another. A natural case to consider here is to assume that this shuffling of behaviors happens according to Brownian motion. We are then led to System [1] presented in Section 3 above that takes into account behavioral heterogeneity and variability. In the SI, we derive this model from more basic considerations. We also describe some mathematical properties of this system. It is enlightening to keep track of the fraction of infected coming from specific strata. To this end, we can introduce a variable I(t, a) representing the number of infected that came from stratus a. It is given by the solution of the equation: It is straightforward that this is consistent with the definition of the total population of infected, that is: In the computations shown later on, we represent the various fractions of I corresponding to the strata they originate from. We may include then here the effect that individuals with high a's are more likely to infect other susceptible. This just amounts to replace I(t) by  https://www.gouvernement.fr/partage/ 11951-conference-de-presse-du-premier-ministre-point-de-situation-sur-la-lutte-contre-la-covid-19.

Supplementary Information
The classical SIR model with time dependent transmission rate. The classical SIR compartmental model (39) 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: Single location growth and behavioral changes.. Changes of behavior driven by phenomena such as awareness of fatalities or fatigue with respect to mobility restrictions can be modeled by a time modulation of the infection rate β.
To start with, we analyze 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 optimization 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 7 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: Berestycki et al. | The effects of heterogeneity and stochastic variability of behaviours on the intrinsic dynamics of epidemics MedRχiv | 9 All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Fig. 8. Rate of infected individuals as a function of days
We developed a toolbox based upon PYGMO (40) particle swarm optimization 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 behavior 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. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  We developed a similar approach with noise added to sampled synthetic data, with satisfactory convergence both in terms of dynamics of infection rate ( Figure 11) and modulation functions (levels and time changes in Figure 12). All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  A model with change of behavior at discrete times in Thau lagoon. 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. Fig. 13. Change of behavior 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 on August 10 th , 2020. The resulting piece-wise constant modulation function in time is represented in Figure 14 below.

| MedRχiv
Berestycki et al. | The effects of heterogeneity and stochastic variability of behaviours on the intrinsic dynamics of epidemics All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  Even though the plateau dynamics seems to be correctly represented in Figure 13, the main shortcoming of this is that it does not really explain the observations of the Thau lagoon experiment. 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 Model [1]
. This section is concerned with the formal derivation of the additional diffusion term in the dynamics of susceptible individuals. As previously emphasized, this diffusion term stems from the combination of the heterogeneity of behaviors and their variability in time. A natural case to consider here is to assume that the shuffling of behaviors 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 In the simple case where σ 2 (a) = d where d is a constant, one ends up with the first equation of [10] by inserting [9] into the heterogeneous SIR model [11] ∂ t S(t, a) = d ∂ 2 aa S(t, a) − β(a)S(t, a) The numerical and theoretical analysis of the case where σ depends on a, and the corresponding extensions of the Fokker-Planck equation with drift terms, which seems relevant from a modeling viewpoint, will be addressed in a forthcoming work. All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted March 29, 2021. ; This paper of Weitz et al. is quite different from our model in that it assumes a homogeneous population: at a given time, all individuals have the same transmission rate. The common behavior simply changes in time by reacting to the outcome of the epidemic and this change reflects in the evolution of the transmission rate. Moreover, we note that the model is calibrated with reported death of various U.S. states and does not involve wastewater concentration measurements. Thus, even though in a different manner from ours, this work also stresses the role of variability for the observed dynamics. Let us also mention another work by Dimarco et al (20) inspired from kinetic theory. Its spirit is akin to ours, but with a different scope. It describes the heterogeneity of individuals in terms of a variable x ≥ 0 corresponding to the number of social contacts. The model then consists of a system of three SIR equations coupled with Boltzmann or Fokker-Planck type equations, local equilibria being expressed as Gamma-type distributions. As a closure of this model, the authors derive a so called S-SIR model where variability of behaviors is expressed as a dependence of mean number of daily contacts upon infection rate. Thus, in the end, this approach is actually quite similar to that of Weitz et al. (2).
Simple properties of Model [10]. For the convenience of the reader, we recall the model introduced in the main paper We also recall the case d = 0 (heterogeneous but without variability): 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, the model [10] may exhibit more sophisticated behaviors such as rebounds in I with multiple local maxima as well as plateaus. This property can be intuitively understood by analyzing the evolution of the growth factor in [12]: 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 characterized by the property S 0 (a)β (a) < 0 such as S 0 (a) = S 1 + (S 0 − S 1 ) exp − a 2 κ 2 for some κ > 0 and 0 < S 1 < S 0 .
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 behavior. To some extent, social contact surveys (6,7,14) justify an initial distribution of susceptible individuals like [16]. However, we are not aware of a rigorous justification of the expression [15] 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) := ∂ a S(t, a) of S with respect to a satisfies the following equation S(0, a).
Since the right hand side of the first equation of [17] 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 [10]. Let us first notice that we can describe the dynamics of susceptible individuals in [10] in terms of its probability density f (t, a) = S(t, a)/S(t) where S(t) = 1 0 S(t, a) da: 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 fueling as it were the epidemic with individuals who had initial lower risk trait. The equivalent of [14] for the average transmission rateβ is In other words, in the absence of behavioral variability (d = 0), the average of the transmission rate decreases under the effect of its variance, the latter being directly linked to behavioral 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 [19] 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.
Theorem 1. Let (t, a) → (S(t, a), I(t), R(t)) be the solution of [11] (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 [8] 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. (18), Dolbeault and Turinici (16) and Almeida et al. (19). 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 [20] 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 [8] with initial conditions (S 0 ,Ĩ 0 , 0) and homogeneous parameterβ 0 , which allows to conclude. An interesting property of System [11] is the conservation of the following entropy-like quantity The property V (t) = V (0) can be deduced from [11]. 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 [10] is non-decreasing in a similar manner as in the framework of viscous perturbations of scalar conservation laws: The following result emphasizes 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 [10]. The following properties hold: a) is not necessarily decreasing in time for a ∈ (0, 1) given. 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 [10] 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.
(ii) Let (S(t),Ī(t),R(t)) be solution of the homogeneous SIR model [8] with initial conditions (S(0) = 1 0 S 0 (a) da,Ī(0) = I 0 ,R(0) = 0) and transmission rate β m = 1 0 β(a) da. Then, denotingS ∞ the large time limit of t →S(t), one has S d ∞ >S ∞ . 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 fueling 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 [10] multiplied by ∂ t S leads to: so that we have (27) 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 [27] generalized 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 [27], 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 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 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 [23] 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 [16] 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 [15] 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 behavior or rebounds and plateau-type dynamics. Figure 15 illustrates the changes of dynamics depending on diffusion parameter d.
The associated dynamics of susceptible individuals are represented in Figure 16 when variable a ∈ (0, 1) is discretized through a finite difference scheme with n g = 50 groups. Colored curves represent each discretized group and the black curve the average of the groups (total susceptible individuals).

| MedRχiv
Berestycki et al. | The effects of heterogeneity and stochastic variability of behaviours on the intrinsic dynamics of epidemics All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  One can observe the effect of diffusion: while individuals with the riskiest behavior 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 [10], it leads to multiple epidemic waves, plateaus, or classical single wave SIR-like dynamics.
Because of the dissipative nature of Model [10], secondary epidemic peaks are of lower amplitude than the first one, as shown in the top left graph in Figure 15. Another simulation derived from the same parameters as the top left graph in Figure 15 consists in modulating in time the β coefficient as a collective homothetic behavior change (β(a) being replaced by β(a)ϕ(t)). Such collective change may be interpreted as reflecting an effect of fatigue with respect to mobility restrictions or a progressive exit from a lockdown period as from Day= 1000. In the simulation below (Figure 17), ϕ(t) = 1 + 3 max(0, t − 1000)/1000, leading to second peak of infectious of higher amplitude than the second one, followed by a plateau. so that oscillations may arise, leading to patterns such as the ones seen on Figure 15. 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). For instance, this is needed to incorporate public health policies that change at a given date. In this case, we can consider for instance a transition between two profiles β 1 (a) and β 2 (a) at a given date. Then, we observe that System [10] 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 behavior 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  β(a, b). These considerations lead us to the following general system.
We observe that our model above [10] is derived from this more general class by assuming that β(a, b) is independent of b. In fact, the same is true for a finite number of states. This type of modeling of heterogeneity can be found in Magal et al. (18) and Arino et al. (15). Still, such works do not embed variability of such heterogeneity, which does not allow for shoulder, plateau and rebound-type patterns. For simplicity we assume that the decay rate γ is uniform but it would not change much if we were to assume that γ = γ(a). 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, hospitalized etc. Likewise, we could consider other variants of the SIR model, such as the SEIR model. It would then be natural to assume that infectious also move around. We are then led to a system of the following kind.
We can also consider a more general diffusion process on the behavioral variables than Brownian diffusion. We already mentioned in the main article a 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 at b jumps to behavior a. Then, the variation of S involves an integral operator rather than a heat operator. We are then led to the following system.
System [10] is supplemented with the same type of initial conditions as above [36]. Note that there are no boundary conditions in this case. It is classical that one can recover System [10] as a certain limit of the more general class of systems [36] (compare e.g. (41)). However, here it is interesting to look at more behavioral 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 (42). 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 (30), Japan (31), Netherlands (28), Detroit (32), Paris (33), Spain (27) and Canada (34). 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 (26, 27, 32). 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 (34), it is not currently used routinely for systematic measurements, mainly for economic reasons.  (46)).

| MedRχiv
Berestycki et al. | The effects of heterogeneity and stochastic variability of behaviours on the intrinsic dynamics of epidemics All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.