Characterizing the binding interactions between P-glycoprotein and eight known cardiovascular transport substrates.

The multidrug efflux pump P-glycoprotein (Pgp) is upregulated in cardiomyocytes following chronic ischemia from infarction and hypoxia caused by sleep apnea. This report summarizes the molecular dynamic studies performed on eight cardiovascular drugs to determine their corresponding binding sites on mouse Pgp. Selected Pgp transport ligands include: Amiodarone, Bepridil, Diltiazem, Dipyridamole, Nicardipine, Nifedipine, Propranolol, and Quinidine. Extensive molecular dynamic equilibration simulations were performed to determine drug docking interactions. Distinct binding sites were not observed, but rather a binding belt was seen with multiple residues playing a role in each studied drug's stable docking. Three key drug-protein interactions were identified: hydrogen bonding, hydrophobic packing, and the formation of a "cage" of aromatic residues around the drug. After drug stabilization, water molecules were observed to leak into the binding belt and condense around the drug. Water influx into the binding domain of Pgp may play a role in catalytic transition and drug expulsion. The cytoplasmic recruitment theory was also tested, and the drugs were observed to interact with conserved loops of residues with a strong affinity. A free energy change of astronomical value is required to recruit the drug from the cytoplasm to the binding belt within the transmembrane domain of Pgp.

Introduction

One of the biggest thorns in modern medicine is acquired drug resistance, which is common to cancer transformed cells. Acquired drug resistance is the primary cause of chemotherapy failure in over 90% of patients with metastatic cancer (Longley and Johnston 2005). A primary factor contributing to cancer drug resistance is the overexpression of the MDR1 gene, which codes for the membrane protein P-glycoprotein (Pgp) (Ueda et al. 1987). Pgp is a 170 kDa transmembrane glycoprotein with approximately 10–15 kDa of glycosylation (Aller et al. 2009). It is a member of the ATP-binding cassette (ABC) group and has two conserved ATP-binding domains on the cytoplasmic side of the protein. Pgp is a 1280-residue, single polypeptide that contains two transmembrane domains (TMDs) and two nucleotide-binding domains (NBDs) with a TMD1–NBD1–TMD2–NBD2 topology (Wise 2012). ATP hydrolysis is stimulated in the presence of drugs that are transported, indicating direct coupling of drug transport and ATP hydrolysis (Urbatsch et al. 1994; Senior et al. 1995). Two ATP molecules bind at Walker A and B motifs that are highly conserved among ABC proteins, which are believed to provide energy for conformational changes in Pgp (Lu et al. 2005). Formation and collapse of the catalytic transition state may be directly coupled to the transport of drugs across the membrane (Ramachandra et al. 1998).

The drug-binding sites (DBS) are formed by the interaction between several transmembrane helixes (Bruggemann et al. 1992; Zhang et al. 1995; Demeule et al. 1998; Loo and Clarke 2002; Loo et al. 2006a,b). It has been demonstrated that Pgp is able to bind multiple ligands, including common cardiovascular drugs, in various sites of the drug-binding domain (Loo et al. 2003). Pgp plays an important role in clinical cardiology. Lazarowski et al. (2005) demonstrated drastic upregulation of Pgp expression in porcine cardiomyocytes exposed to hypoxic conditions as a result of chronic ischemia. In rats subjected to intermittent hypoxia secondary to obstructive sleep apnea, upregulation of Pgp expression in cardiomyocytes is a direct result of the hypoxic conditions (Dopp et al. 2009). Upregulation of Pgp as a result of a common disorder can make effective treatment of cardiac patients rather difficult.

Although the structure of human Pgp has not been determined, mouse (Aller et al. 2009), Caenorhabditis elegans (Jin et al. 2012), and several bacterial homologues exist (Dawson and Locher 2006, 2007; Ward et al. 2007). Several studies have already been conducted using the mouse Pgp model as well as using bacterial homologues. Yet, the exact mechanism of the catalytic activity and drug recruitment of Pgp remains unknown. Two hypotheses for recruitment are; (1) the ligand being drawn into the binding domain of Pgp, while it is passing through the lipid during diffusion into the cell (Raviv et al. 1990; de Graaf et al. 1996; Shapiro and Ling 1997, 1998; Loo and Clarke 1999a,b) (2) the drug being recruited from the cytoplasm and drawn up to the binding domain after complete diffusion into the cell.

