Understanding the impact of numerical solvers on inference for differential equation models.

Most ordinary differential equation (ODE) models used to describe biological or physical systems must be solved approximately using numerical methods. Perniciously, even those solvers that seem sufficiently accurate for the forward problem, i.e. for obtaining an accurate simulation, might not be sufficiently accurate for the inverse problem, i.e. for inferring the model parameters from data. We show that for both fixed step and adaptive step ODE solvers, solving the forward problem with insufficient accuracy can distort likelihood surfaces, which might become jagged, causing inference algorithms to get stuck in local 'phantom' optima. We demonstrate that biases in inference arising from numerical approximation of ODEs are potentially most severe in systems involving low noise and rapid nonlinear dynamics. We reanalyse an ODE change point model previously fit to the COVID-19 outbreak in Germany and show the effect of the step size on simulation and inference results. We then fit a more complicated rainfall run-off model to hydrological data and illustrate the importance of tuning solver tolerances to avoid distorted likelihood surfaces. Our results indicate that, when performing inference for ODE model parameters, adaptive step size solver tolerances must be set cautiously and likelihood surfaces should be inspected for characteristic signs of numerical issues.

Simulating clinical trials for model-informed precision dosing: using warfarin treatment as a use case.

Treatment response variability across patients is a common phenomenon in clinical practice. For many drugs this inter-individual variability does not require much (if any) individualisation of dosing strategies. However, for some drugs, including chemotherapies and some monoclonal antibody treatments, individualisation of dosages are needed to avoid harmful adverse events. Model-informed precision dosing (MIPD) is an emerging approach to guide the individualisation of dosing regimens of otherwise difficult-to-administer drugs. Several MIPD approaches have been suggested to predict dosing strategies, including regression, reinforcement learning (RL) and pharmacokinetic and pharmacodynamic (PKPD) modelling. A unified framework to study the strengths and limitations of these approaches is missing. We develop a framework to simulate clinical MIPD trials, providing a cost and time efficient way to test different MIPD approaches. Central for our framework is a clinical trial model that emulates the complexities in clinical practice that challenge successful treatment individualisation. We demonstrate this framework using warfarin treatment as a use case and investigate three popular MIPD methods: 1. Neural network regression; 2. Deep RL; and 3. PKPD modelling. We find that the PKPD model individualises warfarin dosing regimens with the highest success rate and the highest efficiency: 75.1% of the individuals display INRs inside the therapeutic range at the end of the simulated trial; and the median time in the therapeutic range (TTR) is 74%. In comparison, the regression model and the deep RL model have success rates of 47.0% and 65.8%, and median TTRs of 45% and 68%. We also find that the MIPD models can attain different degrees of individualisation: the Regression model individualises dosing regimens up to variability explained by covariates; the Deep RL model and the PKPD model individualise dosing regimens accounting also for additional variation using monitoring data. However, the Deep RL model focusses on control of the treatment response, while the PKPD model uses the data also to further the individualisation of predictions.

Model-driven optimal experimental design for calibrating cardiac electrophysiology models.

BACKGROUND AND OBJECTIVE: Models of the cardiomyocyte action potential have contributed immensely to the understanding of heart function, pathophysiology, and the origin of heart rhythm disturbances. However, action potential models are highly nonlinear, making them difficult to parameterise and limiting to describing 'average cell' dynamics, when cell-specific models would be ideal to uncover inter-cell variability but are too experimentally challenging to be achieved. Here, we focus on automatically designing experimental protocols that allow us to better identify cell-specific maximum conductance values for each major current type. METHODS AND RESULTS: We developed an approach that applies optimal experimental designs to patch-clamp experiments, including both voltage-clamp and current-clamp experiments. We assessed the models calibrated to these new optimal designs by comparing them to the models calibrated to some of the commonly used designs in the literature. We showed that optimal designs are not only overall shorter in duration but also able to perform better than many of the existing experiment designs in terms of identifying model parameters and hence model predictive power. CONCLUSIONS: For cardiac cellular electrophysiology, this approach will allow researchers to define their hypothesis of the dynamics of the system and automatically design experimental protocols that will result in theoretically optimal designs.

Filter inference: A scalable nonlinear mixed effects inference approach for snapshot time series data.

