Figure 2 summarizes some of the main findings corresponding to the “0D simulations”. The solid line in Figure 2A shows the time course of [ K + ] o during 30 min of simulated progressive ischemia in the virtual cardiomyocyte stimulated at 1 Hz. Data were sampled at 1 s intervals to eliminate the small “micro” fluctuations (in the μmol/L range) in [ K + ] o during the systole-diastole cycle, and highlight the“macro” behaviour of [ K + ] o . It can be seen that [ K + ] o is stable before the onset of ischemia (with a value of 5.4 mmol/L), rises during the first ≈5 min, smoothly reaches a plateau (≈12 mmol/L) that lasts until the 15th minute approximately, and then slowly increases again.

To test the validity of the simulation, results were compared to experimental data published elsewhere (Wilde and Aksnes, 1995). The dots in Figure 2A correspond to adapted experimental data of extracellular potassium accumulation obtained in rabbit hearts subject to acute global ischemia. The data shown in the figure have been shifted +0.9 mmol/L to account for the difference between the normoxic value of [ K + ] o in Wilde and Aksnes. (1995) (4.5 mmol/L) and in our simulations (5.4 mmol/L). After the 16th minute post-occlusion, the experimentally measured [ K + ] o increased faster than in our simulation (not shown). In Supplementary Figure S1, a direct comparison of our simulation with multiple experimental recordings reported in different works, corresponding to a variety of animal species and/or heart rates, is shown without any shift in the [ K + ] o values.

In order to test the soundness and sensitivity of the results to the specific time courses of intracellular ATP and intra/extra-cellular pH, we carried out alternative simulations using different time courses and/or initial normoxic values and/or final ischemic values for said ischemic parameters. The results, which can be found in Supplementary Figure S2, show that, despite of small quantitative differences, the qualitative features of the extracellular K ⁺ behavior remain unaltered. As for the effects of different possible time courses of intracellular ADP, its effects are presented and discussed below.

Finally, the effect of heart rate is shown in Figure 2B. The results show that the rate of rise of extracellular potassium increases with heart rate. In the particular case of a quiescent cell, the rise in [ K + ] o is slow and delayed, following an initial period in which a slight decrease is found.

The increase in [ K + ] o observed in the simulation and the experiments must be due to an imbalance between potassium efflux and influx that causes a net efflux during the [ K + ] o rising phase. To further analyze the behaviour of potassium fluxes, we computed the potassium flux rates using Eqs 10, 11 (see “Materials and Methods”). Figure 2C shows the time course of the different total potassium flux rates (separated into “Efflux” for the unidirectional positive fluxes through potassium channels, “Influx” for the unidirectional negative flux through the NaK pump, and “Net Efflux” for the difference).

During the normoxic period, efflux and influx are equal in magnitude [34.5 (μmol/L)/s], so there is no net potassium efflux. As soon as ischemia begins, an imbalance between the potassium efflux and influx quickly arises, with a rapid decrease in the magnitude of the influx rate being coupled to a lesser increase in the efflux rate. The resulting net efflux rise (Figure 2C) during the first ≈4 min post-occlusion provokes the first rising phase in [ K + ] o Figure 2A.

Just before the fifth minute of ischemia, potassium influx has almost stabilized [at ≈ 13 (μmol/L)/s]. However, the efflux rate suffers an abrupt drop, just prior to showing notorious oscillations [which are due to alternans in the action potential duration (APD), as explained below]. In the following 2 minutes, the efflux keeps decreasing until it almost equals the magnitude of the influx. As a result, the net potassium efflux falls to near zero values, which in turn gives rise to the plateau in [ K + ] o observed in Figure 2A. After several minutes, the efflux rate stabilizes in a lower value [at ≈ 10.5 (μmol/L)/s], while the magnitude of the influx suffers a slow but sustained decrease, provoking a continuous slow increase in the net potassium efflux (Figure 2C). This gives rise to the secondary increase in [ K + ] o observed from roughly the 10th minute onwards in Figure 2A.

To better understand the causes of the dynamic changes in potassium efflux and influx, the individual contributions of the different potassium transmembrane currents were computed using Eq. (10) at selected time instants depicted by the short vertical lines in Figure 3A. The results are depicted in panels B to F in Figure 3. The coloured bars represent efflux rates (when positive) or influx rates (when negative). The continuous data of the time course of the flux rates corresponding to the individual contributions of all the potassium currents is shown in Supplementary Figure S3.

As shown in Figure 3B, which corresponds to the normoxic situation just before the onset of ischemia, the current carrying the highest potassium efflux is I _Kr , with other currents playing a minor role. Unsurprisingly, the NaK pump is the only current that generates a noticeable potassium influx. At this time point, influx and efflux are balanced.

