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

Abstract The L‐type calcium current (ICaL) plays a critical role in cardiac electrophysiology, and models of ICaL are vital tools to predict arrhythmogenicity of drugs and mutations. Five decades of measuring and modeling ICaL 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 ICaL 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 ICaL. This article is categorized under: Cardiovascular Diseases > Computational Models Cardiovascular Diseases > Molecular and Cellular Physiology


| INTRODUCTION
The "long-lasting" or L-type calcium current (I CaL ) is an ionic current found in cardiomyocytes, neurons, endocrine cells, and other cell types throughout the body (Striessnig et al., 2014). It plays a crucial part in critical functions such as hormone secretion, regulation of gene expression, and contraction of cardiac and smooth muscle (Hofmann et al., 2014). Its critical role in cardiac cellular electrophysiology has led to extensive modeling efforts and a large number of models of cardiac I CaL electrophysiology have been published since the 1970s. Incorporated into models of the cardiac action potential (AP), I CaL models have a long and successful history of application in fundamental electrophysiology research (Noble & Rudy, 2001). More recently, they have been used or proposed for use in drug safety assessment (Mirams et al., 2012) and risk stratification in cohorts with ion channel mutations (Hoefen et al., 2012). Choosing an I CaL model for such safety-critical applications is not an easy task, as several models can be found in the literature, which vary widely in their structure and assumptions, and are often published in an article form which does not lend itself to quantitative comparison without considerable effort (Cooper et al., 2011). On a more fundamental level, each model of I CaL represents a testable theory of its physiology, and the existence of so many competing theories reveals both gaps in our knowledge and a need for their comparison and synthesis. In this review, we take the first major steps to facilitate a critical reassessment by the community, by collecting I CaL models, analyzing their qualitative differences, and providing quantitative comparisons based on simulations with freely reusable code. We start with a brief description of I CaL physiology and the common structure of its models.

| Biophysical properties of I CaL
The L-type calcium channels (LCCs) through which I CaL flows consist of a pore-forming α 1 subunit and the auxiliary subunits β and α 2 δ (Dolphin, 2016). The most common α 1 subunit in Purkinje cells and ventricular and atrial cardiomyocytes is Ca V 1:2 (encoded by the gene CACNA1C), while Ca V 1:3 (CACNA1D) is more prevalent in the sinoatrial node (SAN) and atrio-ventricular node (AVN) (Gaborit et al., 2007;Zamponi et al., 2015).
Electrically, I CaL is characterized by fast, membrane voltage V m ð Þ-dependent activation at relatively depolarized potentials (around 0 mV) and slower inactivation that depends on V m (voltage-dependent inactivation, VDI) and, via a separate process (Hadley & Lederer, 1991), on the intracellular calcium concentration (calcium-dependent inactivation, CDI). Although LCCs are highly selective and I CaL is predominantly carried by calcium ions, smaller potassium and sodium components have also been measured (Hess et al., 1986) and are commonly included in models.
An influential study by  on currents through isolated LCCs introduced the idea that they have three different modes of gating: A "mode 0" in which channel openings are rare, a "mode 1" characterized by rapid bursts of brief openings, and finally a "mode 2" (rarely observed without first introducing a channel agonist) that features long channel openings and only brief periods of channel closure. These long openings, together with their large unitary conductance lead to the "L-type" designation (Nilius et al., 1985). Although in this review we focus on wholecell "aggregate" I CaL , many of the models we will discuss are inspired by concepts from single-channel work.
In ventricular and atrial cardiomyocytes, I CaL is the major inward current responsible for maintaining the "plateau" of the AP (Catterall, 2011), making it a major determinant of the AP duration ( Figure 1) and important cell properties such as restitution and refractoriness. The calcium influx during the plateau phase initiates the process of calciuminduced-calcium release (CICR): calcium flowing in through LCCs causes the nearby ryanodine receptor (RyR) channels to open, leading to a further, and much larger, influx of calcium from the sarcoplasmic reticulum (SR, see Figure 3), the major calcium store inside the cell. The resulting increase of calcium concentration in the cytosol allows formation of cross-bridges between myosin heads and actin filaments and leads to contraction (Eisner et al., 2017;Winslow et al., 2016).
The crucial role of calcium channels in maintaining healthy cardiac function is borne out by clinical evidence. For example, reduction of the expression of the Ca V 1:2 α 1 subunit by less than half can lead to heart failure (Goonasekera et al., 2012), and mutations in CACNA1C and in its subunits have been linked to Brugada syndrome, Timothy syndrome, arrhythmia, and structural heart defects (Hofmann et al., 2014). This same vital role makes LCCs an important target for pharmacological antiarrhythmic therapy. For example, nifedipine, verapamil, and diltiazem inhibit I CaL (Ortner & Striessnig, 2016). Unintended pharmacological modulation of I CaL , conversely, can lead to drug-induced arrhythmias and other adverse cardiovascular effects. Initiatives like the comprehensive in vitro proarrhythmia assay (CiPA) therefore take I CaL into account when attempting to predict pro-arrhythmic risk (Li et al., 2019).
Several of the body's regulatory mechanisms target I CaL . Following sympathetic stimulation ("fight-or-flight" response), β1-adrenergic membrane receptors are activated, which triggers intracellular processes resulting in the conversion of ATP to cyclic adenosine monophosphate (cAMP), which activates protein kinase A (PKA), which phosphorylates the LCCs. This increases I CaL leading to a larger Ca 2þ influx and a larger Ca 2þ release through CICR, which ultimately increases contractile strength (Gao et al., 1997;van der Heyden et al., 2005). Parasympathetic stimulation ("rest-and-digest") can reduce cAMP levels, which is known to inhibit the effects of sympathetic stimulation on I CaL but appears to have little effect on the basal current (McDonald et al., 1994). In addition to these mechanisms, Ca 2þ /calmodulin-dependent protein kinase II (CaMKII) is both regulated by Ca 2þ Â Ã , and a regulator of Ca 2þ Â Ã via phosphorylation of both LCCs and RyRs (Maier & Bers, 2007;Winslow et al., 2016). This co-dependence between Ca 2þ Â Ã and CaMKII activity can create a positive feedback mechanism in which peak I CaL increases when successive voltage pulses are applied, in a process called Ca 2þ -dependent facilitation (CDF, Bers & Morotti, 2014) and sometimes equated with "mode 2" channel gating. Unusually high CaMKII levels have been observed in conditions such as heart failure and are thought to be arrhythmogenic (Bers & Morotti, 2014;Soltis & Saucerman, 2010). Finally, LCCs are not spread uniformly throughout the cell membrane, but occur mostly (90%) in specialized regions called dyads (Scriven et al., 2000). In dyads, the outer (t-tubular) cell membrane is in close proximity to the SR membrane (<15 nm, see Fawcett & McNutt, 1969;Eisner et al., 2017) and LCCs are in very close contact with RyRs. This enables CICR (Stern, 1992) and allows fast signaling (Abriel et al., 2015;Harvey & Hell, 2013), but the small size of these "nanodomains" also means that Ca 2þ Â Ã can fluctuate much faster and with a much higher amplitude than in the rest of the cytosol. Models of I CaL are often not created in isolation, but as part of a larger AP model that makes certain assumptions about LCC localization and hence about the Ca 2þ Â Ã affecting I CaL . As a result, the assumed channel localization must be taken into account when comparing models of I CaL .

| Models of cardiac I CaL
The first published models of a cardiac calcium current appeared in the 1970s (Bassingthwaighte & Reuter, 1972;McAllister et al., 1975), only a few years after the first calcium current recordings in cardiac preparations. Although Bassingthwaighte and Reuter (1972) confidently referred to a "calcium current whose role is primarily in the maintenance of the plateau [and in] regenerative depolarization in low-sodium medium", the exact nature of the current and the species it was carried by was controversial at the time, so that the more cautious names "secondary inward" and "slow inward current" were commonly used (Fozzard, 2002). Although the complexity of I CaL models has increased over time (keeping pace with new discoveries about the current and its role in cardiac function), all models that we assess in this manuscript can be expressed in a common form: Here, O is the fraction of open channels (or the open probability of a single channel). The time-varying behavior, that is, the kinetics (or gating), of the channel(s) is captured by O, which is voltage-and time-dependent and may also depend on Ca 2þ Â Ã . The term δ represents the "driving term" accounting for the electrochemical forces driving the ionic movement through the open channels. δ can depend on both V m and the ionic concentrations on either side of the channel. The current is scaled by either a maximal conductance g or a maximum permeability P, with different units depending on the choice of driving force model (see Section 2.5). Conductance and permeability are functions of the single channel permeability and the total number of channels in the membrane (Schulz et al., 2006), which can vary over longer periods (hours to days); however, they are treated as constants in all of the models we surveyed.
We shall discuss a variety of approaches taken to modeling O and δ, as well as factors not included in Equation (1) such as channel localization and regulation by other variables. In general, we shall aim to take a reductionist approach and regard components affecting I CaL (e.g., local calcium dynamics, CaMKII signaling or β-adrenergic stimulation) as separate entities. Examples of the full equations for early (McAllister et al., 1975) and more recent models (Faber & Rudy, 2000) can be viewed in the online repository accompanying this article, at https://github.com/CardiacModelling/ ical-model-review.

| Outline
In this review we identify, discuss, and compare the models that have been used to describe cardiac mammalian I CaL to date. In the next section, we describe the process by which we have identified models to be included in the review and classify the 73 models we identified by their origins, the assumptions made by their authors about LCC localization and the variables that regulate I CaL , and the equations used for I CaL gating and driving force. In Section 3, we compare a subset of 60 models quantitatively, by simulating the application of several experimental protocols. These include voltage-clamp protocols for activation, inactivation, and recovery, as well as an AP-clamp protocol and a protocol in which AP and calcium transient (CaT) are both clamped. We then test whether we can reduce the observed variability in predictions by grouping models according to the qualitative aspects identified earlier. In the final section we discuss possible approaches for choosing an I CaL model for a simulation study and developing and documenting new models of I CaL .

