Emulator-based Bayesian calibration of the CISNET colorectal cancer models

Purpose

To calibrate Cancer Intervention and Surveillance Modeling Network (CISNET) 's SimCRC, MISCAN-Colon, and CRC-SPIN simulation models of the natural history colorectal cancer (CRC) with an emulator-based Bayesian algorithm and internally validate the model-predicted outcomes to calibration targets.

Methods

We used Latin hypercube sampling to sample up to 50,000 parameter sets for each CISNET-CRC model and generated the corresponding outputs. We trained multilayer perceptron artificial neural networks (ANN) as emulators using the input and output samples for each CISNET-CRC model. We selected ANN structures with corresponding hyperparameters (i.e., number of hidden layers, nodes, activation functions, epochs, and optimizer) that minimize the predicted mean square error on the validation sample. We implemented the ANN emulators in a probabilistic programming language and calibrated the input parameters with Hamiltonian Monte Carlo-based algorithms to obtain the joint posterior distributions of the CISNET-CRC models' parameters. We internally validated each calibrated emulator by comparing the model-predicted posterior outputs against the calibration targets.

Results

The optimal ANN for SimCRC had four hidden layers and 360 hidden nodes, MISCAN-Colon had 4 hidden layers and 114 hidden nodes, and CRC-SPIN had one hidden layer and 140 hidden nodes. The total time for training and calibrating the emulators was 7.3, 4.0, and 0.66 hours for SimCRC, MISCAN-Colon, and CRC-SPIN, respectively. The mean of the model-predicted outputs fell within the 95% confidence intervals of the calibration targets in 98 of 110 for SimCRC, 65 of 93 for MISCAN, and 31 of 41 targets for CRC-SPIN.

Conclusions

Using ANN emulators is a practical solution to reduce the computational burden and complexity for Bayesian calibration of individual-level simulation models used for policy analysis, like the CISNET CRC models. In this work, we present a step-by-step guide to constructing emulators for calibrating three realistic CRC individual-level models using a Bayesian approach.

Efficient Bayesian automatic calibration of a functional-structural wheat model using an adaptive design and a metamodelling approach

Functional-structural plant models are increasingly being used by plant scientists to address a wide variety of questions. However, the calibration of these complex models is often challenging, mainly because of their high computational cost, and, as a result, error propagation is usually ignored. Here we applied an automatic method to the calibration of WALTer: a functional-structural wheat model that simulates the plasticity of tillering in response to competition for light. We used a Bayesian calibration method to jointly estimate the values of five parameters and quantify their uncertainty by fitting the model outputs to tillering dynamics data. We made recourse to Gaussian process metamodels in order to alleviate the computational cost of WALTer. These metamodels are built from an adaptive design that consists of successive runs of WALTer chosen by an efficient global optimization algorithm specifically adapted to this particular calibration task. The method presented here performed well on both synthetic and experimental data. It is an efficient approach for the calibration of WALTer and should be of interest for the calibration of other functional-structural plant models.

A Bayesian method with nonlinear noise model to calibrate constitutive parameters of soft tissue

The measured mechanical responses of soft tissue exhibit large variability and errors, especially for the softest brain tissue, while calibrating its constitutive parameters in a deterministic way remains a common practice. Here we implement a Bayesian method considering the nonlinear noise model to calibrate constitutive parameters of brain tissue. A probability model is first developed based on the measured experimental data, likelihood function, and prior function, from which the posterior distributions of model parameters are formulated. The likelihood function considers the nonlinear behaviors of the constitutive response and noise distribution of the experimentally measured data. Meanwhile, the prior predictive distribution is computed to check the probability model preliminarily. Secondly, the Markov Chain Monte Carlo (MCMC) method is used to compute the posterior distributions of model parameters, enabling assessment of parameter uncertainty, correlation, and model calibration error. Finally, the posterior predictive distributions of the overall response, constitutive response, and noise response are computed to validate the probabilistic model, all of which are consistent with the corresponding data. Furthermore, the effect of the prior distribution, experimental data, and noise model on posterior distribution is studied. Our study provides a general approach to calibrating constitutive parameters of soft tissue despite errors and large variability in experimental data.

Modified Taylor Impact Tests with Profiled Copper Cylinders: Experiment and Optimization of Dislocation Plasticity Model

Current progress in numerical simulations and machine learning allows one to apply complex loading conditions for the identification of parameters in plasticity models. This possibility expands the spectrum of examined deformed states and makes the identified model more consistent with engineering practice. A combined experimental-numerical approach to identify the model parameters and study the dynamic plasticity of metals is developed and applied to the case of cold-rolled OFHC copper. In the experimental part, profiled projectiles (reduced cylinders or cones in the head part) are proposed for the Taylor impact problem for the first time for material characterization. These projectiles allow us to reach large plastic deformations with true strains up to 1.3 at strain rates up to 105 s-1 at impact velocities below 130 m/s. The experimental results are used for the optimization of parameters of the dislocation plasticity model implemented in 3D with the numerical scheme of smoothed particle hydrodynamics (SPH). A Bayesian statistical method in combination with a trained artificial neural network as an SPH emulator is applied to optimize the parameters of the dislocation plasticity model. It is shown that classical Taylor cylinders are not enough for a univocal selection of the model parameters, while the profiled cylinders provide better optimization even if used separately. The combination of different shapes and an increase in the number of experiments increase the quality of optimization. The optimized numerical model is successfully validated by the experimental data about the shock wave profiles in flyer plate experiments from the literature. In total, a cheap, simple, but efficient route for optimizing a dynamic plasticity model is proposed. The dislocation plasticity model is extended to estimate grain refinement and volume fractions of weakened areas in comparison with experimental observations.

Demography-environment relationships improve mechanistic understanding of range dynamics under climate change

Species respond to climate change with range and abundance dynamics. To better explain and predict them, we need a mechanistic understanding of how the underlying demographic processes are shaped by climatic conditions. Here, we aim to infer demography-climate relationships from distribution and abundance data. For this, we developed spatially explicit, process-based models for eight Swiss breeding bird populations. These jointly consider dispersal, population dynamics and the climate-dependence of three demographic processes-juvenile survival, adult survival and fecundity. The models were calibrated to 267 nationwide abundance time series in a Bayesian framework. The fitted models showed moderate to excellent goodness-of-fit and discriminatory power. The most influential climatic predictors for population performance were the mean breeding-season temperature and the total winter precipitation. Contemporary climate change benefitted the population trends of typical mountain birds leading to lower population losses or even slight increases, whereas lowland birds were adversely affected. Our results emphasize that generic process-based models embedded in a robust statistical framework can improve our predictions of range dynamics and may allow disentangling of the underlying processes. For future research, we advocate a stronger integration of experimental and empirical studies in order to gain more precise insights into the mechanisms by which climate affects populations. This article is part of the theme issue 'Detecting and attributing the causes of biodiversity change: needs, gaps and solutions'.