Two and a half minutes later (Figure 3C), the most notorious change is the decrease in the magnitude of the influx rate [from 34.5 to 15.0 (μmol/L)/s] caused by the reduction in the current transported by the NaK pump. Concomitantly, the efflux rates through I _K1 and I _K(ATP) increase considerably, counteracting the decrease in the I _Kr efflux rate, with changes in other potassium currents being negligible. Thus, the unidirectional potassium efflux increases (from 34.5 to 37.2 (μmol/L)/s). The overall result is an increase in the net efflux [from 0 to 22.2 (μmol/L)/s], which provokes the primary rise in [ K + ] o shown in Figure 3A.

Five minutes post-occlusion (Figure 3D), the efflux rates through various ion channels have considerably changed. Most noticeably, the efflux generated by I _Kr has strongly diminished and the only currents generating an appreciable potassium efflux are now I _Kr , I _K1 and I _K(ATP). Even if the influx through the NaK pump has very slightly decreased, this important reduction in the efflux provokes a strong decrease in the net efflux rate [now 11.2 (μmol/L)/s]. This is responsible for the decrease in the rate of rise in [ K + ] o observed at that time in Figure 3A, which in turn begins to generate the plateau in potassium accumulation.

Ten minutes post-occlusion (Figure 3E), when [ K + ] o has already reached the plateau, only I _K1 and I _K(ATP) generate potassium efflux, with the efflux of other potassium currents having dropped to zero. This causes a strong decrease in the total efflux compared to Panel D. Meanwhile, the NaK pump has continued decreasing slowly. At this time point, efflux and influx are almost equal [yielding a net efflux rate of 0.9 (μmol/L)/s], which is consistent with the plateau phase in extracellular potassium seen in Figure 3A.

This situation is similar after 20 min of ischemia (Figure 3F), except that the NaK pump influx is further reduced, again generating an imbalance between influx and efflux that gives rise to a net efflux of 4.5 (μmol/L)/s and a secondary slow rise in [ K + ] o shown in Figure 3A.

To further investigate the causes of these changes in efflux and influx as ischemia develops, we plotted the time course of APs and selected potassium currents at the time instants discussed above. The top row in Figure 4 shows how APs change as ischemia progresses. Resting membrane potential becomes increasingly less negative (a direct consequence of the increase in [ K + ] o ), APD decreases (due to hyperkalemia and also to the partial activation of the ATP-sensitive potassium channels), and the cell has become almost non-excitable by the 10th minute. Moreover, AP alternans exist in the fifth minute post-occlusion.

The second row shows the progressive evolution of I _K1. Basically, its diastolic value keeps increasing as ischemia develops due to the hump in its current-voltage relationship that implies a higher current at membrane potentials moderately higher than the potassium Nernst potential (E _K ). This causes the increase in its associated efflux seen in Figures 3D–F.

The third row represents the evolution of I _K(ATP). While it is almost non-existent in normoxia, its amplitude keeps increasing in ischemia (due to the decrease in the intracellular ATP/ADP ratio). Although its duration decreases with time (due to the reduction in APD), the area below the curve (proportional to the efflux generated by the current) is almost constant.

The fourth row shows the progression of I _Kr (which generates the highest potassium efflux in normoxia) with ischemia. Its duration has decreased after 2.5 min of ischemia, thus reducing the efflux, and it has almost vanished 10 min after the onset of ischemia (due to the cell being almost non-excitable at this point). Five minutes post-occlusion, AP alternans cause an alternation in the amplitude and duration of I _Kr , which is very low in the short AP compared to the long one. Further effects of AP alternans in the time course of [ K + ] o will be discussed below.

Finally, the lower row corresponds to the NaK pump current. The decrease in intracellular ATP levels continuously reduces the magnitude of the current, thus reducing the influx of potassium as was apparent in Figures 3B–F.

Once we identified the causes of the triphasic nature of extracellular potassium accumulation in regards to potassium fluxes, we further analyzed the influence of the dynamically changing AP morphology during ischemia on the time course of [ K + ] o .

The data depicted in red in the lower graph in Figure 5A represents the time course of the APD at 90% repolarization (APD ₉₀) during the 30 min period of simulated ischemia at the same time scale as [ K + ] o (upper blue graph). Immediately after the onset of ischemia, APD ₉₀ begins to decline, as expected. To test the validity of these results, the simulated values were compared to experimental results published elsewhere. The four green dots in the lower graph of Figure 5A show Monophasic Action Potential (MAP) duration values (adapted by adding 32 m to account for the basal difference with our APD ₉₀ values) measured in acutely ischemic human hearts in-vivo (Sutton et al., 2000).

