Accessing gap-junction channel structure-function relationships through molecular modeling and simulations
BMC Cell Biology volume 18, Article number: 5 (2017)
Gap junction channels (GJCs) are massive protein channels connecting the cytoplasm of adjacent cells. These channels allow intercellular transfer of molecules up to ~1 kDa, including water, ions and other metabolites. Unveiling structure-function relationships coded into the molecular architecture of these channels is necessary to gain insight on their vast biological function including electrical synapse, inflammation, development and tissular homeostasis. From early works, computational methods have been critical to analyze and interpret experimental observations. Upon the availability of crystallographic structures, molecular modeling and simulations have become a valuable tool to assess structure-function relationships in GJCs. Modeling different connexin isoforms, simulating the transport process, and exploring molecular variants, have provided new hypotheses and out-of-the-box approaches to the study of these important channels.
Here, we review foundational structural studies and recent developments on GJCs using molecular modeling and simulation techniques, highlighting the methods and the cross-talk with experimental evidence.
Results and discussion
By comparing results obtained by molecular modeling and simulations techniques with structural and functional information obtained from both recent literature and structural databases, we provide a critical assesment of structure-function relationships that can be obtained from the junction between theoretical and experimental evidence.
Gap junctions (GJs) are regions of cellular membranes in which transmembrane proteins belonging to adjacent cells are in close contact, thereby forming hydrophilic dual-membrane channels. These channels allow the exchange of nutrients, metabolites, ions and small molecules up to ~1 kDa. GJ channels (GJCs) are formed by the end-to-end docking of the extracellular portion of two hemichannels (HCs) or connexons  each HC being composed of an hexagonal array of connexins (Cx) protomers . GJCs have crucial roles in many processes including differentiation, neuronal activity, development, immune responses and cell synchronization. Moreover, several human diseases are caused by mutations in connexins, including neurodegenerative diseases, skin diseases, deafness and developmental abnormalities .
From rough to fine: the early ages of GJC structure
The discovery of GJs began with the seminal work of Robertson who described them as regular and hexagonal lattices filling the gap between the cellular membranes of adjacent cells [4, 5]. Benedetti and Emmelot  described gap junctions as being particularly abundant in regions related to cellular communication and through electron microscopy of rat liver cells, demonstrated their icosahedral structure and hexameric symmetry (Fig. 1, panels A to E). In 1967 Revel and Karnovsky  presented findings that led to the term “Gap Junctions”. Payton et al.,  conducted experiments demonstrating that GJ transport molecules. Through the next decade, the Goodenough lab made substantial progress by identifying the chemical components of GJs as common lipids and one specific (yet unidentified) protein . Further studies led to the characterization of the constituent protein of gap junctions, named connexin . It was later proposed that one connexon crossed each junctional membrane forming an interconnecting channel from the cytoplasm of one cell to the cytoplasm of the other, spanning the 2 nm gap between the apposed cells . In this way, the term connexon was used to represent the half-channel (hemichannel) contributed by each cell to create the gap junction channel between cells.
By studying the structure of isolated GJs using high-resolution EM up to 18 Å, Unwin and Zampighi  demonstrated that observed hexameric structures formed regular lattices through the cellular membrane. Attention to the electronic density led to their proposal that GJC’s could exist in at least two different states; open and closed. Later, Unwin and Ennis , proposed that the GJC state was sensitive to the concentration of Ca2++, and that these divalent cations promoted the switching between the two states of the channel. Therefore, it was established for the first time that GJCs could not only transport solutes but regulate the permeability between cells.
A topological comparison between the HC of α and β connexins emerged from the work of Yeager and Gilula  who used anti-peptide antibodies directed to different sites of the protein sequence, together with cleavage by an endogenous protease and 2D-EM imaging. In doing so, these authors proposed a 2D-plot denoting the significant difference between the carboxyl terminal region of α (Cx43) and β (Cx32) connexins. Importantly, despite this sequence divergence, both connexin subfamilies share their quaternary structure, denoting a similar HC architecture. In 1994, Zhang and Nicholson , used a similar approach to determine the topological model of Cx26, demonstrating that it is consistent with the previously deduced topology for Cx32 and Cx43. Shortly afterward, in 1997, Unger and colleagues  achieved a structural projection of a Cx43 GJC at 7Å resolution revealing that a ring of transmembrane (TM) alpha-helices flanked a central hydrophilic pore and a second ring of alpha-helices was in close contact with the membrane lipids.
The first 3D structure of a GJC was determined by Unger and colleagues in 1999  (Fig. 2a) using Cx43. Using electron crystallography at 7.5 Å resolution, they confirmed that the channel was formed by the end-to-end docking of two HC, exhibiting 24 internal electron densities that were consistent with an alpha-helical conformation of the four TM domains of each connexin (Fig. 2b). However, it was not possible to assign unequivocally the correspondence of each TM in sequence with the observed densities. This fundamental work was in agreement with most of the accumulated structural and functional studies  and provided a foundational structural source for gap junction biologists.
In the following years, many studies were designed to obtain more detailed structural information about GJCs and also to relate structural features with transport and gating processes. Methods involving electron cryomicroscopy (cryo-EM) , X-ray diffraction , atomic force microscopy (AFM) , computational methods [21, 22] and nuclear magnetic resonance (NMR)  combined with mutational, biochemical, and functional studies (reviewed in ) have provided a wealth of information about GJCs.
The devil is in the details: the race for a high resolution model
The use of computational algorithms in conjunction with the experimental evidence became a promising way of obtaining a more detailed structure of GJC. Fleishman and colleagues  provided the first model of the Cα atoms of a GJC based on an improved cryo-EM map and biochemical experiments (Fig. 2c). In this work, the main challenge was the assignment of the TM helices to the corresponding densities in the EM map. Fleishman used a mixed approach of sequence analysis and experimental data. In detail, the sequence analysis was based both on the conservation and side chain polarity of the residues in each TM, with the following constraints: charged and conserved residues at sites of interhelical interaction; less conserved residues facing the lipids or lumen of the pore. They also employed experimental data on accessibility of each region in GJC, revealed using SCAM (substituted-cysteine accessibility method). This technique was carried out by several groups on several GJC isoforms which not necessarily were in agreement, in fact studies on hemichannels identified the pore lining domain as TM1 [24, 25] while studies of GJC’s reported the pore-lining as TM3 . So, in the Fleishman work , the pore-lining regions were defined as TM1 and TM3 in a rigid-helix structure with a parallel orientation.
A key aspect of GJC structure is related to the interactions between the two connexons of apposing cells. The EM-maps gave little information about this region and the flexible nature of the extracellular domains. In 2007, Kovacs and colleagues  studied the extracellular region of GJC, modeling the interaction as a highly ordered β-barrel-like structure.
Oshima et al.  computed 3D maps from EM of a Cx26 M34A mutant. The M34A mutation was employed to create a more stable channel based on information that amino acid substitutions at position 34 result in channels that appear to be stabilized in a partially closed conformation . Oshima et al.  observed the general structure of Cx26 channel as similar to the Cx43 channel observed at a resolution of 7.5 Å . In addition, they observed a plug-like density in the pore lumen and proposed that the M34A mutation could generate a conformational change around the N-terminal domain, causing it to block transport, posing as a plausible, explanation for how alterations in function relate to structure. To explain the blocking mechanism, Oshima et al.  superimposed the Fleishman et al.  Cα model on their EM-map and argued that the flexibility of the N-terminal domain would allow it to enter to the pore acting as a plug only for M34A mutant due to the smaller side chain of alanine in position 34 . This was consistent with biochemical data showing that sidechain length at position 34 is a key determinant of channel function .
In 2008, Pantano et al. , built an all-atom model of a Cx32 connexin based on TM assignments of Fleishman et al. . In this exhaustive work they rebuilt the side chains using molecular dynamics they fitted the Cα structure obtaining a stable system embedded in a palmitoyl-oleoyl-phosphatidyl-choline (POPC) bilayer. The method introduced in this paper, to reconstruct a protein starting from the protein backbone, was validated on a different channel (the KcsA potassium channel, Pantano et al., ). To date, several protein structure have been modeled using this approach .
A light in the shadow: the crystallographic structure of the human Cx26
Less than a year after the publication of the all-atom model , a 3D structure of human connexin 26 (hCx26) was determined at a resolution of 3.5 Å by X-ray crystallography, providing key structural information about connexins GJCs, including identification of the pore lining regions of hCx26 [19, 32] (Fig. 3). The crystal structure demonstrates that each protomer is composed by four transmembrane helices (namely TM1 to TM4), two extracellular loops (E1, E2) and a N-terminal helix (NTH), forming a typical four-helical bundle in which any pair of adjacent helices are antiparallel. The arrangement of the tertiary structure of hCx26 is as follows: NTH, TM1 and TM2 face the pore; TM3 and TM4 are located within the perimeter of the hemichannel facing the membrane lipids, with E1 and E2 are facing the extracellular environment. Thus, this crystal structure confirmed that TM1, NTH, E1 and, to a lesser extent TM2, are the major pore-lining regions of each connexin subunit, discarding any role of TM3 in this function. Regarding inter-protomer interactions, the Cx26 crystal structure shows that they are mostly located in the extracellular half of the transmembrane helices TM2 and TM4, and in the extracellular loops. Residues included in the inter-protomer interactions, as well as those involved in the intra-protomer interactions, are highly conserved among several members of the connexin family, suggesting that there is a conservation of both the protomer folding and the oligomerization process to form the hemichannel, within the connexin family members . Moreover, this structure shed light on an important unresolved issue: the TM helix assignment. This issue arose because the connecting loops between the transmembrane alpha helices were not revealed in the first map published by Unger et al. . In fact, reference to Fig. 4 of Unger et al.  noted that with an in-plane resolution of 7.5 Å, there was ambiguity in assigning the molecular boundary of the connexin subunits. The helix assignments in the model of Fleishman et al. , which served as a starting point for several other modes [18, 30] are inconsistent with the arrangement of TM helices in the atomic model of Maeda et al. .
Following this work, a deeper analysis of the crystallographic structure  led to heightened focus on the role of the amino terminal helix (NTH) as a key element in channel blockage and gating. New structural information revealed details including a circular hydrogen bond network between the Asp2 side chain and the NTH main chain. Moreover, Trp3 from the NTH and other residues form a hydrophobic patch along the pore lining side of TM1 of the counter-clockwise adjacent subunit. This hydrophobic patch comprises residues M34, V36 and V37. Remarkably, a set of deafness-associated mutations has been mapped to this patch, such as V37I, M34A, M34T and A40V, suggesting that this zone is highly relevant for Cx26 function (reviewed in ).
One template for all: when the models make progress
The crystallographic hCx26 structure not only revealed the detailed structure of a GJC, it opened the door to the development of comparative models using this structure as a template. This is especially important because there are as many as 21 connexin isoforms expressed in different species and different tissues (Harris, ). Moreover, each six-protomer connexon can form homo- or heteromeric GJC creating an almost limitless space for structural and functional specializations [36, 37].
Comparative modeling takes advantage of a known 3D structure to build the structure of a different but related (at least in terms of sequence similarity) protein. Early work in the 1980’s [38, 39] and 1990’s [40, 41] demonstrated that protein structure is more likely conserved than sequence, so sequences with at least 40% identity could have very similar structures. Therefore, one sequence can be modeled using the structure of another if they share at least 40% sequence identity. In the case of membrane proteins, this statement is more general, because in transmembrane regions a high structural correlation occurs even when sequence identity is as low as 20% .
A key step in the comparative modeling procedure is the sequence alignment between the template structure and the target sequence. This alignment guides the progressive building of the main chain for the target sequence. Each software employs a distinct method (or a combinations of these) such as rigid-body assembly, segment matching or satisfaction of spatial restraints . Modern comparative models are developed using alignments that consider not only protein sequence but also some structural features, for instance the position of the secondary structure elements.
However some groups [32, 44] used a more manual approach to build their models of Cx32 GJC. They started from a multiple sequence alignment without considering structural features and manually corrected the mismatched residues on the Cx26 structure with the corresponding residue on the target Cx32. The side chain structure of these mismatched residues were built with the software COOT  and checked for spatial restrictions with the software Procheck . The model of Cx32 was used to generate hetero Cx32-Cx26 GJC and was used to explain their experimental observations related to channel docking. In their work they analyzed inter-connexon interactions in the extracellular space and they identified hydrogen bond networks on the E1-E1 and E2-E2 interfaces, inferring that the latter probably play a pivotal role in docking .
The above procedure was also applied to analyze the parahelix of Cx50 , a region between TM1 and EC1 that forms an imperfect 3.10 helix in the crystallographic structure [19, 48]. This region seems very stable in several simulations and structure-function studies revealed that it is involved in voltage gating [49, 50]. Tong et al.  analyzed the parahelix structure using comparative modelling focusing on charge-changing mutations in this region. They found that local surface electrostatic charges in this region play an important role in the ion permeation.
After main-chain modeling, further refinement can include addition of amino acid side chains. Many programs are capable of using a library of rotamers calculated from experimentally known structures, and selecting according to criteria such as steric clashes, energy considerations, and main chain structure dependency among others . This step is particularly important because the crystallographic structure derived by Maeda et al.,  lacks the start methionine (M1) and sidechains of several residues. The existence of loops not present in the template structure represents a challenge for most of modeling works as well, and it can be solved via searching for segments that fits on fixed-backbone endpoints, or by a conformational search (ab initio) or by a combination of both .
When the dance begins: the role of simulation in revealing molecular motions
During the last 20 years, computational techniques have become an indispensable tool for understanding of processes at the nanoscale level. The number of articles using comparative models of GJC are increasing in number, with many based on the atomic model of Cx26 . However, the intrinsically dynamic nature of GJC function – to transport molecules from one cell to another – presents a challenge to the analysis based on the static models and the static X-ray structure itself.
Among the methodologies belonging to the field of molecular simulation, molecular dynamics (MD) is widely used for the study and characterization of molecules at an atomic resolution . Starting from an experimentally determined structure, or even from a model , MD techniques allow the researcher to “capture” the natural motion of the molecular system and to monitor its dynamic behavior through time. MD deals with the numerical solution of the N-body problem raised by Isaac Newton in the 17th century, by elegantly combining statistical mechanics and classical physics, through the time-dependent integration of Newton’s equation of motion. These equations, even in the simplest systems, are of such complexity that the integration must be done with numerical methods over a huge quantity of discrete time steps, instead of being performed in an analytical, continuous fashion. At each time step, the coordinates of the atoms, are used to calculate the potential energy of the system (V) and the force. The latter is calculated using a molecular mechanics force field (FF) or potential energy function. A long series of these calculations allows the generation of a trajectory through phase space, defined by the three atomic spatial position and momentum arrays, as representative of a statistical mechanical ensemble of the molecular microstates. In this way, the reliability of MD simulations fully depends on two factors; a) the capacity to explore all regions of phase-space, also known as the sampling problem and b) the ability to faithfully reproduce the potential energy surface of the studied system, the scoring problem . The sampling problem has been partially alleviated by consistent software development, together with advances in specialized computer hardware, as evidenced by recent simulations in the order of hundreds of microseconds , moreover special enhanced-sampling techniques can be employed . On the other hand the scoring problem entirely depends on the quality of the FF, which are specifically tailored to a given set of molecular species. Regarding biomolecules, these are commonly termed biomolecular FF, with many “flavors” such as the GROMOS , OPLS , CHARMM , and AMBER  force-fields. Even though their parameterization philosophy differs, they are functionally similar. They generally include terms that describe major bonded (bonds, angles, and dihedral angles) and non-bonded (van der Waals and electrostatic) interactions. Parameters used for these energy terms derive from a combination of experimental data and quantum mechanical calculations, tuned to optimally reproduce the structure and vibrational modes of the molecular systems of interest, as well as their thermodynamic properties .
Normally the configurational part of a MD is of special interest, because it is possible to analyze atom movements and conformational changes. Thus, the trajectory of a MD i.e., the ensemble of snapshots of atomic coordinates in function of time, is critical for the GJC analysis. The trajectories in MD range from picoseconds to hundreds of nanoseconds, but in the case of GJC they should be long enough to allow such a massive system to adopt the conformational equilibrium required and overcome the sampling problem. The crystallographic structure in Maeda et al.  has ~9,800 atoms, the model with completed regions absent in the previous structure and the corresponding hydrogen atoms has nearly 30,000 atoms, moreover the model embedded in bilayer surrounded by solvent (water plus ions) has ~178,000 atoms. So the time needed to get the structure equilibrated from its “frozen” state in the initial structure to the desired temperature structure that mimic the experimental conditions, is near the microseconds time scale .
Over the last few decades, MD has been used in the study of a wide range of biological phenomena, such as protein folding, ion conduction, and muscle elasticity [63–66]. The rapid increase of the computational power and recent developments in simulation software have enabled MD studies of significantly larger systems, such as the ribosome in complex with a protein-conducting channel (2.7 million atoms)  and of much longer processes (hundreds of microseconds) . In particular, our research group has been using MD simulations since 2003, enhancing understanding of applied and pure research areas such as protein engineering , computer based drug design [70, 71] structure and function relationships of GPCRs [72, 73], protein-binding domains in plants , viral recognition and hypersensitivity response in plants , and water permeability properties of archaeal aquaporins . Despite the enormous success of MD protocols, an important gap remains between the time scales currently accessible by MD simulations and the length of manyl biological processes. As a matter of fact, rotational and translational motions of large-scale domains such as those involved in channel gating (microseconds to milliseconds or even longer), are still far from the reach of traditional MD methods. Furthermore, the intrinsic barriers that exist between important microstates hinders proper configurational sampling, indeed statistical mechanics precludes such barrier crossing when these are high above thermal energy (kBT), normally higher than 2.5 KJ/mol (0.6 Kcal/mol) at 300 K. To bridge this gap, special simulation procedures, termed enhanced sampling techniques have been developed to allow the exploration of these slow degrees of freedom [77, 78].
Even though the enhanced sampling techniques have partially alleviated the sampling problem, many biologically relevant processes are above the scope of current atomistic models, consequently simplifications, normally in the form of a drastic reduction of degrees of freedom are employed. In this way, a common approach is to remove unimportant degrees of freedom, normally solvent molecules employing the Brownian Dynamics (BD) formulation . In BD, the dynamics of molecules of interest is controlled by the Potential of Mean Force (PMF), an energy function that includes the averaged effect of the neglected particles, a friction factor and a stochastic term that mimics the effect of random collisions. Proper BD simulations, therefore require the careful parameterization of the PMF term, which can be obtained by rigorous free energy calculations in all-atom MD. BD simulations have been successfully applied in biomolecular simulations of protein crowding and ion-conduction in membrane channels under electro-osmotic gradients, among others [80–82]. In particular for BD simulations of ion-channels, the electrostatic term of potential equation is replaced by a linearized solution of the Poisson-Boltzmann equation. These calculations are numerically costly thus, it is common practice to simulate the channel as a rigid body, impeding the study of gating phenomena due to structural rearrangements. This is the case in the work of Kwon et al.  where BD is used to simulate the ion flux considering a fully rigid protein coupled to an ion bath via the Grand-canonical Monte Carlo (GCMC) technique. In GCMC , apart from the typical randomized translational, rotational and orientational moves typical of standard Canonical Monte Carlo simulations, random insertions (from a virtual reservoir) or deletions of particles (ions, in this case), via the Metropolis-Hastings algorithm are performed [84, 85]. The probability of insertion/deletion will also depend on the chemical potential of the studied species. Upon usage of 2 baths kept at different chemical potential values, a concentration gradient can be established, rendering an ionic flux, thus current, against the gradient. This mixed approximation allowed calculation of current versus voltage relationships. Otherwise the molecular dynamics imposes a methodological problem in the form of the periodic simulation condition.
The periodicity of a system in molecular simulations assures that any atom reaching one side of the simulation box can appear at the opposite side. This property renders a periodic system that can be considered infinite in any direction of the simulation space. Without periodicity, the atoms would continuously hit the “border” and this would influence some macroscopic properties of the system, like volume and pressure and would render simulations far from reality. The problem with transport-through-channel systems is that the solvent molecules (i.e. ions), prefer to move from one side of the bilayer to the other through the energetically free path the periodic box offers instead of through the channel (i.e. connexin pore).
A simple approximation to study the transport of molecules through channels is the application of an external force to push (or pull) the molecule or molecules to go through a specific path. For ion transport, this force could be applied to the entire system, in the form of a constant electric field parallel to the pore axis, giving rise to a voltage difference along the membrane. This external force can be calibrated to provide the necessary work for the molecules to move from one side to the other according to the direction and magnitude of the applied field. This is the case in the work of Zonta et al.,  where potentials of −80 mV and 80 mV were applied to study the binding of calcium ions to the extracellular region.
Another possibility of external force usage, is the so-called steered molecular dynamics (SDM) technique. This technique consists of applying an external force to a specific molecule to pull it into a certain direction. The pulling is not exerted directly over the molecule, instead is applied through a virtual spring with a certain stiffness, attached to the atom (or atoms) of interest. The force applied could be constant and generate variable velocity during the movement or be handled by the software to maintain a constant velocity in the selected direction. The latter is commonly used because the force necessary to maintain a constant velocity must counterbalance the forces exerted by the molecular environment (i.e. the pore internal walls) and this could be registered as a function of the distance of movement, in the form of force exerted along the pathway, in other words the work or PMF. With this information the free-energy landscape along a specific reaction coordinate could be analyzed (reviewed in ).
For GJC, SMD have been employed to explore the transport of chloride and potassium [86, 88], molecules such as calcein  and simple saccharides . In the case of atomic ion transport through gap junctions, a region within the pore opposes their permeation, forming an energetic barrier, the parahelix (PH) region. This region is the second narrowest in the pore, after the NTH region, and is located between TM1 and E1. Interestingly, the NTH region is more flexible and does not provide a high resistance to the passage of solutes, while the PH region is more stable and has a higher number of conserved charged residues, in which induce a very stable electrostatic network that can interact with the studied solutes . In the case of molecules, Luo et al.  argue that there is a difference in the energetics of transport for molecular solutes and atomic ions. Molecules have more degrees of freedom, generating substantial entropic barriers as the molecule potentially adopts a specific conformation in order to pass the narrow regions of the pore. On the other hand, for atomic ions, the main contributor to free energy of permeation should be enthalpic due to the lack of internal degrees of freedom of these species and the presence of charged residues within the channel. In this way, Luo et al.,  analyzed the transport free-energy of two aminopyridyl-labebeld maltosaccharides, one being permeant to GJC (maltose, a disaccharide) the other being impermeant to GJC (maltotriose, a trisaccharide). Briefly, they did not find significant differences in the PMF profiles along the pore pathway (see Fig. 4) and noted configurational entropy as the key to discrimination between solutes. The trisaccharide is bigger, thus more flexible (with more degrees of freedom) and it is therefore less probable that the molecule will adopt the correct conformation to pass through the narrow pore. The disaccharide needs less energy to overcome this conformational barrier and pass through the pore. Nevertheless, proper entropy calculations need to be performed to confirm this claim.
When generating structural models of HCs or GJCs, the explicit representation of lipid bilayers and water molecules provides a more realistic environment for the study of this complex system. Despite the limitations in simulation length due to the number of atoms, MD has been successfully employed to study a wide range of biomolecular systems and phenomena similar in complexity to GJCs . In simulation models employing MD simulations, the system consists of the completed HC embedded in a phospholipid bilayer plus solvent on each side of the bilayer, inside a simulation box of around 1331 nm3 [22, 48, 86] (Fig. 5). In all cases, the initial system is energy minimized to avoid steric clashes and after velocities assignment (from Maxwell-Boltzmann’s distributions) an equilibration dynamics and final production dynamics are performed. Due to the presence of a lipid bilayer and to allow the natural fluctuations of the membrane, this simulations must be performed with the NPT ensemble (constant particle number, pressure and temperature) via coupling the system to a virtual thermostat and barostat with anisotropic cell fluctuation, the latter is achieved by scaling velocities and positions, respectively.
Not all that glitters is gold: revealing the controversies still present after the crystallographic structure
Upon the release of an atomic model of Cx26 using X-ray crystallography , some of the controversies surrounding GJC structure seem resolved. However, this structure is a glimpse of a channel composed of one of over 20 different connexin proteins, captured in one state by the crystallographic process at a resolution which is not accurate enough (3.5 Å). In reality there are more unanswered questions than ever. Some of the key questions relate to dynamics properties, in this way, molecular simulation techniques arise as proper tools to generate hypotheses based on structure-dynamics-function relationships encoded in the molecular architecture of the Cx26 GJC.
One of the key questions related to the Cx26 structure  involves whether the channel was captured in the open or closed state. When the structure was published, it was speculated to represent an open configuration [19, 34] due to the unobstructed path from one side of the pore to the other. A proposed closed state of the channel was believed to have been presented in the low-resolution EM-map derived in Oshima et al.,  using the Cx26M34A mutant. This mutant which displayed a distinctive density in the cytoplasmic side of pore lumen, alleged to be the NTH  which was not observed in the high resolution atomic model. Therefore the configuration in the crystallographic structure was identified as being in an open state.
The above issue was partially addressed by the Bargiello lab using molecular dynamics simulations of the Cx26 hemichannel [22, 48] A goal of the first study involved adding and refining regions not well resolved in the atomic model  including the cytoplasmic loop (residues from 100 to 124), the C-terminus (residues from 218 to 226), the initial methionine residue and the side-chain associated with residues K15, S17 and S19. In their first simulation, Kwon et al.  used ion-flux calculations via grand-canonical monte-carlo brownian dynamics simulations (GCMC/BD) of this structure and showed that the pore region was too narrow to account for the experimentally observed currents. Under a transmembrane potential, these simulations predicted a marked inward current rectification with an almost full anionic selectivity, in full disagreement with experimental evidence . Maeda et al.  reported a minimal pore diameter of 14 Å but Kwon et al.  argued atom diameters, missing residues and sidechains such as Met1 were not considered in the calculation (Fig. 6a). The pore diameter after correction is on average, 10 Å . It is important to mention that the simulations employed herein did not account for any structural re-arrangement due to the protein being simulated as a rigid body. So this first simulation served to assess, the “openness” of the atomic model presented by Maeda et al. . A second simulation was performed on the completed structure, constructed by taking the atomic structure  revised to include side-chains and loops via comparative modelling. This new structure had an internal pore diameter less than 6 Å (Fig. 6b). It was further refined via long MD simulations leading to a slightly increased pore diameter (Fig. 6c) that matched much more closely the single-channel conductance derived from patch clamp experiments .
Another important feature lacking in the x-ray structure is related to co- and post-translational modifications (PTM) identified experimentally that change the net charge of several residues . These are the neutralizing acetylations of the N-terminus and internal lysines K15, K102, K105, K108, K112 and K116. Consequently, Kwon et al.,  carried out MD simulations including these modifications, which remarkably reproduced the experimentally derived current versus voltage (I/V) relations and cation selectivity. The authors concluded that the MD relaxation of the complete structure plus the charge neutralization derived from co- and post-translational modifications were essential for the development of a model of the open conformation that accounts for the cationic-selective nature of Cx26, and closely resembles the experimental single-channel I/V curves .
In a further study, the authors performed extensive MD simulations specifically to explore the structural determinants that stabilize the open conformation . Briefly, they explored the interaction network of the parahelix, which has been suggested as the loop permeability barrier . From energetic analyses of their simulations, it was found that an extensive van der Waals and electrostatic network stabilized the parahelix. Interestingly, this electrostatic network showed a certain level of correlated motion with an adjacent subunit, suggesting a cooperative effect. Consequently, the authors hypothesized that the disruption of this electrostatic network by an external voltage could allow the parahelix to protrude the channel, in a concerted mechanism that would eventually close the pore .
Additional disagreement with the atomic structure of Cx26  arose from molecular dynamics simulations of Zonta et al., . The Cx26 GJC structure was simulated using the atomic model of Maeda et al.  with completed regions, similar to the aforementioned work . The structure of Cx30 GJC was also modeled and analyzed. The authors confirmed the stability of TM regions, but observed a displacement of NTH towards the TM1, widening the pore lumen. They detected the hydrophobic interactions of Trp3 as described by Maeda  but not the hydrogen bonding network around Asp2. The transport process was also analyzed, and the electrostatic interaction networks around the parahelix was described as a potential barrier to ion passage, specifically the residues K41 and E49 in Cx26 and Cx30 respectively .
An old dancer comes to play: the effect of calcium on GJC function
The effect of calcium on GJC function was identified more than 30 years ago  but still presents a challenge for those aiming to understand its complex role in GJC function. Extracellular calcium is capable of blocking the voltage-dependent opening of hemichannels, and therefore it is expected to bind to specific amino acid residues facing the extracellular region of connexons [94–96].
A first approximation to model calcium binding came from Zonta et al. . Four Ca2+ atoms were introduced into the extracellular space, and different membrane potentials were simulated by applying an external electric field. With an electric field mimicking a membrane potential of −80 mV, these ions interacted with residues E42, D46, E47 and E50 located in the PH-EC1 region. The authors hypothesized that the binding of calcium ions in the PH region could cause the physical occlusion of the pore around this narrow sector. At zero field conditions (no electric field applied) the calcium ions lose the interaction and a field of +80 mV completely abolishes the interaction.
In further studies the nature of the binding was analyzed at a molecular level . Specifically the proposed PTM γ-carboxyl-glutamate on residues E42 and E47  was explored as a candidate for calcium ion binding. This PTM have been previously identified as a possible calcium coordinating moiety. This analysis was followed by quantum chemistry (QM) calculations at the density functional theory level. In QM, not only the atomic cores but the electronic degrees of freedom are simulated, thus rendering these type of techniques quite computationally costly, allowing the simulation of a only a few hundred atoms . QM simulations take a completely different theoretical approximation to simulate atomic motions. This technique takes a quantum mechanics approach, simulating not only the atomic cores but the electronic degrees of freedom. These more realistic considerations are computationally costly because they imply a much higher level of calculation, allowing simulation of small systems, or even sub-systems as in the aforementioned case of γ-carboxyl-glutamate binding of calcium ions . The analysis rendered a plausible mode of Ca2+ gating involving structural rearrangements induced via the ion’s interaction with PTM residues, and disruption of a proposed salt bridge network. The analysis rendered a plausible mode of binding of Ca2+.
In 2016, the structures of both a free and calcium-bound Cx26 GJC were determined at resolutions of 3.8 and 3.3 Å respectively . These structures (Fig. 7) bear strikingly similar to each other and also to the crystallographic structure of a Cx26 GJC . The Ca2+ ions are bound to the PH region, with five coordination atoms forming a square pyramidal geometry (Fig. 7c). These atoms belong to the carboxylate moiety of residues E42 and E47 plus the main-chain oxygen of G45 on the next protomer. This structure – with the exceptions of one salt bridge between E42 and R75 – is nearly identical to the calcium-free structure. Any conformational difference between them is localized to the Ca2+ binding site. The pore lining residues and inner diameter of the pore are not altered by Ca2+ binding. This discovery was unexpected because Ca2+ blocks the channel and it was presumed that a significant conformational change would be induced by the binding of calcium ions in this very narrow part of the pore. To uncover some dynamic behavior of the structures, MD simulations were performed confirming the stability of calcium binding, with the exception of a G45 main-chain interaction that was potentially unstable. The authors concluded that the mechanism of calcium block (particularly for potassium conduction) is more related to an electrostatic barrier imposed by the positive charge of calcium (Fig. 7) than a structural occlusion due to conformational changes . Recent work in our lab recently revealed that the charge arrangement within a GJC-generic model is crucial for its cationic selectivity .
Unfortunately the Cx26 structure doesn't resolve the position of all the residues, and the authors didn't discuss the proposed PTM in the same E42 and E47 residues identified as calcium binding sites [93, 97]. However this work is an interesting effort to unveil the structural role of calcium and the mechanism of conductance blockade, and more work in the field of molecular simulation will likely provide more progress in the area.
Finding the lost water: the discovery of the hCx26 IC pocket
About 3 years ago our group initiated molecular modelling studies of GJCs beginning with a study of Cx26. Molecular dynamics simulations of the completed structure of Cx26 HC (with PTMs) offered an excellent source of information about connexons and their transport processes. While studying the dynamic behavior of the NTH and its possible role in the slow gating process, a cavity was found between NTH, TM2 and TM3 in each of the six protomers (Fig. 8) . This cavity was filled with water molecules from the solvent (Fig. 8a), and phylogenetic analysis demonstrated that amino acid residues facing this cavity are highly conserved. Moreover, by conducting comparative molecular modelling the cavity - termed the IC pocket - was detected in all selected representatives of GJC (Pareja et al., unpublished work). The IC pocket is consistent with accessibility of side-chains revealed by substituted cysteine accessibility method (SCAM) that led to the identification of TM3 as the GJC pore lining . In SCAM analysis of Cx32 GJCs, accessible residues were identified in all TM domains including five residues in TM1 and six residues in TM3 (Fig. 8b). Of the five residues in TM1, only one (133C) faces the proposed water pocket. The others lie along the pore-facing region of TM1 (I30, F31, M34 and V35) facing the NT in the crystal structure of Cx26. Interestingly, these were accessible only in a proposed closed state of the channel. Of the six residues in TM3 that were accessible to sulfhydryl reagents, three correspond to residues facing the IC pocket (Y135, S138 and V139). Other accessible residues were located toward the extracellular end of TM3 (F141, L144 and F149) and a proposed mechanism of access to these residues has not been formally proposed.
Protein cavities, in particular those filled by water, play significant biological roles in transmembrane proteins. They are known to influence the activation/deactivation events of GPCRs , the ligand-binding processes of cannabinoid and beta-adrenergic receptors [73, 103, 104], intermolecular recognition events [105–107] and protein folding [107, 108]. By analyzing water dynamics it was determined that the water molecules within the IC pocket are different from bulk water, having a slower dynamic behavior and a longer residence time. Moreover, the number of water molecules in the IC pocket is correlated with the orientation of residues R143 and F29 (Fig. 8c). The former forms an electrostatic network that changes depending on the IC pocket state and could play an important role in channel activity. Interestingly, NTH orientation was related to changes in the electrostatic network involving R143 which influences IC water occupancy. The NTH orientation was measured as distance from protein center of mass, or as angle with pore axis . In vitro functional assessments of hCx26 and charge-altering mutants at position 143 were conducted by evaluating the uptake of molecular tracers through HCs. The results demonstrated that R143 is a key residue regulating conductance of hCx26 HCs to tracer molecules .
Many challenges to better understanding gap junction structure and regulation are related to the complex nature of these large intercellular channels. The complexity of a full GJC embedded in a two membrane system makes these channels difficult to purify and even more difficult to crystallize for x-ray diffraction experiments. This explains the scarcity of crystallographic structures - only three in 50 years. It also explains regions of poor resolution in the existing structural models [19, 99].
In functional analyses, the challenge of working with intercellular channels has been partially circumvented by analysis of HC. These channels are amenable to electrophysiology and reconstitution and have provided a wealth of information about channel gating, regulation and permeability . However, more work is needed to understand gating and regulation of GJC.
Other key questions are related to the mechanism by which disease-related mutations affect the channel, and this can be answered by comparing mutants with the Cx26 structures [19, 99]. Other questions include; what are the mechanisms of GJC gating, is there a structural correlation between cytoplasmic and extracellular portions of the pore and, what are the determinants of transport and permeability. Other important questions are related to the changes in sequence and structure that result in a diverse array of GJCs.
Although the current crystallographic structures are rigid, frozen snapshots of the channel, they have provided an astonishing amount of information and an irreplaceable starting point for simulation studies. These simulations unveil the dynamic behavior of GJC channels, further refine key regions in the structure, and model conformational changes that may occur in response to stimuli such as solute passage or applied voltage. Their use will undoubtedly lead to a better understanding of GJ channel structure and function, particularly when performed in conjunction with experimental techniques.
Molecular simulation, e.g. MD, can be considered as the in silico version of a conventional microscope – the (numerical) molecular microscope – with the advantage of providing a live view of the biomolecule at an atomic resolution. Time-averaged properties that are computed from an MD trajectory can be compared with macroscopic quantities, which are averages over time and multiple copies, that are measured experimentally. Furthermore, molecular simulation techniques offer the opportunity to explore systems under an unlimited number of artificial conditions that are often inaccessible experimentally. For instance, a residue can be changed to mimic the effect of an experimentally studied point mutation, or even neutralize the electrostatic forces generated by any given group of atoms.
The advance in processing power and availability of computational hardware offer more capacity for simulations and extend simulation times offering the ability to explore explicitly the conformational space. The use of advanced sampling techniques, for instance replica exchange with swarms of trajectories, may be the next step in obtaining more detailed models of GJC function.
Goodenough DA. In vitro formation of gap junction vesicles. J Cell Biol. 1976;68:220–31.
Harris AL. Emerging issues of connexin channels: biophysics fills the gap. Q Rev Biophys. 2001;34:325–472.
Saez JC, Berthoud VM, Branes MC, Martinez AD, Beyer EC. Plasma membrane channels formed by connexins: their regulation and functions. Physiol Rev. 2003;83:1359–400.
Robertson JD. The occurrence of a subunit pattern in the unit membranes of club endings in Mauthner cell synapses in goldfish brains. J Cell Biol. 1963;19:201–21.
Robertson JD, Bodenheimer TS, Stage DE. The ultrastructure of Mauthner cell synapses nodes in goldfish. J Cell Biol. 1963;19:159–99.
Benedetti EL, Emmelot P. Electron microscopic observations on negatively stained plasma membranes isolated from rat liver. J Cell Biol. 1965;26:299–305.
Revel JP, Karnovsky MJ. Hexagonal array of subunits in intercellular junctions of the mouse heart and liver. J Cell Biol. 1967;33:C7–12.
Payton BW, Bennett MVL, Pappas GD. Permeability and structure of junctional membranes at an electrotonic Synapse. Science. 1969;166:1641–3.
Goodenough DA, Stoeckenius W. The isolation of mouse hepatocyte gap junctions. J Cell Biol. 1972;54:646–56.
Goodenough DA. Bulk isolation of mouse hepatocyte gap junctions. Characterization of the principal protein, connexin. J Cell Biol. 1974;61:557–63.
Unwin PN, Zampighi G. Structure of the junction between communicating cells. Nature. 1980;283:545–9.
Unwin PN, Ennis PD. Two configurations of a channel-forming membrane protein. Nature. 1984;307:609–13.
Yeager M, Gilula NB. Membrane topology and quaternary structure of cardiac gap junction ion channels. J Mol Biol. 1992;223:929–48.
Zhang JT, Nicholson BJ. The topological structure of connexin 26 and its distribution compared to connexin 32 in hepatic gap junctions. J Membr Biol. 1994;139:15–29.
Unger VM, Kumar NM, Gilula NB, Yeager M. Projection structure of a gap junctions membrane channel at 7A resolution. Nature. 1997;4:686–9.
Unger VM, Kumar NM, Gilula NB, Yeager M. Three-dimensional structure of a recombinant gap junction membrane channel. Science. 1999;283:1176–80.
Yeager M, Unger VM, Falk MM. Synthesis, assembly and structure of gap junction intercellular channels. Curr Opin Struct Biol. 1998;8:517–24.
Oshima A, Tani K, Hiroaki Y, Fujiyoshi Y, Sosinsky GE. Three-dimensional structure of a human connexin26 gap junction channel reveals a plug in the vestibule. Proc Natl Acad Sci U S A. 2007;104:10034–9.
Maeda S, Nakagawa S, Suga M, Yamashita E, Oshima A, Fujiyoshi Y, Tsukihara T. Structure of the connexin 26 gap junction channel at 3.5 A resolution. Nature. 2009;458:597–602.
Müller DJ, Hand GM, Engel A, Sosinsky GE. Conformational changes in surface structures of isolated connexin 26 gap junctions. EMBO J. 2002;21:3598–607.
Fleishman SJ, Unger VM, Yeager M, Ben-Tal N. A Cα model for the transmembrane α helices of gap junction intercellular channels. Mol Cell. 2004;15:879–88.
Kwon T, Harris AL, Rossi A, Bargiello TA. Molecular dynamics simulations of the Cx26 hemichannel: Evaluation of structural models with Brownian dynamics. J Gen Physiol. 2011;138:475–93.
Purnick PE, Benjamin DC, Verselis VK, Bargiello T, Dowd TL. Structure of the amino terminus of a gap junction protein. Arch Biochem Biophys. 2000;381:181–90.
Zhou XW, Pfahnl A, Werner R, Hudder A, Llanes A, Luebke A, Dahl G. Identification of a pore lining segment in gap junction hemichannels. Biophys J. 1997;72:1946–53.
Kronengold J, Trexler EB, Bukauskas FF, Bargiello TA, Verselis VK. Single-channel SCAM identifies pore-lining residues in the first extracellular loop and first transmembrane domains of Cx46 hemichannels. J Gen Physiol. 2003;122:389–405.
Skerrett IM, Aronowitz J, Shin JH, Cymes G, Kasperek E, Cao FL, Nicholson BJ. Identification of amino acid residues lining the pore of a gap junction channel. J Cell Biol. 2002;159:349–59.
Kovacs JA, Baker KA, Altenberg GA, Abagyan R, Yeager M. Molecular modeling and mutagenesis of gap junction channels. Prog Biophys Mol Biol. 2007;94:15–28.
Skerrett IM, Di WL, Kasperek EM, Kelsell DP, Nicholson BJ. Aberrant gating, but a normal expression pattern, underlies the recessive phenotype of the deafness mutant Connexin26M34T. FASEB J. 2004;18:860–2.
Skerrett M, Kasperek E, Cao FL, Shin JH, Aronowitz J, Ahmed S, Nicholson BJ. Application of SCAM (substituted cysteine accessibility method) to gap junction intercellular channels. Cell Commun Adhes. 2001;8:179–85.
Pantano S, Zonta F, Mammano F. A fully atomistic model of the Cx32 connexon. PLoS One. 2008;3:1–11.
Raval A, Piana S, Eastwood MP, Dror RO, Shaw DE. Refinement of protein structure homology models via long, all-atom molecular dynamics simulations. Proteins Struct Funct Bioinforma. 2012;80:2071–9.
Nakagawa S, Maeda S, Tsukihara T. Structural and functional studies of gap junction channels. Curr Opin Struct Biol. 2010;20:423–30.
Abascal F, Zardoya R. Evolutionary analyses of gap junction protein families. Biochim Biophys Acta. 2013;1828:4–14.
Maeda S, Tsukihara T. Structure of the gap junction channel and its implications for its biological functions. Cell Mol Life Sci. 2011;68:1115–29.
Xu J, Nicholson BJ. The role of connexins in ear and skin physiology - Functional insights from disease-associated mutations. Biochim Biophys Acta Biomembr. 2013;1828:167–78.
Sosinsky G. Mixing of connexins in gap junction membrane channels. Proc Natl Acad Sci U S A. 1995;92:9210–4.
Martínez AD, Maripillán J, Acuña R, Minogue PJ, Berthoud VM, Beyer EC. Different domains are critical for oligomerization compatibility of different connexins. Biochem J. 2011;436:35–43.
Chothia C. Principles that determine the structure of proteins. Annu Rev Biochem. 1984;53:537–72.
Lesk AM, Levitt M, Chothia C. Alignment of the amino acid sequences of distantly related proteins using variable gap penalties. Protein Eng. 1986;1:77–8.
Sali A, Blundell TL. Definition of general topological equivalence in protein structures. J Mol Biol. 1990;212:403–28.
Sali A, Blundell TL. Comparative Protein Modelling by Satisfaction of Spatial Restraints. J Mol Biol. 1993;234:779–815.
Olivella M, Gonzalez A, Pardo L, Deupi X. Relation between sequence and structure in membrane proteins. Bioinformatics. 2013;29:1589–92.
Sánchez R, Šali A. Advances in comparative protein-structure modelling. Curr Opin Struct Biol. 1997;7:206–14.
Gong XQ, Nakagawa S, Tsukihara T, Bai D. A mechanism of gap junction docking revealed by functional rescue of a human-disease-linked connexin mutant. J Cell Sci. 2013;126(Pt 14):3113–20.
Emsley P, Cowtan K. Coot: Model-building tools for molecular graphics. Acta Crystallogr Sect D: Biol Crystallogr. 2004;60(12):2126–32.
Laskowski RA, Moss DS, Thornton JM. Main-chain bond lengths and bond angles in protein structures. J Mol Biol. 1993;231:1049–67.
Tong X, Aoyama H, Tsukihara T, Bai D. Charge at the 46th residue of connexin 50 is crucial for the gap-junctional unitary conductance and transjunctional voltage-dependent gating. J Physiol. 2014;592:5187–202.
Kwon T, Roux B, Jo S, Klauda JB, Harris AL, Bargiello TA. Molecular dynamics simulations of the Cx26 hemichannel: Insights into voltage-dependent loop-gating. Biophys J. 2012;102:1341–51.
Verselis VK, Trelles MP, Rubinos C, Bargiello TA, Srinivas M. Loop gating of connexin hemichannels involves movement of pore-lining residues in the first extracellular loop domain. J Biol Chem. 2009;284:4484–93.
Tang Q, Dowd TL, Verselis VK, Bargiello TA. Conformational changes in a pore-forming region underlie voltage-dependent “loop gating” of an unapposed connexin hemichannel. J Gen Physiol. 2009;133:555–70.
Vásquez M. Modeling side-chain conformation. Curr Opin Struct Biol. 1996;6:217–21.
Fiser A, Do RK, Sali A. Modeling of loops in protein structures. Protein Sci. 2000;9:1753–73.
van Gunsteren WF, Berendsen HJC. Computer Simulation of Molecular Dynamics: Methodology, Applications, and Perspectives in Chemistry. Angew Chem Int Ed Engl. 1990;29:992–1023.
Lee MR, Tsai J, Baker D, Kollman PA. Molecular dynamics in the endgame of protein structure prediction. J Mol Biol. 2001;313:417–30.
Shaw DE, Bowers KJ, Chow E, Eastwood MP, Ierardi DJ, Klepeis JL, Kuskin JS, Larson RH, Lindorff-Larsen K, Maragakis P, Moraes MA, Dror RO, Piana S, Shan Y, Towles B, Salmon JK, Grossman JP, Mackenzie KM, Bank JA, Young C, Deneroff MM, Batson B. Millisecond-scale molecular dynamics simulations on Anton. SC 09 Proc Conf High Perform Comput Netw Storage Anal 65. 2009:1. http://dl.acm.org/citation.cfm?id=1654099.
Mori T, Miyashita N, Im W, Feig M, Sugita Y. Molecular dynamics simulations of biological membranes and membrane proteins using enhanced conformational sampling algorithms. Biochim Biophys Acta Biomembr. 2016;1858(7 Pt B):1635–51.
Oostenbrink C, Villa A, Mark AE, Van Gunsteren WF. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem. 2004;25:1656–76.
Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J Chem Phys. 1983;79:926.
MacKerell ADJ, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FT, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher WE, Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiórkiewicz-Kuczera J, Yin D, Karplus M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B. 1998;102:3586–616.
Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc. 1995;117:5179–97.
Guvench O, MacKerell AD. Comparison of protein force fields for molecular dynamics simulations. Methods Mol Biol. 2008;443:63–88.
Grossfield A, Feller SE, Pitman MC. Convergence of molecular dynamics simulations of membrane proteins. Proteins. 2007;67:31–40.
Duan Y, Kollman PA. Pathways to a protein folding intermediate observed in a 1-microsecond simulation in aqueous solution. Science. 1998;282(October):740–4.
Sotomayor M, Schulten K. Single-molecule experiments in vitro and in silico. Science. 2007;316:1144–8.
Becker OM. Protein Folding: Computational Approaches. In: Becker OM, MacKerell Jr AD, Roux B, Watanabe M. Computational Biochemistry and Biophysics. New York: CRC Press; 2001. p. 371–91.
Roux B, Allen T, Bernèche S, Im W. Theoretical and computational models of biological ion channels. Q Rev Biophys. 2004;37:15–103.
Gumbart J, Trabuco LG, Schreiner E, Villa E, Schulten K. Regulation of the protein-conducting channel by a bound ribosome. Structure. 2009;17:1453–64.
Klepeis JL, Lindorff-Larsen K, Dror RO, Shaw DE. Long-timescale molecular dynamics simulations of protein structure and function. Curr Opin Struct Biol. 2009;19:120–7.
Colombres M, Garate JA, Lagos CF, Araya-Secchi R, Norambuena P, Quiroz S, Larrondo L, Pérez-Acle T, Eyzaguirre J. An eleven amino acid residue deletion expands the substrate specificity of acetyl xylan esterase II (AXE II) from Penicillium purpurogenum. J Comput Aided Mol Des. 2008;22:19–28.
Lagos CF, Caballero J, Gonzalez-Nilo FD, David Pessoa-Mahana C, Perez-Acle T. Docking and quantitative structure-activity relationship studies for the bisphenylbenzimidazole family of non-nucleoside inhibitors of HIV-1 reverse transcriptase. Chem Biol Drug Des. 2008;72:360–9.
Lizama C, Lagos CF, Lagos-Cabré R, Cantuarias L, Rivera F, Huenchuñir P, Pérez-Acle T, Carrión F, Moreno RD. Calpain inhibitors prevent p38 MAPK activation and germ cell apoptosis after heat stress in pubertal rat testes. J Cell Physiol. 2009;221:296–305.
Artigas RA, Gonzalez A, Riquelme E, Carvajal CA, Cattani A, Martínez-Aguayo A, Kalergis AM, Pérez-Acle T, Fardella CE. A novel adrenocorticotropin receptor mutation alters its structure and function, causing familial glucocorticoid deficiency. J Clin Endocrinol Metab. 2008;93:3097–105.
Gonzalez A, Duran LS, Araya-Secchi R, Garate JA, Pessoa-Mahana CD, Lagos CF, Perez-Acle T. Computational modeling study of functional microdomains in cannabinoid receptor type 1. Bioorg Med Chem. 2008;16:4378–89.
Stange C, Matus JT, Domínguez C, Perez-Acle T, Arce-Johnson P. The N-homologue LRR domain adopts a folding which explains the TMV-Cg-induced HR-like response in sensitive tobacco plants. J Mol Graph Model. 2008;26:850–60.
Ehrenfeld N, Gonzalez A, Cañón P, Medina C, Perez-Acle T, Arce-Johnson P. Structure-function relationship between the tobamovirus TMV-Cg coat protein and the HR-like response. J Gen Virol. 2008;89:809–17.
Araya-Secchi R, Garate JA, Holmes DS, Perez-Acle T. Molecular dynamics study of the archaeal aquaporin AqpM. BMC Genomics. 2011;12 Suppl 4:S8.
Gilson MK, Given JA, Bush BL, McCammon JA. The statistical-thermodynamic basis for computation of binding affinities: a critical review. Biophys J. 1997;72:1047–69.
Christen M, Van Gunsteren WF. On searching in, sampling of, and dynamically moving through conformational space of biomolecular systems: A review. J Comput Chem. 2008;29:157–66.
Allen MP, Tildesley DJ. Computer simulation of liquids. Liq Oxford Univ Press New York. 1987;18:385.
Graham TGW, Best RB. Force-induced change in protein unfolding mechanism: Discrete or continuous switch? J Phys Chem B. 2011;115:1546–61.
Balbo J, Mereghetti P, Herten DP, Wade RC. The shape of protein crowders is a major determinant of protein diffusion. Biophys J. 2013;104:1576–84.
Turchenkov DA, Bystrov VS. Conductance simulation of the purinergic P2X2, P2X4, and P2X7 ionic channels using a combined brownian dynamics and molecular dynamics approach. J Phys Chem B. 2014;118:9119–27.
Im W, Seefeld S, Roux B. A Grand Canonical Monte Carlo–Brownian Dynamics Algorithm for Simulating Ion Channels. Biophys J. 2000;79:788–801.
Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. J Chem Phys. 1953;21:1087–92.
Hastings WK. Monte carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57:97–109.
Zonta F, Polles G, Zanotti G, Mammano F. Permeation pathway of homomeric Connexin 26 and Connexin 30 channels investigated by molecular dynamics. J Biomol Struct Dyn. 2012;29:985–98.
Roux B. Ion conduction and selectivity in K(+) channels. Annu Rev Biophys Biomol Struct. 2005;34:153–71.
Zonta F, Buratto D, Cassini C, Bortolozzi M, Mammano F. Molecular dynamics simulations highlight structural and functional alterations in deafness-related M34T mutation of connexin 26. Front Physiol. 2014;5(March):85.
Zonta F, Polles G, Sanasi MF, Bortolozzi M, Mammano F. The 3.5 ångström X-ray structure of the human connexin26 gap junction channel is unlikely that of a fully open channel. Cell Commun Signal. 2013;11:15.
Luo Y, Rossi AR, Harris AL. Computational studies of molecular permeation through Connexin26 channels. Biophys J. 2016;110:584–99.
Karplus M. Molecular dynamics simulations of biomolecules. Acc Chem Res. 2002;35:321–3.
Suchyna TM, Nitsche JM, Chilton M, Harris AL, Veenstra RD, Nicholson BJ. Different ionic selectivities for connexins 26 and 32 produce rectifying gap junction channels. Biophys J. 1999;77:2968–87.
Locke D, Bian S, Li H, Harris AL. Post-translational modifications of connexin26 revealed by mass spectrometry. Biochem J. 2009;424:385–98.
Ebihara L, Steiner E. Properties of a nonjunctional current expressed from a rat connexin46 cDNA in Xenopus oocytes. J Gen Physiol. 1993;102:59–74.
Pfahnl A, Dahl G. Gating of cx46 gap junction hemichannels by calcium and voltage. Pflugers Arch Eur J Physiol. 1999;437:345–53.
Gómez-Hernández JM, de Miguel M, Larrosa B, González D, Barrio LC. Molecular basis of calcium regulation in connexin-32 hemichannels. Proc Natl Acad Sci U S A. 2003;100:16030–5.
Zonta F, Mammano F, Torsello M, Fortunati N, Orian L, Polimeno A. Role of gamma carboxylated Glu47 in connexin 26 hemichannel regulation by extracellular Ca2+: insight from a local quantum chemistry study. Biochem Biophys Res Commun. 2014;445:10–5.
Zettili N. Quantum mechanics: concepts and applications. 2009.
Bennett BC, Purdy MD, Baker KA, Acharya C, McIntire WE, Stevens RC, Zhang Q, Harris AL, Abagyan R, Yeager M. An electrostatic mechanism for Ca(2+)-mediated regulation of gap junction channels. Nat Commun. 2016;7:8770.
Escalona Y, Garate JA, Araya-Secchi R, Huynh T, Zhou R, Perez-Acle T. Exploring the membrane potential of simple dual-membrane systems as models for gap-junction channels. Biophys J. 2016, in press.
Araya-Secchi R, Perez-Acle T, Kang S-G, Huynh T, Bernardin A, Escalona Y, Garate J-A, Martínez AD, García IE, Sáez JC, Zhou R. Characterization of a novel water pocket inside the human Cx26 Hemichannel structure. Biophys J. 2014;107:599–612.
Grossfield A, Pitman MC, Feller SE, Soubias O, Gawrisch K. Internal hydration increases during activation of the G-protein-coupled receptor rhodopsin. J Mol Biol. 2008;381:478–86.
Ladbury JE. Just add water! The effect of water on the specificity of protein-ligand binding sites and its potential application to drug design. Chem Biol. 1996;3:973–80.
Gonzalez A, Perez-Acle T, Pardo L, Deupi X. Molecular basis of ligand dissociation in beta-adrenergic receptors. PLoS One. 2011;6:e23815.
Zhou R, Huang X, Margulis CJ, Berne BJ. Hydrophobic collapse in multidomain protein folding. Science. 2004;305:1605–9.
Liu P, Huang X, Zhou R, Berne BJ. Observation of a dewetting transition in the collapse of the melittin tetramer. Nature. 2005;437:159–62.
Levy Y, Onuchic JN. Water mediation in protein folding and molecular recognition. Annu Rev Biophys Biomol Struct. 2006;35:389–415.
Imai T. Roles of water in protein structure and function studied by molecular liquid theory. Front Biosci. 2009;14:1387–402.
Lomize MA, Pogozheva ID, Joo H, Mosberg HI, Lomize AL. OPM database and PPM web server: resources for positioning of proteins in membranes. Nucleic Acids Res. 2012;40:D370–6.
This work was partially supported by Fondecyt #1160574 (to TPA), PFB16 Fundación Ciencia para la Vida and ICM-Economía P09-022-F Centro Interdisciplinario de Neurociencias de Valparaíso (TPA, JAG). Conicyt Scholarship (to CP).
This article has been published as part of BMC Cell Biology Volume 18 Supplement 1, 2017: Proceedings of the International Gap Junction Conference 2015: second issue. The full contents of the supplement are available online at http://bmccellbiol.biomedcentral.com/articles/supplements/volume-18-supplement-1.
Publication charge for this article was funded by grant Fondecyt #1160574 (to TPA).
Availability of data and materials
TPA conceived the original idea. TPA and FV co-wrote and co-edited the final version of the manuscript. JAG, YE, CP and IMS co-edited the text and contributed to the discussion. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
About this article
Cite this article
Villanelo, F., Escalona, Y., Pareja-Barrueto, C. et al. Accessing gap-junction channel structure-function relationships through molecular modeling and simulations. BMC Cell Biol 18 (Suppl 1), 5 (2017). https://doi.org/10.1186/s12860-016-0121-9
- Gap-junction channels
- Structure and function
- Molecular simulation
- Homology modeling