Tumour growth: An approach to calibrate parameters of a multiphase porous media model based on in vitro observations of Neuroblastoma spheroid growth in a hydrogel microenvironment

To unravel processes that lead to the growth of solid tumours, it is necessary to link knowledge of cancer biology with the physical properties of the tumour and its interaction with the surrounding microenvironment. Our understanding of the underlying mechanisms is however still imprecise. We therefore developed computational physics-based models, which incorporate the interaction of the tumour with its surroundings based on the theory of porous media. However, the experimental validation of such models represents a challenge to its clinical use as a prognostic tool. This study combines a physics-based model with in vitro experiments based on microfluidic devices used to mimic a three-dimensional tumour microenvironment. By conducting a global sensitivity analysis, we identify the most influential input parameters and infer their posterior distribution based on Bayesian calibration. The resulting probability density is in agreement with the scattering of the experimental data and thus validates the proposed workflow. This study demonstrates the huge challenges associated with determining precise parameters with usually only limited data for such complex processes and models, but also demonstrates in general how to indirectly characterise the mechanical properties of neuroblastoma spheroids that cannot feasibly be measured experimentally.

PBTK modeled perfluoroalkyl acid kinetics in zebrafish eleutheroembryos suggests impacts on bioconcentrations by chorion porosity dynamics

The zebrafish eleutheroembryo (zfe) is widely used as a model to characterize the toxicity of chemicals. However, analytical methods are still missing to measure organ concentrations. Therefore, physiologically-based toxicokinetic (PBTK) modeling may overcome current limitations to help understand the relationship between toxic effects and internal exposure in various organs. A previous PBTK model has been updated to include the chorionic transport barrier and its permeabilization, hatching dynamics within a zfe population over development, and active mediated transport mechanisms. The zfe PBTK model has been calibrated using measured time-dependent internal concentrations of PFBA, PFHxS, PFOA, and PFOS in a zfe population and evaluated using external datasets from the literature. Calibration was successful with 96% of the predictions falling within a 2-fold range of the observed concentrations. The external dataset was correctly estimated with about 50% of the predictions falling within a factor of 3 of the observed data and 10% of the predictions are out of the 10-fold error. The calibrated model suggested that active mediated transport differs between PFAS with a sulfonic and carboxylic acid functional end groups. This PBTK model predicts well the fate of PFAS with various physicochemical properties in zfe. Therefore, this model may improve the use of zfe as an alternative model in toxicokinetic-toxicodynamic studies and help to refine and reduce zfe-based experiments, while giving insights into the internal kinetics of chemicals.

Usefulness of cultivar-level calibration of AquaCrop for vegetables depends on the crop and data availability

As a result of climate change, climatic extremes are expected to increase. For high-value crops like vegetables, irrigation is a potentially economically viable adaptation measure in western Europe. To optimally schedule irrigation, decision support systems based on crop models like AquaCrop are increasingly used by farmers. High value vegetable crops like cauliflower or spinach are grown in two distinct growth cycles per year and, additionally, have a high turnover rate of new varieties. To successfully deploy the AquaCrop model in a decision support system, it requires a robust calibration. However, it is not known whether parameters can be conserved over both growth periods, nor whether a cultivar dependent model calibration is always required. Furthermore, when data are collected from farmers' fields, there are constraints in data availability and uncertainty. We collected data from commercial cauliflower and spinach fields in Belgium in 2019, 2020 and 2021 during different growing periods and of different cultivars. With the use of a Bayesian calibration, we confirmed the need for a condition or cultivar specific calibration for cauliflower, while for spinach, splitting the data per cultivar or pooling the data together did not improve uncertainty on the model simulations. However, due to uncertainties arising from field specific soil and weather conditions, or measurement errors from calibration data, real time field specific adjustments are advised to simulations when using AquaCrop as decision support tool. Remotely sensed or in situ ground data may be invaluable information to reduce uncertainty on model simulations.

Application of Parameter Optimization Methods Based on Kalman Formula to the Soil-Crop System Model

Soil-crop system models are effective tools for optimizing water and nitrogen application schemes, saving resources and protecting the environment. To guarantee model prediction accuracy, we must apply parameter optimization methods for model calibration. The performance of two different parameter optimization methods based on the Kalman formula are evaluated for a parameter identification of the soil Water Heat Carbon Nitrogen Simulator (WHCNS) model using mean bias error (ME), root-mean-square error (RMSE) and an index of agreement (IA). One is the iterative local updating ensemble smoother (ILUES), and the other is the DiffeRential Evolution Adaptive Metropolis with Kalman-inspired proposal distribution (DREAMkzs). Our main results are as follows: (1) Both ILUES and DREAMkzs algorithms performed well in model parameter calibration with the RMSE_Maximum a posteriori (RMSE_MAP) values were 0.0255 and 0.0253, respectively; (2) ILUES significantly accelerated the process to the reference values in the artificial case, while outperforming in the calibration of multimodal parameter distribution in the practical case; and (3) the DREAMkzs algorithm considerably accelerated the burn-in process compared with the original algorithm without Kalman-formula-based sampling for parameter optimization of the WHCNS model. In conclusion, ILUES and DREAMkzs can be applied to a parameter identification of the WHCNS model for more accurate prediction results and faster simulation efficiency, contributing to the popularization of the model.

Experimentally characterizing the spatially varying anisotropic mechanical property of cancellous bone via a Bayesian calibration method