As shown in the lower graph of Panel A, between the fourth and the sixth minute of ischemia (approximately), APD ₉₀ values begin to oscillate, indicating that the cell is developing electrical alternans. This alternans phase is delimited by the dashed vertical lines. Figure 5B shows the detail of the waveform of two consecutive alternating APs obtained in our simulations (at the 4.6 min mark). At the beginning of this phase, alternans with a 2:2 pattern (i.e., long–short–ong–short) are predominant, but shortly afterwards more complex patterns arise (e.g., long–short–short–long–long–short–long).

As the dashed vertical lines indicate, the period of occurrence of alternans matches the smooth transition from the primary rising phase of [ K + ] o to the plateau, suggesting that the alternans themselves may be responsible for the gradual deceleration of the extracellular potassium rise and the genesis of the plateau. To confirm this hypothesis, APs at selected time intervals spanning four consecutive APs were plotted together with [ K + ] o . At these short time scales, [ K + ] o clearly shows the “micro” fluctuations caused by the systole-diastole cycles (inset in Panel A top graph), as clearly seen in the lower (blue) plots of Panels C through F.

During normoxia (Figure 5C), APs do not alternate at all, and potassium currents are balanced in a way that the diastolic minimum potassium levels (DmKL, indicated by the * signs) reach the same value at the end of each diastolic period. This is consistent with the normoxic stable behaviour of the “macro” [ K + ] o with zero net efflux shown in Figures 2C, Figure 3B.

On the contrary, at the 2.5 min mark in the ischemic cell (Figure 5D), the positive net potassium efflux seen in Figures 2C, Figure 3C is related to a continuous increase in the DmKL from one AP to the following one. Consequently, the “macro” [ K + ] o rapidly increases. However, this trend is attenuated when alternans begin to appear. Five minutes post-occlusion (Figure 5E), when the “macro” [ K + ] o is decelerating (Figure 5A), the difference between consecutive DmKL values diminishes in relation to Panel D. The reason is that the “short” APs originate lower positive shifts in the “micro” [ K + ] o fluctuations due to the fact that I _Kr , the strongest efflux generator at the 2.5 min mark, almost vanishes in the “short” APs (see Figure 4). This is because a) it is a systolic current that is only active during the AP, which is now short, and b) because the driving force for potassium ions is smaller due to the lesser peak of the short AP. Therefore, the average potassium net efflux rate diminishes (Figures 2C, 3D), causing a reduction in the slope of the “macro” [ K + ] o that eventually causes the plateau to appear (upper Panel A).

Finally, 10 min after the onset of ischemia (Figure 5F), only short APs are present in the cell due to its partial loss of excitability provoked by an already elevated value of [ K + ] o (see Panel A). In these circumstances, the low positive shifts in the “micro” [ K + ] o fluctuations are almost equal for consecutive beats, thus the “macro” [ K + ] o remains almost constant (as seen in the upper graph of Panel A).

Figure 6A shows the spatial profiles of [ K + ] o in the tissue preparation from just before and up to 30 min after the onset of ischemia. This profile increases monotonically in the altered tissue for the first 2 minutes post-occlusion, to then adopt a non-monotonic profile with the spatial peak increasing with time and moving toward the center of the altered tissue. These changes in [ K + ] o also affect a small region in the distal part of the “normoxic” tissue. The position of the peak value of [ K + ] o stabilizes after 15 min of the onset of ischemia, defining an ischemic BZ of approximately 1.2 cm, in agreement with experimental results (Hill and Gettes, 1980; Coronel et al., 1991). The slope of the rising phase of the [ K + ] o keeps increasing until minute five post-occlusion before starting to decrease, while the spatial peak of the [ K + ] o profile moves toward the center of the ischemic zone. Around minute 15, the rising phase of the [ K + ] o profile has become almost linear, and the location of the peak has stabilized. From this time onwards, the peak keeps rising together with the concentration in the CIZ until almost reaching a plateau by minute 30. A video with the evolution of the [ K + ] o profile with time is available in the Supplementary Movie SM1.

To test the sensitivity of the results to the time-course of [ADP] _i and the spatial extension of the pH transition zone, simulations considering the ADP time-course proposed by Weiss et al. (1992) were performed for three different pH transition zone sizes (0 cm, 0.5 cm, and 1 cm). Results from the simulations are shown in Supplementary Figure S4. The results show that the ADP data from Weiss et al. (1992) leads to a lower rate of rise of [ K + ] o during the primary rising phase, which delays the potassium plateau with respect to our control simulations based on the ADP data from Terkildsen et al. (2007). However, the results are qualitatively similar. Indeed, a plateau is reached when AP alternans flatten the curve, the same as in our main simulations (compare left and right curves in Supplementary Figure S4B, C). As for the [ K + ] o spatial profiles, they are affected by the time-course of ADP (the size of the [ K + ] o border zone in particular), with changes in the order of a 10%, with the pH transition zone affecting less. However, the time evolution described in the previous paragraphs is observed to be the same for both of the considered ADP time-courses. Furthermore, simulations performed with a mesh size of 0.125 mm showed no difference in the results with respect to those obtained with a mesh size of 0.25 mm (see Supplementary Figure S5).

