Modulation of Asymmetric Flux in Heterotypic Gap Junctions by Pore Shape, Particle Size and Charge

Gap junction channels play a vital role in intercellular communication by connecting cytoplasm of adjoined cells through arrays of channel-pores formed at the common membrane junction. Their structure and properties vary depending on the connexin isoform(s) involved in forming the full gap junction channel. Lack of information on the molecular structure of gap junction channels has limited the development of computational tools for single channel studies. Currently, we rely on cumbersome experimental techniques that have limited capabilities. We have earlier reported a simplified Brownian dynamics gap junction pore model and demonstrated that variations in pore shape at the single channel level can explain some of the differences in permeability of heterotypic channels observed in in vitro experiments. Based on this computational model, we designed simulations to study the influence of pore shape, particle size and charge in homotypic and heterotypic pores. We simulated dye diffusion under whole cell voltage clamping. Our simulation studies with pore shape variations revealed a pore shape with maximal flux asymmetry in a heterotypic pore. We identified pore shape profiles that match the in silico flux asymmetry results to the in vitro results of homotypic and heterotypic gap junction formed out of Cx43 and Cx45. Our simulation results indicate that the channel's pore-shape established flux asymmetry and that flux asymmetry is primarily regulated by the sizes of the conical and/or cylindrical mouths at each end of the pore. Within the set range of particle size and charge, flux asymmetry was found to be independent of particle size and directly proportional to charge magnitude. While particle charge was vital to creating flux asymmetry, charge magnitude only scaled the observed flux asymmetry. Our studies identified the key factors that help predict asymmetry. Finally, we suggest the role of such flux asymmetry in creating concentration imbalances of messenger molecules in cardiomyocytes. We also assess the potency of fibroblasts in aggravating such imbalances through Cx43-Cx45 heterotypic channels in fibrotic heart tissue.

Introduction

Gap junction channels provide pores through cell membranes between two adjoined cells, allowing transport of cytoplasmic molecules amongst them. Each gap junction channel is formed of two coaxially placed connexons (or hemichannels) with one in each cell membrane. Each of these connexons are made of 6 circumferentially arranged connexins (Cx). To date, 21 different Cx types have been identified in humans (Söhl and Willecke, 2004). Gap junction channels are called homotypic when Cx in both hemichannels are made of the same type, and heterotypic if the Cx in each hemichannel is different. They can also be heteromeric or hetero-multimeric where multiple Cx types constitute a hemichannel. Gap junctions play key roles in intercellular communication in different tissues like the skin, nervous system, respiratory epithelium, bones, lens and retina of the eye, inner ear, heart, vasculature, and reproductive systems. Density and regulation of gap junctions are essential for cellular and tissue physiology. How these channels are relevant in the development of several diseases has become clearer as more regulatory mechanisms that involve channel expression, localization and cellular modifications are known (Spray and Dermietzel, 1995; Kelsell et al., 1997; Xia et al., 1998; Jongsma and Wilders, 2000; Uhlenberg et al., 2004).

Currently, x-ray crystallographic structures of only Cx26 and Cx43 are available (Unger et al., 1999; Maeda et al., 2009; Bennett et al., 2016), with data on molecular structure and in-pore charge distribution available only for Cx26. Structural details are a requisite for using existing computational approaches for studying membrane channels like ion channels (Allen et al., 1999; Corry et al., 2001; Chung et al., 2002; Vora et al., 2008) and porins (Schirmer and Phale, 1999; Im and Roux, 2002; Lee et al., 2011). Molecular modeling of Cx26 hemichannels (Kwon et al., 2011; Zonta et al., 2014) and gap junctions (Bennett et al., 2016; Escalona et al., 2016) was performed for studying their structure/function relations. Limitations and difficulties in deriving such structures pose a barrier to studying the structure and permeability of the whole range of gap junction channel types. Thus, most studies rely heavily on in vitro experimental methods like dye diffusion, double whole cell voltage and patch clamping. Fluorescent dyes were a vital tool in such studies with Lucifer yellow (LY) being the most commonly used for dye diffusion studies (Verselis et al., 1986; Cao et al., 1998; Martinez et al., 2002; Valiunas et al., 2002; Beltramello et al., 2003; Dong et al., 2006; Eckert, 2006; Rackauskas et al., 2007; Yum et al., 2007). Gap junctions serve as conduits for intercellular diffusion and transport studies of physiologically relevant molecules and ions were also conducted (Verselis et al., 1986; Veenstra et al., 1994; Goldberg et al., 1999, 2002; Niessen et al., 2000; Beyer et al., 2001; Veenstra, 2001; Martinez et al., 2002; Bedner et al., 2003, 2006; Locke et al., 2004; Valiunas et al., 2005; Ayad et al., 2006; Bukauskas et al., 2006; Harris, 2007). Our research focused on homotypic and heterotypic gap junction channels involving Cx43 and Cx45, which are abundantly present in the heart (van Kempen et al., 1995; Vozzi et al., 1999; Severs et al., 2001). In previous works, we found that Cx43 and Cx45 heterotypic channels do not exhibit electric current rectification beyond rectification caused by voltage dependent gating, which is significantly different for each connexon (Elenes et al., 2001). Calculations of dye diffusion using cell pairs under double whole cell voltage clamp experiments indicate that homotypic Cx43 channels allow higher LY flux than Cx45 channels, while heterotypic Cx43-Cx45 channels show asymmetric fluxes (Martinez et al., 2002; Moreno et al., 2004; Mondal et al., 2017). Interestingly, the recorded fluxes were higher in Cx45 to Cx43 direction than Cx43 to Cx45. Based solely on their respective unitary conductance, it was theorized that Cx43 would be a wider hemichannel than Cx45, which explained the difference in homotypic fluxes, but failed to explain the flux difference in the heterotypic Cx43-Cx45 junctions. We speculated on a structural basis for the mechanism creating this asymmetry, which was also proposed before to explain similar findings (Robinson et al., 1993). The proposed hypothesis was criticized at the time as it implied the presence of an active mechanism (Finkelstein, 1994) and found to be inconsistent with unitary conductance and permeability values from in vitro experiments (Veenstra et al., 1995; Veenstra, 1996). Subsequently, using a simplified Brownian dynamics computational model, we showed that a heterotypic pore can generate such an active mechanism by creating electrical fields across the pore through motion and redistribution of charged particles caused by the pore's shape (Mondal et al., 2017). The combination of a conical and cylindrical hemichannel was found to be a requisite for the underlying mechanism.

