Nucleotide asymmetry and flexible linker dynamics modulate drug efflux cycle of P-glycoprotein, A computational study

TL;DR

This study uses computational methods to uncover how nucleotide asymmetry and flexible linker dynamics influence the drug efflux cycle of P-glycoprotein, a key player in multidrug resistance in cancer.

Contribution

The study reveals the dynamic roles of nucleotide asymmetry and linker flexibility in modulating P-glycoprotein's drug efflux mechanism.

Findings

Asymmetric nucleotide coordination at NBS correlates with TMD restructuring for substrate efflux.

The flexible linker transiently forms α-helices that affect NBD dimerization.

Conformation-dependent substrate pathways and TMD-linker interactions facilitate tunnel formation for drug transport.

Abstract

Despite advancements in oncology, multidrug resistance (MDR) mediated by P-glycoprotein (P-gp/ABCB1) remains a major barrier to chemotherapy. P-gp is an ATP-binding cassette transporter that undergoes nucleotide-driven structural rearrangements to efflux chemotherapeutics, but the mechanistic details of the substrate transport remain poorly resolved. Here, we performed high-throughput multi-replica molecular dynamics to simulate P-gp in a lipid bilayer (totaling ∼110 µs) to dissect nucleotide-dependent conformational changes across the transport cycle. Our adaptive sampling strategy reveals asymmetric nucleotide coordination at nucleotide-binding sites (NBS), which correlates with transmembrane domain (TMD) restructuring for substrate efflux. The experimentally unresolved flexible linker transiently forms up to five turns of α-helix that affects the nucleotide binding domain (NBD) dimerization process. We identified conformation-dependent substrate/allocrite pathways including nucleotide-specific access routes, while TMD-linker interaction facilitates substrate access tunnel formation. Together, these pathways reveal that the concerted interplay of nucleotide occupancy, linker dynamics, and overall protein conformation governs the structural plasticity and broad substrate promiscuity of the substrate binding cavity in P-gp. By integrating these findings, this work bridges static structural data with dynamic functional insights to further our understanding of the P-gp substrate translocation cycle.

Figures

Asymmetry of the nucleotide binding sites and preferential binding toward ATP and flexible linker affects nucleotide affinity. (a) Starting conformation of the four IF P-gp and two OF P-gp human homology models. The two transporter halves are colored in yellow and blue, and the flexible linker in red to highlight the differences between each P-gp starting configuration. Each of the models were incorporated with four nucleotide states; ADP/ADP, ADP/ATP, ATP/ADP, ATP/ATP. (b) Orientation of A-loop in IF P-gp is shown, where Y401 and Y1044 are involved in coordinating the nucleotide adenosine rin

Flexible linker forms a helix and forms stable interaction with TM3/9. (a) ICLs and intracellular portions of TM helices. ICL1 and ICL4 (cyan and purple) interact with the helical segment of the flexible linker via H-bonding interactions between polar side chains. (b) Center-of-mass distance between intracellular portions of TM3–4 and TM9–10 representing the distance between ICL1–4 and ICL2–3, respectively. (c) PCA heatmap with PC1/2 for four nucleotides state. PC1 represents the distance between the N-terminal side of the linker and the ICLs and linker helix. PC2 represents the distance betwe

Ligand accessible tunnels are shaped in a nucleotide dependent manner in IF P-gp for potential substrate ingress into the binding cavity. (a) View from the lipid bilayer cross section highlighting the ligand access tunnels identified by CAVER 3.0. P-gp orientation is set to show the channel orientations in the most optimal way. NBDs may appear to be close, but are separated as the panel shows IF P-gp conformation. Starting position of probes were residues F336 and F983 within the cavity. Tunnels are defined as the following: 1b-intracellular opening near TM4/6 (pink); 1c-intracellular opening

