Computational investigation of cholesterol binding sites on mitochondrial VDAC.

The mitochondrial voltage-dependent anion channel (VDAC) allows passage of ions and metabolites across the mitochondrial outer membrane. Cholesterol binds mammalian VDAC, and we investigated the effects of binding to human VDAC1 with atomistic molecular dynamics simulations that totaled 1.4 μs. We docked cholesterol to specific sites on VDAC that were previously identified with NMR, and we tested the reliability of multiple docking results in each site with simulations. The most favorable binding modes were used to build a VDAC model with cholesterol occupying five unique sites, and during multiple 100 ns simulations, cholesterol stably and reproducibly remained bound to the protein. For comparison, VDAC was simulated in systems with identical components but with cholesterol initially unbound. The dynamics of loops that connect adjacent β-strands were most affected by bound cholesterol, with the averaged root-mean-square fluctuation (RMSF) of multiple residues altered by 20-30%. Cholesterol binding also stabilized charged residues inside the channel and localized the surrounding electrostatic potentials. Despite this, ion diffusion through the channel was not significantly affected by bound cholesterol, as evidenced by multi-ion potential of mean force measurements. Although we observed modest effects of cholesterol on the open channel, our model will be particularly useful in experiments that investigate how cholesterol affects VDAC function under applied electrochemical forces and also how other ligands and proteins interact with the channel.

Introduction

Voltage-dependent anion channels (VDACs) are integral membrane proteins in the mitochondrial outer membrane. VDACs are related to ancient porin proteins commonly seen in bacteria, and the eukaryotic channels regulate ion and metabolite flux between the cytosol and the mitochondrial intermembrane space.[1] Crystallographic and NMR structures of mammalian VDAC revealed a β-barrel with 19 strands connected by loops and an N-terminal helix inside the channel.[2−5] At low membrane potentials, VDAC assumes an open conformation with modest anion selectivity, whereas applied voltages greater than 30–40 mV cause VDAC to gate into a low conducting, cation-selective state.[6]

The VDAC-mediated confluence of the cytosol and mitochondrial intermembrane space suggests low membrane potentials across the mitochondrial outer membrane and predominantly open channels in vivo.[7,8] However, in addition to voltage, VDAC function is modulated by proteins that include Bcl-xL,[9] tBid,[10] hexokinase,[11] and tubulin,[12] and the interactions of these proteins with VDAC critically regulate apoptosis and cellular respiration.[8,13,14] Lipids can also affect channel activity by altering the properties of the surrounding membrane (e.g., lateral pressure) or through direct interactions with the protein.[6,15] In addition to endogenous regulators, VDAC is a target of drugs that include erastin[16] and general anesthetics,[17−20] although the functional influence of these ligands is not entirely clear.

Here, we used atomistic molecular dynamics simulations to investigate VDAC interactions with cholesterol. Cholesterol is known to bind mammalian VDAC in vivo[21,22] and in vitro,[2] and specific cholesterol sites on the protein have been identified with NMR.[2] Functionally, cholesterol enhances the structural integrity of isolated VDAC and aids channel insertion into bilayers,[21,23,24] and cholesterol promotes uniformity of open channel conductances.[21,24] Cholesterol has also been suggested to affect VDAC interactions with other proteins.[25] Together, these effects suggest VDAC function might be modulated when mitochondrial outer membrane cholesterol content increases, such as in cancers.[25−28]

The majority of proteins that bind cholesterol are integral membrane proteins;[22] however, cholesterol interactions with β-barrels, which in mammals consist of VDAC and other mitochondrial proteins,[29,30] have not been extensively characterized. This is in contrast to extensive efforts investigating cholesterol binding to other ion channels that reside on the plasma membrane.[31−37] Therefore, this work represents an early step toward understanding how cholesterol interacts with VDAC, the effects of cholesterol binding on channel structure and dynamics, and the potential interplay associated with ligand and protein binding.

Methods

Docking