Here we use the simplified Brownian dynamics model of a gap junction pore presented in our previous work to mimic LY diffusion experiments under whole cell voltage clamp. We designed our studies to elucidate the roles of pore shape, particle size and charge on flux and flux asymmetry of large molecules. Simulation studies were conducted for both homotypic and heterotypic pores. Our studies revealed a heterotypic pore profile with maximum flux asymmetry and also pore-shapes that matched our previously reported in vitro results with Cx43 and Cx45 pores. We also found that within the set limits, flux asymmetry was almost independent of particle size and scaled by particle charge. While it was not possible to establish a direct relationship between pore-shape, particle charge and size, we identified factors for predicting flux asymmetry. We also speculated on a possible role of such asymmetric fluxes in heart disease.

Methods

Simulation Model

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.

FIGURE 1

Figure 1. In vitro and in silico experiment designs with homotypic and heterotypic gap junction pores. (A–D) Basic setup for in vitro dye diffusion experiments under double whole cell voltage clamp with homotypic (A) Cx43 and (B) Cx45 gap junction pores and (C,D) heterotypic Cx43-Cx45 gap junction pores with LY (yellow arrow) injected from (C) Cx43 side and (D) Cx45 side. Voltage pulses of 10 mV/200 ms/1 s were applied and LY dye diffusion was recorded across the gap junction. (E,F) Simulation experiment design for single (E) homotypic and (F) heterotypic gap junction pore. The cell-pore-cell system was divided in 5 sections: (1) cell 1, (2) left mouth, (3) vestibule (4), right mouth, and (5) cell 2.

Simulation Experiment Design

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.

Flux Simulations in Homotypic Pores (R₁-R₂-R₂-R₁)

Table 1 lists the inner (R₁) and outer (R₂) mouth radii considered for forming homotypic pores. Figure 1E shows the homotypic pore simulation setup. Based on our set criteria for selecting combinations, 56 homotypic pores were considered. Flux simulations were run for each pore 30 times (n = 30).

TABLE 1

Table 1. Sectional radii of homotypic pores.

Flux Simulations in Heterotypic Pores

Our previous work (Mondal et al., 2017) was focused around the hypothesis that flux asymmetry in heterotypic pores can be produced if one of the pore mouths was cylindrical and the other was conical. Simulation were designed to study the effects of R₁, R₂, and R₃ in different simulation sets (rows 1–3 in Table 2). Each simulation case was repeated 50 times (n = 50) on either direction (cone to cylinder and cylinder to cone) of the heterotypic pore. In the 1^st set of simulations (row 1 in Table 2), R₁ was varied while R₂ and R₃ were fixed to 9.2 Å. This basically simulates the effect of varying cone angle of the conical mouth of the pore. In the 2^nd set of simulations, R₂ was varied, with R₁ = 22.5 Å and R₃ = 9.2 Å. Here, the influence of inner mouth size was studied. In the 3^rd set, the cylinder mouth size was varied while the left conical mouth was kept constant, with R₁ = 22.5 Å and R₂ = 9.2 Å. Finally, in our 4^th simulation set, the conical mouth was set to the maximal R₁ and R₂, with R₁ = 25.9 Å and R₂ = 12.2 Å and the effects of varying the cylinder mouth were investigated.

TABLE 2

Table 2. Sectional radii combinations for heterotypic pores.

Flux Simulations with Varying Particle Size

In this set of simulations, the particle size was varied from 0.8 to 6.6 Å. Row 1 in Table 3 lists the heterotypic pore and the particle radii values used. While the particle size was varied, the counter-ion size was kept constant. Particle charge was set to 2e⁻.

TABLE 3

Table 3. Particle sizes and charges for heterotypic pore simulations.

Flux Simulations with Varying Particle Charge

In our previous work, particle charge was an important factor in creating the observed flux asymmetry of large ionic particles. In this set of simulations, we varied the particle charge (Row 2 in Table 3) to study the influence of charge on flux asymmetry. The number of counter-ions was accordingly changed to maintain charge balance. Particle radii of LY and Cs were not changed.

Simulation Output and Analysis

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).

Particle Flux and Flux Asymmetry

Particle flux J was calculated by counting the total number of particles that crossed the pore n_p divided by the total simulation time t.

$\begin{array}{ll}J=\frac{n_p}{t}& (1)\end{array}$

In homotypic pores, flux was only measured in the x direction where the concentration gradient of 2LY + 4Cs to 0 particles was maintained between cell 1 and cell 2 and denoted by J_hom (Figure 1E). In heterotypic pores, flux was measured in the x direction with 2LY + 4Cs to 0 particles in cell 1 to cell 2 (or cone to cylinder direction), and the -x direction with 2LY + 4Cs to 0 particles in cell 2 to cell 1 (or cylinder to cone direction) and denoted by J_het,₁ and J_het,₂, respectively (Figure 1F).

Flux ratio was the measure of flux asymmetry α_het in a heterotypic pore and defined as the ratio of the LY (or particle) fluxes in cone to cylinder (J_het,₁) and cylinder to cone (J_het,₂) directions.

$\begin{array}{ll}{\alpha }_{het}=\frac{J_{het,1}}{J_{het,2}}& (2)\end{array}$

Heterotypic flux ratios were plotted in log scale (base 2) to accurately scale the magnitude of flux asymmetry in reversed cases (e.g., α_het of 0.5 and 2 are equal in magnitude but in opposite directions).

Ratio of the homotypic fluxes of the corresponding hemichannels in a heterotypic pore α_hom was also calculated to identify reversal of flux asymmetry.

$\begin{array}{ll}{\alpha }_{hom}=\frac{J_{hom,1}}{J_{hom,2}}& (3)\end{array}$

where J_hom,₁ and J_hom,₂ are the fluxes in corresponding homotypic pores formed of the left and right hemichannels, respectively.

Particle Probability

Particle probability profiles along the pore's axial profile (along x-axis) were calculated from the simulation data and plotted. Axial particle probability P_dX was the probability of particle presence within a cylindrical section of width dX at x of the pore:

$\begin{array}{ll}{P}_{dx}=\frac{\Delta t}{t_{total}}{\displaystyle \sum _{i=1}^{n_p}{\displaystyle \sum _{n=1}^{n_t}{p}_{dx}}}(i,n)& (4)\\ where{p}_{dX}\left(i,n\right)=\{\begin{array}{c}1,x\left(i,n\right)lies\text{within}dX\\ 0,\text{otherwise}\end{array}\end{array}$

Here, Δt is the time step, t_total is the total simulation time, N_p is the number of particles, N_t is the number of time steps, x is the particle position, and i and n are the particle and time indexes, respectively.

Average Electrical Force per Particle per Pore Section

The average electrical force per particle per pore section was calculated as in our previous work (Mondal et al., 2017). In short, the pore was divided into 5 sections in series: (1) cell 1 (C1), (2) left mouth (LM), (3) vestibule (VT), (4) right mouth (RM), and (5) cell 2 (C2) as shown in Figures 1E,F). The average interparticle electrostatic and electric field forces per particle type F was calculated. The x-component of F F_X was plotted and compared amongst the different simulation conditions.