Potential substrate egress channels identified in OF P-gp showing high structural plasticity upon nucleotide incorporation. (a) Viewed from the lipid bilayer cross section as [Fig. 3](#fig0015). Tunnels are defined as the following: 5a- extracellular side bordering upper leaflet access route via TM1/3 (green); 5b- upper leaflet access route via TM1/3 (green); 4a- extracellular side bordering upper leaflet access route via TM7/9 (blue); 4b- upper leaflet access route via TM7/9 (blue); 6- extracellular space egress toward the solvent near TM6/12. (b) Alternative view of the ligand access tunnels

Full text

Introduction

1

ATP-binding cassette (ABC) transporters are integral membrane proteins that drive active transport of substrates against the thermodynamic gradient by harnessing energy from ATP molecules [1], [2]. This ATP-dependent transport mechanism enables ABC transporters to fulfill critical roles in diverse physiological processes, including detoxification, lipid trafficking, and immune response [3], [4]. In humans, overexpression of ABC transporters like P-glycoprotein (P-gp/ABCB1) compromises chemotherapeutic agents by reducing the intracellular concentration of drugs through active efflux. Understanding the structural dynamics and energetics of ABC transporters is crucial for developing adequate strategies to overcome MDR and improve therapeutic outcomes.

P-glycoprotein (P-gp) is the first drug resistance-related ABC transporter to be identified, characterized as a 170 kDa glycosylated protein capable of effluxing diverse chemotherapeutics [5]. The overexpression of P-gp in tumors such as myelogenous leukemia, liver cancer and lung cancer correlates with poor chemotherapeutic treatment outcomes [6], [7], as P-gp expels hydrophobic drugs (e.g., paclitaxel, doxorubicin) out to the extracellular space [8]. Beyond oncology, P-gp is ubiquitously expressed at blood-tissue barriers (e.g., liver, kidney, brain) and in immune cells, where it regulates the bioavailability of analgesics, antivirals, and other therapeutics, and contributes to endogenous lipid/steroid handling and barrier homeostasis [8]. Polyspecificity of P-gp to transport structurally unrelated substrates stems from a large and flexible substrate-binding cavity within the transmembrane domain (TMD). However, the lack of clear structure-activity relationships for P-gp substrates has hindered the design of selective inhibitors.

P-gp is encoded as a single polypeptide comprising two TMDs that each contains six α-helices and two nucleotide-binding domains (NBDs), with a flexible linker connecting NBD1 and TMD2 that is unresolved in available P-gp structures [2], although possible traces of the linker region were found via cryo-EM [9]. The TMDs form a hydrophobic cavity accessible to substrates from the intracellular milieu in the inward-facing (IF) state [10], [11]. During the transport cycle, the NBDs dimerize as P-gp adopts the outward-facing (OF) state, releasing the substrate into the extracellular space. Recently, a cryo-electron microscopy (cryo-EM) study revealed IF P-gp with an occluded NBD dimer interface, highlighting the conformational plasticity of P-gp undergoing IF-to-OF transition [12]. The NBDs harbor conserved motifs critical for ATP binding and hydrolysis: Walker A binds a phosphate group of the bound nucleotide via a conserved lysine (GxxGxGK(S/T); x = any residue); Walker B (hhhhDE; h=hydrophobic residue) contributes a conserved aspartate that coordinates Mg²⁺; ABC signature motif (LSGGQ) from the opposing NBD stabilizes ATP coordination via glutamate side chain that acts as a general base to deprotonate a water nucleophile for in-line attack on the γ-phosphate while NBDs are dimerized; adenosine-binding loop (A-loop) positions adenine ring of the nucleotide via π-stacking with a conserved tyrosine [13]. ATP binding ultimately leads to NBD dimerization, sandwiching ATP molecules between the Walker A and the signature motif to enable ATP hydrolysis, which is followed by ADP and inorganic phosphate (Pi) release that resets the transporter [2], [14]. Several cryo-EM studies indicate that P-gp predominantly resides in IF conformations with OF states being transient [12], [15], [16], yet the sequence of ATP-driven transitions linking these states remains unresolved. Because nucleotide chemistry is a primary determinant of NBS geometry and NBD dimer stability, our study focuses on the major turnover states by systematically sampling ATP- and ADP-bound combinations at NBS1/NBS2. We did not explicitly model ADP-Pi product states: on the sub-microsecond timescales accessible in our approach, inorganic phosphate exhibits varying residence times and heterogeneous coordination that are difficult to converge, and excluding Pi allows us to isolate γ-phosphate-dependent contacts that cleanly distinguish ATP from ADP. Targeted ADP-Pi simulations are an important future extension but lie beyond the present scope.

There are functional aspects of the cycle that remain unanswered; (i) how the occupancy of the two non-equivalent nucleotide-binding sites (NBSs) alternate during the transport cycle, (ii) how the unresolved flexible linker between NBD1 and TMD2 affects conformational transitions, and (iii) the sequence of structural rearrangements linking ATP hydrolysis to substrate efflux. To study different conformational states of P-gp, we constructed full-length homology models of human P-gp as starting points for the molecular dynamics (MD) simulation; mouse P-gp IF structures (PDB - 4Q9H/4M1M, referred as IF-narrow and IF-wide, respectively), human OF cryo-EM structure (PDB - 6C0V, OF-occluded), and bacterial OF crystal structure from Sav1866 (PDB - 2HYD, OF-open). We avoided bacterial IF templates with extreme NBD separation distance to minimize non-physiological bias while still including a plausible OF-open reference to sample a putative extreme that human P-gp may only transiently visit. In most simulations, human homology OF-open models rapidly compact toward OF-occluded geometries within 100 ns, consistent with OF-open being short-lived; thus, our conclusions do not rely on Sav1866-like persistence.

Combining high-throughput MD with an adaptive, multiple-replica sampling strategy, we explore the interplay between the nucleotide state and linker flexibility to gain structural and mechanistic understanding of the P-gp transport cycle, favoring many independent replicas over a single long trajectory to improve state-space coverage and increase convergence between samples.

Results

2

P-gp exhibits basal ATPase activity in the absence of substrates [17], indicating that nucleotide occupancy alone can bias the transporter along its conformational cycle. To elucidate how nucleotide states modulate structural dynamics across IF and OF states, we conducted an extensive series of multiple-replica molecular dynamics (MD) simulations of human P-gp in various nucleotide states (apo, ADP/ADP, ADP/ATP, ATP/ADP, ATP/ATP) (Table 1, Table 2); we refer to the symmetric nucleotide states (ADP/ADP and ATP/ATP) as “even” and asymmetric states (ADP/ATP and ATP/ADP) as “odd”. To dissect mechanistic links between nucleotide incorporation and conformational shifts in P-gp, we initiated simulations from two IF (IF-wide and IF-narrow; varying in NBD separation distance, see Fig. S1) and two OF (OF-open and OF-closed; varying in cavity access to extracellular solvent) conformations (Fig. 1a).Table 1. Total number and length of initial 10 ns multi-replica MD simulations performed in this study. All replicas were initialed from the same starting configuration with a random seed of initial velocities.Table 1. ConformationIFOFModel PDBwidenarrowopenclosedLinker modelmdlafmdlafmdlmdlNo. of nucleotide system555555No. of replica simulation120120120120120120Run time (ns)101010101010Total run time (μs)666666Table 2Total number and length of extended multi-replica MD simulations performed in this study. All replicas that displayed outlier structural deviations in any of the predefined subdomains of P-gp (see Methods) were extended to 100 ns from the initial 10 ns production run.Table 2. ConformationIFOFModel PDBwidenarrowopenclosedLinker modelmdlafmdlafmdlmdlNo. of nucleotide system555555No. of replica simulation282631272418Run time (ns)100100100100100100Total run time (μs)141315.513.5129Fig. 1Asymmetry of the nucleotide binding sites and preferential binding toward ATP and flexible linker affects nucleotide affinity. (a) Starting conformation of the four IF P-gp and two OF P-gp human homology models. The two transporter halves are colored in yellow and blue, and the flexible linker in red to highlight the differences between each P-gp starting configuration. Each of the models were incorporated with four nucleotide states; ADP/ADP, ADP/ATP, ATP/ADP, ATP/ATP. (b) Orientation of A-loop in IF P-gp is shown, where Y401 and Y1044 are involved in coordinating the nucleotide adenosine ring. Walker A motifs (cyan), signature motif (purple) and the respective A-loop residue pairs in NBS (magenta) are shown. (c) Box plots representing the respective distributions of the residue pairs in IF P-gp. (d) H-bond interaction profile between the bound nucleotide and NBS in IF P-gp. The relevant residues involved in the interaction are labeled, where the ICL1/3 and ICL2/4 are colored in cyan and purple, respectively. (e) Frequency for key polar interactions between bound nucleotides and NBS residues (with distance of <3.0 Å) in OF P-gp normalized over all replica trajectories and displayed as a heatmap.Fig. 1

Asymmetric nucleotide coordination

3

To probe NBS dynamics in P-gp, we surveyed residue-pair distances near the A-loop. These distances reflect conformational stability at NBS1 and NBS2 across IF and OF states under varying nucleotide conditions (Fig. 1c), as utilized in DEER spectroscopy measurements of mouse P-gp in multiple nucleotide states [18].

Inward-facing conformers exhibit NBS1 flexibility

3.1

In IF P-gp, the A-loop residue-pair distance at NBS1 fluctuated extensively (∼25–65 Å) with maximal variability observed in the IF-wide-mdl conformer, while NBS2 exhibited comparatively restricted dynamics (∼25–45 Å) (Fig. 1c). This asymmetry suggests NBS1 adopts greater structural plasticity when NBDs are separated. Free binding energies of nucleotides were calculated using MM/PBSA method [74]. Despite the variability of this protocol in measuring absolute free energy values for highly charged molecules like ATP and ADP, we can qualitative comparison between each condition. The nucleotide binding energy calculations across IF conformers revealed that ATP tended to stabilize both NBSs relative to ADP; however, at NBS1 this difference was marginal except in the ATP/ATP system (Fig. S2). Protein-nucleotide electrostatic interaction analysis highlighted stable coordination of nucleotide phosphate groups by the Walker A residues S434/S1077 in all IF states (Fig. S3). Notably, ATP bound at NBS1 in the IF-narrow conformers formed an additional polar contact with T435, while purine ring of ATP interacting with R905 in intracellular coupling-loop (ICL) 4 correlated with reduced A-loop distances, suggesting nucleotide positioning modulates NBS geometry. The symmetry-related contact at NBS2 involves R262 in ICL2 as the nucleotide purine rings in IF-wide and IF-narrow similarly engage in the hydrogen bonding interaction with the arginine side chain (Fig. S3). Upon transition to OF conformers these ICL-nucleotide interactions (R905 at NBS1 and R262 at NBS2) are lost, indicating state-dependent disengagement of the ICLs from the A-loop region.

Outward-facing conformers reveal nucleotide sensitivity in NBS2

3.2

In OF P-gp, nucleotide chemistry dictated the NBS dynamics. For "odd" nucleotide states (ADP/ATP or ATP/ADP), ATP-bound NBS exhibited closer configurations as displayed in Fig. 1c (median A-loop distances at NBS2: 15.5 Å for ATP bound versus 25 Å for ADP bound). Conversely, "even" nucleotide states (ADP/ADP or ATP/ATP) showed symmetric distances (∼18 Å) at both NBS. This asymmetry was exclusive to the OF-open conformer, with NBS2 displaying greater sensitivity between ADP and ATP in the “odd” nucleotide states (Δdistance for bound ATP and ADP: ∼10 Å in NBS2, versus Δdistance: ∼4 Å in NBS1). Importantly, the A-loop distance metric reports on local NBS tightness rather than global NBD-NBD separation, as the center-of-mass NBD-NBD distances remain indistinguishable between OF-open and OF-closed (Fig. S4), ruling out NBD dimer separation and supporting nucleotide-dependent local aperture changes in OF-open. The observation of asymmetric nucleotide occupancy impacting NBS conformation was further corroborated by larger entropic differences between ATP and ADP in binding energies at NBS2 (Fig. S2). The OF-closed conformers lacked nucleotide-dependent A-loop distance changes, suggesting the sampled conformation is nucleotide insensitive, possibly being a post-hydrolytic state.

Electrostatic interactions underlie nucleotide specificity

3.3

Protein-nucleotide interaction frequency analysis revealed distinct coordination patterns for ADP and ATP (Fig. 1d,e). Both nucleotides formed stable interaction with Walker A residues K433/K1076 and T435/T1078 in all OF states. In contrast, γ-phosphate of ATP formed additional polar contacts with S429/S1072 (Walker A), which was absent in ADP-bound OF states. This ATP-specific interaction was replaced in ADP bound NBS with a stable interaction between β-phosphate of the nucleotide and G430/G1073 backbone. In both OF states, ATP formed more stable electrostatic interaction with NBS residues than ADP, which is in agreement with consistently lower binding energies for ATP distributing −50 kcal/mol and lower, compared to −50 kcal/mol and higher for ADP (Fig. S2). NBS2 nucleotide specificity was further evident in OF-open conformers, where the ribose hydroxyl of ATP formed a hydrogen bond with Q535; the apparent loss of this contact with ADP is likely due to nucleotide reorientation due to missing γ-phosphate mediated contacts. By contrast, in OF-closed conformers the tighter NBS pocket preserves the Q535-ribose interaction with both ATP and ADP (Fig. 1e), consistent with the lower nucleotide sensitivity of OF-closed seen in the A-loop analysis. Notably, ADP bound at NBS2 in OF-open ATP/ADP state exhibited the weakest binding energy at around −110 kcal/mol (Fig. S2), suggesting ATP hydrolysis at NBS2 may initiate conformational transitions within the NBD dimer. However, the absence of significant NBD separation across OF states (Fig. S4) implies the full reversion to IF state requires extended timescales beyond the current simulations as expected.

Structural and functional dynamics of the flexible linker

4

The flexible linker connecting the two homologous halves of P-gp is unresolved in experimental structures due to inherent structural flexibility [10], [11], [19], [20]. To resolve the dynamic behavior of the linker, we modeled the starting configurations of the linker using two distinct approaches: [1] ab initio prediction (MODELLER, referred as “mdl”) and [2] AlphaFold2-derived secondary structure (referred as “af”) (Fig. S1).

Linker architecture and helix-dependent interactions

4.1

In the simulations, both models revealed that the residues 669–690 of the linker form a transient α-helical structure (Fig. 2c/d). In IF-af conformers, the helical propensity of the residues 669–690 was higher than IF-mdl conformers, likely due to the initial helical content in the starting configuration in IF-af models. Notably, we observed the linker residues 669–690 forming electrostatic interactions with TM helices 3, 9, and 10 (Fig. S5). Positively charged arginines in the linker (R669, R673, R680) formed salt bridges with the acidic TM3/10 residues, while negatively charged aspartates/glutamates in the linker (D679, E686, E689) interacted with basic TM9 residue side chains, adhering the linker in between opposing TMDs. This transient helix formation correlated with shorter ICL distances in IF-af conformers (TM3–4: ∼25–30 Å; TM9–10: ∼25–27 Å) compared to IF-mdl conformers (TM3–4: 27–37 Å; TM9–10: 28–35 Å) (Fig. 2b).Fig. 2. Flexible linker forms a helix and forms stable interaction with TM3/9. (a) ICLs and intracellular portions of TM helices. ICL1 and ICL4 (cyan and purple) interact with the helical segment of the flexible linker via H-bonding interactions between polar side chains. (b) Center-of-mass distance between intracellular portions of TM3–4 and TM9–10 representing the distance between ICL1–4 and ICL2–3, respectively. (c) PCA heatmap with PC1/2 for four nucleotides state. PC1 represents the distance between the N-terminal side of the linker and the ICLs and linker helix. PC2 represents the distance between the ICLs and linker helix. A region between −40–0 PC1 and −10 to −40 PC2 is marked with a box to highlight the difference across all PCA results. The representative structures of the linker helix-ICL region are shown for each of the six P-gp conformers in each nucleotide state. (d) Alpha helical propensity of the flexible linker (630−700) normalized over all replica trajectories shown as helicity heatmap for all six P-gp conformers in apo and four different nucleotide states. Red bar with asterisks represents the helical region (679−688) within the AlphaFold2 P-gp model, which has been the template for the linker configuration in IF-wide-af and IF-narrow-af homology models.Fig. 2

Conformational landscape of linker-TMD interface

4.2

Principal component analysis (PCA) of linker-TMD conformations, visualized via PyEMMA-derived kernel density estimation [21], revealed distinct nucleotide- and conformational-dependent populations (Fig. 2c). IF-mdl conformers, characterized by lower helical content (∼10–15 %) within the linker, occupied a broad conformational space in the PCA map (top-left quadrant), representing the disordered linker configuration. In contrast, IF-af conformers were closely distributed at the bottom center region of the heatmap, representing the microstate configuration with closer TMD-linker proximity stabilized by electrostatic interactions. Interestingly, IF-wide-af conformers with NBS1 bound to ADP (ADP/ADP, ADP/ATP) resulted in the linker helix formation closer to the N-terminus (Fig. 2d), which was represented as a unique microstate density in the PCA heatmap (Fig. 2c, annotated box). These observations suggest the global P-gp conformation and possibly nucleotide state can affect conformational variability in the linker.

Conformation-dependent substrate access pathways

5

Having established NBS asymmetry and linker-mediated NBD/TMD coupling, we next explored how these conformations shape drug-translocation tunnels connecting the substrate binding cavity and surfaces of the membrane leaflets and the solvent.

IF state: substrate accessible routes

5.1

CAVER analysis [22] of IF conformers revealed three ligand-accessible pathways connecting the TMD substrate-binding cavity to the intracellular lipid/solvent interface (Fig. 3): (i) tunnels between TM4 and TM6 (channels 2a, 2b) in the N-terminal half of P-gp, (ii) a lateral portal via TM10/12 (channel 2c) in the C-terminal half, and (iii) cytosolic-facing portals near the solvent interface (channels 1b, 1c). In the IF-wide-mdl conformer, the nucleotide state altered the ligand-accessible pathway availability and the separation distance between the two NBDs. Channel 2c appeared exclusively when NBS1 was bound to ADP (Fig. 3c), which correlated with higher NBD center-of-mass (CoM) distances above 50 Å (Fig. S4). Channel 2a was observed in conformers with ADP occupancy at NBS2, which had lower NBD separation distance below 50 Å. Channel 2b remained frequently open (15–44 % of frames), supporting experimental observations that TM4/6 form a common entry portal [11], [12]. In contrast, IF-wide-af formed channels 2a/2b only, independent of nucleotide state (Fig. 3d). The absence of TM10/12-associated tunnels (channel 2c) and consistent occurrence of channels 2a/2b toward the intracellular solvent and lower lipid leaflet interfaces suggest linker-mediated stabilization of binding cavity access via TM4/6. IF-narrow conformers lacked detectable tunnels with adequate frequency (all occurrences were < 10 % of frames) with overall lower NBD CoM distances than IF-wide (Fig. S4). This suggests the conformational microstate observed in IF-narrow conformers represents an IF-occluded state with inaccessible binding cavity and closer NBDs.Fig. 3. Ligand accessible tunnels are shaped in a nucleotide dependent manner in IF P-gp for potential substrate ingress into the binding cavity. (a) View from the lipid bilayer cross section highlighting the ligand access tunnels identified by CAVER 3.0. P-gp orientation is set to show the channel orientations in the most optimal way. NBDs may appear to be close, but are separated as the panel shows IF P-gp conformation. Starting position of probes were residues F336 and F983 within the cavity. Tunnels are defined as the following: 1b-intracellular opening near TM4/6 (pink); 1c-intracellular opening near TM10/12 (purple); 2a-intracellular side bordering lower leaflet access route via TM4/6; 2b-lower leaflet access route via TM4/6; 2c-lower leaflet access route via TM10/12; 3-inner sub-pocket near TM4/8. (b) Alternative view of the ligand access tunnels. (c) Nucleotide dependent ligand access gate formation in IF-wide-mdl. Along the gray lines representing lipid bilayer boundaries to solvent, relevant tunnel paths from all simulation replicas and corresponding tunnel occupancy are shown. Each occupancy value is normalized based on the total number of frames from the respective replica trajectory pool. TM4/6 (pink) forms routes to tunnels 2a and 2b and TM10/12 (purple) routes to tunnel 2c. Last 70 % of frames were used to analyze converged portions of the trajectories. Tunnels with occupancy values higher than 0.10 are depicted. (d) Nucleotide dependent ligand access gate formation in IF-wide-af. No tunnel above 0.10 occupancy were detected in IF-narrow P-gp conformers.Fig. 3

OF state: substrate egress routes

5.2

In OF conformers, three ligand tunnels were identified (Fig. 4): (i) TM1/3-bounded channels [5a, 5b], (ii) TM7/9-associated routes [4a, 4b], and (iii) a solvent-exposed portal above the binding pocket (channel 6). OF-open exhibited nucleotide-dependent formations of the substrate egress tunnels. The apo state showed limited occurrence of channel 5a (Fig. 4c), while ADP/ADP state showed tunnels toward the upper lipid leaflet and extracellular solvent interfaces via TM7/9 (channels 4a/4b) and the extracellular portal (channel 6). Notably, the ADP/ATP state produced the largest ensemble fraction of open egress tunnels in high frequency. The channels 5a/5b toward TM1/3 and channels 4a/4b toward TM7/9 appeared in ∼80 % of the frames, while the solvent-exposed channel 6 was detectable in 86 % of the frames. The extensive opening of the binding cavity in the ADP/ATP state correlated with significant rearrangement in the TMD, as shown by the elevated distances between TM3–11 (33.5 Å), TM6–12 (23 Å), and TM8–9 (31 Å) (Fig. 4e). This correlation hints that loss of the γ-phosphate at NBS1 may relax TMD packing, potentially lowering the barrier to substrate release. Conversely, OF-closed conformers, marked by tightly packed extracellular TM helices (e.g. TM5–11: 32 Å in ADP/ADP state), lacked detectable tunnels with adequate frequency (above 10 % of the frames), which suggests OF-closed microstate being a post-transport configuration (Fig. 4e).Fig. 4. Potential substrate egress channels identified in OF P-gp showing high structural plasticity upon nucleotide incorporation. (a) Viewed from the lipid bilayer cross section as Fig. 3. Tunnels are defined as the following: 5a- extracellular side bordering upper leaflet access route via TM1/3 (green); 5b- upper leaflet access route via TM1/3 (green); 4a- extracellular side bordering upper leaflet access route via TM7/9 (blue); 4b- upper leaflet access route via TM7/9 (blue); 6- extracellular space egress toward the solvent near TM6/12. (b) Alternative view of the ligand access tunnels. (c) Nucleotide dependent ligand access gate formation in OF-open P-gp configurations. Gray lines represent lipid bilayer boundaries to intra- and extracellular space. Relevant tunnel paths from all simulation replicas and corresponding tunnel occupancy are shown. TM1/3 (green) forms routes to tunnels 5a and 5b and TM7/9 (blue) routes to tunnel 4a and 4b. Last 70 % of frames were used to analyze converged portions of the trajectories. Tunnels with occupancy values higher than 0.10 are depicted. (d) Top view of the OF-occluded P-gp configurations. No tunnels above 0.10 occupancy were detected. (e) Extracellular residue pair distances. The two TMD halves, TMD1 and TMD2, are colored in yellow and blue, respectively. The residue pairs used in the box plots are depicted with connected dots in purple.Fig. 4

Discussion

6

High-throughput MD sampling framework

6.1

Recent advances in structural biology, including cryo-EM and deep learning-based structure prediction tools like AlphaFold, have dramatically expanded the availability of static snapshots of ABC transporters [23]. However, capturing transient intermediate conformational states essential for understanding functional dynamics can be guided using computational approaches such as MD simulations, which could bridge temporal and spatial resolution gaps between experimental techniques. Although catalytic turnover of P-gp is millisecond-scale, our aim here is not to traverse the full reaction cycle but to densely sample local equilibria around experimentally observed IF/OF conformations under specified nucleotide states, the microstates that underlie the larger experimentally detected conformations. Although our multi-replica MD framework is in principle well suited to follow species with variable residence time such as inorganic phosphate after ATP hydrolysis, in this work we focused on the ATP-bound and ADP-bound nucleotide conditions that represented the pre-hydrolysis and post-hydrolysis-like states. The explicit ADP and inorganic phosphate intermediate condition was not covered in this study due to the uncertainty in the starting position of these moieties.

The conventional MD (cMD) simulations of membrane proteins like P-gp often suffer from insufficient sampling due to the high conformational degrees-of-freedom as evidenced by poor structural convergence in 3 × 200 ns simulation of mouse P-gp [24] and persistent sampling deficits of membrane bound P-gp even after 500 ns [25]. Here, we implement an adaptive MD workflow combining 120 short, parallel simulations (initial 10 ns scout runs) followed by selective extension (to 100 ns) of replicas exhibiting outlier structural deviations (Table 1, Table 2). All replicas used identical systems, force fields, membrane compositions, and thermodynamic conditions (310 K, 1 atm, NPT), differing only by randomized initial velocities (Table 1). By prioritizing trajectory extension based on structural variance, similar to machine learning-guided latent space sampling approach [26], our methodology displays greater sampling of transient conformations compared to the cMD benchmark (Fig.S2) as this strategy provided broader conformational coverage than a single 1.8 µs trajectory, as quantified by principal component projections. Limited conformational transitions in the cMD benchmark suggests kinetic trapping of P-gp configuration in the initial metastable state due to the rugged free-energy landscape, as evidenced in other MD benchmark studies [27], [28]. Stabilization of PC projections in cMD indicate that the simulation reached local equilibration, while different basins are sampled across replicas in multi-replica approach (Fig. S6). Our strategy partially mitigates the stochastic noise through redundancy with the use of a high number of replica simulations and tailors simulations to explore nearby transient states via targeted extensions of structurally more variable trajectories. This high replica MD protocol could be utilized to investigate different metastable states of other ABC transporters and other complex protein systems. An important limitation of our approach, however, is that the individual replica trajectories reach hundreds of nanoseconds rather than the microsecond to millisecond timescale of slower transitions across the ligand transport cycle, so the present method resolves which conformational intermediates are accessible and how nucleotides engage the protein to affect the protein conformation.

To assess how well the sampled conformations matches existing experimental models, we overlapped the aggregate simulated ensemble and the conformational neighborhoods of multiple cryo-EM/X-ray structures (Fig. S7). In the PCA landscape, PC1 largely reports the NBD-NBD separation, and our ensembles populate the same basin that contains several experimental IF structures (e.g., 4M1M, 8PEE, 7ZK9, 4Q9H, 7OTI). By contrast, intermediate conformations between IF and OF are under-sampled in our runs, which we attribute to the linker configuration transiently positioning between the NBDs and impeding full closure (consistent with our linker analysis). We therefore expect that longer (μs-ms) trajectories and broader linker sampling will further connect these regions; nonetheless, the present data establish that our multi-replica strategy accesses experimentally relevant neighborhoods while revealing nucleotide- and linker-dependent heterogeneity.

To benchmark starting coordinates against more recent human cryo-EM models, we compared our pre-equilibrated starting configurations to structures in PDB accession codes 8GMG (IF) and 8SB8 (OF) [12] by the backbone root mean square deviation (RMSD) (Fig. S8). The IF-narrow model shows similar structural characteristics with the cryo-EM IF state (including the characteristic TM4/TM10 kinks) with RMSD of 2.82 Å, whereas the IF-wide model shows larger RMSD of 3.70 Å driven by wider NBD separation and modest TMD orientational differences. The OF-closed model matches the cryo-EM OF-occluded conformation with 2.02 Å RMSD, while the OF-open model diverges mainly at the extracellular TMD ends, consistent with our observation that OF-open often compacts toward OF-occluded within ∼100 ns and is likely short-lived in human P-gp (Fig. S8).

NBS asymmetry and mechanochemical coupling

6.2

For clarity, we use “NBS asymmetry” in three distinct senses: (i) site asymmetry, NBS1 versus NBS2 within the same nucleotide condition; (ii) occupancy asymmetry, “odd” vs “even” ATP/ADP combinations across the two sites; and (iii) coordination asymmetry, residue-level differences in how ATP vs ADP are engaged.

The asymmetric nucleotide coordination of NBS1 and NBS2, evidenced by divergent A-loop distances (Fig. 1c) and binding energies (Fig. S2), is consistent with the two-stroke ("power stroke") model of ABC transporters, where ATP hydrolysis at one NBS drives conformational transitions [29], [30]. This site/occupancy asymmetry is further corroborated by DEER spectroscopy studies [18], which revealed two distinct populations (45 Å and 60 Å) for distances between spin-labeled residues at NBS1 in mouse P-gp, compared to a single peak (∼45 Å) at NBS2. These experimental observations align with our simulations: In IF conformers, NBS1 exhibits higher flexibility (A-loop distance: 25–65 Å) compared to NBS2 (25–45 Å) (Fig. 1c). Culbertson and Liao [12] reported Walker-A-to-signature motif distances of ∼28–34 Å (NBS1) and ∼21–27 Å (NBS2) for apo and ATP/ATP IF states respectively, indicating a consistently tighter NBS2. Although our A-loop metric reports local pocket tightness rather than interior Walker/LSGGQ spacing, our IF simulations likewise show smaller, less variable A-loop distances at NBS2 (Fig. 1c), in agreement with those cryo-EM-derived trends [12]. For completeness, we also simulated Apo/ATP and Apo/ADP scenarios in IF P-gp (Fig. S7); their global geometries were similar to their two-nucleotide counterparts within our sampling window and thus are not discussed further. In OF conformers, heightened sensitivity to nucleotide anhydride chemistry in NBS2 emerges as ATP and ADP coordinates with uniquely different residues (Fig. 1e). For instance, the dimerized NBDs enable full coordination of γ-phosphate of ATP by the signature motif (Q535) and Walker A (S1072) residues, while these interactions are absent in ADP bound NBS2. This NBS coordination asymmetry in OF P-gp, absent in IF states, may contribute to the mechanochemical coupling between ATP hydrolysis and TMD restructuring during substrate efflux that is unique to OF conformations.

The functional divergence between NBS1 and NBS2 is underscored by mutational studies. For instance, simultaneous mutation of Walker A lysines (K433/K1076) abolishes P-gp activity, while single mutants retain partial ATPase and transport function [31]. These findings align with our observation that K433/K1076 forms a stable electrostatic interaction when γ-phosphate is present (ATP bound NBS), whereas those contacts are weakened or lost in ADP-bound states lacking the γ-phosphate (Fig. S3). How the loss of electrostatic contacts mediated by γ-phosphate upon ATP hydrolysis affects local NBS1/2 geometry, overall NBD dimer and P-gp conformations could be studied incorporating ADP-Pi bound states. The additional electrostatic interactions within NBS established in the presence of ATP is likely related to the thermodynamic preference for ATP over ADP. This idea is supported by free energy calculations of P-gp NBD dimers, where ATP bound states exhibit lower potential of mean force profiles than ADP-Pi bound states [32]. Because hydrolysis yields ADP and Pi, transient ADP-Pi occupancy can reshape NBS electrostatics. We deliberately modeled ATP and ADP states without explicit Pi to isolate the contribution of the γ-phosphate and to avoid poorly converged Pi residence times on the ns timescales sampled. Given our results, we would predict that explicit Pi, before egressing out of NBS, would partially reinstate γ-phosphate-like contacts at NBS2 (e.g., with S429/S1072 and K433/K1076), reduce local A-loop separations in OF-open, and stabilize OF-closed lifetimes, mirroring the states seen in vanadate-trapped preparations. The eventual exit of Pi would destabilize the local chemistry in NBS to facilitate further separation between ADP bound Walker A and opposing signature motif. Testing these predictions with longer simulations that include ADP-Pi will be an important next step but does not alter the conclusion that γ-phosphate-dependent coordination is accentuated at NBS2 in OF states.

The NBS coordination asymmetry in ATP hydrolysis is evidenced by biochemical studies of Walker B human P-gp mutants (D555N/D1200N) [33]. Using vanadate-trapping combined with [α-^32 P^]-8-azido-ATP photolabeling, the authors showed that under turnover conditions, the C-terminal half of P-gp (NBS2) was preferably driven into the post-hydrolysis state unlike the N-terminal site (NBS1) remaining predominantly in an ATP-bound state, suggesting ATP is hydrolyzed preferably at NBS2 than NBS1. This parallels our observation of nucleotide-specific coordination changes for NBS2 in OF-open conformers, and it also rationalizes why product-like states at NBS2 (ADP/ADP or ATP/ADP) show looser local A-loop closure relative to ATP bound states. While ATP binding stabilizes NBS2 via γ-phosphate coordination (Fig. 1e), ADP occupancy correlates with elevated A-loop distances (25 Å and 15.5 Å for ADP and ATP, respectively). The site/coordination asymmetry in NBS is critical for the "alternating access" mechanism [34], where iterative ATP hydrolysis at one NBS drives TMD rearrangements. However, the absence of significant NBD separation in OF states (Fig. S4) suggests the full IF→OF→IF transitions require prolonged timescales or substrate binding. This is consistent with MD studies of Sav1866 showing stable NBD dimerization in ATP bound states but progressive dissociation in ADP bound systems [35], [36].

Taken together with earlier spectroscopy and mutagenesis studies, our simulations add a concrete detail: in the inward-facing ensemble the more flexible NBS1 favors possible nucleotide exchange, whereas full γ-phosphate coordination at NBS2 emerges only after the transporter samples the outward-facing state. The catalytic readiness achieved with full coordination between NBS residues and γ-phosphate is then relieved following hydrolysis/product formation upon ATP hydrolysis. This site-specific difference in ATP versus ADP handling was not captured in the published experiments and refines current descriptions of the alternating-access cycle.

The flexible linker as a potential conformational gatekeeper

6.3

Our findings position the unresolved flexible linker as an important structural component of the P-gp conformational cycle. The partially helical linker in IF-af conformers stabilizes the intracellular interface between opposing TMDs (Fig. 2d) via electrostatic networks. The salt bridge formation between the charged linker residues (R669/D679) and TM3/9/10 (Fig. S5) effectively acts as a "molecular glue", modulating the NBD-NBD distance (Fig. S4) and TMD flexibility as evidenced by linker conformer dependent substrate tunnel formations (Fig. 3c/d). Consistent with this role, cryo-EM studies have shown that drug binding within the central cavity induces a pronounced shift of TM9, described as an “initiator of the peristaltic extrusion” mechanism [38]. Consistent with emergent helicity in our IF ensembles, HDX-MS on murine P-gp reported markedly incomplete deuterium uptake across three linker peptides (residues 619–684), indicating local secondary structure or persistent contacts despite the linker being unresolved in static structures [39]. In another study, persistent helices introduced by insertion within the linker inactivated P-gp [37]. Although the insertion mutation has been recently predicted to form a much longer helix than expected in a recent computational study [40], we view this as consistent with our interpretation that any secondary structure in the native linker is transient and context dependent rather than constitutive. Our simulations further indicate that helix formation within the linker correlates with reduced inter-ICL distances (Fig. 2b). The stable interaction between the linker and TMDs could impede premature NBD dimerization and prolong substrate-binding cavity accessibility, forming an IF occluded microstate. The IF occluded-state with the collapsed binding cavity and close NBD-NBD distance has been previously observed in a bacterial exporter McjD and human P-gp [12], [41]. In ATP/ATP IF ensembles, we observe partial compaction at the NBS surfaces even without bound drug, driven primarily by γ-phosphate-dependent electrostatics that bridge Walker A and the opposing LSGGQ signature motif (Fig. 1e, Fig. S3). Countervailing linker-TMD tethers (R669/D679 to TM3/9/10) and the absence of substrate-induced TMD preorganization limit full dimer closure, rationalizing the modest basal ATPase and the prevalence of IF-wide/ATP-bound cryo-EM states. Thus, ATP binding biases the ensemble toward dimer engagement, whereas completion of tight closure is tempered by linker contacts and the transmission interface until the substrate further reorganizes the TMDs.

Interpretation of the variable helicity in the linker must, however, account for potential bias from the initial model configurations deployed in this study. The linker-mdl construct was introduced without a helical portion in the linker, whereas the linker-af model contained a pre-formed α-helix (residues 679–688) predicted by AlphaFold2 (Fig. S1). Despite these differences in the initial configuration, emergent helicity developed in IF-mdl simulations (∼20–70 % helical propensity in residues 676–682 in IF-wide-mdl; ∼10–70 % helicity in 680–691 in IF-narrow-mdl), and IF-af simulations exhibited dynamic helicity across a broader segment (∼10–100 % helicity in residues 671–684 in IF-wide-af; ∼10–100 % helicity in residues 675–688 in IF-narrow-af). Notably, despite starting with linker with no helicity in the initial configuration, IF-wide-mdl trajectories displayed spontaneous helix formation in residues 679–688 mirroring the AlphaFold2-predicted region (Fig. 2d). IF-narrow-af conformers showed positional variability of the helix, despite the intrinsic α-helix-stabilizing propensity of the AMBER force field [42]. In IF-wide-af conformers, the helical segment shifted toward the N-terminus in ADP/ADP and ADP/ATP states, showing that the nucleotide state can transiently affect the linker conformation. Collectively, these observations indicate that linker helicity is not a modeling artifact, but rather a dynamic feature that reorients in response to the global conformational landscape of P-gp and is modulated by the nucleotide state. A recent report of multi-microsecond simulations that varied starting linker conformations similarly reported secondary structure formation predicted in the central region and persistent α-helical element across the linker in simulations, mirroring our observations [40]. Functionally, helicity appears as short and transient helices spanning a few helical turns that recur across replicas and nucleotide states (Fig. 2d), supporting transient secondary structure in this region.

The regulatory role of the linker has been hypothesized due to the presence of conserved charged residues in this region across ABC transporters [43]. Previous biochemical studies postulated a potential regulatory role of the linker via phosphorylation [37], [44], [45], although the results are contrasting [46]. Phosphorylation of serine residues (S661, S667, S671, S675, S683) within the helical linker region (residues 671–691) of P-gp could modulate its helicity, potentially fine-tuning transporter activity [44], [47]. Interestingly, alanine mutations of these serines do not abolish transport [44], suggesting the change in electrostatic profile of the linker does not critically affect the dimerization kinetics. Functional studies of linker truncations by Hrycyna and colleagues further support the structural role of the linker [37]. Removal of residues 653–686 (Δ34) abolished ATPase activity in human P-gp, while reintroducing a 17 residues long glycine/serine-rich peptide into the Δ34 shortened linker partially restored function [37]. This implies that the linker plays a critical role for the structural stability and conformational dynamics of P-gp. Complementing truncation data, proteolytic cleavage of the linker increased basal and differentially drug-stimulated ATPase activity and was interpreted as partial uncoupling of hydrolysis from transport [48]. These findings reinforce our proposal that an intact, electrostatically engaged linker can act as a gatekeeper that restrains premature NBD closure. From our simulation data, we observe transient linker helicity (residues 673–690) in IF conformers (Fig. 2d) and frequent salt-bridge formation between the charged linker residues in the helical portion and the intracellular segments of TMDs (Fig. S5). We predict that the charged residues in the native linker might introduce a kinetic barrier for NBD dimerization via electrostatic tethering, affecting the conformational transition kinetics which is consistent with the ATPase activation observed upon linker proteolysis [48]. In this framework, proteolysis or charge perturbations may loosen the linker-TMD tethers and lower the barrier to NBD-NBD dimerization, explaining the increased basal ATPase after linker cleavage [48]. MD simulations of murine P-gp corroborate this idea, showing the linker dampens NBD fluctuations and solvent accessibility [49], [50]. However, prior studies failed to observe persistent helical content in the linker due to limited sampling [49], [50]. Despite our efforts to enhance conformational sampling with the use of a multi-replica approach, intrinsic flexibility of the linker and the lack of inherent secondary structure mean our simulations do not exhaustively cover the full dynamic landscape of the linker, as evidenced by PCA heatmaps (Fig. 2c). Nevertheless, the consistent emergence of linker helicity across models and nucleotide states underscores the role of the linker as a potential conformational buffer, transiently modulating TMD/NBD coupling via electrostatic tethering, when localized in between the two wings of P-gp. With this insight, stabilizing the transient helical linker-TM3/9 interaction via allosteric modulator could prolong the IF-occluded intermediate state, which may slow transport activity and enhance chemotherapeutic retention.

If stabilization of the linker helix indeed steadies the inward-facing state, progression through the transport cycle would require this segment to relax or unwind so that the NBDs can finally dimerize. Such short-lived structural shifts lie outside the reach of current structural methods, thus future computational studies can shed light on this potential linker maneuver and point to experiments that test their role.

Dynamic substrate entry and release pathways

6.4

In IF P-gp, we identify formations of ligand accessible tunnels involving TM4/6 and TM10/12 (Fig. 3c), showing the flexible nature of the substrate binding pocket in P-gp to accommodate diverse substrates. For clarity, channel labels denote topology and exit direction: 1b/1c, cytosolic-facing portals near solvent; 2a/2b, lateral tunnels between TM4-TM6 (N-terminal half); 2c, a lateral portal between TM10-TM12 (C-terminal half); 4a/4b, upper leaflet and extracellular space directed routes between TM7-TM9; 5a/b, a TM1-TM3 route toward the extracellular solvent and upper leaflet; and 6, the solvent-exposed extracellular outlet (Fig. 3, Fig. 4).

The involvement of TM4 and TM10 in the formation of ligand access tunnels has been observed in cryo-EM studies of human P-gp bound to substrate/inhibitors, which showed high structural plasticity in TM4/10 that are actively engaged in trapping the drug in the binding cavity [10], [11]. Our simulations add an additional dimension to the ligand binding dynamics, as we observe the partially helical linker conformation stabilizing ligand accessibility to the central cavity via TM4/6 in IF-wide-af conformers (Fig. 3d), suggesting a possible facilitatory mechanism of the linker to aid substrate entry during nucleotide exchange. Although, the linker-TMD interaction and the respective binding cavity rearrangement are likely transient the occlusion of binding cavity near TM4/6 with "molecular plugs" like zosuquidar and elacridar shows high efficacy in P-gp inhibition [51], [52], [53]. Nucleotide-dependent TM10/12 tunnel formation in IF-wide-mdl underscores the dynamic nature of substrate recruitment in the presence of flexible linker configurations. ATP or ADP occupancy in NBS correlates to differential occurrence of ligand accessible channels from the central cavity, including the channel 2c that appears predominantly when NBS1 is bound to ADP. A possible explanation is that loss of the γ-phosphate at NBS1 weakens ATP-specific contacts between the Walker A-signature motif and ICL2/ICL3 transmission surface, biasing local loosening on the TMD2 side and transiently widening the TM10/12 lateral portal. This interpretation is consistent with our A-loop distance shifts and electrostatic interaction patterns (Fig. 1c/e,3c).

As expected from OF rearrangements, the IF tunnels we detect do not persist in OF ensembles: TM4, TM6, TM10, and TM12 are rotated to occlude lateral portals when NBDs are dimerized, in agreement with cryo-EM observations that ligand capture occurs in IF prior to tight dimerization [38]. Asymmetric nucleotide occupancy in OF-open bound to ADP/ATP shows the highest ensemble population of open efflux channels, with solvent-exposed portal (channel 6) and tunnels toward upper lipid leaflet and solvent interfaces (channels 4a, 4b, 5a) observed in most of the frames (Fig. 4c). Recent substrate-tracking cryo-EM with covalently tethered ligands further identified TM1 as a dynamic regulator of transport, including a mid-transport break near Gly72 that modulates gating [54]. Consistent with this, our CAVER analysis detects a TM1/3-directed egress route (channel 5a). TM1/2 loop (TM1-TM2) segments exhibit high positional and secondary-structure variability across our ensembles, features expected to influence egress timing and path selection. ADP/ADP state exhibits transient cavity opening toward the extracellular solvent (channel 6: 14 % occupancy) and upper leaflet via TM7/9 (channels 4a,4b: 13–19 % occupancy), while the TM1/3 side of the cavity remains inaccessible. The occurrence of substrate egress tunnels in partially or fully hydrolyzed nucleotide states (ADP/ATP and ADP/ADP states, respectively) in our simulations aligns with the "ATP switch" model [55], where ATP hydrolysis at NBS drives conformational resetting. ATP/ATP and ATP/ADP states in OF-open results in no detectable egress tunnels, as the distance between the extracellular ends of TM6 and TM12 is lower (TM6–12: ∼16 Å) in the respective nucleotides (Fig. 4e) leading to tighter packing of TM helices and binding cavity. These findings align with MD studies of Sav1866, where ATP/ATP states collapsed the substrate cavity while asymmetric nucleotide states led to opening of the extracellular gate [56], [57].

Unlike the OF-open ensembles, our OF-closed trajectories show no detectable efflux tunnels and retain a tightly packed TM6–TM12 interface (Fig. 4e). The OF-closed model was initiated from the ATP-bound E556Q/E1201Q (“EQ”) mutant, a construct that cryo-EM and biochemical probes have assigned to an outward-occluded, post-hydrolytic state [20], [58]. The persistence of the sealed cavity and dimerized NBDs in our simulations, independent of the nucleotide states, supports the interpretation that our OF-closed conformer represents a post-transport intermediate in which the export pathway has already shut, and the protein awaits nucleotide release to reset. Together with the lack of IF-type lateral portals in OF-open/OF-closed, this complementarity supports a model in which substrate access is favored in IF, while efflux proceeds via OF-only pathways (channels 4,5,6), reconciling our simulations with inhibitor-bound and substrate-bound cryo-EM states [38], [51].

The asymmetric role of NBS1 and NBS2 as observed in tunnel formation from the OF-open conformers is supported by biochemical studies. Hrycyna and colleagues proposed that each hydrolysis event serves distinct functions [33]. ³¹P SSNMR on BmrA trapped in various nucleotide states showed one NBD occupies a high-energy post-hydrolysis conformation, while the other NBD holds a non-hydrolyzed ATP, directly reporting on asymmetric catalytic cycling [59]. H/D exchange data by Vigano et al. [60] showed when E552Q and E1197Q mouse P-gp mutants are trapped in the post-hydrolysis state with ADP-Vi then rescued by adding ATP, the H/D exchange rate was increased from ∼45 to ∼80 % for E552Q P-gp, but no significant change was observed in E1197Q mutant [60]. With this data, the authors suggested hydrolysis at NBS1 drives major membrane-domain rearrangements. Based on the observation of substrate egress tunnels in OF-open conformers, we propose that ATP hydrolysis at NBS1 initiates substrate extrusion, while NBS2 hydrolysis resets the transporter.

Conclusions

7

This study presents a high-throughput MD framework that efficiently characterizes metastable states in P-gp, partially overcoming sampling limitations of conventional simulations. By deploying parallel short replicas and selectively extending trajectories with high structural variance, we capture transient intermediates critical to the substrate transport cycle of P-gp. Our simulations show a stage-wise model in which (i) NBS asymmetry biases inward facing ensembles, (ii) transient linker helicity modulates the inward facing to occluded state transition, and (iii) combined effects reshape substrate tunnels predicted for ligand entry and release. Together, our simulations link nucleotide chemistry, linker flexibility and substrate pathway to provide more detailed picture of how P-gp moves its cargo across the membrane. Beyond P-gp, our workflow provides a paradigm for probing metastable dynamics in other ABC transporters and flexible membrane proteins in general, bridging gaps between experimental insights and functional mechanism elucidation.

Materials and methods

8

Human P-glycoprotein 3D structural modeling

8.1

Crystal structures of IF murine and mouse P-gp (PDB: 4Q9H, 4M1M respectively) [10], [11], OF Sav1866 (PDB: 2HYD) [19] and cryo-EM structure of OF human P-gp (PDB: 6C0V) [20] were selected as templates for building full length human P-gp homology models. All sequences in template structures were mapped to human P-gp sequence (Uniprot ID: P08183) and missing residues were modeled using Modeller v9.23 [61]. For the IF models, a total of 5000 models of flexible linker (residues 630–699) were generated and a conformation with the highest Discrete Optimized Protein Energy (DOPE) was selected as starting configuration (referred as linker-mdl, see Fig. S1). The second starting configuration of the flexible linker was constructed by adopting the AlphaFold2 predicted structure (AF-P08183-F1-v4) [23]. The linker structure for OF P-gp homology models were generated, as described for linker-mdl.

MD simulation set up

8.2

Protein was embedded in 5:1 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC):cholesterol bilayer using CHARMM-GUI [62] based on the protein-membrane orientation predicted by OPM database [63]. Additional minimization was performed during system generation using CHARMM-GUI to relax the membrane lipids [62]. While a coarse-grained pre-equilibration route could further homogenize the bilayer, the present atomistic protocol yielded stable protein–membrane engagement across replicas and was sufficient for the purposes of this study. ATP-Mg^2+^ were docked by aligning nucleotide bound structure of ABCB10 (PDB: 4AYT) [64]. Simulations were performed using AMBER ff14SB force field for protein [65] and LIPID14 FF [66] for POPC-cholesterol membrane, and phosphate and magnesium parameters from Meagher et al. [67] and Allnér et al. [68]. The periodic TIP3P water box was used with Na^+^ and Cl^-^ ion concentrations of 150 mM. System was energy minimized using AMBER20 [69], applying harmonic restraints with a force constant of 1000–0 kcal/mol Å^2^ on the heavy atoms of protein, as described in other membrane simulation protocols [70].

