Lipid-Loving ANTs: Molecular Simulations of Cardiolipin Interactions and the Organization of the Adenine Nucleotide Translocase in Model Mitochondrial Membranes.

The exchange of ADP and ATP across the inner mitochondrial membrane is a fundamental cellular process. This exchange is facilitated by the adenine nucleotide translocase, the structure and function of which are critically dependent on the signature phospholipid of mitochondria, cardiolipin (CL). Here we employ multiscale molecular dynamics simulations to investigate CL interactions within a membrane environment. Using simulations at both coarse-grained and atomistic resolutions, we identify three CL binding sites on the translocase, in agreement with those seen in crystal structures and inferred from nuclear magnetic resonance measurements. Characterization of the free energy landscape for lateral lipid interaction via potential of mean force calculations demonstrates the strength of interaction compared to those of binding sites on other mitochondrial membrane proteins, as well as their selectivity for CL over other phospholipids. Extending the analysis to other members of the family, yeast Aac2p and mouse uncoupling protein 2, suggests a degree of conservation. Simulation of large patches of a model mitochondrial membrane containing multiple copies of the translocase shows that CL interactions persist in the presence of protein-protein interactions and suggests CL may mediate interactions between translocases. This study provides a key example of how computational microscopy may be used to shed light on regulatory lipid-protein interactions.

The mitochondrial adenine nucleotide translocase (ANT, also known as the ADP/ATP carrier, AAC) facilitates export of ATP outward across the inner membrane into the cytoplasm, in exchange for import of ADP back into the matrix.[1] This is driven by the electrochemical gradient across the inner membrane. The translocase is one of the most abundant proteins in the mitochondrial inner membrane.[2] Its abundance and the importance of this transporter have rendered it one of the best characterized members of the mitochondrial carrier (MC) family, which facilitate the movement of a range of metabolites in and out of mitochondria.[3]

The translocase has an intimate relationship with a key lipid in mitochondrial membranes, cardiolipin (CL). The presence of CL is essential for maximal stability of the translocase in vitro, and CL has been suggested to modulate interactions of ANT with other proteins.[4,5] Knockdown of CL synthesis in humans results in Sengers syndrome, a condition in which the carrier is depleted from the inner membrane, leading to impaired oxidative phosphorylation and symptoms ranging in severity from mild difficulties with exercise to neo-natal fatality.[6,7]

To date, atomic-resolution structures of bovine and yeast forms of ANT have been determined in a cytoplasm-facing state.[8,9] A nuclear magnetic resonance (NMR) structure has been determined for another MC family member, uncoupling protein 2 (UCP2).[10] These translocase structures revealed a fold in which six transmembrane helices (H1–H6) form a bundle around a central substrate binding cavity, abutted on the matrix side by three amphipathic helices (MH1–MH3) arranged around the outer rim of the helix bundle (Figure ). There is an internal pseudo-3-fold symmetry, seen also in the UCP2 structure, which sequence comparisons suggest to be conserved across the MC family. The presence of CL was essential for reconstitution and crystallization of ANT, and CL molecules have been resolved in all crystal forms of the translocase, bound between the short matrix helices. This agrees with earlier 31P NMR measurements that indicated tight binding of CL molecules to bovine heart translocase.[11] Although CL molecules were not resolved in the structure of UCP2, the presence of this lipid was essential to obtain interpretable spectra.[10]

Figure 1

ANT1 and CL. (A) Topology and structure of bovine ANT1 colored to indicate the internal three-fold repeat symmetry. The structure shown corresponds to the cytoplasm-facing state of ANT1 (Protein Data Bank entry 1OKC).[8] (B) Coarse-grained (CG) representation of ANT1 embedded in a lipid bilayer. (C) Atomistic and CG models of CL.

Despite the importance of CL in the structure and function of the ADP/ATP translocase and other MC family proteins, characterization of these interactions within a membrane environment is relatively scarce. Multiscale molecular dynamics simulations have proven to be a powerful tool for identifying lipid binding sites[12] and for characterizing the energetics[13−15] and dynamics[16] of the interactions of the lipid with membrane proteins. In this approach, simulations are applied at both coarse-grained (CG) and atomistic resolutions.[17] In CG simulations, the granularity of the system is decreased (Figure B,C), thus easing the limitations on computationally accessible system sizes and time scales inherent in atomistic simulations. This method has been used to predict PIP2 lipid binding sites on Kir channels,[18] subsequently confirmed by crystal structures,[19,20] as well as to predict conserved lipid interaction patterns around all aquaporins of known structure.[21] Coarse-graining a system may also allow very large scale simulations of membrane patches containing multiple copies of a given protein. Such simulations have recently been used to define interactions of CL with the mitochondrial respiratory chain complexes III and IV,[22] to determine interactions of lipids with GPCRs,[23] and to explore bacterial outer membrane protein diffusion and oligomerization.[24]

