Application of stochastic phenomenological modelling to cell-to-cell and beat-to-beat electrophysiological variability in cardiac tissue.

Variability in the action potential of isolated myocytes and tissue samples is observed in experimental studies. Variability is manifested as both differences in the action potential (AP) morphology between cells (extrinsic variability), and also 'intrinsic' or beat-to-beat variability of repolarization (BVR) in the AP duration of each cell. We studied the relative contributions of experimentally recorded intrinsic and extrinsic variability to dispersion of repolarization in tissue. We developed four cell-specific parameterizations of a phenomenological stochastic differential equation AP model exhibiting intrinsic variability using APs recorded from isolated guinea pig ventricular myocytes exhibiting BVR. We performed simulations in tissue using the four different model parameterizations in the presence and the absence of both intrinsic and extrinsic variability. We altered the coupling of the tissue to determine how inter-cellular coupling affected the dispersion of the AP duration in tissue. Both intrinsic and extrinsic variability were gradually revealed by reduction of tissue coupling. However, the recorded extrinsic variability between individual myocytes produced a greater degree of dispersion in repolarization in tissue than the intrinsic variability of each myocyte.

Introduction

Dispersion of repolarization is defined as a spatial heterogeneity in the time of repolarization between different regions of cardiac tissue. Dispersion of repolarization is associated with arrhythmic risk (Chauhan et al., 2006). A potential source of dispersion of repolarization is variability in action potential (AP) morphology and duration in the myocytes within cardiac tissue. Variability in the AP duration (APD) of cardiac myocytes may arise from naturally occurring differences in ion channel density between cells, the spatial location of the cell within the heart, and circadian rhythms (Jeyaraj et al., 2012). Cell isolation procedures may also affect the AP of cells that have been isolated from the myocardium (Yue et al., 1996). We describe these differences in AP morphology between different myocytes as ‘extrinsic variability’.

An alternative source of variability in APD is ‘intrinsic’ beat-to-beat variability of repolarization (BVR) in the same cell. Beat-to-beat variability in the APD of isolated cells (Zaniboni et al., 2000), beat-to-beat dispersion of repolarization in tissue (Hondeghem et al., 2001), and variability in the QT-interval at the ECG level (Thomsen et al., 2004) have all been reported. At a cellular level, BVR is an apparently random variation in the APD (Johnson et al., 2010). This is distinct from alternans, which is a regular short-long oscillation in the APD between subsequent beats. BVR at both a cellular and a whole-heart level has been shown to be a predictor of arrhythmogenic risk, although the mechanisms by which this occurs are unknown (Kääb et al., 2003, Johnson et al., 2010).

Extrinsic variability between different cell and tissue samples is difficult to reproduce experimentally, as it is specific to the sample under investigation. Computational modelling studies may provide a method to understand how this extrinsic variability contributes to dispersion of repolarization in the ventricles. Many simulation studies have focussed on heterogeneities that exist within the ventricle; for example transmurally, or in an apex-base direction (Efimov et al., 1996, Viswanathan et al., 1999, Keller et al., 2012). The first approaches to modelling consequences of variability in ion channel expression in single cells are now beginning to appear (Sarkar and Sobie, 2011, Walmsley et al., 2013, Britton et al., 2013, Gemmell et al., 2014).

Intrinsic beat-to-beat variability at a cellular level has also been suggested as contributing to dispersion of repolarization in tissue. Simulation studies by Pueyo et al. (2011), Lemay et al. (2011), and Heijman et al. (2013) demonstrated that the BVR in isolated cells induced by stochasticity in ion channel gating is dramatically reduced when cells are electrotonically coupled in tissue. These studies all showed that BVR can re-appear as the conductivity of the tissue decreases.

In this study, we investigate how dispersion of repolarization in the ventricular myocardium is affected by the relative contributions of intrinsic beat-to-beat variability, and extrinsic variability between different cells. We also investigate how these contributions are influenced by inter-cellular coupling, as motivated by the findings of Lesh et al. (1989), Lemay et al. (2011) and Pueyo et al. (2011), who observed that the degree of coupling between cells affected dispersion of repolarization.

