Nano-scale solution of the Poisson-Nernst-Planck (PNP) equations in a fraction of two neighboring cells reveals the magnitude of intercellular electrochemical waves

The basic building blocks of the electrophysiology of cardiomyocytes are ion channels integrated in the cell membranes. Close to the ion channels there are very strong electrical and chemical gradients. However, these gradients extend for only a few nano-meters and are therefore commonly ignored in mathematical models. The full complexity of the dynamics is modelled by the Poisson-Nernst-Planck (PNP) equations but these equations must be solved using temporal and spatial scales of nano-seconds and nano-meters. Here we report solutions of the PNP equations in a fraction of two abuttal cells separated by a tiny extracellular space. We show that when only the potassium channels of the two cells are open, a stationary solution is reached with the well-known Debye layer close to the membranes. When the sodium channels of one of the cells are opened, a very strong and brief electrochemical wave emanates from the channels. If the extracellular space is sufficiently small and the number of sodium channels is sufficiently high, the wave extends all the way over to the neighboring cell and may therefore explain cardiac conduction even at very low levels of gap junctional coupling.


S1 Numerical scheme
In this section, we describe the numerical scheme used to solve the PNP model equations in three dimensions.

S1.1 Temporal splitting scheme
For every timestep n, we assume that the concentrations c n−1 k (and consequently ρ n−1 ) are known for t n−1 , and we solve the PNP system in two steps: Step 1: Solve ∇ h · (ε∇ h φ n ) = −ρ n−1 (1) Step 2: Solve where ∇ h refers to a finite difference discretization of ∇.
We consider the solution in a 3D grid of N x × N y × N z points, at N x positions in the x-direction, N y positions in the y-direction and N z positions in the z-direction. The distance between the points in the x-direction are denoted by ∆x 1 , ..., ∆x Nx−1 , and the distance between the points in the y-and z-directions are similarly denoted by ∆y 1 , ..., ∆y Ny−1 and ∆z 1 , ..., ∆z Nz−1 , respectively.
Here, as an example, ε i−1/2,j,q is the value of ε evaluated in the point

S1.3 Finite difference scheme for Step 2
The finite difference scheme follows the same structure for the different ionic species, k. In order to avoid confusion related to the index k for the different ionic species, we therefore describe the scheme for an arbitrary concentration c, with diffusion coefficient D and β-value β, where c, D and β can be replaced by c k , D k and β k for any of the ionic species, k. The scheme reads Written out more completely, and using approximations of the type c i+1/2,j,q ≈ 0.5(c i+1,j,q + c i,j,q ), Here, like for ε i−1/2,j,q , D i−1/2,j,q is the value of D evaluated in the point S2 Temporal development from electroneutrality to the resting state In Figure A, we show how the potential, charge density and ion concentrations close to the membrane change with time in the simulation used to find the resting state of the system for L e = 10 nm. Initially, the potential, φ, and the charge density, ρ are zero, but as the K + channels open, the intracellular potential approaches a value of about −80 mV, and the charge density approaches about −0.2 C/cm 3 on the intracellular side of the membrane and about 0.2 C/cm 3 on the extracellular side. with open K + channels used to find the resting state of the system. We plot the potential, charge density and ion concentrations in the first points outside of the membrane in the x-direction, on the intracellular and extracellular sides. In the y-and z-directions, the plotted points are located in the center of the Na + channel cluster, which is closed and thus acts as a normal membrane. We use L y = L z = 300 nm, L e = 10 nm and a K + channel cluster consisting of 36 channels. In order to speed up the dynamics, we have used a scaling factor of d K + =1 for the K + channel diffusion instead of the default value of d K + =0.5.  Figure B: Illustration of the complex profile of ρ in the K + channel clusters at rest. In the left panel, we show −ρ 0 and F [K + ] along a line in the center of the K + channel cluster of the pre-junctional cell at the end of the resting state simulation for L e = 10 nm. We observe some small deviations between the two curves, giving rise to a non-zero ρ = ρ 0 + F [K + ], illustrated in the right panel. The x-axis is shifted so that the K + channel cluster starts (from the intracellular side) at x = 0.
S3 Profile of ρ in the K + channel cluster at rest In Figure 7 in the main paper, we observed that the charge density, ρ, inside of the K + channel cluster had a rather complex profile. In this section, we wish to illustrate how this profile arises as a consequence of small deviations at steady-state from the completely linear profile defined as initial condition for the K + concentration in the K + channels. Recall from (4) in the main paper, that the charge density, ρ, is defined as where F is Faraday's constant, z k is the charge of the different ion species present, and c k is the concentration of each of these ion species. Furthermore, ρ 0 is defined such that ρ = 0 in the entire domain at the beginning of the simulation. In the K + channel cluster, the concentrations of all ion species except for K + are initially set to zero and remain zero in the entire simulation, because the diffusion coefficient of all other ions that K + is set to zero in the K + channels. This means that in the K + channel cluster, ρ is given by where [K + ] is the concentration of K + ions. In the left panel of Figure B, we plot −ρ 0 and F [K + ]. We observe that the two terms are very similar, making ρ close to zero. However, there are some small differences between the terms at steady-state, caused by small changes in the K + concentration occurred during the simulation, and this gives rise to the profile of ρ (see the right panel).

