Ceramides bind VDAC2 to trigger mitochondrial apoptosis

Ceramides draw wide attention as tumor suppressor lipids that act directly on mitochondria to trigger apoptotic cell death. However, molecular details of the underlying mechanism are largely unknown. Using a photoactivatable ceramide probe, we here identify the voltage-dependent anion channels VDAC1 and VDAC2 as mitochondrial ceramide binding proteins. Coarse-grain molecular dynamics simulations reveal that both channels harbor a ceramide binding site on one side of the barrel wall. This site includes a membrane-buried glutamate that mediates direct contact with the ceramide head group. Substitution or chemical modification of this residue abolishes photolabeling of both channels with the ceramide probe. Unlike VDAC1 removal, loss of VDAC2 or replacing its membrane-facing glutamate with glutamine renders human colon cancer cells largely resistant to ceramide-induced apoptosis. Collectively, our data support a role of VDAC2 as direct effector of ceramide-mediated cell death, providing a molecular framework for how ceramides exert their anti-neoplastic activity.

S phingolipids are essential components of eukaryotic membranes that participate in a broad range of cellular processes by controlling vital physical membrane properties and as signaling molecules in response to physiological cues and stresses [1][2][3] . Notably ceramides, the central intermediates of sphingolipid metabolism, have emerged as key mediators of antiproliferative and tumor suppressive cellular programs such as apoptosis, mitophagy, cell cycle arrest, and senescence 4,5 . Multiple stress stimuli, including tumor necrosis factor α (TNFα) 6,7 , ionizing radiation 8,9 , and chemotherapeutic drugs 10,11 , cause a rise in ceramide levels through stimulation of de novo ceramide synthesis, activation of sphingomyelin hydrolysis, or both. Interventions that suppress ceramide accumulation render cancer cells resistant to these stress-inducing agents, indicating that ceramides are bona fide anti-proliferative and pro-apoptotic signaling molecules. While these findings raised considerable interest in targeting sphingolipid-metabolizing enzymes for cancer therapy 5,12,13 , the mechanisms by which ceramides execute their tumor-suppressive activities are incompletely understood.
Mitochondria serve a central role in apoptosis induced by stress stimuli. The mitochondria from cancer cells are often resistant to induction of mitochondrial outer membrane permeabilization (MOMP), a point of no return in the intrinsic pathway of apoptosis. MOMP allows the passage of intermembrane space proteins such as cytochrome c to activate caspases, a family of cysteine proteases responsible for executing an ordered destruction of the cell 14 . MOMP is controlled by pro-and anti-apoptotic members of the B-cell lymphoma 2 (Bcl-2) protein family, which collectively determine the balance between cell death and survival 15,16 . The main function of the anti-apoptotic Bcl2 proteins is to counter the pro-apoptotic activities of the Bcl-2 proteins Bax and Bak, which directly mediate MOMP by creating proteolipid pores responsible for cytochrome c release 17,18 . Several reports have indicated that ceramides can trigger MOMP by modulating the activity of kinases or phosphatases implicated in controlling Bcl-2 protein function. In cells, elevated ceramide levels have been shown to inhibit phosphoinositide-3-kinase (PI3K) and Akt/PBK signaling, resulting in dephosphorylation and subsequent activation of proapoptotic Bcl-2-family protein Bad 19,20 . Short-chain ceramides can bind and stimulate protein phosphatase 2A (PP2A), which dephosphorylates and inactivates the anti-apoptotic protein BCL2 21,22 .
Other studies revealed that ceramides can also act directly on mitochondria to trigger MOMP and apoptotic cell death 23 . For instance, mitochondrial targeting of a bacterial sphingomyelinase to generate ceramides in mitochondria or directing CERTmediated ceramide transport to mitochondria induces cytochrome c release and apoptosis 24,25 . In addition, ER-like membranes associated with isolated mitochondria appear to produce sufficient amounts of ceramides to enable a transient passage of cytochrome c across the outer membrane 26 . However, the underlying mechanisms remain to be established. Interestingly, ceramides have been shown to form pores in model bilayers as well as in the outer membrane of isolated mitochondria that are large enough to mediate passage of cytochrome c 27,28 . Formation of ceramide channels does not rely on any particular protein but is disrupted by anti-apoptotic Bcl-2 proteins 29 . It has also been suggested that ceramides accumulating in the mitochondrial membrane of mammalian cells upon irradiation form or stabilize microdomains that serve as platforms into which Bax inserts and assembles into an active pore 30 . While ceramides have been proposed to cooperate directly with Bax in the assembly of cytochrome c conducting channels 30,31 , other experiments with isolated mitochondria suggest that metabolic conversion of ceramides into sphingosine-1-phosphate and hexadecenal is necessary to facilitate Bax/Bak activation leading to MOMP 32 .
In this study, we present evidence for an alternative mechanistic view, namely that ceramides mediate their pro-apoptotic activity at least in part by interacting directly and specifically with the voltage-dependent anion channel VDAC2, a mitochondrial platform for Bax/Bak translocation [33][34][35] . Identification of VDAC2 as an effector of ceramide-mediated cell death provides new opportunities for exploiting the therapeutic potential of ceramides as tumor suppressor lipids.