Traditional experimental tests for characterizing bone's mechanical properties usually hypothesize a uniaxial stress condition without quantitatively evaluating the influence of spatially varying principal material orientations, which cannot accurately predict the mechanical properties distribution of bones in vivo environment. In this study, a Bayesian calibrating procedure was developed using quantified multiaxial stress to investigate cancellous bone's local anisotropic elastic performance around joints as the spatial variation of main bearing orientations. First, the bone cube specimens from the distal femur of sheep are prepared using traditional anatomical axes. The multiaxial stress state of each bone specimen is calibrated using the actual principal material orientations derived from fabric tensor at different anatomical locations. Based on the calibrated multiaxial stress state, the process of identifying mechanical properties is described as an inverse problem. Then, a Bayesian calibration procedure based on a surrogate constitutive model was developed via multiaxial stress correction to identify the anisotropic material parameters. Finally, a comparison between the experiment and simulation results is discussed by applying the optimal model parameters obtained from the Bayesian probability distribution. Compared to traditional uniaxial methods, our results prove that the calibration based on the spatial variation of the main bearing orientations can significantly improve the accuracy of characterizing regional anisotropic mechanical responses. Moreover, we determine that the actual mechanical property distribution is influenced by complicated mechanical stimulation. This study provides a novel method to evaluate the spatially varying mechanical properties of bone tissues enduring complex mechanical loading accurately and effectively. It is expected to provide more realistic mechanical design targets in vivo for a personalized artificial bone prosthesis in clinical treatment.

Combining Bayesian Calibration and Copula Models for Age Estimation

Accurately estimating and predicting chronological age from some anthropometric characteristics of an individual without an identity document can be crucial in the context of a growing number of forced migrants. In the related literature, the prediction of chronological age mostly relies upon the use of a single predictor, which is usually represented by a dental/skeletal maturity index, or multiple independent ordinal predictor (stage of maturation). This paper is the first attempt to combine a robust method to predict chronological age, such as Bayesian calibration, and the use of multiple continuous indices as predictors. The combination of these two aspects becomes possible due to the implementation of a complex statistical tool as the copula. Comparing the forecasts from our copula-based method with predictions from an independent model and two single predictor models, we showed that the accuracy increased.

A PBPK model to evaluate zebrafish eleutheroembryos' actual exposure: bisphenol A and analogs' (AF, F, and S) case studies

The zebrafish eleutheroembryo model is increasingly used to assess the toxicity and developmental adverse effects of xenobiotics. However, the actual exposure is seldom measured (poorly accessible), while a predictive model could estimate these concentrations. The predictions with a new eleutheroembryo physiologically based pharmacokinetic (PBPK) model have been evaluated using datasets obtained from literature data for several bisphenols. The model simulated the toxicokinetics of bisphenols A (BPA), AF, F, and S through the eleutheroembryo tissues while considering the body and organ growth. We further improved the predictions by adding dynamic flows through the embryo and/or its chorion, impact of experimental temperature, metabolic clearance, and saturation of the absorption by Bayesian calibration. The model structure was determined using the BPA dataset and generalized to the other bisphenols. This model revealed the central role of the chorion in the compound uptake in the first 48 h post-fertilization. The predictions for the BPA substitutes estimated by our PBPK model were compared to available toxicokinetics data for zebrafish embryos, and 63% and 88% of them were within a twofold and fivefold error intervals of the corresponding experimental values, respectively. This model provides a tool to design new eleutheroembryo assays and evaluate the actual exposure.

Model selection for assessing the effects of doxorubicin on triple-negative breast cancer cell lines

Doxorubicin is a chemotherapy widely used to treat several types of cancer, including triple-negative breast cancer. In this work, we use a Bayesian framework to rigorously assess the ability of ten different mathematical models to describe the dynamics of four TNBC cell lines (SUM-149PT, MDA-MB-231, MDA-MB-453, and MDA-MB-468) in response to treatment with doxorubicin at concentrations ranging from 10 to 2500 nM. Each cell line was plated and serially imaged via fluorescence microscopy for 30 days following 6, 12, or 24 h of in vitro drug exposure. We use the resulting data sets to estimate the parameters of the ten pharmacodynamic models using a Bayesian approach, which accounts for uncertainties in the models, parameters, and observational data. The ten candidate models describe the growth patterns and degree of response to doxorubicin for each cell line by incorporating exponential or logistic tumor growth, and distinct forms of cell death. Cell line and treatment specific model parameters are then estimated from the experimental data for each model. We analyze all competing models using the Bayesian Information Criterion (BIC), and the selection of the best model is made according to the model probabilities (BIC weights). We show that the best model among the candidate set of models depends on the TNBC cell line and the treatment scenario, though, in most cases, there is great uncertainty in choosing the best model. However, we show that the probability of being the best model can be increased by combining treatment data with the same total drug exposure. Our analysis points to the importance of considering multiple models, built on different biological assumptions, to capture the observed variations in tumor growth and treatment response.

Flocculation kinetics and mechanisms of microalgae- and clay-containing suspensions in different microalgal growth phases

Interplays between microalgae and clay minerals enhance biologically mediated flocculation, thereby affecting the sedimentation and transportation of suspended particulate matter (SPM) in water and benthic environments. This interaction forms larger flocs with a higher settling velocity and enhances SPM sinking. The aim of this study was to investigate the flocculation kinetics of microalgae and clay in suspension and to elucidate the mechanisms associated with such interactions. Standard jar test experiments were conducted using various mixtures of kaolinite and microalgal samples from batch cultures (Chlorella vulgaris) to estimate biologically mediated flocculation kinetics. The organic matter (OM) composition secreted by the microalgae was characterized using a liquid chromatography - organic carbon detection system, and quantitative analysis of transparent exopolymer particles was conducted separately. A two-class flocculation kinetic model, based on the interaction between flocculi and flocs, was also adopted to quantitatively analyze the experimental data from flocculation. Results from the flocculation kinetic tests and OM analyses, in association with other data analyses (i.e., floc size distribution and flocculation kinetic model), showed that flocculation increased with OM concentration during the growth phase (10-20 d). However, on day 23 during the early stationary phase, flocculation kinetics started decreasing and substantially declined on day 30, even though the amount of OM (mainly biopolymers) continued to increase. Our results indicate that an adequate quantity of biopolymers produced by the microalgal cells in the growth phase enhanced floc-to-floc attachment and hence flocculation kinetics. In contrast, an excessive quantity of biopolymers and humic substances in the stationary phase enhanced the formation of polymeric backbone structures and flocculation via scavenging particles but simultaneously increased steric stabilization with the production of a large number of fragmented particles.

Development of a Bayesian inference model for assessing ventilation condition based on CO2 meters in primary schools

