Testing the Applicability of Nernst-Planck Theory in Ion Channels: Comparisons with Brownian Dynamics Simulations

The macroscopic Nernst-Planck (NP) theory has often been used for predicting ion channel currents in recent years, but the validity of this theory at the microscopic scale has not been tested. In this study we systematically tested the ability of the NP theory to accurately predict channel currents by combining and comparing the results with those of Brownian dynamics (BD) simulations. To thoroughly test the theory in a range of situations, calculations were made in a series of simplified cylindrical channels with radii ranging from 3 to 15 Å, in a more complex ‘catenary’ channel, and in a realistic model of the mechanosensitive channel MscS. The extensive tests indicate that the NP equation is applicable in narrow ion channels provided that accurate concentrations and potentials can be input as the currents obtained from the combination of BD and NP match well with those obtained directly from BD simulations, although some discrepancies are seen when the ion concentrations are not radially uniform. This finding opens a door to utilising the results of microscopic simulations in continuum theory, something that is likely to be useful in the investigation of a range of biophysical and nano-scale applications and should stimulate further studies in this direction.


S1.1 The calculation of ion concentration from BD trajectory
The calculation of ion concentration is straightforward. As shown in figure S1, we consider 1D situation.
We take a segment of the channel z(1) ≤ z ≤ z(n) and divided it into n − 1 small segments each of which has a length of δ, so we have n sampling points (1, 2, . . . , n). Assuming i is one arbitrary sampling point, then the concentration of the ion species v at point i, C v (i), is calculated from the BD trajectories as follows: N v (i) is the number of ions v that locating at z(i) − δ 2 ≤ z ≤ z(i) + δ 2 in the entire production trajectory, R(i) is the radius of the channel at the sampling point i, and n f is the total amount of the frames in the production trajectories. By this calculation, the concentration of each ion species at all the n sampling points are known.

S1.2 The calculation of potential
In this study we calculated the potential at each sampling point using two steps. In the first we solved Poisson's equation to determine the total drop in potential, V , between the two ends of the segment used for the calculation, i.e., between sampling points 1 and n. For this calculation, the external applied electric field, the dielectric boundary, and the fixed charges on the channel were all considered. While the mobile ions were not. Then, we used the stationary NP equation to determine the potential at the intermediate points by utilising the values of the concentration at each point.
Assuming the potential at the first sampling point then the potential at the last sampling point To fulfil the stationary NP equation, we have Substituting the 1D NP equation into the equation S5, and noting that n v (i) = 10 3 N A C v (i), one gets If we consider the channel cross section area where C v (i) is the concentration of the ion species v, S(i) is the channel cross section area, and Φ v (i) is the potential at the ith sampling point respectively.
Equation S8 gives us n − 2 equations (i = 1, 2, . . . , n − 2). Combining them with the equations S2 and S3, we get a system of n equations. As the ion concentration C v (i), the channel cross section area S(i), and the potential difference V are already known, the system of equations can be solved to get the potential for each ion species Φ v (i). By using this strategy, it is guaranteed that the current through each sampling section is the same and the choice of segment on which the NP calculation is performed does not affect the result.

S1.3 The calculation of ion current
Once we have obtained the ion concentration and potential at each sampling point, the calculation of ion current can be done by: where n v (i) = 10 3 N A C v (i).

S1.4 The selection of the grid spacing
In our calculations, the grid spacing δ was set to 0.5Å, as our tests show that a value of 0.5 ∼ 1.0 gives the most stable prediction about the ion current. Using smaller values requires longer simulations to achieve the same level of sampling, while coarser grids can begin to compromise the results. The table S1 shows an example of how the grid spacing influences the calculated current found in the case of the 'real catenary channel'.

S1.5 The influence of the choice of the segment to perform the NP calculation
In one of our catenary channel calculations (cylindrical part has a radius of 8Å and spans from z=-7.5 A to z=7.5Å. The widest part has a maximum radius of 12Å at z=22.5Å) , we studied if selecting different segments to perform the NP equation would affect the accuracy. The results are shown in table S2. Table S2. The influence of segment selection on NP current (A). The BD results are 1.11×10 −10 A and -5.69×10 −11 A for Na + and Cl − respectively. segment -10∼10 -15∼15 -20∼20 -30∼30 Na + 1.12 ×10 −10 1.17 ×10 −10 1.19 ×10 −10 1.20 ×10 −10 Cl − -7.88 ×10 −11 -7.66 ×10 −11 -7.69 ×10 −11 -9.11 ×10 −11 As can be seen, only after including the ends of the channel for the calculation, there is a relatively larger deviation, especially for Cl − . Therefore, we suggest when doing such kind of calculation, one should try to take part of the channel and avoid including the ends of the channel or water reservoir.

S2 The influence of different dielectric constants
In the paper, all the calculations were performed with the dielectric constant of water set to 60.0. To help see the influence of this choice of dielectric constant on our results, we also performed similar simulations for passive, real and charged cylindrical channels with the dielectric constant of water to be 80 and all the other parameters be the same. The results are shown in figure S2.
Increasing the dielectric constant of water increases the ion conductance through the channel. However, as can be see in the figure  It should be noted that the average value found from the fluctuating currents as shown in figure S3 still do make sense and match well with the BD results for the cylindrical channels. For instance, the average value calculated for the above example using the middle 90 sampling points is 2.42 ×10 −12 A, which is very close to the BD result 2.06×10 −12 A. (The result using the BD-NP approach employed in the manuscript is 2.06 ×10 −12 A.) But for more complex channels, this method fails and can not give stable current values which are dependent on the choice of sampling points. The fact that the results become dependent on the choice of the sampling points in this direct method creates a significant drawback.

S3.1 Potential problem in our strategy
There is one problem about our strategy of calculating the potential from ion concentration: the potential profiles calculated for different ion species do not match with each other.