Structural basis of complex formation between mitochondrial anion channel VDAC1 and Hexokinase-II

Complex formation between hexokinase-II (HKII) and the mitochondrial VDAC1 is crucial to cell growth and survival. We hypothesize that HKII first inserts into the outer membrane of mitochondria (OMM) and then interacts with VDAC1 on the cytosolic leaflet of OMM to form a binary complex. To systematically investigate this process, we devised a hybrid approach. First, we describe membrane binding of HKII with molecular dynamics (MD) simulations employing a membrane mimetic model with enhanced lipid diffusion capturing membrane insertion of its H-anchor. The insertion depth of the H-anchor was then used to derive positional restraints in subsequent millisecond-scale Brownian dynamics (BD) simulations to preserve the membrane-bound pose of HKII during the formation of the HKII/VDAC1 binary complex. Multiple BD-derived structural models for the complex were further refined and their structural stability probed with additional MD simulations, resulting in one stable complex. A major feature in the complex is the partial (not complete) blockade of VDAC1’s permeation pathway, a result supported by our comparative electrophysiological measurements of the channel in the presence and absence of HKII. We also show how VDAC1 phosphorylation disrupts HKII binding, a feature that is verified by our electrophysiology recordings and has implications in mitochondria-mediated cell death.

M itochondria are the primary source of ATP production in eukaryotic cells. Newly generated ATP is transported out of the mitochondrial matrix via the adenine nucleotide transporter at the inner mitochondrial membrane and out to the cytosol through the voltage-dependent anion channel (VDAC) located in the outer membrane of mitochondria (OMM). VDAC1, the most abundantly expressed isoform of VDAC, serves as the main conduit for the large flux of ions, ATP/ ADP, metabolites, organic anions, and various respiratory substrates across the OMM [1][2][3] . In addition to its function as the major gateway in and out of mitochondria, VDAC1 also acts as a scaffold, regulating mitochondrial-triggered apoptotic signaling through interactions with a variety of proteins 2,4-6 .
VDAC1 has been reported to bind to pro-or anti-apoptotic proteins that can alter the permeability of the OMM and promote or prevent cell death [7][8][9][10][11][12][13][14][15] . For example, VDAC1 serves as a receptor for the cytosolic anti-apoptotic protein hexokinase (HK) which enhances cell survival 7,8 . Interaction with VDAC1 not only modulates cell survival but also gives HK preferential access to ATP for its catalytic activity of phosphorylating glucose to glucose-6-phosphate during glycolysis 16 .
HK is known as one of the primary factors of high glycolytic characteristics of rapidly growing tumor cell 17 . In particular, HK isoforms HKI and HKII are overexpressed in many types of cancers 4,[18][19][20] . Elevated levels of HKI and HKII lead to a high rate of glycolysis, known as the Warburg effect, resulting in the enhanced generation of lactic acid, a key component in promoting cell growth [21][22][23] . Coincidentally, HKI and HKII both differ from other HK isoforms with one extra structural element: in addition to two homologous hexose binding domains (N-domain and C-domain, Fig. 1a), HKI and HKII have a short N-terminal hydrophobic helix (termed H-anchor hereafter, Fig. 1a) that is thought to be capable of membrane-binding 24 . Given the importance of HK/VDAC1 interaction in cell growth and survival, disruption of this protein-protein complex has been identified as a potentially effective therapeutic strategy to prevent rapidly growing tumor cells [25][26][27][28][29][30] .
While a number of VDAC1 residues have been identified to be essential for interaction with HKI and HKII isoforms 7,21,[31][32][33][34][35] , structural details of their complex remain unknown. Two protein-protein docking studies have been reported for modeling the complex formation between HK isoforms and VDAC1 36,37 . Both studies proposed direct plugging of the H-anchors of HKI/ HKII into the pore of the barrel-shaped VDAC1. However, these models failed to address several key characteristics of HK/VDAC1 interactions, e.g., how the binding of HKI only partially reduces the maximal conductance of VDAC1 7,21,25,31 , as the plugged model would result in a complete blockade of the channel. Furthermore, at the docked interface, there is a mismatch between the hydrophobic H-anchor from HKI or HKII and the highly charged VDAC1 interior surface (Fig. 1a, b) 38 .
It is known that HKs can bind to membrane even in the absence of VDAC1 [39][40][41] , and the H-anchor is shown to be essential for this membrane interaction 24,39,40 , likely due to its hydrophobic nature. Immunoblotting of mitochondria-bound HKI and HKII has shown that truncating the hydrophobic portions of the H-anchor disrupts HK binding to native OMM (even in the presence of VDAC1) 42,43 . The significance of the H-anchor in OMM binding is further emphasized by the fact that HK homologs lacking the H-anchor show no mitochondrial localization 20,44,45 . Possibly due to the inability to bind the membrane, truncation of the HKI H-anchor also eliminates the effect on the VDAC1 conductance 7 . Thus, the H-anchor in both HKI and HKII is vital for interaction with the OMM, a step that appears to be a prerequisite for the HK-induced VDAC1 blockade. Henceforward, we hypothesized the following scenario in the complex formation between HKs (HKI or HKII) and VDAC1: first, HK binds to the OMM by inserting its H-anchor, a step that aligns the H-anchor for interaction with and binding to the outerrim of VDAC1, in order to form the complex on the cytosolic surface of the OMM (Fig. 1c).
Based on the hypothesis described above, we designed a modeling approach involving two major steps. In the first step, we describe membrane binding of HKII with molecular dynamics (MD) simulations. Next, the resulting membrane-bound HKII is used to study the formation of its complex with VDAC1 using Brownian dynamics (BD) simulations. Our results show that HKII binds to the membrane by inserting its H-anchor, independent of VDAC1. The resulting complex reveals two key functional features of HKII/VDAC1 interaction: (1) a partial, and not complete, blockage of VDAC1 conductance through HKII interaction and (2) how phosphorylation of VDAC1, and the corresponding phosphomimetic mutation, can disrupt the interaction between the two protein. Both of these predictions resulting from the computational model are supported by our designed electrophysiology measurements and mutagenesis experiments. performed 12 independent 200-ns MD simulations using the highly mobile membrane mimetic (HMMM) model 46 (Fig. 2a, more details in SI), during which binding of the N-domain of HKII (referred to as HKII-N) to lipid bilayers representing the OMM was simulated. Each simulation started with HKII-N initially placed in the solution above the membrane with a different initial orientation (Figs. 2a and S1). Binding of HKII-N to the membrane was observed in 10 out of 12 replicas (Fig. 2b), all with the H-anchor inserting into the membrane. Spontaneous membrane binding occurs within the first 50 ns of the simulations in 9 out of the 10 successful replicates. After the first encounter with the membrane, the H-anchor is rapidly inserted into the membrane and the HKII-N remains membrane-bound for the remainder of the simulation (Fig. 2b). The H-anchor reaches a maximum depth of 10 Å below the bilayer's phosphorus plane (as measured by the insertion of residue I2), penetrating well into the membrane's hydrophobic core (Fig. 2b). Membrane partitioning and orientation (with respect to the membrane normal) of the H-anchor were similar in all membrane-bound replicas (Figs. S1 and S2), suggesting convergence in the obtained membranebound configuration of HKII-N.

