­­Mitochondrial ATP synthase dimers spontaneously associate due to a long-range membrane-induced force­

ATP synthases play a central role in cellular bioenergetics. Located in the plasma membrane of bacteria, archaea, and cyanobacteria, or the innermost membranes of chloroplasts and mitochondria, these ubiquitous enzymes harness transmembrane electrochemical gradients of H⁺ or Na⁺, generated from respiration or light harvesting, to produce the majority of ATP used by the cell.

ATP synthases consist of four functionally distinct elements: a transmembrane domain, a catalytic domain, and a central and peripheral stalk. Protons or sodium ions permeating through the transmembrane region, down their electrochemical gradients, cause the rotation of a subcomplex known as the c ring. This rotation is transmitted to the catalytic region via the central stalk, which elicits a cycle of conformational changes that result in the synthesis and release of ATP. These changes in the internal structure of the catalytic domain are possible because the peripheral stalk holds the catalytic region stationary relative to the transmembrane region. Thus, the exergonic process of ion permeation is efficiently coupled to the endergonic process of ATP synthesis.

Although the architecture of the ATP synthase complex is largely conserved across species, the mitochondrial enzyme is distinct in that it exists as a dimer. This dimer is formed through accessory transmembrane subunits not present in prokaryotes or chloroplasts (Arnold et al., 1998; Lapaille et al., 2010). The dimers of mitochondrial ATP synthases are unique compared with other membrane protein oligomers in that the two transmembrane domains are tilted relative to each other, approximately at a right angle (Allegretti et al., 2015; Hahn et al., 2016; Guo et al., 2017; Klusch et al., 2017). This arrangement causes the ATP synthase dimer to have a striking V-like shape, with the two catalytic domains projecting away from one another. In addition, these complexes are also remarkable in that they are not scattered randomly throughout the inner mitochondrial membrane; instead, they are organized into long rows that run along the sharply curved edges of the membrane invaginations called cristae (Strauss et al., 2008; Dudkina et al., 2010; Davies et al., 2011, 2012; Mühleip et al., 2017). These rows extend for hundreds of nanometers and include dozens of dimers, arranged side by side (Fig. 1).

The factors that drive and stabilize these supramolecular structures have been unclear. Using molecular dynamics simulations, we have previously shown that an isolated ATP synthase dimer induces a pronounced, long-range deformation in the surrounding membrane, owing to its V-like shape (Fig. 1; Davies et al., 2012). This curvature deformation occurs both parallel and perpendicular to the dimer interface. In contrast, when multiple dimers were placed side by side and kept a distance similar to that observed in mitochondria, we observed that the membrane flattens in between adjacent dimers, i.e., it restores its natural shape along the direction of the row. Based on a simple elastic model (Helfrich, 1973; Marsh, 2006), we inferred from these different curvatures that the membrane might drive the formation of the ATP synthase dimer rows to reduce its overall deformation (Davies et al., 2012). However, direct evidence that this effective force in fact emerges from the interatomic interactions and dynamics of the molecular system was not obtained. Here, we present free-energy molecular simulations designed to explicitly probe the existence and magnitude of this hypothetical membrane-induced attractive force. We also discuss the implications of our findings and how to test them experimentally.

Two molecular systems were used in this study, both based on the MARTINI 2.1 force field (Marrink et al., 2007). The first system, or system A, consists of a bilayer of 5,832 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipids in a 100-mM NaCl solution. The dimensions of this bilayer are 45 × 45 nm, and the height of this system, including the two solvent layers at either side of the membrane, is 17 nm. In total, system A comprises ∼300,000 coarse-grained (CG) atoms. The second simulation system, or system B, includes two ATP synthase dimers embedded in a much larger lipid bilayer comprising 68,000 POPC molecules. The structural model for the ATP synthase dimer has been described previously (Davies et al., 2012). That is, it is a composite of x-ray structures for the F₁ domain and the c ring, transformed into a CG representation; and for the remainder of the transmembrane domain, a network of hydrophobic CG particles filling up a density envelope derived from single-particle analysis of cryoelectron tomograms of mitochondrial membranes. The shape of this envelope is not identical but closely resembles more recent cryo-electron microscopy (EM) structures of the enzyme. The dimensions of the lipid bilayer for this two-dimer system are 85 × 254 nm, and the height of the system is 33 nm. In total, this second simulation system amounts to ∼6 million CG atoms.

