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 Å²) 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.

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 Å² 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.

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 Table 1). Each residue in Table 1 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 Table 1. 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.

When placed together, the unique binding sites collectively formed a belt of activity that is represented in Figure 1. 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.