Results
A chemical screen for ceramide-binding proteins yields VDACs. To identify proteins involved in ceramide-mediated stress signaling and apoptosis, we used a bifunctional ceramide analog carrying a photoactive diazirine and clickable alkyne group in its N-linked acyl chain (pacCer, Fig. 1a) 36 . Total membranes from human HeLa cells were incubated with pacCercontaining liposomes, subjected to UV crosslinking and click reacted with Alexa Fluor 647 azide (AF647-N 3 ). In-gel fluorescence (IGF) analysis revealed a subset of membrane-bound proteins with affinity for the pacCer probe, which included a mitochondria-associated protein of~33 kDa that was prominently photolabeled (Fig. 1b, Supplementary Fig. 1). As ceramide exerts its apoptogenic activity in mitochondria 8,24,25 , we set out to identify the 33 kDa candidate ceramide-binding protein (CBP). To this end, mitochondria were photolabeled with pacCer and then click reacted with a PEG-based reagent containing an azide, a biotin and a TAMRA fluorophore as functional groups (Fig. 1c). Next, pacCer-crosslinked proteins were isolated using Neu-trAvidin agarose and visualized by IGF (Fig. 1d). The fluorescent 33 kDa protein band was cut from the gel, trypsin-digested, and identified by LC-MS/MS analysis as the voltage-dependent anion channel isoforms VDAC1 and VDAC2 (Supplementary Table 1). In line with the MS data, pretreatment of HeLa cells with VDAC1 and VDAC2-targeting siRNAs effectively depleted the fluorescent 33 kDa protein band from pacCer-labeled and Alexa click-reacted mitochondria (Fig. 1e, f). The photolabeled 33 kDa protein band also cross-reacted with both anti-VDAC1 and anti-VDAC2 antibodies (Fig. 1g, h). In contrast, VDAC3 and TOM40-an outer mitochondrial membrane (OMM) channel protein with a cellular copy number close to that of VDACs 37 -both lacked affinity for pacCer, indicating that VDAC1 and VDAC2 are genuine mitochondrial CBPs.
MD simulations uncover a ceramide-binding site on VDACs.
To search for a putative ceramide-binding site on VDAC1 and VDAC2, we performed coarse-grain molecular dynamics (CG-MD) simulations using the Martini model [38][39][40] . Main simulations were performed with VDAC channels at an aggregate time of 1.23 ms (Supplementary Table 2)-only attainable using CG-MD. A well-resolved structure of mouse VDAC1 (PDB: 4C69) 41 was used as a base template. An available structure of VDAC2 from zebrafish (PDB:4BUM) 42 showed almost perfect structural identity to VDAC1 (1.7 Å barrel backbone RMSD). From the assumption of identical secondary structure we mutated VDAC1 side chains to the mouse sequences of VDAC2 and VDAC3 to obtain all three isoforms for comparison. A bilayer mimicking the OMM 43 was constructed with~630 lipids. To the outer leaflet of this bilayer, 16 molecules of C 16 -ceramide were added. In agreement with the photolabeling data, simulations revealed that VDAC1 and VDAC2, but not VDAC3, have a binding site for ceramide buried in the membrane interior on one side of the barrel wall, comprising β-strands 3-5 (Fig. 2a). This site harbors a uniquely positioned glutamate (Glu) residue in the transmembrane region of β-strand 4-Glu73 in VDAC1 and Glu84 in VDAC2-that faces the bilayer's hydrophobic core (Fig. 2b). In its deprotonated state, this residue seemed to promote direct contact with the ceramide head group ( Fig. 2c; Supplementary Information Videos 1 and 2). In VDAC3, which does not bind ceramide (Figs. 1h and 2a), the bilayer-facing Glu residue is replaced by a glutamine (Gln73; Fig. 2b). Substitution of Gln for Glu73 in VDAC1 or Glu84 in VDAC2 strongly reduced the ceramide occupancy and residence time at the binding sites (Fig. 2d, e). Protonation of the bilayer-facing Glu also greatly diminished ceramide-binding (Fig. 3a, b) while substitution of a deprotonated aspartate for the Glu residue retained ceramide binding ( Supplementary Fig. 2). This indicates that a negative charge on the membrane-buried Glu residue is critical for ceramide binding.
In line with a previous computational study 44 , we also found a number of binding sites for cholesterol. These displayed no overlap with the ceramide-binding site (Fig. 2d). To exclude the possibility that competition with ceramide prevented cholesterol to occupy the ceramide-binding site, additional simulations were performed in the absence ceramide. Also in those simulations no cholesterol binding near the membrane-facing Glu was observed ( Supplementary Fig. 3). We observed zero-specific binding events between phosphatidylcholine (PC) and VDACs in the OMM mimics. Yet when VDAC1 was simulated in a bilayer of 100% dimyristoyl-phosphatidylcholine (DMPC) following the setup of a previous atomic resolution simulation analysis 45 , contacts of the PC head group with the bilayer-facing deprotonated Glu residue could occasionally be observed ( Supplementary Fig. 4a). However, these encounters were much shorter-lived (≤5 ns; Supplementary Fig. 4b) and extremely rare in comparison to those involving ceramide, which displayed average residence times of 0.8 and 1.2 μs for VDAC1 and VDAC2, respectively (Supplementary Fig. 4c). In sum, CG-MD simulations revealed that VDAC1 and VDAC2 each harbor a binding site for ceramide, with a charged Glu residue buried in the membrane interior serving a critical role in ceramide binding.
Ceramide binding relies on a membrane-facing glutamate. We next sought to verify the relevance of the membrane-facing Glu residue in VDACs for ceramide binding. To this end, human VDAC1 and VDAC2 and the mutant channels VDAC1 E73Q and VDAC2 E84Q were produced recombinantly in E. coli and then reconstituted in egg PC liposomes ( Supplementary Fig. 5). Density gradient fractionation analysis revealed that reconstitution efficiencies of wild type and mutant channels were practically indistinguishable. The reconstituted channels were then subjected to photolabeling with pacCer and bifunctional analogs of diacylglycerol (pacDAG), PC (pacPC), phosphatidylethanolamine (pacPE), and cholesterol (pacChol; Supplementary Fig. 6). VDAC1 and VDAC2 could be efficiently and reproducibly photolabeled with pacCer, pacPC, and pacChol, but not with pac-DAG or pacPE (Fig. 4). In agreement with the simulations, replacing the membrane-facing Glu with Gln virtually abolished labeling of both channels with pacCer and pacPC, whereas labeling with pacChol was not or only slightly affected (Figs. 4 and 5a). Moreover, reducing the pH from 7 to 5 caused a significant reduction in E73-dependent photolabeling of VDAC1 with pacCer (Fig. 3c, d). This suggests that ceramide binding is critically dependent on the protonation state of the membraneexposed Glu, as predicted by the simulations. The pronounced labeling of the wild-type channels with pacPC was somewhat unexpected as simulations revealed that specific  c Stills from an MD simulation, showing the approach and binding of a ceramide molecule to VDAC1 in close proximity of the bilayer-facing Glu residue in its deprotonated state. Protein surface colors mark polar (green), apolar (white), cationic (blue), or anionic (red) residues. d Space-filling and wireframe models of VDAC1, VDAC1 E73Q , VDAC2, and VDAC2 E84Q with deprotonated E73/E84. Indicated are the volumes for which there is ceramide occupancy greater than 10% (orange) or cholesterol occupancy greater than 20% (yellow). e Distribution of the durations of ceramide contacts with VDAC1, VDAC1 E73Q , VDAC2, and VDAC2 E84Q at the preferred binding site as in d.
The y-axis indicates the fraction of the total system time spent in binding events of the duration indicated by x. Summing all points' y-values yields the fraction of total simulation time when ceramide was bound encounters of PC molecules with the membrane-facing Glu are relatively rare. This discrepancy could be due to the irreversible nature of photoaffinity labeling in combination with the major difference in timescales between CG-MD simulations (μs range) and photoaffinity labeling (second range). While our data revealed no obvious overlap between the ceramide and cholesterol binding sites on VDACs, a recent photolabeling study mapped the bilayer-facing Glu to a major cholesterol binding pocket on VDAC1 46 . In contrast to the present work, the latter study was performed on VDAC1-containing bicelles with cholesterol probes that carry the photoactive diazirine in the aliphatic tail or at C7. Thus, additional work will be necessary to resolve the discrepancy between the simulations and photolabeling studies of cholesterol binding to VDACs. As aliphatic diazirines display photochemical preference for nucleophilic amino acids 47 , the deprotonated side chain of the membrane-buried Glu in VDACs may provide a site of insertion for the diazirine in pacCer, which may diffuse a short distance to reach such a nucleophile. Therefore, we next labeled VDACs with pacCer in the presence of excess C 16 -ceramide. Labeling by pacCer was progressively reduced by C 16 -ceramide when added in 3-to 27-fold excess (Fig. 5b), arguing against the idea that pacCer labeling of VDACs is primarily driven by affinity of the aliphatic diazirine for the negatively charged Glu. Excess C 16ceramide did not affect labeling of VDACs with pacChol. In line with the simulations, these data indicate that ceramides bind VDACs at the membrane-facing Glu residue. To verify that this concept not only holds for recombinant channel proteins in synthetic bilayers but also for their native counterparts in mitochondrial membranes, we made use of the carboxylmodifying reagent dicyclohexylcarbodiimide (DCCD). This hydrophobic compound irreversibly reacts with a number of integral membrane proteins through covalent modification of membrane-embedded Asp or Glu residues and was previously shown to modify Glu73 in VDAC1 48 . Pretreatment of mitochondria with DCCD selectively abolished photolabeling of the 33 kDa protein band with pacCer ( Supplementary Fig. 7), hence providing complementary proof that the membrane-facing Glu residue in VDAC1 and VDAC2 is part of an authentic ceramide-binding site.
Loss of VDAC2 disrupts ceramide-induced apoptosis. VDAC channels are active participants in the cytosolic release of apoptogenic proteins from mitochondria. VDAC2 serves as a platform for the mitochondrial recruitment of pro-apoptotic Bcl-2 proteins Bak and Bax [33][34][35] , which can commit cells to death by permeabilizing the OMM for cytochrome c 17,18 . In response to various apoptotic stimuli, VDAC1 forms oligomers and participates in the assembly of a cytochrome c-conducting pore 49,50 . Identification of a ceramide-binding site on VDAC1 and VDAC2 raised the question whether ceramides exert their apoptotic activity by interacting with these proteins. To address this, we employed an engineered ceramide transfer protein equipped with an OMM anchor, mitoCERT (Fig. 6a). We previously demonstrated that human HCT116 colon cancer cells expressing mitoCERT undergo Bax-dependent apoptosis by mistargeting newly synthesized ER ceramides to mitochondria (Fig. 6b) 25 . This led us to determine the impact of VDAC removal on mitoCERT-induced apoptosis in HCT116 cells. Loss of VDAC1, VDAC2, or both was verified by immunoblotting and IGF analysis of mitochondria photolabeled with pacCer ( Supplementary Fig. 8). Expression of mitoCERT in wild-type HCT116 cells triggered apoptosis, as indicated by cleavage of caspase substrate PARP1 (Fig. 6c). No PARP1 cleavage was observed in cells expressing a mitoCERT variant that lacked the ceramide transfer or START domain, mitoCERTΔSTART. Loss of VDAC1 had no effect on the ability of mitoCERT to induce PARP1 cleavage. In contrast, VDAC2 removal rendered cells partially resistant to mitoCERT-induced PARP1 cleavage, especially in the absence of VDAC1 (Fig. 6c). c Human VDAC1 and VDAC1 E73Q were produced in E. coli, purified, reconstituted in liposomes, and then photolabeled at the indicated pH with pacCer added from an ethanolic stock. Samples were click-reacted with AF647-N 3 , subjected to SDS-PAGE, and analyzed by IGF and CB staining. d Quantitative analysis of relative pacCer photolabeling efficiencies of reconstituted VDAC1 treated as in c. Data are means ± s.d.; n = 4; **p < 0.01 by two-tailed paired t-test. Source data Glu84 in VDAC2 is critical for ceramide-induced apoptosis. We next analyzed individual wild type and ceramide-binding defective VDAC channels for their ability to support mitoCERTinduced apoptosis. To this end, HCT116 VDAC1/2 double KO cells were stably transduced with haemagglutinin (HA)-tagged VDAC1, VDAC1 E73Q , VDAC2, or VDAC2 E84Q (Fig. 7a). Mitochondrial localization of the tagged channels was confirmed by immunofluorescence microscopy (Supplementary Fig. 9). Heterologous expression of VDAC2, but not VDAC1, restored mitoCERT-induced apoptosis in the double KO cells, as evidenced by PARP1 cleavage and proteolytic activation of caspase-3, the principal executioner caspase in apoptotic cells (Fig. 7a-c). This confirmed that VDAC2, unlike VDAC1, is a key player in ceramide-mediated cell death. Strikingly, replacing Glu84 with Gln in VDAC2 greatly reduced its ability to support mitoCERTinduced apoptosis in the double KO cells (Fig. 7a-c), indicating that cell death triggered by a rise in mitochondrial ceramides requires a ceramide-binding competent form of VDAC2. While VDAC1 or VDAC2 removal had no major impact on cellular Bax levels, cells lacking VDAC2 had strongly reduced Bak levels ( Supplementary Fig. 10a), consistent with the unique role of VDAC2 in stabilizing Bak 33 . Reintroducing VDAC2 or VDA-C2 E84Q in VDAC1/2 double KO cells in each case fully restored Bak levels to those in wild-type cells ( Supplementary Fig. 10b), indicating that an intact ceramide-binding site is dispensable for VDAC2-mediated stabilization of Bak. This notion is supported by the finding that the Bak-stabilizing activity of VDAC2 requires isoform-specific sequence motifs located outside of the region involved in ceramide binding 33 .