All molecular dynamics simulations were performed with GROMACS 4.5.5 (Hess et al., 2008) and PLUMED (Bonomi et al., 2009). The integration time step was 10 fs. The simulations were performed at 300 K and 1 bar. A velocity-rescaling algorithm with time constant 0.3 ps was used to preserve the system temperature. The pressure was sustained semiisotropically using the Berendsen method, with a time constant of 3.0 ps. Nonbonded interactions were calculated using a Coulomb potential, with relative dielectric constant ε = 15, and with a Lennard-Jones potential; a smooth shifting function was used so that both potentials are zero at 1.2 nm.

A molecular dynamics simulation of system A was used to estimate the characteristic bending modulus k_c of a symmetric POPC membrane, as represented by the MARTINI 2.1 force field. The bending modulus can be derived from the undulation spectrum of the membrane, which in turn can be calculated from a molecular dynamics trajectory. Specifically, in the low-frequency regimen, the amplitudes of the different fluctuation modes of the membrane, u, fulfill the relationship (Helfrich, 1973; Marrink et al., 2004):

where < > denotes a time average, A is the surface area of the membrane, q is the mode wave number (inverse of the wavelength), T is the temperature, and k_B is the Boltzmann constant. A trajectory of 275 ns was calculated for this purpose. The undulation spectrum was then obtained from the last 150 ns of this trajectory. Specifically, a Fourier transform of the z positions of the phosphate groups in each leaflet was calculated for each trajectory snapshot, and the resulting amplitudes were averaged over the trajectory. The value of k_c was then estimated by fitting the function shown above to the calculated spectrum in the low-frequency regimen.

To determine whether the membrane alone can drive the association of two ATP synthase dimers, simulations of system B were performed to calculate the change in free energy as a function of the distance between the two dimers, r. This one-dimensional free-energy profile, or potential of mean force, was calculated using the umbrella-sampling method. The total simulation time was 15 μs. Specifically, 200 independent simulations were performed spanning a range of r values from 2.5 to 50 nm, in increments of 0.25 nm (additional simulations were introduced in smaller increments in some cases to guarantee overlap of the calculated probability distributions). In each simulation, the value of r was maintained in the vicinity of the target using a harmonic restraining potential of force constant 95.6 kcal/mol nm⁻². To simplify the calculation, the two dimers were kept at all times approximately parallel to each other. To do so, we monitored the distances between the centers of mass of opposing c rings and applied an additional restraint that minimizes their difference, using a harmonic potential of force constant 28.7 kcal/mol nm⁻². The starting configurations for this series of simulations were produced sequentially; i.e., for a given simulation, the initial snapshot was extracted from the preceding simulation along r. Each simulation comprised a relaxation stage of 37.5 ns, followed by a production phase of 70 ns. Snapshots were collected every 1 ps of simulation time; i.e., >15 million snapshots were analyzed. The free-energy profile was calculated using the DHAM algorithm (Rosta and Hummer, 2015). To estimate the statistical error, the simulation data were split into two sets, and the mean difference of the corresponding potential of mean force profiles was calculated.

The elastic deformation of the membrane induced by the two ATP synthase dimers was quantified as described previously (Davies et al., 2012). In brief, for a given position in the plane of the membrane (x, y) and a given simulation snapshot, we evaluated the instantaneous membrane curvature in the x and y directions by calculating the radius of the circle that best fit the membrane surface at that position, in either direction. The reciprocal of this radius value is the instantaneous local curvature, c_x or c_y. The mean curvature c is the sum of these two local curvatures (Marsh, 2006), time averaged over the simulation length, i.e., < c > = < (c_x + c_y) >. To calculate the membrane curvature for each snapshot of the simulation, the phosphate groups of the lipid molecules were mapped on a two-dimensional grid in the (x, y) plane (0.5-nm grid-point spacing), based on their proximity to each grid point, using a cutoff distance of 9 Å. The circular fitting was then performed at each grid point of the average surface following the method of Forbes (1989).

To evaluate whether ATP synthases are organized by a membrane-mediated attractive force, we simulated the association of two ATP synthase dimers embedded in a phospholipid membrane of dimensions 250 × 85 nm surrounded by solvent (Fig. 2). The number of particles in this simulation model was ∼6 million. A series of 200 molecular dynamics trajectories was then calculated to probe the propensity of the dimers to move closer together or farther apart as the distance between them, r, was gradually varied (Materials and methods). These simulations used the umbrella-sampling method, totaling 15 μs of sampling time. From a global analysis of these trajectories, we deduced the gain or loss in free energy resulting from the change in the dimer-to-dimer distance, i.e., the potential of mean force W(r).