| QUALITATIVE COMPARISON
To identify models of I CaL we searched on PubMed (in March 2021) for publications containing the term "L-type calcium channel". Because many I CaL models are presented as part of a larger modeling effort, we also included the term "cardiac cell model" and scanned lists of models such as those given in Noble et al. (2012) and Heijman (2012). To constrain our study, we only included I CaL models that were presented by the authors as representing I CaL from healthy human or mammalian cells. Similarly, we focused on basal (vs. adrenergically stimulated) I CaL models.
Older models of "calcium current" and "slow" or "second inward" current were included if they modeled atrial, ventricular, or Purkinje cells, in which L-type is the dominant calcium current (Gaborit et al., 2007). The models included this way were by Bassingthwaighte andReuter (1972), McAllister et al. (1975), Beeler and Reuter (1977), DiFrancesco and Noble (1985), Hilgemann and Noble (1987), and Noble et al. (1991). Some older models were excluded if they modeled the SAN, in which the dominant calcium current is of the T-type (Mesirca et al., 2014). Next, we selected modern models where the authors explicitly included an L-type calcium current.
Some studies included several model variants. The models by Pandit et al. (2001), Ten Tusscher et al. (2004), Ten Tusscher and Panfilov (2006, O'Hara et al. (2011), andTomek et al. (2019) have epi-, endo-, or mid-myocardial variants but the I CaL kinetics are the same. In simulations where a whole-cell model was used (Section 3.6), we used the epicardial variant of these models. Similarly, the apical cell model (representing cells at the apex/tip of the ventricles) was chosen for Bondarenko et al. (2004). The I CaL kinetics vary among the different versions presented for Inada et al. (2009), Paci et al. (2013, and Varela et al. (2016), for which we chose the atrio-nodal (a transitional region between atrium and atrioventricular node), atrial-like, and right-atrial model, respectively.
The next selection was based on the model equations. Models were excluded from our overview if they only changed the value of the conductance or permeability parameter from a previous model; models were included if they introduced new equations, made changes to equations, or changed rate or driving force parameter values. Models for which the equations were not fully given in the publication or an online addendum were also excluded in this step (9 models in total), as were models for which the equations were suspected to contain typesetting errors (1 model). Due to its historic importance, an exception was made for Bassingthwaighte and Reuter (1972) for which the published equations did not match the accompanying figures. This model is included in our qualitative, but not quantitative analysis.
Finally, studies that calculated the current by considering the sum of all individual LCCs stochastically were only included for qualitative analysis. These included the studies by Greenstein and Winslow (2002), Restrepo et al. (2008), Hashambhoy et al. (2009), andNivala et al. (2012). The studies by Hinch et al. (2004), Greenstein et al. (2006), Asakura et al. (2014), andHimeno et al. (2015) described I CaL using models in which the gating was linked to the ryanodine receptors, and were also included for qualitative analysis only.
Applying these criteria we found a total of 73 distinct I CaL models, spanning five decades of research.

| Model provenance
Models are not created in isolation but inherit ideas, equations, and parameter values from previous models. Figure 2 presents a tentative "phylogeny" of I CaL models, showing how the models that are in use today were derived from one another over time. To establish these relationships, we searched the model publications for explicit references to parent models. In some cases (e.g. Rasmusson et al., 1990) the parent models did not meet the selection criteria used in this study, which is indicated in the figure by a dashed line. In other cases (e.g. Bondarenko et al., 2004;Demir et al., 1999;Lindblad et al., 1996) the parent model was not mentioned explicitly but could be inferred from the equations. In constructing the figure, we attempted to limit each model to a single parent. This is a simplification and the figure should not be read as a statement on originality. For the models up to around 1995 it is likely that each model "inherited" in some way from all other models published up to that time.
The phylogeny shows that many models have encountered species switches over the course of their development. For example, Grandi et al. (2010) (human) inherits from Shannon 2004 (rabbit), Luo-Rudy 1994 (guinea pig), and Rasmusson 1990 (bullfrog). Similarly, the histories of the O' Hara et al. (2011) and Paci et al. (2013) models also contain several species switches. Note that this does not necessarily suggest re-use of data, as models are frequently reparameterized to new data sets. However, inheritance of parameter values (and thereby indirectly of data sources) is very common in electrophysiology models (Niederer et al., 2009).

| Local Ca 2þ
Â Ã and spatial organization Several studies, both experimental and mathematical, have argued that I CaL , CICR, and I CaL -regulation cannot be understood without treating the Ca 2þ concentration at the intracellular channel mouth as a separate entity from the "global" or "bulk" cytosolic concentration. Perhaps the earliest suggestion is found in Bassingthwaighte and Reuter (1972), who argued for the existence of a local subspace with elevated Ca 2þ Â Ã to explain the observed reversal potential for I CaL , which was much lower than predicted by theory. An influential analysis of CICR by Stern (1992) reviewed the experimental evidence that the amount of Ca 2þ released by the SR is proportional to the "trigger" Ca 2þ Â Ã with a very high amplification factor, and went on to show that such a graded, high-amplification response could not be achieved in models that considered only a single global Ca 2þ Â Ã : A series of experimental studies investigating the mechanism behind (and calmodulin dependence of) CDI suggested that both a "local" and a "global" Ca 2þ Â Ã were involved (Peterson et al., 1999;Tadross et al., 2008) (in this theory, the "local" level refers to very brief spikes in Ca 2þ Â Ã that occur at the channel mouth when an LCC is open, while the "global" level is the lower Ca 2þ Â Ã between spikes, which is dominated by diffusion from the cytosol). Finally, some models of I CaL regulation have argued for highly localized presence of signaling molecules near LCCs, in concentrations very different from the bulk cytosol (Abriel et al., 2015;Harvey & Hell, 2013).
Modelers of the cardiac AP have accounted for these local phenomena by dividing the cell into smaller-volume "subspaces" in which calcium ions, carried in by I CaL or released from the SR, can cause Ca 2þ Â Ã to transiently rise to much higher levels than those seen in the bulk cytosolic space (Colman et al., 2022). flows directly into the bulk cytosol (e.g. DiFrancesco and Noble, 1985)). The second type of cardiac model adds a "dyadic" or "junctional" space, sometimes called the "cleft," between LCCs, located in the t-tubules, and the RyR in the adjacent SR. The illustration also shows LCCs connected directly to the bulk cytosol, as some cardiac models have a fraction ($10%) of I CaL flowing directly into the cytosol (e.g. Tomek et al., 2019). The third type of subspace model retains the dyadic space, but adds a "submembrane" or "subsarcolemmal" space just below the remainder of the cell membrane, so that I CaL flows either into the dyadic or the submembrane space. An example of this configuration, including compartment sizes and experimental references, is given in Shannon et al. (2004). Note that there are no hard boundaries between the subspaces: calcium can diffuse freely between them, but it is assumed that the rate of diffusion is slow enough for subspace Ca 2þ Â Ã elevations to arise. In this review, we use Ca 2þ Â Ã i to denote the bulk cytosolic concentration, Ca 2þ Â Ã s for the submembrane space concentration, and Ca 2þ Â Ã d for the dyadic space concentration. In Figure 3, and in the models it represents, a single t-tubule is used to represent a vast t-tubular network, and a single dyadic space represents all (thousands of) dyads. To allow further spatial heterogeneity, some studies (particularly those interested in calcium "sparks" arising locally within a myocyte) have gone further and discretized the cell into a 3D network of "functional release units", each with their own calcium concentration (see e.g., Nivala et al., 2012;Rice et al., 1999). Alternative, less computationally intense, approaches have been used by for example, Nordin (1993), who define a "superficial, middle, and deep myoplasm", and Koivumäki et al. (2011), who modeled the myocyte as a series of concentric shells.
These different models of cell geometry and Ca 2þ Â Ã concentrations complicate the task of comparing I CaL models in an AP model. Different assumptions about channel localization in an AP model will lead to LCCs experiencing different intracellular calcium concentrations. I CaL models subsequently have different assumptions about CDI and driving force. Furthermore, since I CaL gating models are often calibrated in the context of an AP model, several more parameters in the I CaL model will be sensitive to the chosen localization.
To gain an overview of how the spatial model affects the model of I CaL electrophysiology, we grouped the 73 I CaL models by the calcium concentrations assumed to affect I CaL gating O ð Þ and driving force δ ð Þ, leading to the 12 categories shown in Table 1. Note that not all models shown in this table fit neatly into the three subspace configurations shown in Figure 3: the studies by Asakura et al. (2014) and Himeno et al. (2015) both define an extra intermediate subspace between dyadic and bulk cytosol; the model by Heijman et al. (2011) has a very small subspace within the dyadic space, into which only I CaL flows; and the models by Noble et al. (1998), , and Restrepo et al. (2008) use different concentrations for the gating and driving force, which means they cannot be neatly assigned a position in Figure 3.
Interestingly, only the six models in rows 9 and 10 of Table 1 show a dependency on more than one Ca 2þ Â Ã , and in all six this arises from the assumption that there are two distinct populations of LCCs, located in separate subspaces. The models in rows 11 and 12 assume that gating depends on the Ca 2þ Â Ã in the dyadic nanodomain near the LCCs, but that the driving force is determined by the larger space surrounding it, which is the bulk cytosol in (Noble et al., 1998) and the subsarcolemmal space in  and Restrepo et al. (2008). None of the models we surveyed had the local and global components of CDI suggested by Tadross et al. (2008) in their O term.
The naming of spaces varies between studies, and sometimes overlaps or conflicts. In preparing Table 1 (and Figure 4) we classified a subspace as "dyadic" if its primary role was in calcium handling (e.g. if it contained LCCs and connected to the SR) while we used the term "submembrane space" only if several ion species flowed in and out of it, and it covered all or most of the (outer and t-tubular) membrane. An example where this conflicts is Sato et al. (2006), where we classified as "dyadic" a subspace termed "submembrane" by the authors.