Outdoor fresh air ventilation plays a significant role in reducing airborne transmission of diseases in indoor spaces. School classrooms are considerably challenged during the COVID-19 pandemic because of the increasing need for in-person education, untimely and incompleted vaccinations, high occupancy density, and uncertain ventilation conditions. Many schools started to use CO₂ meters to indicate air quality, but how to interpret the data remains unclear. Many uncertainties are also involved, including manual readings, student numbers and schedules, uncertain CO₂ generation rates, and variable indoor and ambient conditions. This study proposed a Bayesian inference approach with sensitivity analysis to understand CO₂ readings in four primary schools by identifying uncertainties and calibrating key parameters. The outdoor ventilation rate, CO₂ generation rate, and occupancy level were identified as the top sensitive parameters for indoor CO₂ levels. The occupancy schedule becomes critical when the CO₂ data are limited, whereas a 15-min measurement interval could capture dynamic CO₂ profiles well even without the occupancy information. Hourly CO₂ recording should be avoided because it failed to capture peak values and overestimated the ventilation rates. For the four primary school rooms, the calibrated ventilation rate with a 95% confidence level for fall condition is 1.96±0.31 ACH for Room #1 (165 m³ and 20 occupancies) with mechanical ventilation, and for the rest of the naturally ventilated rooms, it is 0.40±0.08 ACH for Room #2 (236 m³ and 21 occupancies), 0.30±0.04 or 0.79±0.06 ACH depending on occupancy schedules for Room #3 (236 m³ and 19 occupancies), 0.40±0.32,0.48±0.37,0.72±0.39 ACH for Room #4 (231 m³ and 8-9 occupancies) for three consecutive days.

Experiments and hyperelastic modeling of porcine meniscus show heterogeneity at high strains

Constitutive modeling of the meniscus is critical in areas like knee surgery and tissue engineering. At low strain rates, the meniscus can be described using a hyperelastic model. Calibration of hyperelastic material models of the meniscus is challenging on many fronts due to material variability and friction. In this study, we present a framework to determine the hyperelastic material parameters of porcine meniscus (and similar soft tissues) using no-slip uniaxial compression experiments. Because of the nonhomogeneous deformation in the specimens, a finite element solution is required at each step of the iterative calibration process. We employ a Bayesian calibration approach to account for the inherent material variability and a Bayesian optimization approach to minimize the resulting cost function in the material parameter space. Cylindrical specimens of porcine meniscus from the anterior, middle and posterior regions are tested up to 30% compressive strain and the Yeoh form of hyperelastic strain energy density function is used to describe the material response. The results show that the Yeoh form is able to accurately describe the compressive response of porcine meniscus and that the Bayesian calibration and optimization approaches are able to calibrate the model in a computationally efficient manner while taking into account the inherent material variability. The results also show that the shear modulus or the initial stiffness is roughly uniform across the different areas of the meniscus, but there is significant spatial heterogeneity in the response at high strains. In particular, the middle region is considerably stiffer at high strains. This heterogeneity is important to consider in modeling the response of the meniscus for clinical applications.

Surrogate-based Bayesian calibration of biomechanical models with isotropic material behavior

This work introduces a computational methodology to calibrate material models in biomechanical applications under uncertainty. We adopt a Bayesian approach, which estimates the probability distributions of hyperelastic material parameters, based on force-strain measurements. We approximate the parametric biomechanical model by combining a reduced order representation of the force response with a Polynomial Chaos expansion. The surrogate model allows to employ sampling-intensive Markov chain Monte Carlo methods and provides an efficient way to estimate (generalized) Sobol coefficients. We use a Sobol sensitivity analysis to assess the influence of material parameters and present an iterative procedure to quantify the accuracy of the surrogate model as additional uncertainty during Bayesian updating. The methodology is illustrated with three cases, tensile experiments on heat-induced whey protein gel, indentation experiments for oocytes and a manufactured example. Real experimental data are used for the calibration.

Modeling nitrous oxide mitigation potential of enhanced efficiency nitrogen fertilizers from agricultural systems

Agriculture soils are responsible for a large proportion of global nitrous oxide (N₂O) emissions-a potent greenhouse gas and ozone depleting substance. Enhanced-efficiency nitrogen (N) fertilizers (EENFs) can reduce N₂O emission from N-fertilized soils, but their effect varies considerably due to a combination of factors, including climatic conditions, edaphic characteristics and management practices. In this study, we further developed the DayCent ecosystem model to simulate two EENFs: controlled-release N fertilizers (CRNFs) and nitrification inhibitors (NIs) and evaluated their N₂O mitigation potentials. We implemented a Bayesian calibration method using the sampling importance resampling (SIR) algorithm to derive a joint posterior distribution of model parameters that was informed by N₂O flux measurements from corn production systems a network of experimental sites within the GRACEnet program. The joint posterior distribution can be applied to estimate predictions of N₂O reduction factors when EENFs are adopted in place of conventional urea-based N fertilizer. The resulting median reduction factors were - 11.9% for CRNFs (ranging from -51.7% and 0.58%) and - 26.7% for NIs (ranging from -61.8% to 3.1%), which is comparable to the measured reduction factors in the dataset. By incorporating EENFs, the DayCent ecosystem model is able to simulate a broader suite of options to identify best management practices for reducing N₂O emissions.

Formulation and probabilistic assessment of reversible biodegradation pathway of Diclofenac in groundwater

We present a conceptual and mathematical framework leading to the development of a biodegradation model capable to interpret the observed reversibility of the Pharmaceutical Sodium Diclofenac along its biological degradation pathway in groundwater. Diclofenac occurrence in water bodies poses major concerns due to its persistent (and bioactive) nature and its detection in surface waters and aquifer systems. Despite some evidences of its biodegradability at given reducing conditions, Diclofenac attenuation is often interpreted with models which are too streamlined, thus potentially hampering appropriate quantification of its fate. In this context, we propose a modeling framework based on the conceptualization of the molecular mechanisms of Diclofenac biodegradation which we then embed in a stochastic context, thus enabling one to quantify predictive uncertainty. We consider reference environmental conditions (biotic and denitrifying) associated with a set of batch experiments that evidence the occurrence of a reversible biotransformation pathway, a feature that is fully captured by our model. The latter is then calibrated in the context of a Bayesian modeling framework through an Acceptance-Rejection Sampling approach. By doing so, we quantify the uncertainty associated with model parameters and predicted Diclofenac concentrations. We discuss the probabilistic nature of uncertain model parameters and the challenges posed by their calibration with the available data. Our results are consistent with the recalcitrant behavior exhibited by Diclofenac in groundwater and documented through experimental data and support the observation that unbiased estimates of the hazard posed by Diclofenac to water resources should be assessed through a modeling strategy which fully embeds uncertainty quantification.