GROMACS [71] 2021.5 package was used for equilibration and production runs. NPT ensemble was used during equilibration with semi-isotropic coupling at the constant pressure of 1 atm and constant temperature of 310 K, using Nosé-Hoover Langevin piston and Langevin thermostat. Time step of 1 fs was used for the initial 2.5 ns equilibration run consisting of 250 ps runs with harmonic restraints of 100, 50, 25, 10, 5, 2.5, 1, 0.5, 0.25, 0 kcal/mol Å^2^, followed by 12.5 ns run without harmonic restraints. Production runs were performed with an increased time step of 2 fs under the same ensemble. The electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method and all bonds to hydrogen atoms were constrained using the SHAKE algorithm. For each P-gp nucleotide system, 120 replicas of 10 ns production runs were performed starting from randomly assigned initial velocities, where the RMSDs of fourteen segments, including twelve TM helices and two NBDs, were computed with respect to the first frame. Among the 120 replicas, domain-specific RMSD outlier detection (Q3 + 1.5 ×IQR) identified replicas exploring under-sampled substates, which were extended for 90 ns. Total production run from all replicas of P-gp conformers and nucleotide states was ∼110 μs in aggregate (Table 1, Table 2).

Analysis of flexible linker interaction

8.3

To conduct principal component (PC) analysis, all replica trajectories were featurized with pyEMMA [21] using the closest heavy atom distances between the residues 159–189, 630–700, 801–831 to include intracellular portions of TM3,9,10 and the flexible linker. The percentage contributions of the top five PCs were 77.15, 6.27, 3.29, 1.73, 1.35, respectively. The first two PCs were used to generate 2D histograms, which were constructed as heatmaps using kernel density estimation. Residue-wise secondary structure analysis was conducted using the DSSP algorithm [72] and CPPTRAJ [73].