S5 Different choice of ρ 0 in the channels
In order to investigate how the choice of the initial conditions and ρ 0 in the channel clusters affect the results, we conducted simulations with a different choice than the default linear transition used in the simulations reported in the main paper. Recall here that ρ 0 is set up such that ρ is zero in the entire domain at t = 0. Moreover, the concentration of all ions are set to zero in the channel clusters except for the ion that is able to move through the channel cluster. Therefore, the profile of ρ 0 in the channel clusters follows directly from the initial conditions in the channel clusters. Figures E-M show the results of simulations using a piecewise constant concentration as initial conditions in the channel clusters. More specifically, the initial concentration is set equal to the intracellular concentration in the half of the channel channel that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space. The results of the simulations using this set up seem to be quite close to the results obtained using the linear profiles used in the simulations shown in the main paper, except for the profile of φ, ρ and the ion concentrations in and close to the K + channel clusters at rest (see Figures  E and F). In addition, the maximum changes in the extracellular potential 7 Figure C: The Ca 2+ concentration in the extracellular space between two cells in simulations with an open Na + channel clusters on the membrane of the pre-junctional cell. The width of the extracellular space, L e , and the size of the Na + channel cluster is varied in the columns and rows of the figure, respectively. The plots show the solution in the (x, y)-plane for the center of the domain in the z-direction at the point in time when the deviation from rest is largest. This time point (defined as the time after the Na + channel is opened) is specified in the upper right corner of each plot. In the y-direction, we focus on the 50 nm closest to the Na + channel clusters. The coordinates on the axes are shifted so that x = 0 marks the end of the membrane of the pre-junctional cell and y = 0 marks the center of the Na + channels. Note that the scaling of the colorbar is different for the different cases. seem to be somewhat smaller for this case (compare, e.g., Figure 12 in the main paper and Figure M).