Our simulations of intrinsic and extrinsic variability are based on experimental recordings from isolated guinea pig ventricular myocytes that exhibit both intrinsic and extrinsic variability. We use a phenomenological model that reproduces both cell-to-cell (extrinsic) and beat-to-beat (intrinsic) variability (Bueno-Orovio et al., 2008, Walmsley et al., 2010). We use four parameter sets that reproduce four specific isolated cell samples so that our simulations are based on physiological levels of single cell variability (Walmsley et al., 2010). We then couple these cells into a simulated tissue using a monodomain approach.

Methods

Reproducing intrinsic beat-to-beat variability

Action potential model

Our phenomenological model of beat-to-beat intrinsic variability is based on the Bueno-Orovio–Cherry–Fenton (BOCF) phenomenological model, described in detail in Appendix A (Bueno-Orovio et al., 2008). We modified the BOCF AP model to incorporate intrinsic beat-to-beat variability as described by Walmsley et al. (2010). A Weiner process is a stochastic process whose increment dW on the interval has mean 0 and variance . We modified Eqs. (A.2), (A.3), (A.4) by including a Weiner increment in each, forming a system of stochastic differential equations (SDEs):where σ, σ, and σ are constants changing the magnitude of the Weiner increment (see, for example, Kloeden and Platen, 2011). The superscript denotes an independent Weiner process for each variable.

Reproduction of experimentally observed variability

Examples of APs simulated using the parameter sets for the SDE model derived by Walmsley et al. (2010) for each of the four cells are shown in Fig. 1A (bottom row). These simulations show good qualitative and quantitative agreement with both the experimentally recorded action potential morphologies, and spread of repolarization times shown in Fig. 1A (top row). There was no variation in the resting potential in the simulations. Fitting of the noise terms was successful in matching the mean and variance of the experimental data (Table 1). APD in the simulations and experiments was quantified by the time from upstroke to 90% repolarization (APD90). Poincaré plots demonstrate the degree of temporal variability by plotting the APD90 of each beat against the APD90 of the preceding beat. Fig. 1B shows Poincaré plots of the temporal variability in APD for apex cell 1 from the experimental data and simulations. The fitted noise terms from Walmsley et al. (2010) are shown in Table 1. For each parameter set we therefore have the stochastic version as shown in Fig. 1, and also a deterministic version where each noise term is set to zero.

Fig. 1

Experimental and simulated temporal variability in repolarization. (A) Experimental data showing repeated 1 Hz stimulations of isolated guinea pig ventricular myocytes from the apex and the base of the heart. Simulated reproductions of the data using four parameterizations of the SDE model are shown below each set of experimental data. Modified from Walmsley et al. (2010). (B) Poincaré plots of the original experimental data (red) and those generated from the guinea pig apex cell 1 parameter set are shown at pacing cycle lengths of 1000 ms, 400 ms, 300 ms, 260 ms, and 220 ms for the stochastic (black) and deterministic (green dot) versions of the model. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Table 1

Experimental and simulated variability in APD90.

Parameter	Apex 1	Apex 2	Base 1	Base 2
Mean (exp) (ms)	219.5	133.7	286.5	185.8
Variance (exp)	37.69	16.80	75.56	21.88
σv	0	0	0	0
σw	0.0012	0.0015	0.0011	0.0011
σs	0.0012	0.002	0.0013	0.0012
Mean (sim) (ms)	223.5	127.3	268.7	182.3
Variance (sim)	39.22	17.64	77.97	23.73

Two-dimensional tissue simulations

Modelling conduction in tissue

All simulations were performed in a two-dimensional tissue sheet of dimension 1 cm ×1 cm. Inter-cellular coupling was simulated using a monodomain approach, where Eq. (A.1) was spatially extended toWe applied zero-flux Neumann boundary conditions to the tissue. In all simulations, the stimulus was applied at the left-hand edge of the tissue ().

We used four different values of the diffusion coefficient D to investigate the effects of reduced coupling upon temporal variability of repolarization. In the original BOCF model D=1.171 cm2 s−1 was used for tissue simulations, based on experimental measurements of human ventricular cell dimensions and cytoplasmic resistivity (Bueno-Orovio et al., 2008), giving a conduction velocity of 72 cm s−1. Simulations were also performed at 50%, 10%, and 5% of this value (D=0.586 cm2 s−1, 0.117 cm2 s−1, and 0.059 cm2 s−1, respectively), giving conduction velocities of 51 cm s−1, 23 cm s−1 and 16 cm s−1, respectively.