One major problem with the current understanding of the mechanism of Pgp is the vast discrepancy in energy required for Pgp conformational change and the energy supplied to Pgp by the hydrolysis of two ATP molecules. In their computational study, St-Pierre et al. (2012) showed that with ADP bound and unbound, 1050.34 and 780.86 kJ/mol energy were required respectively to induce a small conformational change. The only source of energy available to Pgp is from the two ATP molecules binding at the conserved Walker motifs. Given the fact that the amount of energy contributed from hydrolyzing one molecule of ATP is only 30.5 kJ/mol, which is obviously not enough to cause the conformational change in Pgp structure; the presence of a ligand must play a larger role in the process.

Molecular dynamics (MD) simulations show that Pgp readily adopts a closed inward conformation without application of any external force. Aller et al. (2009) states that the competent drug-binding conformation is the open inward conformation. Assuming that the protein should stay in open inward conformation until the drug recruitment, these two results contradict with each other. Although the mouse homology structure has only 3.6 Å resolution, extensive equilibration can allow the side chain residues to be chemically competent, allowing the performance of realistic drug–protein interaction simulations (Wise 2012).

In this study, restraints were placed on the backbone (C – the first carbon atom that attaches to a functional group) atoms of Pgp to allow it to remain in the open inward (toward cytoplasm) conformation so that DBS could be accurately determined. Once binding occurred, the constraint on C atoms were removed to determine if the binding site was successfully located, or if the energy minima location changes. The binding sites of eight cardiovascular drugs were determined and the types of interactions between drug and protein residues were analyzed. The drugs used are: Amiodarone, Bepridil, Diltiazem, Dipyridamole, Nicardipine, Nifedipine, Propranolol, and Quinidine, all of which are known Pgp transport ligands (Abernethy and Flockhart 2000). The cytoplasmic recruitment scenario was also tested to determine its viability. Here we report over 300 ns of combined simulation data discussed in the context of the molecular mechanisms that underlie the drug-binding mechanism of Pgp.

Methods

Materials

The visual molecular dynamics program suite (VMD) (Humphrey et al. 1996) was utilized for visualization and model building. MD simulations were performed using NAMD (Phillips et al. 2005). Local Linux clusters and the Kraken supercomputer at the University of Tennessee were used for calculations.

Model creation for each drug using mouse Pgp structure

Models for Pgp were creating using the solved Mus musculus (Aller et al. 2009) structure (PDB entry 3G61). Backbone and side chain hydrogens were created using psfgen (Humphrey et al. 1996). The Pgp structure was then positioned in a physiologically competent manner into a 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) bilayer. Water (TIP3P model) was added using the solvate function in VMD. Drug PDB and parameter files were obtained from SwissParam (Zoete et al. 2011). Drug structure files (PSF) were created using psfgen. A total of eight different models were created for each drug where the final system contained mouse Pgp, POPC phospholipids, TIP3 water molecules, and one of the eight studied drugs (approximately 169190 atoms). Energy minimizations to eliminate any Van der Waals clashes between atoms and poor geometries of side chain residues and extensive equilibrations of the lipid, the protein, and the entire system were performed using NAMD to ensure the integrity of the system before any simulation experiments were performed.

MD simulations

MD simulations used NAMD 2.6 (Phillips et al. 2005) and the CHARMM27 force field (MacKerell et al. 1998) with constant temperature and pressure (NPT ensemble) in a periodic cell using Langevin temperature and pressure controls, as well as particle-mesh Ewald electrostatic calculations at 310 K. Steered MD simulations (Schlitter et al. 1994) were performed using target coordinates from mouse Pgp. Constraint files were written using the VMD package mdff (Trabuco et al. 2008). Constraining forces with a magnitude of 100.00 (kcal/mol Å2) were applied to C atoms of the mouse Pgp backbone. Drugs were placed equidistant from each TMD within the neutral zone that contains the proposed binding pocket (Bruggemann et al. 1992; Demeule et al. 1998; Loo and Clarke 2002; Longley and Johnston 2005; Loo et al. 2006a,b) of mouse Pgp (refer to Fig.1). The drug was fixed and the system was minimized, allowing just side chain residues to move freely. After 2000 steps of minimization, only C atoms of the protein were constrained, with no forces applied on any other atom in the system. The drug was allowed to move freely within the binding pocket of mouse Pgp. Root mean squared deviation (RMSD) calculations were performed after each simulation using the VMD RMSD Trajectory tool. RMSD calculations were used to determine the stability of the drug–protein interactions and binding site location. After the drug reached the binding site, the constraints on the C atoms were removed and additional 30 ns equilibration simulations were performed. Again, drug–protein interaction stability was analyzed using additional RMSD calculations. Throughout each simulation the RMSD of the protein backbone and secondary structure was tracked to ensure proper folding.

