shortcomings due to carrying out molecular dynamics (MD) simulations in the absence of an explicit membrane environment, our TM model provided important structural insights into the effects of the point mutation methionine to threonine at position 34 of Cx26 (Cx26M34T) within the TM phase of the connexon. The model prediction that this mutation leads to a constriction of the channel pore accords with experimental data obtained by various techniques [17]. A similar model of the wild type Cx26 IGJ channel, combined with use of the Renkin’s equation, predicted a value for the unitary permeability of Cx26 to cAMP that is in agreement with the experimental estimate obtained by a novel method based on FRET microscopy [18]. The application of our modeling procedure to the well-resolved structure of the potassium channel from Streptomyces Lividans (KcsA), taken as a control case (see Supporting Information Text S1), provided root mean fluctuations per residue of nearly 2 A˚ (see Supplementary Information, Figure S1), in line with the results obtained by previous theoretical studies within explicit membrane-like environment [19]. In this article, we report the construction of the first model of a full connexon within a phospholipid bilayer. Our primary aim was obviating limitations present in our original design that, in the absence of modeled extracellular and cytoplasmic loops, lacked connectivity between the TM helices. Prediction of the side chain positions allowed us to re–evaluate the relative orientation of the TM helices. Our calculations provide new molecular insights (i) on the role played by specific amino acids in keeping the correct folding of individual connexins and (ii) on the nature of the protein2protein interactions that hold together the quaternary structure of the connexon. These goals have been achieved by devising a practical protocol that combines modeling techniques with MD to produce interpretable protein models from a set of Ca coordinates corresponding to low resolution structures. Results MD drives the modeled system to a stable state By construction, the overall structural determinants of the connexon (Figure 1a) meet those of recombinant, C–terminal truncated, Cx43 channels obtained by cryo–electron microscopy [10,11]. Molecular mechanics protocols lead to conformations corresponding to local energy minima in the proximity of the starting configuration. We opted instead for MD simulations since temperature effects intrinsic in this technique avoid trapping by shallow local minima, and thus permit the entire complex to explore better the energy landscape. In the course of the MD simulations performed on the Cx32 model, the TM region featured RMSD deviations from the initial conformation ranging from 0.1 up to 0.2 nm for individual connexins, indicative of a very stable behavior. RMSD values for the whole connexin stabilized around 0.3 nm after a simulation time of nearly 5 ns (Figure 1b). The differences between the RMSD calculated over the TM domain of the connexon and those of the separated connexins indicates that the structure of the helix bundles remains stable while global quaternary structure rearrangements take place in the course of the simulation. The observed dynamical behavior is in agreement with a general analysis of the convergence achieved by MD simulations of different membrane proteins within a time window of the order of 10 ns [20]. Pore size The reported maximum constriction point of the gap junction sits within the TM domain [10]. Consistently, the pore diameter predicted by our model has its narrowest point at the level of Tyr151 (Figure 1c), resulting in an average internal diameter, All-Atoms Connexon Model dPore, calculated between opposite OH moieties that fluctuated during the MD simulation from 1.1 to 1.7 nm. A linear fit gives a nearly horizontal line (slope 20.004) at dPore = 1.4 nm, underlin- ing the good convergence for the TM domain (Figure 1d). Considering intrinsic fluctuations, this value agrees well with the range from ,1.1 nm estimated by perfusion of polyethylene glycol probes [21] to ,1.5 nm estimated by [9] based on cryoelectron microscopy data [10]. Our results point to a highly mobile behavior of the channel, originating from quaternary structure oscillations as well as side chain mobility of the single residues. Indeed, if an average dPore size of ,1.4 nm is assumed, fluctuations are needed to account for the passage of larger fluorescent tracers such as Alexa Fluor594 [22] that have a minimum (non solvated) diameter of nearly 1.6 nm. Thus, our simulations suggest that large permeating molecules may slowly creep through gap–junction channels exploiting positional rear- rangements of the connexin subunits powered by thermal energy. In so doing, they unavoidably interact with pore residues. These interactions, together with geometrical constraints, will therefore be the key determinants of the molecular permeability to larger solutes, including key metabolites and second messengers [18]. Comparison with electron microscopy map Aimed to compare the gross geometrical determinants of our model against the raw experimental data we performed a rigid fit of the Ca atoms of the final (averaged, see Methods section) connexon structure, using the same cryoelectron microscopy data of ref. [11]. Small deviations were found respect to the original model in the positions of TM3, which slightly moved towards the pore lumen (Figure 1e). The mismatch in the pitch of helix TM2 is described in the next paragraph. Although the map is too coarse to permit meaningful quantitative comparisons, visual inspection of the best fits indicates that both models match equally well the experimental data (Figure 1e). Orientation and structure of the TM2 helix In the original assignment by Fleishman and coworkers [11], the relative rotation angles of the a helices fitted into the density map were estimated by analysis of evolutionary conservation and hydrophobicity of amino-acid residues. Thus TM2 in ref. [11] occupies the most peripheral position within each of the six connexins, facing the membrane lipid phase, and consequently it is the only helix not involved in inter2connexin interactions. However, 6 out of 13 TM2 residues reported to be involved in CMTX disease (see www.expasy.org/uniprot/P08034 and linked refs. therein) are exposed to the hydrophobic phase, suggesting that an alternative orientation is possible. Therefore, we sought to explore different helical orientations that would increase the contact area between TM2 and the rest of the four helix bundle, and reduce the inter2helical gap volume (Table 1). Using rigid docking techniques, we obtained a different orientation for TM2, rotated nearly 180 degrees respect to ref. [11], which presents a higher degree of surface complementarity. In our model only 3 residues involved in CMTX disease, namely, Leu81, Ser85 and Leu89 face the lipid phase, and thus, do not contribute to the formation of inter helical contacts. Figure 2a,b show a comparison between TM2 in our structure and that of ref. [11]. This new orientation proposed for TM2 is supported by the exhaustive analysis of disease-causing mutations within the segment from Arg75 to Val95, both in Cx32 and in its highly homologous isoform Cx26 (see Discussion). Although the low resolution map does not show any obvious kinks in the TM domain, residues Leu83 through Pro87 defined a ‘‘proline kink’’ (PK) segment in a model of isolated TM2, where PLoS ONE | www.plosone.org 2 July 2008 | Volume 3 | Issue 7 | e2614