Here we apply multiscale MD simulations to investigate the binding of CL to ANT and to UCP2 in model mitochondrial membranes. We assess structural, energetic, and dynamic aspects of lipid interactions with a single copy of each protein, before evaluating the interplay of protein–lipid and protein–protein interactions in large scale simulations of membrane patches containing 25 copies of the ANT molecule.

Methods

Coarse-Grained Simulations

All simulations were performed using the GROMACS 4.6 (www.gromacs.org) simulation package[25] and the MARTINI (http://md.chem.rug.nl/) force field.[26−28] A summary of the simulations performed is provided in Table .

Table 1

Overview of the Simulations Performed

description	membrane composition	duration
1× bovine ANT1,a CG	outer leaflet, PC (100%)	5 × 6 μs
	inner leaflet, PC (90%)/CL (10%)
1× bovine ANT1,a CG	outer leaflet, PC (50%)/PE (50%)	5 × 1 μs
	inner leaflet, PC (42.5%)/PE (42.5%)/CL (15%)
1× yeast Aac2p,b CG	outer leaflet, PC (50%)/PE (50%)	5 × 1 μs
	inner leaflet, PC (42.5%)/PE (42.5%)/CL (15%)
1× UCP2,c CG	outer leaflet, PC (50%)/PE (50%)	5 × 1 μs
	inner leaflet, PC (42.5%)/PE (42.5%)/CL (15%)
PMF calculations, 1× ANT1,a CG	PC (100%)d	32 × 1–2 μs per PMF
25× bovine ANT1,a CG	outer leaflet, PC (60%)/PE (40%)	20 μs
	inner leaflet, PC (50%)/PE (40%)/CL (10%)
25× bovine ANT1,a CG	outer leaflet, PC (60%)/PE (40%)	20 μs
	inner leaflet, PC (60%)/PE (40%)
1× bovine ANT1,a AT	outer leaflet, PC (50%)/PE (50%)	0.5 μs
	inner leaflet, PC (42.5%)/PE (42.5%)/CL (15%)

Simulation based on the crystal structure of bovine ANT1 (PDB entry 1OKC(8)).

Simulation based on the crystal structure of yeast Aac2p (PDB entry 4C9G(29)).

Simulation based on the crystal structure of UCP2 (PDB entry 2LCK(10)).

See Potential of Mean Force Calculations and Methods for further details regarding membrane composition.

Protein Data Bank (PDB) coordinate sets were processed using MemProtMD[30] to remove cocrystallized lipids, waters, and carboxyatractyloside inhibitors. CG simulations used the MARTINI version 2.1 force field.[27] Simulations of ANT utilized an ElNeDyn elastic network.[31] Simulations of UCP2 used an elastic network with a distance cutoff of 1 nm and a force constant of 1000 kJ mol–1 nm–2. In addition, long distance PRE restraints (force constant of 1000 kJ mol–1 nm–2) were applied to residues 68, 105, 202, and 205. The residues were chosen on the basis of the NMR restraint file (PDB entry 2LCK) and were required to prevent collapse of the UCP2 structure, as was also observed by Zoonens et al.[32] Each protein was placed in a box with dimensions of 12 nm × 12 nm × 10 nm, containing 460 randomly oriented PC lipids. The system was solvated using standard MARTINI water particles and neutralized with ∼0.15 M NaCl. A 50 ns self-assembly simulation[33] was performed to allow formation of a lipid bilayer around the protein. PC molecules within the bilayer were subsequently exchanged with other lipid species to form an asymmetric bilayer of specified lipid composition as described previously.[34] For repeat simulations, the lipids were independently exchanged such that each repeat contained a different initial distribution of lipid molecules consistent with the same specified lipid composition. CL parameters were from ref (35). CL was modeled with a net ionization state of −2e (i.e., −1e per phosphate group) based on experimental estimates.[36] The acyl tails of PC and PE lipids were modeled as 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE), while CL was modeled as 1,3-bis(1-palmitoyl-2-oleoyl-sn-3-phosphatidyl)-sn-glycerol.

Extended systems containing 25 copies of the ANT1 molecule were generated using the GROMACS genconf tool to concatenate a system containing a single ANT1 molecule embedded in a PC bilayer onto a 5 × 5 grid (Figure S1). The lipids of the energy-minimized concatenated PC membrane were then exchanged for a model mitochondrial membrane as described above.[34] The resultant system dimensions were 65 nm × 65 nm × 14 nm, containing 25 proteins and ∼11000 lipids, solvated by standard MARTINI water and neutralized by ∼0.15 M NaCl. The systems contained ∼4.7 × 105 particles. An initial 1 μs equilibration period was used to facilitate equilibration of lipid–protein interactions and facilitate randomization of protein orientations. During this period, the backbone particle of a single proline residue in each ANT1 molecule (Pro27) was restrained in the X–Y plane, thus allowing rotation and randomization of protein orientations, but preventing their translational movement. Systems both with and without CL were then run for 20 μs of an unrestrained production run simulation (see Table ).