Figure 1

Structure of P-glycoprotein placed in a POPC phospholipid membrane. Pgp is shown in cartoon representation and colored gray. The phosphate heads are shown in sphere representation and colored purple. Pgp has three distinct regions: the charged extracellular matrix (ECM) region, the neutral lipid domain, and the charged cytoplasmic region. Conserved ATP-binding domains Walker A and Walker B are shown in spheres and colored teal. The drug-binding belt is represented by the residues colored orange that start at the boarder of the neutral lipid domain and extend into the charged ECM region. The drug Diltiazem colored magenta has been placed between the NBDs of Pgp demonstrating the initial drug placement for the cytoplasmic recruitment study. The lipid tails and water molecules are not shown. POPC, 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine.

The cytoplasmic recruitment model was analyzed using SMD. Drugs were placed in a neutral conformation between the two NBDs of Pgp, equidistant from the charged residues (see Fig.1). The C atoms of Pgp were constrained with a force of 100.00 kcal/mol Å2 in the open inward conformation. A continuous pulling velocity equivalent to 5 Å/ns was applied toward to protein interior in a direction normal to the plane of the membrane (Z axis). Each drug was free to move in the X and Y directions. The force applied to each drug was calculated and tracked at each time step until it reached the binding belt.

Statistical analysis

Stability of protein-drug binding was determined using RMSD calculations. Throughout each equilibration simulation the position of drug and protein atoms were tracked and their coordinates were recorded in five picosecond intervals. These coordinates were used to create a plot of atom movement in relation to their individual starting positions. Special attention was given to the portion of the plot that showed suspected binding stability as seen by a large plateau showing minimal additional atom movement. These plots were created for each atom of each residue that was suspected to play a role in drug binding, as well as the individual atoms of the drug. On average the plateau sections of each graph contained 84 data points. The standard deviations of each plateau were determined. A predetermined threshold of 0.5 Å was set to determine significant drug stabilization.

Results

Cardiovascular drug–Pgp interactions involve hydrogen bonding, aromatic cage formation and hydrophobic packing

The structure of Pgp contains distinct large areas that are charged and neutral. The portion of Pgp exposed to both the Extracellular matrix (ECM) and cytoplasm contains a large density of both positively and negatively charged residues while the structure of Pgp embedded in the lipid layer is primarily neutral (see Fig.1). The boundaries of the extracellular charged layer above the TMD and the neutral zone of the TMD do not intersect and are 2.5 Å away from each other. The neutral zone of Pgp is approximately 33 Å thick.

Each drug was placed equidistant from each TMD within the neutral region that contains the binding belt of Pgp. The drug was fixed in place and the system was minimized. The backbone (C) atoms of Pgp were constrained so that Pgp remained in the open-inward conformation, which is the competent drug-binding conformation. The stability of the C atoms was monitored during each simulation, in order to make sure that the protein stays in the open-inward orientation. Five C atoms, randomly selected from start-middle-end regions of each TMD, yielded an average RMSD of 0.259 Å. There was no restriction on the side chain atoms. Constraints were then removed from the drug prior to equilibration simulations. Residues within 5 Å of each drug were tracked throughout the equilibration simulations, and RMSD calculations were performed on the drug at each step. Once an RMSD plateau (defined as 0.5 Å or less for at least 2 ns) was reached, all atoms of each residue (excluding C) within 5 Å were recorded (see Table1). Each residue in Table1 was analyzed for the presence of possible hydrogen bond formations with the drug. The requirement for a hydrogen bond is defined as the existence of an electronegative atom within 3 Å of hydrogen atoms. Both side chain and backbone atoms were analyzed for hydrogen bond presence, with extra attention given to the side chains because they were unrestrained throughout the simulation. All of these atoms and their RMSD calculations are listed in Table1. It was noted that a significant amount of aromatic residues, both phenylalanine, and tyrosine residues were present around the drug both during and at the end of the simulation. F71, F332, F339, F728, F953, F974, F979, and Y110, Y113, Y114, Y303, Y946, and Y949 are members of the binding sites. The rings of each aromatic residue side chain aligned itself so that the region representing the pi electron cloud was facing the drug (see Fig.2). This configuration formed an aromatic cage around the drug. RMSD calculations were performed on one aromatic carbon of each ring and were plotted over time. It was observed that these carbons freely oscillated (at least 1 Å of movement) until getting into contact with the drug. A spatial rearrangement was made, and then the ring was stabilized for the rest of the simulation. This aromatic cage formed around every drug, but the aromatic residues that formed the cage changed slightly, while F71 and Y113 were part of each individual drug-binding site. Hydrophobic packing around the drug was also observed to take part in binding. Glycine, valine, leucine, and isoleucine hydrophobic side chains surround the drug in addition to the aromatic cage. RMSD calculations of hydrophobic side chain atoms revealed the same conclusion as the aromatic carbons in the rings. The hydrophobic atoms freely oscillated, then came in contact with the drug and stabilized in space due to hydrophobic packing interactions. These atoms showed at least 1 Å of movement and then plateaued so that less than 0.5 Å of movement was seen for at least 2 ns, to demonstrate stability. Each drug stabilized in a unique binding site although there was some overlap in terms of residue interaction between the drugs. Residues involved in hydrophobic packing were unique from drug to drug as well as residues involved in hydrogen bond formation. However, M34, M35, M67, M68, and S333, S340 were involved in almost all drug stabilizations, lending their sulfur and hydroxyl side chains, respectively, to provide electrostatic support and stability.