Results and discussion
To ensure that the final membrane-bound model of HKII-N obtained from the HMMM simulations was stable, we converted one of the resulting membrane-bound systems to a conventional, full membrane and simulated it for additional 200 ns. The Hanchor remains buried into the hydrophobic core of the full membrane during the simulation (Fig. 2). In fact, the H-anchor appears to become even more engaged with the membrane during this additional stimulation, as evidenced by deeper penetration of I2 into the membrane, reaching a maximum depth of 17 Å below the phosphorus plane of the lipid head groups' approximately at t = 300 ns (Fig. 2b). These results provide strong support and a model for the direct interaction of HKII with the membrane. Membrane partitioning of the H-anchor obtained from the fullmembrane simulation shows that the first 10 residues remain stably bound and/or insert into the membrane (below the phosphorus plane) (Fig. 2c). This finding substantiates the experimental results where truncation of the first 10 residues from HKII was shown to disrupt its OMM binding 42 . Hence, our results provide the first structural model for membranebound HKII.
HKII/VDAC1 complex formation. Once we established the membrane-binding mode of HKII-N using MD simulations, we extended our study to investigate the molecular interaction between membrane-bound HKII-N with membrane-embedded VDAC1 by performing an aggregate of millisecond-scale atomicresolution BD simulations. BD simulations were performed by placing 100 independent replicas of HKII-N around VDAC1 (Fig. 3a-c); each system was simulated for 20 μs, totaling 2 ms. During the simulations, the membrane-bound pose of HKII-N was preserved using restraints designed based on the MD simulation results. These restraints allowed for the lateral diffusion of HKII-N around VDAC1 but prohibited significant vertical displacement of the system with respect to the membrane plane, thereby approximating the membrane-insertion depth and orientation of HKII-N (Fig. 3a, Supplementary Movie 1). Details for the BD simulations are provided in SI text.
The BD simulations revealed five distinct hot spots for HKII-N to bind VDAC1 while remaining anchored in the membrane (Fig. 3d). To further investigate whether these hot spots represent distinct HKII-N/VDAC1 complexes, we performed cluster analysis of HKII-N's relative position to VDAC1 (using C α root-mean-square deviation as a reference) for BD trajectories when the two are in close contact (i.e., having at least 1 atom within 3.0 Å cutoff distance). The top five distinct clusters resulted from our analysis match the locations of the 5 BD hot spots (Fig. 3d). A representative HKII-N/VDAC1 complex, with the most favorable interaction energy between HKII-N and VDAC1, was selected for each of the five clusters (complexes termed HKV1, HKV2, HKV3, HKV4, and HKV5, respectively, for Clusters 1-5). Each HKII-N/VDAC1 complex model was then used to generate a full-length HKII/VDAC1 complex by extending the HKII C-domain from HKII-N. The resulting complexes were then inserted into a full membrane, creating five independent systems for additional MD refinement simulations (Figs. 3e and S3), during which the proteins and the membrane are allowed to relax conformationally, and the interfacial residues can adjust better to the new environment provided by binding of the partner protein. Each system was simulated in two replicates. Among the five complexes, HKV1 maintained the strongest interaction with VDAC1 during the 650 ns of MD refinement (Figs. S3 and S4, Supplementary Movies 2 and 3). In HKV2, relatively weak HKII/VDAC1 interactions were observed (Figs. S3 and S4). Furthermore, the binding mode observed in HKV2 would interfere with the dimeric interface of VDAC1 47 , reducing the probability of its in vivo relevance (Figs. S3 and S5, Supplementary Movies 4 and 5). In contrast, the binding mode observed in HKV1 does not interfere with the proposed dimeric interface of VDAC1 47 (Fig. S5). For the rest of the complexes (HKV2, HKV3, and HKV5), either a complete dissociation or negligible HKII/VDAC1 interactions were recorded during their respective MD refinement simulations (Figs. S3 and S4, Supplementary Movies 6-11). Moreover, HKV1 corresponds to the highest population cluster (Table S1). Therefore, we consider HKV1 as the most relevant system and used it for further structural analysis.
In HKV1, multiple hydrogen bonds and salt-bridge interactions between VDAC1 and HKII (H-anchor and N-domain) maintained the stability of the complex during the MD simulations (Fig. S6). These interactions promoted slightly deeper membrane insertion and a larger tilt angle (relative to the membrane normal) of H-anchor in HKV1 when compared to the membrane-bound configuration of HKII-N alone (Fig. S7). Though overall interactions seem to reach a steady state (Fig. S3), some individual interactions were only maintained intermittently (Fig. S6). This suggests that much longer equilibration may be needed to obtain a properly relaxed structure, as expected considering the fact that complex formation occurs on a minute time-scale (based on current, indirect experimental observation [see below]). Nevertheless, no contacts were observed between the C-domain of HKII and VDAC1 (Fig. S12). This is consistent with previous experiments where mitochondrial binding of HKII (and HKI) was observed with immunoblotting of H-anchor+Ndomain, but not with the C-domain only 43 .
HKII modulates VDAC1 conductance. Upon complex formation, HKII covers a large fraction of the cytosolic surface of VDAC1 (Fig. 4a), primarily due to the contacts from the Ndomain. This coverage results in an almost halved cytoplasmic opening of VDAC1 (reduced from 13.0 to 7.3 Å, Fig. 4b), suggesting that the VDAC1 conductance might be affected upon HKII binding. To further test this hypothesis and quantify the effect, we performed constant-voltage MD simulations for both HKII-free and HKII-bound states of VDAC1. Three independent  (Table S1). e A snapshot (at t = 350 ns) of membrane-embedded HKII/VDAC1 complex 1 (HKV1) during the full-membrane MD simulation.
simulations were performed for both states, each under a −50 mV membrane potential for 80 ns. Ionic currents obtained from these simulations showed that VDAC1 conductance is indeed partially reduced upon HKII binding (Fig. 4c, 36,37 would virtually block the channel permeability completely, as indicated by the drastic narrowing of the pore to a 2.1 Å effective radius and a negligible current in constant-voltage simulations ( Fig. S13), which cannot account for the partial reduction of current observed experimentally upon HKII/VDAC binding. The hallmarked partial blockage of VDAC1 current upon HKII binding demonstrated in the simulations was further corroborated with a series of electrophysiology measurements. Using the planar lipid bilayer method previously described 48 , the reconstitution of wild-type (wt) VDAC1 and its insertion into the bilayer was evidenced by its unique voltage-dependent gating characteristic in response to a voltage ramp protocol from −80 to +80 mV, where the maximal conductance occurs approximately between −50 and +50 mV and lower conductances at the more hyperpolarizing or depolarizing potentials (Fig. 5a). After washing away the remaining soluble VDAC1, the current was monitored under a −30 mV test potential for 30 s at 2 min intervals. When HKII was added to the cis chamber of the planar lipid bilayer (corresponding to the cytosolic side), there was a marked decrease in the maximal conductance of VDAC1 (Figs. 5b, S14). Steady-state levels were reached in~6 min, and this latency may be due to the time it takes for the membrane insertion of HKII and its diffusion at the bilayer before reaching VDAC1. Upon reaching a steady state, the effect of HKII on VDAC1 was stable for the subsequent 18 min of experiment time.
Note that the addition of HKII resulted in only the reduction of VDAC1 conductance but not a complete blockade. These results strongly support the HKV1 model in MD simulations (Fig. 4) where the N-domain of HKII only partially covers the permeation pathway of VDAC1 and partially reduces the conductance.
Phosphorylation of VDAC1 disrupts HKII binding. The immense HKII/VDAC1 interaction energy shown in MD simulations (Figs. S3 and S4) and the composition of the HKII/ VDAC1-binding interface (Fig. S6) strongly suggest the binding interactions between HKII and VDAC1 are dominated by salt bridges. Therefore, their complex formation might be interfered or disrupted with the presence of additional charged species, such as post-translational modifications to the proteins. Coincidently, an earlier proteomics screening of phosphorylated mitochondrial proteins has identified S215 of VDAC as a phosphorylation site in vivo 49 , which is located at the HKII/VDAC1 binding interface (Fig. S6).
To investigate the impact of S215 phosphorylation on complex formation between HKII and VDAC1, we repeated the BD simulations after introducing a phosphoserine residue at position 215. Whereas BD simulation of wt-VDAC1 resulted in an HKIIbinding hot spot near S215 (Cluster 1), the probability of HKII presence at this region is substantially reduced after S215 phosphorylation (Fig. 6b), suggesting that the formation of HKII/VDAC1 complex can be inhibited when S215 is phosphorylated. The disappearing of population corresponding to Cluster 1 can also be predicted by the greatly reduced HKII/ VDAC1 interaction energy if a phosphoserine is directly placed at S215 position of the BD trajectories of wt-VDAC1: the electrostatic interactions between the two proteins suffered from 2.6-fold reduction on average in Cluster 1, while other clusters showed only minor changes (Fig. S15).
VDAC1 phosphorylation is known to modulate HKII binding in vivo 40 . Therefore, the impact of VDAC1 S215 phosphorylation to HKII/VDAC1 interactions was further corroborated with an electrophysiology experiment using a phosphomimetic mutant of VDAC1, S215E. Comparing to wt-VDAC1, the HKII-induced The average radius profile of the VDAC1 pore calculated using HOLE 95 for both VDAC1 alone, and its complex with HKII (Complex HKV1). c The in silico VDAC1 current at −50 mV for a membrane-equilibrated VDAC1 alone or the VDAC1/HKII HKV1 complex calculated from electric-field MD simulations. The mean and standard deviation of current for each system were calculated from three independent simulations. d The cumulative net number of channel-crossing events by Cl − (red traces) and Na + (blue traces), tracked over the time course of a representative MD simulation for VDAC1 or the HKV1 complex.
partial blockage is rarely observed, if not completely abolished, in membranes containing the S215E phosphomimetic mutant of VDAC1 (Figs. 5b and S14). Additionally, the presence of HKII does not alter the current-voltage characteristics of S215E-VDAC1 significantly (Fig. 5a). Based on these observations, the S215E point mutation appeared to prevent or diminish the interaction of HKII with VDAC1, supporting the predicted outcome of the computational model. To date, it remains unclear which kinase or signaling pathway mediates VDAC1's phosphorylation at S215 position. We used the peptide sequence near S215 for kinase substrate prediction and did not yield a convincing kinase hit. However, since HKII/ VDAC1 binding can be inhibited by the S215 phosphorylation, the kinase responsible for S215 phosphorylation is expected to belong to a pro-apoptotic pathway or at least not in a pathway promoting cell survival. Intriguingly, VDAC1 is a known substrate of glycogen synthase kinase (GSK)-3β 35 , a major promoter of mitochondrial intrinsic apoptotic pathway 50 . VDAC phosphorylation by GSK3β is also linked to a reduced HKII binding 35 , and the disruption of HK binding to VDAC1 has been shown to play an essential role in several diseases including cardiac ischemia-reperfusion injury 51-53 , by promoting mitochondria-mediated cell death 7,35,54 . The matching phenotype leads us to speculate that the kinase responsible for VDAC1 S215 phosphorylation might be GSK3β.
GSK3β is known to highly favor a "primed" substrate where another phosphorylated serine/threonine is located at the +4 residue position downstream of the phosphorylation site 55 . One may argue that in VDAC1, such position is a membraneembedded F219, which cannot be phosphorylated and serve as a priming residue for GSK3β recognition (Fig. 6a). Based on crystal structures of GSK3β with a bound pseudo-substrate phosphopeptide inhibitor 56 , GSK3β seems to recognize substrate peptides with an extended backbone configuration where the priming phosphate sits~13 Å away from the transferring γ-phosphate. Interestingly, the β-barrel VDAC is full of extended backbones and its cytoplasmic rim has a number of serine and threonine residues, among them S215 is located at a β-turn slightly protruding out of the membrane surface and is exposed to the cytosol (Fig. 6a). In addition, residue S240 is located at the immediately following β-turn at the cytoplasmic side and is in line with the downstream residues of S215, spacing on average 14.3 ± 1.7 Å between the two hydroxyl groups in our simulations (Fig. 6a). Most importantly, S240 was also found phosphorylated in vivo in the mitochondrial proteomics study 49 . The proximity of phospho-S240 from S215 and the extended backbone configuration of residues between them make it possible for GSK3β to bind to the cytoplasmic rim of VDAC1, recognizing phospho-S240 and phosphorylating S215.
Although the consensus sequence of GSK3β substrate peptides S/T-X-X-X-S/T(P) 57 would predict subsequent T211 phosphorylation if S215 is phosphorylated by GSK3β, the two residues are situated right at the two ends of the same β-turn (Fig. 6b). It is therefore unlikely that phospho-S215 can be the primer for GSK3β to phosphorylate T211 of VDAC1 due to incompatible backbone configuration. To date, we have yet to discover a literature reporting T211 phosphorylation of VDAC.
Other VDAC1 residues implicated experimentally in complex formation with HK. Several other VDAC1 residues have been previously identified as essential for interaction with HKI and  Fig. 5 Effects of HKII on wt-VDAC1 and S215E phosphomimetic mutant. a Current-voltage relationships for wt-VDAC1 (left) and S215E VDAC1 (right) in the absence (control; black lines) or the presence of HKII (+HKII; red lines), acquired during a −80 to +80 mV linear voltage ramp protocol. b Effects of HKII on VDAC1 conductance monitored over time (purple circles: wt-VDAC1; brown diamonds: S215E phosphomimetic mutant). Change in conductance in the presence of HKII was relative to that prior to the addition of HKII. Summary graphs are shown for both wt-VDAC1 (n = 8) and S215E VDAC1 (n = 6) in the presence of HKII. Representative steady-state current recordings for both wt-VDAC1 and S215E VDAC1 in the presence and absence of HKII are shown in Fig. S14. HKII isoforms 7,21,[31][32][33][34][35] . To examine if our modeled complexes provided any evidence in support of the role for these interactions, we analyzed all possible hydrogen-bond and salt-bridge interactions between VDAC1 and HKII formed at any point during the MD simulations of all five complexes (Figs. S6 and S8-S11). We captured several interaction pairs which involve previously identified VDAC1 residues: E189(VDAC1)-R373(HKII) in HKV1 (Fig. S6), D78(VDAC1)-K346(HKII) in HKV2 (Fig. S8), and D78(VDAC1)-K323(HKII) and T51 (VDAC1)-E317(HKII) in HKV3 (Fig. S9). E189Q and D78N mutations (independently) of VDAC1 have been reported to inhibit the protective effect of HKI overexpression against apoptotic cell death 21,32 . According to our simulations, these mutations might interfere with the complex formation between VDAC1 and HK, for example, by disrupting the interaction of E189(VDAC1) with H373 of HKI, which is likely protonated, or with R373 of HKII. Other interactions affected by these mutations include D78-K346 (K346 is conserved between HKI and HKII) or D78-R323 (K323 of HKII is replaced by an R323 in HKI) bonds. Phosphorylation of VDAC1 at residue T51 by GSK3β has been previously shown to dissociate HKII from the mitochondria 35 . This phosphorylation might destabilize the HKII/VDAC1 complex by disrupting the T51(VDAC1)-E317(HKII) bond due to the electronegative group introduced at T51. The involvement of several other VDAC1 residues in direct interactions with HKII (or HKI) cannot be explained by our modeled complexes. One possible explanation is that these mutations do not disrupt the complex formation by direct inhibition of HK/VDAC1 interactions, rather, through other cellular processes such as VDAC1 oligomerization.
Higher-order HKII and VDAC1 organization. While in this study, we have focused on the complex formation between one VDAC1 and one HKII, reportedly they can also form higherorder organizations. HKII has been reported to exist as a monomer or a tetramer; intriguingly, no intermediates, i.e. dimers and trimers, have been detected 58,59 . The mechanism underlying the formation of tetramers and their exact arrangement are not known. It has been speculated that an initial translocation of an HKII monomer to the OMM may encourage the binding of subsequent HKII monomers to form a tetrameric complex around VDAC1. Confounding the interaction between VDAC1 and a tetrameric complex of HKII is that VDAC1 itself may also oligomerize 60,61 . VDAC1 has been reported to exist as monomers, dimers, trimers, tetramers, and even hexamers 60,61 . Though a tetrameric HKII is speculated to interact with VDAC1 oligomers, the stoichiometry of interactions between HKII tetramers and VDAC1 monomers/oligomers remains unknown. It should also be noted that while translocation of HKII is associated with cell survival, oligomerization of VDAC1 can result in apoptosis 2 . Thus, the dynamic process regulating the interaction of tetrameric HKII with oligomeric VDAC1 still needs to be delineated. Nevertheless, our computational model presents the first steps in the translocation and binding of monomeric HKII to the OMM and its interaction with VDAC1.

Concluding remarks
Atomic structures of protein complexes are usually derived from experimental techniques such as X-ray crystallography and cryo-EM. In the case of HKII/VDAC1 complex formation, a priori requirement of membrane binding of HKII hinders structural determination in typical experimental techniques. Simulation of this process is also challenging because complex formation occurs on a time-scale (in the order of minutes based on current, indirect experimental observation [Fig. 5b]) not feasibly obtained by conventional or enhanced sampling MD simulation techniques. Therefore, we designed a hybrid modeling approach integrating MD and BD simulations and report a stable structural model for the HKII/VDAC1 complex, which we validated using complementary experiments at multiple levels. The complex partially blocks the permeation pathways of the VDAC1 channel, a key feature that was also observed in our electrophysiology measurement, strongly supporting the relevance of the derived complex structure. Furthermore, the binding interface of the complex is highlighted by the presence of a potential phosphorylation site S215, whose phosphorylation/phosphomimeticmutation results in inhibition of HKII/VDAC1 interactions both in our BD simulations and in our electrophysiology recordings. This inhibition appears to be mediated through disruption of the electrostatic attraction between H-anchor and cytosolic outer-rim of VDAC1. Though other phosphorylation sites on VDAC1 have been identified 2,52 and several of these sites have been implicated in cardiac and neurodegenerative diseases, their functional consequences have not been delineated. Our study shows that phosphorylation of a single residue, S215, can disrupt HKII interaction with VDAC1. This provides insight into how phosphorylation of specific sites on VDAC1 can impact its binding with HKII and potentially other proteins that translocate to mitochondria, for example, the pro-apoptotic protein BAX. Characterizing the complex molecular interactions between VDAC1 and anti-apoptotic and pro-apoptotic proteins will help delineate their roles in cell death and survival.

Methods
Membrane binding simulations. We used the X-ray structure of human HKII (PDB ID: 2NZT) as a starting point for our simulations. Glucose and glucose-6phosphate were removed from the crystal structure, and one monomer from the HKII dimer was selected for the simulations. The missing regions (residues 1-16, 98-104, 346, 404-405, 518-525, 546-552, 645-649, 914-917) from the crystal structure were modeled using MODELLER 62 . Among missing regions, residues 1-16 were predicted to be in an α-helical structure by the secondary structure prediction tool YASPIN 63 , while other missing regions were predicted to be unstructured. As residues 1-16 are a part of the H-anchor and known to be crucial for membrane-binding, we carefully modeled this part by using as a template the highly homologous isoform enzyme from rat, HKI (PDB ID: 1BG3), for which the H-anchor is fully resolved as an α-helix. Notably, the H-anchor of rat HKI shares high sequence similarity with human HKII (Fig. S16). Other missing regions were modeled as loops using conjugate gradients and MD with a simulated annealing approach in MODELLER 64 .
To reduce the system size in our membrane-binding simulations, we removed the C-domain and used only the N-domain (residues 1-463) of human HKII for membrane binding simulations. N-domain of human HKII, referred to as HKII-N hereafter is known to be sufficient for binding to the OMM 43 . In order to assign the protonation states of ionizable residues of HKII-N, pK a values were estimated using PROPKA 65,66 . Then, HKII-N was placed above a symmetric HMMM lipid bilayer composed of 108 phosphatidylcholines (PC), 58 phosphoethanolamines (PE), 26 phosphatidylinositols (PI), 4 phosphatidylserines (PS), and 2 phosphatidic acids (PA) lipid molecules in each leaflet. The lipid composition was chosen to approximately replicate the OMM of mammalian cells 67 . HKII-N was placed 5 Å above the membrane, where 5 Å is the minimum distance between any HKII-N atom and any membrane atom along the membrane normal (z-axis). Twelve independent HMMM systems were constructed using CHARMM-GUI 68 , each resulting in a box with dimensions of 125 × 125 × 174 Å 3 , including 240,000 atoms. To minimize possible bias in membrane binding simulations due to a particular initial lipid distribution, the lipid bilayers in all replicas were generated independently using CHARMM-GUI 68 to obtain a different (randomized) initial lipid distribution for each replica. All replicas were neutralized with 150 mM NaCl and solvated with the TIP3 water 69 . All the 12 solvated HMMM systems were then energy-minimized for 2500 steps before used for the subsequent 200 ns membrane binding simulations. A harmonic restraint on the z-position (along the membrane normal) was applied to the carbonyl atoms of short-tailed lipids with a force constant of k = 0.1 kcal mol −1 Å −2 , to mimic the atomic distributions of a fulltailed lipid bilayer and to prevent occasional diffusion into the aqueous phase 70,71 .
Conversion of HMMM lipids to full-tailed lipids. In order to generate a complete model of membrane-bound HKII-N, we used the last frame of one of the HMMM membrane-binding simulations described above and converted it to a fullmembrane construct using CHARMM-GUI 68 . Short-tailed HMMM lipids were transformed into full lipids by removing the DCLE molecules and adding the missing carbons to the lipid tails while preserving the positions of the headgroups and the initial six carbons of the lipid tails. The system was then energy-minimized for 10,000 steps and then equilibrated for 2 ns, while harmonic restraints (force constant, k = 1 kcal mol −1 Å −2 ) were applied to (1) the heavy atoms of the protein and (2) the positions of lipid heavy atoms corresponding to short-tailed lipids only in the cis leaflet of the lipid bilayer (the leaflet facing the protein). These restraints were applied to preserve lipid-protein interactions established during the association of the protein with the lipid membrane in HMMM simulations, while the newly added carbon atoms adjusted to the system 72 . Following this step, the system was simulated without restraints for 200 ns, resulting in a fully relaxed, membranebound protein.
Simulation of membrane-embedded VDAC1. The X-ray structure of mouse VDAC1 (PDB ID: 3EMN) 38 was embedded in a full membrane, with the same lipid composition as described above, using CHARMM-GUI 73 . Mouse and human VDAC1 proteins have a sequence identity of 94%. The system was neutralized with 150 mM NaCl and solvated with TIP3 water 69 resulting in a box size of 120 × 120 × 84 Å 3 including 110,000 atoms. Then the system was energy-minimized for 2500 steps and equilibrated for 1 ns, during which Cα atoms of the protein were harmonically restrained with a force constant k= 1 kcal mol −1 Å −2 . Following this step, the system was simulated without restraints for 200 ns.
MD simulation of HKII/VDAC1 complexes. For each of the complexes (HKV1, HKV2, HKV3, HKV4, and HKV5) derived from the BD simulations, a full-length HKII/VDAC1 model was generated by adding the C-domain of HKII to the coordinates of the HKII-N/VDAC1 complex (based on the known full structure of HKII). Each extension was carried out by adopting coordinates from full-length HKII after superimposing the backbone atoms of residues S449-L463 (an α-helix) into HKII-N. The resulting full-length models were then embedded in a fullmembrane with the same lipid composition as described above using CHARMM-GUI 73 . These five different systems were neutralized with 150 mM NaCl and solvated with TIP3 water 69 . HKV1 simulation system resulted in a box size of 215 × 215 × 160 Å 3 including 700,000 atoms; HKV2 system resulted in a box size of 120 × 120 × 240 Å 3 including 330,000 atoms; HKV3 system resulted in a box size of 160 × 180 × 140 Å 3 including 330,000 atoms; HKV4 system resulted in a box size of 120 × 120 × 230 Å 3 including 320,000 atoms; and HKV5 system resulted in a box size of 120 × 120 × 200 Å 3 including 300,000 atoms. Two independent replicas were created for each of the systems, resulting in a total of 10 independent simulation systems. Each system was energy minimized for 2500 steps and equilibrated for 1 ns while harmonically restraining the Cα atoms of both proteins with a force constant k = 1 kcal mol −1 Å −2 . Each system was then simulated for 650 ns, without restraints.
MD simulation protocols. MD simulations in this study were performed using NAMD 2.13 74,75 utilizing CHARMM36m 76 and CHARMM36 77 force field parameters for proteins and lipids, respectively. Bonded and short-range nonbonded interactions were calculated every 2 fs, and periodic boundary conditions were employed in all three dimensions. The particle mesh Ewald (PME) method 78 was used to calculate long-range electrostatic interactions every 4 fs with a grid density of 1 Å −3 . A force-based smoothing function was employed for pairwise nonbonded interactions at a distance of 10 Å with a cutoff of 12 Å. Pairs of atoms whose interactions were evaluated were searched and updated every 20 fs. A cutoff (13.5 Å) slightly longer than the nonbonded cutoff was applied to search for the interacting atom pairs. Constant pressure was maintained at a target of 1 atm using the Nosé-Hoover Langevin piston method 79,80 . Langevin dynamics maintained a constant temperature of 310 K with a damping coefficient, γ, of 0.5 ps −1 applied to all atoms. HMMM simulations were performed by employing a constant area in the XY dimension (membrane plane). For the full-membrane simulations, a constant ratio was used instead, which keeps the X:Y ratio of the unit cell constant. Simulation trajectories were collected every 10 ps.
Electric field simulations. Ionic current through VDAC1 was calculated by performing simulations of the membrane-embedded form of the channel in the presence of a constant electric field normal to the membrane. Simulations were performed for both HKII-free and HKII-bound states (HKV1) of VDAC1. The starting point for each simulation was the last snapshot of their respective fullmembrane equilibrium MD simulation. For each system, the salt concentration was increased from 150 to 500 mM by using the AUTOIONIZE plugin of VMD 81 , in order to enhance the number of ion permeation events during the electric field simulations. Three independent simulations were performed for each system, under a membrane potential difference of −50 mV, each for 80 ns.
Ionic current (I) was computed by counting the number of ions (Na + and Cl − ) that cross the VDAC1 channel (moving all the way from one side of the membrane to the other) over time, i.e. I = N q τ , where N is the number of crossing events over a time interval τ, and q is the charge of the ion (1.60217662 × 10 −19 Coulombs for Na + and −1.60217662 × 10 −19 Coulombs for Cl − ).
BD simulation setup. We used the membrane-bound HKII-N and membraneembedded VDAC1 to simulate the formation of their complexes using BD simulations. Atomic coordinates of both proteins were extracted from the final snapshot of their respective MD simulations in full-membranes. During the BD simulations, membrane-bound HKII-N was considered as the moving protein, while VDAC1 was the stationary protein, with both proteins modeled as rigid body entities.
Given the rigid representation of the molecular entities in BD simulations, we did not use an explicit membrane in these simulations. In order to take into account membrane placement of HKII-N and VDAC1 (as observed in the MD simulations of the membrane-bound forms of these proteins), in the beginning of each simulation, the z-coordinates of each protein were shifted so that their respective membrane midplanes would match (Fig. 3A). Another set of restraints were then used to maintain both the insertion position and orientation of HKII-N during the BD simulation; we applied restraints on the positions of H-anchor residues along the z-axis in our BD simulations. Restraints were applied on four residues (M1, F10, N15, Q20) of H-anchor using a grid-based potential coupled to each residue. Grid-based potentials were obtained by Boltzmann inversion of the probability distribution function of the z-coordinate of each residue obtained from the last 100 ns of membrane-bound MD simulation of HKII-N in full-membranes. Applied restraints were able to maintain the positioning and the tilt angle distribution of the H-anchor in BD simulations close to those observed in fullmembrane MD simulations (Fig. 3A) while allowing for the lateral diffusion of HKII-N around VDAC1 and their complex formation. To avoid unnecessary sampling of HKII-N during simulation, we employed a circular grid-based potential wall around VDAC1 with a radius of 150 Å.
It is important to note that BD simulations do not allow for internal relaxation of the interacting proteins and the membrane upon their interaction and binding. Therefore, BD here is used to only generate an initial, unrefined candidate structure for the complex purely based on their overall (global) attraction as rigid bodies. These BD-generated structures are then subjected to several hundreds of nanoseconds of MD simulations (details in the "MD simulation of HKII/VDAC1 complexes" section) during which the proteins and their embedding membrane are allowed to relax and detailed interfacial interactions can be refined.
BD simulation protocols. BD simulations were performed using a GPU accelerated BD code, ATOMIC RESOLUTION BROWNIAN DYNAMICS (ARBD) 82 . The masses and moments of inertia of the HKII-N (moving protein) were calculated directly from its atomic coordinates. HYDROPRO program suite 83 was used to estimate the translational and rotational friction coefficients of HKII-N protein which provided Langevin forces and torques at each timestep to keep the system at 310 K. To be noted, our calculated friction coefficients represent the diffusive behavior of HKII-N in aqueous solutions. Though this diffusive property would be significantly different for membrane-bound HKII-N, it will not affect the results of our study, as we are not focusing on either diffusive or kinetic properties of HKII-N binding to VDAC1.
The force and torque acting on HKII-N, due to stationary protein VDAC1, were calculated as follows. First, a grid of electric charge and Lennard-Jones (LJ) particle densities were obtained from the atomic coordinates of membrane-bound HKII-N with a 1 Å resolution. The density in each grid cell experienced a local force due to the corresponding grid-specified stationary VDAC1 potential (electrostatic or LJ terms). These local forces and the corresponding torques were summed over to obtain the total force and torque on the HKII-N molecule. A cutoff of 34 Å was used for the force calculation. Electric charge and LJ particle density of the atoms comprising HKII-N were calculated using the VOLMAP utility in VMD. On the other hand, stationary protein VDAC1 was represented by the electrostatic and LJ potential map of its atomic coordinates. The electrostatic potential map for the stationary protein (VDAC1) was calculated by solving the nonlinear Poisson-Boltzmann equation implemented in the APBS software 84,85 . cglen and fglen options in APBS were chosen to ensure a resolution of 1 Å for computing the cubic electrostatic potential grids. Dielectric constants for the protein interior and solvent were set to 12.0 and 78.0, respectively. Ionic radii were set as per the CHARMM36 force field 76 . A cubic-spline surface 86 model was implemented to model the dielectric interface and ion accessibility. The rate of dielectric transition was set to 0.3 and surface density to 10. The radius of the solvent molecule was set to 1.4 Å, the same as water. These parameters are comparable to previous studies of protein-protein interactions using ARBD simulations 87 .
The LJ parameters of the atoms comprising the stationary protein (VDAC1) were clustered into three categories: one representing all hydrogen atoms, another representing oxygen, and nitrogen atoms, and the final one representing carbon and sulfur atoms. The atoms in each category were assigned an average value for the parameters R min and ϵ. Then, ϵ was scaled down by 0.3, following a similar approach adopted by McGuffee et al. 88 to avoid the stickiness of proteins during the BD simulations. These two parameters, R min and scaled ϵ, were used to generate potential maps for the interaction of the three above-mentioned atom categories with the stationary protein, VDAC1, using the IMPLICIT LIGAND SAMPLING (ILS) implementation of VMD at 1 Å resolution 81,89 . ARBD simulations were performed using a timestep of 200 fs. Simulation coordinates were collected every 0.2 ns.
Generation of an H-anchor plugged model. An alternative model for the HK/ VDAC1 complex has been reported in the literature 36,37 . In this model, which we refer to as the "H-anchor plugged model", the H-anchor completely plugs the VDAC1 pore. To investigate the electrophysiological properties of such a model and how VDAC1 permeability is affected, we generated an approximate model of the H-anchor-plugged HKII/VDAC1 complex by docking the H-anchor into the lumen of VDAC1, which was then used in our electric field simulation for current measurements. Docking was performed using the "easy interface" implemented in HADDOCK2.2 web portal 90 . The X-ray structure of human HKII (PDB ID: 2NZT, missing regions modeled) and the same of mouse VDAC1 (PDB ID: 3EMN) 38 were used for docking. Docking was performed between the residues of H-anchor and the residues of VDAC1 facing the lumen on the cytoplasmic half resulting in multiple structural models. A representative structure was selected based on the best HADDOCK score, which was then embedded in a full-membrane with the same lipid composition as those used for our own model (described above) using CHARMM-GUI 73 . The system was neutralized with 150 mM of NaCl and solvated with TIP3 water 69 , resulting in a box size of 120 × 120 × 200 Å 3 including 270,000 atoms. The newly generated system was energy minimized for 2500 steps and equilibrated for 1 ns while harmonically restraining the Cα atoms of both proteins with a force constant of k = 1 kcal mol −1 Å −2 . Then, the system was simulated without restraints for 100 ns before using it for computational electrophysiological measurements.
Analysis. System visualization and analysis were carried out using VMD 81 . Depth of insertion of HKII-N into the membrane was assessed by calculating the z component of the center of mass (COM) of protein side-chain atoms relative to the average plane of the phosphorous atoms of the membrane. Interaction energy (van der Waals + electrostatic) between HKII-N and VDAC1 was calculated using the NAMDENERGY plugin of VMD. Contacts between HKII and VDAC1 are calculated using a 3 Å distance cutoff between any atoms of the two proteins. A hydrogen bond was counted to be formed between an electronegative atom with a hydrogen atom (H) covalently bound to it (the donor, D), and another electronegative atom (the acceptor, A), provided that the distance D-A is <3 Å and the angle D-H-A is more than 120°. Clustering was performed based on the position of the HKII-N with respect to VDAC1 (calculated using the root-mean-square deviation (RMSD) of HKII-N after superpositioning the BD trajectories using VDAC1). The "measure cluster" module implemented in VMD 81,91 was used for clustering analysis with an RMSD cutoff of 10 Å.
Experimental electrophysiology. Recombinant rat wild-type (wt) or mutant (S215E) VDAC1 proteins were reconstituted into planar lipid bilayers as described previously 48 with some modifications. Briefly, phosphatidylethanolamine (PE) and phosphatidylcholine (PC) (Avanti Polar Lipids) were mixed in a ratio of 7:3 (v/v), dried under nitrogen gas, and resuspended in n-decane (Sigma) for a final lipid concentration of 25 mg/ml. The cis/trans chambers contained symmetrical solutions of 10 mM Trizma base (Sigma), 500 mM KCl (Sigma), and 1 mM CaCl 2 (Sigma), pH 7.4. The cis chamber was held at virtual ground and the trans chamber was held at the command voltage. The pClamp software (version 10, Molecular Devices, San Jose, CA) was used for data acquisition. Currents were digitized at 5 kHz and low-pass filtered at 1 kHz using a voltage clamp amplifier (Axopatch 200B, Molecular Devices) via a digitizer (DigiData 1440A, Molecular Devices). The recombinant VDAC1 proteins were added into the cis chamber. Insertion of VDAC1 into the bilayer membrane and its function was monitored and confirmed by current recordings in response to a ramp protocol from −80 to +80 mV, as VDAC1 is characterized by its uniquely distinct current-voltage relationship. Subsequently, the solution in the cis chamber was replaced with the same initial buffer solution at the speed of 2.5 ml/min to remove non-inserted VDAC1 proteins and prevent additional VDAC1 insertion into the bilayer membrane. Current recordings under control conditions (in the absence of HKII) were initially taken, followed by addition of HKII (human recombinant HKII, 60 kU/ml; Genway Biotech, Slan Diego, CA) into the cis chamber. Currents were monitored during a 30-s recording duration every 2 min up to 24 min, and analyzed using CLAMPFIT (Molecular Devices) and ORIGIN (version 10; OriginLab, Northampton, MA). Mean current from each time point was normalized to its control as a percentage data, which was used for statistics later.
Significant differences between groups were determined with one-way ANOVA (SPSS Statistics 24, IBM) with posthoc significance analysis. P < 0.05 was considered significantly different.
Statistics and reproducibility. The data plotted in Figs. 2c, 3a, and 6b are the mean and standard deviation from more than 100 snapshots extracted from trajectories.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. The raw data used for the data diagrams in the main figures are available in https://doi.org/10.6084/m9.figshare.14511885.v1 92 .

Code availability
Simulation trajectories were collected using NAMD and ARBD. Visualization and analysis were performed using VMD and Python. All of these software packages are publicly available.