Four residues in the mouse crystal structure of VDAC1 (PDB code 3EMN),[3] which was solved at 2.3 Å resolution, were mutated in PyMOL[38] to humanize the protein (Asp55 → Thr55, Val129 → Met129, Ser160 → Ala160, Val227 → Ile227). We specifically used 3EMN because the conformation of the N-terminal helix, which differs among other VDAC1 structures,[2,4] was confirmed with solid-state NMR.[39] The N-terminal methionine,[40] water, and detergent were removed before loading the structure into AutoDockTools.[41,42] To the protein, hydrogens were added, nonpolar hydrogens were merged, and Kollman charges were added. Molecular coordinates for cholesterol were downloaded from the CHARMM small molecule library;[43,44] in AutoDockTools, Geisteiger charges were added, nonpolar hydrogens were merged, and six torsions were allowed (i.e., cholesterol was fully flexible). For docking, grid boxes targeted the specific cholesterol sites and biased the exterior, membrane side of the β-barrel. Residue side chains in the site of interest and projecting outside the barrel were flexible during docking runs with AutoDock Vina.[45] AutoDock was programmed to return six docking results with search exhaustiveness of 100. The highest scoring poses were generally chosen for simulations. The similarity of some docking poses within each site often led to sampling of multiple poses during preliminary simulations used to build the cholesterol-bound VDAC model (see below). For this reason, multiple docking calculations with different grid centers were performed on some sites to increase heterogeneity in the starting poses.

MD Simulation System Setup

After docking, PDB files of VDAC, the docked cholesterol molecules, and corresponding flexed residues were created. The protein pore was oriented along the z-axis, normal to the membrane, using the PPM server of Orientations of Proteins in the Membrane,[46] and the output models were loaded into the CHARMM-GUI Membrane Builder.[47−49] In all simulations, the N-terminal residue (Ala2) was oriented to place it at z < 0. Ala2 was acetylated,[40] and residues were protonated according to their standard states at pH 7.4. The water thickness (minimum height on the top and bottom of the system) was 15.0 Å.

For preliminary simulations used to examine different orientations of cholesterol binding to VDAC, and for production simulations of the final cholesterol-bound model, 5 cholesterols were docked to the protein, and 160 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) molecules (80 in each leaflet) and 11 randomly distributed cholesterols comprised the membrane. Systems without docked cholesterol contained 160 DOPC and 16 cholesterols in the membrane. As a reference, phosphatidylcholines are the most abundant mitochondrial outer membrane phospholipids at 43–50 mol percent, with cholesterol 8–11 mol percent.[50,51] A total of 23 K+ and 25 Cl– ions were randomly placed by the CHARMM-GUI to neutralize the system and provide a salt concentration of 0.15 M KCl. KCl is generally favored over NaCl in VDAC experiments to minimize the effects of ionic size and diffusivity on measurements. Systems were separately generated with the CHARMM-GUI server, hence system details (number of atoms, starting coordinates, etc.) differed between simulations. Systems were approximately 87 Å, 87 Å, and 74 Å in x, y, and z dimensions, respectively, with about 56 000 total atoms that included about 9100 TIP3P waters.

Simulation Details

Atomistic molecular dynamics simulations were run with NAMD v2.9.[52] The CHARMM36 model was used for protein[43,53,54] and phospholipid[55,56] parameters, whereas the modified (CHARMM) 36c model was used for cholesterol;[57] parameters for TIP3P waters[58] and ions[59] are well-established. All simulations used periodic boundary conditions and particle mesh Ewald (PME) electrostatics. Interactions between nonbonded atoms were cut off at 12 Å, and bonds involving hydrogen were constrained using the SHAKE/RATTLE algorithms. A Langevin thermostat and barostat were used to maintain a temperature and pressure of 303.15 K and 1 atm, respectively, and no surface tension was imposed. The simulation time step was 2 fs. Prior to production, 20 000 minimization steps, and a 0.725 ns equilibration protocol detailed elsewhere[48] were performed on each system to gradually release restraints on the protein. In total, 22 separate simulations totaling over 1.4 μs contributed to this work, including 10 simulations of 100 ns.

Trajectory Analyses