While the actual simulations were done in +x and −x directions, the −x data for P_dX and F_X data presented here were flipped horizontally so all flux direction are in the +x direction. This simplified it for the reader to compare the presented data in the different cases.

Results

The results presented in the next sections were derived from ~2.2 TB simulation data resulting from ~28.3 days of computation time per processing core (44 processors). This excludes data handling and postprocessing done on the central workstation (8 processing cores).

Flux Simulations in Homotypic Pores

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).

Outer Mouth Variation in Conical Mouthed Pores: R₁-6.2-6.2-R₁ and R₁-12.2-12.2-R₁

Results for fixed inner mouth radii (R₂) of 6.2 Å (narrow mouth) and 12.2 Å (wide mouth) with varying outer mouth radius are presented in Figures S1, S2, respectively. In pores with narrow inner mouth (R₂ = 6.2 Å), an increasing outer mouth marginally changed the LY flux (Figure S1A). When the inner mouth radius was widened (R₂ = 12.2 Å), increasing the outer mouth radius increased the LY flux which tended to saturate at higher values of R₁ (Figure S2A).

The P_dX profile in cell 1 was of the same order in all cases (Figures S1B, S3B) indicating a constant LY concentration in cell 1. As the outer mouth was varied, the P_dX profile varied across the pore. The vestibule size was kept constant and for wide inner pore mouth (R₂ = 12.2 Å), the P_dX profiles in the vestibule were of the same order. The P_dX profiles tended to be noisier with increasing distance along x from cell 1, and were noisiest at the rightmost region (cell 2). In general, the particle probability profiles were reminiscent of their expected concentration profiles based on the pore's radius profile in simple diffusion, but in the logarithmic scale (base 10). R₂ = 6.2 Å created probability profiles with the most variations along the pore sections, while R₂ = 12.2 Å presented a more gentle probability variation.

Homotypic pores with R₂ = 6.2 Å had the highest F_X magnitudes (Figures S1C, S3C). These were mostly observed in the left and right mouth of the pore. F_X also varied the most in those sections with varying outer mouth sizes. In the left mouth though, the trend was nonlinear. F_X rose in magnitude up to 8.2 Å and then dropped as outer mouth size was further increased to 25.9 Å.

The case R₁ = R₂ = 6.2 Å showed highest F_X magnitudes in the right mouth (Figure S1C). F_X decreased as the outer mouth size was increased. In pores with R₂ = 12.2 Å (Figure S2C), F_X showed less variation than that in R₂ = 6.2 Å. A narrow inner mouth caused more variations in F_X in the left and right mouths (Figure S1C) while a wide mouth created notable F_X variations in the vestibule (Figure S2C). Furthermore, in narrow mouthed pores, F_X in the left mouth rose up to R₁ = 15.8 Å and then dropped with increasing R₁. In the wide mouthed pore, F_X in the vestibule had almost a linear rise with increasing R₁.

Inner Mouth Variation in Conical Mouthed Pores (22.5-R₂-R₂-22.5)

Simulation data were selected and arranged for R₁ = 22.5 Å and varying inner mouth size of R₂ = 6.2 Å to 12.2 Å (Figure S3). LY flux exhibited a nonlinear rise with increasing inner mouth size (Figure S3A). P_dX remained almost the same up to the initial section of the left mouth (Figure S3B). Pronounced variation among the varying inner mouth size conditions was observed closer to the inner mouth and to the right of that region. Variations were amplified at the inner mouth regions on both sides of the vestibule. F_X in the left mouth and vestibule were negative and most distinct. While F_X magnitude dropped with increased R₂, it rose in the vestibule. Increasing R₂ did not generate noticeable variations in F_X in cell 2.

Cylindrical Mouthed Pores (R₁-R₂-R₂-R₁)

Figure S4 summarizes results from flux simulations with cylindrical pores (R₁ = R₂) where the cylinder size was increased from 6.2 to 12.2 Å. The complete pore with R₁ = R₂ value of 12.2 Å was almost a cylinder as the elliptical vestibule almost disappears (bottom right corner of Figure S4A). LY flux rose nonlinearly with increasing pore mouth size (Figure S4A). P_dX profiles of the different pore sizes ran parallel in the left and right mouth sections (Figure S4B). Vestibular sections of the P_dX profile straightened with increase in cylinder radius. As before, the probability profile followed the pore-section in the logarithmic scale. At 12.2 Å, the particle probability was almost a straight line with no elliptical bump at the vestibule region. The left mouth and right mouth had the most pronounced F_X magnitudes (Figure S4C). F_X was negative in all cases, except in the right mouth for cylinder radius of 6.2 Å. In the left mouth, F_X linearly decreased with increasing cylinder radius. F_X magnitudes in cell 1 were marginal. In cell 2, the differences in F_X for different cylinder radii were marginal.

Flux Simulations in Heterotypic Pores

Varying Outer Mouth Size of Conical Section (R₁-9.2-9.2-9.2)

Increasing the outer mouth size of the conical section marginally changed LY flux in the cone to cylinder direction (blue bars in Figure 2A). However, LY flux in the cylinder to cone direction rose (orange bars in Figure 2A). This caused an increase in the difference between the fluxes in the two direction as R₁ was increased. This contributed to an increase in flux ratio, which shows a gentle rise with increasing R₁ (gray bars in Figure 2A). The cone with the widest mouth showed the highest asymmetry (R₁ = 25.9 Å).

FIGURE 2

Figure 2. Properties of heterotypic pores with varying outer mouth size of conical section (R₁-9.2-9.2-9.2). (A) LY fluxes and flux ratios simulated with constant concentration gradient of 2LY + 4Cs in cone to cylinder (blue bars) and cylinder to cone (orange bars) directions. The difference between J_het,₁ and J_het,₂ increased with larger conical mouth size. LY flux asymmetry increased with increasing outer mouth size of the conical section. (B,C) P_dX profiles in simulation (B) cone to cylinder direction and (C) cylinder to cone direction. (D–F) Average force per LY particle in (D) left mouth, (E) vestibule, and (F) right mouth. F_X was in the -x direction in all sections of the pore. (D) F_X magnitudes in the left mouth in either direction were almost equal for different pore mouth sizes. (E) F_X magnitudes in the vestibule rose in the cone to cylinder direction with increase in R₁. (F) The right mouth had high differences in F_X magnitudes between the two directions with F_X magnitude dropping with increasing R₁ in the cylinder to cone direction.

