Regulation of Ion Gradients across Myocardial Ischemic Border Zones: A Biophysical Modelling Analysis

The myocardial ischemic border zone is associated with the initiation and sustenance of arrhythmias. The profile of ionic concentrations across the border zone play a significant role in determining cellular electrophysiology and conductivity, yet their spatial-temporal evolution and regulation are not well understood. To investigate the changes in ion concentrations that regulate cellular electrophysiology, a mathematical model of ion movement in the intra and extracellular space in the presence of ionic, potential and material property heterogeneities was developed. The model simulates the spatial and temporal evolution of concentrations of potassium, sodium, chloride, calcium, hydrogen and bicarbonate ions and carbon dioxide across an ischemic border zone. Ischemia was simulated by sodium-potassium pump inhibition, potassium channel activation and respiratory and metabolic acidosis. The model predicted significant disparities in the width of the border zone for each ionic species, with intracellular sodium and extracellular potassium having discordant gradients, facilitating multiple gradients in cellular properties across the border zone. Extracellular potassium was found to have the largest border zone and this was attributed to the voltage dependence of the potassium channels. The model also predicted the efflux of from the ischemic region due to electrogenic drift and diffusion within the intra and extracellular space, respectively, which contributed to depletion in the ischemic region.


Introduction
Myocardial ischemia is caused by reduced perfusion to regions of the heart leading to a localised reduction in supply of metabolites, limited waste removal and compromised ionic homeostasis. The first 10 minutes of ischemia are associated with an increased risk of arrhythmias peaking after 5-6 minutes [1]. During this period arrhythmias are commonly initiated within the border zone (BZ) separating viable, well perfused, tissue and the ischemic, underperfused, region [2][3][4]. Ischemia causes an increase in extracellular potassium (½K z e ), intra and extracellular proton concentrations (½H z i and ½H z e , respectively), intracellular sodium (½Na z i ) and intracellular calcium (½Ca 2z i ) concentrations [5]. The dominant mechanisms for these changes have been attributed to a shift in the ATP/ADP ratio, which inhibits the Sodium-Potassium ATPase pump (NaK) and increases the conductance of ATP-inactivated K z channels; respiratory acidosis causing an increase in CO 2 ; and metabolic acidosis, where a shift towards anaerobic respiration increases the production of H z in the cell [5]. Inherently, these changes in ionic concentrations in the ischemic region lead to gradients in properties across the BZ, creating electrophysiological heterogeneities that are thought to favour the occurrence of arrhythmias [6][7][8].
Experimentally, the development of gradients of extracellular pH (pH e ) and ½K z e [9,10] have been well characterised using ion sensitive electrodes. Intracellular metabolite gradients have been characterised by fluorescent NADH [11] and biopsy [12] measurements. However, less is known on the gradients of intracellular ions, in particular ½Na z i , ½Ca 2z i and ½H z i , nor are the mechanisms that underpin the spatial and temporal evolution of these ion concentration gradients well characterised or understood. This study aims to investigate the spatial-temporal evolution of ion gradients across ischemic BZ and the primary regulators of the BZ size and rate of development.
Previous measurements of ion concentrations and metabolites across the BZ have either been performed at multiple locations but at a limited number of time points [10,12] or have tracked the time evolution of ion concentrations but only from a limited number of locations [13,14]. Furthermore, these measurements have only been able to characterise a subset of ions of interest across the BZ. The need to track the evolution of multiple ionic species in space and time to understand the gradients of cellular electrophysiology across the BZ motivates the use of biophysical computational modelling. Previous models of electrophysiology during acute regional ischemia have simulated the effects of these spatial gradients but have not simulated their time evolution [15][16][17]. More recent work has simulated the time evolution of ½K z e gradients [18], but have not considered other ionic gradients, the effects of nonlinear interactions between K z and other ions, the effect of K z diffusion in the intracellular space or the effects of potential gradients on ion diffusion.
In this study a new model of cardiac tissue electrophysiology is developed to investigate the spatial-temporal evolution of ionic concentrations across the ischemic BZ, during the first 5 minutes of reduced perfusion. The proposed model extends the conventional bidomain equations to explicitly link membrane potential to ionic concentrations and enforces ionic species conservation. A model of ion regulation across the cell membrane is then developed, parameterized, validated and coupled to the tissue model. This combined model is then used to investigate the spatial and temporal dispersion of ions across the BZ.

Methods
To model the evolution of ionic concentrations in the presence of multiple ionic gradients, electric gradients and heterogeneous tissue properties requires the development of a new set of equations for modelling the myocardium. In the next section the equations to model the BZ are derived and a model of passive cell membrane ion regulation is developed and validated. The changes to the cell membrane model to simulate NaK inhibition, K z channel activation, respiratory acidosis and metabolic acidosis are then described.