Production simulations were analyzed with Visual Molecular Dynamics (VMD).[60] Every 20 ps snapshot was used for all quantitative analyses. Generally, the trajectories were aligned to and centered around the backbone of the crystal structure before performing the following analyses: RMSF: Root-mean-square fluctuation was measured for residue α carbons and describes the fluctuation of the atom relative to its average position throughout each simulation. Ion diffusion: Cl–/K+ permeability ratios were determined after normalizing the number of ions that diffused through the pore to the total number of that type of ion in the system (23 for K+ and 25 for Cl–). Diffusion was defined as traversing the membrane through the channel pore from −20 Å to 20 Å, or vice versa, along z. The protein spanned these coordinates in every frame of every simulation, as detected with the HOLE software.[61]Channel radius: The minimum channel radius for each frame was measured with the HOLE software, which fit consecutive spheres inside the channel with the centers spaced every 0.5 Å along the z axis. Hydrogen bonds: Hydrogen bonds were calculated with the hydrogen bond plugin of VMD with a 3.3 Å and 20° donor–acceptor cutoff. Electrostatics: The protein, membrane, and water contribution to the electrostatic potential was generated with the PME electrostatics plugin (PMEPot) using 88 × 88 × 80 grid position counts (∼1 Å spacing).[62] The generated maps contained an identical number of data points from each simulation (the dimensions of the systems were also essentially identical). Potential of mean force: Averaged multi-ion potential of mean force plots were made for each trajectory. For each frame, we measured both the number of ions and the number of water molecules in a cylinder of radius 15 Å that ran parallel to the z-axis through the center of the channel. This was used to calculate the concentration of ions, C(z), in 1 Å bins, and the relative free energy was estimated with the equation:where R is the gas constant, T is the temperature, and Cbulk is the concentration of ions in the bulk water.[63,64]Ionic density: Ionic density maps for entire trajectories were generated with the VolMap tool plugin of VMD.

Figure Preparation and Statistics

Structural figures were generated with VMD and PyMOL, and graphical figures were generated with GraphPad Prism (GraphPad Software Inc., La Jolla, CA) and KaleidaGraph (Synergy Software, Reading, PA); movies were also generated with VMD. Statistical analyses were performed as described in the text within the GraphPad Prism software. For each analysis, the mean from each simulation was considered a single experimental unit; therefore, where applicable, values are reported as the mean measurement from multiple simulations with standard error. Statistical comparisons were considered significantly different when p < 0.05.

Results and Discussion

Cholesterol-Bound VDAC

Published NMR chemical shift mapping experiments revealed resonances of nine VDAC backbone amides that had a significant chemical shift induced by the presence of cholesterol, with the structure of the protein essentially unchanged (Figure S1).[2] This indicated a change in the chemical environment surrounding the backbone amides and suggested specific binding by cholesterol. In that experiment, cholesterol was dissolved in a hydrophobic micellar phase and interacted with VDAC from outside the barrel, analogous to membrane cholesterol approaching the protein.

With binding demonstrated experimentally, we aimed to identify the effects of cholesterol site occupancy on VDAC. To initiate this, we began a combined docking and molecular dynamics approach targeting cholesterol to the specific VDAC backbone atoms implicated in binding. The topology of the protein and initial docking calculations suggested a model in which five cholesterol molecules with unique orientations were required to simultaneously occupy all residues that bind cholesterol. A 5 to 1 mol ratio of cholesterol to protein also reflected conditions under which chemical shift mapping was performed, and is consistent with the number of cholesterols that bind VDAC in detergent.[2,21] For simplicity, the targeted cholesterol sites are referred to numerically with ascending residue composition (site 1: Lys96, Thr116, Gly117, Asp128; site 2: Phe169, Ala170; site 3: Gly172; site 4: Lys236; site 5: Ser260) (Figure S1).

We docked cholesterol to each site using AutoDock, which gave multiple poses and reasonable docking scores for each site (ΔG ≈ −6.0 to −8.3 kcal/mol). Fifteen atomistic molecular dynamics simulations, averaging 37.3 ns each, were computed in an iterative approach. We sampled multiple starting poses in each site with the goal of optimizing a cholesterol-bound VDAC model (Figure S2). The primary criterion for selecting favorable poses was sustained residency of cholesterol in the sites. Unbinding of cholesterol oriented unfavorably and subsequent diffusion into the membrane was observable on these time scales.

The final cholesterol-bound VDAC model used for subsequent production simulations is shown in Figure 1A. The consistency and fidelity of the model were measured in five 100 ns simulations. As a correlate to chemical shift mapping, we measured the distance between cholesterol and the backbone amide atoms in the sites throughout the trajectories (Figure 1B and Figure S3). In each 100 ns simulation, cholesterol maintained a ∼3–6 Å intermolecular distance to an amide atom from at least three sites. Additionally, each site was essentially occupied for the full 100 ns in at least three simulations (Figure 1B). The unsustained binding we observed in some simulations, however, suggested low affinity interactions, and simultaneous occupation of these five sites cannot be determined at the NMR time scale at which the sites were initially identified.[2]