Binding-free energy calculations

8.4

Binding free energies for nucleotides (ATP/ADP) were estimated by the molecular-mechanics/Poisson-Boltzmann surface-area (MM/PBSA) protocol in 0.15 M salt concentration [74]. Prior to the energy decomposition, membrane lipids, ions and water in the original system were stripped so that only protein atoms were considered as receptor. Fifty uncorrelated snapshots were extracted from each 100 ns trajectory.

CAVER tunnel analysis

8.5

The CAVER 3.0 software [22] was used to identify the possible ligand access tunnels in the simulation trajectories with frames saved at intervals of 200 ps. Initial 30 % of the frames of each of the trajectories were stripped to minimize starting conformation bias. All frames of the trajectories belonging to the unique P-gp conformer and nucleotide state were combined and aligned to an initial frame. For all CAVER analyses, tunnel calculation was performed excluding the residues 1–44, 370–710, 1013–1280 to only include the TMD region. The starting points were defined as the benzene carbons in F336 and F983 to represent the canonical ligand binding site. CAVER tunnel analysis was run using a probe radius of 2.5 Å, shell radius of 6.0 Å, and shell-depth of 4.0 Å.

Author contributions

SP, HF and JW jointly conceived and supervised the study. SBH designed the study, performed calculations, analyzed the data, and wrote the manuscript. SBH, SP, HF and JW jointly edited the manuscript.

CRediT authorship contribution statement

Sungho Bosco Han: Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation. Stephen M. Prince: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. Jim Warwicker: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. Hao Fan: Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization.

Declaration of Competing Interest

There are no conflicts to declare.