Variability is an intrinsic property of biological systems and is often at the heart of their complex behaviour. Examples range from cell-to-cell variability in cell signalling pathways to variability in the response to treatment across patients. A popular approach to model and understand this variability is nonlinear mixed effects (NLME) modelling. However, estimating the parameters of NLME models from measurements quickly becomes computationally expensive as the number of measured individuals grows, making NLME inference intractable for datasets with thousands of measured individuals. This shortcoming is particularly limiting for snapshot datasets, common e.g. in cell biology, where high-throughput measurement techniques provide large numbers of single cell measurements. We introduce a novel approach for the estimation of NLME model parameters from snapshot measurements, which we call filter inference. Filter inference uses measurements of simulated individuals to define an approximate likelihood for the model parameters, avoiding the computational limitations of traditional NLME inference approaches and making efficient inferences from snapshot measurements possible. Filter inference also scales well with the number of model parameters, using state-of-the-art gradient-based MCMC algorithms such as the No-U-Turn Sampler (NUTS). We demonstrate the properties of filter inference using examples from early cancer growth modelling and from epidermal growth factor signalling pathway modelling.

Importance of modelling hERG binding in predicting drug-induced action potential prolongations for drug safety assessment.

Reduction of the rapid delayed rectifier potassium current (I Kr) via drug binding to the human Ether-à-go-go-Related Gene (hERG) channel is a well recognised mechanism that can contribute to an increased risk of Torsades de Pointes. Mathematical models have been created to replicate the effects of channel blockers, such as reducing the ionic conductance of the channel. Here, we study the impact of including state-dependent drug binding in a mathematical model of hERG when translating hERG inhibition to action potential changes. We show that the difference in action potential predictions when modelling drug binding of hERG using a state-dependent model versus a conductance scaling model depends not only on the properties of the drug and whether the experiment achieves steady state, but also on the experimental protocols. Furthermore, through exploring the model parameter space, we demonstrate that the state-dependent model and the conductance scaling model generally predict different action potential prolongations and are not interchangeable, while at high binding and unbinding rates, the conductance scaling model tends to predict shorter action potential prolongations. Finally, we observe that the difference in simulated action potentials between the models is determined by the binding and unbinding rate, rather than the trapping mechanism. This study demonstrates the importance of modelling drug binding and highlights the need for improved understanding of drug trapping which can have implications for the uses in drug safety assessment.

A Bayesian nonparametric method for detecting rapid changes in disease transmission.

Whether an outbreak of infectious disease is likely to grow or dissipate is determined through the time-varying reproduction number, Rt. Real-time or retrospective identification of changes in Rt following the imposition or relaxation of interventions can thus contribute important evidence about disease transmission dynamics which can inform policymaking. Here, we present a method for estimating shifts in Rt within a renewal model framework. Our method, which we call EpiCluster, is a Bayesian nonparametric model based on the Pitman-Yor process. We assume that Rt is piecewise-constant, and the incidence data and priors determine when or whether Rt should change and how many times it should do so throughout the series. We also introduce a prior which induces sparsity over the number of changepoints. Being Bayesian, our approach yields a measure of uncertainty in Rt and its changepoints. EpiCluster is fast, straightforward to use, and we demonstrate that it provides automated detection of rapid changes in transmission, either in real-time or retrospectively, for synthetic data series where the Rt profile is known. We illustrate the practical utility of our method by fitting it to case data of outbreaks of COVID-19 in Australia and Hong Kong, where it finds changepoints coinciding with the imposition of non-pharmaceutical interventions. Bayesian nonparametric methods, such as ours, allow the volume and complexity of the data to dictate the number of parameters required to approximate the process and should find wide application in epidemiology. This manuscript was submitted as part of a theme issue on "Modelling COVID-19 and Preparedness for Future Pandemics".

Models of the cardiac L-type calcium current: A quantitative review.