Tissue Model Derivation
Consistent with previous models of cardiac tissue electrophysiology the myocardium is represented as a two phase medium, with each point in the domain containing a fraction of intra and extracellular space. This gives rise to the scaling variables.
where V x is the volume of the x space in a unit of V myocardium volume, a x is the volume fraction of space x (with x corresponding to the extra (e) or intra (i) cellular space), x is the interspace surface area per unit volume myocardium and C is area of the interspace surface. The intracellular volume fraction (a i ) can be separated into sub volume fractions representing the cytosol (a cyt ), mitochondria (a mito ) or the sarcoplasmic reticulum (a SR ) to model distinct subcellular spaces, as described below. At each point in space there exists an intracellular potential (w i ), an extracellular potential (w e ), a transmembrane potential (V m ) and an intra and extracellular ion concentration for each of the ion species in the model. The movement of each ion species in each domain is driven by diffusion, due to a gradient in ion concentration, and drift, due to a gradient in the electric field. This movement is described by the Nernst-Plank equations where C j x is the concentration of unbound ion j in space x, z j is the charge of ion j, D j x is the effective diffusion of ion j in space x, F is Faraday's constant, R is the gas constant, T is the absolute temperature and w x is the potential in the space x. The conventional Nernst-Plank equations are adapted to represent the movement of ions across the cell membrane between the intra and extracellular spaces. This gives where J j m is the flux of ion j across the cell membrane. J j m is defined in units of mM ms {1 mm. The surface separating the two regions is modelled as a simple capacitor and defining the transmembrane potential to scale with intracellular charge gives where Q i is the intracellular charge per unit cell membrane area and C m is the membrane capacitance per unit cell membrane area. The charge on either side of the membrane is assumed to be equal but opposite, giving Separating the concentration of each ion species in the intra and extracellular space into ions that are membrane bound and make up the membrane charge (Q j x ) and those that are in solution gives.
where C j xsol and C j xQ are the ions in solution and membrane bound ions, respectively. By assuming charge neutrality for the ions within the solute gives: Multiplying Eqn 6 by the ion species charge and Faraday's constant, then summing over all ion species in space x gives.
Converting the concentration of ions per unit volume in domain x to charge per unit cell membrane area and introducing a static charge term (g x ) that characterises all charge not attributable to C j x , gives Combining with Eqn 8 gives This defines the charge on the membrane as equal to the unbalanced charge in space x. Using the charge balance Eqn 5, gives Differentiating Eqn 11 with respect to time, substituting in Eqn 3 and recognising that all transmembrane fluxes are balanced, provides this ensures that there is no net charge accumulation in any unit volume of myocardium. Defining the relationship between the intra and extracellular potentials gives Using Eqn 4 and the definition of charge (Eqn 10), then gives the algebraic definition of the transmembrane potential: Rearranging Eqn 13 and substituting in the definition of V m from Eqn 14 allows w i to be defined in terms of w e and C j x . Combining this definition of w i with Eqns 14, 13, 12, 3 and 2 then represents a closed set of equations. These equations are equivalent to the bidomain equations in the case of a single charge carrier and no gradient in ion concentrations, as shown below.
In cardiac myocytes many important ions, including H z and Ca 2z are heavily buffered, both within the cell and the extracellular space. To account for buffering the free and buffer bound fraction of C j x are calculated. In general, as C j xsol wwC j xQ the effect of ions bound to the cell membrane will not be included in the buffering equations for simplicity. This gives where C j xfree are the unbound ions and C j xbuff are the ions bound to buffers. At this time all buffers will be treated as rapid and to a single representative buffer species, giving where B j x , K j x and n j x are the concentration, binding affinity and Hill coefficient, respectively, for the buffer of ion j in space x. Assuming that ions bound to buffers are immobile, only C j xfree is used to calculate the diffusion and drift of ions in Eqn 3, and similarly in Eqn 12. As binding of ions to a buffer implicitly removes a charged binding site located on a static protein, Eqn 14 remains unchanged.
Due to the complex anatomy of the cardiac myocyte many sub volumes exist within the cell that affect ionic concentrations. The sum of the volume fraction (a x ) values must be less than, but do not have to be equal to, one, allowing the model to represent any volume fractions that are not directly accessible by ions. In particular the a i variable can represent all space in the cell or can be substituted for a cyt representing the volume fraction of the cytosol (a sub volume of a i ). This allows the effects of SR, mitochondrial or other subcellular structure volumes on intracellular ionic concentrations to be accounted for in the model. Equation summary. The modelled equations are given by where it is important to note that for non buffered ions B j x~0 and C j x~C j xfree . Consistency with bidomain equations. Imposing the implicit assumptions of the bidomain equations that charge carriers are not buffered, ion concentrations are homogenous and charge is carried by a single carrier to Eqn 18, the bidomain equations can be derived. Assuming homogenous ion concentrations and considering the case of the intracellular space reduces Eqn 18 to Differentiating Eqn 14 for the intracellular space gives substituting Eqn 18 into Eqn 19 gives As ion concentrations are homogenous conductivity is defined as Then Eqn 20 reduces to Applying the single charge carrier assumption and converting from ionic flux to current gives the first bidomain equation The second bidomain equation is readily derived from applying the homogenous ion concentration assumption to Eqn 12 and multiplying by Faraday's constant (to convert from conserving ion flux to current) giving Substituting in the definition of conductivity from Eqn 21 then gives the second bidomain equation