Determining the effect of intrinsic variability alone

To investigate the effect of intrinsic BVR alone without extrinsic variability between cells, we created a sheet of tissue with a homogeneous AP model. We repeated this for each of the four parameter sets in Section 2.1, giving four different sheets of tissue. 50 stimuli were applied at a pacing cycle length of 1000 ms. We repeated the simulations for both the deterministic version and the stochastic versions of the AP model to observe the effects of intrinsic variability alone. In order to ensure that simulations using the stochastic AP models were comparable, the same pseudo-random number sequence was used for each stochastic simulation. A different sequence was generated for each node in the computational mesh.

Effect of fast pacing rates on dispersion of repolarization

We also investigated the consequences of the alternans-like behaviour induced in the stochastic AP model at short cycle lengths shown in Fig. 1B. In order to determine whether this alternans affected the magnitude of dispersion of APD measured in the intrinsic beat-to-beat variability case we used apex cell 1. The simulations had an initial cycle length of 230 ms, where alternans-like behaviour began to appear in the stochastic cell model using the apex cell 1 parameter set (Fig. 1B). The cycle length was reduced in 2 ms intervals after every 40 beats in the simulation until conduction block occurred. Conduction block was defined as a failure of the stimulus to lead to excitation of the entire tissue. Simulations were performed for both the deterministic and the stochastic models to determine the degree of BVR in tissue due to intrinsic variability.

Comparing the relative effects of intrinsic and extrinsic variability

In order to simulate the effect of extrinsic inter-cell variability, one of the parameter sets corresponding to the four cells in Section 2.1 was applied to the cell model at each node in the computational mesh. The parameter sets were applied randomly with equal probability. Simulations were performed for both the deterministic (extrinsic variability only) and the stochastic version (extrinsic and intrinsic variability) of the resulting tissue model. 50 stimuli were applied at a pacing cycle length of 1000 ms. The same spatial distribution of cell parameter sets was used for both the deterministic and the stochastic simulation to enable direct comparison between results.

Data analysis

Simulation results were analysed by recording APD90 after each simulation. We did not perform these analyses where or to prevent electrotonic effects arising from the Neumann boundary conditions or the stimulus site influencing the APD90 arising from the cell model (Cherry and Fenton, 2011). The maximum dispersion of APD was measured as the maximum difference in APD90 across the region 0.25 cm at each node in the mesh. The difference in APD90 between the stochastic and the deterministic simulation was also computed for every stimulus. Furthermore, when considering BVR across multiple beats in the tissue, we calculated the mean APD90 in the tissue during each activation.

Numerical implementation and simulation software