Simulations were performed under NPT conditions at a temperature of 310 K and a pressure of 1 bar with periodic boundary conditions. The temperature and semi-isotropic pressure were controlled using the Berendsen thermostat and Berendsen barostat[37] with coupling constants of 4 ps. Equations of motion were integrated using the leapfrog algorithm with a time step of 20 fs. Electrostatics were modeled as reaction field with a Coulomb cutoff of 1.2 nm, using a potential shift modifier. van der Waals interactions were smoothly shifted off between 0.9 and 1.2 nm. The LINCS algorithm[38] was employed to constrain covalent bonds to their equilibrium values.

Potential of Mean Force Calculations

PMFs were calculated via umbrella sampling, utilizing a protocol previously described.[14] A 2 μs production run was performed with bovine ANT1 (PDB entry 1OKC)[8] embedded in a PC bilayer containing 10% CL within the inner leaflet. For each of the three sites, a snapshot containing a bound CL/PC molecule was extracted. A new PC bilayer was subsequently assembled around the position-restrained (Fc = 1000 kJ mol–1 nm–2) protein and lipid of interest during a 50 ns self-assembly simulation.[33] This simple composition was chosen to accelerate convergence of calculations and has previously been shown to cause no significant deviation compared to that for protein–lipid PMFs conducted in two-component PC/PE bilayers.[22] A one-dimensional (1D) reaction coordinate ranging from the lipid-bound state to the lipid-unbound state was then generated via a steered MD (SMD) simulation. Each CL/PC/PS molecule was pulled away from the protein over a distance of 3–4 nm at a rate of 0.1 nm/ns (Fc = 1000 kJ mol–1 nm–2). The force was applied to the GL0 headgroup particle of CL, and the PO4 headgroup particle of PC/PS. Weak positional restraints (Fc = 400 kJ mol–1 nm–2) were also applied in the X–Y plane to the backbone particle of each proline residue within each of the three Px[D/E]xx[K/R] motifs of the protein. Application of such restraints acted to prevent rotation of the protein and translational “following” of the tight binding CL molecules as they were pulled away. In addition, weaker positional restraints (Fc = 50 kJ mol–1 nm–2) were applied to the GL0/PO4 particle of the CL/PC lipid in the Y direction. This served to retain the lipid on the 1D reaction coordinate and thus enhance sampling of the chosen pathway. Snapshots spaced by ∼0.1 nm were extracted from the SMD simulation and used as input for umbrella sampling simulations. Each window of the umbrella sampling was run for 1–2 μs, with umbrella biasing potentials (Fc = 1000 kJ mol–1 nm–2) applied between the center of mass of the three restrained proline residues and the GL0/PO4 particle of the CL/PC/PS lipid. Within umbrella sampling simulations, long-range electrostatics were modeled using particle mesh Ewald (PME).[39] The subject lipid was treated separately from bulk lipids for temperature and pressure coupling. Approximately 32 umbrella sampling simulations were run per PMF. PMF profiles were constructed using the GROMACS implementation (g_wham) of the weighted histogram analysis method (WHAM).[40] Bayesian bootstrapping with 200 bootstraps was used to estimate the errors for each profile. Convergence was assessed by comparing profiles calculated from independent 100 ns segments of simulation time (Figure S2), as well as by calculating profiles via replica exchange initiated from different starting configurations. We previously tested the reproducibility of this protocol and found well depths initiated from independent repeats spanned a range of 5 kJ/mol (i.e., 2RT).[14] We tested this for the system presented here by conducting two independent repeats of the entire calculation for site I, along with two repeats utilizing different lateral lipid restraints (Figure S3). The well depths agreed to within 5 kJ/mol, suggesting this range is applicable to the current system and indicating the robustness of the protocol to the choice of initial snapshot and lipid restraint.

Atomistic Simulations

Atomistic simulations employed the GROMOS96 53A6 force field.[41] A final snapshot of a 1 μs CG simulation was converted to atomistic detail using a fragment-based protocol.[42] The system was solvated using the SPC water model and neutralized with ∼0.15 M NaCl. A 1 ns equilibration simulation was performed with protein Cα particles restrained (Fc = 1000 kJ mol–1 nm–2). This was followed by a 0.5 μs unrestrained simulation (see Table ). The temperature was controlled at 310 K using a V-rescale thermostat[43] with a coupling constant τt of 0.1 ps. Pressure was maintained at 1 bar using a Parrinello–Rahman barostat[44] with a compressibility of 4.5 × 10–5 bar–1 and a coupling constant τp of 1 ps. Electrostatic forces were modeled using particle mesh Ewald (PME).[39] van der Waals interactions were cut off at 1.2 nm. A time step of 2 fs was employed, and covalent bond lengths were constrained using the LINCS algorithm.[38]