Application of Distributed Agent-based Modeling to Investigate Opioid Use Outcomes in Justice Involved Populations

Criminal justice involved (CJI) individuals with a history of opioid use disorder (OUD) are at high risk of overdose and death in the weeks following release from jail. We developed the Justice-Community Circulation Model (JCCM) to investigate OUD/CJI dynamics post-release and the effects of interventions on overdose deaths. The JCCM uses a synthetic agent-based model population of approximately 150,000 unique individuals that is generated using demographic information collected from multiple Chicago-area studies and data sets. We use a high-performance computing (HPC) workflow to implement a sequential approximate Bayesian computation algorithm for calibrating the JCCM. The calibration results in the simulated joint posterior distribution of the JCCM input parameters. The calibrated model is used to investigate the effects of a naloxone intervention for a mass jail release. The simulation results show the degree to which a targeted intervention focusing on recently released jail inmates can help reduce the risk of death from opioid overdose.

BayCANN: Streamlining Bayesian Calibration With Artificial Neural Network Metamodeling

Purpose: Bayesian calibration is generally superior to standard direct-search algorithms in that it estimates the full joint posterior distribution of the calibrated parameters. However, there are many barriers to using Bayesian calibration in health decision sciences stemming from the need to program complex models in probabilistic programming languages and the associated computational burden of applying Bayesian calibration. In this paper, we propose to use artificial neural networks (ANN) as one practical solution to these challenges. Methods: Bayesian Calibration using Artificial Neural Networks (BayCANN) involves (1) training an ANN metamodel on a sample of model inputs and outputs, and (2) then calibrating the trained ANN metamodel instead of the full model in a probabilistic programming language to obtain the posterior joint distribution of the calibrated parameters. We illustrate BayCANN using a colorectal cancer natural history model. We conduct a confirmatory simulation analysis by first obtaining parameter estimates from the literature and then using them to generate adenoma prevalence and cancer incidence targets. We compare the performance of BayCANN in recovering these "true" parameter values against performing a Bayesian calibration directly on the simulation model using an incremental mixture importance sampling (IMIS) algorithm. Results: We were able to apply BayCANN using only a dataset of the model inputs and outputs and minor modification of BayCANN's code. In this example, BayCANN was slightly more accurate in recovering the true posterior parameter estimates compared to IMIS. Obtaining the dataset of samples, and running BayCANN took 15 min compared to the IMIS which took 80 min. In applications involving computationally more expensive simulations (e.g., microsimulations), BayCANN may offer higher relative speed gains. Conclusions: BayCANN only uses a dataset of model inputs and outputs to obtain the calibrated joint parameter distributions. Thus, it can be adapted to models of various levels of complexity with minor or no change to its structure. In addition, BayCANN's efficiency can be especially useful in computationally expensive models. To facilitate BayCANN's wider adoption, we provide BayCANN's open-source implementation in R and Stan.

Bayesian versus Empirical Calibration of Microsimulation Models: A Comparative Analysis

Calibration of a microsimulation model (MSM) is a challenging but crucial step for the development of a valid model. Numerous calibration methods for MSMs have been suggested in the literature, most of which are usually adjusted to the specific needs of the model and based on subjective criteria for the selection of optimal parameter values. This article compares 2 general approaches for calibrating MSMs used in medical decision making, a Bayesian and an empirical approach. We use as a tool the MIcrosimulation Lung Cancer (MILC) model, a streamlined, continuous-time, dynamic MSM that describes the natural history of lung cancer and predicts individual trajectories accounting for age, sex, and smoking habits. We apply both methods to calibrate MILC to observed lung cancer incidence rates from the Surveillance, Epidemiology and End Results (SEER) database. We compare the results from the 2 methods in terms of the resulting parameter distributions, model predictions, and efficiency. Although the empirical method proves more practical, producing similar results with smaller computational effort, the Bayesian method resulted in a calibrated model that produced more accurate outputs for rare events and is based on a well-defined theoretical framework for the evaluation and interpretation of the calibration outcomes. A combination of the 2 approaches is an alternative worth considering for calibrating complex predictive models, such as microsimulation models.

Statistical finite elements for misspecified models

We present a statistical finite element method for nonlinear, time-dependent phenomena, illustrated in the context of nonlinear internal waves (solitons). We take a Bayesian approach and leverage the finite element method to cast the statistical problem as a nonlinear Gaussian state-space model, updating the solution, in receipt of data, in a filtering framework. The method is applicable to problems across science and engineering for which finite element methods are appropriate. The Korteweg-de Vries equation for solitons is presented because it reflects the necessary complexity while being suitably familiar and succinct for pedagogical purposes. We present two algorithms to implement this method, based on the extended and ensemble Kalman filters, and demonstrate effectiveness with a simulation study and a case study with experimental data. The generality of our approach is demonstrated in SI Appendix, where we present examples from additional nonlinear, time-dependent partial differential equations (Burgers equation, Kuramoto-Sivashinsky equation).

Development of a fate and transport model for biodegradation of PBDE congeners in sediments

Polybrominated diphenyl ethers (PBDEs) are a family where each congener possesses different physicochemical properties, persistence and/or toxicity. Biodegradation can selectively change the abundance of congeners. These warrant modeling of individual congeners by considering biodegradation pathways together with fate and transport (F&T) mechanisms. Accordingly, this study aims to develop a F&T model (Fate and Transport model for Hydrophobic Pollutants - FTHP) that integrates congener specific biodegradation of PBDEs in sediments. The model is tested using sediment data from a location representing the Lower South Bay of San Francisco. Results demonstrated settling, resuspension, and biodegradation as important mechanisms. FTHP is then used to predict congener concentrations in a period of 20 years for two cases (constant and time-dependent water column concentrations) and four alternative scenarios: no intervention (i.e., natural attenuation, also serves as the base case), no degradation, dredging and biostimulation. The greatest impact on the reduction of total PBDE concentrations was achieved by a reduction in water column concentrations, i.e. source control, and dredging. On the other hand, biostimulation coupled with source control was the most effective in reducing bioaccumulative PBDE congener concentrations and almost as effective as dredging for the rest of congeners. Proposed FTHP model can distinguish between congeners and help devise informed management plans which focus on decreasing risks associated with persistent and bioaccumulative compounds in contaminated sediments.