Figure 1

(A) Optimized cholesterol-bound VDAC model, with cholesterol colored magenta. The N-terminal helix is colored dark blue, and the residues within each site are colored yellow. (B) Minimum distance between cholesterol and an amide atom (N or H) of a residue within that site. A running average of 25 data points was used to reduce noise.

Multiple intermolecular contacts contributed to stable binding in each site (Table 1 and Movies S1–S3). The cholesterol molecules primarily oriented in grooves defined by ridges of hydrophobic and sometimes aromatic side chains that alternately project into the membrane from adjacent β-strands; sites 2 and 3 are defined by a common ridge, as are sites 4 and 5. While we docked cholesterol to the protein in vacuo, we anticipated favorable binding in simulations would be enhanced by the sterol hydroxyl forming hydrogen bonds. The hydroxyls of sites 4 and 5 cholesterols extended into the membrane, contacting water and lipid head groups, and the cholesterol in site 1 formed significant hydrogen bonds with Tyr153. Cholesterol in sites 1 and 3 were also accessible to water funneling from the edge of the protein. Specifically, water was available to site 1 through Tyr153 and the adjacent Ser167, which constitutively contacted water in the channel, and could reach site 3 between the side chains of loop residues His122 and Tyr146, which also contacted bulk water.

Table 1

Residue Contacts in VDAC Cholesterol Sites

site	side chain contacts	backbone contacts	H-bond partner(s)a
1	Leu95, Leu97, Thr116, Tyr118, Leu125, Cys127, Leu142, Met155	Gly117, Gly126, Gly140, Ala141	Tyr153, H2O
2	Ile123, Leu142, Leu144, Ala151, Tyr153, Phe169, Val171, Thr182, Val184, Phe190, Trp210	Gly152, Ala170, Asn183, Gly191	N/A
3	Ile123, Leu144, Tyr146, Trp149, Ala151, Val171	Gly145, Leu150, Gly172	H2O
4	Leu202, Ile221, Ala223, Tyr225, Phe233, Ala235, Val237, Ile243, Leu245	Lys236, Ser234	H2O, DOPC
5	Phe233, Leu245, Tyr247, Leu259, Ala261, Leu263, Leu275	Gly246, Ser260	Asn269, H2O, DOPC

Site 1 cholesterol hydrogen bonded with Tyr153 for 40–50 ns in separate simulations, and site 5 cholesterol hydrogen bonded with Asn269 for 12 ns in one simulation.

In contrast to the above, cholesterol bound the protein parallel to the membrane with the hydroxyl among lipid tails in site 2, stably wedged between side chains, and inaccessible to bulk water and lipid head groups. The initial coordinate of this cholesterol hydroxyl was 1.3 Å from the bilayer midplane. Membrane cholesterol commonly assumes an upright orientation in membranes of saturated phospholipids, with the hydroxyl at the hydrophilic interface ∼16 Å from the bilayer midplane.[65] However, in the presence of polyunsaturated fatty acids, membrane cholesterol can assume a flat orientation in the middle of the bilayer,[65−67] not unlike cholesterol in site 2; therefore, cholesterol could assume this orientation and bind protein without initial placement, especially in more complex membranes.

Our initial simulations were designed to identify favorable orientations of cholesterol in the experimentally identified sites. Validating our approach, we observed replacement of a docked cholesterol by a randomly placed membrane cholesterol, which assumed an identical orientation in site 4 (Movie S4). This should theoretically be observed at all sites; however, measuring multiple unbinding and rebinding events at each site is beyond currently feasible time scales and comes at an uncertain computational cost. While cholesterol potentially binds in other orientations or sites with physiologically relevant affinities, the present model identifies several binding modes that reasonably and reproducibly satisfy experimental binding data.

Cholesterol and VDAC Dynamics