Furthermore, to better understand the mechanisms behind the formation of the BZ, a sensitivity analysis of the effect of the tissue conductivity and the potassium diffusion coefficient was performed. Together with the control simulation, the following combination of parameters were considered: i) The tissue conductivity, D _V , was reduced by 20%, by 60%, and by 80% down to 0.0005 m (conduction velocity of 30 cm/s) while keeping the extracellular potassium diffusion coefficient, D _K , at the control value in Table 1; ii) the extracellular potassium diffusion coefficient, D _K , was multiplied by 0 (no extracellular potassium transport i.e., D _K = 0.0), by 10, and by 100 while keeping the effective tissue conductivity to the value in Table 1. Results of these simulations are show in Supplementary Figure S6, whereas Figure 6B shows only results for four representative cases. Results show that increasing D _K reduces the value of [ K + ] o in the central ischemic zone and enlarges the BZ by increasing the region in the distal part of the normoxic tissue affected by ischemia. This leads to a smoother transition between the normoxic and the CIZ. Conversely, reducing D _V makes the BZ smaller, reduces the maximum [ K + ] o , and leads to a biphasic [ K + ] o spatial profile. The reduction in the maximum [ K + ] o and the biphasic behavior in the spatial profile are caused by an early appearance of AP alternans in the border zone and conduction block in the CIZ as D _V decreases. The early appearance of AP alternans cause a premature ending of the rising phase, whereas the absence of systolic current in the CIZ after conduction block causese a reduction in the extracellular K ⁺ efflux responsible for the biphasic behavior in the spatial profile. Most importantly, results indicate a larger influence of the electrotonic coupling compared to extracellular potassium transport in the formation of the BZ. Indeed, the figure shows that a 60% reduction in tissue conductivity implies a 40% reduction in the size of the BZ (from 1.2 cm to 0.72 cm), whereas a 10 fold increases in D _K only increased the size of the BZ in a 10% (from 1.2 cm to 1.35 cm). In addition, the total suppression of extracellular potassium transport left the BZ almost unaltered (see Figure 6B).

To better understand the role of extracellular potassium transport in the [ K + ] o accumulation, the efflux, influx and net efflux diffusion rates as a function of time are depicted in Figure 7 for two positions in the tissue, namely x = 2.5 cm in the BZ, and x = 3.5 cm in the CIZ. The figure shows that the net diffusion efflux rate is four orders of magnitude smaller than the net efflux associated with the various ion channels shown in Figure 3. The figure also shows that net diffusion efflux rate results non-zero when the [ K + ] o experiences a maximum (see Figures 7A, B), or during the transition from the first rising phase to the second rising phase (see Figures 7C, D). However, even in this case, the net diffusion efflux is much smaller than the transmembrane efflux. From the analysis of Figure 7 it emerges that the evolution in the spatial [ K + ] o profile is associated with a heterogeneous local potassium accumulation in the altered region.

Figure 8 shows the [ K + ] o time profile at different locations within the tissue, together with the evolution of the APs and the electrograms during the first 16 min of ischemia. The time profile of [ K + ] o (Panel A) changes according to the position within the tissue. While in the unaltered region the [ K + ] o remains constant, as expected, in the BZ (position x = 2 cm) the time profile shows a non-monotonic behavior, in good agreement with the experimental results reported by Coronel et al. (1988) in regionally ischemic pig hearts. As moving from the BZ towards the center of the ischemic zone (positions x = 2.9 cm and x = 3.3 cm), the [ K + ] o time profile evolves toward the characteristic triphasic shape shown in Figure 2A for the isolated cell, and in good agreement with experimental data reported by Hill and Gettes (1980). While the local triphasic behavior of [ K + ] o is associated with the ionic mechanisms described in the previous sections, the heterogeneity in the [ K + ] o time profile within the tissue appears to be mainly related to the electrotonic coupling as demonstrated in Figure 6B.

The heterogeneity observed in the [ K + ] o time profile is also observed in the APD (Figure 8B) and in the APs themselves (Figure 8C). Within the CIZ (x = 3.3 cm) the AP shows changes in morphology associated with the alterations due to ischemia, with alternans developing after the 2.5 min mark, followed by total loss of excitability after the 5.5 min mark. Concomitant changes are observed in the proximal side of the CIZ (x = 2.9 cm). Here, a certain degree of alternation is also observed, though the cells remain electrically responsive during the whole simulation period. At the proximal side of the BZ (x = 2.0 cm), AP alternans are hardly noticed and the only changes in the AP morphology happen due to alterations associated with ischemia.

