A physicochemical model for rationalizing SARS-CoV-2 concentration in sewage. Case study: The city of Thessaloniki in Greece.

Detection of SARS-CoV-2 in sewage has been employed by several researchers as an alternative early warning indicator of virus spreading in communities, covering both symptomatic and asymptomatic cases. A factor that can seriously mislead the quantitative measurement of viral copies in sewage is the adsorption of virus fragments onto the highly porous solids suspended in wastewater, making them inaccessible. This depends not only on the available amount of suspended solids, but also on the amount of other dissolved chemicals which may influence the capacity of adsorption. On this account, the present work develops a mathematical framework, at various degrees of spatial complexity, of a physicochemical model that rationalizes the quantitative measurements of total virus fragments in sewage as regards the adsorption of virus onto suspended solids and the effect of dissolved chemicals on it. The city of Thessaloniki in Greece is employed as a convenient case study to determine the values of model variables. The present data indicate the ratio of the specific absorption (UV254/DOC) over the dissolved oxygen (DO) as the parameter with the highest correlation with viral copies. This implies a strong effect on viral inaccessibility in sewage caused (i) by the presence of humic-like substances and (ii) by virus decay due to oxidation and metabolic activity of bacteria. The present results suggest days where many fold corrections in the measurement of viral copies should be applied. As a result, although the detected RNA load in June 2020 is similar to that in April 2020, virus shedding in the city is about 5 times lower in June than in April, in line with the very low SARS-CoV-2 incidence and hospital admissions for COVID-19 in Thessaloniki in June.

Introduction

To assist public health authorities in recommending interventions against the spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) or, conversely, loosening protective measures in large populations, at minimum cost in human lives and the economy, requires reliable large scale early warning methods. Medically confirmed data are by definition not suitable to provide early warning signs because only cases of infected patients are reported, when they are found. In addition, clinical tests, practically targeting small parts of the population, are not sufficient to reveal the true scale of the epidemic and set off alarms, since there are likely many more asymptomatic patients circulating in the community than the identified ones (Ahmed et al., 2020; Randazzo et al., 2020).