Methods for Model Calibration under High Uncertainty: Modeling Cholera in Bangladesh

Background. Published data on a disease do not always correspond directly to the parameters needed to simulate natural history. Several calibration methods have been applied to computer-based disease models to extract needed parameters that make a model's output consistent with available data. Objective. To assess 3 calibration methods and evaluate their performance in a real-world application. Methods. We calibrated a model of cholera natural history in Bangladesh, where a lack of active surveillance biases available data. We built a cohort state-transition cholera natural history model that includes case hospitalization to reflect the passive surveillance data-generating process. We applied 3 calibration techniques: incremental mixture importance sampling, sampling importance resampling, and random search with rejection sampling. We adapted these techniques to the context of wide prior uncertainty and many degrees of freedom. We evaluated the resulting posterior parameter distributions using a range of metrics and compared predicted cholera burden estimates. Results. All 3 calibration techniques produced posterior distributions with a higher likelihood and better fit to calibration targets as compared with prior distributions. Incremental mixture importance sampling resulted in the highest likelihood and largest number of unique parameter sets to better inform joint parameter uncertainty. Compared with naïve uncalibrated parameter sets, calibrated models of cholera in Bangladesh project substantially more cases, many of which are not detected by passive surveillance, and fewer deaths. Limitations. Calibration cannot completely overcome poor data quality, which can leave some parameters less well informed than others. Calibration techniques may perform differently under different circumstances. Conclusions. Incremental mixture importance sampling, when adapted to the context of high uncertainty, performs well. By accounting for biases in data, calibration can improve model projections of disease burden.

Bayesian calibration of a growth-dependent tree mortality model to simulate the dynamics of European temperate forests

Dynamic vegetation models (DVMs) are important tools to understand and predict the functioning and dynamics of terrestrial ecosystems under changing environmental conditions. In these models, uncertainty in the description of demographic processes, in particular tree mortality, is a persistent problem. Current mortality formulations lack realism and are insufficiently constrained by empirical evidence. It has been suggested that empirically estimated mortality submodels would enhance DVM performance, but due to the many processes and interactions within a DVM, the claim has rarely been tested. Here, we compare the performance of three alternative growth-dependent tree mortality submodels in the DVM ForClim: (1) a mortality function with theoretical foundation (ForClim v3.3); (2) a mortality function with parameters directly estimated based on forest inventory data; and (3) the same function, but with parameters estimated using an inverse approach through Bayesian calibration (BC). Time series of inventory data from 30 ecologically distinct Swiss natural forest reserves collected over 35+ yr, including the main tree species of Central Europe, were used for the calibration and subsequent validation of the mortality functions and the DVM. The recalibration resulted in mortality parameters that differed from the direct empirical estimates, particularly for the relationship between tree size and mortality. The calibrated parameters outperformed the direct estimates, and to a lesser extent the original mortality function, for predicting decadal-scale forest dynamics at both calibration and validation sites. The same pattern was observed regarding the plausibility of their long-term projections under contrasting environmental conditions. Our results demonstrate that inverse calibration may be useful even when direct empirical estimates of DVM parameters are available, as structural model deficiencies or data problems can result in discrepancies between direct and inverse estimates. Thus, we interpret the good performance of the inversely calibrated model for long-term projections (which were not a calibration target) as evidence that the calibration did not compensate for model errors. Rather, we surmise that the discrepancy was mainly caused by a lack of representativeness of the mortality data. Our results underline the potential for learning more about elusive processes, such as tree mortality or recruitment, through data integration in DVMs.

Accuracy and biases in predicting the chemical and physical traits of many types of cheeses using different visible and near-infrared spectroscopic techniques and spectrum intervals

Near-infrared spectroscopy (NIRS) has been widely used to determine various composition traits of many dairy products in the industry. In the last few years, near-infrared (NIR) instruments have become more and more accessible, and now, portable devices can be easily used in the field, allowing the direct measurement of important quality traits. However, the comparison of the predictive performances of different NIR instruments is not simple, and the literature is lacking. These instruments may use different wavelength intervals and calibration procedures, making it difficult to establish whether differences are due to the spectral interval, the chemometric approach, or the instrument's technology. Hence, the aims of this study were (1) to evaluate the prediction accuracy of chemical contents (5 traits), pH, texture (2 traits), and color (5 traits) of 37 categories of cheese; (2) to compare 3 instruments [2 benchtop, working in reflectance (R) and transmittance (T) mode (NIRS-R and NIRS-T, respectively) and 1 portable device (VisNIRS-R)], using their entire spectral ranges (1100-2498, 850-1048, and 350-1830 nm, respectively, for NIRS-R, NIRS-T and VisNIRS-R); (3) to examine different wavelength intervals of the spectrum within instrument, comparing also the common intervals among the 3 instruments; and (4) to determine the presence of bias in predicted traits for specific cheese categories. A Bayesian approach was used to develop 8 calibration models for each of 13 traits. This study confirmed that NIR spectroscopy can be used to predict the chemical composition of a large number of different cheeses, whereas pH and texture traits were poorly predicted. Color showed variable predictability, according to the trait considered, the instrument used, and, within instrument, according to the wavelength intervals. The predictive performance of the VisNIRS-R portable device was generally better than the 2 laboratory NIRS instruments, whether with the entire spectrum or selected intervals. The VisNIRS-R was found suitable for analyzing chemical composition in real time, without the need for sample uptake and processing. Our results also indicated that instrument technology is much more important than the NIR spectral range for accurate prediction equations, but the visible range is useful when predicting color traits, other than lightness. Specifically for certain categories (i.e., caprine, moldy, and fresh cheeses), dedicated calibrations seem to be needed to obtain unbiased and more accurate results.

Segmented Bayesian calibration approach for estimating age in forensic science