These electrophysiologic changes are clearly reflected in the electrograms depicted in Figure 8D for different positions in the tissue. The electrograms show the typical depression in the ST segment due to ischemia in the altered region (small black arrows in Figure 8D), whereas the alternans in the AP are reflected in alternans in the T-wave (see Figure 8D right panel). The conduction block in the CIZ at minute 13 is evident by the total suppression in of the T-wave in all electrograms, and the monophasic electrogram at the CIZ (x = 3.3 cm).

The hypothesis that during acute ischemia, a current flowing between the ischemic and the normal tissue, a current of injury, may cause depolarization of the normal cells near the border zone has been proposed in several studies (Janse et al., 1980; Coronel et al., 1991). To explore this hypothesis more in detail by means of numerical simulations, the injury current in the tissue at 3 min and 4 min after the onset of ischemia were computed. The upper panel in Figure 9 show the AP at three position in the tissue located at the normal zone (x = 0.9 mm), the border zone (x = 2.3 mm), and the central ischemic zone (x = 3.2 mm), whereas panel B shows the electrograms at the same positions. The instant t _A has been chosen in correspondence with the end of repolarization of the cell in the normal zone, coincident with the point of deep negative T wave in the electrogram recorded in the ischemic zone. On the contrary, t _B corresponds to an instant of complete repolarization in the three positions. Panel C shows the depolarization current density for the stimulations at 3 and 4 min after the onset of ischemia, at the instant when the propagation front is entering the altered region (x = 2 cm). Panel D depicts the intracellular injury currents at instants t _A and t _B at 3 min and 4 min after the onset of ischemia.

The results in Figure 9 show that the magnitude of the depolarization current, required to generate an ectopic beat, is about eight times larger than the magnitude of the injury current flowing from the ischemic to the normal zone.

Due to the difficulties related to the direct measurement of the individual contributions to the rapid primary rise of [ K + ] o , the classical hypotheses that aim an explanation are basically of speculative nature. In summary, the following factors have mainly been hypothesized to play a role: a) a partial inhibition of the NaK pump, b) an enhancement of the I _K(ATP) current, c) an increase of sodium influx and d) a reduction of the interstitial volume. Also, the discussion as to whether the extent to which this rise is due to an increased efflux and/or a reduced influx is still not resolved (Wilde and Aksnes, 1995; Weiss et al., 2017).

Regarding the partial inhibition of the NaK pump, due to the decrease in the intracellular ATP/ADP ratio, this effect is included in our model. Indeed, the last row in Figure 4 shows the decrease in magnitude of the current carried by the pump. Moreover, the bars in Figure 3 (and also the results shown in Supplementary Figure S3) show that the magnitude of the influx rate through the pump rapidly decreases during the primary rising phase of [ K + ] o , with the decrease in the magnitude of the potassium influx rate seen in Figure 2C being the direct consequence of this partial inhibition. Indeed, our results clearly show that this reduction is pivotal in explaining the primary rise in [ K + ] o , because the reduction in the influx (55% in 2.5 min), almost exclusively generated by the NaK pump depression, is much higher than the increase in the efflux during this phase (10% in 2.5 min). This confirms the experimentally derived hypothesis regarding the importance of the NaK pump inhibition. Also, this is in accordance with the main conclusion of a previous modelling study (Terkildsen et al., 2007). According to their model, the inhibition of the NaK pump is the major driving force for potassium loss. This is difficult to correlate with experimental findings because quantification of the influx rate of the NaK pump relies on indirect methods (Wilde and Aksnes, 1995). Yet, evidence of an ischemia-induced depression of the pump exists (Bersohn et al., 1982; Mitani and Shattock, 1992).

Regarding the enhancement of the ATP-sensitive potassium current, Our results show that, indeed, this factor is not as essential as the depression of the NaK pump in the primary rise of [ K + ] o . On the one hand, the efflux rate increase related to the enhancement of I _K(ATP) 5 min post-occlusion is only 5.8 (μmol/L)/s (see Figure 3D ;Supplementary Figure S3), which is much smaller than the concomitant influx rate decrease (20 (μmol/L)/s). In that same time period, the increase in the efflux rate provoked by the increase in I _K1, for instance, is 5.6 (μmol/L)/s, almost equal to the I _K(ATP) case. The lack of a higher increase in the I _K(ATP)-related efflux seems to be paradoxical with the continuous increase in the peak value of the current shown in the third row of Figure 4. The explanation, also pointed out by Wilde and Aksnes. (1995); Weiss et al. (2017), is that the I _K(ATP) increase is a self-limiting mechanism in terms of potassium efflux, because the continue increase in its peak is counteracted by the decrease in its duration caused by the continuous reduction in APD ₉₀. Thus, our results show that I _K(ATP) does play a role in the primary potassium accumulation phase, but it is quantitatively less important than other currents (mainly the NaK pump).

