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.


a b s t r a c t
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 intercellular 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.
& 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).

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-tobeat variability of repolarization (BVR) in the same cell. Beat-tobeat 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.  , 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 t on the interval ½t; t þΔt has mean 0 and variance Δt. We modified Eqs. (A.2)-(A.4) by including a Weiner increment in each, forming a system of stochastic differential equations (SDEs): where σ v , σ w , and σ s 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 (APD 90 ). Poincaré plots demonstrate the degree of temporal variability by plotting the APD 90 of each beat against the APD 90 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.

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 to We applied zero-flux Neumann boundary conditions to the tissue. In all simulations, the stimulus was applied at the left-hand edge of the tissue (x o0:03 cm). 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 cm 2 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 cm 2 s À 1 , 0.117 cm 2 s À 1 , and 0.059 cm 2 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 pseudorandom 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 alternanslike 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 APD 90 after each simulation. We did not perform these analyses where x r 0:25 cm or x Z0:75 cm to prevent electrotonic effects arising from the Neumann boundary conditions or the stimulus site influencing the APD 90 arising from the cell model (Cherry and Fenton, 2011). The maximum dispersion of APD was measured as the maximum difference in APD 90 across the region 0.25 cm r x r 0:75 cm at each node in the mesh. The difference in APD 90 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 APD 90 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 . 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  .586 cm 2 s À 1 (green), D ¼ 0.117 cm 2 s À 1 (blue), and D ¼ 0.059 cm 2 s À 1 (black). (B) The difference between the stochastic and deterministic APD 90 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 APD 90 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.) 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).

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 APD 90 between beats increased slightly as the diffusion coefficient D was decreased. This was quantified by calculating the variance in mean APD 90 across all 50 beats (see Table 2). The distribution of mean APD 90 for the stochastic simulations with intrinsic variability was centred on the mean APD 90 of the deterministic simulations (no intrinsic variability) in all cases. The beat-to-beat variance in mean APD 90 also increased with the variance in APD 90 of the underlying cell (Table 2, see brackets after cell name for isolated cell variance in APD 90 ).
Spatial dispersion in the APD 90 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 APD 90 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 cm 2 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 cm 2 s À 1 of up to 0.18 ms 2 in the deterministic case and 0.4 ms 2 in the stochastic case (base cell 1).

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 APD 90 recorded at all four values of D over 160 beats. Mean APD 90 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 APD 90 .
We compute the variance within the APD 90 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 Table 2 Summary of intrinsic variability only simulations μ APD : the mean APD 90 in the simulated tissue over one beat. Mean(μ APD ): average of μ APD over all beats in the simulation. Var (μ APD ): variance in μ APD 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   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 APD 90 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 cm 2 s À 1 ) is bounded by 71.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 cm 2 s À 1 ). In the deterministic simulations, there is almost no spatial heterogeneity in the APD 90 , and similarly with D ¼ 1.171 cm 2 s À 1 almost no dispersion is observed in the stochastic model. At D ¼ 0.059 cm 2 s À 1 , the stochastic simulations show a modest spatially heterogeneous alternation in the simulated APD 90 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 cm 2 s À 1 . The maximum amplitude of the oscillations is approximately 2 ms.

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 APD 90 of the deterministic simulations was constant over all beats when including extrinsic variability between cells. There was a beat-to-beat variability in mean APD 90 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 APD 90 behaviour in tissue.
The variability in APD 90 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 APD 90 shown in Fig. 6B. The variance in the APD 90 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 ms 2 , which remains smaller than each simulated cell's variance in APD, which ranges between 17.64 ms 2 for apex cell 2 and 77.97 ms 2 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 cm 2 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 cm 2 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.

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 cellcell 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 . 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 cm 2 s À 1 ). At the lowest value of D (0.059 cm 2 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 beatto-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. 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. The maximum difference in APD between the stochastic and deterministic simulations at each beat, for each simulation. Values of D are 1.171 cm 2 s À 1 (red), D¼ 0.586 cm 2 s À 1 (green), D ¼ 0.117 cm 2 s À 1 (blue), and D ¼ 0.059 cm 2 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.) 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 APD 90 than beat-to-beat variability. The degree of dispersion in APD 90 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 APD 90 variability observed in isolated cardiac myocytes.
τ so ¼ τ so1 þðτ so2 À τ so1 Þð1 þ tanhðk so ðu À u so ÞÞÞ=2; ðA:10Þ The steady state values are given as ( ðA:13Þ w 1 ¼ ð1 ÀHðu À θ o ÞÞð1 À u=τ w1 ÞþHðu À θ o Þw n 1 : ðA:14Þ The currents in the BOCF model represent an approximation to groups of currents observed in cardiac mycocytes. The fast inward current J fi represents the sodium current I Na that generates the upstroke of the AP. The plateau phase of the AP in the BOCF model is generated by a balance between two opposing currents. The slow inward current J si represents inward currents such as the Ltype calcium current I CaL , and the slow outward current J so represents outward currents active in the plateau phase, such as the slow-rectifier and rapid-rectifier potassium currents I Ks and I Kr . The variables v, w, and s represent the availability of channels generating these groups of currents to be activated. These variables recover following the end of the AP, as determined by the Heaviside functions in Eqs. (A.2)-(A.4). The Heaviside functions within the currents in Eqs. (A.5)-(A.7) approximate generic voltage-gated channel opening and closure. The values of the currents and variables making up the BOCF model during one AP are shown in Fig. 8. The parameter sets used in this study can be found in Table A1.

Appendix B. Convergence analysis for the stochastic tissue model
We examined convergence in the conduction velocity and in the maximum difference in the non-dimensional membrane potential u to determine the numerical parameters used in the  Fig. 9. A spatial discretization of 0.01 cm and a PDE time step of 0.01 ms are sufficient for this simulation study.
We use a conserved Brownian path to investigate convergence properties of stochastic models (Gaines and Lyons, 1997). A conserved Brownian path is a path sampled such that at multiple temporal resolutions it remains a Brownian path. To construct such a path, we sample a Brownian path at our lowest temporal resolution of interest (0.01 ms). This Brownian path b(t) is sampled using a Weiner process whose Gaussian increment W h has mean 0 and variance h on the increment ½t tþ h. Hence W t $ Nð0; ffiffiffi h p Þ. We now construct a new Brownian path from the original Brownian path, which preserves the increments of the old path, but has a temporal resolution of h=2. Consider each time increment ½t tþh in the original path in turn. We bisect this interval to get two subintervals ½t tþ h=2 and ½t þ h=2 t þ h. We now wish to construct a random variable for each sub-increment, such that both new variables are Weiner increments on the subintervals ½t tþ h=2 and ½t þ h=2 t þ h. Their sum must be equal to dW t , to conserve the original Brownian path. Let Z $ Nð0; ffiffiffi h p =2Þ be independent of W h . Then consider the increments 1 2 W h À Z and 1 2 W h þ Z. These are then Weiner increments on the intervals ½t tþ Δt=2 and ½t þ Δt=2 t þ Δt, respectively. Using this method, we can iteratively construct a new Brownian path that preserves the previous Brownian path. This results in a nested series of Brownian paths, all with the same characteristics, at different temporal resolutions, as shown in Fig. 10A.
In order to test convergence we pre-generate a conserved Brownian path for the stochastic variables w and s for each node in a 0.5 cm Â 0.5 cm tissue with a spatial discretization of 0.01 cm. Simulations are performed for a duration of 250 ms with the stochastic version of the apex cell 1 model. Simulations are analysed by calculating the difference in APD 90 at each node between the current simulation and the reference solution obtained at a time step of 0.000625 ms (Fig. 10B). An absolute difference in APD 90 of r 0:1 ms between the reference solution and the 0.00125 ms solution is observed in the D¼ 0.059 cm 2 s À 1 case. Simulations are run using a spatial discretization of 0.01 cm, a PDE time step of 0.01 ms and an ODE/SDE time step of 0.00125 ms.