Forensic age estimation is receiving growing attention from researchers in the last few years. Accurate estimates of age are needed both for identifying real age in individuals without any identity document and assessing it for human remains. The methods applied in such context are mostly based on radiological analysis of some anatomical districts and entail the use of a regression model. However, estimating chronological age by regression models leads to overestimated ages in younger subjects and underestimated ages in older ones. We introduced a full Bayesian calibration method combined with a segmented function for age estimation that relied on a Normal distribution as a density model to mitigate this bias. In this way, we were also able to model the decreasing growth rate in juveniles. We compared our new Bayesian-segmented model with other existing approaches. The proposed method helped producing more robust and precise forecasts of age than compared models while exhibited comparable accuracy in terms of forecasting measures. Our method seemed to overcome the estimation bias also when applied to a real data set of South-African juvenile subjects.

Cost-utility analysis of a potential lung cancer screening program for a high-risk population in Germany: A modelling approach

BACKGROUND

Lung cancer is the leading cause of cancer death in Germany. Although several randomized trials in Europe have evaluated the effectiveness of lung cancer screening programs, evidence on the cost-effectiveness of lung cancer screening is scarce.

OBJECTIVE

To evaluate the cost-effectiveness of a population-based lung cancer screening program from the perspective of a German payer.

METHODS

We conducted a cost-effectiveness analysis from the public payer perspective for a high-risk population defined as heavy former and current smokers (≥20 cigarettes per day) between 55 and 75 years of age. The underlying model consisted of two Markov models. We differentiated between a population-based annual screening program and standard clinical care. Depending on stage at diagnosis, simulated patients were assigned to one of five treatment paths according to the German clinical guideline for the diagnosis and treatment of lung cancer. Costs, life years saved, and quality adjusted life years (QALYs) were used as outcomes. Values for input parameters were taken from the literature. The model was run for 60 cycles with a cycle length of three months. Deterministic and probabilistic sensitivity analyses were conducted.

RESULTS

In the base case, annual lung cancer screening led to an increase in incremental costs (€ 1,153 per person) compared to standard clinical care. However, the screening approach was associated with an incremental gain in life years (0.06 per person) and QALYs (0.04 per person). Thus, the incremental cost-effectiveness ratio (ICER) was € 19,302 per life year saved and € 30,291 per QALY. A probabilistic sensitivity analysis with 10,000 draws resulted in average ICERs of € 22,118 per life year and € 34,841 per QALY.

CONCLUSION

We provide evidence that lung cancer screening for a high-risk population may be more effective, but also more costly, than standard clinical care from the perspective of a German payer.

Improving model prediction of soil N2O emissions through Bayesian calibration

The biogeochemical processes that lead to the production of N₂O in arable soils are controlled by temporally and spatially varying drivers. The need for prediction of soil N₂O emissions across scales means that agroecosystem biogeochemistry models are widely used to simulate N₂O emissions. Due to the parameter-dense nature of agroecosystem models their parameters have to be calibrated according to the soil and climatic conditions of the intended area of application. Bayesian calibration is considered one of the most advanced ways to complete this task. In this study, we calibrate nine parameters of the Landscape-DNDC process-based agroecosystem model, which are key to its N₂O prediction. The Metropolis-Hastings algorithm is used at four separate implementations in order to estimate parameter posterior distributions at four arable sites in the UK. The results of this process are visualised, summarised and assessed against measured N₂O data from ten independent arable sites. The study shows that, in many cases, soil N₂O emission peaks that were not predicted with the default model parameters were predicted after calibration. Overall, the prediction of soil N₂O fluxes across all the sites that were considered was improved by 33% when using the calibrated parameters.

Carpals and epiphyses of radius and ulna as age indicators using longitudinal data: a Bayesian approach

The aim of this study is to develop a new formula for age estimation in a longitudinal study of a sample from the radiological collection of wrist bones of growing infants, children, and adolescents recorded at the Burlington Growth Centre. A sample of 82 individuals (43 boys and 39 girls), aged between 3 and 16 years, were analyzed with a total of 623 X-rays of left hand-wrist bones by measuring the area of carpal bones and epiphyses of the ulna and radius (Bo) and carpal area (Ca). The intra-class correlation coefficient (ICC) and its 95% confidence interval were used to evaluate intra-observer agreement. Hierarchical Bayesian calibration has been adopted to exceed the bias deriving from the classical regression approach used for age estimation in forensic disciplines, since it tends to overestimate or underestimate the age of the individuals. Calibration distributions of the dataset obtained by the evaluation of BoCa (the ratio of Bo and Ca) suggested mean absolute errors (MAE) of 1.07 and 1.34 years in boys and girls, respectively. The mean interquartile range (MIQR) was 1.7 and 2.42 years in boys and girls, respectively. The respective bias of the estimates was β_ERR = - 0.025 and - 0.074. Furthermore, a correspondence between different BoCa values and estimated age with its standard deviation (SD) was calculated for boys and girls, respectively. In conclusion, the Bayesian calibration method appears to be suitable for assessing both age and its distribution in subadults, according to hand-wrist maturity. Furthermore, it can easily incorporate other age predictors, obtaining a distribution of the subjects with multivariate predictors.

Statistical emulation of landslide-induced tsunamis at the Rockall Bank, NE Atlantic

Statistical methods constitute a useful approach to understand and quantify the uncertainty that governs complex tsunami mechanisms. Numerical experiments may often have a high computational cost. This forms a limiting factor for performing uncertainty and sensitivity analyses, where numerous simulations are required. Statistical emulators, as surrogates of these simulators, can provide predictions of the physical process in a much faster and computationally inexpensive way. They can form a prominent solution to explore thousands of scenarios that would be otherwise numerically expensive and difficult to achieve. In this work, we build a statistical emulator of the deterministic codes used to simulate submarine sliding and tsunami generation at the Rockall Bank, NE Atlantic Ocean, in two stages. First we calibrate, against observations of the landslide deposits, the parameters used in the landslide simulations. This calibration is performed under a Bayesian framework using Gaussian Process (GP) emulators to approximate the landslide model, and the discrepancy function between model and observations. Distributions of the calibrated input parameters are obtained as a result of the calibration. In a second step, a GP emulator is built to mimic the coupled landslide-tsunami numerical process. The emulator propagates the uncertainties in the distributions of the calibrated input parameters inferred from the first step to the outputs. As a result, a quantification of the uncertainty of the maximum free surface elevation at specified locations is obtained.

Uncertain future soil carbon dynamics under global change predicted by models constrained by total carbon measurements

