Our simplified Brownian dynamic gap junction model (Mondal et al., 2017) was developed using MATLAB (Mathworks Inc., Natick, MA). In short, particles were modeled as spheres with point charges at their center. Stoke's radius values were adjusted for different particles size. The system was charge balanced by adding appropriate counter-ions. The simulation parameters were set to mimick in vitro dye diffusion under double whole cell voltage clamp experiments with coupled cells (Figures 1A–D). We modeled all our pore shape studies with Lucifer yellow (yellow arrows in Figures 1A–D) particles and Cesium (Cs) counter-ions with a 2 mV transjunctional voltage. A 2 mV transjunctional voltage would give a flux equivalent to one from 10 mV/200 ms/1 Hz pulses. A constant concentration gradient of 2 LY and 4 Cs (8 mM LY) was maintained across the pore by replacing any particle that crossed over to cell 2 (central plane of cell 2) back to cell 1. A time step of 1 ps was set. The program saved particles' positions every 10^th time step (10 ps). Data files containing particle positions were saved every 10 μs. The LY particle flux was calculated every 10 μs considering the total number of particles and simulation time. The simulation was terminated when the last 5 calculated LY flux values had a standard deviation <5% or simulation total time was 500 μs.

The simulations were run on 4 Apple Mac Pro computers involving a total of 44 (12 × 3 + 8) processing cores. Our simulation studies were designed to assess the role of pore shape, particle size, and particle charge on flux asymmetry in heterotypic pores. Simulation conditions were defined accordingly as described below.

The averaged gap junction pore was divided into 3 sections defined by 4 radii (R₁, R₂, R₃, and R₄) at the section ends placed in series. A pore was identified by listing these pore-section radii (in Å) in order from left to right as R₁-R₂-R₃-R₄. Homotypic pores were symmetric about the vestibule, R₁ = R₄ and R₂ = R₃ (Figure 1E) and were denoted as R₁-R₂-R₂-R₁. Heterotypic pores studied in this work were always in a cone-cylinder combination of hemichannels, which meant R₃ = R₄. This was a required condition for activating the mechanism that resulted in asymmetric fluxes in heterotypic pores (Mondal et al., 2017). So, in all our simulations, heterotypic pores will be of the form R₁-R₂-R₃-R₃. A range of values were set for the outer (R₁) and inner pore (R₂) radii taking into consideration the particle size and computation time. In our studies, inner mouth radii ranged from 6.2 to 12.2 Å and outer mouth radii ranged from 6.2 to 25.9 Å. In all our pore radii combination, outer mouth radii were never smaller than the inner mouth radii to avoid inverted conical pore sections.

Simulation data were collected and processed for analysis on an 8 core Apple Mac Pro computer using MATLAB R2016a. This included data verification, tabulating particle fluxes, calculating particle probabilities and sectional forces and data organization. The final result tables and plots were created with Excel 2013 (Microsoft Inc., Redmond, WA).

Fluxes across homotypic pores do not present asymmetry. Table S1 list the fluxes for all 56 homotypic pores considered. Figures S1–S4 present LY fluxes, particle probability profiles and average electrical force data arranged to study variations in outer (R₁) and inner (R₂) mouth radii. LY flux values in our simulation with 8 to 0 mM LY gradient ranged from 3,551 molecules per pore per s (m/c/s) to 64,588 m/c/s in homotypic pores (Table S1) and 2,919 m/c/s to 51,409 m/c/s in heterotypic pores (Table S2). Physiological reported values of LY flux range from 4,000 m/c/s (m/c/s) (Valiunas et al., 2002; Ek-Vitorín and Burt, 2005; Ek-Vitorín et al., 2006; Heyman and Burt, 2008; Kanaporis et al., 2013) to 300,000 m/c/s (Weber et al., 2004).

We conducted these studies using the heterotypic pore profile of 25.9-12.2-8.2-8.2. As expected, fluxes dropped as particle size was increased (Figure 6A). Fluxes were higher in the cylinder to cone direction (orange bars in Figure 6A). Flux ratio was maximal for r_p = 4.9 Å, which is the Stokes radius of LY. Even though the flux ratio dropped on either side of r_p = 4.9 Å, the flux asymmetry was still pronounced (~0.7).