Analysis and Visualization

Two-dimensional (2D) density maps for lipids and proteins were calculated as described previously.[24] Three-dimensional (3D) lipid density isosurfaces were computed using the VolMap plugin of VMD.[45] Diffusion coefficients were extracted from the mean square displacement (MSD) using g_msd and a time period of 20 ns. Contact and protein clustering analysis was performed using in-house scripts and MDAnalysis.[46] Lipid contacts were calculated between each residue of the protein and the headgroup particle of CL (defined by GL0) using a cutoff of 0.6 nm to define a contact. A cutoff of 4.4 nm between protein centers of mass was used for protein clustering analysis, while per residue protein–protein contacts were calculated utilizing the approach described in ref (24). The VMD[45] software package was used for visualization of simulation data.

Results

Identification of CL Binding Sites on Bovine ANT1

A CG molecular model of bovine ANT1 was embedded into a simple model membrane containing PC (90%) and CL (10%) in the matrix-facing leaflet, and PC (100%) in the cytoplasm-facing leaflet. This CL content lies within the range estimated for mitochondrial inner membranes from lipidomics.[47−49] Five repeat simulations of bovine ANT1 were conducted, each with a duration of 6 μs and beginning from a different distribution of lipids within the bilayer (see Table and Methods for details).

CL molecules were seen to encounter and bind to three specific sites on the protein [sites I–III (Figure and Movie S1)]. Each site was composed of a groove on the protein surface formed by the N-terminal end of the even-numbered transmembrane helices (H2, H4, and H6) of each repeat, and the short amphipathic matrix helices (MH1–MH3). The location of the sites thus has a symmetry related to the internal three-fold structural repeat. The location of these sites is in excellent agreement with those identified in X-ray crystal structures.[8,9] This indicates that the sites occupied in the crystal structure persist within a membrane environment. Interactions were predominantly mediated by the headgroup moieties of CL, as seen in the time-averaged 2D density maps for CL (Figure B). In contrast, the acyl tails exhibited more dynamic interactions, which resulted in less well-defined density, though a small degree of ordering was seen at site III. This is consistent with a number of crystal structures in which only the headgroup moieties and adjacent atoms of the acyl chains were resolved.[29] The dominance of headgroup interactions seen here is also consistent with functional assays performed in yeast showing Aac2p assembles normally regardless of the acyl chain composition of CL.[50]

Figure 2

CL binding sites on bovine ANT1. (A) In the left panel, the initial simulated system consisted of a single molecule of bovine ANT1 embedded in a mixed PC (lime) and CL (magenta) bilayer, solvated with standard MARTINI water particles (transparent surface) and neutralized with ∼0.15 M NaCl (blue spheres). The right panel is the final snapshot of a 6 μs simulation showing three bound CL molecules (magenta). (B) Time-averaged 2D density maps of CL in the membrane plane for the headgroup (left) and acyl tail moieties (right). (C) View onto the base of the translocase from the matrix showing the positioning of the three CL (magenta) binding sites around the protein within the X-ray structure (left) and CG simulations (right). The CL binding sites within the CG simulations are shown by the time-averaged probability density surface for the CL molecules, contoured to show the three binding sites. The density surface was calculated from 5 × 6 μs simulation repeats, each starting from a different random lipid distribution. See Figure S4 for per simulation density maps.

The number of contacts formed between each residue of the protein and the CL headgroup was calculated over the course of each simulation. Within each site, the CL headgroup showed a preference for interaction with the N-termini of the even-numbered TM helices and matrix helices. In site I, the highest degree of contact was seen with S68, G72, T154, and G155. Residues showing high levels of contact in site II were Q174, G175, G252, and T253, while in site III, they were G52, I53, G272, and W274 (Figure A). A number of the interacting residues form part of the highly conserved [YWF][RK]G and [YF]xG motifs,[2] respectively, located at the N-termini of the even-numbered TM helices and matrix helices (Figure B). In particular, the glycine residues within these motifs were seen to form frequent contacts with CL in the simulations. The colocalization of CL interactions to these conserved sequence motifs is encouraging, given also the suggested role of the glycine residues[2] within the motifs in mediating CL interactions within crystal structures.[9] PC lipids showed comparatively less specificity for these regions (Figure S5). To assess the robustness of the results to the presence of PE lipids (present in vivo(47)) and changes in CL concentration, 5 × 1 μs simulations were conducted in mixed membranes with a composition of PC (50%) and PE (50%) in the cytoplasm-facing leaflet and PC (42.5%), PE (42.5%), and CL (15%) in the matrix-facing leaflet. The simulations yielded very similar patterns of contact of CL with the PC/CL only membranes, suggesting the interaction is not sensitive to the presence of PE (Figure C), as has been seen for the interactions of CL with other mitochondrial proteins.[22]

Figure 3