Pool-based carbon (C) models are widely applied to predict soil C dynamics under global change and infer underlying mechanisms. However, it is unclear about the credibility of model-predicted C pool size, decay rate (k), and/or microbial C use efficiency (e) as only data on bulked total C is usually available for model constraining. Using observing system simulation experiments (OSSE), we constrained a two-pool model using simulated data sets of total soil C dynamics under topical hypotheses on responses of soil C dynamics to warming and elevated CO₂ (i.e., global change scenarios). The results indicated that the model predicted great uncertainties in C pool size, k, and e under all global change scenarios, resulting in the difficulty to correctly infer the presupposed "real" values of those parameters that are used to generate the simulated total soil C for constraining the model. Furthermore, the model using the constrained parameters generated divergent future soil C dynamics. Compared with the predictions using the presupposed real parameters (i.e., the real future C dynamics), the percentage uncertainty in 100-yr predictions using the constrained parameters was up to 45% depending on global change scenarios and data availability for model-constraining. Such great uncertainty was mainly due to the high collinearity among the model parameters. Using pool-based models, we argue that soil C pool size, k, and/or e and their responses to global change have to be estimated explicitly and empirically, rather than through model-fitting, in order to accurately predict C dynamics and infer underlying mechanisms. The OSSE approach provides a powerful way to identify data requirement for the new generation of model development and test model performance.

Using data-driven agent-based models for forecasting emerging infectious diseases

Producing timely, well-informed and reliable forecasts for an ongoing epidemic of an emerging infectious disease is a huge challenge. Epidemiologists and policy makers have to deal with poor data quality, limited understanding of the disease dynamics, rapidly changing social environment and the uncertainty on effects of various interventions in place. Under this setting, detailed computational models provide a comprehensive framework for integrating diverse data sources into a well-defined model of disease dynamics and social behavior, potentially leading to better understanding and actions. In this paper, we describe one such agent-based model framework developed for forecasting the 2014-2015 Ebola epidemic in Liberia, and subsequently used during the Ebola forecasting challenge. We describe the various components of the model, the calibration process and summarize the forecast performance across scenarios of the challenge. We conclude by highlighting how such a data-driven approach can be refined and adapted for future epidemics, and share the lessons learned over the course of the challenge.

A new inverse regression model applied to radiation biodosimetry

Biological dosimetry based on chromosome aberration scoring in peripheral blood lymphocytes enables timely assessment of the ionizing radiation dose absorbed by an individual. Here, new Bayesian-type count data inverse regression methods are introduced for situations where responses are Poisson or two-parameter compound Poisson distributed. Our Poisson models are calculated in a closed form, by means of Hermite and negative binomial (NB) distributions. For compound Poisson responses, complete and simplified models are provided. The simplified models are also expressible in a closed form and involve the use of compound Hermite and compound NB distributions. Three examples of applications are given that demonstrate the usefulness of these methodologies in cytogenetic radiation biodosimetry and in radiotherapy. We provide R and SAS codes which reproduce these examples.

Uncertainty Analysis for Computationally Expensive Models with Multiple Outputs

Bayesian MCMC calibration and uncertainty analysis for computationally expensive models is implemented using the SOARS (Statistical and Optimization Analysis using Response Surfaces) methodology. SOARS uses a radial basis function interpolator as a surrogate, also known as an emulator or meta-model, for the logarithm of the posterior density. To prevent wasteful evaluations of the expensive model, the emulator is built only on a high posterior density region (HPDR), which is located by a global optimization algorithm. The set of points in the HPDR where the expensive model is evaluated is determined sequentially by the GRIMA algorithm described in detail in another paper but outlined here. Enhancements of the GRIMA algorithm were introduced to improve efficiency. A case study uses an eight-parameter SWAT2005 (Soil and Water Assessment Tool) model where daily stream flows and phosphorus concentrations are modeled for the Town Brook watershed which is part of the New York City water supply. A Supplemental Material file available online contains additional technical details and additional analysis of the Town Brook application.

Local Derivative-Free Approximation of Computationally Expensive Posterior Densities

Bayesian inference using Markov chain Monte Carlo (MCMC) is computationally prohibitive when the posterior density of interest, π, is computationally expensive to evaluate. We develop a derivative-free algorithm GRIMA to accurately approximate π by interpolation over its high-probability density (HPD) region, which is initially unknown. Our local approach reduces the waste of computational budget on approximation of π in the low-probability region, which is inherent in global experimental designs. However, estimation of the HPD region is nontrivial when derivatives of π are not available or are not informative about the shape of the HPD region. Without relying on derivatives, GRIMA iterates (a) sequential knot selection over the estimated HPD region of π to refine the surrogate posterior and (b) re-estimation of the HPD region using an MCMC sample from the updated surrogate density, which is inexpensive to obtain. GRIMA is applicable to approximation of general unnormalized posterior densities. To determine the range of tractable problem dimensions, we conduct simulation experiments on test densities with linear and nonlinear component-wise dependence, skewness, kurtosis and multimodality. Subsequently, we use GRIMA in a case study to calibrate a computationally intensive nonlinear regression model to real data from the Town Brook watershed. Supplemental materials for this article are available online.

Efficient Interpolation of Computationally Expensive Posterior Densities With Variable Parameter Costs

Markov chain Monte Carlo (MCMC) is nowadays a standard approach to numerical computation of integrals of the posterior density π of the parameter vector η. Unfortunately, Bayesian inference using MCMC is computationally intractable when the posterior density π is expensive to evaluate. In many such problems, it is possible to identify a minimal subvector β of η responsible for the expensive computation in the evaluation of π. We propose two approaches, DOSKA and INDA, that approximate π by interpolation in ways that exploit this computational structure to mitigate the curse of dimensionality. DOSKA interpolates π directly while INDA interpolates π indirectly by interpolating functions, for example, a regression function, upon which π depends. Our primary contribution is derivation of a Gaussian processes interpolant that provably improves over some of the existing approaches by reducing the effective dimension of the interpolation problem from dim(η) to dim(β). This allows a dramatic reduction of the number of expensive evaluations necessary to construct an accurate approximation of π when dim(η) is high but dim(β) is low. We illustrate the proposed approaches in a case study for a spatio-temporal linear model for air pollution data in the greater Boston area. Supplemental materials include proofs, details, and software implementation of the proposed procedures.