Table 1

Summary of the RMSD calculations used to determine residue stability and drug binding

Drug name	Residue-atom of interaction (RMSD [Å], time to reach stable RMSD [ns])	Drug RMSD (Å) and time to reach stability (ns)
Amiodarone	M74-SD [2.9 A, 3.0 ns], M67-SD [5.2A, 2.5 ns], M67-HE1 [6.6 A, 1.4 ns], L971-CD2 [2.8 A, 4.1 ns], V978-O [3.2A, 2.9 ns], V978-CG1 [2.7A, 2.7 ns], V978-CG2 [3.0A, 2.5 ns], I977-CD [4.4A, 2.5 ns], Y110-OH [2.3 A, 3.6 ns], Y113-OH [4.0 A, 1.2 ns], F953-CD1 [1.6 A, 2.1 ns], F71-CD1 [4.2A, 2.4 ns], F71-HB2 [2.5A, 2.2 ns], F974-CD1 [3.1A, 5.2 ns], F979-CD1 [2.4A, 5.5 ns], Y949-OH [3.1A, 2.7 ns]	9.4A, 5.4 ns
Bepridil	M68-CE [4.8A, 2.5 ns], M68-SD [2.6A, 1.9 ns], M67-SD [3.7A, 2.0 ns], M67-CE [3.9A, 2.4 ns], M945-CE [5.1A, 3.4 ns], M945-SD [2.9A, 4.2 ns], I977-CD [3.6A, 2.1 ns], I977-CG1 [2.4A, 1.4 ns], L64-CD2 [3.1A, 2.5 ns], I117-CD [3.3A, 3.1 ns], Y113-OH [4.8A, 4.0 ns], F953-CD1 [2.3A, 1.3 ns], F974-CD1 [1.9A, 2.9 ns], F332-O [2.9A, 1.4 ns], F332-CD1 [2.4A, 3.3 ns], Y114-OH [1.1A, 3.4 ns], Y946-OH [1.1A, 3.5 ns], Y949-OH [1.5A, 1.7 ns], F71-CD1 [3.5A, 2.1 ns]	5.8A, 4.1 ns
Diltiazem	M68-SD [2.3A, 1.8 ns], M68-CE [3.6A, 3.4 ns], M67-CE [6.3A, 2.5 ns], M67-SD [5.2A, 2.3 ns], S333-OG [1.3A, 3.2 ns], S340-OG [1.0A, 2.2 ns], Q343-OE1 [4.2A, 1.9 ns], Q343-NE2 [3.5A, 2.4 ns], Q343-CD [2.6A, 1.7 ns], I336-O [6.5A, 2.8 ns], L977-CG2 [2.6A, 2.2 ns], L977-CG1 [2.3A, 1.2 ns], P65-CG [1.4A, 2.4 ns], L64-CD2 [4.5A, 2.0 ns], L64-CD1 [3.7A, 1.4 ns], L64-CG [2.8A, 2.8 ns], F728-CD1 [2.4A, 3.4 ns], Y114-CD1 [3.3A, 4.4 ns] F332-CD1 [3.4A, 1.8 ns]	9.8A, 3.2 ns
Dipyridamole	S725-OG [1.8A, 1.4 ns], S333-OG [2.9A, 2.2 ns], M982-CE [3.8A, 2.9 ns], M68-SD [5.0A, 3.2 ns], M68-CE [5.3A, 3.8 ns], A981-CB [1.1A, 2.4 ns], I977-CD [2.6A, 1.4 ns], I977-CG1 [1.8A, 1.1 ns], V978-O [2.4A, 2.6 ns], I336-CD [2.6A, 1.8 ns], L64-CD1 [5.0A, 2.1 ns] L64-CG [3.2A, 1.9 ns], L64-CD2 [4.5A, 2.0 ns], F724-CD1 [1.5A, 3.8 ns], Y303-OH [10.1A, 3 ns], Y303-CD1 [4.2A, 3.2 ns], F332-CD1 [4.1A, 2.2 ns], F332-O [2.8A, 1.2 ns], F979-CD1 [1.2A, 3.2 ns], F974-CD1 [6.0A, 1.8 ns], F974-HB2 [3.1A, 1.7 ns], F728-CD1 [2.2A, 2.3 ns], F71-CD1 [3.5A, 2.9 ns]	7.6A, 5.1 ns
Nicardipine	M945-CE [5.2A, 1.5 ns], M945-SD [3.2A, 1.3 ns], M67-CE [4.7A, 1.9 ns], M68-CE [4.0A, 1.7 ns], I117-CD [1.2A, 2.5 ns], I977-CD [3.1A, 1.5 ns], I977-HG21 [3.3A, 2.1 ns], V121-O [1.1A, 3.1 ns], P65-CG [1.4A, 2.4 ns], L64-CD2 [1.4A, 1.4 ns], Y113-OH [4.5A, 1.9 ns], F71-CD1 [2.8A, 3.1 ns], F953-CD1 [2.0A, 2.6 ns], Y949-OH [3.6A, 3.1 ns], F979-CD1 [1.6A, 2.1 ns], Y946-OH [1.2A, 1.2 ns], F974-CD1 [4.5A, 2.1 ns]	9.3A, 5.9 ns
Nifedipine	M982-CE [1.8A, 2.1 ns], Q721-NE2 [2.4A, 3.1 ns], Q721-OE1 [2.0A, 1.9 ns], P722-CG [1.6A, 2.4 ns], A981-O [1.2A, 2.3 ns], I336-CD [2.2A, 1.9 ns] F979-CD1 [1.5A, 2.5 ns], F724-CE2 [4.2A, 2.5 ns], Y303-OH [2.0A, 2.8 ns], F339-CD1 [3.1A, 2.2 ns]	4.3A, 6.4 ns
Propranolol	M945-CE [5.1A, 5.2 ns], M945-SD [3.1A, 8.2 ns], S975-OG [2.1A, 2.9 ns], M68-CE [3.2A, 6.1 ns], M67-SD [5.2A, 4.9 ns], M67-CE [7.0A, 4.8 ns], I977-CG2 [2.5A, 2.7 ns], I977-CD [3.7A, 6.3 ns], V973-CG1 [2.6A, 1.9 ns], L64-CD1 [4.3A, 7.6 ns], L65-CD2 [4.7A, 8.4 ns], L65-CG [2.7A, 7.4 ns] V973-CG2 [3.0A, 2.1 ns], I117-CD [3.5A, 3.4 ns], V978-CG2-2.2A, 3.9 ns], Y949-OH [2.1A, 3.0 ns], Y946-OH [1.2A, 5.9 ns], F974-CD1 [3.1A, 3.2 ns]	11.1A, 8.1 ns
Quinidine	M74-SD [12.3A, 2.2 ns], M74-SD CG [14.1A, 2.3 ns], I70-CG2 [19.0A, 2.5 ns], V973-CG2 [11.5A, 2.1 ns], V970-CG2 [14.8A, 1.8 ns], V970-CG1 [13.7A, 2.5 ns], F71-CD1 [10.3A, 1.9 ns], F71-CB [11.7A, 2.1 ns], Y113-OH [12.0A, 1.5 ns], Y110-OH [9.2A, 1.9 ns], F974-CD1 [15.4A, 2.6 ns], F953-CD1 [7.8A, 2.8 ns]	6.8A, 2.9 ns