As in homotypic pores, the P_dX profiles varied primarily due to the differences in pore shapes and the differences were most prominent in the section where size was varied. The differences were most prominent in the conical mouth in either direction (Figures 2B,C). In the cone to cylinder direction simulations, P_dX levels in the left mouth were seen to drop in the left mouth as R₁ was reduced (Figure 2B). In the remaining pore (vestibule and right mouth), the P_dX profiles were of the same shape and in the same order of magnitude. In the cylinder to cone direction simulations (Figure 2C), the P_dX levels began to notably vary from the end of the left mouth (end of conical mouth) itself. Particle probabilities were higher in the vestibule for smaller values of R₁ and declined as R₁ was increased. This trend was reversed in the center region of the right mouth (right mouth section in Figure 2C).

As in homotypic pores, F_X in cell 1 was marginal compared to F_X in other pore sections (Figures S5A–E). The left (Figure 2D) and right (Figure 2F) mouths had high F_X values. F_X values were negative throughout the pore in simulations of both direction. The left mouth F_X differences between either directions were marginal for all cases. F_X in the vestibule (Figure 2D), though small in magnitude, showed a widening difference between the 2 directional fluxes as R₁ was increased with F_X in the cone to cylinder direction being larger (Figure 2E). Differences in F_X magnitudes were high in the right mouth (Figure 2F).

Varying Inner Mouth Size of Conical Section (22.5-R₂-9.2-9.2)

Flux values rose with increasing inner mouth size in both directions (blue and orange bars in Figure 3A). Difference in flux in either direction seemed to fluctuate as R₂ was increased. This variation appeared to be multimodal as observed in flux ratios which varied from ~0.7 to 1 (gray bars in Figure 3A). Due to high standard error magnitudes, a clear dependence of R₂ on flux asymmetry was hard to identify with maximum flux asymmetry at R₂ = 9.2 Å and 11.2 Å and no flux asymmetry seen at R₂ = 6.2 Å.

FIGURE 3

Figure 3. Properties of heterotypic pores with varying inner mouth size of conical section (22.5-R₂-9.2-9.2). (A) LY fluxes rose in both directions with increasing inner mouth size (only R₂). Differences between J_het,₁ and J_het,₂ vary amongst the different cases. LY flux ratio variations appeared to be multimodal with increasing R₂ with maximum flux asymmetry at R₂ = 9.2 Å and 11.2 Å. (B) P_dX in cone to cylinder direction dropped near R₂ and maintained almost a uniform difference amongst the different inner pore mouth sizes. (C) P_dX in the cylinder to cone direction decreased in the left mouth and vestibule with increasing R₂. This tendency flipped by the end of the vestibule with P_dX levels increased with increase in R₂ in the right mouth. (D–F) Average force per LY particle in (D) left mouth, (E) vestibule, and (F) right mouth. (D) Increase in R₂ caused F_X magnitudes (-x direction) in the left mouth to gently decrease in the cone to cylinder direction, and gently rise in the cylinder to cone direction. (E) In the vestibule, F_X increased in the cone to cylinder direction, while decreased in the cylinder to cone direction with increase in R₂. (F) F_X magnitudes in the cone to cylinder direction were negative and varied randomly in the right mouth.

Particle probability profiles pointed to the vestibular mouth (at R₂ or end of conical mouth) as the main source of P_dX variation in the cone to cylinder direction (Figure 3B). This order of difference in magnitude due to varying R₂ was maintained in the remaining pore profile. The P_dX profiles in the cylinder to cone direction exhibited different variations (Figure 3C). While P_dX remained in the same order up to the vestibule mouth (or end of the cylinder mouth), it was highest for R₂ = 6.2 Å and lowest for R₂ = 12.2 Å in the vestibule region. By the end of the vestibule (smaller end of the conical mouth), this trend was inversed and P_dX dropped along x with lowest levels at R₂ = 6.2 Å and highest at R₂ = 12.2 Å in the region and the remaining right mouth (conical section).

F_X differences between the two directions were pronounced in the vestibule (Figure 3E) and right mouth (Figure 3F) with higher F_X against the flux (negative) in the cone to cylinder direction. This rose in the cone to cylinder direction with increasing R₂ in the vestibule (blue bars in Figure 3E). The variation in F_X in the right mouth was random. In the left mouth, F_X magnitude (in -x direction) in the cone to cylinder direction fell, while rising in the cylinder to cone direction as R₂ was increased (Figure 3D). In general, the net electrical force acting against the direction of flux was more in the cone to cylinder direction than in the cylinder to cone direction.

Varying Cylindrical Mouth Size (22.5-9.2-R₃-R₃)

Fluxes rose with increased cylinder mouth size in both directions (blue and orange bars in Figure 4A). Flux ratio, on the other hand, followed a bell curve with maximal at R₃ = 8.2 Å. Flux asymmetry appeared to be present only at cylinder mouth sizes between ~7.2 and 9.2 Å.

FIGURE 4

Figure 4. Properties of heterotypic pores with varying cylindrical mouth size (22.5-9.2-R₃-R₃). (A) LY fluxes increased in both directions as the cylinder section was widened. LY flux ratio rose and then dropped with increasing cylinder mouth size with maximal at R₃ = 8.2 Å. (B) P_dX levels remained in the same order up to the mouth of the vestibule in cone to cylinder direction for the different cylinder mouth sizes. Probability profiles almost ran parallel from there on for all cases with probability levels proportional to R₃. (C) P_dX levels dropped at the beginning of the cylindrical mouth (left mouth) and mimicked the pore profile up to the vestibule, after which they gently rose. Profiles for different cylinder sizes ran almost parallel. (D–F) Average force per LY particle in (D) left mouth, (E) vestibule, and (F) right mouth. (E) In the cone to cylinder direction, F_X in the left mouth was marginally affected by increase in the cylinder mouth size. In the cylinder to cone direction though, F_X dropped as the cylindrical mouth was widened. (E) In the vestibule, F_X gently dropped in the cone to cylinder direction while rising in the cylinder to cone direction with widening cylindrical mouth. (F) F_X magnitudes in the right mouth were larger in the cone to cylinder than cylinder to cone direction with magnitude dropping with increasing R₃. F_X was negative for all but R₃ = 6.2 Å.

There were no noticeable differences in P_dX along x up to the vestibule (or end of the conical mouth) with varying cylinder size in the cone to cylinder direction (Figure 4B). In the vestibule, P_dX was inversely proportional to the cylinder radius with smaller cylinders having P_dX profiles an order of magnitude higher than the widest cylinder. In the right mouth, P_dX levels dropped gently and were proportional to the cylinder radius with P_dX profiles parallel to each other. On the other hand, in cylinder to cone direction simulations, P_dX profiles (Figure 4C) varied from the beginning of the cylindrical mouth to the end of the pore with P_dX values almost maintaining contoured profiles to each other with P_dX levels proportional to the cylinder mouth size.