The L-type calcium current ( I CaL ) plays a critical role in cardiac electrophysiology, and models of I CaL are vital tools to predict arrhythmogenicity of drugs and mutations. Five decades of measuring and modeling I CaL have resulted in several competing theories (encoded in mathematical equations). However, the introduction of new models has not typically been accompanied by a data-driven critical comparison with previous work, so that it is unclear which model is best suited for any particular application. In this review, we describe and compare 73 published mammalian I CaL models and use simulated experiments to show that there is a large variability in their predictions, which is not substantially diminished when grouping by species or other categories. We provide model code for 60 models, list major data sources, and discuss experimental and modeling work that will be required to reduce this huge list of competing theories and ultimately develop a community consensus model of I CaL . This article is categorized under: Cardiovascular Diseases > Computational Models Cardiovascular Diseases > Molecular and Cellular Physiology.

A Parameter Representing Missing Charge Should Be Considered when Calibrating Action Potential Models.

Computational models of the electrical potential across a cell membrane are longstanding and vital tools in electrophysiology research and applications. These models describe how ionic currents, internal fluxes, and buffering interact to determine membrane voltage and form action potentials (APs). Although this relationship is usually expressed as a differential equation, previous studies have shown it can be rewritten in an algebraic form, allowing direct calculation of membrane voltage. Rewriting in this form requires the introduction of a new parameter, called Γ0 in this manuscript, which represents the net concentration of all charges that influence membrane voltage but are not considered in the model. Although several studies have examined the impact of Γ0 on long-term stability and drift in model predictions, there has been little examination of its effects on model predictions, particularly when a model is refit to new data. In this study, we illustrate how Γ0 affects important physiological properties such as action potential duration restitution, and examine the effects of (in)correctly specifying Γ0 during model calibration. We show that, although physiologically plausible, the range of concentrations used in popular models leads to orders of magnitude differences in Γ0, which can lead to very different model predictions. In model calibration, we find that using an incorrect value of Γ0 can lead to biased estimates of the inferred parameters, but that the predictive power of these models can be restored by fitting Γ0 as a separate parameter. These results show the value of making Γ0 explicit in model formulations, as it forces modellers and experimenters to consider the effects of uncertainty and potential discrepancy in initial concentrations upon model predictions.

Learning transmission dynamics modelling of COVID-19 using comomodels.

The COVID-19 epidemic continues to rage in many parts of the world. In the UK alone, an array of mathematical models have played a prominent role in guiding policymaking. Whilst considerable pedagogical material exists for understanding the basics of transmission dynamics modelling, there is a substantial gap between the relatively simple models used for exposition of the theory and those used in practice to model the transmission dynamics of COVID-19. Understanding these models requires considerable prerequisite knowledge and presents challenges to those new to the field of epidemiological modelling. In this paper, we introduce an open-source R package, comomodels, which can be used to understand the complexities of modelling the transmission dynamics of COVID-19 through a series of differential equation models. Alongside the base package, we describe a host of learning resources, including detailed tutorials and an interactive web-based interface allowing dynamic investigation of the model properties. We then use comomodels to illustrate three key lessons in the transmission of COVID-19 within R Markdown vignettes.

A Quantitative Systems Pharmacology Perspective on the Importance of Parameter Identifiability.

There is an inherent tension in Quantitative Systems Pharmacology (QSP) between the need to incorporate mathematical descriptions of complex physiology and drug targets with the necessity of developing robust, predictive and well-constrained models. In addition to this, there is no "gold standard" for model development and assessment in QSP. Moreover, there can be confusion over terminology such as model and parameter identifiability; complex and simple models; virtual populations; and other concepts, which leads to potential miscommunication and misapplication of methodologies within modeling communities, both the QSP community and related disciplines. This perspective article highlights the pros and cons of using simple (often identifiable) vs. complex (more physiologically detailed but often non-identifiable) models, as well as aspects of parameter identifiability, sensitivity and inference methodologies for model development and analysis. The paper distills the central themes of the issue of identifiability and optimal model size and discusses open challenges.

Using the Barthel Index and modified Rankin Scale as Outcome Measures for Stroke Rehabilitation Trials; A Comparison of Minimum Sample Size Requirements.