The general P_dX profile remained the same for all particle sizes in their respective direction of flux simulations (Figures 6B,C). The net drop in P_dX across the pore was directly proportional to the particle size in both directions. In cell 1, P_dX remained in the same order of magnitude. The maximum difference in P_dX among different r_p values was seen at the narrowest section of the pore. This includes the mouths on either side of the vestibule and the cylindrical mouth.

In general, F_X magnitudes dropped nonlinearly with increasing particle size (Figures 6D–H). F_X magnitudes in cell 1 in the cone to cylinder direction became more distinct at smaller particle sizes (r_p = 0.8 to 2.7 Å), following a trend of F_X magnitudes dropping with increasing particle size (Figure 6D). F_X magnitudes in the cone to cylinder direction reduced with increasing particle size and saturated around r_p = 4.9 to 6 Å in the left mouth (Figure 6E) and vestibule (Figure 6F). In the left and right mouths, r_p = 6 Å showed small rise in F_X (Figures 6E,G). In the left mouth, F_X magnitude was always higher in the cylinder to cone direction (Figure 6E). On the other hand, F_X magnitudes were always higher in the cone to cylinder direction in the vestibule (Figure 6F) and right mouth (Figure 6G).

In this set of simulation, LY particle radius was fixed to r_p = 4.9 Å. Particle flux reduced with increased particle charge in either direction (Figure 7A). Flux ratio, on the other hand, increases and seems to reach maximal at q_p = 2e⁻ (gray bars in Figure 7A).

P_dX profiles (Figures 7B,C) were similar to the general ones seen before in the previous cases. The higher the particle charge, the lower was the P_dX level.

Since at zero charge there were no electrical forces as interparticle electrostatic forces and applied transjunctional voltage were absent, F_X for 0e⁻ have not been presented in sectional F_X figures (Figures 7D–H). F_X differences between the two directions rose with increasing negative charge. F_X in cell 1 was marginal in magnitude compared to F_X in other sections. In the left mouth, F_X was negative and distinctly higher in the cylinder to cone direction than in the cone to cylinder direction for q_p > 1e⁻ (Figure 7E). In the vestibule (Figure 7F) and right mouth (Figure 7G), it was the other way round with F_X larger in the cone to cylinder direction.

