Our senses of touch and hearing are dependent on the conversion of external mechanical forces into electrical impulses by the opening of mechanosensitive channels in sensory cells. This remarkable feat involves the conversion of a macroscopic mechanical displacement into a subnanoscopic conformational change within the ion channel. The mechanosensitive channel NOMPC, responsible for hearing and touch in flies, is a homotetramer composed of four pore-forming transmembrane domains and four helical chains of 29 ankyrin repeats that extend 150 Å into the cytoplasm. Previous work has shown that the ankyrin chains behave as biological springs under extension and that tethering them to microtubules could be involved in the transmission of external forces to the NOMPC gate. Here we combine normal mode analysis (NMA), full-atom molecular dynamics simulations, and continuum mechanics to characterize the material properties of the chains under extreme compression and extension. NMA reveals that the lowest-frequency modes of motion correspond to fourfold symmetric compression/extension along the channel, and the lowest-frequency symmetric mode for the isolated channel domain involves rotations of the TRP domain, a putative gating element. Finite element modeling reveals that the ankyrin chains behave as a soft spring with a linear, effective spring constantof 22 pN/nm for deflections ≤15 Å. Force–balance analysis shows that the entire channel undergoes rigid body rotation during compression, and more importantly, each chain exerts a positive twisting moment on its respective linker helices and TRP domain. This torque is a model-independent consequence of the bundle geometry and would cause a clockwise rotation of the TRP domain when viewed from the cytoplasm. Force transmission to the channel for compressions >15 Å depends on the nature of helix–helix contact. Our work reveals that compression of the ankyrin chains imparts a rotational torque on the TRP domain, which potentially results in channel opening.
Introduction
NOMPC is a mechanosensitive channel that mediates touch and hearing sensation in Drosophila (Liang et al., 2013; Zanini and Göpfert, 2013; Yan et al., 2013; Zhang et al., 2015), and it is the founding member of the larger transient receptor potential family of ion channels that respond to mechanotransduction (TRPN). The recent NOMPC structure revealed that the channel is composed of four identical subunits, each with a transmembrane (TM) domain and a large cytoplasmic domain formed by 29 ankyrin repeats (ARs; Fig. 1 A; Jin et al., 2017). The TM domain adopts a common fold found in other transient receptor potential (TRP) channels and voltage-gated potassium KV channels, in which the four subunits come together to form a single ion pore. Each TM subunit is composed of six α helices (labeled S1–S6) with the S1–S4 segments surrounding the central pore domain (S5 and S6) and linked via the S4–S5 linker (Fig. 1 A, orange). Directly below, and in contact with the S4–S5 linker, is the TRP domain (Fig. 1 A, dark blue), which is highly conserved in TRP channels and implicated in channel gating (Venkatachalam and Montell, 2007; Cao et al., 2013). In NOMPC, the TRP domain connects to the S6 inner pore helix, and it is sandwiched by the S4–S5 linker from the top and the linker helices (Fig. 1 A, violet) from below. The linker helices are a group of short α helices that includes the preS1 elbow that leads into the S1 TM segment, and the region makes a hydrophobic connection between the AR domain and the TM domain. Each AR chain extends away from the membrane ∼150 Å, and each AR repeat is made up of two antiparallel α helices followed by a β-hairpin loop connecting one repeat to the next (Gaudet, 2008; Fig. 1 A, inset). While ARs are structurally similar, their sequence varies. Due to the geometry of the repeats, when multiple ARs are connected they form a superhelical configuration. In NOMPC, each ankyrin chain forms a left-handed helix tilted ∼20° with respect to the membrane normal (z axis). The N-terminal end of the chain is thought to contact microtubules (MTs), as suggested in Fig. 1 A, and the degree of splay of the four chains is consistent with the diameter of an MT. The bundle formed by the four chains has two contact points: contact point 1 is located toward the N terminus close to the MT, and contact point 2 is located closer to the TM domain.
A long-standing question in sensory biology is how macroscopic, mechanical forces are transduced into electrical signals by mechanosensitive channels. In 2004, Howard and Bechstedt impressively proposed that the 29 ARs in NOMPC form a bundle of helices that constitutes a gating spring that directly couples MTs to the channel domain to sense external forces and directly pull on the pore to open the channel (Howard and Bechstedt, 2004). They based their hypothesis on x-ray data showing that 12 ARs form a regular structure that makes up approximately one half of a turn of a helix (Michaely et al., 2002) and that therefore, 29 ARs should form a nearly full helical turn. There are several positive attributes to this spring hypothesis, as noted by the authors. First, a full helical turn allows forces to be transmitted to the channel without imparting torques, and they hypothesized that this is why 29 ARs are so conserved in nature. Second, they estimated the spring stiffness of an independent chain to be ∼1 pN/nm or ∼4 pN/nm for four independent chains in parallel like the NOMPC helical bundle. Such a soft spring would be a highly compliant object that would allow for large compressions on the order of 100 Å and would provide a graded response to channel opening for small deformations. Both single-molecule experiments and computational simulations have demonstrated that isolated chains formed by 24 ARs behave as biological, linear springs with a spring constant of 2.4–5 pN/nm (Sotomayor et al., 2005; Lee et al., 2006), only slightly larger than the value proposed by Howard and Bechstedt (2004). Nonetheless, the intricate coiling of the ARs in NOMPC, together with their two points of contact, suggests that the chains may not act independently as four parallel springs, at least under some deformations.
Using transmission electron microscopy, the Howard laboratory went on to reveal the ultrastructure of the campaniform sensory receptor in flies (Liang et al., 2013). The images showed that the membrane is connected to an array of parallel MTs via structures termed membrane–MT connectors (MMCs). These connectors share physical properties with the AR domain of NOMPC, and they colocalize with the AR domain. Later it was shown that, in fact, ARs associate with the MT and transmit the forces to open the channel (Zhang et al., 2015). The transmission electron microscopy also showed that the MTs sit on an electron-dense material, which the authors hypothesized to be stiff (Liang et al., 2013). From these images, the authors constructed a finite element model that explored the stress–strain relationship in the system to extract how the membrane, MMC, and MTs interact under load. Their primary result is that the MMCs are likely the most pliant elements, and they account for 50% to 80–90% of the total deformation (depending on the assumed stiffness), followed by the membrane, which accounts for the remaining strain (Liang et al., 2013). While incredibly enlightening, this model was created before the high-resolution determination of the NOMPC channel, and hence, it remains to be shown how the chains and the chain bundle respond to force and how this force eventually opens the channel at the atomic level. The three NOMPC class averaged structures all have closed pore domains, yet differences in the AR chains and linker regions hinted at motions that may be associated with channel opening (Jin et al., 2017). First, comparing class 1 to class 2 revealed a compression–extension motion in the helical bundle with class 2 more extended by ∼3 Å and the greatest displacement occurring over the N-terminal portion of the AR domain. The N-terminal ARs of class 3 are not resolved, suggesting that this region is highly mobile. Second, as the AR chains compress from class 2 to 1, each linker domain, which is adjacent to a TRP domain, undergoes a very subtle clockwise (CW) rotation when viewed from the cytoplasm. The authors hypothesize that this CW linker motion is consistent with the CW TRP domain motion that occurs between the open and closed TRPV1 ion channel structures (Cao et al., 2013; Liao et al., 2013).
Here we employ a number of computational approaches to address questions concerning the mechanical properties and gating motions of the NOMPC channel. Protein function is often probed computationally through the use of fully atomistic MD simulations; however, given the large size of the NOMPC channel—over 1 million atoms when fully solvated in a lipid bilayer—this approach would fail to provide meaningful information about how the channel responds to large deformations over the short, hundred-nanosecond timescales that we could access. As an alternative, we first employed normal mode analysis (NMA) to uncover the low-frequency, large-amplitude, collective motions of the entire protein as well as the motions of the isolated ion channel together with the TRP domain. The lowest-frequency mode recapitulates the compression–extension motion of the helical bundle along the z axis observed in class averages, while showing little to no relative motion in the TM domain. Intramolecular motions in the channel domain were only observed once we computed the modes in the absence of the chains, and in this case, the lowest-frequency mode that preserves fourfold symmetry involved pivoting of the TRP domain, as suggested earlier (Cao et al., 2013; Liao et al., 2013). Since NMA is a linear theory that describes small displacements around a starting structure, we employed a multiscale model to probe large deformations of the helical bundle. We first performed fully atomistic MD simulations on small, isolated stretches of the AR domain to determine the material properties of the chain from the natural fluctuations of its constituent components. We then used this information to construct a finite element model of the entire helical bundle starting from a geometrically faithful representation of the four coiled chains terminating into a fixed plate we call the TRP region (Figs. 1 B and S1). We developed three models based on different assumptions about how the contact points influence the mechanics: (1) no contact; (2) frictionless contact; and (3) rough contact, and we validated the no-contact numerical results analytically. When a load is applied from the N terminus along the z axis with the TRP region fixed in z but free to rotate in the plane of the membrane, all three models behave as soft, linear springs with an effective spring constant for the entire bundle of 22 pN/nm over the first 15 Å of the deformation. Compression imparts a torque on the TRP region causing the entire channel domain to rotate counterclockwise (CCW) in the membrane when viewed from the cytoplasm. A counteracting CW torque develops at the tip of each AR chain as it attaches to the TRP region, which counteracts the rigid body rotation and keeps the system in static equilibrium. This later torque imparts forces to the linker helices and TRP domain that would cause them to rotate CW as observed for the opening of the TRPV1 channel (Cao et al., 2013; Liao et al., 2013). When compressed beyond 15 Å, the rough and frictionless contact models diverge with the rough model “locking up” and increasing its stiffness nearly five times that of the frictionless contact model. In addition, the frictionless model changes the direction of the torque imparted to the end of each AR chain, which would impact the conformation of the TRP domain, potentially closing the channel. Based on our analysis, we suggest additional experiments that would differentiate these models to further elucidate how mechanical forces gate the NOMPC channel.
Materials and methods
NMA
The Hessian matrix was constructed from the Cα backbone atoms of the NOMPC structure (PDB ID 5VKQ; Jin et al., 2017) using the ProDy web server (Bakan et al., 2011, 2014) with default values of 15 Å for the cutoff and 1.0 for gamma. The eigenvectors and eigenvalues for the 10 lowest modes (ignoring the 6 trivial modes corresponding to pure translation and rotation motions) were computed in MATLAB and verified against the results reported by ProDy. For the entire channel we used residues 124–1660 from all four chains, while the Hessian for the channel-linker–only analysis was computed from residues 1140–1660, which includes the linker helices, TRP domain, and TM segments. To mimic ankyrin domain binding to the MT, we modified the Hessian matrix to immobilize the most distal Cα atoms that would contact the MT (residue IDs 145–147, 179–181, and 213–215). Mathematically, we zeroed out all elements of the Hessian matrix in the rows and columns corresponding to these atoms before computing the eigenvectors and eigenvalues. All modes were visualized using the Normal Mode Wizard plugin (Bakan et al., 2011) in Visual Molecular Dynamics (Humphrey et al., 1996).
MD simulations
All-atom MD simulations were performed on the intracellular domain of NOMPC (PDB ID 5VKQ; Jin et al., 2017) using Amber 16 with CUDA-enabled Particle Mesh Ewald MD along with the ff14SB force field (Maier et al., 2015) and TIP4P-Ew (Horn et al., 2005) water model. Specifically, we simulated three stretches of one filament, which we call ANK1 (AR 12–17; blue), ANK2 (AR 17–22; red), and ANK3 (AR 22–27; light green; Figs. 1 A and S2). We did not examine the mechanics of AR 1–11 due to low resolution in this region. Within each stretch, we refer to the ith AR as ARi. For each stretch, we modeled all the missing side chains using Phyre2 (Kelley et al., 2015). Additionally, we use Amber to generate a disulfide bond between C586 and C628 in ANK1. We solvated the systems in 150 mM KCl resulting in box sizes of ∼85 Å in length containing ∼80,000 atoms. Minimization was performed using 500 steps of steepest descent followed by 7,500 steps of conjugate gradient. The system was then initially equilibrated for 200 ps in the NVT ensemble at 298°K using a Langevin thermostat with harmonic restraints applied to all backbone atoms and the side-chain atoms of the modeled residues. Additionally, we applied restraints throughout the entire simulation to the dihedral angles of the flanking ARs (AR1 and AR6) to prevent unfolding of the ends. We further equilibrated the system for a total of 1.4 ns during which the restraints on the middle span of ARs (AR2–AR5) were gradually reduced in five steps starting with spring constants of k = 5 kcal/mol/Å2. Production runs were performed for ∼450 ns in the NPT ensemble, and snapshots were saved every 20 ps for analysis. For the equilibration and production, we used the Langevin thermostat with a friction coefficient of 1 ps−1 to maintain the temperature at 298°K and the Monte Carlo barostat with an isotropic cell to maintain the pressure at 1 bar. Hydrogen atoms were constrained with SHAKE, and a 2-fs integration time step was employed. The nonbonded interactions were cut off at 10 Å, and the electrostatic interactions were calculated using the Particle Mesh Ewald method. All analysis was performed on the internal ARs 2–5 ignoring the restrained flanking repeats. We then analyzed the trajectories using the last 350 ns of each simulation, and the errors represent the 95% confidence interval (CI) obtained using a Monte Carlo bootstrap technique (Rubin, 1981).
The stiffness and torsional spring constants were calculated using the methods outlined in (Matsushita et al., 2010). Briefly, the variances in each relevant degree of freedom were computed directly from the simulation data followed by the use of the equipartition theorem to assign an effective spring constant to each mode of deformation. We calculated the stiffness spring constant by considering the distance di,j(t) between the center of mass (COM) of ARi to ARj (Fig. 3 A):
Then the stiffness spring constant is given by
The torsional spring constant was computed by considering the fluctuations of the dihedral angle θi,j(t) formed between two repeats, ARi and ARj, along the filament. First, a vector was defined for each repeat pointing from its COM to the central Cα atom of the first helix of the repeat (Fig. 3 D). This vector is approximately perpendicular to the long axis of the filament. The dihedral angle is then defined to be the angle between the two planes (one formed by with COM of ARj and the other formed by with COM of ARi). The torsional spring constant is then
We averaged the spring constants from ANK1–ANK3 to compute the material properties used to parameterize the finite element model.
Finite element model
For a linear, elastic, isotropic rod, the material properties are defined fully by the Young’s stretching modulus (E) and the shear modulus (G; Landau and Lifshitz, 1986), which relate the tensile stress (force per unit cross-sectional area) to the stretching strain (change in length) and the shear stress to the shear strain, respectively. For a rod of circular cross section A = πr2 and length L, E and G are related to ks and kθ as (Landau and Lifshitz, 1986; O’Reilly, 2017)
where I = πr4/4 is the area moment of inertia of a circular rod. All finite element calculations were computed using the software Abaqus (2014) with the protein modeled as a linear, elastic, isotropic material defined by E and G. The shape of each AR chain from the cryo-EM structure (Jin et al., 2017) was fit with a fourth-order curvilinear polynomial to reproduce the chain geometry. Cross-sectional slices through each stretch were visualized to determine an effective radius of the chain (Fig. S3). The images revealed roughly circular geometries with radii between 0.75 and 1.2 nm, and we chose a uniform value of r = 1.0 nm for all computations. Calculation of the elastic moduli from the MD simulations required this cross-sectional area and the length of internal AR2–AR5 distance, which was set to L = 30 Å. In Abaqus we used eight-node, linear brick elements (C3D8R) with reduced integration and hourglass control within the standard solver. Additionally, we employed the nonlinear geometry feature and a full Newton–Raphson method with unsymmetric matrix storage. The maximum number of time increments was 1,000,000 with a minimum time increment of 10−8 and a maximum value of 0.1. The total time was 1.0 (no units), and we used the default convergence criteria. The C-terminal end of the chains terminates into a rigid block at the approximate position of the linker helices and TRP domain. This block is free to rotate in the xy plane consistent with TM motion in a fluid lipid bilayer, but it cannot move vertically. The N-terminal ends of the four helices were moved up and down vertically along the z axis to mimic compression and extension of the helical bundle while being fixed in the xy plane as if affixed to an MT. We explored three contact models: (1) no contact between the chains; (2) frictionless contact, in which adjacent chains slide on each other with no resistance when they meet; and (3) rough contact, in which chains are fixed at the point of contact once they meet. These boundary conditions were defined in Abaqus by prescribing rules for the normal and tangential forces at points of contact. In the no-contact model, adjacent chains do not see each other, while the frictionless and rough contact models experience a hard contact force in the normal direction upon contact. In the tangential direction, the frictionless model experiences no force allowing slip without energy loss, while the rough contact model experiences a large force damping out any relative motion. Given that the only external force present in the system is in the vertical direction, all the internal forces in the x and y directions cancel out. This implies that there is no displacement of the COM in the xy plane, yet internal forces in the xy plane will generate rotational moments about the z axis; i.e., the helical bundle is free to rotate about the z axis. Furthermore, under the assumption that the TM and the TRP region are rigid, there is insignificant deformation of the plate, and hence, the top of each chain is radially constrained in the xy plane.
Online supplemental material
The supplemental text describes a model for ankyrin chains as idealized helical rods. We provide a general mathematical framework of helical rod theory and use the theory to describe how the ankyrin chains behave under small and large deformations. The idealized helical model is used to support the finite element results obtained in Abaqus (described in the main text). Fig. S1 shows the geometry of an NOMPC ankyrin chain as an idealized left-handed helix. Fig. S2 shows the structure of the three stretches, called ANK1, ANK2, and ANK3, of six AR repeats used in our molecular simulations. Fig. S3 shows the heterogeneous shape of the cross-sectional area of ANK1, ANK2, and ANK3. Fig. S4 shows the TRP domain motion in the lowest-frequency normal modes (1–10) of the full, free NOMPC channel. Fig. S5 shows the TRP domain motion in the lowest frequency normal modes (1–10) of the full, clamped NOMPC channel. Fig. S6 shows the TRP domain motion of the lowest frequency normal modes (1–10) of the NOMPC TM domain. Fig. S7 shows the stability of ANK1, ANK2, and ANK3 from the molecular dynamic simulations.
Results
NMA suggests key gating motion
We first computed the normal modes of the full, NOMPC channel to determine the structure’s low-frequency, collective motions (Levitt et al., 1985; Bahar et al., 2010). We employed a coarse-grained representation of the protein only relying on the Cα backbone atoms, but it has been well documented that NMA performed on these chemically reduced structures still provides meaningful insight into the functional motions of the protein and serves as a good model for interpreting class averages observed in cryo-EM images (see Bahar and Rader [2005] for a review of the method). The lowest-frequency mode is fourfold symmetric, and Video 1 shows the full range of motion, while Fig. 2 A displays the extremum (cyan and black backbone) of one AR chain superposed on the starting structure after aligning the TRP domain from each snapshot. The largest deviations are localized to the bottom of the chain (N terminus) where we observe a vertical displacement along the z axis of ∼8 Å, but this displacement is transmitted weakly to the top of the chain where there is almost no vertical displacement of the last AR adjacent to the linker helices. Instead, upon compression (cyan) the top of the AR chain moves CCW 1–2 Å when viewed from the cytoplasm. This flexibility predicted by NMA is consistent with the cryo-EM class averages. The largest difference between classes 1 and 2 also occurs in the N-terminal region of the AR domains; however, the deflection is more moderate at 3 Å. It is not surprising that the absolute magnitudes differ as the class averages more likely represent energy minima rather than the extreme deflections reported for the normal mode. Additional support for the flexibility in the N-terminal region comes from the inability to resolve this section in class 3. Nonetheless, the lowest 10 modes fail to show any significant deformation in the TM domain, which we determined by aligning the extremum from each mode on the TRP domain and visualizing the differences, as illustrated in the inset to Fig. 2 A for the lowest mode and in Fig. S4 for all 10 modes. Thus, the lowest-order modes show that the AR domains are the most flexible part of the protein, but they do not lend any insight into how conformational changes in the helical bundle couple to the TM domain.
While these NMA calculations most closely match the cryo-EM conditions, we wanted to explore how clamping the N terminus to an MT altered the modes. To do this, we anchored nine ARs in space as if fixed to an MT by modifying the Hessian matrix as described in the Materials and methods section. We then computed the 10 lowest, nontrivial eigenvectors and eigenvalues, and Video 2 of mode 1 shows that the AR chains undergo compression–extension along the z axis in a fourfold symmetric manner. The two most extreme configurations (cyan and black) from this mode are represented in Fig. 2 B. Unlike the free protein, the motion along the z axis (as viewed by AR 29 adjacent to the channel domain) is severely dampened and almost nonexistent. The motion primarily takes place between the contact point 2 and the linker helices, with either end being largely fixed. We also investigated the motion of the TRP domain in the 10 lowest modes (Fig. S5), and as in the free protein case, the TRP domain moves as a rigid body for all modes (Fig. 2 B, inset). Next, we performed NMA on the isolated TM domain together with the TRP domain and linker helices to attempt to identify intrinsic motions in the channel masked by the soft, helical bundle. In the first nine modes, the TRP domain exhibits asymmetric motions incompatible with the fourfold symmetry in the channel (Fig. S6). However, mode 10 revealed symmetric, torsional displacements of the TRP domains similar to the hypothesized TRP domain gating motions (Cao et al., 2013; Jin et al., 2017). The TRP domain pivots at its C-terminal end, and the two extreme configurations of the mode (cyan and black) are characterized by a rotation of the helix about this pivot point. In the context of the full channel, this mode is certainly higher frequency than the z-axis displacements of the helical bundle, and while this result is suggestive that large compression–extension of the AR domain couples to torques in the TRP domain, it is not possible to say from this analysis.
Material properties of the AR domain determined from MD simulations
NMA analysis is a linear theory that describes motions around a starting configuration; however, it has been suggested that NOMPC may undergo very large deformations (Jin et al., 2017), and single-molecule experiments have shown that AR repeats can undergo large deformations without unfolding (Lee et al., 2006). In an attempt to understand how the helical bundle behaves for large conformational changes, we decided to construct a multiscale model of this portion of the protein using all-atom molecular simulation to determine the material properties of small sections of the chain, which would then be used to parameterize a full, finite element model of the helical bundle. Similar approaches have been used in the past to extract the material properties of proteins from all-atom MD simulations (Matsushita et al., 2010). Additionally, others have constructed finite element models of mechanosensitive channel gating for the stretch-activated MscL channel (Tang et al., 2006). Most similar to our approach, the Sun laboratory constructed multiscale models of protein mechanics using elasticity theory with moduli parameterized from atomistic simulations (Choe and Sun, 2005, 2007).
To start, we performed full-atom MD simulations on three different stretches of one of the AR chains (ANK1 [blue], ANK2 [red], and ANK3 [green] in Figs. 1 A and S2) to extract the stiffness spring constant ks and the torsional spring constant kθ of isolated elements. We examined three distinct stretches to determine if differences along the length of the chain or sequence differences, present in all ARs, significantly affect their mechanical properties. Ideally, we would have chosen stretches evenly spaced from the N terminus to the region just before the linker helices, but the density in the N-terminal region is rather low, and we decided it was best to simulate AR stretches that are well resolved. Thus, the most distal region we felt confident in simulating was ANK1, and there are several buffer repeats between ANK3 and the linker helices. Each stretch was six ARs long to keep the total atom count small allowing us to simulate the system for 450 ns in a reasonable time frame. We calculated the variance in the length and in the twist angle of each stretch between the most distant unconstrained ARs (2 and 5; Fig. 3, A and D). All three stretches exhibited the greatest fluctuations over the first 150–200 ns; however, the fluctuations were most pronounced in ANK1 (blue) with ∼5-Å deviations followed by ANK2 (red) and ANK3 (green), which showed variation between 2 and 2.5 Å (Fig. 3 B). The RMSDs of the entire stretch were 1.5–3.0 Å for all three, but much of this deviation comes from the linker loops as the values computed on the helices in each AR were much smaller (Fig. S7). The mean end-to-end distance for each stretch was between 29 and 30 Å with ANK2 starting out at a larger distance before settling down to a slightly more compact form.
We used the stabilized, absolute magnitude of the fluctuations over the last 350 ns (∼2 Å for all stretches) to calculate the stiffness constants for all three stretches. Using Eq. 2, we computed ks from the natural fluctuations in the protein (for our exact approach, see Matsushita et al. [2010] and the Materials and methods section), and these values ranged from 1,200 to 2,900 pN/nm. Despite the nearly 2.5-fold change in stiffness, the CIs for each stretch encompass the average ks = 2,100 (95% CI, 3,100 to 1,340) pN/nm. Next, we computed the torsional spring constant by monitoring the torsion angle θ2,5 between AR5 and AR2 of each stretch (Fig. 3 E). In the three systems, the values of θ2,5 are negative, indicating the left-handed geometry of the AR chain (Fig. 3 E). While the lengths of the three stretches were all within 1–2 Å, the degree of twist varies with ANK3 being ∼20° overwound compared with ANK2 and ANK1. The θ2,5 value of ANK1 shows greatest fluctuations of ∼40°, but the magnitude is reduced after 100 ns. Nonetheless, kθ computed from Eq. 3 reveals values of 65, 130, and 80 kBT/rad for ANK1–ANK3, respectively (Fig. 3 F), and the mean value across all stretches is kθ = 89 (95% CI, 132 to 53) kBT/rad. While our computations revealed differences in the stiffness values for the individual stretches, the CIs were quite large with a good deal of overlap, especially for the torsional spring constant. Based on this realization, together with the simplification that comes from treating each helical rod as a homogenous material, we chose to parameterize the finite element model with the average ks and kθ values mentioned above.
A finite element model of the helical bundle
We next constructed a continuum, finite element model of the AR four-helix bundle to explore its mechanical response to large deformations up to and beyond 30 Å. The continuum description of the protein is characterized by E and G, and assuming that the chain behaves as an isotropic, linear elastic material, E = 2.0 (95 CI, 2.9 to 1.3) GPa and G = 0.8 (95% CI, 1.0 to 0.4) GPa follow from the torsional and stiffness spring constants determined from our MD simulations according to Eq. 4. For these moduli, we measured the Poisson ratio to be ν ≈ 0.3, the value used by Tang et al. (2006) to study MscL. Next, we constructed a 3-D model of the helical bundle in Abaqus (2014) by parameterizing the shape of each chain mathematically based on the cryo-EM structure. The chains terminate at the linker helices and TRP domain (collectively called the TRP region), which we modeled as a rigid plate (Fig. 1 B). It has been suggested that the NOMPC channel responds to soft touch by compressing the AR domain like a spring (Liang et al., 2013), and we modeled this situation by fixing the N-terminal end of the chains in the xy plane and displacing them vertically in small increments, while maintaining the z coordinate of the TRP region fixed in space (Materials and methods section). Nonetheless, the TRP region was free to rotate in the xy plane.
Each AR chain is a left-handed helix, and the mechanical properties of inextensible, unshearable, constant pitch helical rods was worked out analytically by Love in the 1890s (see Supplemental material; Love, 2013). However, there are two complicating factors in NOMPC that limit the usefulness of this previous work. First, the helices are not constant pitch, as can be seen in Fig. 1, and second, when the helices contact each other, the system becomes four coupled springs that are no longer analytically tractable. As we will see, this first point is minor, but the second point can profoundly influence the results. Thus, using the finite element model we examined three scenarios to probe different possible outcomes: (1) a no-contact model to determine what happens if the helices move past each other; (2) a frictionless contact model that does not allow penetration but does let the helices slide on each other without resistance; and (3) a rough contact model that does not allow helices to slide after contact, but rather the surfaces remain attached upon further compression. A previous finite element model of the MscL stretch-activated ion channel modeled helix–helix contact using a Lennard–Jones potential where the helices influenced each other before direct contact. This model is more physically realistic, but in the limit where the width of the attractive minimum is narrow, it should produce the same result, and we believe that the most important factor is whether the helices are able to slide or not once contact is made.
The final compressed equilibrium configurations (Fig. 4, A–C) and corresponding videos (Videos 3, 4, and 5) reveal how the chains deform and the TRP region moves for the different models. The no-contact model exhibits a linear force–compression profile (reported value is for one chain) over the entire 3-nm deformation (black dashed line in Fig. 4 D) consistent with previous simulations and single-molecule experiments on isolated AR chains (Sotomayor et al., 2005; Lee et al., 2006). The effective spring constant for a single chain obtained from the slope yields a similar value to the analytic result (keff = 5.5 pN/nm versus ktheory = 4.0 pN/nm), indicating that the nonconstant pitch of a single helix is not a crucial determinant of its isolated mechanical properties. Additionally, when we applied our procedure to a 24-AR chain, we obtained a stiffness constant of 4.7 pN/nm in exact agreement with the MD simulations of Sotomayor et al. (2005). Our predicted, effective spring constant for the entire bundle is 22 pN/nm, which is simply four times the single chain value. While the close correspondence between the no-contact model and the analytic results validates our numerics, we can immediately see that this simple model fails as the final configuration (dark gray) shows major overlap of all of the chains (Fig. 4; light, transparent structure is the starting configuration in all panels).
Over the first 15 Å of compression, there is no contact between helices in either the rough or frictionless models (gray zone in Fig. 4, D and E), and thus, the force–displacement and energetics of all three models are identical (Fig. 4, D and E). Helix–helix contact first happens at contact point 2 followed by contact at point 1 at 18 Å compression. Reported force values are per chain, so the total force required to compress the bundle 15 Å is 4 × 9 = 36 pN. The reported energy is for the entire bundle, and it only requires 6.5 kBT to reach the first contact point. Thus, while thermal fluctuations are expected to cause compressions up to 7 Å, they are unlikely to be large enough to cause helix–helix contact. Once contact occurs at point 2, both the rough and frictionless contact models depart from simple linear behavior. The interchain interactions increase the overall effective spring constant of the helical bundle beyond what is expected for four independent AR helical springs. However, the change in the slopes indicates that the frictionless model (blue curve) is only marginally stiffer, while the rough contact model experiences an abrupt change and becomes much more stiff (red curve, Fig. 4 D). Careful examination of the rough contact force–displacement curve actually reveals three zones of increasing stiffness, before contact at point 2 (Δz < 15 Å), after point 2 but before point 1 (15 Å < Δz < 18 Å), and after both contacts (Δz > 18 Å). This means that impeding helix–helix sliding causes the ankyrin chain to lock up, becoming increasingly harder to compress, or in other words, when the chains stick together at the contact points, the effective stiffness of the system increases significantly. This difference in the effective stiffness of the helical bundle has an important consequence for how hard it is to compress the chains by 30 Å, for instance, requiring only 25–30 kBT of energy in the frictionless case but up to 65 kBT for the rough contact model (Fig. 4 E).
One striking feature evident from the final snapshots and the videos is that the channel domain (TRP region) rotates in the membrane during compression. Initially, all three models rotate CCW when viewed from the cytoplasm (Videos 3, 4, and 5), and since the N terminus of each helical chain is fixed in the xy plane, CCW rotation indicates that the chains are overwinding: they are acquiring more helical turns. Overwinding is a direct consequence of imposing compression and end-rotations to a helical rod (see Fig. S1 for geometry), and in general, a helical rod with a fixed radius R subjected to a compression force along its helical axis will naturally tend to overwind as we show analytically in the small deformation limit (Supplemental material). As a consequence, the no-contact model continues to overwind as the compression increases; however, the other two models exhibit more complex behavior and switch directions leading to underwinding. A helical rod under compression can only underwind—decrease the helical turns—if the helical radius R expands. After contact in the frictionless model, the chains start to underwind because chain–chain contact creates a positive radial force on neighboring chains that expands the helical radius R (Video 4). In fact, the in-plane rotation is greatest for this model, and one can readily see the expansion of the entire bundle in the plane where contact site 2 occurs. On the other hand, in the rough model the chains stick together after contact resulting in a radius R that remains approximately constant upon further compression. Hence, the rough model exhibits small rotations of the TRP region that oscillate back and forth between CW and CCW modes (Video 5). Please note that since these rigid rotations simply rotate the channel in the fluid, lipid membrane, they are unlikely to influence the gating of the channel.
Force and torque transmission to the TRP region
Next, we analyzed how the force from the chains is transmitted to the linker helices and TRP domain. To do this, we integrated the stresses over all nodes connecting the helices to the TRP region and decomposed the net force into its radial and tangential components as illustrated in Fig. 5 A (all quantities in Fig. 5 are reported for a single helical chain). Prior to helix–helix contact at contact point 2, which occurs near 15 Å, all three contact models exhibit a positive radial force F|| that would cause the TRP region to radially expand, potentially opening the S6 bundle crossing (Fig. 5 B). This expansion does not occur here because the TRP region is rigid. As the helical bundle continues to compress, the radial forces become smaller in the frictionless and rough contact models, and the rough contact model reverses sign above 20 Å, which may cause the S6 bundle crossing to start to close. Meanwhile, each AR chain exerts a negative tangential force F⊥ on the TRP region before helix–helix contact (gray region in Fig. 5 C), which is expected since the long axis of each helix does not perfectly align with the membrane normal (Fig. 1 B) and vertical displacements of the N terminus impart a force to the TRP region along this direction due to the radial constraint in the xy plane. Thus, each helical chain initially exerts a CCW torque on the TRP region at the center of the chains, which must be counteracted by other torques in the system, discussed next, to remain in static equilibrium for a given compression because the entire TRP region is free to rotate. Once contact occurs in the frictionless model the perpendicular force exerted by the chains switches signs (blue curve in Fig. 5 C). Finally, the magnitude of the in-plane forces are ∼10 times smaller than the out-of-plane forces in Fig. 4 D, but they may evoke gating motions important for opening the pore domain.
Next, we calculated the twisting moment mz on each TRP domain that balances the rotational torque imposed by the tangential forces F⊥ (Fig. 5 D) such that the system remains in static equilibrium. Prior to contact (gray region in Fig. 5 D), all three models exert a positive torque on the TRP domain, which would cause CW rotation of each individual region when viewed from the cytoplasm. For displacements >15 Å, the rough contact model continues to exert a positive torque on the TRP segment, but the frictionless model reverses direction and starts to exert a negative torque on each TRP segment that would cause CCW rotation of the TRP domain. The subtle rotations observed in the linker helices when the chains are compressed from class 2 to class 1 are consistent with an applied CW torque (Jin et al., 2017), and a CW torque is also consistent with the mode of opening suggested by the closed and open TRPV1 channel structures (Cao et al., 2013).
Discussion
We constructed a multiscale model of the NOMPC channel using different levels of theory to investigate how external forces influence the conformation of the channel and transmit forces through the helical bundle to the channel domain. We first employed NMA to explore the natural motions of the starting structure. The lowest-frequency mode of NOMPC shows that the helical bundle compresses and extends along the z axis by several angstroms, consistent with the structure being a soft spring. The room temperature amplitude of the fluctuations predicted from NMA is larger than those observed in the cryo-EM class averages. Calculations reveal that the largest motions are localized to the N-terminal ARs (Fig. 2 A), consistent with poor density in the class 3 averages (Jin et al., 2017). The lowest modes of the entire channel fail to show any significant movements in the channel domain, consistent with the notion that the most accessible, lowest-energy deformations first occur in the helical spring, and much more energy is required to cause conformational changes in the channel domain needed to gate the channel. In this regard, the helical bundle is an actuator that will undergo deformations due to thermal noise, but these room temperature fluctuations are not enough to gate the channel. Clamping the N-terminal ARs revealed no further information about gating. When we examined the normal modes of the TM domain alone, we saw more complex motions in the channel. Most of these modes are asymmetric, but normal mode 10 is the first mode that exhibits a fourfold symmetric motion, wherein the TRP domain rotates in a manner observed between the open and closed TRPV1 channel structures (Cao et al., 2013). Thus, the intrinsic motions of the channel-only construct suggest that TRP domain rotation, while much higher frequency than compression–extension in the helical bundle, may be accessible with enough energy.
In an attempt to access these higher-energy modes of the channel, we created a finite element model of the helical bundle so that we could explore large-scale deformations in the chains. This continuum framework provides a reduced dimensionality, similar to the NMA analysis, making it computationally tractable to calculate both deformations and the energetics of these conformations. We employed a multiscale approach by parameterizing the elastic moduli from all-atom MD simulations on smaller segments of the full AR chain. The biggest uncertainty in our modeling lies in how the helices behave after making contact at points 2 and 1. That said, there are robust conclusions concerning the mechanics of compression before contact that are model independent (stage 1 in Fig. 6). The energy required to compress the helical bundle by 15 Å is only 6.5 kBT (Fig. 4 E), and the corresponding total force is 36 pN (Fig. 4 D). After contact, the effective spring constant dramatically increases in the rough contact model, while it increases only slightly in the frictionless model providing a nearly linear force–displacement curve. Compression of the helical bundle induces overwinding and underwinding of the AR chains, which causes rigid body rotations in the plane of the membrane (stage 1, central panel). While these movements differ from model to model, we do not think that this motion is critical to channel gating. Instead, the most important aspect is that each chain exerts a torque on the linker helices/TRP domain, which would cause it to rotate. As a consequence of the left-handed helical geometry of the NOMPC AR chains alone, this torque is initially in the CW direction for all models when viewed from the cytoplasm (stage 1, right panel), and it would cause the TRP domains to adopt a configuration observed in the open TRPV1 channel (Cao et al., 2013). As the compression continues beyond 15 Å (stage 2 in Fig. 6), the torque on the TRP domain switches directions for the frictionless model but not for the rough contact model. If TRP domain movement is the key to channel gating, then we predict that the frictionless contact model would actually close the channel under high compressions. Additional experiments are needed to determine if the channel deactivates under high compressive loads, but if it does, this finding supports the frictionless contact model. Finally, the helical chains do not adopt a uniform circular cross section with constant radius throughout (Fig. S3), and additional calculations with a 0.75-nm-radius rod revealed that the magnitude of the forces and torques were smaller, but the signs were similar, and all of the deformations and shapes were qualitatively similar.
While our modeling efforts have laid out a number of scenarios with clear, testable predictions, the biggest uncertainty is the nature of the helix–helix contact. It is difficult to predict how these elements will slide on each other by simply inspecting the nature of the amino acid composition, and instead, we believe that the best approach is to attempt cross-linking between the patches at contact points 1 and 2 to determine if there is a strong influence on channel behavior. If cross-linking minimally perturbs channel function, it supports the rough contact model, while the frictionless contact model is favored otherwise. Another potential experiment for discerning which model is most correct is to measure the force–displacement curve using single-molecule methods. We predict that the force–displacement curve for the rough contact model will be nonlinear with a soft and stiff regime, while the frictionless contact model will be linear throughout. Unfortunately, it may be difficult to compress the helical bundle, as previous single molecule experiments on AR chains were only able to extend chains beyond the resting equilibrium length (Lee et al., 2006). Moreover, some studies have suggested that channel opening may involve extension of the bundle (Yan et al., 2013; Zhang et al., 2015). Importantly with respect to extension, we used our finite element model to pull on the helical bundle, and the contact points never touch. Hence, all three models are the same under tension, and they give very similar results to the analytic model for four independent chains in parallel. Thus, up until the contour length of the 29 ARs is reached and they start to unfold, we expect pulling experiments will reveal a linear force–displacement curve with a 22 pN/nm stiffness constant based on the calculations in Fig. 4 D.
There are also additional questions regarding the function of the channel in vivo. For instance, at rest in the hair cell, are the AR repeats under compression, under tension, or at their equilibrium length? The resting open probability of the channel and how changes in the compression–extension of the helical bundle shift this open probability would depend on the resting state of the system. Additionally, it is not known how much energy or force is required to open the channel. If this value was known, it would help interpret where along the energy–displacement curve in Fig. 4 E we expect to see opening, which likely would shed light on the degree of bundle deformation required for opening and on how the mechanical forces are being transmitted to the pore domain at opening (Fig. 5). Our preliminary calculations using our continuum membrane-bending model (Argudo et al., 2017) suggest that the ankyrin chains are much softer than the membrane, and most of the deformations would be localized to the helical bundle, as previously hypothesized by the finite element model of the distal tip constructed by the Howard laboratory (Liang et al., 2013). This result is important, as the greater compliance of the helical bundle is likely essential for generating the required forces needed to gate the channel, rather than simply deforming the membrane. Previous studies have shown that deleting the first 12 ARs leads to higher surface expression of the channel and a higher open probability, while deleting the last 14 or 17 ARs abolishes surface expression (Zhang et al., 2015). Moreover, doubling the AR repeats (29 + 29) or adding additional ARs in between AR29 and TM1 lead to channel expression with spontaneous activity, but only the 29 + 29 AR channel showed mechanosensitivity. Our continuum model is ideal for quickly exploring these additional channel constructs to determine how the forces and torques on the linker helices/TRP domain change with the changes in the linker length.
Finally, we do not know how the forces predicted from our continuum model would actually impact the rearrangements of the helical bundle and TRP domain. Is TRP domain rotation important, or do the large vertical and radial forces exerted by the AR chains gate NOMPC? In the future, we intend to use all-atom MD simulations in the absence of the helical bundle to determine how rearrangements of the TRP domain and linker helices impact ion flow through the inner gate. We will use the forces and torques predicted from the finite element model to guide these simulations.
Acknowledgments
We thank P. Jin, D. Bulkley, Y.-N. Jan, and Y. Cheng for fruitful discussions; O. Choudhary for his support using ProDy for the NMA; and T. Lang for use of Abaqus.
This research was supported by the University of California, San Francisco, Program for Breakthrough Biomedical Research to M. Grabe; an American Heart Association predoctoral fellowship to N.P. Bethel (17PRE33410225), and a National Institutes of Health postdoctoral fellowship to D. Argudo (4T32HL007731-25).
The authors declare no financial competing interests.
Author contributions: D. Argudo, N.P. Bethel, and M. Grabe conceptualized the research, and D. Argudo, S. Capponi, and N.P. Bethel performed the calculations/simulations. All authors analyzed the data and made the figures. S. Capponi, D. Argudo, and M. Grabe wrote the paper.
José D. Faraldo-Gómez served as editor.
References
This work is part of the special collection entitled "Molecular Physiology of the Cell Membrane: An Integrative Perspective from Experiment and Computation."
Author notes
D. Argudo and S. Capponi contributed equally to this paper.