Amino acids identified to interact with each drug in the binding belt of P-Glycoprotein. Single letter amino acid codes are listed along with the atom name, total atom RMSD in Å, and time taken to reach stability in nanoseconds. Green colored entries indicate residues involved in hydrogen bonding. Red colored entries show residues involved in hydrophobic packing. Lastly, purple colored entries demonstrate residues involved in forming an aromatic cage around the drug. The right column of the table shows the drug total RMSD in Å and the time taken to reach stability in nanoseconds. RMSD, root mean squared deviation.

Figure 2

Snapshot of Amiodarone at the binding belt, as a characteristic example of three crucial protein – drug interactions. Residues are color coded according to the type of interaction exhibited with Amiodarone. Residues V978, I977, and L971 are colored red and represent hydrophobic packing interactions characterized by weak electrostatic forces. Residues M74 and M67 are colored green. The dashed lines represent the hydrogen bonding between the sulfur atoms of M74, M67 and hydrogens of Amiodarone as well as between hydrogen atom HE1 of M67 and an iodine atom of the drug. Aromatic cage formation is demonstrated by residues F953, Y110, F979, and Y113. These residues are colored purple and show cage formation by aligning their π electrons with Amiodarone.

When placed together, the unique binding sites collectively formed a belt of activity that is represented in Figure1. This belt is located in the nonpolar lipid region of the protein between the two regions of charged residues (refer to Fig.1), which is also consistent with findings from previous studies. Several experimental studies conducted by Loo et al. show that the TMDs are essential for drug recruitment, binding, and expulsion from Pgp. It is suggested that TM 1, 7, and 11 form part of the drug domain (Loo and Clarke 1999a,b, 2002; Loo et al. 2006a,b). Zhang et al. (1995) propose that TM 12 and the loop between TM 11 and 12 form part of the drug-binding domain as well. A computational study conducted by Wise using Sav1866 bacterial homologues, supports these experimental findings (Wise 2012). Wise observed similar aromatic residue behavior during drug-binding simulations using the Sav1866 model.