Estimating the spread of SARS-CoV-2 through its load in sewage has recently been employed by several researchers due to virus particles and fragments in human excretions, especially in stool, ending up in sewage (e.g., Ahmed et al., 2020; Lodder and de Roda Husman, 2020; La Rosa et al., 2020; Rimoldi et al., 2020; Medema et al., 2020; Randazzo et al., 2020). The notion of fecal-oral transmission route of the virus is supported by findings showing that gastric, duodenal and rectal epithelial cells infected by SARS-CoV-2 release virus to the gastrointestinal tract (Xiao et al., 2020; Wang et al., 2020; Yeo et al., 2020). Wastewater-based epidemiology has been suggested as a promising method for rapid and inexpensive tracking of viruses in the general population thus, allowing more efficient population targeting for clinical diagnostic testing and subsequent prompt advancement of mass protective measures or vaccination campaigns (Daughton, 2020; Mallapaty, 2020; Xagoraraki and O'Brien, 2019, Hellmér et al., 2014).

Several reports, chiefly in online archives yet to be approved by peer review, have presented data of detected coronavirus in wastewater. Even though it has been argued that the virus can be active in sewage for approximately 3 days, its ribonucleic acid (RNA) can still be detected for longer time (Nghiem et al., 2020). Most of these data are based on rather sporadic sampling and are qualitative (presence/absence of virus). Nevertheless, there are also efforts providing quantitative information on the amount of virus found in wastewater (e.g., Ahmed et al., 2020). Exploitation of such data for accurate estimation of infected people in the community still presents a major challenge as it depends on reliable knowledge of person-to-person variability of virus shedding in excretions like sputum, urine, sweat and feces, and gathering such knowledge is by itself a challenging task (Wölfel et al., 2020).

Despite the never-ending need to further increase the accuracy in wastewater molecular analysis, the present level of the RT-qPCR technology appears to provide quite representative estimation of the virus genome concentration in sewage (Randazzo et al., 2020). However, data extrapolation for exact assessment of virus load in a community, is difficult. This is because the coronavirus, like all viruses, has very small size (order of a hundred of nanometers) is extremely surface active (crucial for its survival and spread) and so interacts intensively with its environment (Kuzyakov and Mason-Jones, 2018). Hence, it can penetrate and adsorb strongly onto the highly porous suspended solids (feces, sand, silt, clay, etc) in sewage and remains undetected with current methodologies. As a result, only a small fraction of viral load can be actually detected in sewage samples despite the optimized protocols used in samples preparation to concentrate and recover it from wastewater (Ahmed et al., 2020; La Rosa et al., 2020; Randazzo et al., 2020). A rough average value of suspended solids in sewage is about 1 g/L, this depending enormously on the employed sewer system, (e.g., combined sewer system or separate-storm sewer system). For typical suspended clay particles accommodating a vast micro- and nano-pores network the available surface area is of the order of several hundreds of m2/g (Macht et al., 2011). Small viral particles and empty viral capsids can move into such pores, be adsorbed and so get inaccessible (Kuzyakov and Mason-Jones, 2018). For that reason, virus extraction from wastewaters demonstrates limited and rather specific recovery rates, even when PEG flocculation methods are applied (Ahmed et al., 2020).

Many studies focus on the development of experimental protocols aiming to improve the way of extracting more viral material from wastewater and its suspended solid matter. Apart from these continuous efforts, another way to overcome the above difficulty is by developing appropriate physicochemical models that can account for the strong adsorption of virus particles and fragments onto the surface of pores. Such models should start from the measured small fraction of virus concentration in the sewage and proceed by rationalizing the measurement by taking into account the flow rate of sewage, the suspended solids concentration, the mechanisms of virus adsorption onto the pores surface and the chemical and biological parameters that affect such adsorption.

There is enough evidence in literature that salinity and organic materials widespread in natural soil and aquatic environments determine virus retention, transport and inactivation. For instance, it has been reported that dissolved (liquid-phase) organic matter, e.g., low MW humic acids, competes for adsorption sites with viruses. In addition, bonded organic matter, on one hand occupies sorption sites, thus reducing the available sites for the virus, but on the other hand provides new hydrophobic binding sites. All in all, adsorption of viruses in soil and aquatic environments is a very complex process, which depends greatly on parameters like the form of organic matter (mineral-associated, dissolved or bonded), salinity, and the type of virus. These parameters act in parallel and depending on their relative intensity (e.g., concentration) can dictate different interaction (electrostatic versus hydrophilic) mechanisms (Zhuang and Jin, 2003; Zhao et al., 2008; Cao et al., 2010; Wong et al., 2013). In other studies, it was found that removal and inactivation of viruses in different granular media is non-linear, reflecting chiefly the heterogeneity within soil grains and differences in the physical and chemical structure of viral particles (Betancourt et al., 2019; Schijven and Hassanizadeh, 2000). Overall, the impact of environmental parameters on virus persistence in the environment has received little attention (Waldman et al., 2020). Yet, there is some recent experience on the influence of environmental parameters against the action of biocides (due to adsorption, galvanic deposition, etc) in recycled hygiene and potable water samples under conditions prevailing in the International Space Station (Petala et al., 2017. Giakisikli et al., 2017; Mintsouli et al., 2018; Petala et al., 2020).

This work is about developing a mathematical framework, at various degrees of spatial complexity, of a physicochemical model that can rationalize the quantitative measurements of viral genome (whole non-infectious particles and fragments) in sewage with regards to the effect of environmental parameters on the adsorption of virus onto suspended solids. In order to determine the values of the adsorption model variables, measurements of viral genome concentration and environmental parameters are required over a short period of time (few days) during which virus shedding at the source can be considered as approximately constant with respect to virus viability in humans and the averaging over a large population.

Greece is one of the few European countries that harnessed the spreading of the virus to low levels (National Public Health Organization GR, 2020). This has been attributed to the early implementation of strict quarantine measures, being followed by a gradual release of measures with social distancing, sanitation and personal health protection means (face masks, hand hygiene etc.) still in place. Thessaloniki, a city of around 1 million inhabitants, has experienced an appreciable number of SARS-CoV-2 infections and also many hospital admissions for coronavirus disease 19 (COVID-19) with a peak in April 2020. These numbers became much smaller in June 2020 when the epidemic subsided. This low level of infection in June rate permits the hypothesis that at the last stages of the population infection curve during the quarantine and before the virus becomes eventually non-detectable, virus shedding at the source (wastewater produced by individual households in the city) has been roughly constant for about 10 days. This hypothesis is essential for determining the functional dependence of the adsorption model parameters on the measured physicochemical properties of sewage over a period of time where environmental parameters usually vary.

In the following sections, the collection, and pre-treatment of samples is presented first. Then, the procedure for physicochemical characterization of samples and the molecular analysis for detecting and quantifying SARS-CoV-2 RNA is described. Details of model development follow. Finally, a set of physicochemical and molecular measurements obtained over a certain period of time in the sewage of Thessaloniki is employed to determine the model parameters and then rationalize the measured values of the concentration of virus at the source (in the city) for a subsequent period of time.

Materials and methods

Infected population in Thessaloniki and wastewater samples' collection

Wastewater samples were collected at the entrance of the Thessaloniki's Wastewater Treatment Plant. Samples were collected from a channel where the screened inflow from the city merges with recirculating overflows from sludge treatment units. The flow rate of this recirculating stream is about 25 times smaller than the flow rate of fresh sewage coming from the city. Yet, it has a very high suspended solids and organic matter load that contributes to the overall quality characteristics of influent wastewater. This is an advantage in the present study because both suspended solids and organic matter load enhance the extent of adsorption and so allow more accurate determination of the adsorption model parameters (see below).

The history of measures of the Greek State against the pandemic is as follows. On March 10th schools and universities were closed whereas on March 16th most shops, museums, entertainment and sports places (gyms, athletic centers, etc) were shut down. On March 23rd strict quarantine measures in people's circulation were applied to the entire population. The strict quarantine for people's circulation ended on May 4th. Progressively after that date, small groups of activities and businesses were allowed to resume with sports places being the last where restrictions were released from June 15th onwards. National Public Health Organization (EODY) recorded 51 SARS-CoV-2 infected Thessaloniki prefecture inhabitants in April 2020 and only 3 infected inhabitants in June 2020 (EODY data records). Admission data of AHEPA hospital, the COVID-19 reference hospital in Thessaloniki at that time, showed an appreciable number of patients in April 2020 (28 hospital admissions due to COVID-19) and a steadily declining number thereafter with only 1 hospital admission in June 2020.

Samples were collected three times per week as of April 21st, 2020. Inlet flowrate was registered at the same dates. Samples were obtained using a refrigerated autosampler (6712 Teledyne ISCO) programmed to deliver a 24-h composite sample by mixing consecutive half-hour samples. Samples were transported to the lab on ice. Part of the samples was processed immediately whereas another part was refrigerated at 4 °C until further processing within 24 h.

Chemical analyses and qualitative detection of SARS-CoV-2 RNA

Chemical analysis of the samples was performed by standard lab equipment to determine certain physical and chemical parameters that are suspected to play a role in the adsorption of virus particles or fragments onto suspended solids. These parameters are: Total Suspended Solids (TSS), pH, electrical conductivity/salinity, dissolved oxygen (DO), 5-days Biochemical Oxygen Demand (BOD5), Chemical Oxygen Demand (COD), UV absorbance at 254 nm (UV254), Dissolved Organic Carbon (DOC), Specific UV absorbance (UV254/DOC), ammonium, nitrates, total organic nitrogen, total phosphorus and ortho-phosphates. These parameters were measured according to Standard Methods (APHA et al., 2012).

From the above, the overall organic loading of the samples is described by BOD5, COD, DOC and indirectly by DO, as DO is associated with the activity of cells and biomass inside the sewage piping systems (considering that wastewater residence time along the sewerage network varies from a couple of hours to 24 h). Humic-like substances are depicted by UV254 and Specific UV absorbance, UV254/DOC, represents the fraction of humic-like substances in the total dissolved organic carbon. Ionic strength (in terms of tendency for electrons transfer) is described by electrical conductivity/salinity, pH and redox potential. Nutrients concentration are represented by the concentration of ammonium, nitrates, total organic nitrogen, total phosphorus and ortho-phosphates.

RT-qPCR assays (Thermo Scientific) for the qualitative detection and characterization of SARS-CoV-2 RNA were used in this study. A small portion of the sample, 50 mL, was centrifuged (3000 g, 4 °C, 1 h) to eliminate solids, and the supernatants were used for further processing. In total 40 ml of the solid-eliminated samples were concentrated through filtration using a Centricon® Plus-70 centrifugal ultrafilter with a cut-off of 10 kDa (Millipore, Amsterdam, the Netherlands). For concentration, samples were centrifuged at 3000 ×g for 3 h, using a slightly modified method reported by Medema et al., 2020. Concentrated samples at final volume of up to 300 μL were used for RNA isolation with the Trizol LS reagent (Invitrogen, Cat No 10296028), according to the manufacturer's instructions. To enhance RNA yield, 2 μL of RNA carrier (FG, Carrier RNA, ThermoScientific LSG, Cat No 4382878) were added to each sample before RNA extraction. RNA pellets were resuspended in 50 μL RNAse-free water and the yield of the preparations was estimated spectrophotometrically. The presence of viral genomic RNA was determined through one-step RT-qPCR reactions, using the SuperScript III Platinum One-Step RT-qPCR kit (ThermoScientific LSG, Cat No 11732–088) and TAQMAN chemistry for the detection of viral RNA (TAQMAN 2019 NCOV ASSAY KIT V1, ThermoScientific LSG, Cat No A47532).Three assays that target SARS-CoV-2 genes (N, S, and orf1ab- nucleocapsid, spike protein and open reading frame 1ab, respectively), and one positive control assay that targets the Human RNase P RPPH1 gene, were used. Targeting three different viral genomic regions reduced the risk of false negatives. Assays have undergone bioinformatic selection and analysis to specifically target sequences that are unique to SARS-CoV-2. For each reaction, 5 μL of the extracted RNA (1/10th of the total preparation) was used as template. Each experimental procedure included positive control reactions (TAQMAN 2019-NCOV CONTROL KIT V1, ThermoScientific LSG, Cat No A47533), corresponding to synthetic DNA containing target sequences for each of the assays included in the TAQMAN kit and for the Human RNAse P RPPH1 gene. Non template controls and internal controls (Human RNAse P RPPH1) were also run on each plate of analysis (StepOne Applied Biosystems, Thermo Scientific). Each sample was analyzed in triplicate. RT-qPCR results were interpreted as following: each viral genome assay (orf1ab, S, N) with a Ct value lower or equal to 37 was considered positive. When at least two of the tested targets were evaluated as positive, the corresponding sample was considered positive to SARS-CoV-2.

For SARS-CoV-2 quantification, standard curves generated using 10-fold dilutions (typically from 10 to 106 copies) of a reference material were used. Synthetic DNA (TAQMAN 2019-NCOV CONTROL KIT V1, ThermoScientific LSG, Cat No A47533) containing 1 × 104 copies/ μL of the target sequences for each of the tested assays and for the Human RNAse P RPPH1 gene was used for standard curve preparation. Standard curves were created for each target gene based on known copy numbers and corresponding RT-qPCR results (Ct values) and displayed linear behavior (R2 ranging from 0.88 to 0.99).

For sewage samples' viral load quantification (copy number), extrapolation of corresponding RT-qPCR results (Ct values) to the generated standard curves was used. To determine the viral load per ml of sample, the portion of the tested RNA preparation used (5 μL out of 50 μL) and the original volume of the sample subjected to concentration (40 mL) were taken into account. Viral load estimations were done both independently for each target gene and by averaging the results acquired for each target gene. Data shown in Fig. 5 were determined based on the standard curve for the N target (characterized by the best R2 value, R2 = 0.99).

Fig. 5

SARS-CoV-2 viral copy numbers per mL based on the N protein and Orf1ab standard curves (a), and evolution of physical-chemical parameters throughout the study period: inlet flowrate and total suspended solids concentration (b), dissolved oxygen concentration and specific UV absorbance (c) and electrical conductivity (d).

Method of data analysis

Physicochemical model

The model considers the dynamic of adsorption of virus parts on suspended solids along the complex flow path in a city sewerage system. The physicochemical model (being at the heart of the modeling effort) is designed as the simplest possible that can describe the problem at hand. Nevertheless, several alternatives of variable degree of complexity are presented regarding the topology of the model.

Firstly, the state variables of the physicochemical model are presented. Viral gene copy concentration in the liquid phase is denoted as C, the mass concentration of solids is denoted as m and the adsorbed mass of virus parts per unit mass of solids is denoted as q. The adsorption process is completely described by two equations: The first one is the so-called adsorption isotherm which correlates the values of q and C at equilibrium conditions. Any deviation from this equilibrium creates a driving force for adsorption (or desorption depending on the direction of the deviation). The adsorption isotherm is denoted as (only the concentration dependence appears here) (Ruthven, 1984):

The adsorption kinetics are given from the equation (Tien, 1994):

Eq. (2) is based on first principles. This is unlike the typical pseudo-first and pseudo-second order adsorption kinetic models which are purely empirical and cannot be scaled up to different conditions from which the corresponding parameters should be derived (Ho, 2006). Eq. (2) is the mathematical reduction of a complicated problem that considers convection and diffusion in the exterior of a solid particle and adsorption/desorption, pore diffusion and surface diffusion in its pore structure, Fig. 1 . It is noted that hA is a generalized mass transfer coefficient (units of inverse time) which includes the solids surface area and depends on the external mass transfer to the particle (convective diffusion in turbulent flow field) and to internal mass transfer in the pores (pore and surface diffusion) (Kyzas and Kostoglou, 2014).

Fig. 1

Schematic representation of possible interaction mechanisms of virus particles/fragments with a porous solid particle suspended in sewage. 1: mass transfer (convection and diffusion) from the surrounding bulk liquid, 2: adsorption, 3: desorption, 4: pore diffusion and 5: surface diffusion.

Topological models

The state variables C, m and q of the physicochemical model are, in general, functions of time and spatial location in the city piping system. In order to define the spatial location, a topological description of the piping system is needed. There is a hierarchy of this description according to its degree of detail. Here the two limits of full resolution and completely continuous description are presented. It should be noted that any intermediate situation is possible depending on the available computational resources and on the degree of required spatial details.

Full resolution discrete model

In this limit any pipe section between two junctions can be considered as an independent model unit (unit cell) (Fig. 2). The unit cell is characterized by its length L, flow cross sectional area S, fluid velocity V (computed by dividing the volumetric flow rate by S), and unit cell inlet conditions with respect to viral particles/fragments (bulk and adsorbed) and mass concentration of solids. It is important to realize that a single value of q is not enough to describe the adsorbed virus parts since there are solids introduced at several locations in the flow system having different residence times and exposure to different values of concentration C. This problem is overcome by allowing q to be non-uniformly distributed among the solids. Mathematically, the function f(q) is introduced defined as the mass concentration of solids having on them a weight fraction of viral particles/fragments between q and q + dq. The following condition holds .

Fig. 2

Full resolution discrete topological model accounting for each pipe section and junction in the sewage piping network. In the schematic of the unit pipe (up-right corner) individual virus parts and solid particles as well as solid particles with adsorbed viral particles/fragments are shown to flow in sewage.

The governing equations for C(z) and f(q,z) along a unit cell (distance z from entrance/junction) are:

The mass balances in each junction have the form

Index i refers to the streams met in the junction, the velocities are signed and Eq. (5c) must hold for all q values.

The above system is quite descriptive but rather computationally intractable since a system of several algebraic equations, one non-linear partial differential equation and one ordinary integrodifferential equation must be solved for each unit cell (considering that the piping system of a city can contain thousands of unit cells).

Fully continuous limit

In the limit of slow adsorption in the local residence time scale, a procedure of transforming the discrete problem to continuous can be followed. Here the fully continuous limit is presented (Fig. 3). Only the flow in a central pipe is considered with all the information from the other piping systems entering in the form of continuous variables along the main flow path.

Fig. 3

Fully continuous limit topological model accounting for lumped regions in the sewage piping network.

Let z be the length along the central pipe. The rate of liquid mass, solids volume and virus parts VP mass entering the central pipe per unit length at position z is described by the functions Fin(z), min(z), Pin(z). The governing equations are: where δ is the Dirac delta function and accounts for the fact that solids enter the path with no adsorbed virus parts on them. Apparently, this assumption can be relaxed in case of having more specific information on this quantity.

Generalizations

It is noted that any combination of the discrete and continuous approaches is possible. For example the main flow paths from lumped small regions to the central flow can be treated explicitly in addition to the central pipe, Fig. 4 . The choice is based on the topology knowledge, the degree of information needed and the available computational tools. An additional variable that can guide the topological simplification procedure and the sharing among discrete and continuity description is the ratio of residence to characteristic adsorption time. In case of small values of this ratio the corresponding piping structures can be replaced by continuous ones.

Fig. 4

Combination of the discrete and continuous approaches is possible. For example the main flow paths from lumped small regions to the central flow can be treated explicitly in addition to the central pipe.

Another very important issue is that the parameters of the function I(C) are not constant, but may have a complex dependence on the local physicochemical environment at the specific location along the flow path (e.g. organic matter, salinity, temperature, dissolved oxygen etc). These quantities, denoted by the vector x, should have prescribed values for each unit cell in the discrete model. In case of the continuous model a prescribed function x(z) is needed.

The limit of fast adsorption equilibria

A particularly interesting limit is the one of the adsorption time being very small with respect to the residence time for every flow element in the system (time from entrance source to measurement point). In this case, adsorption equilibria are established, and the mathematical model is defined only at the measurement point i.e. only the entrance source and outflow conditions affect the results. When this holds, all the above models, regardless their topological structure, degenerate to the following one (for completeness a finite volume fraction of solids is handled):

A mass balance at the measurement point leads to mq + (1-φ)C = constant = (1-φ)Co where Co is the inlet virus parts VP concentration in the entrance source and φ = m/ρ is the solids volume fraction. The dominant equation takes the form (Langmuir isotherm is considered according to which q = QmKC/(1 + KC) where Qm is the corresponding maximum value and K is a parameter):

In the particular case of small concentrations of virus parts VP, the Langmuir isotherm can be approximated by a linear one and the above equation is simplified and solved analytically to give:

The two unknown parameters K and Qm can be unified to one as A = KQm. In case of small values of φ its contribution to the result is relatively small so only a rough guess for ρ is adequate for further computations. The new working equation is

As discussed before, in principle A depends on a vector x of variables at the measurement point (organic matter, salinity, DO etc). This is actually a fast simplification arising at the fast equilibrium limit. The state variables at the measurement point depends only at the local value of x at this point and not to the complete history of x along the flow path as it is the general case. It follows how the relation (11) can be exploited in practice.

Let us say that at two different days the measurements for C and m, x are C1, C2, m1, m2 and x 1, x 2 respectively. Assuming a similar initial viral particles/fragments VP concentration, Co, leads to the relation

In this equation the only unknowns are the values of A. In case of N measurements, N-1 relations of the above type are derived. The scope is to estimate the multivariable function A(x) from N-1 scattered measurements. There is no general approach to this and the procedure which will be followed should be based on the particular actual values of vectors x i.

When assuming a similar shedding rate of viral particles/fragments (Ro) among the sampling days, instead of concentration, the procedure has to be changed as follows: Eq. (11) is tautotically equivalent to:where R is the mass flow rate of virus parts at the measurement point and the subscript “o” denotes again the inlet virus parts at the entrance source. Let us denote F the total flow rate. The equation corresponding to Eq. (12) but for constant Ro is

The aim is (i) to estimate the function A(x) from data and (ii) to use it in order to relate the measured quantities to Ro. This will be done as follows: Let us say that we have found an estimation of function A(x). The measurements at two different days are C1,x1,m1,F1 and C2,x2, m2, F2, respectively. The scope here is to identify the meaning of these measurements for the spread of the disease in the community. This means that we want to find the ratio of the viral particles/fragments shedding rates Ro1 and Ro2. This is given as

So selecting a particular day as a reference, one can check the relative variation of the mass of produced viral particles/fragments.

Results

Biological and physicochemical measurements

Fig. 5a-d presents SARS-CoV-2 gene copies number per mL of the tested wastewater samples (n = 29) based on the N protein standard curve along with various physicochemical parameters measured in the collected wastewater samples. The detectable viral copies in the samples collected between April 21st and May 3rd exhibit an overall declining trend along the sampling days. For the rest of May 2020, the virus genome targets are undetectable in the analyzed samples. The latter is in line with the high number of SARS-CoV-2 infections in Thessaloniki prefecture (n = 51) and hospital admissions (n = 28) in April 2020 and the very low number of infections (n = 3) and admissions (n = 1) in AHEPA hospital in Thessaloniki during strict quarantine in June 2020 (National Public Health Organization GR, 2020). Yet, as of June 3rd viral copies are detected again at an approximately constant level until June 22nd, which is the last day of measurements in this study.

The sample-by-sample variability of all physicochemical parameters, Fig. 5b-d, is rather usual for the urban wastewater of Thessaloniki which utilizes in part a combined sewer system. The concentration of TSS is high due to leachates recirculation at the sampling point. For the same reason, the organic load is also high. The observed small increase of humic-like substances over the samples might influence the detection of viral copies and this is explored further below. Thessaloniki has since its Roman times a well-documented rich aquifer which was utilized continually, until the early 20th century, when it has been abandoned as a water resource. However, the western part of the urban area, where the sewage treatment plant is, lies at or below sea level, leading to interference of saline water with the sewage, which explains to the unexpected high values of electrical conductivity (Stavridis et al., 2017).

Model application

The physicochemical interactions in adsorption are too complicated so no quantitative information about the form of the function A(x) can be theoretically invoked to facilitate the parameter estimation process. On the other hand, some qualitative information can be extracted in order to support the quantitative analysis. For example, it is expected that salinity (i.e. conductivity) mainly affects the adsorption process through the adjustment of electrical double layer interactions. However, it is well known that for the measured high salinity levels the double layer is already quite compressed so salinity is not expected to influence A(x).

Due to the unknown A(x) structure, the problem is undetermined (it is non-linear having N-1 relations and N unknowns). Clearly the virus shedding rate Ro is required to close the problem. However, in the specific case of the majority of viral particles/fragments being in adsorbed state (Kuzyakov and Mason-Jones, 2018) the problem can be linearized overcoming the requirement for knowledge of an additional variable. In this case, Eq. (14) takes the form:

This relation can be applied to the actual data for N days to find the relative values of A(x). Without loss of generality, the first day can be considered as the reference one with A1 = 1. So a vector of length N with values Ai (i = 1,2,3…,N) is generated. As explained in the Introduction, the time period of the N measurements (days) must by definition be much smaller than the viability of the virus in humans, to permit the hypothesis of an approximately constant virus shedding rate, and large enough to give meaningful results. So, a typical range of values N between 6 and 10 must be considered. In the present case, the choice is N = 6 (dictated by the sampling performed every two days). The best possible choice would have been N = 10 corresponding to everyday measurements.

The problem here is that a function of M parameters (M is the length of vector x) has to be fitted using N data points. The number M can be initially reduced on a physical basis. Apart from removing the conductivity/salinity, the ambient temperature can be also removed from x since the temperature in sewage pipes is pretty constant within a period of a few days -because of the thermal inertia of the ground- and independent from the ambient temperature.

As a systematic procedure, one must first check if a correlation between A and each measured environmental parameter exists (which henceforth will be termed as variable to ascribe better the algebraic sense of variability). If this is the case, only a single variable fully determines adsorption. If, however, this is not the case linear expansions must be fitted to the data points for sets up to N-1 variables. Then a combination of a single quadratic and N-3 linear expansions must be tested. In this respect, all the possible combinations of variables (per two, per three etc) up to total number M must be tested employing all possible combinations of polynomial expansions. It is stressed that for a given confidence level (based on the estimated accuracy of the data) the simplest fitting (in terms of the number of independent variables) that achieves this accuracy level must be selected. Increasing the number of variables is only adding noise to the problem.

In case of no acceptable results from the above systematic procedure, a non-linear generalized fitting for pairs or triplets of variables must be pursued. There is no rigorous way to do this so there are elements of art in this procedure and it is important how much one insists in achieving the goal before arguing that there is no intrinsic correlation of A to x.

Let us apply the above procedure to the present data having N = 6. The vector of Ais takes the values 1, 1.05, 0.6, 0.99, 1.38, 1.32 (the units can be ignored since only the relative values are essential). There are two important observations: a) for a confidence interval of 5% there are actually three levels of A with relative values 1, 0.6 and 1.32. The differences are significant and cannot be characterized as random noise or experimental error. b) The fact that the value of A is reduced by 40% from one sample to the next confirms the significance of the present approach. Given the several days life expectancy of virus in the human body and the averaging performed over the city's population, the observed sudden decrease of A value in just two days (sampling interval), during which the virus production rate is assumed approximately constant, confirms that adsorption is heavily affected by the composition of the liquid.