OBJECTIVES: Underpowered trials risk inaccurate results. Recruitment to stroke rehabilitation randomised controlled trials (RCTs) is often a challenge. Statistical simulations offer an important opportunity to explore the adequacy of sample sizes in the context of specific outcome measures. We aimed to examine and compare the adequacy of stroke rehabilitation RCT sample sizes using the Barthel Index (BI) or modified Rankin Scale (mRS) as primary outcomes. METHODS: We conducted computer simulations using typical experimental event rates (EER) and control event rates (CER) based on individual participant data (IPD) from stroke rehabilitation RCTs. Event rates are the proportion of participants who experienced clinically relevant improvements in the RCT experimental and control groups. We examined minimum sample size requirements and estimated the number of participants required to achieve a number needed to treat within clinically acceptable boundaries for the BI and mRS. RESULTS: We secured 2350 IPD (18 RCTs). For a 90% chance of statistical accuracy on the BI a rehabilitation RCT would require 273 participants per randomised group. Accurate interpretation of effect sizes would require 1000s of participants per group. Simulations for the mRS were not possible as a clinically relevant improvement was not detected when using this outcome measure. CONCLUSIONS: Stroke rehabilitation RCTs with large sample sizes are required for accurate interpretation of effect sizes based on the BI. The mRS lacked sensitivity to detect change and thus may be unsuitable as a primary outcome in stroke rehabilitation trials.

A Voltammetric Perspective of Multi-Electron and Proton Transfer in Protein Redox Chemistry: Insights From Computational Analysis of Escherichia coli HypD Fourier Transformed Alternating Current Voltammetry.

This paper explores the impact of pH on the mechanism of reversible disulfide bond (CysS-SCys) reductive breaking and oxidative formation in Escherichia coli hydrogenase maturation factor HypD, a protein which forms a highly stable adsorbed film on a graphite electrode. To achieve this, low frequency (8.96 Hz) Fourier transformed alternating current voltammetric (FTACV) experimental data was used in combination with modelling approaches based on Butler-Volmer theory with a dual polynomial capacitance model, utilizing an automated two-step fitting process conducted within a Bayesian framework. We previously showed that at pH 6.0 the protein data is best modelled by a redox reaction of two separate, stepwise one-electron, one-proton transfers with slightly "crossed" apparent reduction potentials that incorporate electron and proton transfer terms ( E app 2 0 > E app 1 0 ). Remarkably, rather than collapsing to a concerted two-electron redox reaction at more extreme pH, the same two-stepwise one-electron transfer model with E app 2 0 > E app 1 0 continues to provide the best fit to FTACV data measured across a proton concentration range from pH 4.0 to pH 9.0. A similar, small level of crossover in reversible potentials is also displayed in overall two-electron transitions in other proteins and enzymes, and this provides access to a small but finite amount of the one electron reduced intermediate state.

Barely sufficient practices in scientific computing.

The importance of software to modern research is well understood, as is the way in which software developed for research can support or undermine important research principles of findability, accessibility, interoperability, and reusability (FAIR). We propose a minimal subset of common software engineering principles that enable FAIRness of computational research and can be used as a baseline for software engineering in any research discipline.

Recent advances and future perspectives for automated parameterisation, Bayesian inference and machine learning in voltammetry.

Advanced data analysis tools such as mathematical optimisation, Bayesian inference and machine learning have the capability to revolutionise the field of quantitative voltammetry. Nowadays such approaches can be implemented routinely with widely available, user-friendly modern computing languages, algorithms and high speed computing to provide accurate and robust methods for quantitative comparison of experimental data with extensive simulated data sets derived from models proposed to describe complex electrochemical reactions. While the methodology is generic to all forms of dynamic electrochemistry, including the widely used direct current cyclic voltammetry, this review highlights advances achievable in the parameterisation of large amplitude alternating current voltammetry. One significant advantage this technique offers in terms of data analysis is that Fourier transformation provides access to the higher order harmonics that are almost devoid of background current. Perspectives on the technical advances needed to develop intelligent data analysis strategies and make them generally available to users of voltammetry are provided.

Using Purely Sinusoidal Voltammetry for Rapid Inference of Surface-Confined Electrochemical Reaction Parameters.