We used the cholesterol-bound VDAC model to investigate the effects of the sterol on the protein. Five 100 ns simulations with an identical mole fraction of cholesterol in the membrane, but not docked to the protein, were computed (the protein in these simulations is referred to as “apo-VDAC” as opposed to “cholesterol-bound VDAC” using our model). Cholesterol did not achieve binding poses in apo-VDAC simulations that were equivalent to cholesterol-bound VDAC. Our experimental approach therefore provided comparable systems that decoupled the effects of cholesterol as a ligand, which was investigated here, from its influence on the membrane. Five 100 ns replicates of VDAC with and without bound cholesterol allowed for a measure of reproducibility and biased our sampling toward the intended apo- or cholesterol-bound states.

The protein backbone equilibrated rapidly and remained stable in the presence or absence of bound cholesterol (Figure S4), and the RMSD of backbone α carbons from averaged structures of apo- or cholesterol-bound VDAC was 0.6 Å. We calculated the average RMSF of the residue α carbons throughout the trajectories (Figure 2A,B). With cholesterol bound, 64 ± 7% of residues per simulation had decreased RMSF compared to the average from apo-VDAC simulations. Residue stabilization by cholesterol was generally global but not absolute; for example, residues 130–138 (between β-strands 8 and 9) were consistently more dynamic (Figure 2C). For any individual residue, however, the absolute difference in averaged RMSF (i.e., the average RMSF from cholesterol-bound simulations minus that from apo-VDAC simulations) was quantitatively small and between −0.63 Å and +0.55 Å.

Figure 2

(A) The average RMSF value for each residue from five simulations of apo-VDAC and (B) cholesterol-bound VDAC systems. Standard error (n = 5) is indicated, and when not visible, the error was smaller than the size of the point. (C) The number of simulations (out of 5) that each residue from cholesterol-bound VDAC simulations had increased α carbon RMSF relative to the average from apo-VDAC simulations. (D) Percent change in RMSF averaged across five cholesterol-bound and five apo-VDAC simulations. Points below the dotted lines indicate that the residue was less dynamic in cholesterol-bound VDAC simulations relative to apo-VDAC simulations.

Generally, the loops that connect β-strands were most dynamic (Figure 2A,B), as previously observed.[68,69] In an aqueous environment, these loops lack the structural support provided by the membrane and inter-residue electrostatic interactions. Even without direct cholesterol binding, the dynamics of the loops were most susceptible to change. The ability of distant loops to be affected by changes in global protein dynamics has been well-characterized in soluble proteins.[70−72] Here, VDAC residues that bind cholesterol might form networks with adjacent amino acids that relay the presence of cholesterol to the loops. Across these loops, the largest percent increases in averaged RMSF were seen in residues Ala134 and Lys274, and the largest percent decreases were seen in residues Pro105, Gly187, and Asn269 (Figure 2D).

We investigated further the effects of cholesterol on the N-terminus, as the helix and cholesterol bind the same wall of the β-barrel, but inside and outside of the pore, respectively (Figure S1). Bonding pairs between helix residues 2–25 and the barrel wall were variable both within and between the two groups (with or without bound cholesterol) (Table S1), and this variability extended to residues that hydrogen bond in the crystal structure.[3] However, the helix maintained 3.1 ± 0.1 and 2.9 ± 0.1 hydrogen bonds per frame with the barrel wall throughout apo-VDAC and cholesterol-bound VDAC simulations, suggesting a dynamic network maintained helical stability. Separate from hydrogen bonding, channel stability is also both dependent on and indicated by a hydrophobic contact between helix residue Leu10 and barrel residue Val143.[39,73] The average distance (∼3 Å) and maximum distance (∼5 Å) between these residues throughout simulations were further consistent with a constitutively open channel for both apo- and cholesterol-bound VDAC. Interestingly, we detected a cholesterol-induced increase in the average minimum channel radius from 5.5 ± 0.1 Å (apo-VDAC) to 5.8 ± 0.1 Å (cholesterol-bound VDAC) using an algorithm that fits spheres inside the channel along the z axis[61] (comparison of minimum channel radii: p < 0.05, two-tailed t test, n = 5 simulations each group). The increased channel radius with cholesterol bound does not necessarily imply greater area for water and ion accessibility; however, it does suggest larger circular or less elliptical constriction sites.

VDAC Electrostatics and Ion Permeability

