Two heads are better than one: current landscape of integrating QSP and machine learning : An ISoP QSP SIG white paper by the working group on the integration of quantitative systems pharmacology and machine learning.

Quantitative systems pharmacology (QSP) modeling is applied to address essential questions in drug development, such as the mechanism of action of a therapeutic agent and the progression of disease. Meanwhile, machine learning (ML) approaches also contribute to answering these questions via the analysis of multi-layer 'omics' data such as gene expression, proteomics, metabolomics, and high-throughput imaging. Furthermore, ML approaches can also be applied to aspects of QSP modeling. Both approaches are powerful tools and there is considerable interest in integrating QSP modeling and ML. So far, a few successful implementations have been carried out from which we have learned about how each approach can overcome unique limitations of the other. The QSP + ML working group of the International Society of Pharmacometrics QSP Special Interest Group was convened in September, 2019 to identify and begin realizing new opportunities in QSP and ML integration. The working group, which comprises 21 members representing 18 academic and industry organizations, has identified four categories of current research activity which will be described herein together with case studies of applications to drug development decision making. The working group also concluded that the integration of QSP and ML is still in its early stages of moving from evaluating available technical tools to building case studies. This paper reports on this fast-moving field and serves as a foundation for future codification of best practices.

Generative adversarial networks for construction of virtual populations of mechanistic models: simulations to study Omecamtiv Mecarbil action.

Biophysical models are increasingly used to gain mechanistic insights by fitting and reproducing experimental and clinical data. The inherent variability in the recorded datasets, however, presents a key challenge. In this study, we present a novel approach, which integrates mechanistic modeling and machine learning to analyze in vitro cardiac mechanics data and solve the inverse problem of model parameter inference. We designed a novel generative adversarial network (GAN) and employed it to construct virtual populations of cardiac ventricular myocyte models in order to study the action of Omecamtiv Mecarbil (OM), a positive cardiac inotrope. Populations of models were calibrated from mechanically unloaded myocyte shortening recordings obtained in experiments on rat myocytes in the presence and absence of OM. The GAN was able to infer model parameters while incorporating prior information about which model parameters OM targets. The generated populations of models reproduced variations in myocyte contraction recorded during in vitro experiments and provided improved understanding of OM's mechanism of action. Inverse mapping of the experimental data using our approach suggests a novel action of OM, whereby it modifies interactions between myosin and tropomyosin proteins. To validate our approach, the inferred model parameters were used to replicate other in vitro experimental protocols, such as skinned preparations demonstrating an increase in calcium sensitivity and a decrease in the Hill coefficient of the force-calcium (F-Ca) curve under OM action. Our approach thereby facilitated the identification of the mechanistic underpinnings of experimental observations and the exploration of different hypotheses regarding variability in this complex biological system.

Can endothelial hemoglobin-α regulate nitric oxide vasodilatory signaling?

We used mathematical modeling to investigate nitric oxide (NO)-dependent vasodilatory signaling in the arteriolar wall. Detailed continuum cellular models of calcium (Ca2+) dynamics and membrane electrophysiology in smooth muscle and endothelial cells (EC) were coupled with models of NO signaling and biotransport in an arteriole. We used this theoretical approach to examine the role of endothelial hemoglobin-α (Hbα) as a modulator of NO-mediated myoendothelial feedback, as previously suggested in Straub et al. (Nature 491: 473-477, 2012). The model considers enriched expression of inositol 1,4,5-triphosphate receptors (IP3Rs), endothelial nitric oxide synthase (eNOS) enzyme, Ca2+-activated potassium (KCa) channels and Hbα in myoendothelial projections (MPs) between the two cell layers. The model suggests that NO-mediated myoendothelial feedback is plausible if a significant percentage of eNOS is localized within or near the myoendothelial projection. Model results show that the ability of Hbα to regulate the myoendothelial feedback is conditional to its colocalization with eNOS near MPs at concentrations in the high nanomolar range (>0.2 µM or 24,000 molecules). Simulations also show that the effect of Hbα observed in in vitro experimental studies may overestimate its contribution in vivo, in the presence of blood perfusion. Thus, additional experimentation is required to quantify the presence and spatial distribution of Hbα in the EC, as well as to test that the strong effect of Hbα on NO signaling seen in vitro, translates also into a physiologically relevant response in vivo.NEW & NOTEWORTHY Mathematical modeling shows that although regulation of nitric oxide signaling by hemoglobin-α (Hbα) is plausible, it is conditional to its presence in significant concentrations colocalized with endothelial nitric oxide synthase in myoendothelial projections. Additional experimentation is required to test that the strong effect of Hbα seen in vitro translates into a physiologically relevant response in vivo.

Mathematical model for shear stress dependent NO and adenine nucleotide production from endothelial cells.

We developed a mass transport model for a parallel-plate flow chamber apparatus to predict the concentrations of nitric oxide (NO) and adenine nucleotides (ATP, ADP) produced by cultured endothelial cells (ECs) and investigated how the net rates of production, degradation, and mass transport for these three chemical species vary with changes in wall shear stress (τw). These simulations provide an improved understanding of experimental results obtained with parallel-plate flow chambers and allows quantitative analysis of the relationship between τw, adenine nucleotide concentrations, and NO produced by ECs. Experimental data obtained after altering ATP and ADP concentrations with apyrase were analyzed to quantify changes in the rate of NO production (RNO). The effects of different isoforms of apyrase on ATP and ADP concentrations and nucleotide-dependent changes in RNO could be predicted with the model. A decrease in ATP was predicted with apyrase, but an increase in ADP was simulated due to degradation of ATP. We found that a simple proportional relationship relating a component of RNO to the sum of ATP and ADP provided a close match to the fitted curve for experimentally measured changes in RNO with apyrase. Estimates for the proportionality constant ranged from 0.0067 to 0.0321 µM/s increase in RNO per nM nucleotide concentration, depending on which isoform of apyrase was modeled, with the largest effect of nucleotides on RNO at low τw (<6 dyn/cm(2)).

Stochastic model of endothelial TRPV4 calcium sparklets: effect of bursting and cooperativity on EDH.

We examined the endothelial transient receptor vanilloid 4 (TRPV4) channel's vasodilatory signaling using mathematical modeling. The model analyzes experimental data by Sonkusare and coworkers on TRPV4-induced endothelial Ca(2+) events (sparklets). A previously developed continuum model of an endothelial and a smooth muscle cell coupled through microprojections was extended to account for the activity of a TRPV4 channel cluster. Different stochastic descriptions for the TRPV4 channel flux were examined using finite-state Markov chains. The model also took into consideration recent evidence for the colocalization of intermediate-conductance calcium-activated potassium channels (IKCa) and TRPV4 channels near the microprojections. A single TRPV4 channel opening resulted in a stochastic localized Ca(2+) increase in a small region (i.e., few µm(2) area) close to the channel. We predict micromolar Ca(2+) increases lasting for the open duration of the channel sufficient for the activation of low-affinity endothelial KCa channels. Simulations of a cluster of four TRPV4 channels incorporating burst and cooperative gating kinetics provided quantal Ca(2+) increases (i.e., steps of fixed amplitude), similar to the experimentally observed Ca(2+) sparklets. These localized Ca(2+) events result in endothelium-derived hyperpolarization (and SMC relaxation), with magnitude that depends on event frequency. The gating characteristics (bursting, cooperativity) of the TRPV4 cluster enhance Ca(2+) spread and the distance of KCa channel activation. This may amplify the EDH response by the additional recruitment of distant KCa channels.

Novel and flexible parameter estimation methods for data-consistent inversion in mechanistic modelling.

Predictions for physical systems often rely upon knowledge acquired from ensembles of entities, e.g. ensembles of cells in biological sciences. For qualitative and quantitative analysis, these ensembles are simulated with parametric families of mechanistic models (MMs). Two classes of methodologies, based on Bayesian inference and population of models, currently prevail in parameter estimation for physical systems. However, in Bayesian analysis, uninformative priors for MM parameters introduce undesirable bias. Here, we propose how to infer parameters within the framework of stochastic inverse problems (SIPs), also termed data-consistent inversion, wherein the prior targets only uncertainties that arise due to MM non-invertibility. To demonstrate, we introduce new methods to solve SIPs based on rejection sampling, Markov chain Monte Carlo, and generative adversarial networks (GANs). In addition, to overcome limitations of SIPs, we reformulate SIPs based on constrained optimization and present a novel GAN to solve the constrained optimization problem.

Multiple factors influence calcium synchronization in arterial vasomotion.

The intercellular synchronization of spontaneous calcium (Ca(2+)) oscillations in individual smooth muscle cells is a prerequisite for vasomotion. A detailed mathematical model of Ca(2+) dynamics in rat mesenteric arteries shows that a number of synchronizing and desynchronizing pathways may be involved. In particular, Ca(2+)-dependent phospholipase C, the intercellular diffusion of inositol trisphosphate (IP(3), and to a lesser extent Ca(2+)), IP(3) receptors, diacylglycerol-activated nonselective cation channels, and Ca(2+)-activated chloride channels can contribute to synchronization, whereas large-conductance Ca(2+)-activated potassium channels have a desynchronizing effect. Depending on the contractile state and agonist concentrations, different pathways become predominant, and can be revealed by carefully inhibiting the oscillatory component of their total activity. The phase shift between the Ca(2+) and membrane potential oscillations can change, and thus electrical coupling through gap junctions can mediate either synchronization or desynchronization. The effect of the endothelium is highly variable because it can simultaneously enhance the intercellular coupling and affect multiple smooth muscle cell components. Here, we outline a system of increased complexity and propose potential synchronization mechanisms that need to be experimentally tested.

Modeling Ca2+ signaling in the microcirculation: intercellular communication and vasoreactivity.

A network of intracellular signaling pathways and complex intercellular interactions regulate calcium mobilization in vascular cells, arteriolar tone, and blood flow. Different endothelium-derived vasoreactive factors have been identified and the importance of myoendothelial communication in vasoreactivity is now well appreciated. The ability of many vascular networks to conduct signals upstream also is established. This phenomenon is critical for both short-term changes in blood perfusion as well as long-term adaptations of a vascular network. In addition, in a phenomenon termed vasomotion, arterioles often exhibit spontaneous oscillations in diameter. This is thought to improve tissue oxygenation and enhance blood flow. Experimentation has begun to reveal important aspects of the regulatory machinery and the significance of these phenomena for the regulation of local perfusion and oxygenation. Mathematical modeling can assist in elucidating the complex signaling mechanisms that participate in these phenomena. This review highlights some of the important experimental studies and relevant mathematical models that provide the current understanding of these mechanisms in vasoreactivity.

Neuronal population models reveal specific linear conductance controllers sufficient to rescue preclinical disease phenotypes.

Preclinical drug candidates are screened for their ability to ameliorate in vitro neuronal electrophysiology, and go/no-go decisions progress drugs to clinical trials based on population means across cells and animals. However, these measures do not mitigate clinical endpoint risk. Population-based modeling captures variability across multiple electrophysiological measures from healthy, disease, and drug phenotypes. We pursued optimizing therapeutic targets by identifying coherent sets of ion channel target modulations for recovering heterogeneous wild-type (WT) population excitability profiles from a heterogeneous Huntington's disease (HD) population. Our approach combines mechanistic simulations with population modeling of striatal neurons using evolutionary optimization algorithms to design 'virtual drugs'. We introduce efficacy metrics to score populations and rank virtual drug candidates. We found virtual drugs using heuristic approaches that performed better than single target modulators and standard classification methods. We compare a real drug to virtual candidates and demonstrate a novel in silico triaging method.

Model order reduction for left ventricular mechanics via congruency training.

Computational models of the cardiovascular system and specifically heart function are currently being investigated as analytic tools to assist medical practice and clinical trials. To achieve clinical utility, models should be able to assimilate the diagnostic multi-modality data available for each patient and generate consistent representations of the underlying cardiovascular physiology. While finite element models of the heart can naturally account for patient-specific anatomies reconstructed from medical images, optimizing the many other parameters driving simulated cardiac functions is challenging due to computational complexity. With the goal of streamlining parameter adaptation, in this paper we present a novel, multifidelity strategy for model order reduction of 3-D finite element models of ventricular mechanics. Our approach is centered around well established findings on the similarity between contraction of an isolated muscle and the whole ventricle. Specifically, we demonstrate that simple linear transformations between sarcomere strain (tension) and ventricular volume (pressure) are sufficient to reproduce global pressure-volume outputs of 3-D finite element models even by a reduced model with just a single myocyte unit. We further develop a procedure for congruency training of a surrogate low-order model from multi-scale finite elements, and we construct an example of parameter optimization based on medical images. We discuss how the presented approach might be employed to process large datasets of medical images as well as databases of echocardiographic reports, paving the way towards application of heart mechanics models in the clinical practice.