Interactions of CL with bovine ANT1. (A) Normalized average frequency of CL headgroup contacts per residue. Residues within each of the even-numbered helices (H2, H4, and H6) and the matrix helices (MH1–MH3) exhibiting the highest degree of contact are labeled (see Methods for full details of contact analysis). (B) Lateral view into each CL binding site with the backbone trace of the protein colored by the relative extent of contacts formed with CL, from white (no contact) to magenta (high contact). Insets show the corresponding view onto the whole translocase colored according to a three-fold repeat topology. The positions of the conserved motifs are labeled. (C) Normalized average frequency of phospholipid headgroup contacts per residue, for the lipid species within three-component membranes.

To further assess interactions within the identified binding sites, a CG snapshot of ANT1 was converted[42] to atomistic detail and simulated for 0.5 μs. Sites I and II of the chosen snapshot were occupied by CL molecules, while site III was unoccupied, thus allowing us both to test the stability of bound CL molecules and to explore possible atomistic binding events within the same extended simulation. Calculation of the Cα root-mean-square deviation (RMSD) during the simulation revealed a stable value of ∼0.3 nm for the core structural elements (i.e., helices) and ∼0.4 nm for all residues, indicating that, despite the absence of the carboxyatractyloside inhibitor, the apo-cytoplasm-facing structure is conformationally stable on the time scale simulated (Figure S6). The CL molecules within sites I and II remained bound (Figure ), with the phosphate moieties contacting those conserved motifs seen to form frequent contacts in the CG simulations. In contrast, the acyl chain tails exhibited more dynamic interactions and adopted a range of conformations (Figure S6,B). Of interest within site I, the side chain of R71 remained pointing downward, consistent with its position within crystal structures (Figure B).[9] No lipid binding events were observed at site III, despite one CL being positioned in the vicinity of the site. This suggests longer time scales may be required to capture the full process of lipid binding at atomistic resolution.

Figure 4

Atomistic simulation of bovine ANT1. (A) The left panel is a view from the matrix side, with the phosphate particles of the bound CL molecules shown as spheres. For each phosphate particle, 1000 evenly distributed snapshots are shown over the 0.5 μs simulation, colored according to simulation time. This may be compared to the phosphate distribution (right) for a noninteracting bulk CL. (B) Final snapshot of site I showing the arrangement of conserved motifs (gray), N-terminal ends of the helices (green and cyan), and CL (magenta).

Energetics of CL Binding

To assess the strength and selectivity of CL interactions, we calculated a potential of mean force (PMF) for the CL–ANT1 interaction via CG simulation and umbrella sampling. The PMF (or free energy profile) between two species describes the change in free energy along a particular reaction coordinate and is calculated from the probability distribution of the two species along that coordinate.[51] This approach has recently been used to define the binding of CL to mitochondrial respiratory chain complexes III and IV[22,52] and to characterize lipid interactions of the epidermal growth factor receptor (EGFR) transmembrane domain.[14] For each of the three sites, a steered MD simulation was conducted in which a force was applied to pull the CL molecule out of its binding site into the bulk membrane along a 1D reaction coordinate (r) perpendicular to the protein surface in the plane of the bilayer (Figure A). The free energy profile along this coordinate was then calculated via umbrella sampling, with r defined as the distance between the center of mass of three conserved prolines within ANT1 and the β-glycerol moiety of CL.

Figure 5

Potentials of mean force for binding of the lipid to ANT1. (A) View onto the base of ANT1 illustrating the reaction coordinates used in free energy calculations. ANT1 is colored according to its three-fold repeat topology, with the position of each binding site indicated by the time-averaged probability density isosurfaces for CL phosphates (magenta). The arrows indicate the approximate 1D reaction coordinates along which the PMF was calculated. (B) PMF profiles for interactions of CL (magenta) and PC (black) with each of the three binding sites, and a control non-CL binding region. For site I, a PMF profile was also calculated for a further anionic lipid species, PS (green). The standard deviations estimated from bootstrapping[40] are shown as the shaded area behind each curve.

The profiles revealed comparable tight binding of CL molecules to each of the three sites, with well depths of approximately −20 kJ/mol (Figure B). In contrast, pulling a CL molecule away from a region not predicted to bind CL revealed a well depth of approximately −4 kJ/mol, indicating no significant interaction in this region. Conducting the same calculations for zwitterionic PC lipids showed no significant interaction of this species beyond the value of RT (2.5 kJ/mol), demonstrating the marked selectivity of these sites for CL over zwitterionic phospholipids. These values may be compared to those obtained for cytochrome c oxidase (CIV) via the same method,[15] revealing that the binding sites on bovine ANT1 have an affinity for CL higher than those of most of the sites on CIV, other than for two high-affinity sites on CIV that are estimated to bind with an interaction free energy of approximately −32 kJ/mol. Interestingly, calculation of the PMF profile for a second anionic lipid species, phosphatidylserine (PS), at site I revealed a well depth of approximately −15 kJ/mol, i.e., intermediate between those of CL and PC. Although the abundance of PS[47] within the inner mitochondrial membrane is comparatively low, these calculations suggest the capacity of the binding sites to interact with other anionic lipid species, albeit less favorably than with CL.