Modelling the Membrane Fluxes
Cardiac electrophysiology is predominantly determined by the movement of Ca 2z , Na z and K z . For charge neutrality Cl { must also be included in the model. To simulate the evolution of acidosis requires the inclusion of H z , HCO { 3 and CO 2 in the model. All of these ions (and CO 2 ) were modelled in the intra and extracellular space, with H z and Ca 2z being buffered. The goal of the model, in this study, was not to track the propagation of the action potential but to simulate the gradients of ions that exist over the BZ. These ion gradients were modelled based on the diastolic properties of the cell. This assumption was also a requirement to enable the simulation of minutes, while remaining computationally tractable.
The membrane ion transport pathways are described first for Na z , H z and HCO { 3 . Ion specific channels are then described that balance the flux of each ion species. The channel and transporter densities were determined by imposing zero net flux for each ion species, using the relative densities of Na z and H z transporters recorded experimentally, and intra and extracellular ionic concentrations and membrane potential values derived from the literature. Where possible experimental data was taken preferentially from rabbit or guinea pig data at body temperature. This limited number of constraints then allowed the model transporters and channel densities to be uniquely determined.
Sodium regulation. The model of Na z regulation included representations of the NaK, sodium calcium exchanger (NCX ), sodium hydrogen exchanger (NHE), sodium bicarbonate cotransporter (NBC) and a lumped sodium channel (I Na ). NaK was modelled using the thermodynamically consistent equation set proposed by Smith and Crampin [19]. This model was subsequently revised by Terkildsen et al., [20] and this parameter set that was used here. The model for NaK was fitted to guinea pig data as limited rabbit data was available. However, the maximum flux was rescaled to match rabbit data, as described below.
The NCX model was taken from Weber et al., [21]. The model has been fitted to rabbit experimental data at 37uC. The NHE model was based on the model developed by Crampin and Smith [22] and reparameterized by Niederer and Smith [23]. In this study extracellular Na z and H z regulation of NHE were included. This model was fitted predominantly to sheep Purkinje data [24], although the H z dependence of NHE remains relatively consistent between species [25]. The NBC model was taken from Crampin and Smith [22]. The model assumes NBC is electro neutral, which is true for only part of the NBC population [25]. There was not sufficient data to fully characterise the electrogenic and electro neutral forms of NBC, hence the electro neutral model was used. Background Na z flux across residual open fast Na z and persistent Na z channels is limited when the cell is quiescent. However, some flux is still present [26] and a simple lumped background ionic flux equation was used, given by to model the residual Na z flux across any open Na z channels. The same equation form was used for modelling all background ion channels. Proton regulation. In the proposed model ½H z i was regulated by NHE, described above, chloride-hydroxide exchanger (CHE), hydrolysis and buffering. Background H z leak or other H z exchangers were not considered in the general model of H z regulation, described here, but do include models of H z -lactate exchange and intracellular metabolism derived H z sources in the model of ischemia, described below. The CHE model comes from Niederer et al., [27] and was fitted to guinea pig data at 37 o C. The hydrolysis of CO 2 into HCO { 3 and H z was governed by where k f and k r are the forward and reverse rates of hydrolysis. Hydrolysis occurs in both the intra and extracellular space and the rate constants were assumed to be the same in both domains. H z buffering in the intra and extracellular space is due to mobile and static H z buffers and HCO { 3 [28]. The intra and extracellular buffering of H z were assumed to be instantaneous and represented by a single population of buffers. To reduce the model size, partial differential equations were only solved for the total concentration of ½H z i or ½H z e . The concentration of free ions were then calculated by where B H and K H represent the concentration of the buffer and the binding affinity, respectively. This buffer model was used for both intra and extracellular H z , with a separate set of parameters for each ion. Calcium regulation. A simplified model of intracellular cardiac Ca 2z was developed assuming that Ca 2z in the intracellular space reaches an approximate equilibrium over the time scales of interest. Furthermore, SERCA ATPase function was modelled with a Hill coefficient of one as opposed to two, to allow the definition of Ca 2z to remain deterministic. The intracellular space was assumed to consist of a sarcoplasmic reticulum (SR) and a cytosolic space. The subsarcolemmal and dyadic space are small and are also likely to be in equilibrium with the cytosolic Ca 2z , so were not included in the model. The Ca 2z dynamics were described by where ½Ca 2z SR is the SR Ca 2z , J leak is the flux of calcium out of the SR, J SERCA is the uptake of Ca 2z by SERCA, V SERCA is the maximum SERCA flux, P leak is the diffusion permeability of the SR membrane, K SERCA is the binding coefficient of Ca 2z to SERCA, a SR and a cyt are the volume fractions of the SR and cytosol, respectively, I CaL is the L-type calcium channel and I Cab is the background Ca 2z channel. In this model the background and L-type Ca 2z channels were modelled as a single lumped generic Ca 2z channel. Introducing cytosolic buffering, ignoring the effects of SR buffering and setting nH as one and assuming that the cytosol and the SR are in equilibrium then gives and collecting terms gives As C is always negative the cubic always has at least one positive real root for possible values of ½Ca 2z tot . The value of ½Ca 2z tot was then found using the root finding method first proposed by Francois Viete in 1600 and reused more recently by Faber and Rudy [29]: The model of Ca 2z dynamics assumes that all Ca 2z buffers were static and that the transport of ions via mobile buffers was accounted for in the effective diffusion parameters of free Ca 2z . It is possible to extend the proposed model to include mobile buffers but they were assumed to play a secondary role in the current model. ½Ca 2z e was assumed to be buffered by a single species and was modelled using the same framework described above for H z (Eqn 28).
Chloride regulation. In this model Cl { homeostasis was maintained by the CHE and the Cl { -HCO { 3 exchanger (AE), which bring Cl { into the cell, and a Cl { channel that allows Cl { to flow out of the cell. The CHE model is described above. The AE model was taken from Crampin and Smith [22] and was developed using guinea pig data at 37 o C. The Cl { channel uses the conventional background channel formulation. A linear H z dependence of the background Cl { current was added to the model based on observations from Komukai et al., [30].
Potassium regulation. In dynamic action potential models of cardiac electrophysiology there are a large number of K z channels [31] that bring K z into the cell. This influx was balanced by the NaK pump, described above, that extrudes K z . For the passive membrane model, all of the K z channels were lumped into a single background current (I Kb ) formulation that was set to balance the flux of K z on NaK. It was assumed that the ð38Þ membrane potential and K z reversal potential are the dominant factors affecting this channel and other forms of regulation have not been considered.
Bicarbonate regulation. HCO { 3 was assumed to be regulated principally by hydrolysis and through NBC and AE. All of these components have been described above and it was assumed that there are no other HCO { 3 pathways across the membrane. Carbon Dioxide. CO 2 is regulated primarily through hydrolysis and can diffuse relatively freely across the membrane. The model of hydrolysis is described above and CO 2 diffusion was assumed to obey Fick's law.

Model Parameters
For each transmembrane ion pathway described above, all kinetic, binding affinity and membrane potential dependencies were taken from the original models. Here the definition of geometrical parameters, ionic concentrations, buffering parameters and the density/scaling of each transmembrane pathway are motivated from data in the literature.
Geometrical parameters. The extracellular space is estimated to be between 17{30% [32-35], 17:7+0:6% and 24:6+0:6% [33,36] of the volume of the heart in rabbit, cat, and rat hearts, respectively. This gave an a e value of 0:2 leading to an a i value of 0:8. The surface to volume ratio of a cell is reported as 0:31{0:35m m {1 [37,38] in rat and rabbit myocytes, corresponding to a x value of 264mm {1 . The relative SR volume was set to 3% of intracellular volume, giving an a SR value of 0:024, based on reported values of 1:4{3:5% of cell volume [37,39,40] in mouse, rat and swine. The relative mitochondrial volume (a mito ) was set to 0:24 based on an estimated mitochondrial cell volume fraction of 30% [37,40]. The cytosol volume fraction was set to 67% of the intracellular space, resulting in an a cyt value of 0:536. To account for the effects of subscellular domains on intracellular ionic concentrations in the model, all references to a i in Eqn 18 were replaced by a cyt . A summary of geometrical parameters is given in Table 1.
Intracellular ionic concentrations. ½Na z i has been measured using SBF1 fluorescence and Na z sensitive electrodes. A significant range of values have been reported from 2{10mM [26,[41][42][43][44] to 20{30mM [34,45,46]. Early measurements of intracellular ionic concentrations were performed using ion sensitive electrodes. These experiments measure ion activity and not ionic concentration and are often performed in multi-cellular preparations, confounding measurements. For these reasons ½Na z i in quiescent myocytes was set to 4mM, consistent with recent calibrated fluorescent measurements in isolated rabbit myocytes [26,44].
No fluorescence dye is routinely used for measuring ½K z i . Using ion sensitive electrodes Lee et al., [45] were able to calibrate their measurements of ion activity in rabbit myocytes using an estimated K z activity coefficient of 0:613, giving a value of 135mM. This compares with a range of 147{236mM calculated by applying the Lee et al., K z ion activity coefficient to ion activity measurements in rabbit, cat and guinea pig [34,43,47,48]. Alternate measurement using flame emission spectrometry by Powell et al., [49] measured ½K z i in rat myocytes, giving a concentration of 120:8+1:7 mM. Given the lower values of the two calibrated measurements, ½K z i was set to 135 mM.
No dye is routinely used for measuring Cl { concentration in cardiac cells, however, ½Cl { i can be measured using ion sensitive electrodes. Cl { activity has been reported as 14{20 mM [35,48,50,51] in sheep, rabbit and guinea pig heart cells. Estimations of ½Cl { i from total tissue Cl { concentrations have resulted in values of 38:9+2:5 mM [35] and 16:2 mM [34] in rabbit cells and 19:4+1:2m mol/g dry wt (or 9:7 mM using the 0:49 mM per mmol/g dry wt scaling factor from Bers [52]) in rat. The higher value of 39 mM may be attributed to the higher extracellular space used in these calculations (30% as compared to 20%). Considering the relative convergence of values ½Cl { i was set to 18 mM.
Resting free Ca 2z is measured using calibrated fluorescence measurements. These measurements range from 44{100 nM in rabbit and guinea pig preparations [53][54][55][56][57]. Given this consistency the ½Ca 2z i will be set to 80 nM. SR Ca 2z concentration is calculated from integrating the current across the cell membrane following the release of Ca 2z from the SR in response to caffeine. These measurements show two populations with high values in rat (73{140 mM [58][59][60][61]), canine (80 mM [62]), rabbit (87-106 mM [60,61]) and ferret (76 mM [63]), compared to lower measurements in guinea pig (17:9{27 mM [58,64]). Given that the majority of species have a higher reported concentration, including rabbit, simulations were run with SR Ca 2z load set to 100 mM. The buffering of Ca 2z can be described by Hill equation(s), mass action equation(s) or a constant buffering power. To compromise between biophysics and complexity, buffering was modelled by a single Hill equation. Hove-Madsen and Bers [65] fitted Ca 2z buffering in rabbit myocytes using two Hill curves; however, the lower affinity buffer will not play a significant role at passive diastolic Ca 2z concentrations. For this reason the high affinity site, with cooperativity reduced from 1:27 to unity, was used to model Ca 2z buffering, giving a buffer concentration of 208.98 mM (converted using a scaling factor of 2:43 from Bers [52]) and an affinity of 0:42 mM. was modelled explicitly and the intrinsic buffers were assumed to be in rapid equilibrium. Leem et al., [70] measured (and modelled) H z buffering by two populations of buffers with binding affinity pK values of 6:03 and 7:57 and concentrations of 84:22mM and 29:38mM. Zanbioni et al., [25] differentiated between mobile and fixed buffers and found that the fixed buffers had a consistent concentration of 50{70mM and binding affinity (pK) value of 6:1{6:2 across rat, rabbit and guinea pig, while the mobile buffers had a constant b value of approx10. Fitting a single buffering curve to these two models over a pH range of 6:5{7:5 gives concentrations of 64{68mM and pK values of 6:45{6:65, which gave a value of 65mM and pK value of 6:5 for this model. Extracellular ionic concentrations. The concentration of the majority of ions in the extracellular space have been measured in canine [46,72], rat [36,73] and cat [74] hearts. These measurements provide a consistent range of ion concentrations for Na z , K z and Cl { , giving ½Na z e as 134{159mM, ½K z e as 3:3{5:8mM and ½Cl { e as 105{120mM. In the model ½Na z e was set to 140mM, ½K z e was set to 4mM, consistent with measurements in guinea pig hearts [75] and rabbit atrium [76] and ½Cl { e was set to 110mM. The ½Ca 2z e values reported range from 1:03{5:16mM, however, these values do not differentiate between buffered and ionized Ca 2z . The properties of Ca 2z buffering in the extracellular space are not well characterised and are generally ignored in previous cardiac cell models. To approximate the buffering properties of extracellular Ca 2z using a simple single species steady state mass action model, with no cooperative binding, requires two parameters, the concentration of the buffer and the binding affinity. Assuming the ratio of free Ca 2z in the extracellular space to bound ions is similar to serum [77] and assuming that the primary buffer of Ca 2z in the extracellular space are phospholipids, then extracellular Ca 2z buffering will have a binding affinity of 1:1mM [78], within the range observed across multiple species [79], and a buffer concentration of 2:3mM based on a free Ca 2z concentration of 1:2mM and assuming 50% of Ca 2z are bound to buffers [77]. pH e was set to 7:4, to be consistent with the majority of experimental studies [57,[69][70][71] and measurements across a range of species [80]. Limited measurements were available to model ½H z e buffering, however, Yan and Kleber [81] reported 39 mM of buffered H z in the extracellular space. By assuming similar binding affinities for the intra and extracellular buffers and that pH e is 7:4, then gave a concentration of extracellular proton buffers of 350mM. A summary of ion concentrations and buffering parameters are given in Table 2 and Table 3, respectively.
The concentration of CO 2 was calculated using the parameters proposed and measured for guinea pig ventricular myocytes at 37 o C by Leem and Vaughan-Jones [82]. The concentration of CO 2 in the extracellular solution was calculated using where a CO 2 is the solubility of CO 2 , set to 0:0307 mM mmHg {1 from human measurements at pH 7:4 and 37 o C [83], %CO 2 is the fraction of air that is CO 2 , set to 5% at baseline and P atm is atmospheric pressure, set to 760 mmHg {1 . This gives a partial pressure of CO 2 (P CO 2 ) of 38 mmHg, consistent with although slightly higher than the 34+3 mmHg measured in rabbit hearts [84]. The hydration of CO 2 was modelled by a mass action reaction (Eqn 27  [82], based on measurements in red blood cells. This value was reused despite its lack of species and cell type consistency, as there were no recent studies characterising the permeability in cardiac myocytes and the high permeability leads CO 2 to be close to equilibrium between the intra and extracellular spaces. Current, exchanger and pump densities. The model of each transmembrane ion pathway (y) was separated into a kinetic regulatory component (c y ) dependent on transmembrane potential and ionic concentrations, and a scalar (V y ) representing either the maximum flux or channel conduction of the pathway. The flux across a pathway (G y ) is then given by The V y values were determind by calculating all of the c y values (excluding c NBC , c hyd and c CO2 as ½CO 2 i and ½HCO { 3 i were unknown) using the ionic concentrations in Table 2 and assuming a membrane potential of 280 mV. The remaining unknown V y and ionic concentrations were then determined from a limited number of measurements and by enforcing a zero net flux condition described by  As the parameters for hydrolysis and permeability of the cell membrane to CO 2 (P CO2 , k f and k r ) were derived from the literature, the constraints on V CO2 and V hyd were used to determine the concentrations of ½CO  [85]. The G NaK Na z flux was assumed to be equal to the sum of all Na z influx. As the model does not include Na z -K z -2-Cl { co-transporter (NaK2Cl) or Na z -Mg 2z exchanger (NaMg), primarily due to the limited data to constrain the kinetics, G NaK was reduced from the flux measured by Despa et al., [85] to 1.13|10 {5 mM ms {1 . Setting these G y values defines V Na , V NCX , V NHE , V NaK , V kb and V Cab . In the experimental and modelling work by the Vaughan-Jones group [70,86] the relative size of G CHE to G NHE is 0:12{0:19 at pH i~7 :1, this gave an estimated ratio of 0:16. Scaling G NHE determined from [85] by 0.16 provided an estimate of G CHE of 0:032|10 {5 mM ms {1 . This allows V CHE and G hyd to be defined. G hyd must be balanced by G CO2 and as P CO2 and ½CO 2 e were known, this flux was used to set ½CO 2 i . Combining the pH i and ½CO 2 i concentrations with the defined hydrolysis parameters and the known value of G hyd then gave ½HCO { 3 i . Knowing G hyd and G NBC gave G AE and hence V AE . Similarly, G Cl was calculated from G AE and G CHE ,which gave V Cl . Finally, V NBC was calculated using G NBC and c NBC , calculated using the derived ½HCO { 3 i value. By automating this parameter derivation process the model parameters could be updated to ensure static ion concentrations, and hence membrane potential, for any perturbation in model parameters or ionic concentrations. Parameters for all simulations in this study were derived following this process.
Intracellular calcium dynamics. The intracellular regulation of Ca 2z was treated as an equilibrium system. This resulted in the SR effectively acting as an additional buffer on ½Ca 2z i . The parameters V leak , V SERCA , K SERCA , a SR and a cyt were defined from enforcing a zero net flux constraint and experimental measurements. Measurements of K SERCA are similar in rat and  Table 4. Black arrows indicate where the model matches experimental data, white arrows indicate where there is no or inconsistent data, gray arrows indicate where no change was observed experimentally and striped arrows indicate where the model does not match experimental data. Data was considered inconsistent if different studies reported opposite changes in intracellular ion concentration or membrane potential in response to a change in extracellular ionic concentrations (for example the change in membrane potential following a decrease in ½Na z e ). In cases where one study reported no change in intracellular ion concentration or membrane potential and another found a change, it was assumed that the change was correct (for example the change in ½K z i in response to an increase or a decrease in ½K z e ). For example in panel A) corresponding to a decrease in ½Na z e (as indicated by the panel label), the experimental data comes from row 1 of Tableô 4. The model predicts that ½Na z i decreases and that ½Ca 2z i and ½H z i increase, consistent with experimental measurements and the arrows are shaded black. There is no consistent observed change in the transmembrane potential so the V m arrow remains white. No data is available for the change in ½K z i in response to a decrease in ½Na z e so the K z arrow remains white. The model predicts a decrease in ½Cl { i yet experimental measurements found an increase, hence the Cl { arrow is striped. doi:10.1371/journal.pone.0060323.g001 rabbit [65] and have been reported as 0.260-0.350mM in rabbit [65,87] and 250-280mM in rat [65,87,88]. A value of 0.3mM was used in the model. J Leak is reported as 0.001-0.012 mMms {1 in rabbit, mouse and rat [54,[89][90][91]. A value of 0.01mMms {1 was used in the model. Enforcing a zero net flux constraint on the SR gave J SERCA as 0.038 mMms {1 , comparable with values of 0.03-0.08mMms {1 measured in rabbit [65,87,92] with units converted using scale factors from Bers [52].

Membrane Model Validation
To validate the passive membrane model, a series of tests based on the response of the cell models intracellular ionic concentrations and transmembrane potential to changes in extracellular ionic concentrations or inhibition of major ion transporters was performed. In order to maximise the number of tests the model was compared against data from cardiac cells, regardless of temperature, species or preparation type. This maximised the number of tests but meant that a quantitative comparison was not valid and so only a qualitative comparison was performed. A summary of experimental data used in the validation is provided in Table 4. Fig. 1 shows the comparison between the model and data. From the 72 simulations performed 43 matched the experimental data, data was not available or was inconsistent for 19, changes were too small to be measured in 5 and the model did not match experiments in 5 cases.

Diffusion Parameters
The tissue model required the definition of diffusion parameters for each of the ionic species in the intra and extracellular space. The diffusion parameters used in the model were all for the apparent diffusivity of free ions, lumping the effects of any mobile buffers, tortuosity and gap junctions into a single diffusion parameter. Using Eqn 21 the conductivity parameters were shown to be directly related to the conductivity parameters in the bidomain equations. From the review of bidomain conductivities and relative values in the intra and extracellular space by Roth [93] a conductivity value of 0.25 Sm {1 was taken for both the intra and extracellular space. This conductivity, along with the ionic concentrations from  [94] was used and the effect of diffusion in the SR was not included, as it was assumed to be non contiguous between cells. It was assumed that Ca 2z diffusion is limited by buffering, resulting in the same value in the intra and extracellular space. Na z , K z and Cl { were all assumed to have a common diffusion coefficient in either the intra or extracellular space as these ions are only nominally buffered and are only affected by gap junctions and tortuosity. Diffusion of CO 2 was estimated from previous modelling/experimental studies of CO 2 in tissue as 11:3|10 {7 mm 2 ms {1 [97], and was assumed to be the same in the intra and extracellular space. Solving Eqn 21 for the diffusion of intra and extracellular Na z , K z and Cl { gives the diffusion parameters summarised in Table 5.
Finally, the effects of ½H z i on gap junction conductivity were included. Measurements in cell pairs by Swietach et al., [98] have demonstrated that the permeability of gap junctions has a biphasic dependence on the ½H z i concentration. To introduce these effects into the model the intracellular diffusion constants were scaled by where n1/K d1 and n2/K d2 are the cooperativity and binding affinity for the activation and deactivation of gap junctions by protons, respectively and c diff scales all intracellular diffusion constants. Using the measured parameters from end-to-end cell pairs in Swietach et al., [98] n1~2:062, n2~1:81, K d1~1 0 {3:959 mM and K d2~1 0 {3:934 mM. In this model of proton effects on gap junction permeability there is an implicit assumption that the gap junction permeability plays a dominant role in defining diffusion. At this stage no other regulators of permeability were included in the model, notably Ca 2z regulation is absent but this can readily be included in the modelling framework as required.

Simulating Ischemia
In this study the effect of NaK inhibition, I Kb activation and respiratory and metabolic acidosis during ischemia on cell ionic homeostasis were considered. This list is not exhaustive and absent factors are discussed below. Ischemia was modelled by respiratory acidosis, metabolic acidosis, increased K z channel conductance and decreased NaK function. The specific time course of each of these changes is poorly characterised. Previous modelling studies have assumed that changes in cell function with ischemia have evolved linearly [99] or as a nonlinear function of prescribed metabolite concentrations [20]. To avoid any undue bias from the arbitrary selection of a time course, all changes are initially considered instantaneous.
NaK inhibition and K z channel activation were modelled by scaling the respective fluxes. Respiratory acidosis was assumed to result from an imbalance of production and washout of CO 2 . To simulate this, an (implicitly electro neutral) intracellular source of CO 2 was introduced into the ischemic region. There was no flux of CO 2 out of the extracellular space and this resulted in a build up of CO 2 in both the intra and extracellular space in the ischemic region.
To model metabolic acidosis required the introduction of an additional intracellular H z source. However, introducing a source of cations into the cell compromised the conservation of charge constraint. To provide an electro neutral source of H z required the concurrent introduction of a source of anions that match the production of H z inside the cell. It was implicitly assumed that the anions entered the cell in an electro neutral form with an H z bound, and the anions and H z separate within the cell due to metabolic processes. The source of H z in ischemia is likely to be due to increased ATP production through glycolysis [100]; in the absence of oxygen this also results in increased lactate production. As lactate readily dissociates from H z it was assumed that the anion source that matches H z flux has similar characteristics to lactate including transmembrane regulation via the lactate-H z membrane exchanger MCT1 [101]. Given the ambiguity in structure, MCT1 was modelled as an ordered exchanger, using the kinetic parameters from Vinnakota and Beard [102] and the maximum influx value of 4:8|10 {5 mMms {1 recorded in guinea pig myocytes [103]. In the absence of metabolic acidosis, lactate concentration was assumed to be nominal, consistent with the small flux of lactate observed in rabbit hearts under normal conditions [104]. Model simulations were performed on a 1D strand with length 32 mm, oriented in the preferred conduction or fibre direction, with a transition between ischemic and viable tissue at 16 mm to ensure the simulation captured the BZ width with nominal boundary condition artefacts. The transition between viable tissue and ischemic tissue was approximated by a Hill equation with a Hill coefficient set to ensure a steep transition between the ischemic region and viable tissue over 440mm, consistent with the rapid drop in oxygen pressure across the BZ in swine [105].

Numerical Methods for Tissue Model
The nonlinear equations were solved using a fully implicit finite difference scheme with a line search Newton-Raphson method. The transmembrane flux component of the Jacobian was calculated using finite differencing with the remainder calculated analytically. The Jacobian was inverted directly using MatLab and was only recalculated if the residual fails to decrease or convergence was not reached within 10 Newton-Raphson iterations.
The dependence of the BZ width on the spatial and temporal discretizations was determined to test for numerical convergence. The width of the transition for each ion concentration between the viable and ischemic region, referred to as the BZ width for each ion, was calculated by fitting a Hill curve to each ionic profile. The width of the BZ for each ion was then calculated as the distance between 10 and 90% of the change in concentration.
A convergence analysis was performed and a mesh discretization of 50mm and a time step of 1000ms was used. Increasing the spatial discretization by a factor of 2 results in a maximum change in BZ width of any ion of 0:001mm. Decreasing the time step to 100ms increases the maximum BZ width of any ion by 0:001mm.

Results
The individual effects of each of the four components of ischemia were first demonstrated. A combined model of ischemia was then developed and the width and magnitude of the changes in ionic concentrations across the ischemic BZ were predicted. The effect of movement of K z , Na z and Cl { within and between intra and extracellular spaces across the BZ were calculated to show the net movement of ions across the BZ.

Simulating Individual Components of Ischemia
The level of inhibition of NaK, the activation of I Kb and the level of CO 2 and H z production in the cell during ischemia is not known. To investigate the effects of each of these aspects of ischemia, they were each individually introduced into the model at five levels of severity, for xw16mm along a 32mm strand. Fig. 2 shows the effect of each component on the change in ½Na z i , ½K z e , ½Ca 2z i and pH i across the ischemic BZ. The range of alterations in ion concentrations is shown by the shaded regions. The minimum and maximum changes in each component are shown by dashed lines and the mid change is shown by the solid line. NaK (yellow) was inhibited by up to 100% (no flux) in 20% increments, I Kb (red) was scaled by up to 250% in 50% increments, metabolic acidosis was simulated by introducing an H z flux in five increments up to a maximum value of 0:1mMms {1 and respiratory acidosis (blue) was simulated by a intracellular CO 2 flux increased in five increments up to a maximum value of 0:12mMms {1 .

Modelling Ischemia
Inherently, there is no single mode of ischemia and the relative contribution of acidosis, NaK inhibition or I Kb activation will depend on the residual flow, age, gender, disease state, location of ischemic region and species under study. To provide a representative case to study, 5 minutes of ischemia were simulated in the presence of all four ischemic mechanisms that match representative results from the literature.
Partial pressure measurements of CO 2 in canine ischemic models show a 100{300% increase in CO 2 partial pressure following occlusion, depending on the level of flow inhibition, after 10{15 minutes. The elevation in CO 2 was approximately linear over time and hence in the model it was assumed CO 2 concentration increases by 150% increase in the first 5 minutes of ischemia [106,107]. This corresponded to an increase in the concentration of CO 2 in the model from 1:1mM to 2:8mM. In the model respiratory acidosis was caused by an increase in ½CO 2 i flux that cannot be vented from the extracellular space. Although it is recognised that the decrease in pH e due to metabolic acidosis will also contribute to elevated ½CO 2 e , initially the CO 2 flux was set at a 60% level to achieve an increase in CO 2 concentration to 2:2mM, which resulted in a decrease of pH i to 7:0. During ischemia pH i decreases rapidly before plateauing after approximately 15 minutes. The decrease in the initial 5 minutes has been reported to fall between 0:2 to 0:5 pH units in rat [108], ferret, [109] and guinea pig [110] preparations. Respiratory and metabolic acidosis will both contribute to this drop in pH i . In the model setting the level of intracellular H z flux to 30% caused a 0:17 pH unit drop and an increase in ½CO 2 e to 1:6mM. Combined metabolic and respiratory acidosis caused a 0:28 unit drop in pH and an increase in ½CO 2 e to 2:85mM.
Given the inhibition of NaK and levels of respiratory and metabolic acidosis a sweep of I Kb activation values was performed in the presence of these changes to achieve the desired level of extracellular K z accumulation. During ischemia, ½K z e increases over three characteristic phases, with the first phase occurring during the initial 5 minutes of ischemia prior to reaching a plateau from minutes 5 to 10, before increasing again. In guinea pigs, ischemia caused ½K z e to increase from 4:2 to 12:6mM over 10 minutes [122] or 4 to 5:8mM over 6 minutes [123]. In swine, ischemia caused ½K z e to increase from 3:3{4:2mM to 10:1{11:5mM after 7{8 minutes [124,125]. In rat, ischemia caused an increase from 4:2mM to 8mM over 5 minutes [110], although other groups have seen a biphasic change in ½K z e in rat with an increase from 5mM to 8:3mM before falling back to 6mM then continuing to rise again, observed over the first 8:2 minutes [126]. In rabbit, ischemia caused ½K z e to increase from 4{4:5mM to 6{9:4mM over 5{10 minutes [10,110,125]. In canine, ½K z e increases from 3:2{4:2mM to 8:1{14:5mM after 6{12 minutes [127,128]. These results indicate an increase in ½K z e from 4mM to 10{14mM after 10 minutes and assuming 80% of this rise occurs in the first 5 minutes [127], and combined with the 5 minute data gives an estimate of ½K z e after 5 minutes of ischemia as 5:5{11mM. In the model an increase in I Kb of 120% was used to achieve an elevation of ½K z e to 6:5mM, resulting in an increase in the membrane potential of 6:6mV to {73:4mV. This is consistent with measurements in cat (5:4{8:1mV over 10 minutes) [129,130], sheep (5mV over 5 minutes) [131], guinea pigs (13:12mV over 20 minutes) [132] and mice (16mV over 10 minutes) [133], but less than measurements in rabbit (33mV over 12 minutes [10]) and guinea pig (32:5mV over 15 minutes [120] or 22{29mV over 30 minutes [134][135][136]). The broad variation in membrane potential changes was not unexpected given the range of changes in ½Na z i and ½K z e reported and in simulations an intermediate value has been achieved.
The individual and combined effects of the four components of ischemia on the temporal evolution of the maximum change in ionic concentrations and the spatial concentration and potential distribution profile after 5 minutes of simulated ischemia are shown in Fig. 3. As ischemia progressed, ½Na z i continued to increase. The early rise was attributed to NaK inhibition; the later increases were due to metabolic and respiratory acidosis. The early rise in ½K z e was significantly affected by NaK inhibition with I Kb activation playing a greater role as ischemia progresses and the membrane potential diverges from the K z reversal potential. The decrease in pH i was solely due to acidosis with no impact from I Kb activation or NaK inhibition. ½Ca 2z i elevation was contributed to principally by NaK inhibition, while increased I Kb decreased ½Ca 2z i . Ischemia caused an initial drop followed by a sustained rise in V m . I Kb activation caused an early hyperpolarisation of the membrane potential, which was subsequently countered by the depolarising effects of NaK inhibition, with acidosis having a limited effect. Figure 4 plots the width of the BZ for each ion with striped bars corresponding to extracellular space and the darker the bar the more significant the concentration gradient relative to the initial concentration. This plot shows that ½K z e had a significantly wider BZ with greater magnitude than ½Na z i . For pH regulation, H z had a narrower BZ compared to the significantly wider HCO { 3 BZ, which may indicate the facilitation of proton transport via HCO { 3 diffusion.

Extracellular Potassium Gradients
The ½K z e BZ width was significantly wider than other BZ ion widths and notably larger than the ½Na z i BZ width. To determine the cause of this extended ½K z e BZ the source of the cumulative changes in ½K z e due to transmembrane flux, drift or diffusion over the 5 minutes of ischemia were calculated and plotted in Fig. 5. This showed that only the transmembrane flux had a significant gradient across the ischemic region. Separating the transmembrane flux into the NaK and I Kb components then identified the K z channel, which includes the ATP-inactivated K z current, as the cause of this gradient. The gradient of I Kb was due to the extensive membrane potential gradient into the ischemic region (see Fig. 3). To confirm that this was the cause  of the ½K z e BZ width, the membrane potential was calculated as normal, but an additional clamped membrane potential was calculated at each time step. The clamped membrane potential had the same maximum and minimum values as the correct membrane potential but instead of a smooth gradient across the BZ it had a sharp transition over the BZ. A comparison of this clamped and the control membrane potentials is shown in the Fig. 5C. The effects of using a clamped membrane potential to calculate the NaK flux or the I Kb current on the cumulative changes in ½K z e after 5 minutes of simulated ischemia are plotted in Fig. 5B, demonstrating that by removing the gradient in the membrane potential experienced by I Kb there is a significant narrowing in the ½K z e BZ.

Drift and Diffusion
To investigate the relative contribution of drift and diffusion to intra region fluxes (inter or extracellular) the drift and diffusion fluxes in each region were plotted, alongside the transmembrane flux, over the length of the strand after 5 minutes of ischemia. Fig. 6 shows the differences in drift and diffusion between Cl { , Na z and K z . As expected from the intra and extracellular gradients, intracellular K z and Cl { diffused into and out of the ischemic region in the intra and extracellular space, respectively. The converse was the case for Na z . Due to the decrease in transmembrane potential, characteristic of ischemic regions, there was a convergence of intra and extracellular potentials with the extracellular potential decreasing in the ischemic region and the intracellular potential increasing. This gradient caused positive ions to drift into the ischemic region in the extracellular space and drift out of the ischemic region in the intracellular space. The converse was true for negatively charged ions. As drift is proportional to the ionic concentration, Na z drift was significant in the extracellular space and K z drift was significant in the intracellular space. For Na z , drift and diffusion operated in the same direction in the intra and extracellular space. The result was a cyclical movement of Na z moving into the ischemic region in the extracellular space, while moving out of the ischemic region in the intracellular space. The Cl { movement was also circular but in the opposite direction (Fig. 6P). However, for K z drift and diffusion were in opposite directions. In the extracellular space where the K z concentration was low, diffusion dominated and K z moved out of the ischemic region. In the intracellular space, where there was a higher concentration of K z , drift dominated, also causing K z ions to move out of the ischemic region. Thus ischemia caused a depletion of K z in the ischemic region through both the intra and extracellular space and, contrary to previous hypothesis [13], the model suggests that intracellular K z movement is the dominant path for K z to leave the ischemic region.

Discussion
In this study a new model of cardiac tissue electrophysiology was developed. The model predicted that the width of the ½K z e gradient across an ischemic BZ would be significantly wider the ½Na z i BZ. The cause of this difference was attributed to the voltage dependence of the I Kb channel. The model also demonstrated that, due to electrogenic drift, K z moved out of the ischemic region in both the intra and extracellular space which will lead to K z depletion.
The model of ionic movement and tissue electrophysiology was developed by combining the Nernst-Plank equations with the bidomain framework. No attempt was made to explicitly validate the proposed tissue model equations due to the paucity of experimental data. However, applying simplifying assumptions with regards to ionic or voltage gradients reduces the proposed equations to the well validated bidomain equations [137,138] or coupled reaction-diffusion equations [28,86], respectively, providing support for the validity of the proposed modelling framework. The limited attempts at simulating the spatial temporal evolution of ionic gradients across the ischemic BZ have largely uncoupled the movement of ions and the electric field. Potse et al., [18] demonstrated that measured K z gradients across an ischemic BZ could be simulated using a model of K z diffusion coupled to a source term. The spatially varying ½K z e gradient could then be included as a boundary condition to models of transmembrane current in the bidomain equations. This model did not include any effect of electric gradients on K z movement, ½K z i movement, inter ionic species interactions or the effect of the ischemic region on any other ion gradient. Similar sets of equations to those proposed here have been used for simulating the potential gradient surrounding cells, including the Debye layer [139], electrical propagation along strands of cardiac cells [140,141] and for modelling ion diffusion in the cable equation [142]. These previous models have either explicitly represented the intra and extracellular domains or only considered the intracellular domain or NaK (blue dashed line) are exposed to a clamped membrane potential. The reference BZ and BZ when I Kb is exposed to a clamped membrane potential are indicated by the yellow and purple shaded regions, respectively. doi:10.1371/journal.pone.0060323.g005 but have not modelled the tissue within the bidomain framework, as derived and implemented here.
The proposed equations can be applied generally in three dimensions as opposed to the one dimensional simulations presented here. The current study does not consider the effects of anisotropy on ion or membrane potential gradients across the BZ, however, if implemented in two or three dimensions the model is capable of representing tissue anisotropy and any effects this may have on BZ gradients. By conserving ionic species the proposed equations provide a more biophysical representation of cardiac electrophysiology than the bidomain equations and can appropriately be applied to simulate a broader range of conditions. However, these benefits come at a cost. Unlike the bidomain equations, with two partial differential equations that can readily be uncoupled and solved as two sets of linear equations [143], the proposed model is nonlinear and contains two parabolic partial differential equations for each ionic species and one elliptic partial differential equation to model the electric potential. This results in a significant increase both in the complexity and number of equations that must be solved and hence comes at a significant increase in computational cost. The proposed framework can be used to simulate electrically active tissue by introducing a full action potential cell model [31]. This would require identifying and separating out the transmembrane pathways for each ionic species present in the cell model to calculate the net transmebrane flux for each ion or J j m in Eqn 18. The internal cell model state variables, including gating variables, Markov states and intracellular Ca 2z dynamics would also need to be solved, introducing a system of nonlinear ordinary differential equations at each grid point. The small time steps and increased number of degrees of freedom required to simulate electrically active tissue would further increase the computational cost of the proposed model.
A new model of ion movement across the un-stimulated cell membrane was developed. This model has been published online and is available at cellml.org. The model creation approach demonstrated two novel methods. Firstly, the ion transporter densities were uniquely constrained by a small number of experimental data sets and a zero net flux constraint. This provided a repeatable and unique method for determining model parameters under quiescent conditions and could readily be applied to any cardiac electrophysiology model to constrain model parameters. Secondly, in the development of this model a grid of 53 experimental observations were created to provide a comprehensive validation of the model response to changes in ion concentrations and in the presence of blockers of major transporters. Although these observations only provide a qualitative comparison, they do provide a benchmark for quantifying the generality of models. Furthermore, Table 4 identified 19 experiments that do not appear to have been performed or remain controversial, highlighting the potential for additional experiments.
In this study the process of validating the model against 53 experimental observations demonstrated the general capacity of the model to replicate the majority of experimental results. However, the model was unable to match ten of the 53 observations. Five of the experimental observations found no change in a measurement. Due to the numerical nature of the model, even very small changes in a concentration can be observed and without adding in a semi arbitrary threshold that should reflect the variability and confidence of each experimental measurement it was not possible for the model to return no change in a value. In general the five remaining failed observations can be attributed to absent mechanisms in the model or inconsistencies with the model and specific experimental setups. The model did not predict the ½Cl { i response to depressed ½Na z e and the ½K z i response to depressed ½Cl { e . This is potentially due to the absence of NaK2Cl from the model, largely due to the lack of data characterising the transporter. The model was unable to replicate the decrease in V m due to elevated ½Na z e . This may be due to the specifics of the experimental protocol. Decreasing ½Na z e , which should have the opposite effect on V m , caused both increases and decreases in V m in different studies indicating that the V m dependence on ½Na z e is sensitive to the specifics of the experimental setup. The model did not predict the change in V m or H z with a decrease in ½Ca 2z e . Experimental observations report an increase in V m with both a decrease or an increase in ½Ca 2z e . This may indicate that the V m is at some minima with respect to ½Ca 2z e , although this seems unlikely. It more likely reflects inconsistent data due to differences in experimental setups that were impossible for the model to replicate. In cardiac myocytes it is postulated that Ca 2z and H z compete for common buffering sites [144,145]. The decrease in ½Ca 2z e may drain the cell of Ca 2z , reducing the Ca 2z bound to buffers that could then be occupied by H z , resulting in a decrease in ½H z i . A common pool of Ca 2z and H z buffering was not present in the model and this may explain the disparity between model predictions and experimental results. Despite the inability of the model to replicate 5 of the 53 observations, the majority of the absent mechanisms are expected to play a secondary role during ischemia. NaK2Cl is a potential contributor to the elevation in ½K z e [146], although this has been questioned [131], and its absence is unlikely to affect the general model conclusions.
Models simulations found a BZ width for the different ions of between 0:3mm (½Na z e ) and 2:4mm (½K z e ). The variation in the width of each ion across the BZ means that there was no single BZ width predicted by the model. The majority of ions transitioned from viable to ischemic concentrations over 0:4{0:8mm, with HCO { 3 and ½K z e notable outliers. The model predictions are consistent with previous measurements of the BZ between 0:6 and 2mm [10,11,[147][148][149]. However, the model did not include the sharp gradients in metabolites over 0:1mm [11] or the effects of regions of reduced perfusion affecting mechanics over 20 mm distance from the ischemic region [150].
The model predicted that ½K z e and ½Na z i have significantly different BZ widths (Fig. 4). It was expected that the coupling of ½Na z i and ½K z e via NaK would result in concordant gradients in these two ions across the BZ. The extended ½K z e gradient was attributed to the voltage dependence of the I Kb (Fig 5). Unlike ½Na z i , ½K z e is regulated by two electrogenic transmembrane pathways I Kb and NaK, whereas ½Na z i is regulated by multiple exchangers that are electro neutral, or have attenuated voltage dependence, compared to ion channels. This discordance in ½Na z i and ½K z e gradients will have a significant effect on the gradient of action potential morphology across the BZ. Elevation of ½K z e depolarises the cell and shortens the action potential duration, whereas elevation of ½Na z i causes a reduction in action potential duration [151,152]. The combined effect of the two gradients will be a slower change in resting membrane potential on the length scale of the ½K z e gradient and a much more rapid transition in action potential duration due to the faster change in ½Na z i in conjunction with the change in ½K z e over the BZ. The temporal evolution of the two gradients are also distinct, with a sustained constant increase in ½Na z i over the first 5 minutes, while ½K z e increases rapidly for 3 minutes before plateauing (Fig. 3). These simulation results are consistent with experimental measurements of ionic concentrations [111,152] and changes in ECG morphology, which report an early elevation of resting membrane potential, followed by a decrease in action potential duration [153] during early ischemia. These spatial-temporal increases in ion concentrations and secondary effects on electrophysiology will increase tissue heterogeneity and have the capacity to play an important role in the BZ arrhythmogenic substrate.
Experimental measurements have reported a decrease in ½K z i during ischemia [152]. Loss of ½K z i has been attributed to a combination of increased I Kb conduction and NaK inhibition [152,154]. However, depletion of ½K z i due to transmembrane K z movement, in the absence of a potential gradient, would cause an intracellular flux of K z into the ischemic zone, due to diffusion, to replenish ½K z i , mitigating the effects of changes of K z transmembrane flux on ½K z i . The model proposed here demonstrated that this is not necessarily the case. The model predicted a significant electrogenic K z flux in the intracellular space out of the ischemic zone (Fig. 6), resulting in a net efflux of K z out of the ischemic zone both in the intra and extracellular space. The movement of K z in the intracellular space out of the ischemic zone would further deplete ½K z i in the ischemic region and exacerbate ½K z i loss, but may limit K z accumulation in the extracellular space.

Limitations
The model is inherently an approximation and hence represents a finite set of known cellular properties and changes that occur during ischemia. In particular, the model treated all buffers as static and rapid, the effects of protons on channel, exchanger and co-transports were not considered, ischemic changes were instantaneous and the model did not include all possible changes or pathways that may affect ionic homeostasis during ischemia.
The tissue model assumed that H z and Ca 2z buffers are static, rapid and made up of a single population of binding sites. It is known that some of the buffers for both H z and Ca 2z are mobile, these could be introduced into the model framework as an additional concentration but these effects were approximated, without the additional computational cost of adding a additional ionic concentration, by using effective diffusion constants. The equilibrium assumption was likely to be valid in the current model due to the long time scales of interest; however, simulation of cardiac action potentials would require the re-evaluation of this assumption. It is also known that H z [28] and Ca 2z [155] are buffered by multiple proteins with distinct binding kinetics but over the range of concentrations simulated these multiple buffer species were unlikely to have significant effects.
The model of intracellular Ca 2z dynamics assumes that the SR and cytosolic Ca 2z concentrations remain in equilibrium, which is clearly not the case during an action potential. The model of intracellular Ca 2z dynamics provided a numerically efficient representation of intracellular Ca 2z buffering and SR Ca 2z uptake. The use of a Hill coefficient of one for SERCA as opposed to the more common and biophysical value of two removed the need for an additional differential equation to model ½Ca 2z SR or the solution of a set of nonlinear equations to model Ca 2z regulation. Furthermore, over the range of Ca 2z values studied it was possible to adjust the maximum SERCA flux to minimize discrepancies between a model with a Hill coefficient of one or two.
It is well recognised that H z play an important role in regulating cellular electrophysiology [156]. Experimental and modelling studies have demonstrated the effects of H z on ryanodine receptor opening probabilities, NCX , Ca 2z channels, Ca 2z buffering and SERCA [22,157]. The majority of the effects of H z on ½Ca 2z i regulation are unlikely to play a significant role in determining the spatial and temporal K z and Na z . However, inhibition of NCX potentially contributes to ½Na z i gradients but this will be secondary to the effects of NaK inhibition.
The model treats all changes for ischemia as instantaneous. The time dependence of inhibition of NaK, activation of I Kb , increased ½CO 2 i flux or increased H z flux are not known and can only be approximated. A linear ramp in ischemic changes has been used previously, but this fails to consider the possibility that changes occur over different time scales. In order to minimise ambiguity in model simulations an instantaneous change in NaK, I Kb , ½CO 2 i flux and H z flux was chosen.
To limit the scope of this study, the effects of cell swelling and the effects of this on changes in ion concentrations [122] were not included in the model. However, previous, modelling studies have found these effects to not significantly alter ½K z e accumulation and may not fundamentally alter the study conclusions [20]. Mitani and Shattock [146] identified the Na dependent potassium channel as a potential contributor to the elevation of ½K z e . However, the effects of an increase in any K z current are captured by the elevation in I Kb that does not necessarily need to be the sole result of an increase in conduction in the ATP inactivated K z channel. Furthermore, the model does not include an increased Na z channel conductance during ischemia. This has been reported during ischemia in the form of increased permeability of Na z through the persistent Na z current [116] and the ATP activated K z channel [158]. However, other groups have found a limited impact of the persistent Na z channel in Na z accumulation during ischemia [159] and previous modelling studies have shown that its inclusion is not required to capture the salient features of ischemia in the single cell [20]. For these reasons a potential increase in the Na z channel conductance was not included in this study.

Summary
A new mathematical framework was derived for simulating cardiac tissue electrophysiology with ion species conservation. The model was used to simulate the movement of ions due to transmembrane channels, pumps and transporters, diffusion and drift in the intra and extracellular space. The model predicted that 1) the sodium BZ is approximately a quarter of the length of the potassium BZ and this is due to the effects of the membrane potential gradient on I Kb and 2) that during ischemia there is a gross movement of potassium ions out of the ischemic region in both the intra and extracellular space due to the effects of drift, which will lead to a depletion of K z from the ischemic region.

Author Contributions
Conceived and designed the experiments: SAN. Performed the experiments: SAN. Analyzed the data: SAN. Contributed reagents/materials/ analysis tools: SAN. Wrote the paper: SAN.