The fitting procedure is performed for M = 4: specific absorption (UV254/DOC), COD, Conductivity, Dissolved Oxygen. Conductivity is included just to test the theoretical argument of no influence. The systematic procedure has not led to acceptable results so the generalized non-linear fit for pairs of independent variables is pursued. The only pair that leads to reasonably correlated results is the one of specific absorption (SA) and dissolved oxygen (DO). The fitting correlation is.

For the most part, the adsorption of virus onto suspended solid matter increases with specific absorption and decreases with dissolved oxygen. An increased specific absorption reflects the higher amount of humic-like substances, represented by UV254 values, adsorbed onto the solids surface. These substances henceforth provide new binding sites for viruses. On the other hand, dissolved oxygen has been suggested to cause not only virus inactivation, but also capsid oxidation and further disassociation (Pinon and Vialette, 2018). What's more, sufficient levels of dissolved oxygen increase the metabolic activity of bacteria and bacterial enzymes, factors that are linked to viral destruction. Decayed or damaged viruses may lose their capacity to get adsorbed onto the solids surface or interact with humic-like substances leading eventually to a lower extent of virus adsorption. It is stressed that empirical correlations like Eq. (17) can be used only for the range of variable values for which it is derived (i.e. no extrapolation is possible).

The above expression can be used to rationalize (so, make them directly comparable) the measurements of C, F, m, SA, DO for two different days (denoted by subscripts “1” and “2”, respectively) giving the ratio of the actual viral particles/fragments shedding rates.