| Regulating variables
In addition to the membrane potential V m ð Þ and (local) intracellular Ca 2þ Â Ã , I CaL models have been developed that are sensitive to several other variables. The most common of these is the CaMKII concentration, which is tracked in the model by Hund and Rudy (2004) and its descendants, and is used to estimate the fraction of CaMKII-phosphorylated channels in order to model CDF. The models by Matsuoka et al. (2003), Asakura et al. (2014), and Himeno et al. (2015) multiplied O by a term depending on the intracellular ATP concentration to reproduce an effect observed by Noma and Shibasaki (1985), while the model by Michailova et al. (2005) based itself on experimental work by O' Rourke et al. (1992) showing that these effects were due to [MgATP] rather than the free ATP concentration. Parasympathetic (ACh) inhibition of I CaL was included in the model by Pohl et al. (2016). No. Models O δ 1 Bassingthwaighte and Reuter (1972), McAllister et al. (1975), Liu et al. (1993), Demir et al. (1994), Lindblad et al. (1996), Zhang et al. (2000),  DiFrancesco and Noble (1985), Hilgemann and Noble (1987), Luo and Rudy (1994), Zeng et al. (1995), Dokos et al. (1996), Priebe and Beuckelmann (1998), Fox et al. (2002), Cabo and Boyden (2003), Matsuoka et al. (2003),  . All models shown here depend on voltage. I CaL models inside the green, blue, and red boxes depend on bulk cytosolic calcium, dyadic space calcium, and submembrane space calcium, respectively, and models inside the black-edged central box also depend on the concentration of CaMKII. Models marked with an asterisk depend on additional variables (as described in the main text) model to describe the LCCs and ryanodine receptors, leading to a dependency on the state of the SR (interestingly, this model uses the Ca 2þ Â Ã i to estimate the concentration in the dyadic space). Not many models surveyed in this review include effects of β-adrenergic stimulation on I CaL , but this is partly due to our selection criteria (see Section 2). An overview of studies focusing specifically on β-adrenergic stimulation (which almost always included effects on I CaL ) is given in the supplement to Heijman et al. (2011).
An overview of the dependencies of each I CaL model in this study is provided in Figure 4.

| Gating mechanisms
The feature that varies most between I CaL models is the set of equations used to describe the fraction of open channels, O. This part of the model determines how the current changes over time, that is, it describes the current kinetics or gating. In this section we present a classification of I CaL models into 15 distinct groups based on their gating equations. I CaL gating is commonly modeled using either Hodgkin-Huxley models (HHMs) or Markov models (MMs). In HHMs (Hodgkin & Huxley, 1952), the open probability is calculated as the product of several independent "gates", for example, one gate representing voltage-dependent activation, one representing VDI, and one representing CDI. The opening and closing of each gate is modeled using a single ordinary differential equation (ODE), although CDI gates are frequently assumed to be very fast so that their ODE can be replaced by an analytic equation for their steady state (e.g. the Hill equation). In MMs, the channel is postulated to be in one of a finite number of states, and transition rates are defined between the states that can depend on V m or Ca 2þ Â Ã : In contrast to HHMs, where each process (activation, VDI, CDI) is independent, MMs allow a more general structure where for example, activation and VDI are coupled (Rudy & Silva, 2006).
In the original and most common formulation, where the open probability of a HHM is the product of a number of independent gates, an equivalent MM can be written for each HHM (Keener & Sneyd, 1998). However, many I CaL models use extensions to the classical Hodgkin-Huxley (HH) scheme, for example, they introduce a non-inactivating fraction of channels (McAllister et al., 1975, see B in Figure 5), or they model the current as if arising from two separate channel families (Nygren et al., 1998). Such models have no obvious single MM equivalent. For examples of MM with and without a HMM counterpart and a graphical explanation of how to convert between the formalisms, see Figure 4 in Rudy and Silva (2006). Taking these equivalences into account, we identified 15 distinct gating models, as shown schematically in Figure 5. Where possible, we have shown models in an MM representation.
Each gating type is described briefly below. Some models define more than one copy of the same gating model but with changes in parameters or Ca 2þ Â Ã to account for LCC localization or phosphorylation; we do not focus on these properties in this section and are only concerned with the gating model structure.
Type A is the earliest and most straightforward gating mechanism, consisting of voltage-dependent activation d ð Þ and inactivation f ð Þ, modeled as independent gates that are not affected by Ca 2þ Â Ã (Bassingthwaighte & Reuter, 1972). Wilders et al. (1991) is a unique model within this group because its time constant of inactivation is not voltage-dependent but rather depends on the fraction of inactivation f . Type B extends type A by adding a noninactivating fraction of channels d 0 , leading to an open probability McAllister et al., 1975).
The type C gating mechanism extends type A by adding a CDI gate, usually written as f Ca . In type C1, this is modeled as an instantaneous process, so that the fraction of open f Ca gates is given by a Hill equation with Hill coefficient 2 in Luo and Rudy (1994) but 1 in the other C1-type models. Types C2 and C3 use an ODE to model the evolution of f Ca . C2 has CDI but calcium-independent recovery (DiFrancesco & Noble, 1985), while both inactivation and recovery are calcium-dependent in C3. Type C4 has an instantaneous f Ca gate, but also incorporates calcium-sensitivity in its gate for VDI, making it a dual VDI-CDI gate.
Gating type D is similar to type C in its HH form, except that the activation gate d is cubed (Nordin, 1993). In the MM structure, this equates to adding several closed states, as shown in Figure 5D.
An even larger Markov scheme is used by models of type E1, in which five steps are required to fully activate (similar to a fifth power in HH terms). In the original implementation this model was written as a combination between a HHM (for VDI) and a MM for activation and CDI (Jafri et al., 1998). The mechanism for CDI in this model is a switch to a "mode Ca" in which transitions to the open state are extremely slow (Imredy & Yue, 1994). A slight variation, which we still classified as E1, is given by Iyer et al. (2004) who removed the O Ca and OI Ca states. Several simplifications of the E1 scheme have been proposed, reducing it from 24 to 20 states (type E2, Michailova et al., 2005), 10 states (type E3, Greenstein & Winslow, 2002), and even three states (type E4, Hinch et al., 2004). Note that, to arrive at the representation of the E2 type shown here, we have reinterpreted the model by Michailova et al. (2005) as an MM and omitted the term accounting for [MgATP] effects and a constant factor x. The study by Li et al. (2010) extended the E4 type by F I G U R E 5 Models of I CaL gating. Each label in this figure (A, C1, C2, …) indicates a distinct gating mechanism, as described in detail in the text. Where possible, models are shown in a MM representation. MM state labels include Closed, Inactivated, Open, Activated, Uncovered, Resting, and Primed. Subscripts s and f are used to indicate states involved in "slow" and "fast" transitions, while the Ca subscripts indicate states involved in calcium-dependent transitions. The conducting states in each Markov structure are shown in bold. Where models have more than one conducting state, the total open fraction is found by adding the occupancy of all conducting states together. The open probability equation for type G is shown within the figure. In the HH equations, d and f denote activation and inactivation gates, respectively, and f Ca gates are calcium-dependent. re-adding VDI, creating the six-state type E5. Similarly, Asakura et al. (2014) extended E4 with a fourth state to create type E6.
Type F gating mechanisms again assume that I CaL has fast and slow modes (citing evidence from Imredy & Yue, 1994;You et al., 1995). In type F1, the fraction of channels in slow mode is determined by an instantaneous calcium binding process (Nygren et al., 1998) or fixed to a constant (Kettlewell et al., 2019), while in type F2 this process is modeled with an ODE (Pandit et al., 2001). Type F3 separates fast and slow kinetics from calcium by fixing the fraction of channels with fast VDI to 73% and introducing a separate CDI process (Pohl et al., 2016). A similar two-component model but without CDI was introduced by Inada et al. (2009), which we have labeled as type F4.
The type G gating mechanism (based on earlier work by Shirokov et al., 1993, for which we could not find the full equations) consists of three independent MMs for activation, VDI, and CDI (Matsuoka et al., 2003). In the MM for CDI, the binding of a Ca 2þ ion brings the system into a state where recovery from inactivation is much slower.
A new type H gating mechanism was introduced by Lovell et al. (2004), consisting of a four state MM with activation, VDI, and CDI. In this model, CDI can only occur when the channel is already in a nonconducting state. More complex MMs were introduced by Bondarenko et al. (2004), Faber et al. (2007), and  to create types I, J, and K respectively. The scheme in type I was inspired by structure-function relationships in a homology model based on Shaker B voltage-gated K þ models. Type J includes a fast and a slow VDI, to reflect gating current experiments by Ferreira et al. (2003). The Type K model incorporates ideas from Soldatov (2003) and Cens et al. (2006), who argued that CDI is a voltage-dependent process that is greatly accelerated when Ca 2þ binds to the channel (see also Pitt et al., 2001). Its Markov scheme is reduced by combining two open states into one, and reducing the five-state activation to a three-state process (Rose et al., 1992).
The type L1 gating mechanism by Hund and Rudy (2004) consists of activation, fast and slow VDI, and fast and slow CDI. In this model, CDI depends on Ca 2þ Â Ã d but also has an I CaL -dependence, to approximate the raised Ca 2þ Â Ã at the cytosolic channel mouth (this was based on the model by Hirano & Hiraoka, 2003, for which we could not find the full equations). The activation variable is raised to a time-and voltage-dependent power n (modeled as an ODE) which accounts for a "fast voltage-dependent facilitation" described by Kamp et al. (2000). In type L2, this is simplified to use a constant n ¼ 1 (Aslanidi, Stewart, et al., 2009), and type L3 simplifies this further by combining fast and slow VDI into a single gate (Hund et al., 2008). All type L models have a fast CDI gate f Ca f that accounts for CDF by CaMKII-dependent phosphorylation.
The study by Ten Tusscher and Panfilov (2006) extended their previous C3 type model with an extra gate, leading to type M. Note that this can also be seen as a simplification of type L2.
The MM scheme for Type N is somewhat similar to C2, but with the important difference that it has two conducting states (so that this MM can no longer be reduced to an equivalent HHM). Like type K, this gating model is based on the idea that Ca 2þ binding in CDI "removes a brake" on a voltage-dependent process (Pitt et al., 2001). This is very similar to the "mode-switching" idea underlying E1. Heijman et al. (2011) uses a variation of this scheme, in which the rates of CDI and VDI are modulated by [CaMKII]. In addition, this model uses a second copy of the same MM structure but with modified rates to represent PKA-phosphorylated (β-adrenergically stimulated) LCCs. Another variation is used by Bartolucci et al. (2020), who borrow two state variables from the model by O'Hara et al. (called n and j in the publication, but shown as f Ca and r, respectively in Figure 5O) which are used to determine the transition rate into the "Ca" states. As in the model by Heijman et al., two copies of the same model are maintained, but this time the phosphorylated population is CaMKII-rather than PKA-phosphorylated.
Type O gating, introduced in O'Hara et al. (2011), splits the channels into several fractions. First, a Ca 2þ Â Ã dependent variable (called the "n-gate" in the original publication but f Ca in our notation) splits the channels into a fraction in "VDI mode" and a fraction in "CDI mode." As in types K and N, the underlying assumption is that Ca 2þ binding (modeled by the "n-gate") brings the LCCs into a state where a voltage-dependent "CDI" process can occur Pitt et al., 2001). In turn, each fraction (VDI mode and CDI mode) is split into a fast and slow fraction, with a fixed ratio (x in Figure 5) between fast and slow "VDI-mode inactivation" and a voltage-dependent ratio (y in Figure 5) between fast and slow CDI-mode inactivation. All three type O models in this review include CaMKII-phosphorylation via a second copy of the model with altered rate equations. Table 2 shows how each model fits into this classification. In addition to the gating type, it lists each model's cell type, species, and the simulated temperature (where stated). The final column provides an indication of the major data sources used in the model's construction. For the earlier generation of models this data source is often hard to establish, as parameters were commonly set by hand to create a model that was acceptably close to a wide range of observations T A B L E 2 Models classified by their gating type, with additional information of cell type, species, temperature (T), and the major data sources used in each model's construction.