F_X was negative or marginal in all cases except for R₃ = 6.2 Å in the cylinder mouth (right mouth, Figure 4F) in the cone to cylinder direction. As R₃ was increased, F_X dropped in the cylinder to cone direction while remaining almost constant in the cone to cylinder direction. In the vestibule, F_X remained constant up to R₃ = 8.2 Å and then decreased with widening cylinder mouth in the cone to cylinder direction (blue bars in Figure 4E). A similar trend was seen in the cylinder to cone direction where F_X remained almost the same up to R₃ = 8.2 Å, but then rose with widening of the cylindrical mouth (orange bars in Figure 4E). In the right mouth, while F_X remained marginal in the cylinder to cone direction, it dropped with increasing R₃ from 7.2 to 12.2 Å (Figure 4F). We looked at differences in F_X amongst the two flux directions and found that the F_X difference reversed in the left mouth and vestibule after R₃ = 9.2 Å.

Varying Cylindrical Mouth Size in Pore with Wide Conical Mouth (25.9-12.2-R₃-R₃)

Results were similar to those with pore profile of 22.5-9.2-R₃-R₃ (Figure 5A). Fluxes rose with wider cylinder mouth and flux ratio was seen to be maximal at R₃ = 8.2 Å. Maximum LY flux asymmetry (0.49) was slightly higher than observed in 22.5-9.2-R₃-R₃ simulations (0.55).

FIGURE 5

Figure 5. Properties of heterotypic pores with wide conical mouth and varying cylindrical mouth size (25.9-12.2-R₃-R₃). (A) LY fluxes rose in both directions with increasing cylinder size. Fluxes in cone to cylinder direction were smaller than in cylinder to cone direction for all cylindrical section sizes. LY flux ratios varied nonlinearly with cylinder mouth size with maximal at R₃ = 8.2 Å. (B,C) P_dX profiles in simulations with concentration gradient in (B) cone to cylinder and (C) cylinder to cone direction. P_dX profile variations were similar to those seen in the 22.5-9.2-R₃-R₃ simulations. (D–F) Average force per LY particle in (D) left mouth, (E) vestibule, and (F) right mouth. F_X in the left mouth and vestibule in the cone to cylinder direction simulations remained constant with increasing cylinder size. In the left mouth (D), F_X magnitude in cylinder to cone direction varied nonlinearly with a maximal at R₃ = 7.2 Å. (F) F_X magnitudes in the right mouth were larger in the cone to cylinder than cylinder to cone direction with magnitude dropping with increasing R₃. F_X was negative for all but R₃ = 6.2 Å.

P_dX profiles (Figures 5B,C) and variations amongst different cylinder sizes were also similar to those for 22.5-9.2-R₃-R₃, with the profiles showing pronounced differences starting from the end of the vestibule in the cone to cylinder direction and beginning of the cylinder mouth in the cylinder to cone direction. Since a wider conical mouth was used, the drop in P_dX along x at the beginning and end of the conical mouth was smaller than seen in 22.5-9.2-R₃-R₃. (Figure 5C vs. Figure 4C).

While P_dX variations were similar to those in the previous pore-shape cases (22.5-9.2-R₃-R₃), F_X variations were notably different. F_X in the cone to cylinder direction remained almost constant in the left mouth (Figure 5D) and vestibule (Figure 5E) with increasing cylinder mouth size. In the right mouth, F_X was positive at R₃ = 6.2 Å, then dropped to negative at R₃ = 7.2 Å, reaching minimal at R₃ = 8.2 Å and then rising again up to R₃ = 12.2 Å (Figure 5F). On the other hand, F_X magnitudes in the cylinder to cone direction were marginal in the vestibule and right mouth. In the left mouth, it increased up to R₃ = 7.2 Å and then dropped, while being greater in magnitude to its corresponding cone to cylinder direction up to R₃ = 11.2 Å (Figure 5D). While in the left mouth, the force against the flux (negative F_X) was higher in the cylinder to cone direction (R₃ = 6.2 to 10.2 Å), it was higher in the cone to cylinder direction in the vestibule and right mouth (R₃ = 7.2 to 12.2 Å).

From all different heterotypic shapes studied, a pore-profile of 25.9-12.2-8.2-8.2 yielded maximum LY flux asymmetry with flux ratio of ~0.49. Flux ratio of the corresponding homotypic pores of the left and right hemichannels was ~9.18.

Influence of Particle Size on Flux Asymmetry in Heterotypic Pores

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).

FIGURE 6

Figure 6. Role of particle size on flux asymmetry in heterotypic pores. Simulation results in heterotypic pore with profile 25.9-12.2-8.2-8.2 with constant concentration gradient, 2 mV transjunctional voltage and q_p = 2e⁻. (A) Particle fluxes dropped nonlinearly as particle radius was increased from 0.8 to 6.6 Å. Flux ratio values (gray bars) ranged from ~0.5 to 0.8. Maximum asymmetry of 0.51 was for particle size of 4.9 Å. (B,C) P_dX levels remain in the same order in cell 1. P_dX profiles were similar to the pore profiles (log 10 scale). (D–H) F_X in the different sections of the cell-pore-cell system. (D) F_X values in cell 1 were marginal compared to other sections. In general, F_X magnitudes dropped in all sections of the pore with increasing particle size. (E) In the left mouth, F_X dropped with increasing particle sizes with F_X in cylinder to cone direction were higher than the cone to cylinder direction for all cases. F_X for r_p = 6 Å in the cylinder to cone direction was unexpectedly high. (F,G) In the vestibule and right mouth, F_X was negative and distinctly higher in the cone to cylinder direction than in the cylinder to cone direction. (H) F_X magnitudes in cell 2 were larger in the cylinder to cone direction in all cases.

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).

Influence of Particle Charge on Flux Asymmetry in Heterotypic Pores

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).

FIGURE 7

Figure 7. Role of particle charge on flux asymmetry in heterotypic pores. Simulation results in heterotypic pore with profile 25.9-12.2-8.2-8.2 with constant concentration gradient, 2 mV transjunctional voltage and r_p = 4.9 Å. (A) Particle fluxes decreased with increasing particle charge. Flux ratio (gray bars) also rose up to q_p = 2e⁻. (B,C) The P_dX profiles were similar in shape but scaled down for high particle charge. (D–H) F_X variation patterns for different q_p were similar in the different sections, but the differences between the 2 directions were magnified with increasing q_p. (D) F_X magnitudes in cell 1 were marginal compared to F_X in other sections. (E) In the left mouth, F_X was negative in cylinder to cone direction and larger in magnitude. (F–G) In the vestibule and right mouth, F_X were larger in the cone to cylinder direction. (H) F_X in cell 2 were small and differences amongst the two directions were marginal.

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.

Discussion

Flux Asymmetry and Heterotypic Pores That Best Fit Cx43-Cx45 In vitro Data

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).