Conservation of Binding Sites

The strong binding of CL to bovine ANT1 together with the structural conservation[53] of the internal three-fold repeat of mitochondrial carriers suggests that these sites are likely to be conserved over other members of the family. We tested this by running CG simulations of yeast Aac2p and mouse UCP2. Multiple (5 × 1 μs) simulations were run for Aac2p and UCP2 in the same manner as they were for ANT1 (Table ). Analysis of the average probability distribution of CL over the course of the simulations showed preferential localization to the three equivalent sites seen for ANT1 (Figure ). These positions compare well with the multiple crystal structures of bovine ANT1[8,9] and yeast Aac2p and Aac3p.[29] The simulations are supportive of direct interactions of CL with UCP2 and predict that these interactions occur at sites equivalent to those seen in ANT1 and Aac2p.

Figure 6

Conservation of CL binding sites. (A) Location of cocrystallized lipids (magenta) in bovine ANT1 (PDB entry 1OKC)[8] and yeast Aac2p (PDB entry 4C9G).[29] (B) CL binding sites identified by CG–MD simulations. The magenta surface shows the probability density for the CL headgroup phosphate particles around a snapshot of each protein, contoured to indicate the positions of the phosphate moieties of CL within each site. (C) Normalized average number of CL headgroup contacts per residue within simulations of bovine ANT1, yeast Aac2p, and mouse UCP2.

Analysis of lipid contacts formed in the yeast Aac2p and UCP2 simulations showed patterns spatially similar to those of bovine ANT1 (Figure C), despite differences in sequence and the somewhat “looser” packing of the transmembrane helices in the UCP2 structure. This suggests the groovelike architecture of the sites is also an important contributor to CL binding, and that a degree of sequence substitution may be tolerated. The conservation of these sites within both crystals (obtained in the presence of detergent) and a simulated membrane environment along with their presence in three MC proteins is supportive of both their biological significance and their likely conservation in other members of the family.

Lipid Interactions and Organization of Bovine ANT1 in Large Membrane Patches

Large scale simulations of membrane patches containing multiple protein copies have recently proven to be useful in understanding the dynamic organization of lipids and proteins.[54] In light of this, a single molecule of bovine ANT1 embedded in a small membrane patch was replicated (Figure S1) to create a large membrane patch with 25 ANT1 molecules in an ∼65 nm × 65 nm area of membrane containing ∼11000 lipid molecules (see Table ), with a cytoplasm-facing leaflet containing PC (60%)/PE (40%) and a matrix leaflet containing PC (50%), PE (40%), and CL (10%). This composition lies within the range estimated for the inner mitochondrial membrane from lipidomics.[47,55] This membrane was used as the basis of 20 μs of unrestrained molecular dynamics.

Examination of the final snapshot showed a degree of protein oligomerization (Figure A), with a mixture of monomers, dimers, and higher-order oligomers observed. Analysis of the dynamics of protein clustering showed the self-assembly of oligomers began on a submicrosecond time scale (Figure B and Movie 2), with ∼50% of the proteins participating in an oligomeric state by 8 μs. Calculation of the averaged 2D protein density maps for ANT1 and CL showed that this interaction occurred at three specific interaction interfaces, which colocalized with the three CL binding sites identified previously (Figure C). Indeed, visual inspection of the simulations showed that CL molecules were often found at these interfaces, loosely bridging proteins. However, in contrast to the crystallographic dimers,[29] we did not observe a stable arrangement in which two CL molecules were present between the two proteins. Such an arrangement may be unfavorable because of the proximity of the negatively charged (−2e) CL headgroups, as well as a degree of steric conflict. This is suggestive of possible competition between protein–protein and protein–lipid interactions. Calculation of 2D density maps and ANT-CL contacts for the first 1 μs of the simulation (during which period the ANT molecules are predominantly monomeric) and the final 1 μs (during which period the ANT molecules are predominantly oligomeric) revealed similar profiles, indicating the broad patterns of CL interaction remain similar under the two regimes (Figure S7). However, further simulations with improved sampling and resolution would be needed to unpack any finer details of CL behavior within oligomers.

Figure 7

Organization of bovine ANT1 in large membrane patches. (A) Final snapshot (t = 20 μs) of the matrix leaflet showing the organization of 25 ANT1 proteins in a PC/PE/CL mixed membrane. (B) Clustering dynamics shown as the percentage of ANT1 cluster size over the simulated time course, along with a snapshot of an ANT1 dimer with CL mediating the interaction. Other lipid molecules have been omitted for the sake of clarity. (C) Density plots showing the frequency of occurrence of ANT around ANT, and CL around ANT. The protein has been omitted from the CL plot for the sake of clarity.