The results obtained very clearly demonstrate that the two dimers are driven together by a strong, long-range force induced by the membrane (Fig. 3). Specifically, the calculated potential of mean force shows that this effective interaction begins to act at distances smaller than 40 nm (center to center), when the membrane deformations induced by each dimer begin to coalesce, gradually forming a ridge between the dimers (Fig. 3 A). The free energy of the molecular system decreases steadily as the dimers come closer together, reaching a minimum at 11–13 nm (Fig. 3 B). Up to this point, the dimers are never in contact, and direct physical interactions between them are negligible. Beyond this point, however, the free energy rises quickly as the dimers begin to interact, and ultimately their steric repulsion becomes dominant. Thus, these large-scale free-energy calculations show that the membrane alone brings the ATP synthase dimers together, up to an optimal noncontact distance where the two dimers form a membrane ridge, flat in cross section (Fig. 1).

Analysis of electron cryotomograms of mitochondrial membranes from Saccharomyces cerevisiae (Davies et al., 2012) shows that the typical distances between ATP synthase dimers observed in situ are highly consistent with those predicted by the theoretical model. Specifically, the most frequently observed distance between adjacent dimers in a row is 12 nm (Fig. 3 C). Thus, the membrane-mediated attractive force quantified in the computational model appears to optimally explain the organization observed physiologically.

The magnitude of the calculated free-energy minimum is also worth noting. To translate this value into a “bound-state” probability, one must integrate the one-dimensional profile W(r) to deduce a two-dimensional dissociation constant:

where k_B is the Boltzmann constant, T is the temperature, and R is a threshold in the dimer-to-dimer distance differentiating bound and unbound states. Note that the units of K_d in the above expression are molecule/nanometers squared. Although this expression is general, the quantity C is specific to our case, as it is a correction for the fact that throughout the calculation of W(r), the dimers are kept approximately parallel for computational convenience (Materials and methods). To a first approximation, C can be estimated as Δθ/2π, where Δθ is the rotational fluctuation of one dimer relative to the other, around the membrane perpendicular, in the associated state (note the rotational freedom of the unbound dimer is 2π, and not 8π², as in solution, on account of the orientational restrictions imposed on the dimer by the membrane itself). The calculated values of K_d and C from the simulated data are ∼3 × 10⁻¹⁵ nm⁻² and ∼0.03, respectively. Finally, the bound-state probability is

where D is the area density of dimers in the membrane.$$For example, the system depicted in Fig. 2 approximately corresponds to a mole-fraction of 5 × 10⁻⁵ (one dimer per 20,000 lipids), which translates into D ∼7 × 10⁻⁵ nm⁻² (for an area per lipid of ∼70 Ų). At this density, the population of the dissociated state would be virtually zero—specifically, 1 – P_b ∼10⁻¹⁰.

Several factors might modulate this association process in a physiological setting. Biological densities might be significantly lower than what Fig. 2 implies, and changes in the lipid composition of each of the monolayers of the inner mitochondrial membrane (Horvath and Daum, 2013) might influence the magnitude of the attractive force between dimers. It should also be noted that the multidimer free-energy profile describing the assembly of a row would be naturally shallower than that describing the association of two dimers only. These considerations notwithstanding, our data indicate that isolated dimers will be extremely rare in any type of membrane.

At the present time, the question posed in this study cannot be quantitatively evaluated using a fully atomistic representation of the molecular system. Hence, we have resorted to a CG representation, specifically that implemented in the MARTINI force field. Whether this widely used approach is sufficiently accurate to study close-range interactions between membrane proteins has been recently questioned (Javanainen et al., 2017). The process we described does not entail such protein–protein interactions, and therefore MARTINI appears to be a valid approximation in our view. Nonetheless, and given that the process we describe is driven by the membrane’s resistance to deformation, we asked whether MARTINI provides a reasonably accurate description of the elastic properties of the specific model lipid membrane considered in this study. We therefore performed a molecular dynamics simulation of a protein-free POPC bilayer (∼5,800 lipid molecules), and calculated its bending modulus from analysis of its undulatory spectrum (Materials and methods). As shown in Fig. 4, undulations of amplitude u have an intensity that is inversely proportional to the fourth power of the wavelength, q, in the long-wavelength regimen, i.e., q < 0.1 Å⁻¹, as predicted from elastic theory (Helfrich, 1973; Marrink et al., 2004). The bending modulus k_c can be derived from the slope of the line fitting the calculated intensity < u² > as a function of q in a logarithmic scale. The calculated value of k_c from this analysis is 1.2 × 10⁻¹⁹ J, which is in good agreement with the value experimentally determined for POPC, namely ∼10⁻¹⁹ J (Marsh, 2006). The CG MARTINI representation therefore appears to be a valid approach to evaluate the association process considered in this study, notwithstanding the simplifications inherent to this force field and its reported limitations for other types of processes controlled by direct protein–protein interactions (Javanainen et al., 2017).

