The recently published crystal structure of the Cx26 gap junction channel provides a unique opportunity for elucidation of the structure of the conductive connexin pore and the molecular determinants of its ion permeation properties (conductance, current–voltage [I-V] relations, and charge selectivity). However, the crystal structure was incomplete, most notably lacking the coordinates of the N-terminal methionine residue, which resides within the pore, and also lacking two cytosolic domains. To allow computational studies for comparison with the known channel properties, we completed the structure. Grand canonical Monte Carlo Brownian dynamics (GCMC/BD) simulations of the completed and the published Cx26 hemichannel crystal structure indicate that the pore is too narrow to permit significant ion flux. The GCMC/BD simulations predict marked inward current rectification and almost perfect anion selectivity, both inconsistent with known channel properties. The completed structure was refined by all-atom molecular dynamics (MD) simulations (220 ns total) in an explicit solvent and POPC membrane system. These MD simulations produced an equilibrated structure with a larger minimal pore diameter, which decreased the height of the permeation barrier formed by the N terminus. GCMC/BD simulations of the MD-equilibrated structure yielded more appropriate single-channel conductance and less anion/cation selectivity. However, the simulations much more closely matched experimentally determined I-V relations when the charge effects of specific co- and posttranslational modifications of Cx26 previously identified by mass spectrometry were incorporated. We conclude that the average equilibrated structure obtained after MD simulations more closely represents the open Cx26 hemichannel structure than does the crystal structure, and that co- and posttranslational modifications of Cx26 hemichannels are likely to play an important physiological role by defining the conductance and ion selectivity of Cx26 channels. Furthermore, the simulations and data suggest that experimentally observed heterogeneity in Cx26 I-V relations can be accounted for by variation in co- and posttranslational modifications.
Proteins encoded by the 21 members of the connexin gene family oligomerize to form intercellular channels (gap junction channels) by the head-to-head union of two hemichannels that span the plasma membranes of adjacent cells. A hemichannel is comprised of six connexin subunits, each of which crosses the membrane four times, which form a large central aqueous pore estimated from a variety of experimental approaches to have a limiting diameter of 12–15 Å (Oh et al., 1997; Gong and Nicholson, 2001; Harris and Locke, 2009). The recent publication of the crystal structure of a Cx26 gap junction channel at 3.5-Å resolution (Maeda et al., 2009) confirmed the overall architecture of connexin channels derived from earlier lower resolution structural work (Unger et al., 1999; Fleishman et al., 2004; Oshima et al., 2007) and resolved the controversy concerning which domains line the channel pore (Verselis, 2009). The crystal structure demonstrates that, starting at its extracellular end, the pore is lined by residues in the N-terminal half of the first extracellular loop (E1), the C-terminal half of the first transmembrane domain (TM1), the first half of the N terminus (NT), which contains the N-terminal helix (NTH), and the C-terminal half of the second transmembrane domain (TM2). This topology is congruent with the contribution of regions of the NTH, TM1, and E1 derived from substituted cysteine accessibility method (Karlin and Akabas, 1998) studies of connexin hemichannels (Zhou et al., 1997; Kronengold et al., 2003; Oh et al., 2008; Tang et al., 2009; Verselis et al., 2009; Sánchez et al., 2010).
The crystal structure provides a unique and long-awaited starting point for quantitative investigation of connexin channel structure–function. Because the voltage sensors for connexin channels are most likely located within the channel pore (Harris et al., 1981; Bargiello and Brink, 2009) rather than in discrete voltage-sensing domains outside the channel pore, as in K+ channels, knowledge of the structure and dynamics of the open state and the pore-forming region is a critical first step to understanding the molecular mechanisms of connexin voltage-dependent gating at the atomic level, as well as determinants of ionic and molecular selectivity.
The published Cx26 crystal structure was surmised to represent the conformation of an open channel (Maeda et al., 2009). This assertion was based on a continuous unobstructed path through the center of the crystal structure and use of crystallization conditions that favor occupancy of an open state. The minimum pore diameter of the crystal structure was reported to be ∼14 Å and located at the constriction formed by the tip of the NTH, which folds into the pore lumen. However, the crystal structure did not include the coordinates of the N-terminal methionine (Met1), whose inclusion would reduce the diameter of the channel pore, nor of the side chains of residues K15, S17, and S19 contained in the N terminus (NT). In addition, residues forming most of the cytoplasmic loop (CL; residues 110–124) and C terminus (residues 218–226) domains were not defined in the crystal, owing to the flexibility of these domains in Cx26 channels (Müller et al., 2002; Maeda et al., 2009). Although these domains may not be pore lining, they may be important in structure–function relations and must be included to provide a structure that can be refined by molecular dynamics (MD) simulation. Thus, the completion of the crystal structure, its refinement by MD simulation, and validation of its functional properties are critical to the elucidation of the structure–function relations of voltage-dependent gating and permeation of connexin channels.
A powerful approach to determine whether a structure represents an open ion channel is application of the grand canonical Monte Carlo Brownian dynamics (GCMC/BD) algorithm (Im et al., 2000; Im and Roux, 2001; Noskov et al., 2004; Egwolf et al., 2010). The resulting I-V relation and charge selectivity can be compared with those experimentally determined. Correspondence between the single-channel properties derived from GCMC/BD simulations and experimental results provides a sensitive test of the validity of a structure (Egwolf et al., 2010; Lee et al., 2011; Rui et al., 2011). Recent studies of ion permeation in large pore voltage-dependent anion channels (VDAC) demonstrate that the application of GCMC/BD and MD methods provides nearly equivalent results that closely correspond to experimental values and have sufficient sensitivity to discriminate between alternative models of voltage-dependent anion channel structure (Lee et al., 2011; Rui et al., 2011). In addition, the multi-ion potential of mean forces (PMFs) determined from the simulations provide information concerning the mechanism of ion permeation by defining the positions and heights of barriers to ion flux along the channel pore and their voltage dependence. GCMC/BD has several intrinsic limitations, including a rigid channel, implicit solvent-based protein–ion interactions, and an implicit membrane bilayer (Lee et al., 2011). Although these limitations are significant in the case of narrow channels, ion flux through large pore ion channels such as those formed by connexins are accurately described by GCMC/BD simulations (Noskov et al., 2004; Egwolf et al., 2010; Lee et al., 2011). A specific advantage of GCMC/BD over all-atom MD is that it is currently the only method able to produce I-V curves under asymmetric ion concentrations. There are currently no methods to simulate ion flux through a channel with asymmetric ionic concentration under periodic boundary conditions.
Intercellular channels formed by Cx26 display linear I-V relations in single-channel recordings from voltage-clamped pairs of transfected Neuro-2a cells. The slope conductance is ∼140 pS in symmetric 140-mM CsCl recording solutions (Oh et al., 1999) and ∼135 pS in 120-mM KCl solutions (Suchyna et al., 1999). The channels display moderate cation preference, 2.6:1, determined by conductance ratios (Suchyna et al., 1999). Cx26 gap junctions have been reported to be preferentially permeable to large cationic dyes over large anionic dyes (Harris and Locke, 2009), although a preferential permeability to anionic dyes correlates with Cx26 expression in specific regions of cochlear sensory epithelium (Zhao, 2005). The molecular basis of this apparent difference in Cx26 charge selectivity is unknown, although as discussed in this paper, it may result from specific co- and posttranslational charge modifications.
The linear I-V of Cx26 intercellular channels does not indicate that the I-Vs of the component hemichannels are linear, as demonstrated by Trexler et al. (2000) for Cx46 channels, where undocked hemichannels rectify inwardly while the I-V of the intercellular channel is linear. The difference is a consequence of asymmetry in the distribution of fixed charge in the connexin hemichannel, which becomes symmetric when two hemichannels are paired. In one case, the I-V of undocked Cx26 hemichannels expressed in Xenopus laevis oocytes was reported to have a slight inward rectification and a slope conductance of ∼385 pS in 100 mM KCl (González et al., 2006), whereas another report shows a slight outward rectification and a slope conductance of ∼340 pS in 140 mM KCl (Sánchez et al., 2010).
The charge selectivity of undocked Cx26 hemichannels among atomic ions has not been determined. However, the molecular determinants and mechanisms of charge selectivity have been reported for other connexins, including Cx46 (Trexler et al., 2000) and a chimeric Cx32 hemichannel in which the first extracellular loop of Cx32 is replaced with that of Cx43 (Cx32*Cx43E1; Oh et al., 2008). These studies demonstrate that the effect of a given charged residue on single-channel conductance, direction of current rectification, and charge selectivity depends on its sign, its position within the pore, the diameter of the pore in the vicinity of the charged residue, and its relative position to other charged residues. The determinants for Cx26 hemichannels are expected to be similar.
The electrostatic surface potential of the Cx26 hemichannel crystal structure has been calculated. As reported by Maeda et al. (2009), the intracellular entrance to pore has a net positive potential ≥40 kT/e, whereas the extracellular half of the pore has a net negative potential less than or equal to −40 kT/e. The corresponding segregation of charges along the length of the pore would result in formation of a P/N junction expected to display marked inward rectification and complex charge selectivity depending on the ionic strength of the recording solutions and the positions of the fixed charges in the pore (Kienker and Lear, 1995; Oh et al., 2008).
Several charge-changing co- and posttranslational modifications of human Cx26 channels expressed in stably transfected HeLa cells have been identified by mass spectrometry (Locke et al., 2009). Acetylation of the α amine of the N-terminal methionine residue (Met1 acetylation) was shown by tandem MS/MS. Other modifications shown by MS using stringent criteria include acetylation of up to seven internal lysine residues at the intracellular entrance of the channel and γ-carboxyglutamation of three glutamate residues in the channel pore. Acetylation of lysine and Met1 residues neutralizes the positive charge of these residues, whereas γ-carboxyglutamation increases the negative charge of glutamate residues. Because they alter charge, these modifications could have a substantial effect on conductance, current rectification, and charge selectivity of Cx26 channels by altering the electrostatic potential within the pore.
N-terminal acetylation, one of the most common protein modifications in eukaryotes, is the irreversible cotranslational transfer of the acetyl group of acetyl-CoA to the α-amino group of the N-terminal residue of a nascent polypeptide chain (Arnesen et al., 2009; Arnesen, 2011). N-terminal acetylation is involved in cell cycle regulation, cancer development, entry into protein degradation pathways, and protein targeting to the endoplasmic reticulum (Starheim et al., 2009; Hwang et al., 2010; Forte et al., 2011). Modification of functional properties of proteins by N-terminal acetylation has been reported in only a few cases (Hwang et al., 2010; Arnesen, 2011). N-terminal acetylation proceeds by two principle paths (Polevoda et al., 2009; Arnesen, 2011). In one, the N-terminal Met is cleaved by one of two methionine amino peptidases, and subsequently the N-terminal amine of the exposed second residue is acetylated by either the NatA or NatD enzyme complexes. In the second, the N-terminal Met is not cleaved and is acetylated by the NatB, NatC, NatE, or NatF complexes, depending on the identity of the second amino acid (i.e., the second codon rule). The N-terminal sequence of Cx26 (Met-Asp), is predicted to be acetylated by the ribosome-associated NatB complex and to occur without removal of Met1, in agreement with the MS/MS studies of Locke et al (2009). Both the catalytic (naa20p) and ancillary subunit (naa25p) proteins of NatB are expressed in Xenopus oocytes (http://www.xenbase.org) and share 99 and 83.4% sequence identity, respectively, with their human counterparts. Berns et al. (1972) reported that the N-terminal methionine of the A2 chain of calf α-crystallin (Met-Asp, as in Cx26) is acetylated when exogenously expressed in Xenopus oocytes.
Acetylation of internal lysine residues is a reversible posttranslational process governed by the opposing actions of lysine acetyltransferases (KATs) and deacetylases (Yang, 2004; Yang and Seto, 2008; Sadoul et al., 2011). The sequence specificity of targeted lysine residues is being actively investigated (Basu et al., 2009), but prediction of acetylated lysines is less precise than that of N-terminal acetylation, owing in part to the large numbers of KATs and their localization in different cellular compartments. Proteomic studies show that the homologous proteins from bacteria to humans (the “lysine acetylome”) have a high likelihood of acetylation (Choudhary et al., 2009). Consequently, both N-terminal and internal lysine Cx26 acetylations reported in HeLa cells will very likely also occur in all eukaryotes, including Xenopus oocytes. Furthermore, depending on their accessibility to relevant enzymes that add or remove the modification, lysine acetylations can be dynamically regulated in situ (Yang and Seto, 2008). Such dynamic regulation could result in variation of channel properties within a cell and among different cell types, and provide a means for physiological regulation of intercellular communication in native tissue.
There is substantially less information concerning protein modification by γ-carboxyglutamation by the vitamin K–dependent enzyme γ-glutamyl carboxylase. The enzyme is an integral membrane protein of the endoplasmic reticulum whose targets include proteins involved in blood coagulation, bone mineralization, signal transduction, and snail conotoxins (Furie et al., 1999; Brown et al., 2005). There have been no large proteomic surveys to identify its targets, nor has its expression been examined in Xenopus oocytes.
In this paper, we complete and equilibrate the atomic structure of the Cx26 hemichannel with all-atom MD simulations in an explicit solvent and 1-palmytoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) membrane system. We perform GCMC/BD simulations to determine the I-V relations and PMFs of three structures: (1) the crystal structure of the Cx26 hemichannel, which lacks Met1 (Maeda et al., 2009); (2) the completed crystal structure that includes the coordinates of residues not resolved in the crystal structure; and (3) a structure that is closest to the average structure obtained by MD simulations (defined as the “average equilibrated structure”). The results of these analyses strongly suggest that the pore diameter of the crystal and completed structures are too narrow to account for experimentally observed I-V relations and conductance. The I-V relations of the average equilibrated structure closely match those of Cx26 hemichannels experimentally determined in this study only when several of the charge-changing protein modifications reported for Cx26 channels by Locke et al. (2009) are considered. We suggest that the average equilibrated structure more closely represents the structure of the open Cx26 hemichannel than does the crystal structure, and that protein modifications are important determinants of Cx26 permeability.
MATERIALS AND METHODS
The construction of the Cx26 hemichannel was based on the Protein Data Bank coordinates (accession no. 2ZW3) of the 3.5-Å x-ray crystal structure of the Cx26 docked channel (Maeda et al., 2009). The published structure was incomplete. The atomic coordinates of Met1, all residues in the CL, the C-terminal domain, the side chains of residues K15, S17, and S19, and all the hydrogen atoms were not defined. We completed the structure in the following manner. (a) The coordinates of Met1 were taken from the nuclear magnetic resonance (NMR) structure of an N-terminal Cx26 peptide (Purnick et al., 2000). (b) The Cα coordinates of the 15-residue CL were initially assigned by linear interpolation of the distance between residues 109 (TM2) and 125 (TM3). (c) The nine residues of the C-terminal domain were built as a linear extension of the TM4 helix. (d) Subsequently, the coordinates of undefined atoms in the CL and C-terminal sequences, and those of residues K15, S17, and S19 in the NT, were designated by their CHARMM protein parameters (MacKerell et al., 1998, 2004). (e) Disulfide bonds, as identified in the crystal structure (Maeda et al., 2009), were created between C53–C180, C60–C174, and C64–C169 in each of the six subunits. (f) Protonation of histidine residues was determined by the hydrogen bond pattern in the crystal structure and assigned as either HSD or HSE. The assignments inferred from the crystal structure were in good agreement with the pKa’s predicted with PROPKA 3.0 (http://propka.ki.ku.dk/; Li et al., 2005; Bas et al., 2008; Olsson et al., 2011). (g) The original coordinates of the atoms in the x-ray structure were fixed and the completed structure was minimized in vacuum using alternate cycles of two algorithms, steepest descent (SD) and adopted basis Newton-Raphson (ABNR), until the energy difference between steps was <10−7 kcal/mol (3,570 SD and 800,000 ABNR total steps). The channel was then integrated into a POPC lipid membrane. This membrane system was chosen because the experimentally determined thickness, 38 Å (Kucerka et al., 2008), is identical to the thickness of the TM region in the crystal structure (Maeda et al., 2009), and the POPC membrane system maintains a disordered liquid state at 310 K (Buffy et al., 2004). A tetragonal periodic boundary box, 150 × 150 × 110 Å including the protein, lipid membrane, TIP3 waters, and 150 mM KCl, was constructed with a program contained in CHARMM-GUI (Membrane Builder; Jo et al., 2007, 2008, 2009). All protein atoms in this system used the CHARMM 22 force field with CMAP corrections (MacKerell et al., 1998, 2004). The CHARMM 27r force field was used for lipid molecules (Feller and MacKerell, 2000). The system was initially equilibrated for 375 ps in CHARMM (Brooks et al., 1983, 2009) at 303.15 K, using a method developed by Woolf and Roux (1994, 1996). Force constraints were reduced, first with NVT dynamics using the Langevin temperature control, and then with NPnAT dynamics using the Nosé-Hoover thermostat (Hoover, 1985). The system was then equilibrated for 140 ns at 310 K using NPnAT dynamics in NAMD2.7b2 (Phillips et al., 2005) with no constraints. NVT is the ensemble name for constant number of particles (N), volume (V), and temperature (T). NPnAT is the ensemble name for constant number of particles (N), pressure in normal direction (Pn), cross-sectional area (A), and temperature (T). Trajectories were calculated with a 2-fs time step (except the initial 75 ps where they were calculated with a 1-fs time step) with the SHAKE algorithm at 1 atm and 310 K. Electrostatic interactions were calculated using particle mesh Ewald method (Essmann et al., 1995). After a 140-ns equilibration, four independent 20-ns production-stage MD simulations were performed.
GCMC/BD simulations were performed with a computer program provided by B. Roux (University of Chicago, Chicago, IL) using methods described previously (Im et al., 2000; Noskov et al., 2004; Roux et al., 2004; Egwolf et al., 2010; Lee et al., 2011). The ion concentrations in the upper and lower buffer compartments were maintained at 100 mM KCl to correspond to the conditions used in single-channel recordings. The position-dependent ion diffusion coefficients along the z axis of the channel pore that were used in the simulation were calculated by the procedure described by Noskov et al. (2004) using pore radii determined with HOLE (Smart et al., 1996). An initial 2-µs simulation was performed to determine the time required to obtain steady-state ion flux in the system. The trajectory of ion flux reached steady state within 45 ns. I-V relations were determined by performing 20 replicate 450-ns simulations at seven voltages: 0, ±50, ±100, and ±150 mV. The PMF was calculated from the average density of K+ and Cl− along the pore axis as described in Roux et al. (2004). The values of PMF were scaled such that the PMF of ions was zero in bulk solution. Protein modifications corresponding to acetylation and γ-carboxyglutamation were mimicked by changing the charge of the modified amino acid in the input Protein Data Bank file.
RNA synthesis and oocyte injection
Two different Cx26 clones were examined, each reflecting a naturally occurring polymorphism in Cx26 sequence. Cx26S86 sequence was cloned into the pGEM-7zf+ vector (Promega) after PCR amplification of commercially obtained human DNA. The Cx26T86 clone, which corresponds to the protein sequence of the Cx26 crystal, was provided by J. Contreras (University of Medicine and Dentistry of New Jersey, Newark, NJ). RNA was transcribed from linear DNA templates by using the Message-Machine kit with T7 RNA polymerase (Applied Biosystems/Life Technologies) according to the manufacturers’ instructions. 100 nl of ∼1 ng/nl RNA was injected into each oocyte obtained from Xenopus (Xenopus 1, Inc.) by standard procedures. After RNA injection, oocytes were cultured at 16°C in pH 7.6 media containing (in mM): 100 NaCl, 2 KCl, 1.8 CaCl2, 1 MgCl2, 5 HEPES, 0.01 EDTA, 0.1 DTT, and 1 µg/ml gentamicin sulfate (Invitrogen).
Single-channel records were acquired using pCLAMP 7.0 software, an Axopatch 200B integrating patch amplifier, and a Digidata 1200A interface (Axon Instruments). I-V relations were obtained from excised patches with voltage ramps of 1.6, 3.2, or 8 s in duration. Data were acquired at 5 kHz and filtered at 1 kHz with a four-pole Bessel filter. Experiments were performed in an RC11 recording chamber (Warner Instruments). In all experiments, the pipette and bath solution contained (in mM): 100 KCl, 1 CaCl2, 3 EGTA, and 5 HEPES, adjusted to pH 8.0 with KOH. The bath solution was connected to the ground with a 3-M KCl agar bridge.
Online supplemental material
Figs. S1–S9 and Tables S1 and S2 further describe the Cx26 structure that results from MD simulations and comparison to the crystal structure. Additional GCMC/BD simulations of non-equilibrated models of Cx26 and surveys of the effects of charge modifications of the average equilibrated structure are included. Table S1 is a list of potential electrostatic interactions among amino acids and lipids in the crystal- and production-stage MD simulations. The cutoff value for hydrogen bonds was 2.4 Å. Interactions listed occur in two or more of the six subunits. Table S2 is a list of potential van der Waals interactions among amino acids in the crystal- and production-stage MD simulations. A cutoff distance of 3.0 Å was used to identify possible interactions.
Completion of the Cx26 crystal structure
The crystal structure of the Cx26 gap junction channel at 3.5-Å resolution (Fig. 1, A and B), was proposed to be that of an open channel (Maeda et al., 2009). The dimensions of the pore of the Cx26 hemichannel derived from the crystal structure, determined by HOLE (Smart et al., 1996), are shown in Fig. 1 C. The narrowest region of the pore, ∼10 Å in diameter, is formed by the D2 residue, at the tip of the NT, located well within the pore (z = −8 Å). Additional pore diameter minima are in the vicinities of residues K41, E42, A49, and D50. The 10-Å minimal value determined by HOLE, at D2, is less than the 14 Å reported (Maeda et al., 2009), but the published pore diameter was based on minimal center-to-center distances of opposed heavy atoms and did not consider atom diameters. The contribution of the N-terminal methionine (Met1) to pore diameter was not considered, as it was not resolved in the crystal. Inclusion of this residue would further narrow the pore.
The omission of the coordinates of Met1 in the crystal structure was attributed to its disordered structure (Maeda et al., 2009), although the possibility that Met1 was removed by cleavage with methionine aminopeptidase (MAP) in the insect Sf9 cells used to obtain Cx26 channels cannot a priori be ruled out. However, excision of the N-terminal methionine of proteins with N termini Met-Asp by MAP is unlikely (Sherman et al., 1985; Helbig et al., 2010). Mass spectroscopy shows that Met1 is not cleaved in Cx26 expressed in HeLa cells (Locke et al., 2009), and protein sequencing with Edman degradation has shown that it is not cleaved in native Cx32 (Met-Asn) purified form rat liver (Nicholson et al., 1981; Zimmer et al., 1987; Hertzberg et al., 1988).
In addition to Met1, the inability to resolve the coordinates of 15 residues forming part of the CL (residues 110–124), all the nine residues forming the C-terminal domain (residues 218–226), and the side chains of residues K15, S17, and S19 contained in the NT (residues 1–22) were attributed to their inherent flexibility (Maeda et al., 2009). This is in agreement with studies of these domains with atomic force microscopy in Cx26, Cx43, and Cx40 channels (Müller et al., 2002; Liu et al., 2006; Allen et al., 2011), and NMR studies of Cx26 and Cx32 NT peptides (Purnick et al., 2000; Kalmatsky et al., 2009).
The Cx26 crystal structure was completed in the manner described in Materials and methods. The 15-residue CL segment not defined in the crystal structure has high probability of disordered structure, ranging from 0.67 to 0.78, as assessed by metaPrDos (http://prdos.hgc.jp/meta/; Ishida and Kinoshita, 2007) using seven different predictors. Consequently, this segment was added as a random structure connecting the defined helical portions of TM2 with TM3. Similarly, metaPrDos predicted that the nine-residue intracellular C-terminal domain would adopt a random structure.
The resulting completed structure is shown in Fig. 1 (D and E). As expected, inclusion of Met1 reduced the minimal pore diameter substantially (Fig. 1 E), from ∼10 Å in the published structure to <6 Å, a value less than the diameter of hydrated K+ (6.62 Å) and Cl− (6.64 Å) ions (Nightingale, 1959), and close to the van der Waals diameters of dehydrated K+ (5.5 Å) and Cl− (3.5 Å). For this reason, it seems unlikely that the crystal structure would permit conduction of hydrated ions, or molecular tracers whose permeabilities define the open state of these channels. The minimal pore diameter is substantially less than the ∼12 Å determined experimentally for Cx26 (Gong and Nicholson, 2001). Consequently, we performed all-atom MD simulations in an explicit solvent and membrane system to refine the “completed” crystal structure.
All system parameters, including total system energy, periodic box size in the z dimension, and surface tension of the lipid, equilibrated within the initial 140-ns simulation (Fig. S1, A–C). The order parameter of the lipid system (Fig. S1 D) confirmed that the POPC membrane continued to be in the “disordered liquid state” (Klauda et al., 2010). The trajectory of root mean square deviations (RMSDs) of the Cx26 hemichannel equilibrated within 140 ns (Fig. S1 E), as did the trajectories of all channel parameters examined. The relaxation time of all trajectories examined, other than total RMSD, varied from 1 to 10 ns as determined by autocorrelation analyses (not depicted). Based on these relaxation times, four independent 20-ns simulations were performed on the Cx26–POPC membrane system starting from the endpoint of the 140-ns simulation to increase the probability that the system sampled an adequate conformational space (see Caves et al., 1998; Grossfield et al., 2007). MD simulation performed with the CHARMM22 force field with CMAP correction improves the prediction of side-chain structure and dynamics (Buck et al., 2006), although it has not been validated fully. NMR studies have shown that side-chain dynamics can occur in the ps to ns time range (i.e., well within the timeframe of our 140-ns equilibration). However, comparative NMR and MD studies of tryptophan side chains of gramicidin A using modified CHARMM27 force fields with CMAP and dCMAP correction suggest that simulations much longer than 20–50 ns are required to equilibrate the side-chain dynamics of some tryptophan residues (Ingólfsson et al., 2011). Notably, however, the backbone dynamics of gramicidin appeared to be adequately sampled by shorter simulations. Similarly, a comparative analysis of side-chain dynamics of the oncoprotein MDM2 with a modified AMBER force field and NMR relaxation (Li et al., 2010) suggested that 32 of 350 dihedral angles sampled may require simulation times longer than 10 ns to equilibrate. Consequently, it is possible that simulations longer than those performed in this study (220 ns total) may be necessary to ensure equilibration of all side chains. However, it is reasonable to expect that the side chains of solvent exposed residues, such as those exposed to the pore lumen, equilibrate faster than buried residues, and that the backbone of the Cx26 structure is well equilibrated by the 140-ns equilibration-stage simulation.
The atomic coordinates of the four simulations were averaged. The structure shown in Fig. 1 G was selected as representative of the average of all structures sampled during the four 20-ns production-phase simulations, based on three criteria. (1) It has the smallest RMSD from the average structure (2.76 Å), (2) the smallest deviation from the average pore dimension (RMSD = 4.69 Å), and (3) comparable probability that any given residue lines the channel pore. We define this structure as the “average equilibrated structure.”
The gross structural features of the crystal structure do not change substantially after equilibration with MD simulations.
The dimensions of the average equilibrated structure after MD simulations are compared with those of the Cx26 hemichannel crystal in Fig. 1. Although the average equilibrated structure displays a marked loss of sixfold symmetry (Fig. 1, G and H), the overall shape of the equilibrated hemichannel, a “cutoff cone” or tsuzumi (Maeda et al., 2009), is preserved. The outer diameters of both the extracellular and intracellular region of the channel at the membrane interfaces increase, from 62 to 67 Å and from 86 to 93 Å, respectively, whereas the length of the channel increased from 89 to 96 Å after MD equilibration. The increase in channel length was caused by extension of the CL from the initial completed structure of Cx26. The increase in the outer diameter of the channel is a result of relaxation in the packing of all four TM helices of each connexin monomer and a small change in the tilt angle of all helices. This is illustrated in Fig. 2, which schematically compares the “skeleton” of a portion of the crystal and the average of all equilibrated structures. The largest change in pore diameter occurred in the region of the pore formed by the NT (Figs. 1, C, F, and I, and 2). The dominating constriction of ∼6 Å at Met1 is eliminated, with the narrowest region now being distributed over the extracellular half of the pore (z = −5 to +40 Å) at ∼12–14-Å diameter.
The overall topology and orientation of the TM helices with respect to each other was unchanged with equilibration, but packing was relaxed. The channel pore remained lined by residues contained in the N-terminal half of E1, the second half of TM1, the first half of N-terminal domain, the C-terminal part of TM2, and portions of the CL.
The relaxation of the crystal structure after equilibration results in numerous changes in the positions of amino acids located within the interfaces of the four TM domains. Tables S1 and S2 provide a list of electrostatic and van der Waals interactions for each residue in the crystal and equilibrated hemichannel structures. The structure and dynamics of the TM1/E1 and N-terminal region and their relations to functional properties of connexin channels will be described in subsequent papers.
Cx26 equilibration by MD simulations substantially increases pore diameter but does not significantly change the pore-lining residues.
The pore dimensions of the completed crystal and equilibrated channel are compared in Fig. 1 (F and I). The range of fluctuations in the pore dimension observed during the four replicate 20-ns production-stage simulations is depicted by the red lines in Fig. 1 I. Notably, the pore-forming region in the segment spanning residues 42–51 (E1 near the TM1/E1 border) has the smallest standard deviation and consequently is the most stable region of the channel. Paradoxically, this region of the pore undergoes large conformational changes in forming the voltage-dependent loop–gate permeability barrier (Tang et al., 2009; Verselis et al., 2009). The region of the channel pore formed by the NT and cytoplasmic entrance (z = −45 to −25 Å) has the greatest dynamic variation.
The largest increase in pore dimension occurs in the region of Met1 (Figs. 1, F and I, and 2). At this position, the average diameter of the channel pore increases from ∼6 to ∼12.8 Å with equilibration. This value is in close agreement with the minimum pore diameter (∼12 Å) inferred from PEG exclusion studies of Cx26 intercellular channels (Gong and Nicholson, 2001). In addition to the absence of Met1, the smaller pore diameter observed in the crystal structure is partly a consequence of the inferred presence of an intersubunit hydrogen bond between the side chain of D2 and the primary amine of the backbone of T5 in adjacent subunits (Maeda et al., 2009). However, this hydrogen bond was lost within 5 ns of MD simulation and did not reform for the duration of the simulations. An intersubunit van der Waals interaction between residue W3 and M34 was also proposed to stabilize the position and structure of the NT. The MD simulations result in substantial changes in the average position of the NT and identify additional interactions between residues located in the NT and other domains. The W3–M34 interaction is not preserved in many of these alternate conformations. The new interactions result in the loss of a cavity observed in the crystal structure (Fig. 4 of Maeda et al., 2009) that separates the helical segment of the NT (residues 2–9) from the “body” of the channel formed by TM1 and TM2. Changes in interactions of the NT with residues in TM1 (including W3) and TM2 have the net average result of repositioning the NT away from the cytoplasmic entrance of the channel (Fig. 2); however, the position of Met1 varies dynamically and can also move toward the cytoplasmic entrance of the pore.
Comparisons of the average water accessibility of side chains and pore-lining probability of the crystal and equilibrated structures are presented and discussed in terms of available experimental results in Figs. S2 and S3. There is good overall agreement between experimental datasets for all substituted cysteine residues that have been reported to be modified by thiol reagents and their water accessibility and pore-lining probability in structural models. The only exception is the reported modification of V38C in Cx32*43E1 (Tang et al., 2009) and its homologue A39C in Cx46 (Kronengold et al., 2003) by MTS reagents in hemichannels expressed in Xenopus oocytes, which neither the crystal structure nor our studies show as pore lining.
Additional features of the equilibrated structure and comparisons to the crystal structure are presented in the online supplemental material. These include membrane–protein interactions, positions of extracellular domains, and comparisons of the flexibility of channel domains in the crystal and dynamic structures.
Single-channel recordings of human Cx26 hemichannels in Xenopus oocytes
We determined the I-V relations of the human Cx26 undocked hemichannels with single-channel recordings in symmetric 100 mM KCl in the outside-out patch configuration. Cx26 undocked hemichannels displayed a range of near linear I-V relations, ranging from slight inward to slight outward rectification (Fig. 3). The average slope conductance ± SEM of these channels at 0 mV are (in pS): inward/sigmoidal (Fig. 3 A), 218 ± 5; slight inward (Fig. 3 B), 270 ± 10; linear (Fig. 3 C), 240 ± 8; slight outward (Fig. 3 D), 180 ± 1; and more outward (Fig. 3 E), 205 ± 5. The I-V relations shown are representative of 17 single-channel records obtained from four batches of oocytes. Similar results were obtained for 18 cell-attached patches (not depicted). The records shown were obtained from a Cx26 clone containing a naturally occurring polymorphism at the 86th locus, S86. Similar records were obtained from a wild-type Cx26 clone, T86, which corresponds to the primary sequence of the crystallized protein. The single-channel record depicted in Fig. 3 E, which displays some outward rectification, is comparable to that reported for Cx26S86 (Sánchez et al., 2010) for inside-out patches recorded in symmetric 140 mM KCl. In contrast, González et al. (2006) reported that the I-V relation of Cx26 hemichannels displays slight inward rectification in cell-attached oocyte patches, similar to Fig. 3 B, but the reported slope conductance (∼385 pS in 100 mM KCl) is larger than the 260-pS conductance that we observe in any records. These results confirm that Cx26 hemichannels, even expressed in the same cell type, can show variation in I-V relations. The source of this heterogeneity is unknown.
Assessment of structural models with BD simulations
A central objective of this study is to identify a structural model that represents a functional, biologically relevant open state of the Cx26 hemichannel. We address this goal by comparing the single-channel I-V relations described above with those computed with the GCMC/BD algorithm developed by Roux and colleagues (Im et al., 2000; Im and Roux, 2001). The position-dependent ion diffusion coefficients used in GCMC/BD simulations are calculated from physical–chemical principles as described by Noskov et al. (2004), using the pore radius profile of the channel and a hydrodynamic approximation. The diffusion coefficients are not adjusted to fit the conductance of the experimental data. GCMC/BD simulations of wide pore ion channels provide an excellent computational description of the complex relation among parameters, such as the position and sign of fixed charges in the permeation pathway, pore diameter, and the mobility and concentrations of permeating ions. Collectively, these parameters determine the conductance and degree of current rectification. The procedure is not curve fitting; all parameters are derived from electrochemical principles calculated in terms of the atomic coordinates of the modeled channel and the dielectric properties of the implicit membrane. The resulting calculated ion fluxes are directly compared with experimental data.
Oh et al. (2008) reported that experimentally determined I-V relations, conductance, and charge selectivity were very sensitive to the position and charge of amino acid substitutions along the z axis of the Cx32*43E1 hemichannel pore. For example, the placement of a negative charge at either the second or fifth positions could be readily distinguished from the placement of a negative charge at the eighth position by differences in conductance, charge selectivity, and single-channel current rectification. The close relation among these three parameters is expected because all are determined by the relative positions of fixed charges in the channel pore, which determine the distribution and concentration of mobile ions. Thus, comparisons of experimental and computationally derived I-V provide a sensitive test of the validity of different structural models of Cx26 channels and are informative as to the molecular determinants of the ion permeation properties of Cx26 channels.
The I-V relations and associated predictions of current rectification, conductance, and charge selectivity calculated with GCMC/BD simulations are best understood in terms of the determination of the PMF of cations and anions. In simple terms, PMF describes the number of ions and their distribution along the permeation pathway as determined by fixed charge in the channel, and how ion distribution within the channel changes with the application of voltage. A region of the channel having a lower concentration of ions than the bulk phase reflects an energy barrier of defined height (kcal/mol), whereas a region with a greater concentration of ions than the bulk phase reflects an energy well. A channel will be anion selective if the summation of the PMF of anions over the length of the channel pore (i.e., the total energy of anions in the pore) is less than the corresponding PMF of cations. The degree of current rectification in the I-V relation is reflected by the voltage dependence of PMF, i.e., how the strength of the electrostatic interaction between a mobile ion and a fixed charge changes with voltage. For example, inward rectification occurs when PMF decreases with application of negative voltage; i.e., the current calculated at negative potentials will be larger than that calculated at positive potentials. Finally, the PMFs calculated in the multiple ion case include the effects of screening of fixed charge by mobile counterions. Thus, for example, the height of an energy barrier for anion in a region of negative fixed charge will be decreased by the presence of cations attracted to this region (charge screening). The Cx26 GCMC/BD simulation system is shown in Fig. 4.
The I-V relations and slope conductances of all Cx26 non-equilibrated structures do not match experimental data.
The permeability of the completed but non-equilibrated crystal structure to K+ and Cl− was calculated for the multiple ion case with GCMC/BD simulations. The low conductance and marked inward rectification of the I-V relation do not correspond to those of any open Cx26 channel determined experimentally (compare Fig. S5 A with Fig. 3). The simulations show this structure to be almost perfectly anion selective (PK/PCl = 0.010 at −100 mV), also in contrast to experimental findings, which show a moderate preference for cations (PK/PCl = 2.6) in intercellular Cx26 channels (Suchyna et al., 1999; see also Harris and Locke, 2009). As described in Fig. S5, a large cation barrier formed by the positive charge associated with Met1 appears to play a major role in determining the I-V relation and charge selectivity of this channel.
The contribution of Met1 and its associated positive charge to the barrier was further examined by performing GCMC/BD simulations of channels in which: (a) Met1 was absent (ΔMet1), as in the published crystal structure; (b) the positive charge of the N-terminal amine of Met1 was neutralized; and (c) Met1 was removed and the positive charge of the N-terminal amine of the second residue (D2) was neutralized. The first case is the closest representation of the published crystal structure, in that Met1 is absent and the positive charge of the primary amine of the D2 residue was not neutralized, but considers all other atoms not defined in the crystal structure. The second case explores the effect of the reduced pore volume resulting from Met1 presence, without the effect of its positive primary amine (as if it was acetylated). The third case removes Met1 and neutralizes the positive charge of the N-terminal amine of D2. This mimics the case where Met1 is cleaved by MAP, followed by acetylation of the penultimate D2 residue. We consider this case for completeness to explore the bases of the deviations of the crystal structure from experimental results. We reiterate that the sequence Met-Asp is not a target for efficient cleavage by MAP in eukaryotes (Sherman et al., 1985; Helbig et al., 2010) nor for acetylation by NatA acetyltransferase (Polevoda et al., 2009).
None of the simulated I-V relations for these structural variants resemble those observed experimentally. Details are provided in Fig. S5. Collectively, the results support the view that the position of Met1 and its associated positive charge are major determinants of ion permeation in the completed crystal structure. The narrow pore diameter of the channel at this position accentuates the electrostatic contribution of the positively charged Met1 residue, resulting in almost perfect anion selectivity, whereas neutralization of the positive charge of Met1 eliminates permeability to Cl− and increases permeability to K+, by unmasking the neighboring negative charge, D2. Significantly, the increase in pore diameter, by removal of Met1, which corresponds to the published crystal structure, does not resolve the deviations of the simulated currents from those observed experimentally, even when the positive charge of the N-terminal amine of D2 is neutralized. We conclude that the crystal structure is unlikely to represent the structure of an open Cx26 hemichannel.
The I-V relation and conductance of the average equilibrated structure more closely approximate, but still deviate from, experimental data.
The GCMC/BD-derived I-V relation of the average equilibrated structure is presented in Fig. 5. The properties of this channel more closely approach those of the experimental data but deviate in the degree of inward rectification, conductance, and predicted anion selectivity. The conductance of the average equilibrated channel compared with the completed crystal structure increases from 78 to 250 pS with an applied potential of −150 mV, and from 8 to 48 pS at +150 mV. These conductances reflect a large change in current rectification, from ∼1,000:1 in the completed crystal structure and ∼10:1 in ΔMet1 crystal structure to ∼5:1 in the average equilibrated channel. The near perfect anion selectivity of the completed crystal structure is substantially reduced in the average equilibrated structure: PK/PCl = 0.45 at −100 mV compared with 0.010 for the completed crystal structure.
The PMF of K+ in the multi-ion case is shown in Fig. 5 C. Notably, the K+ barrier formed in the vicinity of Met1 (Fig. S5 C) in the completed crystal structure is substantially reduced to ∼3.0 kcal/mol at 0 mV, as is the barrier formed in the vicinity of K41 (compare Fig. S5 C with Fig. 5 C). The largest cation barrier (∼3.5 kcal/mol) is now formed in the vicinity of K41. The height of this barrier is decreased at −150 mV to ∼2.2 kcal/mol and increased at 150 mV to ∼4.2 kcal/mol, reflecting the inward rectification of the current carried by K+. Overall, the changes in K+ PMF reflect the increase in pore diameter resulting from channel equilibration by the MD simulations. The pore diameter at Met1 is increased and almost the same as that in the region bounded by residues K41–D50 (Fig. 5 B).
The PMF of Cl− is presented in Fig. 5 D. Again, the barrier associated with D2 in the crystal structure (Fig. S5 D) is substantially reduced, consistent with the increase in pore diameter, whereas the barriers formed in the vicinity of E42 and A49 persist but increase slightly, from 2.2 kcal/mol in the completed crystal to 2.4 kcal/mol in the average equilibrated structure. The channel remains anion selective because the free energy of Cl− in the channel pore, as determined by the PMF of Cl−, is less than that of K+. The rectification of anion currents is primarily a result of the voltage dependence of the Cl− barrier formed in the vicinity of E42, which decreases substantially at negative potentials.
Consideration of these PMFs suggests that the deviation between the simulated and experimental I-V relations results from the incorrect distribution of charges in the channel pore.
Effects of charge changes corresponding to protein modifications
As noted previously, charge-changing co- and posttranslational modifications of human Cx26 expressed in transfected HeLa cells were identified by mass spectrometry (Locke et al., 2009). These include neutralization of the positive charge of the N-terminal amine of Met1 by N-acetylation and neutralization of the positive charge of acetylation of seven internal lysine residues at the cytoplasmic entrance, including the 15th, 102nd, 103rd, 105th, 108th, 112th, and 116th loci as well as increases of the negative charge of glutamic acid by γ-carboxyglutamation at 42nd, 47th, and 114th loci. The positions of these residues in the structure of the completed crystal structure channel are shown in Fig. 4. These posttranslational modifications would change the electrostatic potential within the channel pore and are expected to simultaneously alter the degree and direction of current rectification, charge selectivity, and channel conductance.
Charge modifications simultaneously alter the conductance, charge selectivity, and current rectification of the average equilibrated structure.
We explored the consequences of charge modifications at these positions on the I-V relations of the average equilibrated channel with GCMC/BD simulations by altering the charge of the identified amino acids in the Protein Data Bank file of the average equilibrated structure. This strategy was adopted because performing separate MD simulations of each of the 11 described charge-changing protein modifications explicitly is computationally expensive, and we had no a priori knowledge of which charge modifications, if any, would produce significant changes in I-V relations.
We explored the potential charge effects of several combinations of acetylations of N-terminal residues and internal lysines, and γ-carboxyglutamations by GCMC/BD simulations. The results are presented in Figs. S6 and S7. The PMFs of K+ and Cl− for the unmodified channel (top row in Figs. S6 and S7) indicate that the dominant K+ barriers are centrally located between z coordinates −20 to 20 Å, whereas the dominant Cl− barriers are located toward the extracellular half of the channel between z coordinates 0 to 20 Å. The general pattern resulting from neutralization of the positive charges is reduction in height of K+ barriers in the central region of the channel and a slight increase and broadening of the Cl− barriers. This effect is primarily a result of neutralization of positive charge of Met1 (N-acetylation), which is sufficient to change the charge selectivity of the channel from anion (Cl−) preferring to cation (K+) preferring (Fig. S6, column A). Neutralization of positive charge of lysine residues located close to the intracellular entrance of the channel has little effect on height of the centrally located cation barrier and anion barriers and consequently does not play a major role in determining charge selectivity. Rather, positive charges of lysine residues located at the intracellular entrance modulate the degree of charge selectivity; their neutralization by acetylation increases cation preference.
Increasing the negative charge of centrally located glutamate residues E42 and E47 by γ-carboxyglutamation has analogous effects on the PMFs of K+ and Cl− (Fig. S7). Increasing negative charge in the central region of the channel pore decreases the height and narrows the width of the K+ barrier and increases the height and width of the Cl− barriers. The net result is a decrease in the free energy of K+ in the channel pore, such that it is less than the free energy of Cl−; i.e., it becomes a cation-selective channel. All combinations of protein modifications in this central region of the channel pore that neutralize positive charge and increase negative charge result in channels predicted to have large conductance and to be almost perfectly cation selective. In the limiting case, where all residues are modified, the channel is predicted to display marked outward rectification (Fig. S7, bottom row).
These results demonstrate the sensitivity of GCMC/BD simulations to changes in the distribution of charge in the channel pore. In these simulations, each alteration of charge has a distinguishable, intuitively understandable effect that can be compared with experimental data; charge changes at different positions/regions have different effects. Notably, charge modifications at intracellular positions have much less effect on I-V relations than those more centrally located. This is a consequence of the substantially larger pore diameter at the intracellular entrance. The diameter of the channel pore is narrowest in the central region between z coordinates −20 to 20 Å. Not surprisingly, charge modifications in this region have the largest effect on the electrostatic profile of the channel and ion permeation.
Neutralization of positive charge is sufficient to replicate experimental I-V relations.
Several combinations of acetylation-mimicking charge neutralizations of Met1 and internal lysines produced simulated I-V relations that closely matched those observed experimentally (Fig. 6). These include neutralizations of Met1 and all six lysine residues located at the TM2/CL border (Fig. 6, A–D), as well as Met1 and all seven internal lysine residues (Fig. 6, E–H). The location of these charged residues is shown in Fig. 4 B. The GCMC/BD calculations reproduce both the conductance and reduced current rectification observed experimentally. Note that the number of lysine residues neutralized modulates both the slope conductance and degree of outward current rectification.
Both of the modified channel forms in Fig. 6 (A–H) are predicted to be cation selective with PK/PCl = 4.6 and 3.6 at −100 mV, respectively. These values are close to the charge selectivity PK/PCl of ∼2.6 of intercellular Cx26 channels (Suchyna et al., 1999). The basis for the increased cation selectivity of these channels is evident in the PMFs of K+ and Cl− ions shown in Fig. 6 (C, D, G, and H). In both cases, the heights of all cation barriers are substantially reduced compared with those of the unmodified average equilibrated structure. The slight outward rectification of the K+ currents is explained by the voltage dependence of K+ PMF. The peak heights of the Cl− barriers do not change substantially, but as indicated in Fig. 6 (D and E), the Cl− barriers are broadened, extending toward the cytoplasmic entrance of the channel. This is most likely a result of “unmasking” the negative charge of D2 by neutralization of the positive charge of Met1, noted previously in the non-equilibrated structure (Fig. S5, E–H). The net effect is a substantial decrease in the contribution of Cl− current and increase in the contribution of K+ current to total current.
The negative charge of D2 is also required for Cx26 cation selectivity.
We further explored the basis for the anion to cation switch in charge selectivity resulting from Met1 neutralization by simulations of average equilibrated channels in which (a) both the positive charge Met1 and negative charge of D2 were neutralized and (b) only the negative charge of D2 residue was neutralized (Fig. S8). In both cases, neutralization of the negative charge of D2 results in anion-selective channels. Thus, neutralization of the positive charge of Met1 switches charge selectivity from anion to cation only when the second locus is negatively charged (compare Fig. S8, B with D).
The reason for this is illustrated by comparing the PMFs of K+ and Cl− (Fig. S8). In the anion-selective unmodified channel, the height of the K+ barrier formed in the vicinity of K41 dominates but its contribution to the electrostatic profile of the channels is decreased when the second residue is negatively charged and Met1 is neutralized (compare Fig. S8, I with J). The width of the Cl− barrier is increased when the positive charge of Met1 is neutralized (compare Fig. S8, M with N). Thus, neutralization of Met1 accentuates the contribution of the negative charge of the second residue to the electrostatic profile of the channel. This causes the free energy of K+ to be less than that of Cl− in the channel, resulting in cation selectivity. Note that the PMFs of both K+ and Cl− are similar for the unmodified channel, the D2-neutral channel, and the D2-neutral Met1-neutral channel as shown in Fig. S8, accounting for their anion selectivity.
To summarize, the GCMC/BD simulations described in the above sections suggest that the cation selectivity of the average equilibrated structure of Cx26 requires both neutralization of Met1 and the presence of a negative charge at the second position. Neutralizations of internal lysine residues serve only to modulate the degree of Cx26 cation selectivity but cannot induce it on their own.
Increasing negative charge by γ-carboxyglutamation alone does not replicate experimental I-V relations.
As illustrated in Fig. S7, increasing the negative charge of centrally located glutamate residues alone by γ-carboxyglutamation does not result in computed I-V relations that correspond to any experimental I-V relations. However, neutralization of the six TM2/CL lysines coupled with γ-carboxyglutamation of glutamate residues results in a predicted I-V relation (Fig. 6, I–L) that corresponds closely to the inward/sigmoidal I-V relation observed experimentally (Fig. 3 A). Similarly, increasing the negative charge of the three glutamate residues combined with the neutralization of Met1, K15, and K116 produces an I-V relation with GCMC/BD simulations (Fig. 6, M–P) that resembles the experimentally observed I-V relations that display slight inward rectification (Fig. 3 B). In this case, we cannot ascertain if the experimental I-V relation would be sigmoidal at larger negative potentials, as predicted by GCMC/BD. It is difficult in experimental studies to determine the I-V relation over the full ±150-mV range because of voltage-dependent channel closure and seal instability at large voltages.
Although increasing negative charge by γ-carboxyglutamation is by itself not sufficient, it can in combination with neutralization of positive charge by acetylation produce I-V relations that resemble some I-V relations that are observed experimentally. But the near perfect cation selectivity of γ-carboxyglutamated channels contrasts experimental observations. Overall, the results predict that acetylation of the N-terminal Met and internal lysines are the key protein modifications that determine Cx26 ionic permeation.
Surprisingly, only the experimentally observed linear I-V relation (Fig. 3 C) was not matched by any of the predicted I-V relations for the tested combinations of acetylations and γ-carboxyglutamations (Figs. S6 and S7). However, it is likely that some combinations that neutralize fewer than all six lysine residues in TM2/CL would produce a linear I-V relation, as these would fall between the slight inwardly rectifying channel (Met/K15) and the slightly outwardly rectifying channel (Met/CL) shown in Fig. S6. Another possibility is the formation of “heteromeric” channels formed by combinations of unmodified and modified subunits or combinations of differentially modified subunits. The possibility of such channels arises from the likelihood that not all sites of acetylation or γ-carboxyglutamation are so modified at a given time. Differences in pKa of lysine residues, leading to differences in protonation, may also play a role, although determination of pKa of these residues with PROPKA 3.0 suggests that all will be deprotonated at neutral pH (all predicted pKa’s were in the range of 10).
Protein modifications of completed and incomplete, but non-equilibrated, crystal structures do not result in I-V relations that match experimental results.
In contrast to the results based on the average equilibrated structure, the I-V relations produced by GCMC/BD simulations of the crystal structure and the completed but non-equilibrated crystal structure that incorporated all combinations of charged protein modifications described in Figs. S6 and S7 deviated strikingly from all experimentally observed I-V relations. This is illustrated in Fig. S9 for the limiting cases of the reported modifications: (a) neutralization of the N-terminal amine only (considered previously; Fig. S5); (b) neutralization of the N-terminal amine and all internal lysines; and (c) neutralization of N-terminal amine, all internal lysines, and increasing the negative charge of the three glutamate residues by carboxylation. The deviations from the experimentally derived I-V relations are primarily the result of a dominating barrier to either or both K+ and Cl− that are formed in the vicinity of the constriction at the tip of the NT. We conclude that neither the completed nor ΔMet1 crystal structures represent the structure of the biological open Cx26 hemichannel, in that no experimentally determined protein modifications can resolve the deviations in the computationally derived I-V relations and those observed experimentally.
The GCMC/BD simulations strongly suggest that the average equilibrated structure resulting from MD simulations closely approximates the structure of the biological open Cx26 hemichannel. Although the MD simulations do not change the overall structural features evident in the crystal, they refine the crystal structure by relaxation of the packing of all four TM helices of each connexin monomer and a small change in the tilt angle of all helices. The net effect is a substantial increase in the channel pore diameter, primarily in the region formed by the NT, and the repositioning of many amino acids located at the interfaces of the four TM domains, resulting in more stable networks of electrostatic and van der Waals interactions (unpublished data).
The I-V relations of the average MD-equilibrated structure determined by GCMC/BD closely match experimentally observed I-V relations when the charge effects of specific combinations of co- and posttranslational modifications identified by mass spectroscopy of Cx26 channels exogenously expressed in mammalian cells (Locke et al., 2009) were incorporated into the model. In contrast, the pore dimensions of the crystal and completed crystal structures, without MD equilibration, are too narrow to allow significant ion flux. In these structures, the tip of the NT forms a constriction that results in a large barrier to ion flux. GCMC/BD simulations showed that none of the modifications of the NT made to the crystal structure or a completed crystal structure that includes Met1 could account for any of the experimentally observed I-V relations.
A potential limitation of this study is that the structures of the charge modifications were not explicitly defined in GCMC/BD simulations. Instead, the charge effect of a given modification was mimicked by altering the charge of the relevant group in the Protein Data Bank file, and the resulting structures were not further equilibrated by MD. Also, simple neutralization does not exactly reproduce the effect of acetylation of a positively charged residue, in that it does not result in zero charge distribution of the modified residue; acetylation will produce a residue with some dipole character. However, the calculated dipole moment of acetylated methionine and lysine does not differ substantially from that of methionine and lysine neutralized by modification of the parameter file, and thus mimicking the charge modification in this manner may be adequate to represent, to a first order, the behavior of the system in GCMC/BD simulations. We used this simple charge-neutralization strategy because of the high computational cost of performing MD simulations of all 11 charge-modifying protein modifications individually. It is possible that the structure of the acetylated residue affects protein conformation. However, six of seven modified lysines and one of three sites of glutamate carboxylation occur at sites located in the highly flexible TM2/CL region and do not form charge interactions with residues outside this region (Table S1). They are therefore unlikely to have a direct effect on the structure or dynamic motions of other Cx26 channel domains. Of particular interest are structural changes that may be caused by N-terminal acetylation and by carboxylation of E42 and E47. We are currently performing additional MD simulations to determine if the Cx26 structure changes substantially when these particular modifications are included explicitly, and if they affect permeability determined by GCMC/BD simulation.
A principal prediction of the GCMC/BD simulations is that neutralization of positive charges by acetylation of the N-terminal amine of Met1 and seven lysine residues located at the cytoplasmic entrance to the pore are critical determinants of the I-V relations and charge selectivity of Cx26 hemichannels. The computationally derived I-V relations of channels with charge neutralizations at these positions closely matched the single-channel I-V relations of Cx26 channels expressed in Xenopus oocytes reported both in this study and in that of Sánchez et al. (2010). Of particular note is the predicted role of N-terminal acetylation in determining the charge selectivity of Cx26 hemichannels. All channels formed by an unmodified (positively charged) N-terminal residue are predicted to be anion selective, whereas all channels with an acetylated (uncharged) N-terminal residue are predicted to be cation selective. The switch from anion to cation selectivity arising from Met1 neutralization is a consequence of the accentuation of the contribution of the negative charge of D2 in determining the electrostatic profile of the channel. The acetylated (uncharged) state of the seven internal lysine residues is predicted to modulate the degree of cation charge selectivity, single-channel conductance, and the amount and direction of current rectification. Notably, the charge of the second residue was also shown to determine the charge selectivity, conductance, and current rectification of the closely related Cx32*43E1 hemichannel (Oh et al., 2008). There, the substitution of a negative charge at the second locus, N2E, more than doubled single-channel conductance, increased cation selectivity severalfold, and reversed the direction of current rectification from inward to slightly outward. The positive charge associated with Met1 in the N2E Cx32*43E1 channel is expected to be neutralized by N-terminal acetylation as its sequence (Met-Glu) is also an established target for the NatB acetylating enzyme complex (Starheim et al., 2009). Consequently, the marked effects of N2E in determining ion permeation of the Cx32*43E1 hemichannel appear to result from both the addition of negative charge at the second locus and neutralization of positive charge by acetylation of Met1.
The omission of the N-terminal Met in the crystal structure was attributed to disordered structure reflecting the flexibility of the NT rather than cleavage by MAP in insect Sf9 cells that were the source of the Cx26 used in the crystallization. Although the absence of Met1 in the crystal structure should be examined further, the view that it is present but disordered is supported by proteomic studies that define the cleavage rules for N termini by MAPs and established the evolutionary conservation of protein target sequences (Frottin et al., 2006; Helbig et al., 2010). Met1 in proteins containing Asp or Asn at the second position, like Cx26 and other members of Group 1 connexins (Beyer and Berthoud, 2009), are rarely excised in eukaryotes including Drosophila melanogaster (Sherman et al., 1985; Bradshaw et al., 1998; Frottin et al., 2006; Goetze et al., 2009; Drag et al., 2010; Helbig et al., 2010). The MS/MS studies of Locke et al. (2009) demonstrated conclusively that the N-terminal Met residue is not cleaved in Cx26 exogenously expressed in transfected HeLa cells, and that the N-terminal amine is acetylated. The N-terminal Met1 is cleaved and the N-terminal amine of the second residue is acetylated in Cx46 and Cx50 in bovine lens (Shearer et al., 2008). These are Group 2 connexins, characterized by Met-Gly N termini, which are predicted to be cleaved with MAP and acetylated by the highly conserved NatA complex. Thus, all three connexin N termini examined to date by MS/MS conform to established rules governing cleavage and acetylation of N termini.
The evolutionary conservation of the mechanisms of N-terminal acetylation strongly suggests that the acetylated state of the N-terminal Met is the same in Xenopus oocytes and mammalian cells. Indeed, González et al. (2006) report that the I-V relations of Cx26 undocked hemichannels expressed in oocytes are similar, if not identical, to those expressed in transfected human Neuro-2a cells. Both display linear/slight inward current rectification. Similarly, Oh et al. (1999) report that the degree of inward rectification of the I-V relations of single Cx26/Cx32 heterotypic intercellular channels expressed in Neuro-2A cells is identical to the voltage dependence of initial currents (open-channel rectification) of Cx26/Cx32 intercellular channels expressed in Xenopus oocytes. The sensitivity of the GCMC/BD-calculated I-Vs to each acetylation suggests that the consistency seen in mammalian cells and Xenopus oocytes would only occur if the connexins had the same protein modifications in both expression systems.
The heterogeneity in the Cx26 I-V relations that we report may be related to the observation that the acetylated state of protein targets is sensitive to the cytoplasmic concentration of acetyl-CoA and hence to the metabolic status of cells (Wellen et al., 2009). Rathmell and Newgard (2009) suggested that differences in the affinity of KATs for acetyl-CoA might account for the differences in sensitivity of various protein targets to metabolic stress. In our study, oocytes were cultured in minimal salt media lacking glucose and pyruvate, conditions that may cause metabolic stress and reduced concentration of acetyl Co-A.
The potential physiological relevance of Cx26 acetylation is apparent by studies of the permeation of charged dyes through Cx26-junctional channels in mammalian cochlea to fluorescent dyes. Cx26 and the closely related Cx30 are widely expressed in supporting cells of cochlea (Zhao and Yu, 2006), yet there are marked differences in junctional charge selectivity in different regions of cochlea (Zhao, 2005). Permeability to anionic dyes only, in specific regions of the cochlea, is correlated with the presence of cells that express only Cx26, whereas cells that express both Cx26 and Cx30 in other regions of the cochlea favor permeability to cationic dyes. However, spread of cationic dyes is observed among cochlear cells in Cx30 knockout mice, which continue to express Cx26 (Chang et al., 2008). Collectively, these results strongly suggest that the charge selectivity of Cx26 in the cochlea is cell specific and appears to be spatially distinct within the cochlear sensory epithelium. Our results, which indicate that charge selectivity is likely determined by the acetylated state of the N-terminal methionine and further modulated by acetylation of internal lysines, provide an explanation for the apparent regional differences in Cx26 perm-selectivity. N-acetylated Cx26 channels are predicted to be cation selective, whereas Cx26 channels with an unmodified Met1 are predicted to be anion selective. Thus, regulated Cx26 acetylation may play a fundamental role in cochlear health, physiology, and in other tissues as well by fine-tuning Cx26 charge selectivity in a cell- and tissue-specific manner.
We thank Dr. Benoit Roux (University of Chicago); Dr. Jeff Klauda (University of Maryland); and Dr. Wonpil Im, Dr. Kyu Il Lee, and Dr. Sunhwan Jo (University of Kansas) for their assistance and advice in constructing the computational systems used in this study. We thank Dr. Jorge Contreras (University of Medicine and Dentistry of New Jersey [UMDNJ]) for providing the human Cx26 clone and Dr. Qingxiu Tang (Albert Einstein College of Medicine) for advice in single-channel recording and for providing Xenopus oocytes. We thank Dr. Bogdan Polevoda (University of Rochester Medical Center) and Dr. Thomas Arnesen (University of Bergen) for discussion of N-terminal acetylation, and the reviewers and editors for their constructive criticisms of the manuscript. The computational aspects of this work were performed on the IBM iDataPlex Cluster managed by the Division of High Performance and Research Computing at UMDNJ and is gratefully acknowledged. Additional computational resources were provided by the National Science Foundation through TeraGrid resources provided by National Center for Supercomputing Applications (grant TG-IBN090012) and by the High Performance Computing Core Service at Albert Einstein College of Medicine.
This work was supported in part by grants from the National Institutes of Health (RO1 NS056509 to A.L. Harris and RO1 GM064889 to T.A. Bargiello). Additional financial support from Albert Einstein College of Medicine to T.A. Bargiello is gratefully acknowledged.
Edward N. Pugh Jr. served as editor.
grand canonical Monte Carlo Brownian dynamics
nuclear magnetic resonance
potential of mean force
root mean square deviation