Alternating current (AC) voltammetric techniques are experimentally powerful as they enable Faradaic current to be isolated from non-Faradaic contributions. Finding the best global fit between experimental voltammetric data and simulations based on reaction models requires searching a substantial parameter space at high resolution. In this paper, we estimate parameters from purely sinusoidal voltammetry (PSV) experiments, investigating the redox reactions of a surface-confined ferrocene derivative. The advantage of PSV is that a complete experiment can be simulated relatively rapidly, compared to other AC voltammetric techniques. In one example involving thermodynamic dispersion, a PSV parameter inference effort requiring 7,500,000 simulations was completed in 7 h, whereas the same process for our previously used technique, ramped Fourier transform AC voltammetry (ramped FTACV), would have taken 4 days. Using both synthetic and experimental data with a surface confined diazonium substituted ferrocene derivative, it is shown that the PSV technique can be used to recover the key chemical and physical parameters. By applying techniques from Bayesian inference and Markov chain Monte Carlo methods, the confidence, distribution, and degree of correlation of the recovered parameters was visualized and quantified.

A Monte Carlo method to estimate cell population heterogeneity from cell snapshot data.

Variation is characteristic of all living systems. Laboratory techniques such as flow cytometry can probe individual cells, and, after decades of experimentation, it is clear that even members of genetically identical cell populations can exhibit differences. To understand whether variation is biologically meaningful, it is essential to discern its source. Mathematical models of biological systems are tools that can be used to investigate causes of cell-to-cell variation. From mathematical analysis and simulation of these models, biological hypotheses can be posed and investigated, then parameter inference can determine which of these is compatible with experimental data. Data from laboratory experiments often consist of "snapshots" representing distributions of cellular properties at different points in time, rather than individual cell trajectories. These data are not straightforward to fit using hierarchical Bayesian methods, which require the number of cell population clusters to be chosen a priori. Nor are they amenable to standard nonlinear mixed effect methods, since a single observation per cell is typically too few to estimate parameter variability. Here, we introduce a computational sampling method named "Contour Monte Carlo" (CMC) for estimating mathematical model parameters from snapshot distributions, which is straightforward to implement and does not require that cells be assigned to predefined categories. The CMC algorithm fits to snapshot probability distributions rather than raw data, which means its computational burden does not, like existing approaches, increase with the number of cells observed. Our method is appropriate for underdetermined systems, where there are fewer distinct types of observations than parameters to be determined, and where observed variation is mostly due to variability in cellular processes rather than experimental measurement error. This may be the case for many systems due to continued improvements in resolution of laboratory techniques. In this paper, we apply our method to quantify cellular variation for three biological systems of interest and provide Julia code enabling others to use this method.

Accounting for variability in ion current recordings using a mathematical model of artefacts in voltage-clamp experiments.

Mathematical models of ion channels, which constitute indispensable components of action potential models, are commonly constructed by fitting to whole-cell patch-clamp data. In a previous study, we fitted cell-specific models to hERG1a (Kv11.1) recordings simultaneously measured using an automated high-throughput system, and studied cell-cell variability by inspecting the resulting model parameters. However, the origin of the observed variability was not identified. Here, we study the source of variability by constructing a model that describes not just ion current dynamics, but the entire voltage-clamp experiment. The experimental artefact components of the model include: series resistance, membrane and pipette capacitance, voltage offsets, imperfect compensations made by the amplifier for these phenomena, and leak current. In this model, variability in the observations can be explained by either cell properties, measurement artefacts, or both. Remarkably, by assuming that variability arises exclusively from measurement artefacts, it is possible to explain a larger amount of the observed variability than when assuming cell-specific ion current kinetics. This assumption also leads to a smaller number of model parameters. This result suggests that most of the observed variability in patch-clamp data measured under the same conditions is caused by experimental artefacts, and hence can be compensated for in post-processing by using our model for the patch-clamp experiment. This study has implications for the question of the extent to which cell-cell variability in ion channel kinetics exists, and opens up routes for better correction of artefacts in patch-clamp data. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

Chaste: Cancer, Heart and Soft Tissue Environment.