Flux Asymmetry Established By Pore Shape While Scaled by Particle Charge and Size

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.

Predicting Flux and Flux Asymmetry in Heterotypic Pores from Sectional Particle Probabilities and Electrical Forces

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.

Possible Role of Heterotypic Cx43-Cx45 Gap Junctions in Cardiac Disease

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.

Limitations

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.

Conclusions

Our investigations shed light on the role of pore shape and particle charge in observed flux asymmetry of large molecules in inactive heterotypic pores. In particular, we identified pore shapes that match in vitro data to a considerable degree. We found the pore shape of radius profile 25.9-12.2-8.2-8.2 produced the maximal flux asymmetry of 0.49. The pore shape came out as the key regulator of flux asymmetry. Flux asymmetry existed irrespective of particle size. Though particle charge was vital, it merely scaled the asymmetry up to a certain magnitude. It is important to note that such asymmetry is possible only in this nano-scale and it will disappear on enlarging the pore. The studies did not reveal a simple relationship for explaining flux asymmetry based on pore-shape, particle charge and size, but we identified factors that aid in predicting such events.

Our work provides a solid foundation to investigate flux symmetry of biologically relevant molecules like ATP and IP₃, which have major metabolic consequences to cell function. We envision that the conclusions from our in silico results will inspire future in vitro studies to demonstrate the role of such preferential fluxes in physiological and pathological conditions of cardiac or other tissues.

Author Contributions

AM developed the simulation model, managed and processed the data and prepared the results. AM prepared the initial draft of the manuscript and implemented corrections and suggestions toward the final draft of the manuscript. APM was the principal investigator and supervised in designing and executing the research. FBS and APM provided advice during data analysis and reviewed the manuscript. FBS also provided computational resources for this research.

Conflict of Interest Statement

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.

Acknowledgments

We acknowledge funding for this work from the Nora Eccles Treadwell Foundation (APM) and NIH R01 R01-63969 (APM). We would like to thank Haonan Yang, Nancy Allen and Dr. Phil Ershler for their computer support at CVRTI.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2017.00206/full#supplementary-material

References

Allen, T. W., Hoyles, M., Kuyucak, S., and Chung, S. H. (1999). Molecular and Brownian dynamics study of ion selectivity and conductivity in the potassium channel. Chem. Phys. Lett. 313, 358–365. doi: 10.1016/S0009-2614(99)01004-0

CrossRef Full Text | Google Scholar

Ayad, W. A., Locke, D., Koreen, I. V., and Harris, A. L. (2006). Heteromeric, but not homomeric, connexin channels are selectively permeable to inositol phosphates. J. Biol. Chem. 281, 16727–16739. doi: 10.1074/jbc.M600136200

CrossRef Full Text | Google Scholar

Bedner, P., Niessen, H., Odermatt, B., Kretz, M., Willecke, K., and Harz, H. (2006). Selective permeability of different connexin channels to the second messenger cyclic AMP. J. Biol. Chem. 281, 6673–6681. doi: 10.1074/jbc.M511235200

PubMed Abstract | CrossRef Full Text | Google Scholar

Bedner, P., Niessen, H., Odermatt, B., Willecke, K., and Harz, H. (2003). A method to determine the relative cAMP permeability of connexin channels. Exp. Cell Res. 291, 25–35. doi: 10.1016/S0014-4827(03)00323-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Beltramello, M., Bicego, M., Piazza, V., Ciubotaru, C. D., Mammano, F., and D'Andrea, P. (2003). Permeability and gating properties of human connexins 26 and 30 expressed in HeLa cells. Biochem. Biophys. Res. Commun. 305, 1024–1033. doi: 10.1016/S0006-291X(03)00868-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennett, B. C., Purdy, M. D., Baker, K. A., Acharya, C., McIntire, W. E., Stevens, R. C., et al. (2016). An electrostatic mechanism for Ca(2+)-mediated regulation of gap junction channels. Nat. Commun. 7:8770. doi: 10.1038/ncomms9770

PubMed Abstract | CrossRef Full Text | Google Scholar

Bevans, C. G., Kordel, M., Rhee, S. K., and Harris, A. L. (1998). Isoform composition of connexin channels determines selectivity among second messengers and uncharged molecules. J. Biol. Chem. 273, 2808–2816. doi: 10.1074/jbc.273.5.2808

PubMed Abstract | CrossRef Full Text | Google Scholar

Beyer, E. C., Gemel, J., Martínez, A. D., Berthoud, V. M., Valiunas, V., Moreno, A. P., et al. (2001). Heteromeric mixing of connexins: compatibility of partners and functional consequences. Cell Commun. Adhes. 8, 199–204. doi: 10.3109/15419060109080723

PubMed Abstract | CrossRef Full Text | Google Scholar

Bukauskas, F. F., Kreuzberg, M. M., Rackauskas, M., Bukauskiene, A., Bennett, M. V. L., Verselis, V. K., et al. (2006). Properties of mouse connexin 30.2 and human connexin 31.9 hemichannels: implications for atrioventricular conduction in the heart. Proc. Natl. Acad. Sci. U.S.A. 103, 9726–9731. doi: 10.1073/pnas.0603372103

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, F., Eckert, R., Elfgang, C., Nitsche, J. M., Snyder, S. A., Hülser, D. F., et al. (1998). A quantitative analysis of connexin-specific permeability differences of gap junctions expressed in HeLa transfectants and Xenopus oocytes. J. Cell Sci. 111, 31–43. Available online at: http://jcs.biologists.org/content/111/1/31.article-info

PubMed Abstract | Google Scholar

Chung, S. H., Allen, T. W., and Kuyucak, S. (2002). Modeling diverse range of potassium channels with Brownian dynamics. Biophys. J. 83, 263–277. doi: 10.1016/S0006-3495(02)75167-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Corry, B., Allen, T. W., Kuyucak, S., and Chung, S. H. (2001). Mechanisms of permeation and selectivity in calcium channels. Biophys. J. 80, 195–214. doi: 10.1016/S0006-3495(01)76007-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, L., Liu, X., Li, H., Vertel, B. M., and Ebihara, L. (2006). Role of the N-terminus in permeability of chicken connexin45.6 gap junctional channels. J. Physiol. 576, 787–799. doi: 10.1113/jphysiol.2006.113837

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckert, R. (2006). Gap-junctional single-channel permeability for fluorescent tracers in mammalian cell cultures. Biophys. J. 91, 565–579. doi: 10.1529/biophysj.105.072306

PubMed Abstract | CrossRef Full Text | Google Scholar