| Driving force
Two types of driving force are common in I CaL models. An Ohmic driving force for a current carried by a single species takes the form where E rev is the reversal potential, at which the current reverses direction. For currents carried by a single species, the reversal potential can be calculated using the Nernst equation where R is the universal gas constant, T is temperature, F is the Faraday constant, z x is the valence of a single ion of species X (2 for Ca 2þ Â Ã Þ, and X ½ o and X ½ i are the external and internal concentrations of X, respectively. Although commonly used in the form given above, a more accurate version uses γ o X ½ o =γ i X ½ i , where γ o and γ i are dimensionless activity coefficients that account for nonideal behavior of the solutes. The units of the gas constant are frequently chosen as mJ=K=mol to yield a reversal potential in mV. A model with an Ohmic driving term is written A more complex model of the electrochemical driving force is given by the Goldman-Hodgkin-Katz (GHK) flux equation (Goldman, 1943;Hodgkin & Katz, 1949): Because its derivation involves assuming a constant electrical field throughout the channel, the GHK equation and its derivatives are also known as "constant field theory." Assuming concentrations in mM ¼ mol=m 3 and using matching units for V and RT=F, a model with a GHK driving term is written as I ¼ P Á O Á δ GHK , where P is in L=ms (¼ m 3 =s) for a current in A, in L=F=ms for a current in A=F, or in cm=s for a current in μA=cm 2 :

Models
Cell Species T Major data sources A number of early studies used a modified GHK equation, which accounts for the (hypothesized) presence of a charged particle near the channel mouth, causing a change, V 0 , in the voltage across the channel (Frankenhaeuser, 1960): So far, we have assumed a current carried exclusively by Ca 2þ through a perfectly selective channel. For this situation, δ Ohmic , δ GHK , and δ GHK 0 all predict a current that reverses at the potential given by Equation (3). For typical concentrations, this would lead to a current that reverses in the range 50-150 mV, which is significantly higher than experimental estimates of 40-70 mV. If, however, we assume that the channel is even very slightly permeable to Na þ and K þ , then the much higher internal concentrations of those ions compared to Ca 2þ Â Ã have a significant effect on the reversal potential.
A simple model for currents with multiple carriers can be made by assuming that each species travels through the channel without affecting the others, so that the total current is simply a sum of single-species currents (independent flux assumption). Equations for E rev can still be derived with this assumption (see e.g., Campbell et al., 1988;Keener & Sneyd, 1998), but the predicted reversal potentials from Ohmic and GHK models are no longer the same. In both equations, the reversal potential of a multi-species current depends strongly on the ratio between the species' permeabilities P X (or, if an Ohmic term is used, their conductances g X ). A study by Campbell et al. (1988) measured the reversal potential of I CaL in bull-frog atrial myocytes and compared this to predictions using different ratios of P K : P Ca . They found that good agreement with experiment could be found using a GHK equation (Equation 4), provided the channel was highly selective to Ca 2þ Â Ã (>95% of the I CaL is carried by Ca 2þ Â Ã at 0 mV). To compare different models' driving terms, we set internal and external concentrations to fixed levels, and plotted the driving term as a function of voltage in Figure 6. We used Na (Note that these values may differ from the concentrations that the models were designed for, the concentrations used in the experiments they were calibrated with, or the physiological concentrations in real cells.) In models where I CaL is carried by multiple species, we calculated the "net driving term" by summing up the contributions of each component, and weighting by the component's permeability relative to the calcium component (e.g. the calcium component has weight P Ca =P Ca ¼ 1, the potassium component has weight P K =P Ca , etc.).
The left panel of Figure 6 shows δ for the models in this study that used an Ohmic driving term (models for which we could not find the full equations are omitted). Rather than using the Nernst equation, most of these models set a constant reversal potential in the range of 45-70 mV. As expected, models that do use the single-species Nernst equation without any modifications (Corrias et al., 2011;Priebe & Beuckelmann, 1998) show a reversal potential well outside of the physiological range. The right panel in Figure 6 shows δ for models with a GHK (or GHK 0 ) driving term. Most of these models do not show a true reversal of current, with the driving force approaching zero for high voltages. Only models with a high permeability for sodium and potassium ions (DiFrancesco & Noble, 1985;Hilgemann & Noble, 1987) show a positive driving force for high voltages. Finally, the model by Liu et al. (1993) uses neither an Ohmic nor a GHK driving term, but an exponentially increasing flux and is shown in the right panel.
A list of models with Ohmic driving terms is given in Table 3. For models using a fixed E rev the numerical value is shown, while the equation used for E rev is shown for the remainder. None of the models with Ohmic driving terms in this study explicitly modeled potassium or sodium components of I CaL (although their existence may be implicitly acknowledged in the models with a fixed reversal potential).
A similar list for models with GHK and GHK 0 driving terms is given in Table 4, along with the offset V 0 (if used), the activity coefficients for internal and external calcium, and the ratio between the permeability coefficients of potassium and sodium with respect to calcium.
The earliest studies in this table all use a voltage shift parameter V 0 ¼ 50 mV, but this form is rarely seen in later models. Models with a large V 0 also assume a large contribution of sodium and potassium ions, consistent with the analysis by Campbell et al. (1988).
The model by Luo and Rudy (1994) abandons the use of V 0 and introduces the values γ i ¼ 1 and γ o ¼ 0:341 for the activity coefficients of internal and external calcium. These values can be calculated using the Pitzer equations (Pitzer & Mayorga, 1973), while the frequently used value γ i ¼ 1 appears to arise from assuming that low internal Ca 2þ Â Ã leads to near-ideal behavior. Instead of using constant values, the model by Tomek et al. (2019) calculates the activity coefficients using the Davies equation (Davies & Malpass, 1964), which is a precursor to Pitzer's work (although the first version of this model erroneously used a natural logarithm instead of a base-10 logarithm for the calculation, see https://bit.ly/3tqD4gP). Activity coefficients for sodium and potassium are usually taken to be 0.75 internally and externally, and are not shown in the table.
Although there is a common consensus that I CaL is carried by multiple ion species, the ratio of permeabilities between the different species varies greatly between models, and several studies have cited the low permeability to potassium and sodium as a reason to omit these components from the model altogether (avoiding their impact on the reversal potential by using a fixed rather than a calculated value). Estimates of LCC selectivity have gone up over time: in general, older models assume larger contributions from sodium and potassium, while more recent models assume a much more selective channel. The high selectivity of LCCs is especially evident when comparing values from Table 4 to estimated selectivity ratios of other channels, for example, P Na : P K ≈ 1:140 for I Kr (Sanguinetti et al., 1995), while P K : P Na has been estimated as 1:10 for I Na (Amin et al., 2005).
Modelers starting with Jafri et al. (1998) have attempted to incorporate evidence that monovalent permeation goes down when there is a significant Ca 2þ influx-in other words, that the independent flux assumption does not hold. In these models, which consider only a calcium and a potassium component, the permeability to potassium is a function of the calcium driving force, so that P K goes down when P Ca is high (Cortassa et al., 2006;Fox et al., 2002;Iyer et al., 2004;Jafri et al., 1998;Winslow et al., 1999). This is indicated with f V, Ca ð Þin Table 4. Ca 2þ influx, on the other hand, seems able to continue even when net I CaL is outward (Zhou & Bers, 2000). (It is also interesting to speculate whether Ca 2þ flux might reverse locally, for example, when a high Ca 2þ Â Ã is reached in the dyad.) Finally, some of the models in Table 4 use a slight variation of the (normal or modified) GHK equation. The model by Jafri et al. (1998) and many of its descendants, (Cortassa et al., 2006;Iyer et al., 2004;Michailova et al., 2005;Rice et al., 1999;Winslow et al., 1999) use a constant internal calcium activity γ i Ca 2þ Â Ã i ¼ 0:001 mM (Smith, 1996). The model by Ten Tusscher and Panfilov (2006)

| QUANTITATIVE COMPARISON
In this section, 60 models are reviewed and compared functionally (Cooper et al., 2011) by simulating voltage-clamp protocols and studying the predicted I CaL response. Voltage-step experiments are simulated to compare activation, VDI Ohmic driving force (mV) GHK driving force (mC/μl) F I G U R E 6 Driving term as a function of voltage for models with an Ohmic driving term left: and a (normal or modified) GHK driving term right. In the left panel, the majority of models are plotted in gray and are shown to form a tight cluster, reversing at around 60 mV. The lower and upper bounds for this cluster are marked by McAllister et al. (1975) and Kurata et al. (2002), which are indicated in color, along with three outliers (Beeler & Reuter, 1977;Corrias et al., 2011;Priebe & Beuckelmann, 1998) that do not show a reversal potential within the physiological range. In the right panel, most currents are again shown in gray, but a selection of models are highlighted in color to show the range of behaviors observable by varying the activity coefficients and voltage offset V 0 (DiFrancesco & Noble, 1985;Hilgemann & Noble, 1987;Hund & Rudy, 2004;Luo & Rudy, 1994;Ten Tusscher & Panfilov, 2006;Tomek et al., 2019). An interactive view of these results is available at https://chaste.cs.ox.ac.uk/q/2022/ical/fig6. and recovery, and CDI. To compare the predicted I CaL during an AP, we perform simulations where either the AP (APclamp) or AP and CaT (AP-CaT-clamp) are "clamped" to a predetermined time-course.  Bassingthwaighte and Reuter (1972) 49.5

| Ionic concentrations
Note: The final two columns provide representative E rev values during the AP ("systole") and in the refractory phase ("diastole"), for models in which E rev is given by an equation and for which a full AP model was available.
T A B L E 4 GHK and GHK 0 parameters in models using GHK driving terms (Equation 5).

| Simulation methods
As the basis for our simulations, we used CellML 1.0 or 1.1 files (Hedley et al., 2001) containing either a single I CaL model, or an I CaL model embedded into a larger model of the AP. When simulating basic voltage-clamp protocols (activation, VDI, recovery, CDI) AP models were reduced to I CaL -only models by fixing the values of V m and all external and internal ionic concentrations to predetermined levels (see Sections 3.2.2 and 3.1). In the AP-and AP-CaT-clamp simulations, some non-I CaL model features are preserved (e.g. the internal calcium and CaMKII dynamics), so that these simulations require a full AP model (see Section 3.6). Where possible, model implementations were obtained from the Physiome Model Repository (PMR, Yu et al., 2011). Thirty-five models obtained this way could be used directly, while a further thirteen needed corrections to either units or equations. To limit the scope of this correction work, five models were reduced to I CaL -only models in the process. Eight models were corrected in AP form, after which the new versions were uploaded to the PMR. A file for the model by P asek et al. (2006) was also obtained from the PMR. This model has a separate variable for the voltage across the Ttubular membrane, distinct from the "main" transmembrane potential. To compare the channel kinetics consistently, we modified this model manually to set this secondary voltage equal to V m .
For 11 models no CellML files were available, so new implementations were created based on published equations or code. An This resulted in a set of 60 model files that could be used to simulate I CaL in isolation, and a subset of 44 files that provided a full model of the AP.

| Software
Simulations were run using the Cardiac Electrophysiology Web Lab (Cooper et al., 2016;Daly et al., 2018). This is a resource developed to allow the characterization and comparison of cardiac electrophysiology models in a wide range of experimental scenarios. It uses a novel protocol language that allows the user to separate the details of the mathematical model from the experimental protocol being simulated. We also performed simulations for calcium sensitivity (see Section 3.5) by solving the model equations using CVODE version 3.1.2 as incorporated in Myokit (Clerx et al., 2016) version 1.28.4, using Python version 3.7.3. All model and simulation files are available under an open source BSD 3-Clause license, and can be obtained from https://github.com/CardiacModelling/ical-model-review (doi: 10.5281/ zenodo.6653898).

| Model annotations and modifications
To run simulations in the Web Lab, models are encoded in CellML format and annotated with metadata terms. These terms allow variables in different models with the same physiological meaning to be identified, read out in units of choice, and manipulated. This is a useful feature for this study as the 60 models do not use consistent naming. The following variables were annotated in all models: time, membrane potential V m ð Þ, the total I CaL , driving term δ ð Þ, and the open probability O ð Þ. When annotating driving term variables, distinct metadata terms were used for Ohmic and GHK formulations. When annotating I CaL in models that split the current into multiple components (by ionic species, phosphorylated/non-phosphorylated fraction, or biological subcompartment) we selected the variable that represents the sum over all these aspects. For this purpose, a new variable was often introduced in the CellML files. Similarly, O is sometimes split into components with or without phosphorylation, δ can be split by ionic species, and both O and δ are occasionally split by the biological compartment. For all such cases, we introduced and annotated new O and δ variables as the weighted average over all these components (see e.g. the annotations in our file Bartolucci et al., 2020). This weighting was made based on the fraction of phosphorylated channels, the relative permeability of each ionic species with respect to calcium, and/or the fraction of channels in each biological compartment. For example, where I CaL was calculated as a linear combination of a calcium, potassium, and sodium component a weighted term for a GHK driving force and current carried by three ion species was calculated as: where δ Ca , δ K , and δ Na are the driving force terms calculated for each species separately. Finally, we annotated the variables representing intra-and extra-cellular concentrations of Ca 2þ , Na þ , K þ , and Cl À . For external concentrations, most models use a single constant parameter, but some make a distinction between a constant "bath concentration" and a time-varying "cleft concentration" near the cell. All three types of external concentration ("bath", "cleft", or simply "external") were annotated and subsequently manipulated by the Web Lab during simulation. Similarly, we annotated (and later manipulated) intracellular concentrations for bulk cytosol, submembrane space, and dyadic space concentration.

| Voltage dependence of activation
To show model predictions of voltage-dependence of activation, a voltage-clamp protocol was adapted from O' Hara et al. (2011). This is a typical example of the type of protocol commonly used to study I CaL activation experimentally, but there are some notable differences between these experiments and our simulation. First, the simulated voltage is set exactly to the values dictated by the protocol (corresponding to an ideal rather than a realistic voltage clamp, see . Second, internal and external concentrations are held constant during the simulation (see Section 3.1), so that CDI is kept at a fixed (but model-dependent) level. When using AP models, this "clamping" of voltage and concentrations also serves to remove all interaction with other ionic currents, so that this is a simulation of I CaL in isolation.
The activation protocol consists of repeating units called sweeps. Each sweep starts with a 20 s segment where V m is kept at the holding potential of À90 mV, followed by a 120 ms step to the test potential, denoted by P1 in Figure 7 (left). The test potential P1 is À60 mV for the first sweep and increases by 2 mV per sweep up to þ60 mV in the final sweep (for clarity, only a subset of these fine-increment sweeps is shown in the figure). In our simulations, current was recorded during P1 to measure the peak I CaL (denoted I) at each sweep (Figure 7, center) which was then normalized as I ¼ I= j min I j. Next, the normalized peak current, I, was plotted against the test potential for all models, leading to the normalized I-V curves shown in Figure 7 (right). We also calculated the midpoint of activation (V 0:5 ) for each model, by finding the point where the normalized I-V curve first crosses À0.5.
Most I CaL models predict that channels begin to activate at membrane voltages around À40 mV and that I CaL reaches its maximum magnitude between À20 and þ 20 mV. The V 0:5 varies between À48 and þ15:3 mV for the models in this study while the median V 0:5 is À13:5 mV. By comparison, examples from the experimental (patch-clamp) literature show that LCCs start to activate at around À40 mV, that peak I CaL is observed between about À5 to 10 mV, and that V 0:5 lies between À25 and À10mV (Magyar et al., 2000;Pelzmann et al., 1998;Van Wagoner et al., 1999).
In Figure 8, we take a closer look at the outlier models not included in Figure 7. The left panel shows the I-V curve for all 60 models, including outliers, which are highlighted in color. This includes the two oldest models in our quantitative review (Beeler & Reuter, 1977;McAllister et al., 1975), both of which are based on tissue rather than isolated-cell measurements and so a difference in their activation kinetics is not surprising. However, some more interesting effects are observed when we dissect the currents into open probability and driving force, as is done in the right panel which shows normalized open probability versus voltage (defined as O ¼ O= j max O j). Of the four outliers, only the McAllister et al. (1975) model has an unusual O-V curve, causing peak I CaL to be reached at a lower voltage than other models. The models by Beeler and Reuter (1977) and Corrias et al. (2011) have more typical O-V curves, but exhibit unusually large currents at higher voltages (left panel) due to their unusually large driving force (see Figure 6 and Table 3). Interestingly, the model by Priebe and Beuckelmann (1998) has a similar driving term (Figure 6), but this is compensated by an outlier O-V curve, leading to a regular I-V curve (Figure 8, but see also Figure 17). Finally, the model by Liu et al. (1993) has an unusual I-V curve not due to its O-V curve but due to its exponentially increasing driving term. From this dissection of I CaL into O and δ we can see how irregularities in one aspect of a model (e.g. O) can sometimes be compensated by irregularities in another (e.g. δ).

| Voltage-dependent inactivation
Next, we illustrate VDI using simulations with a protocol adapted from Magyar et al. (2000), shown in Figure 9 (left). As with activation, we show predictions for ideal voltage clamp and perfectly buffered ionic concentrations leading to a fixed but model-dependent level of CDI, so that there are important caveats when comparing these simulations to experimental data. Each sweep in the protocol starts with 20 s at a holding potential of À90 mV, followed by a 500 ms step (P1) to a variable conditioning potential and a 120 ms step (P2) to a fixed test potential of 0 mV. The conditioning potential at P1 is À90 mV for the first sweep and increases by 2 mV per sweep up to 30 mV for the final sweep (a subset of which is shown in the figure). Current was recorded during P2 to measure peak current (I, Figure 9, center), which was normalized (I ¼ I= min I ) and plotted against conditioning potential to construct the I-V curve shown in Figure 9 (right).
This figure shows that most models predict inactivation to set in at around À60 mV or greater. The conditioning potential after which the peak I CaL is halved V 0:5 ð Þ varied between À59:3 and þ10 mV for the models in this study while the median À29:9 mV. By comparison, examples from ventricular cells at physiological temperature show LCCs Peak current (normalized) F I G U R E 7 Left: Voltage-clamp protocol used to characterize activation (see main text for details). Center: The current predicted by the Grandi et al. (2010) model, simulated with fixed internal and external concentrations and a constant level of CDI. Traces from successive sweeps are overlaid and color-coded to match the P1 steps in the first panel. The peak I CaL at each voltage step is marked with a gray square. These values are used to construct the I-V curves on the right, as shown in the animation online at https://github.com/CardiacModelling/ ical-model-review. Right: Normalized peak I CaL versus step voltage for 56 (out of 60) models (the direction of the current was preserved during normalization). The black dashed line shows the I-V curve from Faber et al. (2007), which is close to the median response. The models by McAllister et al. (1975), Beeler and Reuter (1970), Liu et al. (1993), and Corrias et al. (2011) showed outlier responses, and are plotted separately in Figure 8. An interactive view of these results is available at https://chaste.cs.ox.ac.uk/q/2022/ical/fig7. beginning to inactivate after a conditioning voltage of around À60 mV, with V 0:5 lying between À30 and À60 mV (Kim et al., 2016;Li et al., 1999).

| Time-course of recovery
Model predictions of the time I CaL takes to recover from inactivation are shown using a protocol adapted from Fülöp et al. (2004) (Figure 10, left). Each sweep in this protocol starts with 10 s at a holding potential of À90 mV, followed by a 500 ms step (P1) to a conditioning potential of 0 mV. The potential is then set back to À90 mV for a duration that varies from sweep to sweep, followed by a 120 ms second step (P2) to 0 mV. The interval duration is 10 ms for the first sweep and increases by 10 ms per sweep up to 50 ms, after which it increases by 25 ms per sweep up to 1000 ms (a subset of the sweeps are shown in the figure).
Peak I CaL was recorded during P1 (I conditioning ) and P2 (I test ) for all sweeps (Figure 10, center). This was normalized as I ¼ I test =I conditioning and the normalized peak I CaL (I ) was plotted against interval length to construct the current-time curves shown in figure 10 (right).
An exponential function 1 À exp Àt=τ ð Þ was fitted to each model's current-time curve to yield the time constant of recovery from inactivation τ. This τ varies between 1.4 and 995.7 ms for the models in this study while the median τ is 49.2 ms. This is much faster than examples from ventricular cells at physiological temperature, in which τ ranged from 50 to 700 ms (Kim et al., 2016;Li et al., 1999).

| Calcium-dependent inactivation
To study the effects of CDI, we repeated the inactivation protocol of Figure 9 at constant internal (bulk, dyadic, and submembrane) Ca 2þ Â Ã of 1 pM, 0.1 μM, and 300 μM. The lowest value represents a situation with near to no Ca 2þ , while the other two values represent the normal range of calcium concentrations in the dyadic space during an AP (Fearnley et al., 2011). Postprocessing proceeded similarly to before: peak I CaL was identified at each sweep, but normalization was performed using the maximum peak I CaL measured at 1 pM, so that the magnitude of the low, moderate, and high Ca 2þ Â Ã curves can be compared. Figure 11 shows the resulting I-V curves for the three experimental conditions. For ease of presentation, results have been split into three categories: insensitive or weakly sensitive to Ca 2þ Â Ã ; mildly sensitive; and strongly sensitive to Ca 2þ Â Ã : To assign a model to a sensitivity class, we defined a "root mean squared difference" measure between the I-V curves for extreme concentrations as Peak current (normalized) Peak probability (normalized) F I G U R E 8 Left: I-V curve and right: Normalized peak-open probability plotted against test potential for all 60 models. The four outliers and the model by Priebe and Beuckelmann (1998) are highlighted in color and discussed in the main text; the others are shown in gray. An interactive view of these results is available at https://chaste.cs.ox.ac.uk/q/2022/ical/fig8.
where n is the number of points in each I-V curve. We classified a model as weakly sensitive if RMSD < 0:1 (9 models), mildly sensitive if 0:1 ≤ RMSD ≤ 0:5 (17 models), and strongly sensitive if RMSD > 0:5 (34 models). Of the 9 models in the lowest sensitivity category, only Noble et al. (1991) and Pandit et al. (2001) account for Ca 2þ Â Ã in their equations at all; the remaining 7 have no Ca 2þ Â Ã terms in either O or δ (see Table 1, Row 1).
3.5.1 | Ca 2þ Â Ã sensitivity and localization To further investigate the differences in predicted sensitivity to Ca 2þ Â Ã we used a simple protocol consisting of a À90 mV holding potential and a single brief step to 0 mV (i.e. the first sweep shown in Figure 9). Starting from an internal (bulk, dyadic, and submembrane) Ca 2þ Â Ã of 1 pM (i.e., almost no Ca 2þ ), we repeated the experiments with increasing concentrations until a 50% drop in the peak I CaL was observed. Note that each repeat was started from the same initial conditions: no build-up of CDI (or CDF) occurred between repeats. We termed the resulting concentration the CDI 50 Ca 2þ Â Ã . Using this procedure, the CDI 50 Ca 2þ Â Ã was obtained for 49 out of the 51 mildly and strongly Ca 2þ Â Ã sensitive models, as shown in Figure 12. The models by Wilders et al. (1991) and Nygren et al. (1998) did not produce a 50% drop at any of the levels we tested. The calculated CDI 50 is shown in Figure 12 (left) and varies across different models from 0.16 μM (10 th percentile) to 338 μM (90 th percentile) while the median CDI 50 is 83.6 μM.
In a more realistic situation, where Ca 2þ Â Ã is not held constant throughout the cell, Ca 2þ flowing into a subspace will induce a change in local Ca 2þ Â Ã with a rate that depends inversely on the subspace volume V local (assuming there is no buffering or diffusion of Ca 2þ ). We may therefore expect models in which I CaL flows into large subspaces to compensate for the slower local Ca 2þ Â Ã -change with a more sensitive CDI, and vice versa. This is explored in Figure 12 (right), in which the CDI 50 concentration is plotted against V local . For models in which LCCs are present in multiple biological subspaces, the volume was calculated as a weighted average: For example, 80% of the LCCs in the model by Tomek et al. (2019) lie in the dyadic space, whereas only 20% lie in the submembrane space, leading to a weighted average volume V local ¼ 0:8V d þ 0:2V s . The models by Beeler and Reuter (1977) and  were omitted as the volume of their localization space could not be determined from the manuscripts. Figure 12 (right) shows a weak correlation between the predicted CDI 50 and the corresponding V local , suggesting that modeled CDI is to some limited extent influenced by the size of the compartment that I CaL flows into. This is another example where we observe that one aspect of the model can compensate for another.

Peak current (normalized)
F I G U R E 9 Left: Voltage-clamp protocol used to characterize voltage-dependent inactivation (see main text for details). Center: Current response during P2 as predicted by the Grandi et al. (2010) model, with constant ionic concentrations and a fixed level of CDI. Traces from successive sweeps are overlaid and color-coded to match the P1 steps in the left panel. The peak I CaL at each voltage step (marked with gray squares) is used to construct the I-V curves on the right, as shown in an animation online at https://github.com/CardiacModelling/icalmodel-review. Right: Normalized peak I CaL versus step voltage curve for 60 models. The I-V curve from Kurata et al. (2002), which is close to the median response over all models, is highlighted with a black dashed line. An interactive view of these results is available at https:// chaste.cs.ox.ac.uk/q/2022/ical/fig9.

| AP and AP-CaT clamp predictions
The experiments simulated so far use voltage-step protocols, which bring out a wide range of I CaL kinetics but are very different from the voltage signals experienced by LCCs in situ. The final aspect we review is the models' prediction of the I CaL current elicited by a predefined AP. As before, we simulate an idealized voltage-clamp, but this time we allow internal Ca 2þ Â Ã to vary. Because models including CDI differ in their geometrical layout, there is no obvious "best" way to do this. Assumptions about internal calcium dynamics will impact the (estimated or calibrated) parameters governing I CaL kinetics, and to bring these to the foreground we may choose to let the internal calcium concentrations in these simulations vary as defined by the model equations. If, however, we are interested only in the models' I CaL equations, it is preferable to force an identical CaT at the internal channel mouth in each simulation.
To address these conflicting aims, we performed a series of three complementary simulations. First, we simulated an AP-clamp in which we allowed the internal calcium concentration(s) to vary according to model equations (AP clamp 1). Second, we performed a dual AP-CaT-clamp, in which all internal Ca 2þ Â Ã (including the subspaces) were clamped to the same precalculated bulk cytosolic CaT (AP clamp 2). Third, we performed a combined AP-CaT-clamp in which bulk cytosol, submembrane space, and dyadic space Ca 2þ Â Ã were clamped to precalculated transients appropriate to the subspace (AP clamp 3). We chose to use the AP and CaT measurements from the Grandi et al. (2010) human ventricle model during 1 Hz steady pacing (Figures 3 and 13), as this is the earliest human model in which all three subspaces are defined. (Note that these strategies are meant to illustrate model characteristics and are not necessarily experimentally feasible.) Because AP clamp 1 requires an AP model, this simulation was only performed for the 44 models for which a CellML file for the full AP model was available. For AP models with epi-, endo-, and midmyocardial variants, the epicardial version was used. AP clamps 2 and 3 were performed for all 60 models. In all three simulations, I CaL was normalized with respect to the total charge carried: I t ð Þ ¼ I t ð Þ= j R t I t ð Þdt j : The resulting normalized currents are shown in Figure 14. 3.6.1 | Variation in gating across models is not easily explained All three AP clamps in Figure 14 show considerable variation in predicted (normalized) I CaL between models. Some of this variation might be explained as a result of differences in for example, species and cell-type, or more mundane factors such as the year of publication. To test this assumption, we divided the 60 models into groups according to several criteria and investigated whether this resulted in groups with reduced variation. We chose the "compromise" protocol AP clamp 3 for this purpose. To remove effects of varying driving terms and maximal conductances, we compared the Peak current (normalized) F I G U R E 1 0 Left: Voltage-clamp protocol to characterize recovery from inactivation, with a variable time interval between P1 and P2 of 10, 20, 30, 40, 50 ms and increment steps of 25 ms thereafter until 1000 ms (see main text for details). Center: Example of I CaL predictions (constant Ca 2þ Â Ã ) during P2, from the model by Grandi et al. (2010). The colors of each line match the color shown above the interval length in the left panel. Traces from successive sweeps are overlaid and the peak I CaL is marked at each sweep. This is then used to construct the curves on the right, as shown in an animation online at https://github.com/CardiacModelling/ical-model-review. Right: Normalized peak I CaL versus time curves for all 60 models. The prediction from Fink et al. (2008) is close to the median response across all models and is highlighted with a black dashed line. An interactive view of these results is available at https://chaste.cs.ox.ac.uk/q/2022/ical/fig10 predicted open probability (O) rather than I CaL . To give an indication of the variation within a group, we used the "normalized root mean square difference" (NRMSD) of the observed open probabilities of all models in the group, defined as:
Here, O ij represents the open probability measured in response to the "AP-CaT clamp 3" protocol at the i th time step for the j th model, m is the total number of time points, and n is the total number of models within the group.
If grouping by a particular category reduced variability, we would expect the NRMSD within each group to be lower than the NRMSD of the original group-of-all-models NRMSD total ¼ 0:139 ð Þ . This can be quantified by defining a relative measure NRMSD relative ¼ NRMSD group =NRMSD total . A grouping can be considered to be better than baseline if its NRMSD relative is less than one.
First, we investigated the effects of grouping models chronologically by placing each model into one of four categories depending on the year of publication. This is shown in Figure 15, where we can observe some convergence in model predictions over time, indicated by a decreasing NRMSD relative . Although this is encouraging, some degree of convergence should be expected as more recent models are often refinements of previous models in the same age category. Large variability should also expected from the models published before the 1990s, which depend on data obtained from a wider range of experimental techniques.

Bulk cytosol Dyadic space
Volume of localization subspace (pl) F I G U R E 1 2 Left: CDI 50 for 47 out of the 49 mildly or strongly Ca 2þ Â Ã sensitive models (see main text for details). Right: CDI 50 plotted against the volume of the subspace assumed to contain the LCCs, for the same 47 models. Colors indicate the subspace containing the channels, or the subspace containing the majority of channels in situations where LCCs are assigned to multiple spaces.
F I G U R E 1 3 The AP waveform used in AP-and AP-CaT-clamp simulations in Figures 14-16. This action potential was derived by simulating steady 1 Hz stimulation with the model by Grandi et al. (2010).
Next, we grouped models according to the approach used to model how the calcium concentrations affect the open probability, depicted by the rows in Table 1. On this basis, Figure 16 shows the models grouped into six classes: no Â Ã near the LCC have at least $20% less disagreement among the predictions across models than the baseline. This indicates that the modeled localization of the LCCs plays an important role in determining model predictions. Note that class 5 contains only a single model (Tomek et al., 2019) so that no NRMSD relative can be calculated. Similarly, all three models in class 6 (Grandi et al., 2010(Grandi et al., , 2011Shannon et al., 2004) were developed by the same research group, which explains the very low NRMSD in this group.
Of all tested criteria, grouping by local Ca 2þ Â Ã near the LCCs exhibited the most consistent predictions.

| WHAT NEXT? OPEN QUESTIONS AND CHALLENGES
In the previous sections of this review, we identified 73 models of mammalian I CaL , mapped their historical relationships, and analyzed the qualitative differences in their constituent parts. Driving-force models differed in their choice of Ohmic, GHK, or modified GHK equations as well as their choice of charge carrier(s). Fifteen major gating schemes were identified, with a further thirteen subtypes bringing the total to twenty-eight. Modeling of LCC localization was much less diverse, with nearly all recent (adult myocyte) models using a cell layout where LCCs emerge into a dyadic subspace. Using simulation, we reviewed the predictions of the current from 60 models using three common voltageclamp protocols, as well as more complex AP and AP-CaT-clamp protocols that illustrated CDI. Great variability was observed between model predictions, which was somewhat diminished when grouping the models by their assumptions about LCC localization, but no grouping could be found that significantly reduced (or explained) this variability.  Grandi et al. (2010). Results for AP clamp 1 are limited to 44 cases where an AP model CellML file was available, whereas all 60 models are shown for the remaining AP clamps. Predictions from Grandi et al. (2010) are highlighted in red. An interactive view of these results is available at https://chaste.cs.ox.ac.uk/q/2022/ical/fig14 terms were canceled out by unusual gating, and calcium-sensitivity of CDI correlated with the volume of the subspace LCCs were assumed to occupy.
These 73 models, along with the further models we could not cover in this review, represent a huge effort and a great achievement by the cell electrophysiology community. It is perhaps surprising, however, that no consensus or Open probability F I G U R E 1 5 Open probability in response to AP clamp 3, grouped according to date of publication, shows a degree of convergence over time. Blue curves in each panel correspond to I CaL predictions of models within the corresponding period, while the remaining predictions are shown in gray. (Note how some models allow O to be greater than 1, so that it is not strictly a probability in these models.)

Open probability
Open probability

Bulk cytosol
Submembrane space F I G U R E 1 6 Open probability (O) in response to AP clamp 3, grouped by local Ca 2þ Â Ã , as shown by the rows in Table 1 where each row represents an O or driving term dependence on different local Ca 2þ Â Ã : The first panel "no Ca 2þ Â Ã " shows the models from rows 1 and 2, followed by "bulk cytosol Ca 2þ Â Ã " (rows 3 and 6); "submembrane space Ca 2þ Â Ã " (rows 4 and 7); "dyadic space Ca 2þ Â Ã " (rows 5, 8, 11, and 12); " Ca 2þ Â Ã i and Ca 2þ Â Ã d " (row 9); and " Ca 2þ Â Ã s and Ca 2þ Â Ã d " (row 10). In each panel, blue curves show I CaL predictions from models within the corresponding group, while the remaining predictions are shown in gray. even "gold-standard" model has emerged for I CaL , or even subcomponents δ and O. This situation presents a challenge for future modelers who need to choose an I CaL model for their studies, or may even wish to add a 74th model. In this section of the review we will focus on two open problems: (1) How can we select the best I CaL model for any particular study-and how can we make this task easier in the future? (2) What new experimental or analytical work is needed to help us choose between different models of the electrochemical driving force and (voltage and calcium-dependent) gating?

| Choosing I CaL models
A popular quote from statistician George Box runs "all models are wrong, but some are useful" (Box, 1979). This has occasionally been applied to models of electrophysiology, and indeed if we consider the size of an LCC alpha subunit ($170 kDa), let alone of the full macromolecular complex formed by an LCC and its various subunits and regulators, it is easy to see how a 2, 8, or even 24-state model could never hope to capture its full complexity. Instead, Box (1976) argues, we should (1) apply Occam's razor to select the simplest, most parsimonious, model that can accommodate the relevant data, and (2) choose a model that is wrong in ways that are both known and acceptable. Note that these criteria depend on the model's intended context of use (Pathmanathan et al., 2017;Whittaker et al., 2020), so that there is room for more than one "best" I CaL model in cell electrophysiology.
A good example of this philosophy is the Ohmic driving term, which is used in many models of ionic currentsdespite the common knowledge that the GHK equation is more accurate, and that this in turn is a clear oversimplification of the complex process of ion permeation (see e.g., Roux et al., 2004). For many currents, the Ohmic and GHK models both provide a good enough fit in the situations of interest and so the Ohmic model can be chosen based on its simplicity. The model's limitations (e.g. poor predictions at extreme "non-physiological" voltages or ionic concentrations) are known to the modeler and understood to be acceptable within the model's intended context of use.
Can we use these principles to choose between our 73 models? If our goal is to match a particular data set, a reasonably straightforward (but very time consuming) approach could be to simply try fitting each gating scheme and each driving term model, choosing the simplest combination that fits, and then verifying its predictive power on a separate data set (which should tell us something about its limitations). Sections 2.4 and 2.5 and the CellML code supplied with this review provide an excellent starting point for such studies.
If our goal is to pick an already parameterized model from the literature, it is much more difficult to see if these principles apply. I CaL model development is relatively well described in the literature: AP modeling articles like Luo and Rudy (1994), Jafri et al. (1998), and  devote several paragraphs and figures to it. In some cases, there is evidence of model complexity increasing when new biophysical detail was added, for example, the addition of new gates in gating schemes A-C-L or A-B-F-O. But it is interesting to see that many of the gating schemes in Figure 5 represent the same biophysical ideas of mode-switching  or the removal and application of a "brake" (Pitt et al., 2001). We also found very little evidence of existing I CaL models being tried and found inadequate before new ones were introduced (a notable exception being Luo & Rudy, 1994). It is unclear if such comparisons were not performed, for example, due to unavailability of previous model code and data sets, or if they were merely not published. Methods for model selection have been proposed, and applied in ion channel modeling as early as Horn and Vandenberg (1984). Recent progress includes enumerating, fitting, and comparing a large number of possible MM structures (Mangold et al., 2021), and such model selection techniques should be explored more widely in this field . Limitations of I CaL models are better documented, although these are usually given exclusively in terms of potential detail that could be added in future work.
Looking forward, there are several things we can do to facilitate future model comparison and selection. The field of cardiac cell electrophysiology has an unusually long history of sharing model and simulation code, dating back at least to the 1980s , and considerable effort has been expended to retroactively encode and share models in CellML (Yu et al., 2011). There is some distance to go yet-just over half the models in this review could be downloaded and used without corrections-but model code sharing on the whole appears to be getting both easier and more common.
Data sharing has a more limited history. Although the need for sharing data (and sufficient meta data) is widely knowledged (Quinn et al., 2011) and technical barriers have all but disappeared in the last decade, it is still uncommon for electrophysiological recordings to be published in full digital form. As a result, modelers have frequently resorted to digitizing published figures and then fitting to "summary curves" (e.g., I-V curves) that represent averaged data from several cells in a condensed form. Both averaging (Golowasch et al., 2002) and fitting to summary curves  can lead to inaccuracies, so there is an evident need for more widespread sharing of experimental data, perhaps in a systematic form such as the Physionet database for ECG signals (Goldberger et al., 2000). This situation will need to be drastically improved before we can answer simple questions like "which model best fits the currently known data" (for a particular species and cell type), or set up public "benchmarks" that models can be tested against, for example, using the Cardiac Electrophysiology Web Lab (Cooper et al., 2016;Daly et al., 2018).
In the articles we reviewed, there is a notable shift in focus from the model itself (with new insights gained from the model in a short section at the end of the article) to the model application and clinical/biophysical interpretation of its predictions. This can result in vital information about the model or data being moved into (less thoroughly examined) supplements, or even omitted. An interesting new approach to making models and data useful beyond a single article, is to supply a detailed description (but also comparisons to previous models and more information about model limitations) in a secondary publication, for example in The Physiome Journal for models or Scientific Data, F1000Research, or PLOS ONE for data.

| Choosing I CaL model parts
A major difficulty when comparing I CaL models is that, although they share a common structure, their driving term δ ð Þ, voltage-dependent kinetics, and calcium sensitivity interact in ways that allow irregularities in one aspect to be compensated in another. For example, the unusual open probability O ð Þ of the model by Priebe and Beuckelmann (1998) (Figure 8) was shown to be compensated by a large reversal potential in its driving term ( Figure 6) to give a typical I CaL (https://chaste.cs.ox.ac.uk/q/2022/ical/priebe). Such compensation allows the model to predict a current similar to other models for some protocols (Figure 17, middle) but different for other protocols (Figure 17, right). Similarly, the models by Demir et al. (1994), Lindblad et al. (1996), and  allow O to be greater than one (so that the usual interpretation of open probability does not apply), which can be compensated for in either δ or the maximal conductance. The model by Beeler and Reuter (1970) has a strong sensitivity to Ca 2þ Â Ã ( Figure 11) but does not include CDI, deriving its Ca 2þ Â Ã -dependence from an unusually strong driving force instead ( Figure 6). Compensation also comes into play when creating or modifying AP models, which often require calibration of the sub-models (e.g. the I CaL model) to cell-level observations . For example, Figure 12 shows how a high sensitivity to calcium can be compensated by postulating a bigger subspace near the LCCs, leading to a lower local Ca 2þ Â Ã : : More generally, an I CaL model's CDI gating may have to be tuned to the specific CaT in the AP model it was designed for (which suggests that some of the variability observed in Section 3.6 may be attributed to compensation for differences in CaT). Similarly, irregularities in an I CaL model can, to a degree, be compensated by irregularities in the outward K þ currents and in the other components involved in CICR.
A model with subtly compensated defects may reproduce many electrophysiological phenomena seen in the experimental data, but through mechanisms that do not reflect those in the real cell. As a result, its predictions are likely to be unreliable when extrapolating to new (perhaps pathophysiological) situations. In this section, we once again break down an I CaL model into its constituent parts, and discuss possible experiments that could let us study each in isolation.

| Driving term
The study of the electrochemical force driving I CaL flux through the open channels is complicated by the wide range of calcium concentrations in a myocyte and the much higher concentrations of Na þ and K þ ions near the LCCs. Different models have been proposed as a result, starting with the Ohmic and GHK forms discussed in this review, but quickly moving on to more elaborate equations Roux et al., 2004) that (perhaps due to their computational cost) were never incorporated into AP models. This has left the question of how best to approximate δ in a computationally tractable I CaL model unanswered, and none of the studies we surveyed showed a comparison of different drivingterm models to motivate their choice.
To decide between driving force models, specialized experiments may be required that can separate driving force from kinetics and CDI. Single-channel "excised-patch" experiments have shown promise in this respect (Greenstein & Winslow, 2002;Hess et al., 1986), and have the potential to allow rapid variation of internal (inside-out) or external calcium concentrations (outside-out), to levels that a whole cell may not tolerate, and in removing obfuscating effects from intracellular buffers and transporters (Gradogna et al., 2009). Single-channel results could then be scaled up to wholecell level, with stochastic channel simulations used to check if a simple scaling approximation holds well enough when only small numbers of channels are open, or if a more complex method is needed.
Another interesting approach could be to fit Ohmic and/or GHK models to predictions from more sophisticated models of ionic permeation (Roux et al., 2004).
The same approaches could be used to try and resolve differences in relative permeability to different ionic species (Table 4). In the models reviewed in this study, we found that the activity coefficients were highly similar. It is interesting to note that their values depend on experimental conditions and so they could be calculated as part of the model (although within physiological ranges of temperature and ionic concentrations they change by less than 0.1% during an AP at 1 Hz pacing).

| CDI and calcium dynamics
CDI is one of the most challenging aspects of modeling I CaL and the degree of CDI predicted by different models varies widely ( Figure 11). CDI depends on the (fluctuating) calcium concentration at the intracellular channel side, and so model assumptions about channel localization and cell geometry have a strong impact on modeled CDI, as evidenced by the relationship between subspace volume and CDI 50 in Figure 12 and by the reduction in variability between models when grouping by channel localization.
Experimentally, CDI can be separated from voltage-dependent kinetics and driving term to some degree by using voltage protocols with a single repeating step, so that voltage effects are consistent-if still unknown. However, the local calcium concentration is difficult to control, as the very low Ca 2þ Â Ã maintained in living cells means the calcium influx through I CaL causes a significant change. In myocytes, a bewildering array of processes affect and are affected by local Ca 2þ Â Ã (Bers, 2008). In heterologous expression systems, calcium-buffers are required to keep Ca 2þ Â Ã at levels the cells can tolerate (You et al., 1997). Bulk cytosolic Ca 2þ Â Ã can be measured optically using fluorescent (Herron et al., 2012) but measuring local calcium elevations in subspaces (or "microdomains") is technically challenging (Acsai et al., 2011). As a result, there is no simple set-up in which CDI can be studied, and analyzing data from experiments targeting CDI requires consideration of intracellular Ca 2þ Â Ã gradients, diffusion, and buffering. Although analytical Peak current (normalized) Current (normalized) Current (normalized) F I G U R E 1 7 Compensation can hide model differences. Left: The I-V curves for Inada et al. (2009) and Priebe and Beuckelmann (1998), normalized to the peak I CaL . These models were chosen because their I-V curves almost overlap, suggesting similar activation characteristics. However, the model by Priebe and Beuckelmann is known to have an unusually low open probability at higher voltages ( Figure 8) which is compensated by a stronger driving term ( Figure 6). Middle: A step protocol showing quantitative, but not qualitative model differences. The membrane potential is held at À90 mV (not shown), then stepped up to 10 mV (0-150 ms), and then down to 5 mV. Currents are shown using the same normalization as before. In both models, the step to 10 mV elicits the maximal response, followed by inactivation, which continues after the step to 5 mV. Right: A step to 60 mV elicits a much weaker response in both models (notice the scale on the y-axes). At this potential, we expect the models to fully activate and inactivate. However, the subsequent step to 5 mV reveals further activation in the model by Priebe and Beuckelmann, leading to an unexpected biphasic response, and a qualitative difference between the models' activation characteristics.

| Voltage-dependent kinetics
The equations for O vary more from model to model than any other aspect of I CaL (Figure 5), and we found large variability in the predicted response to voltage-clamp protocols, which cannot be easily reduced by grouping models by cell type or other factors. While most studies of I CaL kinetics have used conventional voltage-step protocols (such as shown in Figures 7, 9, and 10), recent advances in protocol design for I Kr Clerx et al., 2019) may be translatable to I CaL . Application of such protocols in high-throughput settings (Lei et al., 2019) with careful consideration of experimental artifacts Lei, Fabbri, et al., 2020) may help produce rich data sets that will allow us to distinguish between the different proposed gating mechanisms and study (biological or experimental) variability. Such an approach, however, would require that voltage and calcium effects are independent (or at least can be studied independently) and may rely on accurate CDI and driving term models already being available. Alternatively, it may be possible to derive gating models as simplifications of more detailed molecular dynamics models (Ramasubramanian & Rudy, 2018;Silva, 2018). In either case, care must be taken to create independent training and validation data sets, so that the most crucial part of the gating model, its predictive power, can be assessed .