The interfaces at sites I and II showed a stronger preference for protein–protein interaction, with comparatively weaker density seen at site III. However, we did not observe a clear preference for certain permutations of protein–protein interaction interfaces, but rather a mixture of interfaces were seen. For example, the interface at site I was able to interact with all other interfaces. For all regions, the contact area between interacting proteins was relatively small, with contacts to neighboring protein molecules seen almost exclusively at the cytoplasm-facing ends of the TM helices, and the short matrix helices (Figure S8). Residues exposed to the core of the membrane on the concave surface of the translocase formed comparatively few direct protein–protein contacts, in keeping with the “hourglass” shape of ANT1. Contacts of this hydrophobic surface were instead mediated by lipid tails.

Interestingly, simulating a comparably sized membrane patch without CL in the bilayer (Table ) resulted in a greater degree of oligomerization of the ANT population over the 20 μs time course (Figure S9). This suggests that CL may have the capacity to modulate the oligomerization process. Calculation of diffusion coefficients for the ANT1 protein molecules showed values comparable to those seen for GPCRs[23] and did not exhibit a significant difference in the presence and absence of CL (Figure S9). This suggests the difference in oligomerization is unlikely to be due to differences in diffusive behavior in the two lipid mixtures and may instead be indicative of modulation via direct interaction.

Discussion

Cardiolipin molecules have been observed bound to several crystal forms of the bovine ANT1.[8,9] Three CL molecules are bound per ANT1 monomer, with the phosphate groups positioned within grooves formed by the short matrix helices. The sites identified within our simulations are in complete agreement with the experimentally determined sites, demonstrating that the latter persist within a membrane environment and thus are likely to be present in vivo. Calculation of free energy profiles for lateral lipid interaction with these sites shows moderately tight binding of CL molecules compared to other regions on the protein, along with clear selectivity for CL over zwitterionic phospholipids, in agreement with indications from 31P NMR measurements.[11] The well depths for these profiles are −20 to −25 kJ/mol and may be compared to profiles calculated via the same method for other protein–lipid interactions, including interactions of glycolipid and phosphoinositide with the EGFR,[14] binding of phosphoinositide to PH domains,[13] and in particular interactions of CL with mitochondrial CIII and CIV.[22,52] This overall suggests the sites on bovine ANT1 have a clear affinity for CL, which is greater than that of the majority of CL binding sites on CIV, and is superseded by only two particularly high-affinity sites on CIV that exhibit interaction free energies of approximately −32 kJ/mol.[15]

Extending the simulations to the two other members of the mitochondrial carrier family of known structure, yeast Aac2p and mouse UCP2, revealed comparable direct interaction of CL with both proteins and demonstrated conservation of binding site location. For yeast Aac2p, this agrees well with the location of cocrystallized CL molecules.[29] At present, there are no structural data for binding of CL to UCP2. However, the presence of CL was required to obtain usable NMR spectra, suggestive of direct interaction.[10] Our simulations support this and predict CL to bind directly to three sites on UCP2 equivalent to those seen in ANT1 and Aac2p. Interestingly, recent biochemical work[56] on the related protein, UCP1, suggested a binding stoichiometry of three CL molecules per monomer, in support of our simulation-based prediction for UCP2. The conservation in binding site location over several proteins suggested by simulation, structural and biochemical studies, and the predicted conservation of the internal three-fold repeat is strongly indicative of likely conservation of the CL binding sites in other members of the mitochondrial carrier family. It will be interesting to see how far this conservation extends as further structures emerge.

Advances in simulation methodology[12,57] allow access to extended time and length scales. We were thus able to investigate the larger (∼0.1 μm) scale organization and protein–lipid interactions of bovine ANT1. These simulations suggest the cytoplasm-facing conformation of bovine ANT1 has the capacity to form self-interactions at specific interfaces and thus to form oligomers within model membranes. The protein–protein interfaces colocalize with the CL binding sites, and CL is often seen to bridge between adjacent proteins. However, in contrast to crystallographic dimers,[9] we did not observe a stable arrangement in which two CL molecules may be present at an interface, with the interfaces rather being mediated by a single CL. This may be related to the unfavorable electrostatics of having two negatively charged (−2e each[58]) CL headgroups in the proximity, in a membrane environment that allows for protein–protein and protein–lipid assemblies more dynamic and transient than what may be observed in a crystal. It would therefore be attractive to analyze the relative free energies of the different modes of protein–protein and protein–CL–protein interactions observed within these dynamic complexes. However, this would be likely to require the development of more effective simulation-based sampling regimes before such free energies may be reliably estimated.[59] Interestingly, simulations of membrane patches in the absence of CL showed a greater degree of oligomerization over the 20 μs simulation compared to the membrane with CL present. This suggests a more complex relationship between CL and ANT–ANT interactions, with a single CL molecule bound between interfaces that can mediate protein interactions as a molecular glue,[60] but the presence of two CL molecules between interfaces possibly acting as a block to disfavor oligomerization.

