There have been major advances recently on the therapeutic approaches of hepatitis B virus (HBV)-related hepatocellular carcinoma (HCC). Surgical treatments are the key curative treatments of HCC, whereas local ablative treatments may also achieve clinical remission in selected cases. Trans-arterial locoregional therapies are regarded as palliative but still lead to improved survival. There have been major breakthroughs in the systemic therapies for HCC. The first marketed targeted therapy, sorafenib, was shown to improve survival in patients with advanced HCC. Studies on other targeted therapies also showed promising results. Suppressing HBV with effective antiviral treatment would also benefit HCC patients by reducing recurrence and improving liver function.
1.
Introduction
The COVID-19 pandemic has exerted huge and unprecedented pressure on public health resources globally. Cross-sectional surveys to establish disease prevalence are likely to be financially unsustainable in the long term and rely heavily on continued cooperation from the public [1]. Wastewater monitoring to detect and quantify SARS-CoV-2 viral RNA shed by infected individuals in the population, and to indicate infection prevalence, was adopted relatively early in the course of the pandemic across a number of countries [2,3], expanding to 67 countries by mid-2022 [4]. Although the demographic coverage and utility of wastewater monitoring varies across adopters of this approach, the method is generally less intrusive, relatively unbiased in terms of its demographic and epidemiological coverage, and costs significantly less per capita than clinical testing programmes (e.g. fanxiexian_myfh148 per individual PCR test cf. fanxiexian_myfh300 per wastewater sample representing larger populations [5]). Wastewater surveillance is thus, arguably, an alternative, or at least complementary, approach to clinical testing programmes.
Wastewater-based epidemiology has been used in some areas of public health for decades [6], but is a relatively novel tool for emerging pathogens. Applications include tracking viral dynamics to monitoring chemical exposures and prescription drug consumption [7,8]. As wastewater sampling for the detection and monitoring of SARS-CoV-2 has been developed and applied at an unprecedented pace, uncertainty remains when interpreting the measured viral RNA signals and their spatiotemporal variation. Variation in the underlying sample, sampling method, and testing, due in part to lack of standardisation, as well as systematic variability in space and time in the measurement environment (e.g. sewersheds), can result in a large degree of noise in the observed signal [9]. Sampling frequency is typically dependent on cost constraints, resulting in sparse and irregularly sampled data. Furthermore, wastewater measurements are typically left-censored if they fall below certain analytical thresholds, such as the limit of detection (LOD), the lowest concentration at which viral RNA is detectable with a given probability (typically 95%); and the limit of quantification (LOQ), the lowest concentration at which viral RNA can be reliably measured with a predefined accuracy. Methods to handle measurements that fall below these limits (e.g. statistical methods, imputation, and scalar or zero replacement) are not standardised and depend on the interpretation of the data (for example, if low values do not impact interpretation they may be omitted from downstream analysis), and information available to the analysts [10,11].
There are several approaches to infer wastewater concentration from noisy, censored, incomplete time series measurements. One common and straightforward approach for denoising time series is to calculate a moving average (MA). However, MA are sensitive to outliers and missing values. There is ambiguity about the most appropriate window-size, and whether to calculate a weighted or ordinary average. Uncentred MAs also operate with a lag, where larger windows create larger lags, delaying reactivity of surveillance in time-dependent operations. Lastly, an MA estimate can never be smaller than the censoring threshold, which leads to biased estimates.
State-space methods model observed data as functions of latent, unobserved stochastic processes and can better account for missing data, observational noise, and censoring. Recently, others have proposed state-space methods to infer viral concentrations from wastewater time series. The underlying "true" viral concentration at time $ X_t $ is modelled as a first-order auto-regressive (AR1) process [12,13]. To account for measurement noise and outliers in the observations, measurements $ Y_t $ are assumed to be equal to $ X_t $ plus an independent mean-zero Gaussian observation error. To account for outliers, from time to time $ Y_t $ is assumed to be replaced by an independent Uniformly distributed random variable that is unrelated to $ X_t $. Left-censoring is accounted for by capping $ Y_t $ at the (known) limit of quantification. Using a Kalman Filter and numerical approximations, the state variable $ X_t $ is inferred from observation data $ Y_t $ to produce a smoothed estimate of viral concentration, with outliers removed, that can extend below the known limit of quantification [12].
In this paper, we propose and test a simpler, more realistic, and more flexible state-space model. Our latent variable is modelled by a first order random walk (RW1) instead of an AR1 process, which reduces the number of model parameters. Instead of randomly replacing observations by a random number, our model generates outliers by assuming observation errors from a heavy-tailed t-distribution. This has the benefit that observations classified as "outliers" can still be informative about viral concentrations.
Our model is implemented in the Stan modelling language [14], which allows for fast Bayesian inference and straightforward extensions of the model.
2.
Methods
2.1. Wastewater data source
Untreated influent samples were collected from sewage treatment plant sites across England at a frequency of four days a week by the Environmental Monitoring for Health Protection (EMHP) programme, led by the UKHSA. The sampling strategy provides coverage of approximately 40 million people across England. Samples were analysed for SARS-CoV-2 RNA by quantifying the number of copies of the nucleocapsid gene (N1) using RT-qPCR. Concentrations under the limit of detection were assigned a value of -4, to be handled during the data processing pipeline depending on the use case. Only sites sampled 30 times or more (around seven weeks' worth of data) were included; median sample count across sites was 145, ranging from 31 to 323.
2.2. Flow normalisation
Extraneous sources of flow, such as heavy rainfall, snow melt, or groundwater ingress into sewers, may dilute wastewater and impact estimates of SARS-CoV-2 RNA concentration. Studies have indicated that the effect of dilution in most cases are minor, but in periods of high dilution events, normalisation is critical [15]. The normalisation approach applied by the English wastewater surveillance programme mitigates this by adjusting measured SARS-CoV-2 concentrations to consider flow. The model is based on the assumption that flow $ F_t $ at time $ t $ is not directly observable. Instead, information about flow is obtained by observing the correlation of concentrations $ \rho_{ti} $ of different markers $ i $ (orthophosphate and ammonia nitrogen). The model assumes:
where λ2 is the flow variance, xti is the load of marker $ i $ at time $ t $, μi and σ2i are the mean and variance of the load of marker $ i $ (all in log space). ⟨logFt⟩ is fixed at 0 to identify the model. Using multiple markers jointly to estimate flow variability can improve the accuracy of estimates [9,16].
2.3. A left-censored dynamic linear model
In our model, the (unobserved) viral concentration signal $ X_t $ is modelled as a first-order random walk (RW1) process
where $ \epsilon_t \sim N(0, 1) $ is an independent and identically distributed normal random variable for $ t = 2, ..., n $. The measured concentrations $ Y^*_t $ are modelled by adding independent measurement noise to $ X_t $:
where the independent measurement error $ \epsilon'_{t} \sim t_{\nu} $ has a Student t-distribution with $ \nu $ degrees of freedom. The actually observed, censored data, are modelled by truncating $ Y^*_t $ at the known censoring threshold $ \ell_t $:
As samples are taken only four times a week, the vector of measurements $ {\bf{Y}} $ contains data observed at a subset $ \mathcal{T} $ of all $ n $ available time points. We infer the viral concentration $ X_t $ from $ {\bf{Y}} $ by Bayesian inference [17], i.e. by calculating the posterior distributions of the latent state $ X_1, \dots, X_n $ and hyperparameters $ \sigma $, $ \tau $ and $ \nu $, conditional on $ {\bf{Y}} $. The posterior distribution is given by:
The first line on the right hand side of Eq 2.7 is determined by the distribution of independent measurement errors, and left-censoring, of the $ Y_t $. The second term is determined by the RW1 time series model for the $ X_t $. The distributions $ p(\tau) $, $ p(\sigma) $, and $ p(\nu) $ in the last line are prior hyperparameter distributions: we specify uninformative uniform prior distributions for $ \tau > 0 $ and $ \sigma > 0 $, and a left-truncated Normal prior for $ \nu $, with prior expectation 3, prior variance 1, and truncated at 2. The parameters of the truncated Normal prior for $ \nu $ were selected by simulation and based on subjective judgements about the likely magnitude of measurement errors. The (multiplicative) proportionality constant in Eq 2.7 is inferred by Markov-Chain Monte-Carlo (MCMC) using the Stan software [14].
The hyperparameters $ \tau $ and $ \nu $ of the measurement process can be interpreted as measurement error variance (larger $ \tau $'s correspond to noisier measurements), and the tendency to generate outliers (smaller $ \nu $'s generate greater deviations from measured viral concentrations). Posterior estimates of these parameters are thus interesting for diagnostic purposes, e.g. to identify anomalous sites.
2.4. Dynamic linear model (DLM)
The DLM was implemented in the open-source programming language Stan [14], which provides efficient sampling of probabilistic models via MCMC and other inference algorithms. Code specifying the model is provided in Supplementary Information Figure S17. MCMC convergence statistics for the fit examples shown in Figure 1 can also be found in the SI (Figure S4–S8, Tables S1–S4).
2.5. 10-fold cross-validation
10-fold cross-validation was performed on the data across the 286 sites that had at least 30 samples. For each iteration:
● Raw SARS-CoV-2 N1 gc/L (with no normalisation for flow) was used with a $ \log_{10} $ transformation.
● Data were randomly split (90%/10%) into training and test sets
● Pre-existing missing values (days when samples were expected but were not collected) were included in the training set but not in the test set.
● The censoring threshold was set to a single value $ \log_{10} {\rm{(133.0 \ gc/L)}} $ for simplicity. In reality the limit of quantification will vary across samples.
● Fit DLM, KS and MA models to training data (details below)
● generate estimates (MA) and posterior samples (DLM, KS) of $ Y_{\text{pred, }}{t} $ at times $ t $ that were left out during training
$ Y_{\text{test, }}{t} $ (the left-out observation data) are then compared to $ Y_{\text{pred, }}{t} $ inferred with the three methods, via mean squared error (MSE) and interval coverage.
$ Y_{\text{pred, }}{t} $ for the DLM were generated by using Stan to sample from the joint posterior distribution of $ X_1, \dots, X_n $ and hyperparameters $ \sigma $, $ \nu $, $ \tau $, inferred from the training data. We then inferred posterior predictive samples $ Y_{\text{pred, }}{t} $ at times $ t $ left out during training by adding t-distributed measurement errors to posterior samples of $ X_t $, and applying censoring if the sampled observation was below the censoring threshold.
$ Y_{\text{pred, }}{t} $ for KS was generated by taking the fitted parameters $ \tau $, $ p_{\text{outlier}} $, $ \mu_{\text{X-test}} $, $ \sigma_{\text{X-test}} $. To get a posterior distribution on $ X_t $ 4000 samples were generated from a normal distribution with
To get $ O_t $ 4000 samples were generated from a binomial with
To get uncensored observations $ Y_{\text{outliers}} $ 4000 samples were taken from a uniform distribution
To get $ Y_{\text{pred, }}{t} $ $ X_t $ is passed to a Normal distribution with scale $ \tau $
Simulating outliers in $ Y_{\text{pred}} $ was done by
Finally, $ Y_{\text{pred, }}{t} $ is censored at some limit $ l $
2.6. Forward prediction
We then further validated model performance by testing how well it is able to predict 10 samples (2.5 weeks) ahead. We refit the model on all samples for all sites minus the final 10 samples, and then predict the left-out samples.
2.7. Exploratory analysis of $ \nu $ and $ \tau $ Outputs from the DLM
Analyses were performed using ${\textsf{R}}$ statistical software (Version 4.1.3) to establish whether the DLM is more likely to observe data variability - characterised by $ \nu $ and $ \tau $ outputs - at sites that show greater concentrations of SARS-CoV-2 RNA in wastewater. For this purpose we regressed the median gene copies per litre (gc/L) obtained over all time ($ \log_{10} $ transformed) against mean $ \nu $ or $ \tau $, controlling for the standard deviation of $ \nu $ and $ \tau $, respectively. We obtained the residuals from these linear regression models to identify sites where the $ \nu $ and $ \tau $ outputs from the model vary in excess of what is accountable to median gc/L and the posterior standard deviation of $ \nu $ or $ \tau $. Residuals for both models were then mapped to the Lower Layer Super Output Area (LSOA) for a given site using the Simple Features (sf) package in ${\textsf{R}}$ [18] (Figure 4).
3.
Results and discussion
3.1. Model fit: Simulated data
To test the output of the DLM, we first simulated data by generating a random walk with variance parameter $ \sigma^2 $ to model the underlying state, which was then sampled with measurement error parameters $ \tau $ and $ \nu $. Exact values are provided in the code. Any values below a predefined limit $ l $ are set to the value of $ l $. These synthetic data were then fit with the DLM. Supplementary Figure S1 shows that the underlying state $ X $-true is tracked rather well by the $ X $-smoothed estimate and lies within the inferred credible intervals, demonstrating that the model can reliably recover the underlying state from noisy observations in a synthetic dataset.
3.2. Model fit: Real data
We fit the DLM to data from 286 sewage treatment works across England, restricted to sites with greater than 30 samples present. Each site is sampled four times a week. Figure 1 shows a range of fitted sites selected, based upon their estimated parameters $ \tau $ and $ \nu $, to illustrate model behaviour at the extreme ends of the spectrum, i.e. low $ \nu $ and low $ \tau $ (Figures 1b and 1d) or high $ \nu $ and high $ \tau $ (Figures 1a and 1c). Sites with high parameter values typically show low levels of SARS-CoV-2 (N1 gene gc/L, the target used to approximate viral concentration in the sample) recovery and more frequent censoring. More censoring leads to more estimation uncertainty (wider credible intervals) as less information is available to constrain viral concentration estimates. Conversely, sites with low parameter values generally correspond to high levels of SARS-CoV-2 recovery and less censoring, therefore providing more information and tighter credible intervals. Supplementary Figure S2 shows a strong positive correlation between $ \nu $ and $ \tau $, and Figures 4b and 4b show a strong negative correlation between $ \nu $ and $ \tau $ in relation to the site's median viral RNA concentration $ (\log_{10} {\rm{(N1 \ gc/L)}}) $, respectively. We note that our model seems to produce realistic estimates of viral concentration during long periods of censoring, and on days where observations are missing entirely.
3.3. Model validation
Model performance was assessed by comparing the MSE produced by the DLM, KS and a seven-day centered MA over 10 folds of cross-validation (see Methods). The MA represents a simple way to remove noise from data, and is used here as a benchmark for comparison. All three models generated comparable MSE per site (Figure 2). However, the DLM and KS can estimate viral concentrations below the censoring threshold and, therefore, provide additional information on value for applications, such as case prevalence estimation (see Applications section). In addition both the DLM and KS provide useful parameters for quantifying uncertainty and outliers within the data (DLM: $ \sigma $/$ \tau $/$ \nu $, KS: $ \sigma $/$ \tau $/$ p_{\text{outlier}} $). This is particularly useful to identify sites generating unexpected data. So, while an MA scores equally well in terms of the MSE, the smoothing methods still confer additional advantages. A boxplot of the pairwise MSE differences, shown in Supplementary Figure S6, shows that the differences are not consistently better or worse for the DLM when compared to the KS or MA models.
As MSE assesses the accuracy of a single point estimate of the predictive distribution, it cannot inform on the reliability of the whole model distribution. In Figure 2, the coverage frequency of prediction intervals was used to characterise the reliability of the predictive distribution. Coverage frequency assesses how well the fitted model represents the variability of the data by analysing to what extent the observations could pass as a random sample from the predictive distribution. If observations and samples from the predictive distribution are statistically indistinguishable, we should expect a 90$ % $ chance that the observation is included in a 90% prediction interval derived from the predictive distribution. See Methods for information on calculating coverage. Figure 2b shows mean coverage frequencies across all sites. For nominal interval widths between 0.8 and 0.95, the KS coverage frequencies lie above the dashed line indicating that the model intervals are slightly wider than the true interval and are thus slightly under-confident. For the DLM, the coverage is too wide below nominal values of 90% and appears more reliable between 0.90 and 0.95 than KS. Both models appear over-confident at nominal values above 95%. For additional information on the distributions coverage values see Figures S4 and S5.
Cross-validation was also performed for forward prediction by removing the last 10 samples and predicting them with either the DLM or KS. Figures S9 and S10 show that both the DLM and KS perform equally well at forecasting up to 10 days of samples.
The DLM performed equally as well as the KS in cross-validation, but with greater parsimony: we removed the Bernoulli outlier functionality, and autoregressive and offset parameters ($ \eta $ and $ \delta $), to specify a simpler model. By providing full Bayesian posterior information, the DLM offers more information on the distributions of all the parameters in the model, thereby facilitating greater quantification of model uncertainty. Furthermore, the Stan framework offers flexibility for modification of the underlying state model (e.g. AR(2) random walk) or the addition of autoregressive parameters, if desired. The MCMC inference algorithm provided in Stan also allowed the model to be estimated more than 10x faster than the Kalman Smoother: the mean runtime of the DLM for each fold in 10-fold cross-validation of 10 sites was 14 seconds compared to 155 seconds with the KS, although with known parameters the prediction speed by the KS is much improved. Results of the test are provided in Supplementary Table S5, however in both cases the run times are small enough that we believe the difference is of little practical significance. The speed difference that is of more practical relevance (although difficult to quantify) is that our model was written in a general purpose modelling framework and so is easier to maintain, modify and adapt than the handcrafted R code of the Kalman Smoother. On the other hand, only the Kalman Smoother is able to quantify the probability of a given sample being an outlier and, therefore, this model will be more desirable for specific use cases. The DLM can only inform on whether a given sample lies outside of a predefined interval of the estimated underlying state, as shown in Figure 1.
4.
Applications
4.1. Comparison of smoothed estimate with Coronavirus (COVID-19) infection survey
Work by multiple groups has shown that SARS-CoV-2 gc/l concentrations in wastewater measurements can track case prevalence ('positivity rate', the percentage of people who have tested positive for COVID-19 on a polymerase chain reaction (PCR) test at a point in time) [19,20,21]. In England, the latter has been measured by the Office for National Statistics' Coronavirus (COVID-19) Infection Survey (CIS), a randomised household survey that provides an estimate of disease prevalence at sub-regional, regional and national levels [22]. Therefore, smoothed estimates of $ \log_{10} {\rm{(N1\; gc/L)}} $ from a DLM or KS can be compared with flow-normalised raw estimates to establish which correlates more strongly with $ \log_{10} {\rm{(CIS\; prevalence}} %) $ over time. Figure 3 compares correlations of CIS with (i) flow-normalised $ \log_{10} {\rm{(N1 \ gc/L)}} $, (ii) flow-normalised $ \log_{10} {\rm{(N1 \ gc/L)}} $ with a 7-day centered MA, and (iii) flow normalised $ \log_{10} {\rm{(N1 \ gc/L)}} $ smoothed estimate of $ X $ for all nine English regions between 1st September 2020 and 1st March 2022. This time range includes a period in which wastewater RNA concentration rates decoupled from clinical measures of disease prevalence, of which the cause is unknown [23]. It is worth noting that this relationship is likely not deterministic, i.e. they are not equivalent and are subject to their own spatiotemporal variation and uncertainty that would manifest in significant changes in the ratio of the measures. Such observations have not been limited to England, and the cause is likely to have multiple factors, both epidemiological (i.e., changes in viral shedding distribution as circulating virus variants emerge and evolve) and metrological (e.g., degree of clinical testing coverage can be demographically biased; laboratory sensitivity can vary significantly with virus concentration method employed for wastewater analysis) [24,25,26].
Smoothed wastewater concentration rates using a DLM or KS correlate more strongly with CIS positivity rate than raw or averaged rates (Figure 3). The enhanced correlation performance of the DLM and KS is likely due to both models' ability to generate data from below the censoring limit. This assertion is supported by the comparison between the smoothed estimates improvement in correlations verses the averaged $ \log_{10} {\rm{(N1 \ gc/L)}} $, which according to the MSE cross-validation should perform equally well. The key difference being that the DLM and KS infer values below the censored limit, thus we attribute at least part of the increase in correlation to this aspect of the models. Figure S9 provides a time series comparison of $ \log_{10} {\rm{(CIS\; prevalence}} %) $, $ \log_{10} {\rm{(N1 \ gc/L)}} $, and smoothed estimates. Smoothed estimates show a specific advantage over raw $ \log_{10} {\rm{(N1 \ gc/L)}} $ during times of low case prevalence. Using a simple sensitivity analysis to exclude the period in which wastewater concentration rates diverged from case rates to train the models, we find the same results (Figures S10–S11). DLM-smoothed rates therefore better complement CIS data, providing additional useful information for public health decision makers.
4.2. Exploration of fitted $ \nu $ and $ \tau $ parameters
The DLM provides two useful parameters for a given site: the extent to which outliers are observed ($ \nu $), with smaller values indicating greater frequency and size of outlier values, and the amount of measurement noise at a given site ($ \tau $), with larger values indicating noisier measurements. Figure 4a shows the geographical distribution of $ \nu $ values for fitted sites mapped to each Lower Layer Super Output Area (LSOA) in England. There is some evidence of localised behaviour, with areas of large $ \nu $ in the North and East, and low values found in the West and London regions. However, interpretation of this map is challenging as $ \nu $ is strongly related to $ {\mathrm{median}}(\log_{10} {\rm{(N1 \ gc/L)}}) $ and the quality of fit, quantified here as the posterior standard deviation of $ \nu $ ($ SD_{\nu} $ (Figure 4b). To account for these relationships and draw more insight from the $ \nu $ parameter we performed a multivariable linear regression analysis where $ {\mathrm{median}}(\log_{10} {\rm{(N1 \ gc/L)}}) $ was regressed onto the mean of the posterior of $ \nu $, controlling for the standard deviation of the posterior of $ \nu $ (see Methods). Figure 4c plots the absolute values of the regression model residuals (see Figure S12 for distribution of residuals); sites with the highest absolute residuals (i.e., the most variance not explained by either $ {\mathrm{median}}(\log_{10} {\rm{(N1 \ gc/L)}}) $ or quality of fit) are clustered in the North West. We repeat this analysis for $ \tau $ in Figure S3; again sites with the largest absolute residuals are concentrated in the North West, with additional large residuals seen in the South and East of England. Observed non-linearity is potentially attributable to high levels of censorship at low levels of $ {\mathrm{median}}(\log_{10} {\rm{(N1 \ gc/L)}}) $. Future analyses should explore this suggestion, potentially with a censored regression model.
Further examples of sites with high and low parameter values are provided in Figures S13–S16.
5.
Conclusions
We show that use of a Bayesian Dynamic Linear Model is a viable method for smoothing left-censored wastewater SARS-CoV-2 measurement data. Handling outliers through a t-distribution, rather than through an independent Bernoulli distribution, as applied in a previously published Kalman Smoother [12], is likely to more directly relate to the underlying state to be recovered. While the DLM and KS perform equivalently with mean squared error under cross-validation, the proposed DLM is more parsimonious (fewer model parameters), has a faster computational time, and is implemented in a more flexible modelling framework, allowing for easier modifications. Additionally, the DLM produces two site-specific parameters, $ \nu $ and $ \tau $, which are able to highlight sites with variable performance. This can be useful when assessing sampling strategies applied at scale (e.g. national or regional surveillance). Sites identified as providing inconsistent, noisy, or low information data may be removed from multi-site monitoring campaigns, for example.
The smoothed data, using our method, more closely correlate with regional infection survey data (CIS) than untransformed raw measurements. Wastewater data, smoothed in this fashion, are therefore more robust, capable of better complementing traditional surveillance, and providing additional confidence and utility for public health decision making.
Nevertheless, our approach has some limitations. The limit of censorship was set to a single value during cross-validation $ \log_{10}({\rm{133.0 \ gc/L)}} $, for simplicity. In reality this limit can vary across samples. From September 2021 SARS-CoV-2 RNA measurements from English wastewater diverged from reported clinical data where it had been previously tracking it. The reason why has still not been established but is potentially attributable to differential shedding rates between variants. Our sensitivity analyses reported in the Supplementary material found this does not impact the performance of our model.
Acknowledgments
The United Kingdom Government (Department of Health and Social Care) funded the sampling, testing, and data analysis of wastewater in England. Obépine funded the work of Marie Courbariaux and provided the R code of the modified Kalman Smoother and support to use it.
Conflict of Interest
The authors declare no conflicting interests in this paper.
Supplementary Information
DownLoad:
