A Heart for Diversity: Simulating Variability in Cardiac Arrhythmia Research.

In cardiac electrophysiology, there exist many sources of inter- and intra-personal variability. These include variability in conditions and environment, and genotypic and molecular diversity, including differences in expression and behavior of ion channels and transporters, which lead to phenotypic diversity (e.g., variable integrated responses at the cell, tissue, and organ levels). These variabilities play an important role in progression of heart disease and arrhythmia syndromes and outcomes of therapeutic interventions. Yet, the traditional in silico framework for investigating cardiac arrhythmias is built upon a parameter/property-averaging approach that typically overlooks the physiological diversity. Inspired by work done in genetics and neuroscience, new modeling frameworks of cardiac electrophysiology have been recently developed that take advantage of modern computational capabilities and approaches, and account for the variance in the biological data they are intended to illuminate. In this review, we outline the recent advances in statistical and computational techniques that take into account physiological variability, and move beyond the traditional cardiac model-building scheme that involves averaging over samples from many individuals in the construction of a highly tuned composite model. We discuss how these advanced methods have harnessed the power of big (simulated) data to study the mechanisms of cardiac arrhythmias, with a special emphasis on atrial fibrillation, and improve the assessment of proarrhythmic risk and drug response. The challenges of using in silico approaches with variability are also addressed and future directions are proposed.

Introduction

Beginning with the seminal paper by Hodgkin and Huxley, 1952, mathematical models of electrophysiology have proven to be valuable tools for better understanding many physiological processes, especially in cardiac arrhythmia research (Noble et al., 2012; Dibb et al., 2015). Fifty-six years after publication of the first cardiac model (Noble, 1962), there is currently a computational model for almost every cell type of the heart, including nodal, atrial, ventricular, and Purkinje cells (Beeler and Reuter, 1977; Difrancesco and Noble, 1985; Luo and Rudy, 1991; Inada et al., 2009; Maltsev and Lakatta, 2009; Sampson et al., 2010; Grandi et al., 2011; O'Hara et al., 2011), for numerous species, and for various levels of complexity across multiple spatial scales (e.g., inclusion of disease-associated remodeling, drug action, contractile, energetics, and signaling modules) (Fink et al., 2011). Most of these models use average data from voltage-clamp experiments of individual ionic membrane currents, and while they have led to many important advances in studies of cardiac electrophysiology and pathology, especially cardiac arrhythmias (Sepulveda et al., 1989; Courtemanche and Winfree, 1991; Panfilov and Holden, 1991; Gray et al., 1995; Krogh-Madsen and Christini, 2012; Roberts et al., 2012; Bueno-Orovio et al., 2014), they typically represent the average behavior of a particular cell type. Because these models ignore evident inter- and intra-personal variability, they can fail to capture the properties of any given individual cell or cells in a given patient (Golowasch et al., 2002; Dokos and Lovell, 2004; Sarkar and Sobie, 2010; Marder, 2011; Zaniboni, 2011; Groenendaal et al., 2015; Pathmanathan et al., 2015). This is in part because incorporating variance that reflects biological data into cardiac models requires significant computational capacity, particularly as compared to what was available when mathematical modeling of electrophysiology first emerged. Given ever-increasing computational capabilities, new modeling approaches have been developed to reproduce and analyze the immense physiological diversity observed in electrophysiology.

In this review, we discuss the importance of accounting for variability when building models of cardiac electrophysiology in both physiological and diseased conditions, and describe new tools being developed to address the limitations of traditional modeling approaches. In particular, we focus on two computational approaches that have emerged as leading methodologies for studying variability in cardiac electrophysiology, which we will refer to as (1) population-based modeling and (2) sample-specific modeling. Both methodologies have provided valuable insights in the fields of neuroscience, cardiology, and pharmacology. Here we review how they have advanced our understanding of the basic mechanisms of cardiac arrhythmias, and particularly atrial fibrillation (AF). As these in silico approaches have led to important insights into arrhythmia risks, mechanisms of arrhythmogenesis, and variable response to drugs, they should be considered when determining the regulatory requirements for the proarrhythmia assessment and drug efficacy and safety evaluation of candidate compounds.

Importance of taking into account variability in cardiac electrophysiology

Sources of variability in cardiac electrophysiology encompass multiple spatial and temporal scales, and include inter-species (Sham et al., 1995; Su et al., 2003), inter-ethnic (Lau et al., 2012; Fender et al., 2014), inter-subject (Taneja et al., 2001; Batchvarov et al., 2002), regional (Feng et al., 1998; Yan et al., 1998), and temporal (Jeyaraj et al., 2012) differences. Variability in experimental electrophysiological data does not only represent natural physiological diversity, but also reflects, in part, differences in the experimental conditions in which electrophysiological measurements are performed (Groenendaal et al., 2015). These conditions can vary from laboratory to laboratory (Niederer et al., 2009; Fink et al., 2011), experiment to experiment, or during the same experiment, e.g., due to deterioration of the experimental preparations over time (Fink et al., 2011). There is also instrument noise (Mirams et al., 2016), artifacts, and use of non-physiological temperatures and solutions (Groenendaal et al., 2015), all of which impact the structure and function of the cellular or tissue components being studied. These sources of variation are difficult to control even for experienced electrophysiologists and are equally as challenging to account for by modelers.

Mathematical cardiomyocyte models have remained useful tools for unraveling physiological and pathophysiological mechanisms, including mechanisms of arrhythmia, and identifying antiarrhythmic strategies without accounting for variability (Sepulveda et al., 1989; Courtemanche and Winfree, 1991; Clancy and Rudy, 1999, 2002; Clancy et al., 2002; Rivolta et al., 2002; Gong et al., 2007; Noble et al., 2007; Tsujimae et al., 2007; Zhang et al., 2007; Zhu and Clancy, 2007; Campbell et al., 2008; Comtois et al., 2008; Kharche et al., 2008; Sale et al., 2008; Ahrens-Nicklas et al., 2009; Butters et al., 2010; Adeniran et al., 2012; Edwards et al., 2014; Grandi and Maleckar, 2016; Morotti et al., 2016; Ni et al., 2017). Although average models have also been successfully applied to the study of sources of variability, such as sexual and hormonal factors (Yang et al., 2010, 2017; Yang and Clancy, 2012), age (Behar and Yaniv, 2017), and circadian regulation (Fotiadis and Forger, 2013), the rationale for developing novel computational approaches that specifically account for electrophysiological variability can be summarized by two main reasons.

Average data may not accurately represent any specific individual or behavior well

The traditional cardiac model-building scheme involves averaging over samples from multiple experiments from many individuals, both to parameterize the model and validate it, which may not represent any specific measured physiological behavior very well. This “failure of averaging” has been demonstrated in many fields, most recently in neuroscience (Golowasch et al., 2002; Marder, 2011), and was particularly well-documented in 1952, the same year that the seminal Hodgkin and Huxley paper was published, when Lt. Gilbert S. Daniels of the U.S. Air Force published a Technical Note that highlighted the fundamental problem with fitting data to the mean (Daniels, 1952; Rose, 2017). Using data from 4,063 pilots, Lt. Daniels calculated the average of 10 physical dimensions believed to be most relevant for design of the cockpit on a plane, including a pilot's height, chest circumference, and sleeve length. Surprisingly, he found that a total of zero individuals fit within the middle 30% of the range of values for each dimension, and less than 3.5% of pilots would be average sized on any three dimensions. After this finding, the U.S. Air Force completely moved from standardizing all dimensions to an “average pilot,” to making all the dimensions adjustable to each individual pilot, which immediately and drastically improved performance and was soon adopted by all branches of the American military. Modeling of electrophysiology is undergoing a similar evolution, which will likely improve the translational significance of the models.

Variability has implications on genesis and treatment of arrhythmia

Variability plays an important role in arrhythmia generation and treatment, as exemplified by AF. The atria are characterized by a high degree of phenotypic variability in physiological properties, with broad and diverging distributions of biomarkers in patients in normal sinus rhythm (nSR) or chronic AF (cAF, Figure 1A) (Ravens et al., 2015), likely due to innate variability of the ionic currents (perhaps due to stochastic gating) that can affect whole cell and/or tissue proarrhythmic behavior (Pueyo et al., 2011; Heijman et al., 2013). This phenotypic variability can be captured by adding variability in the conductance parameters of a mathematical model of the action potential (AP, Figures 1B,C). In some circumstances, physiological variability itself can be the substrate for arrhythmia. For example, increased heterogeneity of refractoriness is important for the maintenance of AF (Moe et al., 1964; Boutjdir et al., 1986; Misier et al., 1992; Sato et al., 1992; Wang et al., 1995, 1996; Gaspo et al., 1997; Liu and Nattel, 1997; Ramirez et al., 2000), and regional differences in atrial ionic currents play a significant role in atrial arrhythmia initiation (Feng et al., 1998; Gaborit et al., 2007; Colman et al., 2013). Consequently, pharmacotherapy that increases dispersion of refractoriness is an adverse side effect of drugs for the treatment of AF (Ramanna et al., 2001; Soylu et al., 2003).

Figure 1

Variability in cardiac electrophysiology.(A) Histograms of AP duration at 90% repolarization (APD90) for patients in nSR (black) and cAF (red) show substantial variability. Reproduced from Ravens et al. (2015) with permission. (B) Example of APs and (C) histogram of APD90 produced using models incorporating variability in conductances of ionic currents; some models (in blue) are rejected due to non-physiological behaviors.

It is well-known that individuals may present largely different responses to same pharmacological interventions. As an example, it has been shown that drugs that block the hERG (human ether-à-go-go-related gene) channel are generally responsible for drug-induced long QT syndrome (diLQTS), but in a population this adverse response is highly variable, from minimum changes in the electrocardiogram (ECG) to induction of lethal ventricular arrhythmias (Singh et al., 2000; Kannankeril et al., 2011). Accounting for physiological variability may help better understand why some individuals display adverse side effects, while others do not. Given the different etiologies of many cardiac arrhythmias, such as AF, computational approaches that take into account variability may help us identify subpopulations in which a particular antiarrhythmic therapy will be effective and safe, or toxic. Furthermore, when assessing the efficacy and safety of a drug administration for heart conditions, it is important to take into account physiological and pathological variabilities to make sure that results are quantified and valid at the population level. Such approaches will potentially be more clinically useful in simulating the effects of drugs and aiding the design of safer and more effective therapies (Britton et al., 2013, 2017a; Passini et al., 2016; Yang et al., 2016).

Approaches and insight on the impact of variability on cardiac electrophysiology

Although many methods have been developed, two families of approaches have emerged as leading methodologies to account for variability in cardiac electrophysiology: (1) population-based and (2) sample-specific modeling (Figure 2). Both methods generally start with the building or use of a baseline cardiac cell model, which has been parameterized and validated to average data. Population-based approaches generate model variants of the baseline model that fit given experimental distributions of electrophysiological outcomes or biomarkers, while sample-specific modeling approaches re-parameterize the baseline model based on cell- or patient-specific datasets (Figure 2). Because their implementation requires computational power, the advancements in computing capabilities and techniques (Pitt-Francis et al., 2006; Abramson et al., 2010), especially in parallel computing (Wang et al., 2011), have helped these methods gain traction in the last decade.

Figure 2

Flowchart connecting traditional cardiac modeling approach to the new methodologies that account for variability.

Population-based modeling

Population-based modeling approaches have been developed and employed to obtain results at the population level, which led to many novel insights into physiological and pathophysiological variabilities, and variable responses to drug administration. We refer the readers to a recent review from the Rodriguez group (Muszkiewicz et al., 2016), where this methodology is described in detail. Here, we briefly describe the general approach of population-based modeling, and summarize how it has contributed to advancing the field as exemplified by some important studies.

Creating populations of AP models

Populations of models are generally created by modifying sets of parameters in a baseline model (Figure 2). This process involves determination of the parameters to be varied, over what range, and a sampling method to select the parameter values. Frequently, maximal conductances or maximum transport rates of ion channels, pumps and transporters in AP models are selected to vary. The parameter space over which these model parameters vary can be chosen either to reflect the experimental range, when available, or theoretical upper and lower bounds. Then, populations of parameter sets are created by sampling the parameters within the predefined parameter spaces. Typically, four types of sampling methods have been applied to obtain the populations of parameter sets: uniform-interval sampling (Romero et al., 2009, 2011; Corrias et al., 2011), log-normal sampling (Sobie, 2009; Sadrieh et al., 2013; Ellinwood et al., 2017a; Morotti and Grandi, 2017), Latin hypercube sampling (LHS) (Britton et al., 2013) and its variants such as orthogonal sampling (Burrage et al., 2015; Donovan et al., 2018), and sequential Monte Carlo sampling (Muszkiewicz et al., 2017).

After generating hundreds or thousands of model variants, calibration can be performed to exclude models that display non-physiological behaviors (Figure 2). This can be done, for example, by removing models showing repolarization failure (Sobie, 2009), or exhibiting AP duration (APD) more than three standard deviations from the population mean (Devenyi and Sobie, 2016). Population of models are also calibrated to measured data from patient samples, whereby model variants are selected based on simulated electrophysiological properties, such as APD, upstroke velocity, resting membrane potential and/or Ca2+ transient (CaT) (Britton et al., 2013, 2017a; Sanchez et al., 2014; Passini et al., 2016; Rees et al., 2018). Other studies use additional information such as ionic current data (Muszkiewicz et al., 2017), or ECG data (Mann et al., 2016). This calibration step is meant to ensure that (1) variants displaying non-physiological properties are discarded before analysis, and (2) the simulated electrophysiological properties of models in a given population are in the same range as experimental data, or, more recently, correspond to the same distribution of observed experimental biomarkers (Lawson et al., 2018), thus possibly making inferences from in silico experiments more physiologically relevant.

Analyzing populations of AP models

Once a population of cardiac AP models is generated, and electrophysiological simulations have been performed corresponding to the scientific question at hand, mechanistic insights can be obtained using various analysis techniques. These analyses have contributed to our understanding of the relative role of the underlying parameters in modulating the physiological properties of interest (i.e., sensitivity analysis), or revealing association of certain parameter ranges or properties with specific physiological behaviors (e.g., repolarization abnormalities, ectopic activity, drug response). Many relevant examples have recently been reviewed (Muszkiewicz et al., 2016). Here we highlight new recent developments and discuss details of parameter sensitivity analysis.

Performing parameter sensitivity analysis

A common systematic analysis of populations of models is sensitivity analysis. It assesses how model outputs, which typically represent whole cell behavior (e.g., APD), are sensitive to changes in model parameters, (e.g., conductances and maximum transport rates). Because many parameters are often varied to generate the populations of models, multivariable linear regression (Hair et al., 2010; Draper and Smith, 2014) has emerged as a frequently utilized tool to perform sensitivity analysis in cardiac electrophysiology. Moreover, as the number of independent parameters varied is used to predict a smaller set of dependent variables, sensitivity analysis is typically performed using partial least squares regression (Geladi and Kowalski, 1986; Sobie, 2009), as compared to standard multivariable regression. The result of linear regression is a set of coefficients (forming a “regression model”) describing how perturbing a particular parameter influences an output of interest. This method has been successfully utilized in other fields such as molecular biology (Janes et al., 2005) and neuroscience (Weaver and Wearne, 2008), and was first used in cardiac electrophysiology by Sobie (2009), who applied it to study sensitivities of properties such as APD and pacing rate threshold for inducing AP alternans. Since the regression model represents a linear approximation of a highly non-linear system, it is important to always check the goodness of fit. Several papers by the Sobie's group have indeed shown that the linear approximation actually provides a very good fit of the AP biomarkers, which was not trivially predictable (Sarkar et al., 2012).

The approach of varying multiple ionic conductances at once in a systematic fashion (as opposed to one at a time) and employing sensitivity analysis using multivariable regression has led to many important insights in cardiac electrophysiology (Sarkar and Sobie, 2011; Mann et al., 2012; Heijman et al., 2013; Walmsley et al., 2013), some of which have been confirmed experimentally (Lee et al., 2013; Devenyi and Sobie, 2016; Devenyi et al., 2017). For example, it has been used to study how different diseased conditions affect the sensitivities of given electrophysiological properties (Sadrieh et al., 2013; Ellinwood et al., 2017a; Vagos et al., 2017; Koivumaki et al., 2018), mechanisms of physiological phenomena (Lee et al., 2013), and for constraining free parameters (Sarkar and Sobie, 2010). Through sensitivity analysis, Cummins et al. identified multiple potential ionic targets mediating forward rate dependence (FRD) of AP, and demonstrated that modulation of the inward rectifier K+ current (IK1) or the Na+/K+ pump current was more likely to produce FRD (Cummins et al., 2014). Devenyi and Sobie performed sensitivity analysis of rat ventricular myocyte models, and quantitatively compared the modulatory role of transient outward K+ current (Ito) and sarco/endoplasmic reticulum Ca2+-ATPase (SERCA) in CaT amplitude. They found that in rat epicardial cells Ito plays a more important role than SERCA in regulating CaT amplitude, and this was analogous to human atrial cells, where both Ito and ultra-rapid delayed-rectifier K+ current (IKur) had greater impacts on CaT amplitude than did SERCA (Devenyi and Sobie, 2016). These studies highlight how sensitivity analysis can be applied to compare and contrast roles of different ionic processes and Ca2+ handling in regulating physiological properties and behaviors between cell types and species. Sensitivity analysis has also been used to compare the dependence of AP biomarkers on ionic processes in healthy and diseased conditions. For example, Lee et al. compared the impact of ionic processes on APD in control and AF-remodeled cells and found that the Na+/Ca2+ exchanger (NCX) current has little influence on APD in control cells but more markedly impacts AF cells; the analysis also revealed that IK1 upregulation plays a dominant role in APD shortening in AF, and that the L-type Ca2+ current (ICaL) significantly contributes to rate-dependent APD changes in both control and AF myocytes (Lee Y. S. et al., 2016). Most recently, Gong and Sobie described a novel use of regression models, cross-cell regression, to predict adult myocyte drug responses from induced pluripotent stem-cell-derived cardiomyocytes (iPSC-CMs) behaviors (Gong and Sobie, 2018).

Multivariable linear regression is used if the physiological output of interest is continuous, but, for the study of arrhythmia mechanisms, another particularly useful and efficient regression technique is logistic regression, which is used when the outcome of interest is Boolean (i.e., yes/no, true/false). Applying logistic regression in studies of physiology is well-described by Lee et al. who examined Ca2+ spark triggering (an all-or-none event), and demonstrated the accuracy of logistic regression using receiver operator characteristic curves (Lee et al., 2013). This method has since been used to study the probability that a certain arrhythmic event such as early afterdepolarizations (EADs) will occur and suggest underlying factors (Morotti and Grandi, 2017).

The main limitation of regression (both linear and logistic) analysis is that it only highlights how inputs are correlated to outputs, so the conclusions drawn from the analysis can be misleading if only a few outputs are considered. For example, it has been shown that completely different parameter combinations could produce essentially identical AP shapes but substantially different CaT amplitudes (Figure 3) (Sarkar and Sobie, 2010). However, sensitivity analysis can still help determine whether the relationships between inputs and outputs in computational models match experimental findings and assumptions, and whether there are particularly influential parameters that can be exploited therapeutically or targeted to better understand a given physiological phenomena (e.g., arrhythmia mechanism).

Figure 3

Different subcellular parameter combinations can lead to same AP shape. Example of how different model parameter combinations (e.g., ion channel conductances and maximum transport rates) can produce nearly identical atrial AP morphologies, but notably different CaTs.

Comparing subpopulations of models

Comparing subgroups in a population of models (often using statistical difference tests of parameters of interest) can help identify underlying determinants of different phenotypes, behaviors, and pathological conditions (Sanchez et al., 2014; Zhou et al., 2016; Britton et al., 2017b; Muszkiewicz et al., 2017; Vagos et al., 2017; Lawson et al., 2018). For example, through characterizing ionic parameters of models that are prone to repolarization abnormalities, Britton et al. found that the electrogenic Na+/K+ pump is a key determinant of susceptibility to repolarization abnormalities in human ventricular cardiomyocytes by applying arrhythmia-provoking conditions to a population of experimentally-calibrated cardiac cells (Britton et al., 2017b). A population-based approach has also been used to tease out the ionic mechanisms underlying variability in iPSC-CMs (Paci et al., 2017). By calibrating generated subpopulations of human atrial myocyte models to ranges of experimental data from a large number of patients with nSR or cAF, Sanchez et al. characterized potential ionic determinants of inter-subject variability in AP biomarkers, and identified similar changes in IK1, IKur, and Ito in cAF vs. nSR subpopulations that were consistent with experimentally reported AF-induced remodeling effects (Sanchez et al., 2014). In a more recent study, instead of calibrating population of models to the range of experimental dataset, Lawson et al. proposed a novel method to calibrate these models to the distributions of multiple experimentally measured biomarkers (Lawson et al., 2018), which led to an improved characterization of ionic differences between nSR and cAF. These studies focused on AP biomarkers at a fixed pacing rate. In a different study, Vagos et al. expanded the use of population of models to compare the steady state and dynamic restitution behaviors of AP in nSR and cAF populations (Vagos et al., 2017). By combining population-based modeling and experiments, Muszkiewicz et al. characterized variability in AP and ionic densities and their impact on CaT in atrial cells from right atrial appendage of patients exhibiting nSR (Muszkiewicz et al., 2017). In addition to calibrating model outputs to measured AP biomarkers, they also extended the experimental calibration of population of human atrial models to model parameter (inputs) by using experimental data of major ionic currents.

Quantifying drug modulatory effects, understanding variability in drug response, and identifying phenotype-specific therapy

By using a population of models that incorporate variabilities, drug modulatory effects on electrophysiological properties can be interpreted at a whole population level, which also contributes to limiting potential model-dependent results. For example, Yang et al. used a population-based approach to simulate effects of late Na+ current (INaL) and hERG block and found that the selective INaL blocker GS-458967 could suppress proarrhythmic markers after hERG block in ventricular myocytes (Yang et al., 2016). Population-based modeling has also allowed for more rigorous quantitative comparison of modulatory effects between multiple drugs. A recent study by Britton et al. calibrated populations of ventricular models to specific individuals using data from human trabeculae (Britton et al., 2017a). They then assessed the effects of four different (selective and non-selective) blockers of the rapid delayed-rectifier K+ current (IKr), dofetilide, sotalol, quinidine, and verapamil, to quantitatively compare changes in AP biomarkers, and demonstrated good agreement with experiments for the selective IKr blockers (dofetilide and sotalol) but not for the non-selective IKr inhibitors (quinidine and verapamil). Paci et al. utilized populations of in silico iPSC-CMs to evaluate antiarrhythmic effects of mexiletine and ranolazine to treat iPSC-CM long QT syndrome type 3 (LQT3) mutants and demonstrated that mexilitine stops spontaneous APs in more LQT3 models than ranolazine due to its stronger effects on the fast Na+ current (INa) (Paci et al., 2017). In contrast to the traditional modeling approach using a single model, the population-based modeling can gain insights into the physiologically relevant variability of predictions made in silico, as demonstrated in these studies.

By taking a step further, simulations using populations of models incorporating variabilities can also help recognize the contributing factors underlying the variability observed in response to drugs. One relevant example is the variable outcomes of hERG inhibition, which is frequently linked with diLQTS. Population-based modeling has offered insights into the mechanisms underlying the fact that individuals will not exhibit the same degree of QT interval prolongation after hERG block (Singh et al., 2000; Kannankeril et al., 2011; Weeke et al., 2014). Employing a population of models of ventricular myocytes, Sobie and Sarkar attributed the variable outcomes to the different ionic properties of the cells (Sarkar and Sobie, 2011). In another interesting application, Passini et al. implemented an in silico drug trial using experimentally-calibrated populations of AP models to investigate the risk of drug-induced arrhythmias, and to identify specific subpopulations at higher risk for proarrhythmic cardiotoxicity (Passini et al., 2017). Their methodology not only demonstrated higher accuracy than animal models in predicting arrhythmia risk, but also provided mechanistic insight into the underlying ionic contributors to repolarization/depolarization abnormalities.

Understanding the bases of variability in electrophysiological behavior and arrhythmia proclivity may also allow developing specific antiarrhythmic therapies for different disease phenotypes. For example, Liberos et al. compared AF models that had self-sustained vs. self-terminating reentries (Liberos et al., 2016). They found that AF maintenance was correlated with high ICaL and INa, and that ICaL block could be an effective treatment depending on the basal availability of Na+ and Ca2+ ion channel conductivities (INa depression increased efficacy). Mayourian et al. employed a comprehensive integrated approach to study the mechanisms of cardiac contractility and arrhythmogenicity using experimentally-calibrated human mesenchymal stem cells (hMSCs) (Mayourian et al., 2017). In simulations testing proarrhythmic effects, they found that hMSCs paracrine signaling protected such adverse effects of heterocellular coupling at various levels of engraftment. This work highlights that antiarrhythmic strategies can move beyond simply considering repolarization abnormalities.

Sample-specific modeling

Instead of taking a baseline cardiac model and introducing variability by randomly varying the ionic conductances, optimization and statistical techniques can also be used to tailor the baseline model to describe a specific experimental sample. Depending on the characteristics of the dataset at hand, sample-specific models can be representative of either individual myocytes or a particular group of cells. The former approach, cell-specific modeling, can be helpful when integrating mathematical modeling into an experimental setup. For example, Ravagli et al. characterized the role of the “funny” current If in sinoatrial myocytes using the dynamic clamp technique by adapting the extent of injected If in a cell-specific fashion, i.e., based on the basal firing rate measured in each individual cell (Ravagli et al., 2016). Despite the use of average data, sample-specific models built from a group of cells (e.g., a cell line developed in a certain laboratory, myocytes isolated from disease models, iPSC-CMs derived from a single patient) can allow for specific characterization of conditions that are far from the average, or even of personalized physiology (Barichello et al., 2018). For example, monophasic AP data recorded in AF patients undergoing ablation procedures have been used to construct atrial cell models that capture patient-specific electrophysiological properties (Lombardo et al., 2016). This approach has the promise of making patient-specific predictions given interventions such as arrhythmia-provoking protocols or drug application. Here we summarize methodologies for building and improving sample-specific cell models. For more detail, we refer the readers to a previous review on the topic (Krogh-Madsen et al., 2016).

Fitting sample-specific models

Sample-specific models can be constructed by fitting the parameters of a baseline model so that the model outputs match the corresponding physiological behaviors seen in a single patient or myocyte (Figure 2). Cardiac electrophysiology models can be optimized using many different fitness functions (Druckmann et al., 2007; Tomaiuolo et al., 2012), such as global search heuristics (Vanier and Bower, 1999; Dokos and Lovell, 2004; Bueno-Orovio et al., 2008; Guo et al., 2010). Recently, many sample-specific models are generated using the genetic algorithm (GA), which traces its beginnings to evolutionary biology (Fraser and Burnell, 1970; Crosby, 1982), but is still being applied in new ways today (Chen and Guan, 2004; Hussein and El-Ghazaly, 2004; Leung et al., 2004; Vieira et al., 2004). Its use for optimization of ionic models is relatively new in both neuroscience (Achard and De Schutter, 2006; Gurkiewicz and Korngreen, 2007; Hobbs and Hooper, 2008; Ben-Shalom et al., 2012) and cardiac electrophysiology (Syed et al., 2005; Bot et al., 2012; Kaur et al., 2014; Groenendaal et al., 2015). Syed et al. demonstrated its feasibility for atrial cell models when they proved they could fit two different cell models (Courtemanche et al., 1998; Nygren et al., 1998) to any given atrial AP (Syed et al., 2005). Essentially, the GA optimization procedure is initialized in the same way as for the population-based approach (varying maximal conductance and/or transport rates), and then it iteratively evolves toward better solutions in generations, while the underlying parameters can be varied, swapped, or discarded. Sensitivity analysis can be used in conjunction with generating sample-specific models as it can inform the design of the error function (i.e., weights) by revealing the conductances that more significantly impact the electrophysiological outputs used for fitting. For example, Krogh-Madsen et al. recently combined sensitivity analysis and global optimization (using a GA) of a ventricular myocyte model to clinical long QT data and intracellular Ca2+ and Na+ concentrations, to better constrain the model parameters (Krogh-Madsen et al., 2017). They found that this improved prediction of drug-induced torsades de pointes (TdP), especially in eliminating false-positive outcomes generated by the baseline model parameters.

Improving fidelity of sample-specific models

The final solution of an optimization procedure using some fitness function may not match experimental data well if only fitting to a single electrophysiological output such as a single AP, because multiple parameter combinations can potentially produce the same AP (Figure 3) (Syed et al., 2005; Druckmann et al., 2007; Sarkar and Sobie, 2010; Guo et al., 2013; Kaur et al., 2014; Groenendaal et al., 2015). In this case, although fitness function itself can be improved, for example, by increasing the population size or diversity for GAs can improve the fit of a sample-specific model, it may not necessarily guarantee that the final solution relates to the global minimum. To address this issue, many methods have been developed using (1) additional electrophysiological properties for fitting and/or (2) more complex electrophysiological protocols to improve model fidelity. It has been shown that model faithfulness can be improved by adding membrane resistance as an objective (Kaur et al., 2014), by optimizing to Ca2+ handling (Dokos and Lovell, 2004; Sarkar and Sobie, 2010; Rees et al., 2018), or by accounting for experimental data generated from multiple pacing frequencies (Syed et al., 2005; Lombardo et al., 2016) or irregular pacing protocols (Guo et al., 2013; Groenendaal et al., 2015).

In addition to using multiple electrophysiological properties to improve the fit of sample-specific models, more intricate voltage-clamp protocols that capture complex and rich electrophysiological dynamics have been employed, as first demonstrated to improve the fit of Markov models of ionic currents with many parameters (Dokos and Lovell, 2004; Zhou et al., 2009; Beattie et al., 2018). These can be implemented over a short time frame and may be used to emphasize certain currents over others. In the absence of pharmacological isolation, Groenendaal et al. used only a 6-s protocol that effectively isolated IK1, ICaL, and slow delayed-rectifier K+ current (IKs) given their disproportionately large contribution at voltage steps of −120, +20, +40, and −30 mV, respectively (Groenendaal et al., 2015). They found that using this protocol alone cannot isolate all ionic currents, and when used in combination with a stochastic pacing protocol there was a considerable improvement in parameter estimation. Developing short but information-rich protocols is useful especially when trying to improve the results of an optimization procedure for cell-specific modeling, because longer protocols take longer to implement experimentally and thus are difficult to perform in a single cell. In a recent study, Beattie et al. proposed an innovative experimental and mathematical modeling method that allows to concisely measure the dominant processes involved in hERG channel gating by applying a short (8-s long) “sum of sinusoids” voltage-clamp protocol (Beattie et al., 2018). The sinusoidal waves were able to provoke a wider range of non-equilibrium behavior than traditional square voltage steps, thus allowing rich and complete characterization of hERG channel kinetics in the same cell and efficient model fitting (Figure 4A).

Figure 4

Improving fit of sample-specific models. (A) Experimental and simulated IKr time courses (bottom) evoked in response to an efficient, information-rich sum-of-sinusoid voltage protocol (top) that allows rapid characterization of IKr behavior. (B) Steady-state peak IKr-voltage curves comparing cell-specific model predictions (bold, colored) to cell-specific experimental recordings (dashed, colored). The black lines in each plot are from the model calibrated to averaged sinusoidal data from all the cells (light gray). Reproduced from Beattie et al. (2018) with permission.

The final step in improving fidelity of sample-specific models is to directly experimentally test the predictions of the model given new protocols (Groenendaal et al., 2015; Devenyi et al., 2017; Beattie et al., 2018). Figure 4B reports an example of such validation experiments, where predictions obtained with cell-specific IKr models (identified applying the sinusoidal protocol in Figure 4A in 9 cells) are compared to the IKr-voltage relationships experimentally determined in each cell (Beattie et al., 2018). The order of the panels in Figure 4B is based on an index of recording stability (lowest to highest difference in leak resistance between the vehicle and dofetilide recordings) that is associated to “data quality”. Cell-specific predictions are excellent for cells 1–5 (higher data quality), but less accurate for cells 6–9 (lower data quality). The analysis also shows that cell-specific models provide better predictions than the average model for the cells with the highest data quality (cells 1–5). Experimental validation is an important last step in improving cell-specific models, as generating cell-specific models is potentially more susceptible to observational error. Devenyi et al. used a GA to re-parameterize the Livshitz-Rudy model of the guinea pig ventricular cardiomyocyte (Livshitz and Rudy, 2009), and predicted an increase in IKr and a drastic decrease in IKs given a dynamic clamp protocol as compared to the original model, and this was validated experimentally (Devenyi et al., 2017). Their adjusted model predicted that IKs can stabilize the AP and EADs better as compared to IKr—which improved the ability to assess arrhythmia risk, given the baseline model did not produce physiological behaviors that were quantitatively similar to their experiments.

Models of patient-specific anatomy

While a detailed discussion of patient-specific anatomical models is beyond the scope of our review, recent studies have begun investigating how inter-patient differences in myocardial structure affects atrial arrhythmia, as reviewed by Barichello et al. (2018). For example, Zhao et al. developed a 3D human heart-specific atrial computer model integrating 3D high resolution structural and functional mapping data to test the impact of wall thickness, fibrosis, and myofiber orientation on AF induction, maintenance, and ablation strategies (Zhao et al., 2017). Deng et al. demonstrated that reentrant driver localization dynamics are influenced by inter-patient variability in the spatial distribution of atrial fibrosis, APD, and conduction velocity (Deng et al., 2017). This suggests that incorporating patient-specific electrophysiological models in patient-specific geometries might enhance their predictive value. We discuss the computational challenges associated to this task in the section entitled “Arrhythmia Research Requires Understanding Variability at Larger Spatial Scales”. Furthermore, obtaining patient-specific electrophysiological data might constitute another logistical roadblock.

Overall, methods that incorporate variability are particularly useful for (1) analyzing variability in cardiac electrophysiology, (2) assessing proarrhythmic risk, (3) determining the underlying factors contributing to variable drug response, and (4) identifying phenotype-specific (and in the future patient-specific) antiarrhythmic targets. Table 1 summarizes applications of these approaches and new insights provided by the studies (shaded areas indicate atrial studies).

Table 1

Applications and main findings of computational methods incorporating cardiac electrophysiological variability (**shaded areas indicate atrial studies).

Author, Year	Baseline Model	Model Generation	Approach	Insight
ANALYZING VARIABILITY
Sobie, 2009	MVMMs: Luo and Rudy, 1991; Fox et al., 2002; Kurata et al., 2005	Sampled from LND	MLR	New method to rapidly identify ionic mechanisms shaping AP properties, CaT, and alternans
Sanchez et al., 2014	HAMMs: Courtemanche et al., 1998; Maleckar et al., 2009; Grandi et al., 2011	Sampled over a ±100% variation range around their baseline values as described by Marino et al. (2008); calibrated to AP recordings in atrial trabecula	Statistical difference tests	Ionic determinants of variability in human AP in nSR vs. cAF
Lee Y. S. et al., 2016	HAMM: Courtemanche et al., 1998	Sampled from LND	MLR	Comparison of parameter sensitivity between nSR and AF condition. Ionic contributions to rate-dependence and spiral wave dynamics in AF
Devenyi and Sobie, 2016	HAMM: Grandi et al., 2011; HVMM: Grandi et al., 2010; RVMM: Pandit et al., 2001	Sampled from LND	MLR	In human atrial myocytes, both Ito and IKur had greater impacts on CaT amplitude than did SERCA. This was similar in rat left ventricular epicardial cells, where Ito played a more important role than SERCA
Vagos et al., 2017	HAMM: Skibsbye et al., 2016	Sampled from Gaussian distribution	MLR	Ionic determinants of unstable behaviors in nSR vs. cAF
Ellinwood et al., 2017b	HAMM: Grandi et al., 2011	Sampled from LND	MLR	IKur impacts APD and effective refractory period more in cAF (even though it is downregulated) vs. nSR
Muszkiewicz et al., 2017	HAMMs: Courtemanche et al., 1998; Maleckar et al., 2009; Grandi et al., 2011	Sampled using LHS and sequential MC; calibrated to experimental recordings	PCCs, statistical difference tests	Ionic determinants of electrophysiological and CaT properties
Lawson et al., 2018	HAMM: Courtemanche et al., 1998	Sampled over a ±100% variation range around their baseline values; sequential MC; model calibrated based on distributions of biomarkers estimated from multivariate kernel density estimation	Statistical difference tests	Accurate identification of inherent variability within the experimental population and improved characterization of ionic differences between nSR and cAF
ASSESSING ARRHYTHMIA RISK
Walmsley et al., 2013	HVMM: O'Hara et al., 2011	MC sampling from a uniform distribution (±30%); calibrated to mRNA expression data in failing and non-failing hearts	MLR	Combination of low SERCA activity and high ICaL conductance impacted the formation of alternans the most in the non-failing heart population, but low hERG conductance was the main contributor to alternans in the failing heart population
Zhou et al., 2016	HVMM: O'Hara et al., 2011	Sampled using LHS; calibrated to in vivo recordings	PCCs, statistical difference tests	ICaL and NCX current determine the cell-to-cell differences in repolarization alternans through intracellular and sarcoplasmic Ca2+ regulation
Britton et al., 2017b	HVMM: O'Hara et al., 2011	Sampled using LHS; calibrated to data in human ventricular trabeculae	Logistic regression, PCCs, statistical difference tests	Na+/K+ pump is a key determinant of repolarization abnormality susceptibility
Morotti and Grandi, 2017	HAMM: Grandi et al., 2011	Sampled from LND	Logistic regression	EADs are particularly sensitive to conductances of INa, acetylcholine-sensitive and ultra-rapid K+ channels, and NCX transport rate
Devenyi et al., 2017	Guinea pig left ventricular myocyte model Livshitz and Rudy, 2009	Fit using GA; sampled from LND	Dynamic clamp data for fitting, MLR	IKs is more capable to stabilize AP and EADs as compared to IKr
IDENTIFYING ANTIARRHYTHMIC TARGETS
Cummins et al., 2014	MVMMs: Luo and Rudy, 1991; Fox et al., 2002; Hund and Rudy, 2004; Ten Tusscher et al., 2004; Ten Tusscher and Panfilov, 2006; Livshitz and Rudy, 2009; O'Hara et al., 2011	Sampled from LND	MLR	IK1 and Na+/K+ pump currents favor forward rate dependence
Liberos et al., 2016	HAMM: Skibsbye et al., 2016	Sampled using LHS; calibrated to AP recordings in atrial trabeculae in patients with AF	PCCs, statistical difference tests	AF maintenance was correlated to high ICaL and INa, and ICaL block could be an effective treatment depending on the basal availability of Na+ and Ca2+ channel conductivities
Passini et al., 2016	HVMM: O'Hara et al., 2011	Sampled using LHS; calibrated to non-diseased and HCM myocytes AP recordings	Analysis of repolarization properties	ICaL re-activation is the key mechanism for repolarization abnormalities in HCM myocytes, and combined NCX, INaL and ICaL block is effective to partially reverse the HCM phenotype
Yang et al., 2016	Rabbit ventricular myocyte model Soltis and Saucerman, 2010	Randomly selected within ±10% of nominal value (uniform distribution)	Analysis of TRIaD pro-arrhythmic markers	GS-458967 suppressed proarrhythmic markers following hERG block
ASSESSING THE VARIABLE RESPONSE TO DRUGS
Sarkar and Sobie, 2011	MVMMs: Fox et al., 2002; Hund and Rudy, 2004; Ten Tusscher et al., 2004; Kurata et al., 2005; Grandi et al., 2010	Sampled from LND	MLR	Individuals do not exhibit the same degree of QT interval prolongation due to different ionic ensembles
Britton et al., 2013	Adapted rabbit Purkinje cell model Corrias et al., 2011	Sampled using LHS; calibrated to experimental data	PCCs	Quantitatively predicted the arrhythmia risk of four concentrations of the K+ channel blocker dofetilide; baseline IKr conductance is the primary determinant of APD prolongation caused by dofetilide
Lancaster and Sobie, 2016	HVMMs: Ten Tusscher et al., 2004; Grandi et al., 2010; O'Hara et al., 2011	Sampled from LND; calibrated to experimental data	PCA, ROC curves, TdP risk scores	TdP risk assessment could be improved by quantifying the impact of multiple cardiac ion channels (even those not typically considered to affect risk)
Britton et al., 2017a	HVMM: Adapted O'Hara et al., 2011	Sampled using LHS, calibrated to heart-specific ex vivo measurements	Coefficients of variation	Good agreement with experiments for selective IKr blockers, but notable differences for the non-selective IKr inhibitors
Passini et al., 2017	HVMM: O'Hara et al., 2011	Sampled using LHS; calibrated to experimental data	TdP scoring system	In silico drug trial demonstrated higher accuracy than animal models in predicting arrhythmia risk (89%); underlying ionic contributions to repolarization/depolarization abnormalities
Krogh-Madsen et al., 2017	HVMM: O'Hara et al., 2011	Fit using GA; optimized to clinical data	TdP risk prediction	TdP risk assessment could be improved by using global optimization methods and multi-variable objectives
Gong and Sobie, 2018	HVMM O'Hara et al., 2011 and human iPSC-CM models Paci et al., 2017	Sampled from LND	Cross-cell MLR	Cross-cell regression predicted adult ventricular myocyte drug responses from the behaviors of an iPSC-CM in silico population

HAMM, human atrial myocyte model; HCM, hypertrophic cardiomyopathy; HVMM, human ventricular myocyte model; LHS, Latin hypercube sampling; LND, log-normal distribution; MLR, multivariable linear regression; MVMM, mammalian ventricular myocyte model; RVMM, rat ventricular myocyte model; PCA, principal component analysis; PCC, partial correlation coefficient; ROC, receiver operator characteristic; MC, Monte Carlo; TRIaD, Triangulation, Reverse use dependence, beat-to-beat Instability of action potential duration, and temporal and spatial action potential duration Dispersion.

Challenges and future directions

We reviewed the most common methods used to account for variability in cardiac electrophysiology, which largely fall into the two categories of (1) population-based modeling and (2) sample-specific modeling. These methods complement each other well, as population-based methods can help characterize behavior in a particular patient group (healthy, diseased, stressed, etc.), and sample-specific modeling shows significant promise to develop personalized medical approaches for individual patients. Both methods have led to many important insights into the mechanisms of arrhythmogenesis and antiarrhythmic strategies. However, there are several important limitations to consider, which suggest potential future developments in modeling of cardiac electrophysiology.

Analysis of electrophysiology from populations of models may require different statistical methods

As opposed to the traditional approach of producing a single value from a single baseline model, models that incorporate variability have allowed statistical methods to be applied that can either ask new scientific questions or quantify the impact of variability on electrophysiological outputs, as performed in experimental studies. While the statistical analysis methods used in experimental studies can be directly applied in the in silico population-based studies, differences in the nature of experimental and simulation studies may need to be considered. For example, some population-based techniques generate model population sizes (often sample sizes in the 1000s) that are much greater than could be achieved by experiments (often sample sizes of 3–12) or the traditional cardiac modeling approach alone. Therefore, given similar effects, results produced in the population-based simulations have greater statistical power to detect differences. Furthermore, because even very small effects can reach statistical significance with large samples, physiological significance should be assessed (White et al., 2014). Additionally, when evaluating drug effects on electrophysiology, in simulations of the same virtual cell (a single model out of the population-models) can be used to perform both control and with-drug studies, allowing for paired comparisons, which is often not practical in experimental studies. The methodologies for analyzing and interpreting the “big data” produced by the population of models should be carefully considered and standards should be established going forward.

Variability does not fully account for uncertainty

Physiologic variability should be thought of in the context of the broader umbrella of uncertainty, which is the confidence (or precision) with which a quantity, such as an electrophysiological output, can be given a value (Mirams et al., 2016). While here we reviewed how cardiac electrophysiology models have begun to account for physiological and experimental variability, uncertainty analysis should determine whether the baseline model itself is a valid representation of its physical system. The extent and rigor of validation during model development affects uncertainty, whereby the broader the set of constraints, e.g., the model recapitulates both voltage and Ca2+ responses, their pacing rate-dependence, short- and long-term behavior, the lesser the uncertainty in the obtained parameters. Uncertainty analysis should also verify that the experiments used to construct the model are appropriate. For example, in experiments, voltage-clamp protocols used to characterize ionic currents are often done using non-selective pharmacological block which may have unidentified effects, over a range narrower than the entire physiological range, or on larger cells that are easier to patch-clamp with intrinsically greater than normal maximal conductances (Courtemanche et al., 1998). All of these would lead to uncertainties in the initial parameters and conditions due to experimental error and lack of knowledge. Likewise, the choice of the computational methods or resources used to perform the model parameterization and simulations can produce uncertainty in model results. This is because cardiac models may use different mathematical equations to describe the same physiological process, perhaps based on different analyses or assumptions on the physical-world process. Using more than one (e.g., cell) model to confirm predictions or validate the mechanistic understanding of a process is therefore a useful strategy (see for example Sarkar and Sobie, 2011; Sanchez et al., 2014; Lancaster and Sobie, 2016; Muszkiewicz et al., 2017). Additionally, even the choice of the numerical solver used by the software can lead to variability in model outputs, i.e., simulator uncertainty (Pathmanathan et al., 2012). Moreover, uncertainty in model outputs may arise if the code has not been verified to truly represent the mathematical equations in the computational model (Niederer et al., 2011a). Finally, optimization procedures can also introduce uncertainty, whereby the choice of whether to optimize simplified models with few parameters (Bueno-Orovio et al., 2008; Al Abed et al., 2013; Guo et al., 2013) or detailed models but only a few properties (e.g., focusing on specific currents) (Zhou et al., 2009; Fink et al., 2011) can lead to multiple distinct models given the same experimental data.

Used in conjunction with the approaches discussed in this review that take into account electrophysiological variability, new methods have been developed that try to quantify uncertainty more generally (Marino et al., 2008). Uncertainty quantification methods aim to quantify uncertainties in model inputs and propagation through the model to see how they affect model predictions (Pathmanathan et al., 2015; Mirams et al., 2016). This is typically done by assigning probability distribution functions, rather than fixed values to model parameters, as done for example by Pathmanathan et al. and applied to the study of INa steady-state inactivation (Pathmanathan et al., 2015). However, this process can be slow and tedious (requiring lots of simulations), especially if using a Monte Carlo sampling method that chooses input values from a statistical distribution. Also, in some cases, this statistical distribution of input parameters can be difficult to obtain or justify experimentally. To solve this issue, uncertainty quantification analysis has developed surrogate models or emulators (e.g., polynomial chaos expansions, and Gaussian process emulators; Chang et al., 2015), which are fast-running statistical approximations of the computational models and are quite powerful when fit to carefully constructed training data. Formal studies using uncertainty quantification in cardiac models are still limited, given the huge number of parameters in cardiac models, and may require the development of new methods or computational techniques (Johnstone et al., 2016).

Potential covariance in ionic conductances challenges the current method of incorporating variability

Currently, populations of cardiac models and sample-specific models typically calibrate or fit to maximal conductance values or transport rates of channels or receptors, based on the observation that changes in expression levels of ion channels and transport proteins are the primary contributors to (inter-species) variability (Rosati et al., 2008). However, this approach does not take into account that the expression of ion channels will vary over relatively short time scales given changes in transcription, translation, degradation, or even circadian rhythm. Moreover, with a few exceptions (Sarkar and Sobie, 2011; Cummins et al., 2014) these methods do not typically account for variability in ion channel kinetics, which is known to change especially in response to drugs (Clancy et al., 2007). The methods discussed in this review can attempt to account for these properties using additional parameters.

Although the correlation between parameters (i.e., maximal conductances) is assessed sometimes (Britton et al., 2013), neither population-based nor the sample-specific approaches account for possible covariance in ion channel conductances, despite the fact it has been observed in neurons (Schulz et al., 2006, 2007; Tobin et al., 2009) and cardiac tissue (Schram et al., 2002; Rosati and Mckinnon, 2004; Deschenes et al., 2008; Xiao et al., 2008; Banyasz et al., 2011; Milstein et al., 2012). The exact mechanisms responsible for these covariances are still being explored. Xiao et al. found that sustained reductions in IKr may lead to compensatory upregulation of IKs through post-transcriptional upregulation of underlying subunits (Xiao et al., 2008), which potentially underlie the observed phenomenon of repolarization reserve (Roden, 2008). Macromolecular complexes or post-transcriptional modifications could also facilitate coregulation of ionic conductances, as demonstrated by the observed structural or functional complex between Ito and INa (Deschenes et al., 2008). Rees et al. recently argued that sensing of aggregate CaT may be sufficient in itself to regulate ionic conductances (of K+ and inward Ca2+) to maintain normal Ca2+ handling (Rees et al., 2018). Moreover, knockout and knockdown studies are consistent with the idea that cardiac cells have compensatory mechanisms to maintain AP or CaTs (perhaps to prevent arrhythmias) (Guo et al., 1999; Zhou et al., 2003). The covariance of ionic conductance can have significant implications for both calibrating populations of models or fitting sample-specific models, because it could propose additional constraints for how the underlying parameters of the computational model can be varied. Thus, new methods have begun to be developed that measure all ionic conductances at once, and can not only tease out how ionic conductances are correlated, but the extent to which they vary between cells (Banyasz et al., 2011; Groenendaal et al., 2015).

Arrhythmia research requires understanding variability at larger spatial scales

Accounting for variability at tissue and organ-level scales is a logical, but not trivial (Elshrif and Cherry, 2014), next step. A thorough investigation of variability would first require including differences among the cells in the same tissue, and evaluating the impact of diverse geometrical distributions. One should also account for patient-specific structural differences, based on measures of tissue conductivity and anatomic properties, including heterogeneity in signaling due to non-uniform innervation. This last step can be particularly problematic when investigating diseased conditions affected by pronounced structural remodeling, such as fibrosis, organ dilation, and alterations in gap junction coupling. Where there has been meaningful progress in accounting for variability is in developing personalized atrial model structures based on medical images (Dossel et al., 2012; Trayanova, 2014). These methods have shown some promise in developing personalized ablation strategies (Mcdowell et al., 2015). For example, recently, Soor et al. implemented a modeling approach to optimize ablation times based on patient-specific atrial geometries to create lesions for a given atrial wall thickness (Soor et al., 2016). Combining these methods that utilize medical images with the methods described here, could significantly improve the clinical value of both methods alone (Hansen et al., 2017; Zhao et al., 2017). For example, Arevalo et al. developed personalized heart models based on cardiac imaging and published patch-clamp data to better predict arrhythmic events and possibly avoid unnecessary implantable cardioverter defibrillators (Arevalo et al., 2016). Developing multi-scale frameworks that account for variability is the next frontier in cardiac modeling that will greatly benefit from further advancements in computing capabilities. Beyond the availability of large experimental and clinical datasets, the development of novel techniques to speed model derivations and to integrate automation will be crucial to capture variability for different cell types and conditions at various scales.

Safety pharmacology requires complementary electrophysiological experimental methods

The in silico approaches described here are being combined with other state-of-the-art tools to improve the evaluation of drug safety. Of significance, these approaches can help further the mission of the CiPA (Comprehensive in vitro Proarrhythmia Assay) initiative, which aspires to develop better methods to predict TdP. Beyond exclusively using steady-state hERG block as the main predictor of arrhythmia and not at all using QT interval prolongation, the CiPA initiative attempts to gain a more comprehensive understanding of proarrhythmic risk by combining (1) mechanistically-based in vitro assays, (2) in silico reconstructions of cardiac electrophysiology, and (3) confirmation using human iPSC-CMs (Colatsky et al., 2016). The methods described in this review are being utilized to help meet the mission of the CiPA initiative (Cummins et al., 2014; Lancaster and Sobie, 2016; Britton et al., 2017a; Passini et al., 2017). Most of the methods described here that assess the effects of drugs on populations of cardiac myocytes use simple pore block schemes. However, it is also clear that the sole use of steady-state hERG block assays is insufficient to predict arrhythmia risk, and thus studies are beginning to simulate the effects drug-binding kinetics and state-specific binding, which have been shown to affect electrophysiological outcomes (Lee W. et al., 2016; Dutta et al., 2017; Ellinwood et al., 2017b; Li et al., 2017). Incorporating more detailed drug-binding effects may allow studying the effects of populations of drugs characteristics (e.g., state-dependent block and kinetics) on populations of cardiomyocytes.

Conclusion

Computational approaches that have been developed over the past decade to account for variability in cardiac electrophysiology have led to important insights into mechanisms of arrhythmogenesis, etiology of disease, and variable response to drugs. The approaches outlined in this review are used in basic research studies, i.e., quite separately from actual clinical workflows, where decisions are made sometimes for a particular patient within minutes. Advanced computing facilities now allow near real-time simulations of anatomically realistic, biophysically detailed models of human cardiac electrophysiology (Niederer et al., 2011b; Okada et al., 2015, 2017). Such massively parallel processes could be optimized to run personalized cardiac simulations pre-determined to have clinical value. However, implementing these approaches more comprehensively into clinical workflows still presents challenges and simulation of variability may not find immediate application.

Author contributions

HN, SM, and EG reviewed the literature and wrote the manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.