The oligomeric state of ANT1 has been much debated.[4,61,62] While the functional unit of the translocase is widely accepted to be monomeric, this does not preclude oligomerization under certain conditions. It is important to remember that our observations pertain to the cytoplasm-facing state of the translocase within simple model membranes: both conformational changes and biologically complex lipid mixtures may also influence oligomerization, as observed for GPCRs.[23,63] In particular, other anionic lipids, e.g., PS, PA, and PI, at low levels in the inner mitochondrial membrane[47−49] may interact with ANT within living cells. Other in vivo factors such as the high degree of membrane curvature of the mitochondrial inner membrane, along with a more heterogeneous protein population, may also influence ANT–ANT interactions. We note that a number of in vitro model membrane studies have suggested CL may regulate the activity of certain proteins via modulation of curvature and other biophysical properties of the inner membrane.[64−67] Larger scale simulations capable of assessing such macroscopic properties would be of interest in this regard. Likewise, we note that key inner membrane-resident compounds such as electron-transferring quinones have been suggested to influence membrane properties.[68] Structural data for a number of other mitochondrial membrane proteins are now available, along with EM data on their in vivo arrangement and membrane curvature.[69] As the number of lipid parameters and time and length scales accessible to simulations continue to increase,[54] it may be possible to address these factors in large scale simulations. With regard to mechanistic insights into the influence of CL on the ANT transport cycle, as information about the possible nature of ANT conformational states continues to emerge,[2,70] a key extension of the current simulations would be to examine the differential interactions of CL with multiple conformational states of ANT to reveal the detailed mechanism of the influence of CL on ANT function.

Methodological Considerations

It is important to be aware of possible limitations in the simulation methods employed. Use of CG models results in an inherent loss of chemical detail, which may be expected to influence the accuracy of protein–lipid interactions. Although the ionization state of the CL headgroup has been estimated to be −2e in solution,[36] it is possible this may change depending on the microenvironment, as may protein ionization states. Capturing such effects would require constant-pH simulations, with a concomitant increase in computational load.[71] Despite the absence of such effects, good agreement is seen between our CG and atomistic simulations, together with previously resolved crystal structures. In addition, we note the method has previously been shown to correctly predict lipid binding sites on an increasing number of other membrane proteins,[19,21,52] in agreement with experiment,[12] and is capable of predicting the effects of mutation on lipid binding.[13,14]

A further key limitation of current simulations is the accessible time scale, which remains a challenge even when using CG models. This is particularly pertinent for large systems containing multiple copies of proteins. These limitations may impede the investigation of slowly converging system properties such as protein–protein interactions.[22,23] In the current study of a 20 μs time scale, we capture only part of the self-assembly process, and an equilibrated state for protein–protein interactions has likely not been reached. Thus, insights into the strength of protein–protein interactions are circumscribed because of limited statistics for association and dissociation events. Nonetheless, important information can be gleaned in identifying protein interaction interfaces and the influence of lipids on the process, as seen for, e.g., respiratory complexes III and IV,[22] and for the S1P1 receptor.[23]

With regard to protein–protein interactions, the MARTINI model has been suggested to overestimate interaction strengths in aqueous solutions.[72] However, within a membrane environment, the model performs well, with the free energy profiles for glycophorin A transmembrane helix dimerization estimated to be approximately −30 to −40 kJ/mol in MARTINI CG,[73,74] and −45 to −60 kJ/mol in all-atom simulations,[75,76] which agree reasonably with available experimental data.[77] The application of MARTINI dihedral restraints together with an ElNeDyn elastic network[31] prohibits conformational transitions, and thus, we consider only the cytoplasm-facing state of the translocase, for which crystal structures are known. One might anticipate that potential conformational transitions[70] within the protein could influence lateral protein–protein interactions as well as protein–CL interactions.

Conclusions

In summary, we have characterized three CL binding sites on the cytoplasm-facing conformation of bovine ANT1 within model mitochondrial membranes. Free energy calculations reveal the ∼100 μM affinity of these sites for CL, and their selectivity over zwitterionic phospholipids. Equivalent sites are present for mouse UCP2 and yeast Aac2p, indicative of both their biological significance and conservation within the transporter family, as further supported by recent biochemical studies of UCP1.[56] Simulations in extended membrane patches demonstrate that CL interactions at these sites persist in the presence of protein–protein interactions, and indeed, CL molecules may mediate interactions between adjacent ANT1 molecules, permitting formation of a range of oligomeric states. More generally, our simulations demonstrate further the ability of computer simulations to identify lipid binding sites on membrane protein, providing near-atomic-resolution characterization of protein–lipid interactions to complement biochemical and structural studies.