External force is required for drug recruitment from cytoplasmic side of lipid bilayer

In order to fully understand the mechanism of Pgp, it is necessary to determine the drug recruitment process. The exact mechanism of drug recruitment is not known, with the uncertainty residing in whether recruitment occurs during or after drug diffusion into the cell. In order to further understand the recruitment process, the alternative recruitment method and the cytoplasmic approach hypothesis was investigated to demonstrate its validity. Simulations were run with drugs placed inside the cell to mimic conditions similar to those found in vivo just after drug diffusion occurs.

Drugs were placed in a neutral conformation in the cytoplasmic side of the model within the space between the NBDs of Pgp (refer to Fig.1 for drug placement). The drugs were placed approximately 8 Å away from the Walker A motifs on the NBDs. During simulations with Pgp C atoms constrained, the drug did not cross the plane of lipid phosphates on the cytoplasmic side of the lipid bilayer. The drug migrated towards the lipid phosphate plane and seemed to electrostatically interact with protein charged residues at the border between the TMD and NBD. However, near the end of the simulation, the drug was released from the charged residues and fell slightly back into the cytoplasm. Following the initial simulation, an additional 20 ns simulation with Pgp C atoms relaxed was run. The drug remained in the cytoplasmic side of the lipid bilayer and was unable to readily enter the protein cavity, and reach its respective binding site. The drug interacted with the side chain atoms of the NBDs of Pgp with no appreciable specificity (results not shown). From the simulation, it appears that in order for the drug to travel up to the binding site from the cytoplasmic side of the membrane, in a timescale of tens of nanoseconds, an external force is needed.

Drugs were placed back in their original cytoplasmic position between the NBDs and constraints were applied to the C atoms of Pgp. SMD simulations were performed where a continuous pulling velocity equivalent to 5 Å/ns was applied to push the drug into the protein cavity. The drug was free to move on the plane parallel to the lipid phosphate plane. During the first 5 ns of simulation, drugs were observed to immediately interact with the Walker A motif and the conserved Q loop region. A physiologically unrealistic, 6,000 kcal/mol free-energy change was needed to free the drugs from those residues and resume movement towards the binding belt. Approximately 7 ns, and 5 Å, after removal from the Walker A and Q loop regions, drugs irreversibly attached to another mostly conserved region of the NBD. This region is made up of the following residues: G890, K891, I892, T894, E895, I897, E898, N899, and F900. When the required free-energy change to separate the drug from this site had reached 8000 kcal/mol, we stopped the Steered MD simulation, and concluded to the unlikelihood of drug recruitment from the cytoplasm. This is consistent with previous findings as experimental evidence suggests recruitment of drug molecules by Pgp directly from the lipid bilayer (Raviv et al. 1990; de Graaf et al. 1996; Shapiro and Ling 1997, 1998) and an active involvement of the TMD regions in the recruitment process (Loo and Clarke 1999a,b).

Water diffuses into the binding pocket postdrug stabilization