Regarding the modification in sodium influx and its role in the primary phase of potassium accumulation, the classical hypothesis is that the increase in potassium efflux is a passive consequence of an increased sodium influx to maintain electroneutrality (Wilde and Aksnes, 1995; Shivkumar et al., 1997). Further, there is evidence that [ N a + ] i increases during acute ischemia (see Wilde and Aksnes (1995) for a review), although the exact extent is unknown. [ N a + ] i does indeed increase in our simulations (from 7.3 mmol/L in normoxia to 9.3 mmol/L in 5 min, as seen in Supplementary Figure S7). Although we did not analyze the causes of this increase, the partial inhibition of the NaK pump seems to be the main factor responsible. This increase in intracellular sodium is necessarily coupled to the increase in extracellular potassium (Wilde and Aksnes, 1995; Shivkumar et al., 1997). According to our results, however, the change of one of them is not the driving cause of the other, but rather they are both caused by their common underlying mechanism: the reduction of the NaK pump itself. As pointed out in Terkildsen et al. (2007), the partial inhibition of the pump accounts for both concentration changes.

Finally, a shrinkage of the interstitial extracellular space provoked by the ischemia-induced alteration of osmotic pressure has also been hypothesized to play a role in the increase in [ K + ] o (Wilde and Aksnes, 1995). We preferred not to include this mechanism in our main set of simulations due to its intrinsically different nature compared to the efflux-influx interplay. However, we conducted separate simulations to assess the quantitative importance of the extracellular volume regulation. The main result is shown in Supplementary Figure S8, where we compare the basic “control” simulation of Figure 2A with a situation in which extracellular volume was reduced by 9% in 30 min (data extrapolated from Terkildsen et al. (2007)). The curves show that the effect of extracellular shrinkage was almost negligible during the primary [ K + ] o rising phase, and its effects were only notorious (but very limited) in the long term. Indeed, the contribution of the extracellular volume reduction was only +0.6 mmol/L after 20 min. These results are quantitatively similar to a previous modeling study (Terkildsen et al., 2007) and experimental results (Yan et al., 1996).

One of the most intriguing and novel result from our simulations refers to the causes of the smooth transition to the extracellular potassium plateau and its relationship with AP alternans. As explained in the Results section, the model predicts the appearance of electrical alternans between the fourth and the sixth minute of ischemia in the isolated cell model (Figures 5B,E). Moreover, the tissue simulations also show AP alternans in the distal altered zone (Figure 8B). This is in correspondence with experimental results in animals. Indeed, these type of alternans, occurring at normal/low stimulation frequencies (1 Hz in our case) around 5 min post-occlusion (or even earlier), are a well-known feature of acute ischemia and have been observed in AP or MAP recordings and/or in electrograms (Kléber et al., 1978; Nakashima et al., 1978; Janse, 1990; Martišienė et al., 2015; Watanabe and Gettes, 2016). Moreover, the experimentally-obtained waveforms of the alternating APs are very similar to our simulated results shown in Figure 5B. The “long” AP has the upstroke divided in two phases, whereas the “short” one lacks the secondary phase. These exact features were found in the APs directly recorded in pig hearts 5 min after clamping the LAD coronary artery in pig hearts [see Figure 2 in Kléber et al. (1978); Figure 4 in Janse (1990)]. The morphology of the AP alternans seen in Figure 1 in Downar et al. (1977) are very similar to those found in our tissue simulations (Figure 8).

In our “0D simulations” (isolated cell), the onset of the [ K + ] o plateau coincides with the alternans phase (Figure 5E). Similarly, in our “1D simulations” (heterogeneous tissue), the alternating period is also simultaneous with (x = 3.3 cm in Figure 8) or precedes (x = 2.0 cm and x = 2.9 cm in Figure 8) the onset of the plateau. In view of this result, we hypothesize that AP alternans are the main cause of the stabilization in the [ K + ] o rise. To confirm this hypothesis, we conducted a separate simulation in which the alternans were inhibited by artificially enhancing the I _CaL current at critical instants of the simulation in which the alternans were just about to appear. The results, shown in Supplementary Figure S9, indicate that the absence of AP alternans abolishes the potassium plateau.