Our simulations with heterotypic pores predicted the maximal possible flux asymmetry due to pore shape. We found that flux asymmetry was directly proportional to the width of the conical section from our R₁-9.2-9.2-9.2 results in Figure 2A. Also we found that a constant cylindrical mouth section of R₃ = 8.2 Å yielded the highest flux asymmetry, as shown in our flux simulation on heterotypic pores for 2 different cone mouth sizes of 22.5-9.2-R₃-R₃ (Figure 4A) and 25.9-12.2-R₃-R₃ (Figure 5A). Our previously published in vitro data with Cx43-Cx45 reported α_het = 0.34 and corresponding α_hom = 3.88 with Cx43 as the left hemichannel (red in Figures 1C,D) and Cx45 as the right hemichannel (white in Figures 1C,D) in the heterotypic junctions. The in silico flux data for heterotypic pore profile of 22.5-9.2-8.2-8.2 (S# 18 in Table S2) was identified to best fit these data points with α_het = 0.55 and α_hom = 3.38. This certainly makes a strong case where in vitro data could be explained by pore shape variation to a very significant degree. Perhaps a finer optimization might yield exactly matching data. Another close candidate from our simulations was 25.9-12.2-9.2-9.2 (S# 26 in Table S2) with α_het = 0.54 and α_hom = 4.98.

Heterotypic pore profile of 25.9-12.2-8.2-8.2 gave maximal flux asymmetry with α_het = 0.49, but the corresponding α_hom = 9.18. This was one of the highest magnitudes of preferential flux reversals (i.e., α_het < 1 and α_hom > 1 or vice versa) between the heterotypic pore and the homotypic pores of the corresponding hemichannels. Considering α_het < 0.75 to be a threshold for significant flux asymmetry, from the 28 heterotypic pores, 14 exhibited flux asymmetry (Table S2). Out of those 14, 13 had preferential flux reversal between the heterotypic and their homotypic counterparts (comparing α_het and α_hom columns in Table S2).

We conducted separate simulations for studying the effect of pore sectional radius in heterotypic pores, particle size and charge on flux asymmetry. We first identified the radii combination that yielded the maximum flux asymmetry (25.9-12.2-8.2-8.2) and performed the particle size and charge studies using this heterotypic pore to maximize sensitivity. In our previous work, we established that particle charge was a vital element toward creating asymmetric fluxes and this asymmetry disappeared when charge was removed. Our results in this work indicate that while charge is vital, it acts as a scaling factor to flux asymmetry, where it increases with rising charge magnitude to about q_p = 2e⁻ (Figure 7A). Particle size seems to play a similar role, though there was no clear relationship between particle size and flux asymmetry. In fact, flux asymmetry was present for all particle sizes considered with maximal asymmetry at r_p = 4.9 Å (Figure 6A). We note that even though charge was a vital factor, the flux asymmetry was purely regulated by the pore shape. Variation of different aspects of the pore-shape in heterotypic pores had strong effects on flux asymmetry. The pore shape features governed the direction and magnitude of asymmetry. Flux asymmetry was directly proportional to the outer mouth size of conical section R₁ (Figure 2A), nonlinearly dependent on the cylinder mouth size R₃ (Figures 4A, 5A) and varied in a multimodal pattern with inner pore mouth size R₂ (Figure 3A). As noted earlier, only half the heterotypic pores exhibited flux asymmetry with charged particles (Table S2). From our results, we conclude that flux asymmetry was primarily regulated by the heterotypic pore's conical and cylinder mouth sizes. Particle charge and size were secondary factors.

In our previous work we suggested a mechanism in which the pore shape controlled the distance between particles in a pore and the charge cloud at the end of the pore (cell 1 in our simulations). A wider or conical mouth pore section allowed this charge cloud to invade into the conical section of the pore, reducing the distance between the charged cloud and any particle that got into the pore, thus increasing the pull-out (in -x direction) electrical force. When the charge cloud was placed on the narrower mouth side, it maintained a higher distance and thus a lower pull-out electrical force. While the underlying mechanism remains the same, particle probability distribution also plays a key role and flux regulation cannot be solely explained based on sectional force differences. Particle probability profiles across the pore were analogous to concentration gradient. P_dX profiles in all cases tend to resemble concentration gradients we would expect in simple diffusion through those pore shapes, but in exponential scale. The narrowness of the pore amplified the drop in probability (or concentration) along x. The smaller the pore section size, the higher was the drop in P_dX across the section and vice versa. This trend in P_dX was particularly consequential in the narrowest sectional pore radii (R₁, R₂, and R₃ = 6.2 Å) in homotypic (Figure S1B, R₂ = 6.2 Å in Figure S3B and R₁ = R₂ = 6.2 Å in Figure S4B) and heterotypic (R₂ = 6.2 Å in Figures 3B,C, R₃ = 6.2 Å in Figures 4B,C, 5B,C) pores. Such drops caused major variations in net charge distribution, leading to the electrical force distribution across the pore. An extreme example of this can be seen in right mouth F_X data for R₃ = 6.2 Å in 22.5-9.2-R₃-R₃ (Figure 4F) and 25.9-12.2-R₃-R₃ (Figure 5F) where F_X were positive with very high magnitudes in the cone to cylinder direction (negative for all other cases in the same direction). This implied large flux as LY particles were pushed out of the pore in the direction of flux. But, because the right mouth has the lowest P_dX level (R₃ = 6.2 Å in the right mouth section in Figures 4B, 5B), there were few particles to push out. So if any particles ended up in the right mouth, they moved out from the pore to cell 2 quicker than in other pore shapes. This eluded a simple relationship to explain or predict accurately why maximum flux asymmetry occurs at 25.9-12.2-8.2-8.2. Even heterotypic pores with no flux asymmetry showed pronounced F_X differences amongst each other (R₂ = 6.2 Å in Figures 3A,D–F and R₃ = 10.2 to 11.2 Å in Figures 4A,D–F). As of now, the best way to select sections for F_X values consequential to flux is to consider them in order of their distance from cell 1 where the particle probability was the highest. Following this criterion, the left mouth and vestibule F_X were most consequential. The right mouth might or might not be significant based on relative P_dX levels in the right mouth and sections to the right of it. We presented only the F_X values for the left mouth, vestibule and right mouth for our heterotypic pore studies here. F_X in cell 1 and cell 2 (along with others) have been presented in the supplemental section (Figure S5). F_X magnitudes in cell 1 were marginal compared to other sections. Though F_X magnitudes were pronounced in cell 2, and they had minimal effect on final flux value as they had the minimal P_dX levels in the pore.

Molecular permeability of biological molecules in gap junction has been extensively studied for key functional cytoplasmic molecules, including secondary messengers like IP₃, cAMP, cGMP, and Ca²⁺ (Bevans et al., 1998; Niessen et al., 2000; Bedner et al., 2003, 2006; Locke et al., 2004; Ayad et al., 2006), nucleotides like RNA and siRNA (Valiunas et al., 2005), and glucose and its metabolites like ATP, ADP and glutamate (Goldberg et al., 1999, 2002). These studies reported 3 to 33-fold differences in permeabilities between different molecules in different connexin types. Our in vitro studies focused on homotypic and heterotypic pores formed out of Cx43 and Cx45: connexin pores present in the heart. While Cx43 predominates in the ventricular myocytes, in diseased or failing human ventricles Cx43 expression decreases while that of Cx45 increases (Yamada et al., 2003), potentially leading to a rise in heterotypic Cx43-Cx45 junctions. Such heterotypic junctions not only increase impedance to conduction in the ventricles, but can also cause functional effect by creating concentration imbalances of secondary messengers like Ca²⁺ and IP₃ through asymmetric fluxes amongst neighboring myocytes. Such imbalances can create disruptions in excitation contraction coupling. Moreover, such Cx combinations can cause variations in glucose availability among neighboring cells, which can hinder ATP production. Furthermore, ventricular fibrosis (Krenning et al., 2010) may cause coupling of myocytes to Cx45 expressing non-myocytes: fibroblast (Gaudesius et al., 2003) and myofibroblasts (Rohr, 2009). Fibroblasts being 20 to 40 times smaller in volume than myocytes would be more easily able to create changes in the transjunctional (Cx43-Cx45) concentration gradient, leading to higher concentration of certain molecules in the myocyte (assuming Cx43 is conical and Cx45 is cylindrical as per our simulation results). This remains to be demonstrated in further studies.

Our model was primarily optimized for studying LY fluxes in gap junction pores mimicking in vitro LY dye diffusion under double whole cell voltage clamp experiments. Most of our simulation data reported here were for LY particles with Cs counter-ions. Our simulations on particle size and charge variation deviated from this and faced limits due to computational and statistical constraints. We allowed a maximal simulation time of 500 μs or the time when the fluxes in the last 5 values had less than 5% standard deviation. These were set taking into consideration the computational resources available, computation time per simulation and the number of different simulation conditions in our studies. This provided limitations to the smallest measurable flux in our model and restricted the minimum pore mouth size to 6.2 Å, maximum particle size to 6.6 Å and maximum particle charge to 2.5e⁻. Going beyond these limits either gave us a considerable number of simulations not reaching the 5% standard deviation criteria or just 0 flux values. Therefore, for parameters leading to lower fluxes, the standard deviation limit will need to be reduced to a lower threshold or the simulation time gap between the points considered for standard deviation would have to be increased.

While the presented results were focused around a neutral pore, pore electrostatics due to wall surface charges can have major significance in regulating gap junction permeability of charged molecules. Our in vitro experimental results with Cx43-Cx45 heterotypic channels exhibited particle charge independence (Mondal et al., 2017), but Cx43-Cx40 heterotypic channels have been shown to be charge selective (Lin et al., 2014). A lack of structural details of Cx isoforms is a major hurdle for accurate modeling at this stage. Simplified models implementing pore wall charges in in silico studies with Cx26 channels have been reported (Escalona et al., 2016). Such simplified application of pore electrostatics can be envisioned to study the role of different combinations of charges and in-pore charge-locations.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.