The Pgp model was equilibrated extensively prior to each run of simulations for 30 ns. Then during simulation of each drug constraints were applied to backbone C atoms of Pgp to allow an open-inward conformation. During equilibration simulations, water molecules were noted to be flowing in through a small opening in the TMD of Pgp just above the lipid membrane plane on the ECM side of the model (see Fig.3). The opening occurred during the drug docking equilibration simulation just after drug stabilization was reached. Residues that constitute the opening are as follows: I58, H60, G61, P65, F190, M193, A194, F197, I201, S333, I336, G337, and Q343. Water molecules initially began to surround the drug and then aggregated at the bottom of the drug (see Fig.3). The water molecules had inserted themselves between the drug and protein residues, which terminated drug-protein stabilization. This behavior was noted for all of the drugs. After displacement from stable binding, each drug was simulated for an additional 20 ns with its trajectory tracked using RMSD calculations. On average, each drug moved an additional 4.1 Å in the direction of the ECM after the water molecules displaced it from its stable binding site. It was noted that no drug was able to reach stable binding within the simulation timeframe following its displacement from its original binding site. To confirm the validity of this observation, drug-free control simulations were performed with identical conditions as in the binding equilibration study. During the control drug-free simulations, no water molecules were noted to flow into the binding region showing that the presence of the drug is required for water flow. In either simulation condition, no appreciable amount of water was seen to flow up from the cytoplasmic side of the protein.

Figure 3

Timescale representation of water influx into the binding belt of Pgp, after the drug reached stabilization. (A) 4.5 ns, (B) 5.5 ns, (C) 6.5 ns, and (D) 7.5 ns. The drug Amiodarone is shown in sticks representation and colored teal. Pgp is shown in cartoon representation and colored gray.

Discussion

The drug-binding domain of Pgp is not well characterized. Specific TMD helices have been identified experimentally to play a role in drug recruitment and binding. To further characterize the drug-binding domain of Pgp equilibration simulations were run using eight cardiovascular drugs that are known to bind and be transported by Pgp. Each of these ligands has stabilized in their corresponding binding sites within the first 10 ns of the equilibration. Additional 30 ns long simulations with C atoms relaxed were performed to confirm the stability of drug binding. Although these simulation times may appear truncated, it was noted that, after the drug free energy converged to local minima within the first 10 ns, additional simulation time, up to 50 ns, did not produce new information. Therefore no additional equilibration was performed, however we accept the possibility that additional changes in binding orientations might very well occur on a much longer timescale. After each simulation, residues within 5 Å of the drug were tracked and analyzed for potential roles in binding.

It is found that, the binding site of Pgp is not defined by a small number of residues, but rather can be deemed a belt of binding activity (see Fig.1). Each drug studied had a distinct location that was most energetically favorable. While there were unique residues to each drug, there was a significant overlap in binding locations between the drugs, supporting the idea of the presence of the binding belt. The interactions between the drug and the TMD of Pgp can be classified into three main categories: hydrogen bonding, hydrophobic packing, and aromatic interaction. Throughout the movement of each drug, various numbers of hydrogen bonds with the protein were observed. Each drug had at least two hydrogen bonds present and many other weaker electrostatic interactions upon binding. The sulfur atom from the side chain of M67 and M68 were common among the drugs as well as the OH group from several different Serine residues. Extensive hydrophobic packing was also observed during binding. The hydrophobic residues involved were not conserved across the different drugs, yet V978 was present at almost every binding site.

The most striking observation made during the drug-binding simulations was the presence of distinct aromatic cages that formed around the drug just after binding. Both Phenylalanine and Tyrosine residues were involved in the cage formation. Y113 and F71 were highly conserved across the different binding sites with Y113 present five and F71 six occasions out of eight drugs simulated. In all cases, the aromatic rings of the residues aligned their site representing the π electron cloud to face the drug to provide electrostatic stability for binding. It should be noted that MD simulations do not allow pi electron clouds to be explicitly present. The rings were present on all sides of the drug and ranged in number from 3 as seen in Diltiazem and Propranolol, to as many as 8 and 10 as seen in Bepridil and Nifedipine, respectively. It was noted that based off of the RMSD calculations of the aromatic carbons, the atoms were free to oscillate until they were stabilized in place when in contact with the drug. The RMSD of the CD1 atoms of each aromatic ring before drug binding was between 1 Å and 2 Å of harmonic movement. When the aromatic rings came in contact with the drug, the RMSD showed a large single rapid movement and then a plateau with less than 0.5 Å of oscillation, showing that the aromatic ring was stabilized in place forming part of the cage around the drug. This observation was seen for every drug studied, demonstrating its importance in drug binding and stabilization. These results can be applied to new developments to help combat acquired drug resistance by creating an inhibitor specific to these types of drug–protein interactions.