It is noted that the above analysis holds in the case that the total viral particles/fragments shedding rate in the city is of interest. In the case that information on the spatial distribution of the production rate at different sources across the city is needed, the fast adsorption (equilibrium) hypothesis may not be valid (since local residence times must be set to the corresponding criterion) so the complete kinetic model developed here must be solved. Apparently, many measurements at different locations along the sewerage piping system are needed in order to set-up a very complex inverse engineering problem for the unknown sources.

With the values of the A(x) variables at hand, it is interesting to apply the aforementioned procedure to the data collected for the period after June 3rd (between May 5th and June 3rd viral RNA copies are not detected in the samples). Fig. 6 displays the corrected by the model evolution of the viral particles/fragments shedding rate (in the city) after June 3rd, R, with respect to the reference value, Ro, for the period from April 21st to May 3rd. Apparently, there is a small increase (R/Ro ~ 1.2) in viral RNA copies shedding rate in the first days of June most probably as a result of the crowding of young people at beaches and bars after the end of the quarantine on May 4th. However, after June 10th the viral particles/fragments shedding rate decreases towards a low range of values approximately between 0.2 and 0.3. This reduction is largely attributed to the low amount of suspended solids after June 10th (Fig. 5b), which drives the model to lessen the estimated shedding rate because virus inaccessibility due to adsorption is limited. This is clearly in line with the calm conditions at the city hospitals during the whole June 2020. This is clearly in line with the very low number of COVID-19 patients admitted to AHEPA hospital in Thessaloniki during June 2020 (1 patient) as compared to 28 patients in April 2020.