Global Sensitivity Analysis of Ventricular Myocyte Model-Derived Metrics for Proarrhythmic Risk Assessment.

Multiscale computational models of the heart are being extensively investigated for improved assessment of drug-induced torsades de pointes (TdP) risk, a fatal side effect of many drugs. Model-derived metrics such as action potential duration and net charge carried by ionic currents (qNet) have been proposed as potential candidates for TdP risk stratification after being tested on small datasets. Unlike purely statistical approaches, model-derived metrics are thought to provide mechanism-based classification. In particular, qNet has been recently proposed as a surrogate metric for early afterdepolarizations (EADs), which are known to be cellular triggers of TdP. Analysis of critical model components and of the ion channels that have major impact on model-derived metrics can lead to improvements in the confidence of the prediction. In this paper, we analyze large populations of virtual drugs to systematically examine the influence of different ion channels on model-derived metrics that have been proposed for proarrhythmic risk assessment. We demonstrate via global sensitivity analysis (GSA) that model-derived metrics are most sensitive to different sets of input parameters. Similarly, important differences in sensitivity to specific channel blocks are highlighted when classifying drugs into different risk categories by either qNet or a metric directly based on simulated EADs. In particular, the higher sensitivity of qNet to the block of the late sodium channel might explain why its classification accuracy is better than that of the EAD-based metric, as shown for a small set of known drugs. Our results highlight the need for a better mechanistic interpretation of promising metrics like qNet based on a formal analysis of models. GSA should, therefore, constitute an essential component of the in silico workflow for proarrhythmic risk assessment, as an improved understanding of the structure of model-derived metrics could increase confidence in model-predicted risk.

Effect of spatial heterogeneity and colocalization of eNOS and capacitative calcium entry channels on shear stress-induced NO production by endothelial cells: A modeling approach.

INTRODUCTION: Colocalization of endothelial nitric oxide synthase (eNOS) and capacitative Ca2+ entry (CCE) channels in microdomains such as cavaeolae in endothelial cells (ECs) has been shown to significantly affect intracellular Ca2+ dynamics and NO production, but the effect has not been well quantified. METHODS: We developed a two-dimensional continuum model of an EC integrating shear stress-mediated ATP production, intracellular Ca2+ mobilization, and eNOS activation to investigate the effects of spatial colocalization of plasma membrane eNOS and CCE channels on Ca2+ dynamics and NO production in response to flow-induced shear stress. Our model examines the hypothesis that subcellular colocalization of cellular components can be critical for optimal coupling of NO production to blood flow. RESULTS: Our simulations predict that heterogeneity of CCE can result in formation of microdomains with significantly higher Ca2+ compared to the average cytosolic Ca2+. Ca2+ buffers with lower or no mobility further enhanced Ca2+ gradients relative to mobile buffers. Colocalization of eNOS to CCE channels significantly increased NO production. CONCLUSIONS: Our results provide quantitative understanding for the role of spatial heterogeneity and the compartmentalization of signals in regulation of shear stress-induced NO production.

Novel Two-Step Classifier for Torsades de Pointes Risk Stratification from Direct Features.

While pre-clinical Torsades de Pointes (TdP) risk classifiers had initially been based on drug-induced block of hERG potassium channels, it is now well established that improved risk prediction can be achieved by considering block of non-hERG ion channels. The current multi-channel TdP classifiers can be categorized into two classes. First, the classifiers that take as input the values of drug-induced block of ion channels (direct features). Second, the classifiers that are built on features extracted from output of the drug-induced multi-channel blockage simulations in the in-silico models (derived features). The classifiers built on derived features have thus far not consistently provided increased prediction accuracies, and hence casts doubt on the value of such approaches given the cost of including biophysical detail. Here, we propose a new two-step method for TdP risk classification, referred to as Multi-Channel Blockage at Early After Depolarization (MCB@EAD). In the first step, we classified the compound that produced insufficient hERG block as non-torsadogenic. In the second step, the role of non-hERG channels to modulate TdP risk are considered by constructing classifiers based on direct or derived features at critical hERG block concentrations that generates EADs in the computational cardiac cell models. MCB@EAD provides comparable or superior TdP risk classification of the drugs from the direct features in tests against published methods. TdP risk for the drugs highly correlated to the propensity to generate EADs in the model. However, the derived features of the biophysical models did not improve the predictive capability for TdP risk assessment.