According to our simulations, the flattening of the [ K + ] o curve occurs mainly because of a drastic reduction in potassium efflux rate (see Figure 2C) caused by the AP alternans. Comparing Panels C and D from Figure 3, indicates that this reduction is mostly due to the depression of the efflux carried by I _Kr and, to a lesser extent, by I _Kb . These currents are active during the AP plateau (i.e., during systole). When alternans appear, the duration of the plateau is reduced at least in one of every two APs, which provokes a reduction in the potassium carried by said currents. This reduces potassium efflux and decelerates the “macro” [ K + ] o curve (Figure 5A) due to the behaviour of the “micro” [ K + ] o fluctuations (Figure 5E).

Different experimental results have been published in which APD and [ K + ] o (or potassium activity) are measured simultaneously during no-flow acute ischemia (Weiss and Shine, 1982b; Kléber, 1983; Weiss and Shine, 1986; Weiss et al., 1992). In one case (Weiss et al., 1992), the potassium plateau was not reached during the duration of the experiment. In another (Weiss and Shine, 1982b), the existence of electrical alternans cannot be inferred (see their Figure 2). This could be due to the low resolution with which they measure APD (one measurement per minute), which makes it impossible to detect short periods of AP alternation. However, in one publication (Kléber, 1983), APD is measured with enough resolution as to detect alternans. In their Figure 2A, alternans are found between 4.5 and 5.5 min post-occlusion, just as in our case. The alternans occurred before the onset of the extracellular potassium plateau. In our “1D simulations” (heterogeneous tissue), similar results have been found. Moreover, Watanabe and Gettes. (2016) found AP alternans in pig hearts subject to regional ischemia before and during the plateau phase (see their Figure 5).

The mechanisms that give rise to the [ K + ] o plateau phase are not yet clarified. According to Weiss and Shine. (1986), the reduction of APD as ischemia progresses contributes to the existence of the plateau. They found that the onset of the plateau occurred when APD had reduced to about two-thirds. In our “0D simulations”, that was the case when the curve begins to flatten just when the alternans began to appear: the APD ₉₀ value transitioned from its normoxic value of 289 milliseconds–192 milliseconds just before the onset of alternans, which corresponds to a 66% percentage reduction of the control value. They hypothesized that the decrease in APD reduced the potassium efflux rate until it equaled the influx rate, giving rise to the plateau. This is consistent with our results (Figures 2C, Figure 5).

Indeed, in our “0D simulations”, the net potassium efflux rate is almost zero in the plateau around the 10 min mark. It stays below 2.2 (μmol/L)/s (10% of its peak value) from the 6.8 to the 13.4 min mark. During this period, influx and efflux rates are almost balanced with a value of ≈10 (μmol/L)/s, which is much smaller than in normoxia (≈34.5 (μmol/L)/s). On the one hand, influx rate is smaller than in normoxia because of the progressive partial inhibition of the NaK pump. On the other hand, efflux rate is also smaller due to the short duration of the AP plateau, which practically eliminates all potassium currents except for I _K1 and I _K(ATP). The former is not AP plateau-dependent (due to its inward rectification nature), while the latter has a high peak due to the low value of the ATP/ADP ratio at this point in time. The overall result is a new stable [ K + ] o at ≈ 11.7 mmol/L, the value up to which the initial 6.8 min lifted extracellular potassium. This can be better understood when comparing the “micro” [ K + ] o fluctuations in Panels C and F from Figure 5: while both correspond to stable situations (normoxia in Panel C, potassium plateau in panel F), the smaller fluctuations in the latter case are the consequence of smaller (but approximately equal) efflux and influx compared to the normoxic case.

A further proof of the influence on APD on the plateau phase is given in Supplementary Figure S9. As explained above, artificial interventions were made to eliminate AP alternans. The lack of a plateau phase in this separate simulation is because such interventions did not only abolish alternans but also maintained APD ₉₀ levels high compared to the control simulation.

Another classical explanation of the potassium plateau, different to considerations on APD, is related to the behaviour of the NaK pump. This has been argued both from experimental observations (see Wilde and Aksnes. (1995) for a review) and from previous computational modeling studies (Rodríguez et al., 2002; Terkildsen et al., 2007). Simulations carried out by Terkildsen et al. (2007) showed that the time course of I _NaK is triphasic and, according to them, the [ K + ] o plateau is a reflection of the plateau in the activity of the NaK pump, which, in turn, is due to the biphasic nature of the time course of intracellular ADP. In our simulations, we also found a deceleration in the decrease of the magnitude of the NaK-related potassium influx (see the bars in Figure 3 an also Supplementary Figure S3) which was concomitant to the potassium plateau, but the elimination of the alternans (Supplementary Figure S9), while maintaining the NaK pump plateau, abolished the potassium plateau. It must be noted that the model used by Terkildsen et al. is different from ours, as it corresponds to the previous generation of AP models (Luo and Rudy, 1994) and is not as comprehensiveness as the model used in this study (O’Hara et al., 2011). Also, their model of ischemia is also less comprehensive. Using a similar AP and ischemia model, Rodríguez et al. (2002) explained the plateau by the combined effect of the above discussed self-limiting effect of the activation of I _K(ATP) and an enhancement of the NaK pump which was not found with our model. In the model used by Rodriguez et al., this enhancement was due to the regulatory effect exerted by the increase in both [ K + ] o and [ N a + ] i on the pump activity. Although this regulatory effect was included in our formulation of the NaK pump following Cortassa et al. (2006); Terkildsen et al. (2007), its effect was more modest than in the former model. Our model is consistent with the fact that the increase produced by [ K + ] o elevation between 4 and 12 mmol/L is mild (Gadsby, 1980; Glitsch et al., 1981; Weiss and Shine, 1986) and surpassed by the opposite effect caused by the decrease in the intracellular ATP/ADP ratio.