It has long been assumed that mitochondrial ATP synthase dimers assemble into rows as a result of direct protein–protein interactions between dimers (Minauro-Sanmiguel et al., 2005; Dudkina et al., 2006; Thomas et al., 2008; Couoh-Cardel et al., 2010). This notion stems in part from analysis of two-dimensional EM images of detergent-solubilized dimers, which appeared to reveal two populations differing in the relative angle between monomers (Dudkina et al., 2006). This result seemed to indicate the possibility of conformational adjustment upon direct interaction between dimers. However, subsequent tomographic analyses of mitochondrial cristae have shown that the conformation of the dimer in situ is invariant for a given species (Davies et al., 2011, 2012). Moreover, these tomographic data are in agreement with three-dimensional structures determined through single-particle cryo-EM (Hahn et al., 2016). A likely explanation for this discrepancy is that the reported structural heterogeneity of the dimer (Minauro-Sanmiguel et al., 2005; Dudkina et al., 2006; Thomas et al., 2008; Couoh-Cardel et al., 2010) results from different two-dimensional projections of the same conformation. Accordingly, in our simulations, we study an invariant structure of the dimer throughout the range of distances considered.

The concept that dimer-to-dimer contacts occur also originates in biochemical and mutational studies indicating that subunit 4 (or b in bovine) in the peripheral stalk and subunits e and g in the transmembrane region are involved in oligomerization. These subunits were proposed to mediate two different dimerization interfaces, thereby generating the arrangement of the dimers observed in mitochondria (Everard-Gigot et al., 2005). Subunit 4 was specifically proposed to form dimer-to-dimer contacts through its soluble region (Spannagel et al., 1998). However, the more recently obtained high-resolution structures of the yeast ATP synthase dimer make very clear that subunits e and g and the N-terminal transmembrane helix of subunit 4 are located at the interface between the two monomers of a dimer (Hahn et al., 2016; Guo et al., 2017). Indeed, the interaction between these subunits is the reason for the striking V shape of the complex. The cryo-EM structures of the ATP synthase dimer also make clear that the peripheral stalks project away from each other, and thus the proposed interaction between the soluble regions of two subunits 4 cannot take place either within a dimer or between dimers. A probable explanation for this discrepancy is that the abovementioned cross-linking experiments were performed with samples lacking subunit e and g, i.e., for monomeric ATP synthases. It is thus questionable whether these interactions are biologically relevant.

Subunit e has also been suggested to be involved in stabilizing the ATP synthase dimer rows through the formation of a coiled-coil by the C-terminal helix (Spannagel et al., 1998; Brunner et al., 2002). However, this putative interaction also appears to be incompatible with the cryo-EM structures of the yeast enzyme, which clearly show this C-terminal helix projecting straight down into the intercristae space (Hahn et al., 2016; Guo et al., 2017). The fact that this element could be discerned even at a 7-Å resolution suggests that it is relatively rigid and thus unlikely to adopt an entirely different conformation. In summary, although the notion of direct dimer-to-dimer interactions is intuitive, the supporting biochemical evidence is, in our view, unconvincing, as it is incompatible with structural data more recently obtained through both high-resolution single-particle analyses and cryotomography of mitochondrial membranes. We therefore posit that subunits e, g, and 4(b) are important for the formation of the ATP synthase rows not because they mediate interactions between adjacent dimers, but simply because they are required for the dimerization of the enzyme, and more specifically, its V-shaped configuration. Accordingly, the preferred distance between dimers observed in both experiment and computation (around 12 nm) is much too large for a direct contact, or even a direct physical interaction (e.g., electrostatics).