Understanding the recruitment method of Pgp is pivotal to understanding its overall mechanism and developing methods to combat its drug resistance properties. The exact mechanism remains unknown, however there is much experimental evidence supporting drug recruitment from the lipid bilayer as opposed to from the cytoplasmic side of the membrane. When placed in the cytoplasmic side of the membrane, the drug was unable to diffuse into the phosphate layer to reach the binding belt of Pgp. When placed within the cytoplasmic structure of Pgp, the drugs did not migrate up towards the binding belt, but instead interacted with the side chain atoms of the NBD regions of Pgp. Most drugs are weak organic acids or bases, existing in unionized and ionized forms in an aqueous environment. The unionized form is usually lipid soluble (lipophilic) and diffuses readily across cell membranes. In physiological conditions diffusion, driven by a concentration gradient, is often enough for drugs to cross the lipid layer. Since Pgp binds and transports ligands effectively, the concentration of the drug would be higher in the ECM than in the cytoplasm. In the case of the cytoplasmic recruitment mechanism, the drug would not only have to flow against its concentration gradient, but overcome any electrostatic interactions that would naturally occur with the NBDs of Pgp. During pulling, drugs were free to move in the either the X or Y directions, and at least once during the simulation was recruited by a section of the NBDs. The Walker A motif, Q loop, and another mostly conserve loop of residues were the sites of recruitment and caused in a physiological sense, irreversible binding of the drug to the NBD due to the high amount of force required to move the drug away from the NBDs and towards the binding belt. During the SMD simulations, an excessive free energy change was required to free each drug from its binding position on the Walker A motif and Q loop. After removal from the Walker A motif and Q loop, the drugs were immediately recruited by another mostly conserved loop of residues and again another astronomical free energy change was required to free the drugs. Since two ATP molecules are the only known source of energy in the system, and they are present after drug binding, this scenario shows little promise.

Besides potential drugs to be transported and ATP, the only molecules in contact with Pgp are water and lipids. The availability of water molecules may play a role in drug expulsion and conformational changes of Pgp. During each constrained simulation with drugs present, water molecules were observed to be flowing into the binding region of Pgp from the ECM through a pore in the TMD which only occurs when the drugs dock to the binding belt. In the absence of the drug, no water was observed in the binding region coming neither from the cytoplasm nor the ECM. The presence of the drug clearly plays a role with water influx into the binding pocket of Pgp. The exact consequences of the existence of water are yet to be determined. Although this report does not investigate the full drug expulsion mechanism by Pgp, one might speculate that water influx may play a role in helping Pgp to adopt an open outward conformation and dispel the drug back into the ECM. Since the amount of brute force required to open Pgp is much higher than what two ATP can provide, another process must supply the energy required for the conformational change. Water may function to provide a positive change in internal pressure of Pgp's binding pocket. This would cause the hydrostatic pressure outside Pgp's binding pocket to be close to the internal hydrostatic pressure, allowing less energy to be required to make the conformational change and expel the drug. This finding should encourage future studies to determine the role of water in Pgp mechanism.

We acknowledge that very recently there has been an updated mouse Pgp structure published (Li et al. 2013). While there have been some improvements to the structure, we believe that these changes have not affected our results since the updates did not occur within our proposed binding belt. Our model compares well to the revised structure with an average backbone C RMSD of 1.45 Å and average side chain RMSD of 1.63 Å. Our adequate equilibration procedure has placed the side chain residues in chemically competent locations despite the experimental uncertainty suggested by the resolution of the solved structure. Using MolProbity (Davis et al. 2007; Chen et al. 2010), a postequilibration analysis was performed. A single MolProbity score was reported following analysis. The MolProbity analysis combines individual parameters including the clashscore, rotamer, and Ramachandran evaluations into a single score, which is normalized to be on the same scale as X-ray resolution. After 5 ns of equilibration, a MolProbity score of 1.69 Å was observed. After 25 ns of equilibration, a MolProbity score of 1.53 Å was observed. A resolution of 1.53 Å is sufficient for determining activity of side chain residues and allows for accurate analysis of drug–protein interactions. However, all models have small inherent flaws whether they are solved crystal structures or self-built homology models. At this current time, the mouse homology model is the best available for simulation study. Our model is chemically competent and behaves similar to its in vivo counterpart (Li et al. 2013).