Chaste (Cancer, Heart And Soft Tissue Environment) is an open source simulation package for the numerical solution of mathematical models arising in physiology and biology. To date, Chaste development has been driven primarily by applications that include continuum modelling of cardiac electrophysiology ('Cardiac Chaste'), discrete cell-based modelling of soft tissues ('Cell-based Chaste'), and modelling of ventilation in lungs ('Lung Chaste'). Cardiac Chaste addresses the need for a high-performance, generic, and verified simulation framework for cardiac electrophysiology that is freely available to the scientific community. Cardiac chaste provides a software package capable of realistic heart simulations that is efficient, rigorously tested, and runs on HPC platforms. Cell-based Chaste addresses the need for efficient and verified implementations of cell-based modelling frameworks, providing a set of extensible tools for simulating biological tissues. Computational modelling, along with live imaging techniques, plays an important role in understanding the processes of tissue growth and repair. A wide range of cell-based modelling frameworks have been developed that have each been successfully applied in a range of biological applications. Cell-based Chaste includes implementations of the cellular automaton model, the cellular Potts model, cell-centre models with cell representations as overlapping spheres or Voronoi tessellations, and the vertex model. Lung Chaste addresses the need for a novel, generic and efficient lung modelling software package that is both tested and verified. It aims to couple biophysically-detailed models of airway mechanics with organ-scale ventilation models in a package that is freely available to the scientific community. Chaste is designed to be modular and extensible, providing libraries for common scientific computing infrastructure such as linear algebra operations, finite element meshes, and ordinary and partial differential equation solvers. This infrastructure is used by libraries for specific applications, such as continuum mechanics, cardiac models, and cell-based models. The software engineering techniques used to develop Chaste are intended to ensure code quality, re-usability and reliability. Primary applications of the software include cardiac and respiratory physiology, cancer and developmental biology.

Four Ways to Fit an Ion Channel Model.

Mathematical models of ionic currents are used to study the electrophysiology of the heart, brain, gut, and several other organs. Increasingly, these models are being used predictively in the clinic, for example, to predict the risks and results of genetic mutations, pharmacological treatments, or surgical procedures. These safety-critical applications depend on accurate characterization of the underlying ionic currents. Four different methods can be found in the literature to fit voltage-sensitive ion channel models to whole-cell current measurements: method 1, fitting model equations directly to time-constant, steady-state, and I-V summary curves; method 2, fitting by comparing simulated versions of these summary curves to their experimental counterparts; method 3, fitting to the current traces themselves from a range of protocols; and method 4, fitting to a single current trace from a short and rapidly fluctuating voltage-clamp protocol. We compare these methods using a set of experiments in which hERG1a current was measured in nine Chinese hamster ovary cells. In each cell, the same sequence of fitting protocols was applied, as well as an independent validation protocol. We show that methods 3 and 4 provide the best predictions on the independent validation set and that short, rapidly fluctuating protocols like that used in method 4 can replace much longer conventional protocols without loss of predictive ability. Although data for method 2 are most readily available from the literature, we find it performs poorly compared to methods 3 and 4 both in accuracy of predictions and computational efficiency. Our results demonstrate how novel experimental and computational approaches can improve the quality of model predictions in safety-critical applications.

Rapid Characterization of hERG Channel Kinetics I: Using an Automated High-Throughput System.

Predicting how pharmaceuticals may affect heart rhythm is a crucial step in drug development and requires a deep understanding of a compound's action on ion channels. In vitro hERG channel current recordings are an important step in evaluating the proarrhythmic potential of small molecules and are now routinely performed using automated high-throughput patch-clamp platforms. These machines can execute traditional voltage-clamp protocols aimed at specific gating processes, but the array of protocols needed to fully characterize a current is typically too long to be applied in a single cell. Shorter high-information protocols have recently been introduced that have this capability, but they are not typically compatible with high-throughput platforms. We present a new 15 second protocol to characterize hERG (Kv11.1) kinetics, suitable for both manual and high-throughput systems. We demonstrate its use on the Nanion SyncroPatch 384PE, a 384-well automated patch-clamp platform, by applying it to Chinese hamster ovary cells stably expressing hERG1a. From these recordings, we construct 124 cell-specific variants/parameterizations of a hERG model at 25°C. A further eight independent protocols are run in each cell and are used to validate the model predictions. We then combine the experimental recordings using a hierarchical Bayesian model, which we use to quantify the uncertainty in the model parameters, and their variability from cell-to-cell; we use this model to suggest reasons for the variability. This study demonstrates a robust method to measure and quantify uncertainty and shows that it is possible and practical to use high-throughput systems to capture full hERG channel kinetics quantitatively and rapidly.