Fig. 6

Evolution of relative shedding rate of viral RNA copies after June 3rd, 2020, as corrected by the present model. The overall low values after June 10th, is in line with the calm conditions observed in city hospitals during this period.

Restrictions on international flights to most Greek airports were released on July 1st (with a few exceptions from countries where the pandemic is at a peak) and the effect of travelers from abroad on virus spreading in the community remains to be seen.

It will take extensive testing with data from different sources and with different viruses to certify the success of the present model. However, even at this point the model is still useful paving the way for proper rationalization of measurements that will transform sewage surveillance from just a useful indication to hard evidence and, eventually, to a reliable tool for public health monitoring. The most prominent advantage of monitoring pathogens in public wastewater is the detection of asymptomatic carriers within the healthy population virus shedding, providing to the authorities a robust early diagnosis indicator that will allow applying or withdrawing measures to prevent spreading. Let's not forget that wastewater-based epidemiology is tremendously cheaper and faster than clinical screening although, admittedly, it cannot replace it.

Conclusions

A mathematical framework is developed for the rationalization of quantitative measurements of SARS-CoV-2 RNA copies in sewage with respect to wastewater environmental parameters at a varying level of spatial distribution, e.g., samples obtained at different spots across a city, across a suburb or just from a single neighborhood. The simplest case is treated herein where sewage samples from just one spot is employed. The heart of the framework is a physicochemical model which incorporates established knowledge on possible effects of environmental parameters such as the concentration of suspended solids and the amount/type of dissolved organic matter on the adsorption of viruses on soil and other solids. The model is applied to Thessaloniki, a city of about one million inhabitants in Greece, which represents a convenient case study because of the observed low level of infection rate along the course of several days which allows making the hypothesis of constant virus shedding during these days. By performing a non-linear generalized fitting procedure for pairs or triplets of variables (i.e. measured environmental parameters) against the sample-by-sample variation of the measured viral copies concentration, an effective functional dependence of the linearized Langmuir adsorption isotherm parameter, A(x), of the model is achieved. The results of this procedure indicate the ratio of the specific absorption (UV254/DOC) over the dissolved oxygen (DO) as the parameter with the highest correlation with the viral copies. This implies a strong effect on viral inaccessibility in sewage caused by the presence of humic-like substances combined with virus decay induced by oxidation and increased metabolic activity of bacteria and bacterial enzymes. More work over broader sewage conditions and at different locations is needed before conclusive statements can be made. Even so, applying the model to the measured viral copies concentrations in Thessaloniki sewage after June 3rd 2020 leads to an overall low range of values of virus production rate, in agreement with the calm conditions in the city hospitals this period.