Ek-Vitorín, J. F., and Burt, J. M. (2005). Quantification of gap junction selectivity. Am. J. Physiol. Cell Physiol. 289, C1535–C1546. doi: 10.1152/ajpcell.00182.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ek-Vitorín, J. F., King, T. J., Heyman, N. S., Lampe, P. D., and Burt, J. M. (2006). Selectivity of connexin 43 channels is regulated through protein kinase C – dependent phosphorylation. Circ Res 98, 1498–1505. doi: 10.1161/01.RES.0000227572.45891.2c

PubMed Abstract | CrossRef Full Text | Google Scholar

Elenes, S., Martínez, A. D., Delmar, M., Beyer, E. C., and Moreno, A. P. (2001). Heterotypic docking of Cx43 and Cx45 connexons blocks fast voltage gating of Cx43. Biophys. J. 81, 1406–1418. doi: 10.1016/S0006-3495(01)75796-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Escalona, Y., Garate, J. A., Araya-Secchi, R., Huynh, T., Zhou, R., and Perez-Acle, T. (2016). Exploring the membrane potential of simple dual-membrane systems as models for gap-junction channels. Biophys. J. 110, 2678–2688. doi: 10.1016/j.bpj.2016.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Finkelstein, A. (1994). Gap junctions and intercellular communications. Science 265, 1017–1018. doi: 10.1126/science.265.5175.1017-b

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaudesius, G., Miragoli, M., Thomas, S. P., and Rohr, S. (2003). Coupling of cardiac electrical activity over extended distances by fibroblasts of cardiac origin. Circ. Res. 93, 421–428. doi: 10.1161/01.RES.0000089258.40661.0C

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberg, G. S., Lampe, P. D., and Nicholson, B. J. (1999). Selective transfer of endogenous metabolites through gap junctions composed of different connexins. Nat. Cell Biol. 1, 457–459. doi: 10.1038/15693

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberg, G. S., Moreno, A. P., and Lampe, P. D. (2002). Gap junctions between cells expressing connexin 43 or 32 show inverse permselectivity to adenosine and ATP. J. Biol. Chem. 277, 36725–36730. doi: 10.1074/jbc.M109797200

PubMed Abstract | CrossRef Full Text | Google Scholar

Söhl, G., and Willecke, K. (2004). Gap junctions and the connexin protein family. Cardiovasc. Res. 62, 228–232. doi: 10.1016/j.cardiores.2003.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, A. L. (2007). Connexin channel permeability to cytoplasmic molecules. Prog. Biophys. Mol. Biol. 94, 120–143. doi: 10.1016/j.pbiomolbio.2007.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Heyman, N. S., and Burt, J. M. (2008). Hindered diffusion through an aqueous pore describes invariant dye selectivity of Cx43 junctions. Biophys. J. 94, 840–854. doi: 10.1529/biophysj.107.115634

PubMed Abstract | CrossRef Full Text | Google Scholar

Im, W., and Roux, B. (2002). Ion permeation and selectivity of OmpF porin: a theoretical study based on molecular dynamics, Brownian dynamics, and continuum electrodiffusion theory. J. Mol. Biol. 322, 851–869. doi: 10.1016/S0022-2836(02)00778-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Jongsma, H. J., and Wilders, R. (2000). Gap junctions in cardiovascular disease. Circ. Res. 86, 1193–1197. doi: 10.1161/01.RES.86.12.1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanaporis, G., Brink, P. R., and Valiunas, V. (2013). Gap junction permeability : selectivity for anionic and cationic probes. Am. J. Physiol. Cell Physiol. 300, C600–C609. doi: 10.1152/ajpcell.00316.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelsell, D. P., Dunlop, J., Stevens, H. P., Lench, N. J., Liang, J. N., Parry, G., et al. (1997). Connexin 26 mutations in hereditary non-syndromic sensorineural deafness. Nature 387, 80–83. doi: 10.1038/387080a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Krenning, G., Zeisberg, E. M., and Kalluri, R. (2010). The origin of fibroblasts and mechanism of cardiac fibrosis. J. Cell. Physiol. 225, 631–637. doi: 10.1002/jcp.22322

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwon, T., Harris, A. L., Rossi, A., and Bargiello, T. A. (2011). Molecular dynamics simulations of the Cx26 hemichannel: evaluation of structural models with Brownian dynamics. J. Gen. Physiol. 138, 475–493. doi: 10.1085/jgp.201110679

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, K. Il., Rui, H., Pastor, R. W., and Im, W. (2011). Brownian dynamics simulations of ion transport through the VDAC. Biophys. J. 100, 611–619. doi: 10.1016/j.bpj.2010.12.3708

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, X., Xu, Q., and Veenstra, R. D. (2014). Functional formation of heterotypic gap junction channels by connexins-40 and -43. Channels 8, 433–443. doi: 10.4161/19336950.2014.949188

PubMed Abstract | CrossRef Full Text | Google Scholar

Locke, D., Stein, T., Davies, C., Morris, J., Harris, A. L., Evans, W. H., et al. (2004). Altered permeability and modulatory character of connexin channels during mammary gland development. Exp. Cell Res. 298, 643–660. doi: 10.1016/j.yexcr.2004.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeda, S., Nakagawa, S., Suga, M., Yamashita, E., Oshima, A., Fujiyoshi, Y., et al. (2009). Structure of the connexin 26 gap junction channel at 3.5 A resolution. Nature 458, 597–602. doi: 10.1038/nature07869

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez, A. D., Hayrapetyan, V., Moreno, A. P., and Beyer, E. C. (2002). Connexin43 and Connexin45 form heteromeric gap junction channels in which individual components determine permeability and regulation. Circulation 90, 1100–1107. doi: 10.1161/01.RES.0000019580.64013.31

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondal, A., Appadurai, D. A., Akoum, N. W., Sachse, F. B., and Moreno, A. P. (2017). Computational simulations of asymmetric fluxes of large molecules through gap junction channel pores. J. Theor. Biol. 412, 61–73. doi: 10.1016/j.jtbi.2016.08.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreno, A. P., Hayrapetyan, V., Zhong, G., Martínez, A. D., and Beyer, E. C. (2004). “Homomeric and heteromeric gap junctions,” in Cardiac Electrophysiology, eds D. P. Zipes and J. Jalife (Philadelphia, PA: Saunders), 120–126. doi: 10.1016/B0-7216-0323-8/50017-8

CrossRef Full Text

Niessen, H., Harz, H., Bedner, P., Kramer, K., and Willecke, K. (2000). Selective permeability of different connexin channels to the second messenger inositol 1,4,5-trisphosphate. J. Cell Sci. 113(Pt 8), 1365–1372. Available online at: http://jcs.biologists.org/content/113/8/1365.short

PubMed Abstract | Google Scholar