S6 A different set of boundary conditions
In order to investigate how the choice of boundary conditions affect the results, Figures N-V show the results of simulations using a different set of boundary conditions than those used in the remaining simulations of this study. These boundary conditions are described in Section 2.3.3 in the main paper.
More specifically, we apply a Dirichlet boundary condition of on the rightmost boundary of the post-junctional cell. Here, E K is the Nernst equilibrium potential for K + ions, i.e., where k B , T , z K + and e are parameters specified in the paper, and c K + ,i is the intracellular K + concentration. Moreover, c K + ,e is the extracellular K + concentration measured in the extracellular point (L i +L m +L e /2, 0, 0) at the end of the simulation used to find the resting state. On the remaining part of the outer boundary of the computational domain, we apply the Neumann boundary condition ε∇φ · n = 0.
In addition, we apply no-flux Neumann boundary conditions for the ionic concentrations on the entire boundary of the computational domain D k ∇c k · n = 0, for k = {Na + , K + , Ca 2+ and Cl − }.
Comparing the results of Figures N-V to the corresponding figures for the default boundary conditions described in the paper, we observe that the results seem to be very similar for the two sets of boundary conditions. The perhaps largest difference between the two cases is that the extracellular K + concentration at rest for a narrow cleft of L e = 5 nm is about 1 mM higher in the simulations with Neumann boundary conditions everywhere for the concentrations than for the default boundary conditions (compare Figure 6 in the main paper and Figure N).
10 Figure E: Stationary solution for the potential, φ, the concentration of Na + , K + , Ca 2+ , and Cl − ions, and the charge density, ρ, in the extracellular space between the two cells in simulations with open K + channels, but closed Na + channels. The figure corresponds to Figure 6 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space. Figure F: Stationary solution of the potential, φ, the concentration of Na + , K + , Ca 2+ , and Cl − ions, and the charge density, ρ, along lines in the xdirection for open K + channels and closed Na + channels. The full green line represents the solution along a line crossing through the K + channels and the dotted black line represents the solution along a line about 100 nm below the K + channel cluster. The light green areas mark the membrane. The figure corresponds to Figure 7 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space.  Figure 8 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space. Figure H: The potential, φ, in the extracellular space between the two cells in simulations with open Na + channel clusters on the membrane of the prejunctional cell at the point in time when the most negative potential occurs. The figure corresponds to Figure 9 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space.  Figure 10 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space.  Figure 11 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space.  Figure C, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space.  Figure D, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space.  Figure 12 in the paper, but with a different profile for the initial conditions and ρ 0 in the ion channel clusters. Specifically, the initial condition in the channel clusters is a equal to the intracellular concentration in the half of the channel cluster that is closest to the intracellular space and equal to the extracellular concentration in the half that is closest to the extracellular space. Figure N: Stationary solution of the potential, φ, the concentration of Na + , K + , Ca 2+ , and Cl − ions, and the charge density, ρ, in the extracellular space between the two cells in simulations with open K + channels, but closed Na + channels. The figure corresponds to Figure 6 in the paper, but with a different set of boundary conditions applied. Specifically, no-flux Neumann boundary conditions are applied for the ionic concentrations on all outer boundaries of the computational domain.

S7 Different choice of intracellular concentrations
In Figures W-AD, we compare the results for two sets of simulations with different initial conditions for the intracellular concentrations. In the first set of simulations, we use the default initial conditions reported in the main paper, and in the second set of initial conditions, the intracellular Na + concentration is set to 8 mM and the intracellular Cl − concentration is set to 133.0002 mM. Otherwise, the set-up used in the two sets of simulations is the same. In Figures W-AD, we observe that the change in initial conditions appears to have no significant effect on the results.
29 Figure W: Stationary solution of the potential, φ, the concentration of Na + , K + , Ca 2+ , and Cl − ions, and the charge density, ρ, in the extracellular space between the two cells in simulations with open K + channels, but closed Na + channels. The figure corresponds to Figure 6 in the paper, but only considers L e = 10 nm, and compares the default initial conditions ([Na + ] 0 i = 12 mM) to a case with [Na + ] 0 i = 8 mM.      Figure AD: The membrane potential, v, the extracellular potential, φ, the extracellular Na + and K + concentrations, the charge density, ρ, and the total I Na current of the pre-and post-junctional cells as functions of time.
The figure corresponds to Figure 13 in the main paper, but only considers L y = L z = 300 nm, and compares the default initial conditions ([Na + ] 0 i = 12 mM) to a case with [Na + ] 0 i = 8 mM.