The prior hypotheses regarding the causes of the secondary slower rise in potassium accumulation are scarce. Classically, this phase has been related to cell irreversible damage (Weiss and Shine, 1982b). It is not completely clear if the cause is related to an increase in potassium efflux provoked by the loss of sarcolemmal integrity or by a further depletion of NaK pump activity (Weiss and Shine, 1986). According to a previous modeling study, the latter could be the main cause (Terkildsen et al., 2007). From our results, a complementary hypothesis may be established. Cell membrane damage was not included in our model, so our results discard the relevance of this factor. According to our simulations, potassium efflux rate almost plateaus during the secondary rise in [ K + ] o while the magnitude of the influx keeps decreasing, giving rise to a secondary increase in the net potassium efflux rate (see Figure 2C). The fact that the efflux is maintained is mainly due to plateauing of both I _K1 and I _K(ATP) (see the bars in Figure 3, the ionic currents in Figure 4 and also Supplementary Figure S3). As for the secondary decrease in the influx rate, it is due to the secondary depression of the NaK pump (see the same figures), which is consistent with the fact that intracellular ATP levels keep decreasing, though at a smaller rate compared to the primary phase (Figure 1).

This phenomenon may be alternatively explained from a different perspective. The rate of rise in the secondary [ K + ] o increase (from ≈14 min onwards, see Figure 2A) is in the same order of magnitude as the primary rise when the cell is quiescent (lower curve in Figure 2B). This can also be seen in experimental results (e.g., Wilde and Aksnes. (1995); Figure 1). The reason is that both situations are very similar in terms of the electrical activity. On the one hand, in the quiescent case, there are no APs, so the efflux and influx completely depends on the diastolic currents. On the other hand, in the secondary rise when pacing, the isolated cell (Figure 4, last row at the 20 min mark) and the ischemic tissue (Figure 8, last row at the 13 min mark) are unresponsive, having lost excitability due to the high [ K + ] o value.

The potassium gradient that is established during acute regional ischemia near the metabolic border, the BZ, is believed to be related with the extracellular potassium transport (Coronel et al., 1988; Coronel et al., 1989). However, our simulations indicate the size of the BZ to be related to the electrotonic coupling between cells close to the metabolic border, with a minor role of the extracellular potassium transport in modulating the BZ (see Figure 6B). Our results indicate that, an effective modulation of the BZ by potassium transport requires a ten-fold increase in the potassium diffusion coefficient (see Figure 6B). These findings are in agreement with the work from Potse et al. (2007b) who conclude that the mechanism that leads to a smooth distribution of [ K + ] o cannot be physical diffusion alone. They suggest pulsating flow in the arterial and venous bed as possibly responsible for the increase in the potassium diffusion coefficient, a hypothesis also proposed by Janse and Wit. (1989). However, it is retained very unlikely that pulsating flow can increase the effective diffusion coefficient ten times with respect the measured diffusion constant of potassium in water, as to be able to sustain the smooth distribution of extracellular potassium observed in acute ischemia. However, the results in Figures 7B, D, indicate a transport of potassium from the ischemic toward the normal myocardium where the excess of potassium is removed by the washout. Our results clearly evidenced the need of systolic currents i.e., a propagating front in the tissue, for the development of the BZ. In fact, considering only diastolic currents as in the work by Niederer. (2013) leads to a BZ of only few millimeters wide and a limited rise of the extracellular potassium. Considering a BZ equally ischemic, the mechanism revealed by our model says that electrotonic coupling modifies the net potassium efflux in the BZ by modulating the current through ionic channels responsible for potassium efflux, with cellular potassium loss being similar in closely adjacent cells. This explains the small invasion of ischemic conditions within the healthy metabolic area, limited by the presence of potassium washout in the unaltered tissue. Moving into the ischemic zone, the effect of the metabolic border reduces until reaching a situation like what is seen in an isolated cell. In this regard, enhancing potassium diffusion makes this transition smoother as shown in Figure 6B.