The Thau lagoon data
The measurement campaign concerned four wastewater treatment plants (WWTP) in the Thau lagoon area in France, serving the cities of Sète, Pradel-Marseillan, Frontignan and Mèze. The measurements were obtained by using digital PCR20 (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.
Concentrations of SARS-CoV-2 (genome units per litre in logarithmic scale) from four WWTPs in Thau lagoon, measured weekly with dPCR technology from May 12th 2020 to January 12th, 2021. Note that there are some missing points.
Figure 1 shows the outcomes in a logarithmic scale over a nine months period. We summarise now their main features.
- 1.
An exponential phase starts simultaneously in Mèze and Frontignan WWTPs in early June.
- 2.
The genome 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 have stabilised ; then, after two weeks, they too stabilised 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.
The epidemiological model with heterogeneity and natural variability of population behaviour
The appearance of such plateaus and shoulders need not be explained either by psychological reactions or by public health policy effects. Indeed, the 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. To understand this phenomena we propose a new model.
Given the complexity and multiplicity of behavioural factors favouring the spread of the epidemic, we assume that the transmission rate involves a normalised variable (a in (0,1)) that defines an aggregated indicator of risky behaviour within the susceptible population. Thus, we represent the heterogeneity of individual behaviours with this variable. 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.
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 mapsto beta (a)) to be an increasing function of a. In the Supplementary Information (SI) Appendix, the 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 behaviour of individuals usually changes from one day to another21. 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 behaviours: variations of the population density for a given a do not only come from individuals becoming infected and leaving that compartment but also results from individuals moving from one state a to another21. 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:
$$begin{aligned} frac{{partial S(t,a)}}{{partial t}} & = d{mkern 1mu} frac{{partial ^{2} S(t,a)}}{{partial a^{2} }} – beta (a)S(t,a)frac{{I(t)}}{N} frac{{{text{d}}I(t)}}{{{text{d}}t}} & = frac{{I(t)}}{N}{mkern 1mu} intlimits_{0}^{1} beta (a)S(t,a);da – gamma I(t), frac{{{text{d}}R(t)}}{{{text{d}}t}} = & gamma I(t), end{aligned}$$
(1)
where (gamma) denotes the inverse of typical duration (in days) of the disease and d a positive diffusion coefficient. System (1) is supplemented with initial conditions
$$begin{aligned} S(0,a) = S_0(a), quad I(0) = I_0, quad hbox {and} quad R(0) = 0, end{aligned}$$
(2)
and with zero flux condition in a at (a=0, 1). In the “Methods” section below, we discuss the relation of this model with other current works.
A more general model
In a more general version of our model, we can consider the population of infected as also structured by the parameter a. The equations are as before but now we keep track of the class a in the infected population. The mechanism here is that a susceptible individual from class a can be infected by infectious from any class I(t, b) but then gives rise to an individual I(t, a) of the same parent class. We also assume that there is a diffusion of the infected behaviours. We denote by ({mathfrak {B}}(a,b)) the transmission rate of S(t, a) by I(t, b). For simplicity and because it is natural, we will assume that it is of the form
$$begin{aligned} {mathfrak {B}}(a,b)= beta (a) beta (b) end{aligned}$$
where (beta) is as before. For full generality, we can also envision multi-dimensional parameters (ain {mathbb {R}}^d), with (a_iin (0,1)). We are then led to the system:
$$begin{aligned} frac{{partial S(t,a)}}{{partial t}} & = d;Delta _{a} S(t,a) – S(t,a)frac{{beta (a)}}{N}intlimits_{0}^{1} beta (b)I(t,b);db frac{{partial I(t,a)}}{{partial t}} & = d^{prime}Delta _{a} I(t,a) + S(t,a)frac{{beta (a)}}{N}intlimits_{0}^{1} beta (b)I(t,b)db – gamma I(t,a), frac{{{text{d}}R(t)}}{{{text{d}}t}} & = gamma intlimits_{0}^{1} I (t,b){mkern 1mu} db, end{aligned}$$
(3)
In the SI we write further, more general, forms of this model, with ({mathfrak {B}}(a,b)) and more general diffusion of behaviours, that can include jumps or non-local variations. The type of models we discuss here may also shed light on the initial phase of the epidemic. We plan to investigate these questions in future work.
Patterns generated by the model
In the next section, we will discuss how the model fits the data observed in the Thau lagoon measurements. But before that, we start by showing that the above model (1) can generate the different patterns we mentioned. For this we rely on numerical simulations without fitting real data. And indeed we obtain plateaus, shoulders, and oscillations. The latter can be interpreted as epidemic rebounds.
The key parameter here is the diffusion coefficient d, which controls the amplitude of behavioural variability (see Fig. 2). Large values of d rapidly yield homogenised behaviours, leading to classical SIR-like dynamics of infectious individuals. For very small values of d, the system also has a simple dynamics, in the sense that I(t) has a unique maximum, and therefore has no rebounds. We derive this in the limit (d=0) for which we show in the SI that there are neither plateaus nor rebounds.
For some intermediate range of the parameter d, plateaus may appear after an exponential growth, like in the initial phase of the SIR model. A small amplitude oscillation, called “shoulder”, precedes a temporary stabilisation on a plateau, followed by a large time convergence to zero of infectious population. We also show that for small enough d, time oscillations of the infectious population curve, i.e. epidemic rebounds, may be generated by Model (1). Such oscillations also appear after a plateau, in a similar way to what one can see in observations.
Simulations in Fig. 2 illustrate the various patterns obtained on the dynamics of infected population as a function of the diffusion parameter. For small enough d, in the top left graph of Fig. 2, one can see oscillations of the fraction of infectious individuals. These oscillations cannot be achieved in the classical SIR model. In fact, the two lower graphs of that figure, for somewhat larger values of d, exhibit the SIR model outcomes. Indeed, for sufficiently large d, the system becomes rapidly homogeneous (i.e. constant with respect to a). Yet, such oscillations are standard in the dynamics of actual epidemics, like the current Covid-19 pandemic. The intermediate value of d, represented in the upper right corner of Fig. 2 shows the typical onset of a plateau at a rather high value of I. Note that this plateau is preceded by a first small dip and then a characteristic “shoulder-like” oscillation.
Model behaviour depending on diffusion parameter values: infected rate dynamics in logarithmic scale. From left to right and then top to bottom, graphs are associated with (d=10^{-5}), (d=5times 10^{-5}), (d=10^{-3}) and (d=5times 10^{-3}) (in (day^{-1}) unit).
Secondary epidemic peaks are of lower amplitude than the first one, as shown in the top graphs of Fig. 2. This empirical observation leads us to conjecture that, at least in many cases, it is a general property of this model (with (beta) independent of time). This property would then reflect a kind of dissipative nature of Model (1). It is natural to surmise that a change of behaviours in time may generate oscillations with higher secondary peaks. Such changes result for instance from lifting social distancing measures or from fatigue effects in the population.
We illustrate this with numerical simulations in Fig. 3. We assume a collective time modulation of the (beta (a)) transmission profile. That is, we replace (beta (a)) by (beta (a)varphi (t)) for some time dependent function (varphi), the other parameters are the same as in the simulations shown in Fig. 2. We look at the effect of a “lockdown exit” type effect. Then, (varphi (t)) takes two constant values, 1 from (t=0) to (t={1000}) and 1.2 after (t={1100}). In between, that is, for (tin ({1000}, {1100})), (varphi (t)) changes linearly from the value 1 to 1.2.
Multiple epidemic rebounds: susceptible individuals are divided into 50 discrete groups in the case where relaxation of social distancing measures starts on Day (t=1000) and ends up on Day (t=1100). The fraction of infected individuals in the population is represented in the left graph in logarithmic scale and in linear scale in the right graph.
One can clearly see a secondary peak with higher amplitude than the first one. Note that after this peak, a third one occurs, with a lower amplitude than the second one. This third peak happens in the regime when (beta) is again constant in time.
The effect of variants
Another important factor that yields secondary peaks with higher amplitudes is the appearance of variants. Consider the situation with two variants. We denote by (I_1(t)) and (I_2(t)) the corresponding infected individuals. The first variant, which we call the historical strain, is associated with (beta _1) and (I_1(0)) and starts at (t=0). The variant strain corresponds to (beta _2) and (I_2) and starts at Day (t=1000). In this situation, the system (1) is extended by the following system:
$$begin{aligned} frac{{partial S(t,a)}}{{partial t}} & = d{mkern 1mu} frac{{partial ^{2} S(t,a)}}{{partial a^{2} }} – left( {beta _{1} (a)I_{1} (t) + beta _{2} (a)I_{2} (t)} right)frac{{S(t,a)}}{N}, frac{{{text{d}}I_{2} (t)}}{{{text{d}}t}} & = frac{{I_{2} (t)}}{N}{mkern 1mu} intlimits_{0}^{1} {beta _{2} } (a)S(t,a){mkern 1mu} da – gamma _{2} I_{2} (t), frac{{{text{d}}I_{1} (t)}}{{{text{d}}t}} & = frac{{I_{1} (t)}}{N}{mkern 1mu} intlimits_{0}^{1} {beta _{1} } (a)S(t,a){mkern 1mu} da – gamma _{1} I_{1} (t) frac{{{text{d}}R(t)}}{{{text{d}}t}} & = gamma _{1} I_{2} (t) + gamma _{1} I_{2} (t), end{aligned}$$
(4)
The total infected population is (I(t)=I_1(t)+I_2(t)). Figure 4 shows a simulation of this system. Before the onset of the second variant, i.e. for (t< 1000), we observe a peak, followed by a small shoulder and a downward tilted plateau. The second variant corresponds to a higher transmission coefficient: namely, we take here (beta _2(a) equiv frac{3}{2} beta _1(a)). When it appears at time (t=1000), initially there is no effect, because the initial number of infectious with variant 2 is very small. Then, there is an exponential growth caused by this second variant gaining strength. The secondary peak is then higher than the first one. A very small shoulder precedes another stabilisation on a downward plateau.
Figure 4 also shows the dynamics of fractions of infected with each one of the variants. Note that the infectious with variant 1 very rapidly all but disappear at the onset of the second exponential growth phase. One might have expected that the historical strain would be gradually replaced by the new strain, merely tilting further downward the plateau. But that does not happen. Thus, it is remarkable that the historical strain gets nearly wiped out at the very beginning of the second exponential growth.
Multiple epidemic rebounds due to a variant virus: susceptible individuals are divided into 50 discrete groups in the case where a new variant appears at Day (t=1000). The transmission rate (beta _2) is taken as (beta _2(a) = 1.5 , beta _1(a)), (d=0.0002), (gamma _1=0.1) and (gamma _2= 0.05). The fraction of infected individuals in the population is represented in the left graph in logarithmic scale. The total infected population is represented in linear scale in the right graph (black curve), variant 1 in red and variant 2 in green.
Application to the Thau lagoon measurements
Model (1) describes the dynamics of the fraction of infectious in the population, that is (t mapsto I(t)/N). Therefore, we need to derive this fraction from the wastewater measurements. To this end, we use an “effective proportionality coefficient” between the two quantities. This coefficient itself is derived from the measurements (compare Section “SARS-CoV-2 concentration measurement from wastewater with digital PCR” in the “Methods” part below). Calibration of model (1) also requires fitting the values of (gamma), the profiles (a mapsto beta (a)) and the initial distribution of susceptible individuals in terms of a.
We carried this procedure and the resulting fitted curve is displayed in Fig. 5. Note that the outcome correctly captures the shoulder and plateau patterns.
Calibrated model on Sète area: blue dots are measures of SARS-CoV-2 genome units and black curve represents the total infected individuals as an output of the model discretized into (n_g=20) groups in a. Initial distribution of susceptible individuals and (beta) function are taken as described in supplementary information. Parameters d and (gamma) are taken as follows: (d=2.5 times 10^{-4}) (day^{-1}), and (gamma =0.1) (day^{-1}).
The underlying dynamics of the rate of susceptible individuals is given in Fig. 6 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 state.
Calibrated model on Sète WWTP: density of susceptible individuals of each group for (n_g=20). The densities of susceptible of each group is represented in colour curves as functions of time. The curves are ordered from top to bottom according to increasing a. The resulting average total susceptible population is represented in black. Susceptible individuals of highest a trait, which are represented in the bottom light blue curve exhibit a non monotonic behaviour.
Source: Ecology - nature.com