| Determining the impact of I CaL model choice on electrophysiological properties
I CaL is an important determinant of key electrophysiological properties at the cell, tissue, and organ level, including APD restitution (Guzman et al., 2010;, development of cellular, discordant, and T-wave alternans (Zhu & Clancy, 2007), EAD formation (Heijman et al., 2014;Kettlewell et al., 2019;Weiss et al., 2010), and cell-to-cell conduction (Joyner et al., 1996;Rohr & Kucera, 1997;Shaw & Rudy, 1997). It is therefore of interest to compare how the 73 models described in this study affect such emergent properties. However, a comprehensive study along these lines requires not just an isolated I CaL model, but rather an I CaL model in the context of a full AP model (or tissue model).
An obvious way to start such a study would be to pick an AP model, replace its I CaL formulation and test its predictions (repeating the process for all or a subset of the models in this review). This is a nontrivial exercise for several reasons, including: (1) The shape of the AP could be significantly affected, even after tuning the I CaL conductance. (2) The CaT and detailed calcium dynamics could be similarly affected, and it may not be possible to set a new maximum conductance that simultaneously corrects the AP shape and CaT. (3) The I CaL model may be incompatible with the CaT and AP that the cell model provides (e.g. it may not activate or inactivate to the right degree at the voltages and Ca 2þ Â Ã that the AP model predicts). (4) Different assumptions about LCC localization would need to be reconciled. (5) The new I CaL model could predict different Na þ and K þ fluxes, which could affect the model's long-term stability.
One way in which a preliminary study can be set up is by only changing the O from one I CaL model to the other, and subsequently tuning the g Ca =P Ca (along with P Na =P K if applicable) so that the resultant AP and/or CaT most closely resemble their original shapes (thereby treating obstacles 1 and 2). Assuming this can be done, this leaves obstacles 3-5, but these might be avoided if a subset of "similar enough" (but interestingly different) I CaL models can be found. By changing the I CaL kinetics through O, a study could evaluate for example, how changes in the time constant of recovery affect restitution  or how changes in a model's I CaL "window-current" affect the onset of EADs (Kettlewell et al., 2019). Note however, that such properties are not always directly encoded in a single model parameter, so that more complex multi-parameter adjustments may need to be made.
A final important complication arises since emergent cell, tissue, and organ properties depend on more than a single "molecular factor" (Weiss et al., 2015); nonlinear interactions of I CaL with other currents that are also modeled differently in various AP models will influence emergent behavior. For instance, an I CaL model "correctly" producing alternans in one AP model may fail to do so in another (Hopenfeld, 2006). So any study inserting different I CaL models into an AP model would also need to consider the influence of the AP model choice.

| CONCLUSION
Five decades of careful research has resulted in an impressive list of 73 mammalian I CaL models, 15 models of LCC gating, and a smaller but similar lack of consensus on how to model I caL 's electrochemical driving force. Accordingly, the 60 models studied quantitatively made very different predictions about activation, inactivation, recovery, and CDI. However, as practical limitations on model sharing, data sharing, and extended publication forms are rapidly disappearing, the field is poised to take advantage of new technological possibilities (and publishing norms) that might allow these decades-old questions to be systematically investigated and resolved. We believe that this review is the first to systematically survey and compare I CaL models and hope that it serves as a starting point for a critical re-assessment of L-type calcium current modeling from which a synthesized, community consensus model may emerge.