Rackauskas, M., Verselis, V. K., and Bukauskas, F. F. (2007). Permeability of homotypic and heterotypic gap junction channels formed of cardiac connexins mCx30.2, Cx40, Cx43, and Cx45. Am. J. Physiol. Heart Circ. Physiol. 293, H1729–H1736. doi: 10.1152/ajpheart.00234.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, S. R., Hampson, E. C., Munro, M. N., and Vaney, D. I. (1993). Unidirectional coupling of gap junctions between neuroglia. Science 262, 1072–1074. doi: 10.1126/science.8093125

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohr, S. (2009). Myofibroblasts in diseased hearts: new players in cardiac arrhythmias? Hear. Rhythm 6, 848–856. doi: 10.1016/j.hrthm.2009.02.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Schirmer, T., and Phale, P. S. (1999). Brownian dynamics simulation of ion flow through porin channels. J. Mol. Biol. 294, 1159–1167. doi: 10.1006/jmbi.1999.3326

PubMed Abstract | CrossRef Full Text | Google Scholar

Severs, N. J., Rothery, S., Dupont, E., Coppen, S. R., Yeh, H. I., Ko, Y. S., et al. (2001). Immunocytochemical analysis of connexin expression in the healthy and diseased cardiovascular system. Microsc. Res. Tech. 52, 301–322. doi: 10.1002/1097-0029(20010201)52:3<301::AID-JEMT1015>3.0.CO;2-Q

PubMed Abstract | CrossRef Full Text | Google Scholar

Spray, D. C., and Dermietzel, R. (1995). X-linked dominant Charcot—Marie—Tooth disease and other potential gap-junction diseases of the nervous system. Trends Neurosci. 18, 256–262. doi: 10.1016/0166-2236(95)80004-L

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlenberg, B., Schuelke, M., Rüschendorf, F., Ruf, N., Kaindl, A. M., Henneke, M., et al. (2004). Mutations in the gene encoding gap junction protein α12 (connexin 46.6) cause Pelizaeus-Merzbacher-like disease. Am. J. Hum. Genet. 75, 251–260. doi: 10.1086/422763

CrossRef Full Text | Google Scholar

Unger, V. M., Kumar, N. M., Gilula, N. B., and Yeager, M. (1999). Three-dimensional structure of a recombinant gap junction membrane channel. Science 283, 1176–1180. doi: 10.1126/science.283.5405.1176

PubMed Abstract | CrossRef Full Text | Google Scholar

Valiunas, V., Beyer, E. C., and Brink, P. R. (2002). Cardiac gap junction channels show quantitative differences in selectivity. Circ. Res. 91, 104–111. doi: 10.1161/01.RES.0000025638.24255.AA

PubMed Abstract | CrossRef Full Text | Google Scholar

Valiunas, V., Polosina, Y. Y., Miller, H., Potapova, I. A., Valiuniene, L., Doronin, S., et al. (2005). Connexin-specific cell-to-cell transfer of short interfering RNA by gap junctions. J. Physiol. 568, 459–468. doi: 10.1113/jphysiol.2005.090985

PubMed Abstract | CrossRef Full Text | Google Scholar

van Kempen, M. J. A., Velde, I. T., Wessels, A., Oosthoek, P. W., Gros, D., Jongsma, H. J., et al. (1995). Differential connexin distribution accommodates cardiac function in different species. Microsc. Res. Tech. 31, 420–436. doi: 10.1002/jemt.1070310511

PubMed Abstract | CrossRef Full Text | Google Scholar

Veenstra, R. D. (1996). Size and selectivity of gap junction channels formed from different connexins. J. Bioenerg. Biomembr. 28, 327–337. doi: 10.1007/BF02110109

PubMed Abstract | CrossRef Full Text | Google Scholar

Veenstra, R. D. (2001). “Determining ionic permeabilities of gap junction channels,” in Connexin Methods and Protocols, eds R. Bruzzone and C. Giaume (New Jersey, NJ: Humana Press), 293–311. doi: 10.1385/1-59259-043-8:293

CrossRef Full Text | Google Scholar

Veenstra, R. D., Wang, H. Z., Beblo, D. A., Chilton, M. G., Harris, A. L., Beyer, E. C., et al. (1995). Selectivity of connexin-specific gap junctions does not correlate with channel conductance. Circ. Res. 77, 1156–1165.

Google Scholar

Veenstra, R. D., Wang, H. Z., Beyer, E. C., and Brink, P. R. (1994). Selective dye and ionic permeability of gap junction channels formed by connexin45. Circ. Res. 75, 483–490. doi: 10.1161/01.RES.75.3.483

PubMed Abstract | CrossRef Full Text | Google Scholar

Verselis, V., White, R. L., Spray, D. C., and Bennett, M. V. (1986). Gap junctional conductance and permeability are linearly related. Science 234, 461–464. doi: 10.1126/science.3489990

PubMed Abstract | CrossRef Full Text | Google Scholar

Vora, T., Corry, B., and Chung, S. H. (2008). Brownian dynamics study of flux ratios in sodium channels. Eur. Biophys. J. 38, 45–52. doi: 10.1007/s00249-008-0353-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Vozzi, C., Dupont, E., Coppen, S. R., Yeh, H.-I., and Severs, N. J. (1999). Chamber-related differences in connexin expression in the human heart. J. Mol. Cell. Cardiol. 31, 991–1003. doi: 10.1006/jmcc.1999.0937

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, P. A., Chang, H. C., Spaeth, K. E., Nitsche, J. M., and Nicholson, B. J. (2004). The permeability of gap junction channels to probes of different size is dependent on connexin composition and permeant-pore affinities. Biophys. J. 87, 958–973. doi: 10.1529/biophysj.103.036350

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, J. H., Liu, C. Y., Tang, B. S., Pan, Q., Huang, L., Dai, H. P., et al. (1998). Mutations in the gene encoding gap junction protein beta-3 associated with autosomal dominant hearing impairment. Nat. Genet. 20, 370–373. doi: 10.1038/3845

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamada, K. A., Rogers, J. G., Sundset, R., Steinberg, T. H., and Saffitz, J. E. (2003). Up-regulation of Connexin45 in heart failure. J. Cardiovasc. Electrophysiol. 14, 1205–1212. doi: 10.1046/j.1540-8167.2003.03276.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yum, S. W., Zhang, J., Valiunas, V., Kanaporis, G., Brink, P. R., White, T. W., et al. (2007). Human connexin26 and connexin30 form functional heteromeric and heterotypic channels. Am. J. Physiol. Cell Physiol. 293, C1032–C1048. doi: 10.1152/ajpcell.00011.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Zonta, F., Buratto, D., Cassini, C., Bortolozzi, M., and Mammano, F. (2014). Molecular dynamics simulations highlight structural and functional alterations in deafness related M34T mutation of connexin 26. Front. Physiol. 5:85. doi: 10.3389/fphys.2014.00085

PubMed Abstract | CrossRef Full Text | Google Scholar