The computational analysis presented here generates two clear-cut experimentally testable predictions: first, that ATP synthase dimers will self-associate into rows without requiring any other protein or a specific membrane component, ultimately creating a membrane ridge along the length of the row; second, that the driving force for the formation of dimer rows originates in the membrane deformation imposed by the V-like shape of the dimers. To test the first hypothesis, purified ATP synthase dimers could be reconstituted into tethered lipid bilayers or synthetic liposomes; the distribution of these dimers in these membranes could be then analyzed using high-resolution atomic-force microscopy or electron cryotomography. To test the second prediction, the same experiment could be conducted for ATP synthase monomers, e.g., constructs lacking the dimer-specific subunits e and g (Davies et al., 2012). In the first experiment, we expect that the dimers will self-associate into long rows; indeed, the magnitude of the attractive force computed here suggests that isolated dimers will be seen rarely or not at all. The second experiment, in contrast, should show no evidence of association or membrane organization.

All cellular membranes undergo constant shape changes. Gentle curvatures can occur through the asymmetric distribution of lipids, whereas acute curvatures require specialized proteins, such as clathrin, which coat the membrane surface and alter its morphology through the formation of a large oligomeric structures (McMahon and Boucrot, 2015; Jarsch et al., 2016). Mitochondrial cristae are formed from the repeated infolding of the inner membrane. This process is driven by lipid biosynthesis and the resulting increase in the surface area of the membrane. As this growth occurs in a confined space, i.e., within the outer membrane, the inner membrane buckles inward. However, membrane growth alone does not explain the morphology of cristae. Mathematical modeling of this growth process shows that although multiple buckling events might occur initially, fusion events would eventually lead to several large vesicle-like compartments (Renken et al., 2002; Kahraman et al., 2012). Mitochondrial ATP synthases are unlike most other membrane-bending proteins in that they are integral transmembrane complexes whose primary biological function, individually, is not membrane remodeling. As dimers, however, we propose that their biological role is also to prevent this kind of vesiculation and to foster instead flattened lamellar-shaped cristae, or other types of geometry that compartmentalize the increasing surface area of the inner mitochondrial membrane. We envisage that early on in the development of the mitochondrion, rows of ATP synthase dimers form at multiple locations in the inner membrane, priming these sites for invagination. As the membrane grows, the inner membrane preferentially buckles at these primed locations, carrying the rows of dimers into the matrix. Because of the strong membrane curvature generated by the row, the invaginations would adopt a flattened-disk geometry, and the presence of ATP synthase dimer rows along the edges would prevent unregulated fusion or division of the inner mitochondrial membrane (Fig. 5).

We have presented a molecular-simulation analysis designed to probe the factors that govern the membrane organization of mitochondrial ATP synthase dimers in vivo. The results presented explicitly demonstrate that there exists a force that drives the dimers together, side by side, to form or extend a row. This attractive force originates in the membrane deformation created by the structure of the isolated ATP synthase dimer. This deformation is energetically costly, but is gradually reduced as the dimers approach each other to form a ridge. The range of this effective attractive force is tens of nanometers, where no direct physical interactions between dimers occur. At the most favorable dimer-to-dimer distance, shown to be around 12 nm by both computations and experiments, protein–protein contacts cannot occur either. We therefore conclude that the formation of dimer rows is a process of spontaneous self-association, driven by the membrane itself. Importantly, this hypothesis is experimentally verifiable, for example through imaging of purified ATP synthase dimers and monomers reconstituted in liposomes or tethered bilayers. The mechanism we propose does not require a specific lipid composition, because all membranes resist deformation. It is essential, however, that the V shape of the dimer is intact. We further posit that rows of ATP synthase dimers prime the inner mitochondrial membrane to develop cristae and contribute to determining their shape and stability. Thus, these remarkable enzymes would not only produce most of the ATP consumed by the cell, but also help define the macroscopic morphology of the organelle in which they operate.

Computational resources were in part provided by the Biowulf Supercomputer at the National Institutes of Health, the Jülich Supercomputer Center, and the Partnership for Advanced Computing in Europe. We thank L.R. Forrest, V. Leone, and F. Marinelli for their comments on this manuscript.

This work was supported by the Max-Planck-Gesellschaft (C. Anselmi, K.M. Davies, and J.D. Faraldo-Gómez), the Cluster of Excellence “Macromolecular Complexes” of the Deutsche Forschungsgemeinschaft (K.M. Davies and J.D. Faraldo-Gómez), and the Intramural Research Program of the National Heart, Lung, and Blood Institute, National Institutes of Health (C. Anselmi and J.D. Faraldo-Gómez).

The authors declare no competing financial interests.

Author contributions: C. Anselmi: Investigation, formal analysis, and writing. K.M. Davies: Conceptualization and writing. J.D. Faraldo-Gómez: Conceptualization, supervision, writing, funding acquisition, and resources.

Lesley C. Anson served as editor.