CRediT authorship contribution statement

M. Petala: Methodology, Conceptualization, Investigation, Data Curation, Formal analysis, Visualization, Writing - Original Draft, Writing - Review & Editing, D. Dafou: Methodology, Conceptualization Investigation, Data Curation, Formal analysis, Writing - Original Draft, Writing - Review & Editing, M. Kostoglou: Methodology, Conceptualization, Data Curation, Visualization, Software, Validation, Formal analysis, Writing - Original Draft & Editing, Th. Karapantsios: Supervision, Project administration, Funding acquisition, Methodology, Conceptualization, Writing - Original Draft, Writing - Review & Editing, E. Kanata: Investigation, A. Chatziefstathiou: Investigation, F. Sakaveli: Investigation, K. Kotoulas: Resources, M. Arsenakis: Methodology, Conceptualization, Writing - Original Draft, Writing - Review & Editing, E. Roilides: Conceptualization, Data curation, Writing - Original Draft, Writing - Review & Editing, T. Sklaviadis: Methodology, Conceptualization, Writing - Original Draft, Writing - Review & Editing, S. Metallidis: Resources, A. Papa: Resources, E. Stylianidis: Funding acquisition, Conceptualization, A. Papadopoulos: Resources, Writing - Original Draft, Conceptualization, N. Papaioannou: Supervision, Funding acquisition, Methodology, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.