Simulations in tissue were performed using the finite element method, as implemented in the open source cardiac simulation package Chaste (Mirams et al., 2013). The deterministic or the stochastic cell model was applied to each node of the tissue mesh as appropriate. The solution of the deterministic ODE at each node in the mesh was determined using a forward Euler scheme, and the stochastic SDE models were solved using an Euler–Marayama scheme (Kloeden and Platen, 2011). Independent Weiner processes were generated for each node in the tissue and each stochastic variable at that node. Numerical convergence for both the deterministic and stochastic schemes was established and a spatial discretization of 0.01 cm, PDE time step of 0.01 ms, and an ODE/SDE time step of 0.00125 ms were used in this study, as described in Appendix B. A user project containing all of the code used in this study is available to download from the Chaste website (http://www.cs.ox.ac.uk/chaste).

Results

Intrinsic variability alone has a small effect on AP dispersion

Tissue simulations of intrinsic variability only showed temporal variability in the mean repolarization time for each of the four parameter sets, shown in Fig. 2. Variability in mean APD90 between beats increased slightly as the diffusion coefficient D was decreased. This was quantified by calculating the variance in mean APD90 across all 50 beats (see Table 2). The distribution of mean APD90 for the stochastic simulations with intrinsic variability was centred on the mean APD90 of the deterministic simulations (no intrinsic variability) in all cases. The beat-to-beat variance in mean APD90 also increased with the variance in APD90 of the underlying cell (Table 2, see brackets after cell name for isolated cell variance in APD90).

Fig. 2

Reducing coupling in tissue simulations with intrinsic variability increases dispersion in APD90. (A) The mean APD90 across all nodes in the central 0.5 cm of the tissue is plotted for both the deterministic (left) and stochastic (right) simulations for each parameter set and each value of D. The values of D shown are 1.171 cm2 s−1 (red), D=0.586 cm2 s−1 (green), D= 0.117 cm2 s−1 (blue), and D= 0.059 cm2 s−1 (black). (B) The difference between the stochastic and deterministic APD90 maps is calculated at each node and each AP in the central 0.5 cm of the tissue. The variance across these nodes is then calculated for each beat. The variance in the difference between the deterministic and stochastic APD90 maps is plotted as a function of the number of beats. (C) The absolute difference between the deterministic and stochastic simulations is calculated at each node in the central 0.5 cm of the tissue. The maximum of these absolute differences for each beat are then plotted as a function of the number of beats. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Table 2

Summary of intrinsic variability only simulations : the mean APD90 in the simulated tissue over one beat. Mean(): average of over all beats in the simulation. Var(): variance in over all beats in the simulation. Det.: deterministic simulations. Stoch.: stochastic simulations. Variance in brackets below each cell is from the isolated cell simulations shown in Fig. 1, shown for comparison.

Cell		D (cm2 s−1)
		1.171	0.586	0.117	0.059
Apex cell 1	Det. Mean(μ_APD) (ms)	219.8	219.8	219.7	219.7
(var. 39.22 ms2)	Stoch. Mean(μ_APD) (ms)	219.8	219.8	219.7	219.7
	Stoch. Var(μ_APD) (ms2)	0.020	0.023	0.033	0.037
Apex cell 2	Det. Mean(μ_APD) (ms)	134.3	134.3	134.3	134.3
(var. 7.64 ms2)	Stoch. Mean(μ_APD) (ms)	134.3	134.3	134.3	134.3
	Stoch. Var(μ_APD) (ms2)	0.012	0.014	0.020	0.021
Base cell 1	Det. Mean(μ_APD) (ms)	288.4	288.4	288.3	288.3
(var. 77.97 ms2)	Stoch. Mean(μ_APD) (ms)	288.4	288.4	288.4	288.3
	Stoch. Var(μ_APD) (ms2)	0.041	0.047	0.065	0.067
Base cell 2	Det. Mean(μ_APD) (ms)	190.8	190.8	190.8	190.8
(var. 23.73 ms2)	Stoch. Mean(μ_APD) (ms)	190.8	190.8	190.8	190.8
	Stoch. Var(μ_APD) (ms2)	0.013	0.016	0.022	0.025

Spatial dispersion in the APD90 resulting from intrinsic variability, as measured by the variance in the difference between the stochastic and the deterministic simulation for each beat is shown in Fig. 2B. Fig. 2B shows that the variance increased as D decreased for each parameter set, demonstrating an increase in dispersion of APD90 resulting from intrinsic variability.

The difference between the deterministic and stochastic simulations did not increase with the number of beats, as shown in Fig. 2C, showing that intrinsic variability does not cause increased tissue level BVR over time in tissue. The maximum difference between the stochastic and deterministic simulations increased as the diffusion coefficient decreased. The largest difference observed was in base cell 1 with D= 0.059 cm2 s −1, with a difference of 2.20 ms from the deterministic simulation.

Heat maps of the APD in the stochastic case for base cell 1 are shown in Fig. 3, demonstrating the increase in dispersion in APD resulting from intrinsic temporal variability in the cell model at low conductivities. A slight dispersion of APD arising from boundary effects can also be seen at higher conductivities, explained by the dependence of the size of the boundary effect with the conductivity of the tissue (Cherry and Fenton, 2011). This was also present in the deterministic simulations. The dispersion created a variance in both the deterministic and stochastic models at D= 1.171 cm2 s−1 of up to 0.18 ms2 in the deterministic case and 0.4 ms2 in the stochastic case (base cell 1).

Fig. 3

Heat maps of spatial dispersion in APD90 resulting from intrinsic variability for base cell 1. APD90 maps shown are recorded from base cell 1 with (A) D = 1.171 cm2 s−1, (B) D = 0.586 cm2 s−1, (C) D = 0.117 cm2 s−1, and (D) D = 0.059 cm2 s−1. The heat maps shown are generated using the stochastic cell model and are taken from the final beat in the series. APD90 maps are plotted from the central 0.5 cm to reduce distortion due to boundary effects.

Intrinsic variability can produce tissue level alternans at fast pacing rates

We ran simulations starting at a cycle length of 230 ms, and reduced the cycle length by 2 ms after every 40 beats until a conduction block occurred at a pacing cycle length of 222 ms. Fig. 4A shows a plot of mean APD90 recorded at all four values of D over 160 beats. Mean APD90 can be seen to oscillate dramatically in response to a change in cycle length, but oscillations are damped over time in both the deterministic and the stochastic cases, leaving a small beat-to-beat alternans in mean APD90.

Fig. 4

Summary statistics on maximum APD dispersion for the fast pacing protocol. Results are shown for D=1.171 cm2 s−1 (red), D= 0.586 cm2 s−1 (green), D= 0.117 cm2 s−1 (blue), and D= 0.059 cm2 s−1 (black). Vertical dashed lines denote a change in the pacing cycle length. (A) The mean APD90 is calculated across the nodes in the central 0.5 cm in the tissue at each beat in the deterministic (left) and stochastic (right) models. The variance of the APD is calculated over the nodes in the central 0.5 cm of the tissue for each beat. (B) The variance in APD90 is plotted for the deterministic and stochastic models. The insets show an enlarged view of beats 80–120. (C) The absolute difference in APD90 between the deterministic and stochastic simulations is calculated at each node in the central 0.5 cm of the tissue. The maximum absolute difference at each beat is shown. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

We compute the variance within the APD90 map at each beat in order to evaluate the degree of dispersion in the APD, as shown in Fig. 4B. The dispersion is comparable between the stochastic and deterministic simulations following a change in pacing cycle length. In fact, the stochastic model displayed a reduced change in the mean value in the beats immediately following the change in cycle length as compared to the deterministic model, as is also seen in the isolated cell simulation in Fig. 1B. The transient increase in variance was due to the change in cycle length inducing a global alternans in both models that was modulated by boundary effects. To remove this effect, we also examined the mean and variance of the APD once the oscillations have stabilized in the final ten beats of the AP at each cycle length. For these beats, the mean value of the APD shows a moderately increased variability at low D as compared to the deterministic model. However, there is an increase in variance in the stochastic model, but not in the deterministic model, as the diffusion coefficient of the tissue is reduced, showing that the dispersion of the APD90 increases when including intrinsic variability at low values of D.

The overall effect of intrinsic variability on the simulation can be measured by calculating the maximum difference between the stochastic and deterministic simulations, as shown in Fig. 4C, for all values of D. The difference between the stochastic and deterministic simulations increases as D is reduced. The maximum difference for the lowest value of D (0.059 cm2 s−1) is bounded by ±1.5 ms, again smaller than the APD range in the stochastic version of apex cell 1 (31.4 ms).

In the simulations with intrinsic variability, we observed phenomena that resembled a spatially discordant beat-to-beat alternans. Fig. 5 shows a sequence of six APD maps in the stochastic model from the end of the 228 ms pacing cycle length section of the pacing protocol at the lowest value of D (0.059 cm2 s−1). In the deterministic simulations, there is almost no spatial heterogeneity in the APD90, and similarly with D= 1.171 cm2 s−1 almost no dispersion is observed in the stochastic model. At D= 0.059 cm2 s−1, the stochastic simulations show a modest spatially heterogeneous alternation in the simulated APD90 that is not present in either the equivalent deterministic recording at the same value of D, or in either the stochastic or deterministic recording for D=1.171 cm2 s−1. The maximum amplitude of the oscillations is approximately 2 ms.

Fig. 5

Spatially discordant alternans-like behaviour recorded in the action potential duration at 228 ms pacing cycle length in the stochastic model. The APD90 maps shown are simulated at a cycle length of 228 ms (beats 75–80) for the stochastic model for D= 1.171 cm2 s−1 and 0.059 cm2 s−1. APD90 maps are plotted from the central 0.5 cm only to remove boundary effects.

Extrinsic variability between cells has a greater effect on dispersion of repolarization than intrinsic variability in APD

As in the intrinsic beat-to-beat variability-only simulations, the mean APD90 of the deterministic simulations was constant over all beats when including extrinsic variability between cells. There was a beat-to-beat variability in mean APD90 in the stochastic simulations that also included beat-to-beat intrinsic variability in the cells (Fig. 6A). Thus, given extrinsic variability in tissue generated from our four cells using a uniform random variable at each node, including intrinsic temporal variability, induces beat-to-beat variability in mean APD90 behaviour in tissue.

Fig. 6

Mean and variance in APD90 from tissue simulations incorporating extrinsic variability. (A) The mean of the APD90 is calculated over all nodes in the central 0.5 cm of the tissue at each beat and plotted for both the deterministic simulations (left) and stochastic simulations (right) for each beat. (B) The variance in the APD90 for both the stochastic simulations (solid lines) and deterministic simulations (dashed lines) for each beat. (C) The maximum difference in APD between the stochastic and deterministic simulations at each beat, for each simulation. Values of D are 1.171 cm2 s−1 (red), D=0.586 cm2 s−1 (green), D= 0.117 cm2 s−1 (blue), and D= 0.059 cm2 s−1 (black). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

The variability in APD90 grows as D decreases due to the unmasking of extrinsic variability in the AP between cells, as shown by the increase in variance in the APD90 shown in Fig. 6B. The variance in the APD90 for the simulations of extrinsic variability between cells-only (deterministic) was slightly lower than when both extrinsic variability between cells and intrinsic variability from beat-to-beat in each cell were simulated (stochastic) as shown in Fig. 6B. However, there is still a much larger increase in the variance in the extrinsic variability-only case than is observed for any of the intrinsic variability only cases from Section 3.1. At the lowest value of D, the variance reaches a maximum of 6.1 ms2, which remains smaller than each simulated cell׳s variance in APD, which ranges between 17.64 ms2 for apex cell 2 and 77.97 ms2 for base cell 1 (Table 2). The difference between the deterministic and stochastic simulations remained at similar values to those seen for the intrinsic variability only simulations, being bounded between −1 ms and +2 ms at the lowest value of D (0.059 cm2 s−1). This is again dramatically smaller than the dispersion in the isolated cell simulations, which ranged between 22.3 ms for apex cell 2 and 47.1 ms for base cell 1.

Heat maps of both deterministic and stochastic simulations from the final beat are shown in Fig. 7B, together with the underlying cell distribution in Fig. 7A. This figure demonstrates the increase in dispersion of repolarization as D is reduced. The maximum dispersion observed is 11 ms with D= 0.059 cm2 s−1. As the differences arising from the extrinsic variability between the simulated cells are much larger than those resulting from intrinsic temporal variability in the simulated cells, the heat maps in Fig. 7 are almost identical.

Fig. 7

Distribution of APD90 from tissue simulations incorporating extrinsic variability. APD90 maps shown are recorded are for (A) the underlying cell distribution, then in (B) at D= 1.171 cm2 s−1, 0.586 cm2 s−1, 0.117 cm2 s−1, and 0.059 cm2 s−1 for the deterministic and stochastic simulation results. Simulation results shown are for the 50th beat in the series, and from the central 0.5 cm of the tissue only.

Discussion

In this study we have applied cell-specific parameterizations of a phenomenological AP model representing intrinsic temporal variability and extrinsic inter-cell variability to investigate consequences of variability in tissue. We have previously shown that a simple phenomenological model can reproduce cellular recordings of BVR (Walmsley et al., 2010). Using this model, we have now investigated the extent to which the intrinsic cellular APD impacts tissue level dispersion of APD, as well as the role of extrinsic cell–cell variability in APD. We have further investigated how this spatial dispersion in APD is affected by the coupling of the tissue.

In tissue simulations with intrinsic temporal variability only, we observe spatial dispersion of APD resulting from temporal variability in the APD of each of the four cells investigated. The dispersion of APD increases as the diffusion coefficient D is decreased, in accordance with previous studies. The dispersion resulting from temporal variability in the APD alone is relatively low. In absolute values, the maximum difference between simulations with and without intrinsic variability observed was around 2 ms in the case of base cell 1. We observed similar results at fast pacing rates, with the alternans observed both with and without intrinsic variability being damped in tissue. Intrinsic temporal variability in the APD produced a similarly low contribution to dispersion in the APD at fast pacing rates. BVR as observed in tissue is therefore unlikely to arise from cellular BVR alone. This does not exclude a role for ectopic activity such as early after depolarizations (Sato et al., 2009), or delayed after depolarizations (Johnson et al., 2013). We did not observe either of these phenomena in the four cells we investigated.

We additionally investigated whether the alternans-like behaviour observed in an intrinsic variability-only cell model at short pacing rates could cause increased dispersion of repolarization. Some simulations resemble a spatially discordant alternans in the simulations, which has been shown in animal models to be arrhythmogenic (Myles et al., 2008). The amplitude of this oscillation is small (2 ms) The amplitude of both the beat-to-beat and the spatially discordant alternans was dramatically lower than those in the single cell simulations (see Fig. 1B). One explanation for this phenomenon is that the stochastic nature of beat-to-beat variability in the model has resulted in some cells entering alternans with a long APD, and others entering alternans with a short APD, so that the overall behaviour in different regions is out of phase, but the difference is heavily damped by diffusion.

When including extrinsic variability between cells, we observed a spatial dispersion in APD that was visible even at the highest value of D (1.171 cm2 s−1). At the lowest value of D (0.059 cm2 s−1) this dispersion reached a value of 11 ms (Fig. 7), an order of magnitude larger than that observed in the temporal variability simulations in the intrinsic variability only cases (Fig. 3). This increase in dispersion resulting from inter-cell extrinsic variability alone can also be measured quantitatively by observing the large increase in variance in the extrinsic variability-only simulations as compared to the change in variance observed in the intrinsic variability-only simulations. Including intrinsic temporal variability in the extrinsic variability simulations resulted in a relatively small further increase in the dispersion as measured using variance. We conclude that the extrinsic cell-to-cell variability between the four cells we investigated provided a greater contribution to dispersion of repolarization in tissue than their intrinsic beat-to-beat variability.

Note that the degree of dispersion resulting from intrinsic beat-to-beat variability alone is too small to be detected even by microelectrode recordings. The dispersion observed in the signal when extrinsic variability is included is large enough that it could be detectable using microelectrode arrays. With a more spatially organized cell distribution, the dispersion, and hence its detectability, would increase; as seen in previous studies (Viswanathan et al., 1999).

When attempting to parameterize a model based on the AP morphology recorded from an individual cell, information on the conductances of all currents making up the observed AP are not available. Even if such information were available, ionic currents are modified by many factors, rendering them a ‘moving target’ (Roden, 2008, Carusi et al., 2012). Zaniboni et al. (2010) have shown that, due to poor parameter identifiability from the AP morphology alone, multiple behaviours in tissue are possible from apparently identical biophysically detailed models. Simplified phenomenological models such as the BOCF model therefore offer an alternative method for reproducing inter-cell variability based upon reproducing the APs recorded from individual cells (Fenton and Karma, 1998, Mitchell and Schaeffer, 2003, Bueno-Orovio et al., 2008, Bueno-Orovio et al., 2012). This method may avoid some of the risks associated with over-fitting biophysically detailed models to an individual AP, whilst still allowing investigation of the consequences of the observed behaviour. An alternative approach would be to use a calibrated population of biophysically detailed models based upon AP recordings from a large number of cells, as proposed by Britton et al. (2013), and then using the resulting parameter distributions to parameterize a tissue model containing extrinsic variability between cells.

Variability in resting potential was observed in the experiments (Fig. 1). As our work was focussed on intrinsic beat-to-beat variability, and the variability in the resting potential of these four cells was not correlated with their action potential durations, we did not include this effect. Future work using a modified phenomenological model, or a biophysically detailed model could consider the impact of this effect on dispersion of repolarization in tissue.

Limitations

The consequences of some cellular-level phenomena such as raised resting potential and accommodation of the AP in response to a change in pacing frequency cannot be investigated using the BOCF model. However, the more recent Bueno-Orovio et al. (2012) model does allow the inclusion of APD accommodation. In this study, our primary objective was to investigate the consequences of BVR in tissue. Therefore, we did not take into account mechanisms such as stochastic ion channel gating that may explain the source of fluctuations in membrane potential that lead to BVR. Exploring how stochastic gating phenomena contribute to BVR requires a biophysically detailed model for which a Markov formulation could be used within the framework of the Langevin equations (Dangerfield et al., 2012). We instead used a constant additive noise term within a phenomenological model as this approach was capable of reproducing the observed behaviour in a simple manner.

The nodes of a finite element mesh do not represent cells as such, but represent the homogenized behaviour of a small volume of cells. Thus, the level of extrinsic variability at each node will not correspond exactly to that in the isolated cells used for parameterization. Our intrinsic variability data arise from the number of ion channels (or, assuming uniform ion channel density, the membrane surface area) of a single cell. This may differ from the number of ion channels/surface area of membrane in the tissue volume represented by the AP model at a single mesh node, and so a future refinement of the variability to match mesh element size may be in order. At very low levels of cellular coupling, it may be necessary to use a discrete modelling approach as assumptions used in the monodomain model may break down. The approach presented here also inherits a limitation from the data set used for its parameterization. We only use four cells, and they come from different areas of two different hearts. Without further experiments we cannot say whether this variability arises as a result of abnormal cells, or in some way embodies the variability across one heart. Our conclusions are relevant to the coupling together of these four cells in tissue, rather than directly to the hearts in question.

In the real heart, the diffusion coefficient D is an anisotropic tensor whose principal direction follows the fibrous arrangement of myocytes within the myocardium (Hooks et al., 2007). In our simulations we have considered only isotropic diffusion coefficients. The interaction between fibre orientation and dispersion of repolarization arising from extrinsic and intrinsic variability may merit future study.

Conclusions

We have developed a phenomenological model capable of reproducing both intrinsic beat-to-beat variability in APD from an isolated guinea pig ventricular myocyte, and extrinsic variability in the AP between different cells. In this study, extrinsic and intrinsic variability were incorporated into tissue simulations. We conclude that, for the guinea pig APs used in this study, the degree of inter-cell variability makes a larger contribution to spatial dispersion in APD90 than beat-to-beat variability. The degree of dispersion in APD90 increases as tissue coupling is decreased. Furthermore, if BVR manifests in tissue, then it is likely to be as a result of ectopic activity, instead of the naturally occurring APD90 variability observed in isolated cardiac myocytes.

Table A1

Guinea pig cell-specific parameter sets for the BOCF model.

Parameter	Apex 1	Apex 2	Base 1	Base 2
θv	0.2985	0.3209	0.2749	0.3005
θw	0.1299	0.1302	0.1254	0.1202
θv−	0.2016	0.1943	0.1908	0.1714
θo	0.006	0.0056	0.0058	0.0062
us	0.9502	0.9802	0.9212	0.9391
ww−	0.77	0.77	0.77	0.77
uso	0.6456	0.6475	0.6645	0.7072
uu	1.546	1.617	1.592	1.575
uo	0.0001	−0.0008	0.0004	0.0003
τ_w∞	0.0273	0.041	0.036	0.029
τo1	462.4	416.8	476.7	456.3
τo2	5.951	6.1	5.579	6.042
τs1	2.689	3.057	2.482	2.813
τs2	1.977	2.192	1.764	1.949
τso1	39.32	30.91	39.13	28.81
τso2	1.193	0.2569	1.107	1.287
τw1−	80	80	80	80
τw2−	100	100	100	100
τv1−	75.39	76.26	72.99	66.59
τv2−	10.03	10.2	9.637	10.03
τsi	2.934	2.727	2.902	2.784
τv+	1.462	1.353	1.525	1.333
τw+	275.9	276.3	352.8	273.8
τfi	0.0957	0.1094	0.095	0.0944
w∞⁎	0.95	1	0.95	0.95
kso	1.946	2.006	2.023	1.942
kw−	198.2	204	188	207
ks	2.203	2.293	2.653	2.812
θwo	1.1	1.1	1.1	1.1
V (mV)	94.1	103.8	94.6	90.5
vR (mV)	−73.9	−77.8	−78.1	−76.1