Discussion
The current study identified a role of VDAC2 as a direct and specific effector of ceramide-mediated cell death. This function critically relies on a uniquely positioned charged Glu residue that mediates direct contact with the ceramide head group in the bilayer interior, potentially driven by electrostatic attraction. We speculate that an amide-containing backbone combined with a small polar head group renders ceramide the preferred lipid binding partner of VDAC isoforms containing the membraneburied Glu residue. Although it is energetically unfavorable for charged residues to face the bilayer's hydrophobic core, recent work revealed that the pKa of Glu73 in VDAC1 is closely tuned to the physiological pH of the cytosol (pKa~7.4) 51 . This implies that under stress-free conditions, the membrane-buried Glu is in its deprotonated fully charged state at least a significant amount of time. Interestingly, ceramide binding to the late endosomal protein LAPTM4B depends on a membrane-embedded aspartate 52 . We anticipate that also other proteins with charged acidic residues in their membrane spans may have affinity for ceramide and potentially participate in ceramide-operated signaling pathways. While ceramide-induced apoptosis requires Bax 8,24,25 , a recent study revealed that VDAC2 specifies Bax recruitment to    Fig. 4 The bilayer-facing Glu is a critical determinant of pacCer photolabeling of VDACs. a Human VDAC1, VDAC1 E73Q , VDAC2, and VDAC2 E84Q were produced in E. coli, purified, reconstituted in liposomes, and incubated for 30 min at 37°C with liposomes containing 1 mol% of pacCer or photoactive and clickable analogs of diacylglycerol (pacDAG), phosphatidylcholine (pacPC), phosphatidylethanolamine (pacPE), or cholesterol (pacChol). Samples were UV irradiated and then click-reacted with AF647-N 3 , subjected to SDS-PAGE, and analyzed by IGF and CB staining. b Quantitative analysis of relative labeling efficiencies of reconstituted VDAC1, VDAC1 E73Q , VDAC2, and VDAC2 E84Q with pacLipids as indicated in a. Data are means ± s.d.; n ≥ 3. Source data mitochondria and concomitantly ensures Bax inhibition by mediating its retrotranslocation into the cytosol 34 . By establishing a dynamic equilibrium between mitochondrial and cytosolic Bax pools, this VDAC2-dependent shuttling is ideally suited for regulation by pro-and anti-apoptotic cues. Our present findings suggest that ceramide binding to VDAC2 may commit cells to death by blocking Bax retrotranslocation. This concept is distinct from previous models postulating that ceramides accumulating in the OMM: (i) self-assemble into cytochrome c-conducting channels 27,28 ; (ii) form lipid macrodomains into which Bax inserts and functionalizes as a pore 30,31 ; (iii) affect mitochondrial shape to facilitate Bax recruitment and apoptosis 53 ; (iv) require metabolic conversion to gain apoptogenic activity 32 .
How ceramide binding tips the balance in VDAC2-mediated shuttling of Bax to trigger mitochondrial apoptosis remains to be established. Our simulations and photoaffinity experiments indicate that ceramide binding to VDACs is pH sensitive and controlled by the protonation state of the membrane-buried Glu. Interestingly, acidification has been shown to promote association of two VDAC1 monomers into a dimer 51 . Assembly of a highaffinity dimer relies on protonation of the membrane-facing Glu (E73) and likely involves hydrogen bonding between this residue and a serine (Ser43) at the dimer interface 51 . VDAC oligomerization has been implicated as a component of the mitochondrial pathway of apoptosis 49,50 and may be part of the mechanism by which VDAC2 stabilizes the mitochondrial pool of Bax in response to apoptotic stimulation 34 . In view of our present findings, a prospect that merits further investigation is whether binding of ceramide to the charged Glu residue lowers the threshold for VDAC oligomerization at neutral pH. Ceramide binding to VDACs may also influence interactions with other proteins, as the membrane-facing Glu residue is critical for association of VDAC1 with hexokinase I 54,55 . VDAC-bound hexokinases are thought to play a pivotal role in promoting cell growth and survival in rapidly growing, hyperglycolytic tumors 56,57 . Consequently, our current study establishes a molecular framework to unravel how ceramides execute their tumor suppressor functions.
MD simulations. For VDAC MD simulations, a well-resolved structure of mouse VDAC1 (PDB: 4C69 at 2.8 Å resolution) 41 was selected as a main template. This structure shares a β-barrel backbone RMSD of 1.7 Å with a structure of VDAC2 from zebrafish (PDB: 4BUM) 42 . Based on the assumption of identical secondary structures, VDAC1 side chains were mutated to the mouse sequences of VDAC2 (NCBI ID:NP_035825.1) and VDAC3 (NCBI ID:NP_035826.1) to obtain structures of all three isoforms for comparison using the PyMOL software. The mutation process was semi-automated and leveraged the used CG representation (see below), in which side chains are represented by up to four particles and are therefore simple to replace. The N-terminal α-helix of the VDAC2 was truncated to match VDAC1 sequence length (with no impact to lipid interactions because the helix sits inside the channel's β-barrel). VDACs were embedded in an outer mitochondrial membrane mimic with about 630 lipids using the insane script 59 16 -ceramide in the outer leaflet. About 150 mM NaCl was added to the system, plus an excess of Na + ions to reach charge neutrality. The Martini coarse-grain (CG) force field 60 was used to model the simulated systems, together with the ElNeDyn elastic-network approach to restrain the protein secondary structure 61 . Simulations were run with GROMACS version 5.0 62 in the isothermal-isobaric (NpT) ensemble, at 300 K and 1 bar. Temperature was controlled using the vrescale thermostat with a coupling constant of 1.0 ps. Pressure was coupled semiisotropically, independently in xy and z, using the Parrinello-Rahman barostat with a coupling constant of 12.0 ps and a compressibility of 4.5 × 10 −5 bar −1 . Standard Martini parameters were used to represent interparticle interactions 60,63 : electrostatic interactions were modeled with a Coulombic potential cutoff at 1.2 nm, shifted to zero along the entire range, in conjunction with an implicit screening constant of 15. The Van der Waals interactions were calculated using a shifted Lennard-Jones potential, with a cutoff of 1.2 nm and a shift to zero from 0.9 nm. A simulation integration time step of 20 fs was used, and particle neighborsearching was performed over a 1.4 nm radius once every 10 simulation steps. As the apparent pKa of the membrane-facing glutamate E73 in VDAC1 is closely tuned to the physiological pH of the cytosol (pKa~7.4) 51 , main simulations of VDAC1 and VDAC2 were initially done with deprotonated E73/E84 in six replicates (0.3 ms total) each. Other simulated systems include VDACs with protonated E73/E84, or with E73Q/E84Q or E73D/E84D mutations. VDAC1 with deprotonated E73 was also simulated in 100% DMPC. The total simulation time for all systems was 1.23 ms (Supplementary Table 2). Trajectories were saved and analyzed every 600 ps.
Simulation analysis. Two main types of analysis were performed: contact analysis and space occupation analysis. Contact analyses were performed either as residuediscriminated analyses or as residence-time distributions. Residue-discriminated contacts were calculated by finding all frames of ceramide (particles AM1 or AM2) or cholesterol (particle ROH) head groups within 7 Å of any relevant residue particle. To allow comparison between simulations of different lengths, the number of in-contact frames was normalized by the total number of frames. Contact residence-time distributions were plotted by logarithmically histogramming the duration of all individual continuous binding events of each single ligand to the VDAC-binding site. To correctly convey the impact of long-term binding events, yvalues indicate the sum of binding durations in that bin. To compare between simulations of different duration, the y-values were further normalized for the respective total simulation time. Analysis of continuous contact times is very sensitive to the use of a sharp cutoff, and long binding events will be underestimated if the ligand briefly drifts past the threshold, even if not actually leaving the binding site. To correct for this, for binding duration analysis a contact was defined when a head group particle (AM1 or AM2 in ceramide; PO4 or NC3 in DMPC) came within 8 Å of at least two groups of residues, out of three defined on the β-strands that compose the binding site (VDAC1 residues 58-60, 73-75, and 81-83). Spuriously brief contacts, or spuriously brief losses of contact, were further filtered by applying a smoothing window five-frames wide; if there were contacts on three or more of those frames the center frame was set to be in-contact, otherwise it was set to be not-in-contact. Binding events that were still bound at the end of a simulation run were disregarded. Space occupation analyses were performed by gridding the simulation box at 5 Å spacing and comparing ceramide versus cholesterol occupancy. The number of occupied frames for each cell was calculated for the entire simulation time, and then normalized by the total number of frames, yielding a grid of occupancy values ranging from zero (never present in that cell) to one (always present). Occupancy volumes were drawn as isooccupancy surfaces at given thresholds (10% for ceramide, 20% for cholesterol), adjusted to compensate for the difference in the number of lipid molecules. To allow meaningful occupancy surfaces to be drawn relative to the protein, its position and rotation in the xy plane were aligned beforehand across all analyzed trajectories.
Algorithms and statistical analysis. Analysis algorithms were programmed in Python, using MDAnalysis (https://www.mdanalysis.org) 64,65 and NumPy 66 packages. The MDreader package (https://github.com/mnmelo/MDreader) was used to allow analysis parallelization, which was essential in tackling the multigigabyte trajectory analysis in a short time. Two types of statistical analyses were performed on the distribution of contact residence times: data are presented as time-weighted average values, with 95% confidence intervals estimated by a 10kresample bootstrapping of the time-weighted averaging, and the significance level for the difference between binding durations was obtained by a continuitycorrected one-tailed Mann-Whitney-Wilcoxon U test. The Visual Molecular Dynamics (VMD) software 67 was used to create images and movies.
Reporting summary. Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Data supporting the findings of this manuscript are available from the corresponding authors upon reasonable request. A Reporting Summary for this Article is available as a Supplementary Information file. Full scans of gels and blots are provided in Supplementary Information