We generated averaged electrostatic potential maps, which revealed entirely positive net potentials inside the anion channel. Cholesterol-induced protein rigidity resulted in larger regions of strong positive potential, with the strongest potentials in the channel arising from residues with decreased positional fluctuation in cholesterol-bound VDAC (Figure 3). Mutagenesis of multiple residues in these positions modulates ion permeability under low voltage in silico[74,75] and in vitro.[76] Similarly, negatively charged residues that were also stabilized by cholesterol are located along the barrel wall opposite the N-terminal helix, providing a path for K+ ions through cholesterol-bound VDAC. This was seen adjacent to Glu84 (Figure S5), previously shown to influence K+ permeability through VDAC in silico.[64,69]

Figure 3

Electrostatic potential maps showing positive potential inside apo-VDAC and cholesterol-bound VDAC. Blue electrostatic potential maps correspond to (A) +64 kT/e and (B) +38 kT/e. The N-terminal helix is colored green, and the protein backbones correspond to the averaged α carbon structures from the simulations with or without bound cholesterol. The N- and C-termini are at the left entrance (in lateral structures) or on the near side of the channel (when looking into the barrel). In lateral structures, the orange spheres correspond to the α carbons of Leu69 and Ser101 which, for reference, are located at approximately z ≈ −8 and z ≈ 8, respectively. In (A), Lys12, Lys20, and Lys236 are shown in blue sticks, and they are the charged residues that predominantly contribute to the electrostatic field. Arg15, Lys32, Lys119, Lys174, and Lys224 are also shown in (B). All of the mentioned residues are stabilized by cholesterol (i.e., have a decreased average RMSF).

Despite changes in charge distribution, there was no significant change in the average number of Cl– or K+ ions that diffused through the channel in each simulation (number of ions that translocate per simulation — Cl–: 26 ± 3 (cholesterol-bound) and 23 ± 1 (apo-VDAC); K+: 10 ± 2 (cholesterol-bound) and 7 ± 1 (apo-VDAC)). The Cl–/K+ permeability ratios from apo- and cholesterol-bound VDAC were 3.3 ± 0.4 and 2.8 ± 0.7, in reasonable agreement with the experimental value of ∼2 for the open channel.[77] We estimated the potential of mean force for Cl– and K+ as a function of the ion height z in the pore (Figure 4). As described previously,[69,74] free energy wells for Cl– were observed inside the channel on either side of the N-terminal helix, which is the source of the strongest positive potentials (Figure 3). The free energy well at −10 Å < z < −2.5 Å coincided with the location of the minimum pore radius, which was also the location of the free energy maxima for K+. The averaged magnitude and shape of the K+ profiles are similar to previous computational work[69,74] despite variability between simulations at 0 Å < z < 10 Å, possibly due to under-sampling of K+ inside the anion channel.

Figure 4

Potential of mean force for (A) Cl– ions and (B) K+ ions along the z axis, with and without bound cholesterol, estimated using eq 1.

Conclusions

By examining the interactions of VDAC with cholesterol, which is both a membrane component and endogenous ligand,[31] we demonstrate stable binding and alteration of the dynamics of the mitochondrial anion channel. Guided by experimental data, we identified reasonable binding orientations for cholesterol in specific VDAC sites. Although we observed replicable cholesterol binding, it is unclear from our simulations whether the five investigated sites would be simultaneously occupied in vivo. However, occupation of all five sites does not appear necessary to affect VDAC dynamics. Experimentally, in addition to modifying the starting positions of molecules, the stoichiometry of site occupancy can be modulated by the replacement rate of vacated sites by cholesterol from the membrane, which in turn is affected by membrane content. Mitochondrial outer membrane cholesterol content varies in disease states to potentially modulate VDAC function, with 2–10-fold increases in cholesterol per mg of protein observed,[26−28] suggesting higher site occupancy under these conditions. However, varying membrane cholesterol content might also affect VDAC through membrane-mediated mechanisms in addition to direct effects of specific protein binding,[6] which was investigated here. Regulation of membrane protein function by lipid composition can likely occur by exerting the physical properties of the membrane and its components through specific protein sites. Studies of β-barrels, whether computational or experimental, should therefore occur in the presence of at least physiologically relevant cholesterol concentrations. This will be particularly critical as investigation proceeds toward clinically relevant and structurally similar sterols that also bind these proteins,[17] potentially at similar sites.