Theoretical studies on 1,4-dihydropyridine derivatives as P-glycoprotein allosteric inhibitors: insights on symmetry and stereochemistry

Abstract P-glycoprotein (P-gp) is a key efflux pump involved in cellular multidrug resistance (MDR), lowering the concentration of many anticancer drugs in tumor cells by pumping them into the extracellular milieu. While previous studies identified 1,4-dihydropyridines (DHP) as putative P-gp allosteric inhibitors, none reported the effect of stereochemistry on the ability of DHPs to bind P-gp. In the present study both symmetric (1) and asymmetric (2 R,S) DHPs were designed as P-gp inhibitors and, after biological evaluation, molecular docking and molecular dynamics simulation (MD) studies were performed to gain insights on how both scaffolds interact with P-gp. The results were thoroughly analyzed i) to evaluate the role of the different substituents and ii) to assess how stereochemistry may affect binding of DHPs to P-gp. Our results suggest that both the substitution pattern and stereochemistry may have a significant impact not only in drug binding but also on membrane permeation/accumulation, thus compromising in which site the DHPs may exert their effect as P-gp efflux inhibitors. Therefore, it is our conclusion that the stereochemistry cannot be neglected during the development of novel 1,4-dihydropyridine derivatives. Communicated by Ramaswamy H. Sarma


Introduction
P-glycoprotein (P-gp, ABCB1) is an ATP Binding Cassette (ABC) transporter that uses the energy of ATP hydrolysis to efflux structurally unrelated molecules and drugs into the extracellular milieu. P-gp is commonly related to multidrug resistance (MDR) phenomenon because this protein is often over-expressed in tumor cells, lowering the intracellular concentration of anticancer agents and thus reducing the efficacy of anticancer drugs. Hence, the expression and activity of P-gp in tumor cells became a serious concern in current chemotherapy regimens (Eckford & Sharom, 2009;Sharom, 2008). P-gp is composed by two transmembrane domains (TMD1 and TMD2) and two nucleotide-binding domains (NBD1 and NBD2), with each TMD containing six transmembrane (TM) a-helices, enclosing a large drug-binding pocket (DBP) deeply buried within the membrane (Aller et al., 2009;Loo & Clarke, 2001) (Figure 1). Binding of substrates (or inhibitors) to the different drug-binding sites (DBSs) within the internal DBP is thought to be the driving force for conformational changes, which are then propagated through intracellular coupling helices (ICHs) into the NBDs. Following, ATP binding and consequent hydrolysis shifts P-gp into an outward-facing conformation, and ADP release resets the transporter into its initial state (Futamata et al., 2020;Gyimesi et al., 2011;Siarheyeva et al., 2010;van Wonderen et al., 2014). Only recently, both inward-(PDB ID: 6QEX) and outward-facing (PDB ID: 6C0V) conformations for human Pgp became available (Alam et al., 2019;Kim & Chen, 2018). Still, the releases of the early murine crystallographic structures (PDB ID: 3G5U; 4M1M) enabled studies focused on the mechanism by which substrates are effluxed from the cell into the extracellular matrix (Pan & Aller, 2015;Prajapati & Sangamwar, 2014;Wise, 2012). However, the efflux activity of membrane proteins is known to also be dependent on the composition of the surrounding lipid bilayer, namely on the lipid type and cholesterol content (Ferreira et al., 2015). More recently, studies in both mouse and human P-gp models were published in an attempt to shed additional insights about the mechanism by which P-gp is able to transport drugs (Alam et al., 2019;Furuta & Sakurai, 2018;McCormick et al., 2015;Sz€ oll} osi et al., 2020).
As P-gp inhibitors have been found to be effective in overcoming MDR, the design and synthesis of efflux pump inhibitors remains an important goal (Ferreira et al., 2015b). An increasing number of molecules, divided into three major generations, have been found to possess inhibitory properties against P-gp (Mollazadeh, Sahebkar, et al., 2018;Palmeira et al., 2012). Herein, verapamil was one of the firstly discovered and is still considered a reliable positive control in biological evaluations (however, due to the cardiovascular side effects its use as an inhibitor is limited) (Tsuruo et al., 1981). Currently, tariquidar analogs (derived from third-generation P-gp modulators) are the only ones still being investigated in clinical trials (Weidner et al., 2016). Another alternative, however, relies on natural products to modulate P-gp (Palmeira et al., 2012), but because of the dual inhibitory effect against Cytochrome P450 and P-gp a wide range of the identified P-gp inhibitors could not reach clinical trials (Ondieki et al., 2017).
In silico methodologies are increasingly becoming the first step in drug discovery, and Molecular Dynamics (MD) simulation techniques are a particular form of these methods (Lill, 2015). Herein, MD studies have been employed to gain important mechanistic insights from such targets and to guide the development of novel efflux modulators (Baptista et al., 2016;Ferreira et al., 2012Ferreira et al., , 2018Isca et al., 2020). More particularly, and due to the high P-gp flexibility, such methods are also able to provide additional insights about conformational changes and other structural alterations related to ligand binding (Ferreira et al., 2015a;Subramanian, 2015).
Recently, it has been proposed the existence of allosteric sites, in close vicinity to the ICHs, that upon ligand binding may impair signal transmission between the transmembraneand nucleotide-binding (TMD/NBD) domains responsible for driving conformational changes leading to efflux (Ferreira et al., 2017). Initially identified from photoaffinity labeling studies, a putative binding site was proposed to exist within TMD4/TMD5, ranging aminoacids 468-527 (Borchers et al., 2002;Bruggemann et al., 1989;Ferry et al., 1992;Isenberg et al., 2001). However, the first experimental evidence was only obtained in 2010 using a combined approach involving auxin transport, binding studies and molecular docking (Kim et al., 2010), identifying E502 (NBD1) and F792 (NBD2) as key residues in modulating the efflux activity of an A. thaliana Pgp homologue (corresponding to E526 and F804 in human P-gp). As these studies were recently corroborated by molecular docking and MD studies (Ferreira et al., 2015a;Shahraki et al., 2018;Subramanian et al., 2015), both approaches can be considered as valuable tools for the development of a new generation of P-gp efflux inhibitors.
Two recent papers, by Shahraki et al. (2018) and Ranjbar et al. (2017), used most of the above information to gain insights on the binding modes of 1,4-dihydropyridines (DHP) and tetrahydroquinolinones at the NBDs of human P-gp, upholding the importance of several residues for DHP binding (V472, S474, Q475, R905, S909, T911 and Q912). Herein, and based on our first-generation DHPs, new derivatives containing a thiophenyl substitution at positions C-3 and/or C-5 were designed and studied as MDR reversers (Mollazadeh et al., 2017(Mollazadeh et al., , 2019Mollazadeh, Moosavi, et al., 2018) (Scheme 1). However, and unlike all previous studies, we are interested to gain further insights on the potential impact that stereochemistry and symmetry may have on MDR-reversal activity reported for the dihydropyridine scaffold.
Both symmetric (1) and asymmetric compounds (2R,S) were evaluated as P-gp inhibitors and compared with verapamil through in silico and in vitro approaches. The first test molecule, a previously synthesized achiral compound (1), was designed with a 4-chlorothiophenyl group at both positions C-3 and C-5. The other two test compounds, the (S)and (R)-enantiomers of 2, have a distinct phenylcarbamoyl group replacing one of the thiophenyl moieties and were used to assess the effect of stereochemistry on DHP binding. Herein, and taken into account that (i) F804 is located at the C-terminal NBD (Ferreira et al., 2017) and (ii) flavonoids were also described to bind at a vicinal ATP-binding site located at this particular NBD (Conseil et al., 1998), we additionally evaluated the interaction of all test compounds at the NBD2.

Synthesis of compounds
The synthesis of compounds 1 and 2 was already published elsewhere (Mollazadeh, Moosavi, et al., 2018), but nonetheless and due to their capability to quickly oxidize, an alternative approach was used to synthesize and protect them from oxidation. In this way, we used a heating bath at 72 C in which a sealed tube containing all reactants (3-nitrobenzaldehyde, phenyl thioacetone and ammonium carbonate) was put under constant stirring. Time of reaction and amounts of reactants remained similar to that previously published (NMR data available in the Supplementary material file, Figures S1-S2).

Citotoxicity and doxorubicin combination assays
The cell viability for normal, parental and resistant cell lines was determined using the 3-(4,5-dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide (MTT) assay, in which the MTT reagent is reduced to formazan via mitochondrial dehydrogenases only in viable cells. Cell viability of compounds 1 and 2 was evaluated in human carcinoma gastric (EPG85-257P) cells, its resistant counterpart (EPG85-257RDB) and also in normal HUVE cells. All cell lines were harvested and seeded at a density of 5 Â 10 3 cells/well in flat-bottom 96well plates and incubated for 24 h at 37 C. Following, compounds 1, 2 and verapamil (VER) were added to each well at a final concentration of 15 lM.
For the combination studies, doxorubicin was additionally added to a final concentration of 0.7 lM. After 72 h, the medium was removed and cells were washed with PBS. MTT reagent (5 mg/ ml) was added to each well and, after 4 h incubation, cell viability was measured by a spectrophotometric method. Excitation and emission wavelengths were 550 and 630 nm, respectively.
2.3. Flow cytometric rhodamine-123 (R123) efflux assay EPG85-257RDB cells were harvested, seeded at a density of 2 Â 10 5 cells/well in 12-well micro plate and incubated for 24 h at 37 C, after which R123 (10 lM) was added to each well and the cells were incubated for 45 min at 37 C. Following, verapamil (positive control, 15 lM), compounds 1 (15 lM) and 2 (as a racemic mixture, 15 lM) were added to special wells and the cells were further incubated for 30 min at 37 C. Cells were washed two times with PBS and harvested by trypsinization.
The cells were pelleted at 2000Âg for 5 min at 4 C, supernatants were removed and the cells were re-suspended in ice-cold PBS. Fluorescence intensity of R123 accumulated in the cells was measured by a BD FACS Calibur Flow Cytometry System (BD Biosciences, San Jose, USA). Excitation and emission wavelengths were 488 and 530 nm, respectively. Mean fluorescence intensities (MFI) were then calculated as

MFI ¼
FLOW of test group FLOW of control group

Molecular docking
A human P-gp model, previously refined and equilibrated in POPC lipid bilayer during 500-ns and ready for in silico studies in GROMACS was downloaded from http://chemistrybits.com/ downloads/systems/ (Bonito et al., 2020). This structure was used to explore the direct interaction between protein and ligand comprising the whole DBP. As dihydropyridines are known to bind at a specific location close to NBD1 (Isenberg et al., 2001) and flavonoids are found to bind to NBD2 (Conseil et al., 1998), all compounds but verapamil were additionally evaluated in both locations. Chemical structures of these compounds were drawn and minimized in Discovery Studio Accelrys software, exported into AutodockTools-1.5.6 (Morris et al., 2009) and saved as pdbqt files for docking. Autodock VINA 1.1.2 (Trott & Olson, 2010) was used for molecular docking protocol, following previous studies (Ferreira et al., 2013;Shahraki et al., 2018). A docking box, comprising the whole internal cavity identified by Aller et al. (2009) and Li et al. (2014), was defined with dimensions xyz of 33.75 Â 26.25 Â 40.50 Å 3 , and centered at the DBP (xy corresponds to the membrane plane). Two other docking boxes with dimensions xyz of 28.50 Â 28.50 Â 28.50 Å 3 were defined and centered at NBD1 (42.42 Â 56.44 Â 44.94 Å) or at Scheme 1. Structure of symmetric compound 1 and of both enantiomers for asymmetric compound 2. NBD2 (78.82 Â 63.54 Â 44.94 Å) ( Figure S3, Supplementary material). Herein, E526 and F804 (corresponding to E502 and F792 in A. thaliana) (Kim et al., 2010) were considered as structural references while defining both boxes. Due to the large search volume in all three locations, VINA's 'exhaustiveness' parameter was manually set to 50 and twenty conformations were generated.

Molecular dynamics simulations
The top-ranked pose for all compounds was saved in pdb format for further usage in molecular dynamics simulations. Ligand parameterization was done by using PRODRG server (Sch€ uttelkopf & van Aalten, 2004) while Antechamber program was used to calculate AM1-BCC partial charges, as previously suggested (Lemkul et al., 2010). Both human P-gp model and lipid topology were parameterized according to the GROMOS96 force field with the 54a7 parameter set, with the latter based on previous studies by  and . Energy minimizations were carried out using the steepest-descent algorithm. Temperature equilibration employed the NVT ensemble with leap-frog integrator and V-rescale (Bussi et al., 2007) temperature coupling to 303 K (10 ps). Four independent 40 ns (first MD run) or 50 ns (three MD replicates) were performed for each test compound, using an NpT ensemble with Nos e-Hoover as thermostat (303 K) and Parrinello-Rahman as barostat (1 bar) (Hoover, 1985;Nos e, 2002;Nos e & Klein, 1983). Only one 50 ns MD run was performed for verapamil (for comparison purposes). Simulations were performed with full periodic boundary conditions (PBC). LINCS (Hess et al., 1997;Hess, 2008) or SETTLE (for water molecules) (Miyamoto & Kollman, 1992) algorithms were utilized to constrain the lengths of all bonds. Short-range electrostatic and short-range van der Waals interactions used a cut-off of 1.2 Å. Long-range electrostatic interactions were computed using the particle mesh Ewald (PME) method (Darden et al., 1993;Essmann et al., 1995). Verlet (P all & Hess, 2013) cut-of schemes were applied to allow calculation of non-bonded interactions on GPUs. The results report the analysis of all four replicates for each MD run (after discarding the first 20 ns as equilibration for each run).

Data analysis
In molecular docking, the number of poses for each drugbinding site and their corresponding binding energies were obtained with Autodock Vina default functions. Residue interactions and cross interaction capabilities were obtained from LigPlot (Wallace et al., 1995) and in-house python scripts. The calculation of relative energies of binding through MD simulations used a Molecular Mechanics-Poisson Boltzman surface area (MM-PBSA) approach to estimate interaction free energies, a technique that has been increasingly used in the study of biomolecular interactions (Homeyer & Gohlke, 2012). However, and as P-gp is a membrane protein, a correction to the polar solvation energies must be taken to account the presence of the low-dielectric environment. Herein, an implicit membrane approach was used to generate ion-accessibility and dielectric maps that specifically encompass the presence of the lipid bilayer (Ferreira et al., 2015a). Hydrogen-and hydrophobic bonding and residue contact frequencies were estimated with gmx hbond (Luzar, 2000;van der Spoel et al., 2006) and g_contacts (Blau & Grubmuller, 2013), respectively. To allow a better comparison between molecules, ligand efficiency (LE) was calculated as the ratio between the number of protein-ligand contacts with frequencies (from g_contacts) greater than 0.50 and the total number of registered protein-ligand contacts. The cross-interaction capability (CIc) of the molecules was estimated as a ratio between the number of hydrogen-bond and non-bonded interactions observed between each molecule and P-gp's N-terminal and C-terminal halves, exclusively at the M-site (Ferreira et al., 2015a).

Synthesis of compounds 1 and 2
All compounds were synthesized as previously published (Mollazadeh, Moosavi, et al., 2018), with yields between 43% and 62%, and the correspondent NMR data is available in the Supplementary material file (Figures S1-S2).

Biological evaluation of compounds 1 and 2
We have reported in our previous papers (Mollazadeh et al., 2019) that DHP derivatives were cytotoxic to cancer cell lines on both EPG85-257P and its drug-resistant EPG85-257RDB cell lines, at 15 mM, via MTT assay and after 72 h incubation. Oppositely, human umbilical vein endothelial (HUVE) cells showed negligible toxic effects (Figure 2).
Herein, compound 2 is more cytotoxic on EPG85-257RDB cells when compared with 1, being this effect potentiated when each compound is combined with the anticancer drug doxorubicin. Regarding the flow cytometry efflux assay, the fluorescent dye R123 is a well-known P-gp substrate used to assess the P-gp inhibitory potential of tested compounds.
Our results showed that all tested compounds increased R123 accumulation inside cells when compared to control. Interestingly, while 1 (symmetric) was qualitatively similar to verapamil, the asymmetric DHP 2 substantially increased R123 accumulation in EPG85-257RDB cells, when compared with verapamil. The MFI (mean fluorescent intensity) of cells incubated in the presence of with compounds 1 and 2 at 15 mM were 1.7-fold and 6.1-fold higher than control, respectively, while the MFI for verapamil-treated cells was only 1.9fold higher (Figure 3).

Molecular docking studies
Molecular docking was used to obtain an initial assessment on which locations the test compounds are found within P-gp. We individually tested each compound (1, 2R and 2S) as a way to unambiguously characterize their interaction with P-gp and to evaluate any possible influence of stereochemistry and/or substitution pattern on protein-ligand binding. Our results showed that both the top-ranked docking pose ( Figure S4, Supplementary material) and most of docked conformations ( Figure S5) were found at a previously characterized M-site (Ferreira et al., 2013) and with more favorable binding energies than verapamil (Table 1). Interestingly, only symmetric compound 1 and the (R)-enantiomer of compound 2 were able to show a moderate cross-interaction capability (CIc), which seems to be in agreement with a higher potency described for other efflux inhibitors as dexverapamil or dexniguldipine (Ferreira et al., 2015b). Nonetheless, as the (S)-enantiomer has a higher number of poses in each of the three considered sites (H, R and M) with lower binding energies, we cannot fully exclude that it may also exert its activity simultaneously as a non-(when at Msite) or competitive inhibitor (at substrate-binding H-and Rsites). Regarding the interaction of DHPs and verapamil with Pgp, no significant difference was found.
When comparing the top-ranked binding energies to those obtained at both NBDs, our results show that for all cases DHPs appear to have slightly lower binding energies by $0.5 kcal.mol À1 , but well within the margin of error reported for VINA (2.85 kcal.mol À1 ) (Trott & Olson, 2010). Therefore, we proceeded to more thorough evaluations through molecular dynamics (MD) simulations to better clarify these findings.

Molecular mechanics/Poisson-Boltzmann binding energies
Protein-ligand binding can be more accurately calculated by means of a molecular mechanics-Poisson-Boltzman surface area (MM-PBSA) approach. Although neglecting the entropic term, the herein obtained relative free energies of binding can still be used to compare how compounds bind to P-gp. However, and as previously referred, P-gp is a membrane protein that requires a correction to the polar solvation energies (Ferreira et al., 2015a;Kumari et al., 2014). From the results in Table 2 it is possible to verify that the results only partially confirmed the molecular docking data. While verapamil displayed the lowest DG bind at the M-site (À69.7 ± 5.6 kcal mol À1 ) and a strong CIc (90), all other compounds did not have their highest binding affinity at the nucleotide-binding domains as suggested by the docking experiments but instead at the modulator M-site. Accordingly, the reported binding affinities were similar, ranging from À55.2 kcal.mol À1 up to À49.6 kcal.mol À1 for compounds 2R and 1, respectively (average DG bind of À52.0 ± 2.88 kcal.mol À1 ). This was accompanied by a significant increase in the CIc for the (S)-enantiomer of compound 2, achieving an almost equivalent CIc to that calculated for verapamil (80 ± 9 against 90, respectively). This seems to be in more in agreement with the lower IC 50 's reported for several (S)-enantiomers as QZ59 (Aller et al., 2009), methadone (St€ ormer et al., 2001) or mefloquine (Pham et al., 2000;Riffkin et al., 1996), when compared with their counterparts. Compounds 1 and 2R maintained or slightly reinforced a moderate CIc, in agreement to that obtained from the docking experiments.
Despite the results above, it is worth noticing that the free energies of binding for each compound at the NBDs were still calculated as favorable (Table 2). Herein, while the (R)-enantiomer reported the lowest binding affinities (i.e. highest free-energies of binding) in both NBDs, the (S)enantiomer of compound 2 revealed slightly better binding energies at NBD2 (À31.4 ± 5.3 kcal mol À1 against À26.3 ± 5.4 kcal mol À1 in NBD1). Quite interestingly, compound 1 was the only one with similar relative DG bind at both NBDs (average DG bind À35.2 ± 1.98 kcal mol À1 ).
Visual evaluation of the MD trajectories revealed that reorientation of the thiophenyl and phenylcarbamoyl moieties could explain most of the differences observed between molecular docking and MD studies.
Concerning 1, while at the M-site an initial intramolecular p-p interaction was replaced by an CH-p with little influence on the protein-ligand interactions, at both NBD1 and NBD2 a shift from its initial position was observed, bringing the molecule closer to ICH1/ICH4 and establishing HBs with D164/R905 (NBD1) or R262/T811 (NBD2).
For the (R)-and (S)-enantiomer, and when at the M-site, while for the former an intramolecular NH-p interaction between the amide group of the phenylcarbamoyl moiety and the thiophenyl group was observed (allowing additional HBs between the carbonyl group and Q725/S979), in the latter the aromatic rings are kept parallel to each other, slightly tilted ($45 ) regarding the dihydropyridine z axis and promoting p-p interactions with the protein. As in 1, both 2R and 2S shift from their initial docked position at the NBDs to new positions closer to the intracellular helices, but while the phenylcarbamoyl moiety of (R)-enantiomer was found to establish HBs with E902/R905 in NBD1, at NBD2 both aromatic moieities were found closer to R262 (HBs with the carbonyl group) and F804 (p-p interactions). Regardind 2S, interactions with ICH2 were also observed, with the thiophenyl (NBD1) or phenylcarbamoyl (NBD2) groups reorientated towards the 'linker' moiety in order to avoid the aqueous surrounding environment. Therefore, the above results suggest that while the (S)enantiomer of compound 2 may exert its inhibitory activity directly at the M-site (as a non-competitive inhibitor), the higher affinity of compound 1 towards NBDs points out a possible allosteric mechanism for this particular compound. The (R)-enantiomer of compound 2, however, as it does not seem to bind with enough affinity to both NBDs or to have strong CIc when at the M-site, may only act as a competitive inhibitor at the substrate-binding sites.

Hydrogen bonding and ligand efficiency index
The ability to establish a stable hydrogen-bond network strengthens the interaction of any given ligand to its target (Chufan et al., 2016;Desai et al., 2012). Therefore, and while in symmetric compound 1 its functional groups are mainly hydrogen-bond acceptors, is crucial to evaluate the importance of the additional amide group, able to act as hydrogen-bond donor. Simultaneously the definition of a ligand efficiency (LE) index, herein calculated as the ratio between the number of protein-ligand contacts with frequencies above 0.5 and total number of protein-ligand contacts, is also an important metric to evaluate how the compounds interact with P-gp.
From Figure 4 some interesting results came forth. Regarding the average number of HBs, hNi, while the (R)enantiomer of 2 had the highest number at NBD1, its (S)enantiomer registered the highest number at both M-site and NBD2. Interestingly, and despite the lower hNi registered for 1, the calculated HB lifetime (s) was high in both M-site (3793 ps) and NBD1 (904 ps), suggesting that 1 is capable of inducing the formation of a stable hydrogen bond networks. A closer inspection of the MD trajectory allowed us to identify that such feature was due to the interaction of the 3nitrophenyl group with I340 (M-site) or D164/A640/R649 (NBD1). Oppositely, in both compounds 2R and 2S only Y310 seem to be involved in hydrogen bonds (at the M-site), although with distinct moieties (NO 2 group in 2R and the N-H of the dihydropyridine ring in 2S). Finally, the increased ability of 2S to establish a strong HB network at the NBD2 is due to interactions with K645 and T811. Herein, both D164 (ICH1) and T811 (ICH3) are located at the intracellular coupling helices (ICHs), though to participate in signal transmission between TMD and NBD during the efflux cycle.
Regarding the compounds' ligand efficiency index, higher LEs were registered at the M-site for all compounds when compared with verapamil (0.17). At the nucleotide-binding domains, i) both (R)-and (S)-enantiomers of compound 2 had higher LEs when compared with the symmetric compound 1 and ii) LEs were in general higher at the NBD1, except for the (S)-enantiomer of compound 2 (similar to compound 2R at NBD1, but higher at NBD2).
We additionally aimed to characterize the residues that contribute the most for these values. At the M-site, LE is almost exclusively due to the interaction with aromatic residues as F72, Y310, F336, F343, F732 and F983, with S979 being the only polar residue with significant contributions for the referred index. At NBD1, the only common residue to all three compounds was R905, experimentally implied in protein maturation (Loo et al., 2008), and while compound 1 had higher interaction frequencies with residues belonging to the NBD (G430, S434 and Q438), both (R)-and (S)-enantiomers interacted with D164 and F904, identified to contribute for both the assembly and activity of human P-glycoprotein (Kwan & Gros, 1998;Loo & Clarke, 2015). Regarding NBD2, again the only maintained residue across all compounds was R262, also implied in signal transmission by mutational studies (Loo et al., 2013). Other involved residues were L258/A259 (ICH2), F804/ D805 and F811 (ICH3) and Q1118/E1119 (NBD). As experimental studies that locked ICH2/ICH3 together also inactivated human P-gp (Loo & Clarke, 2014), the simultaneous interaction with residues in both ICH2 and ICH3 by dihydropyridines is interesting because it may reveal a similar mechanism underlying allosteric inhibition when bound at the NBDs.
Therefore, it appears that either at the drug-binding Msite or at both NBDs, the tested DHPs are able to interact with key residues experimentally involved with drug efflux. However, the data also suggests that the P-gp inhibitory activity is also intimately related with the substitution pattern and the stereochemistry of the dihydropyridine scaffold. While the (S)-enantiomer is expected to have a major impact on impairing confirmation changes through binding at the M-site (due to its strong cross-interaction capability), the symmetric compound 1 is expected to bind the TMD/NBD more efficiently, thus impairing the signal transmission created by substrate binding at the drug-binding pocket.

Passive membrane permeation
We additionally aimed to evaluate how the compounds interact with the membrane, namely if (i) they are able to fully permeate into the cytoplasm, or (ii) if they accumulate inside the lipid bilayer. Numerous studies pointed out the surrounding lipid bilayer as an important counterpart in drug efflux (Ferreira et al., 2015). In order to bind at the intracellular ICHs (TMD/NBD interfaces), DHPs must be able to efficiently cross the lipid bilayer and accumulate inside the cell. Oppositely, increased retention at the plasma membrane may favor diffusion into the P-gp's drug-binding pocket (Ferreira et al. 2015c), thus enabling them to act as competitive or non-competitive efflux inhibitors (Ferreira et al. 2015b). Herein, we predicted passive membrane permeation through a dioleoyl-phosphatidilcholine (DOPC) bilayer using the PerMM online server , based on the inhomogeneous solubilitydiffusion model. Through such approach, not only the potential of mean force (PMF) profile across the membrane can be estimated ( Figure 5) but also their permeability coefficient across several membranes (Table 3).
From Figure 5 is clear that the symmetric compound 1 has a much lower energetic barrier (-2.26 kcal mol À1 ), enabling faster 'flip-flop' motions between leaflets and, theoretically, achieving larger concentrations at the cytoplasmic side. Oppositely, the energetic penalty for crossing the hydrophobic core of the membrane is larger and similar for both (R)-and (S)-enantiomers of compound 2 (À5.75 and À4.69 kcal.mol À1 , respectively), meaning that such compounds are expected to build up within the membrane and thus increasing the probability of 'sliding' into the P-gp's internal cavity where the drug-binding sites are located (Ferreira et al., 2015c;Subramanian et al., 2015). Quite interestingly, the small difference between both enantiomers regarding its the energy minima within the membrane ($0.4 nm, 1.2 < z < 1.6 nm) also agrees with the localization of the putative 'entrance' gate between TM helices 10/12 (Ferreira et al., 2012;Ferreira et al., 2015c), which may favor (S)-enantiomer access to the drug-binding sites in detriment of the (R)-enantiomer.
The logarithms of permeability coefficients estimated by the PerMM server also supports the above assumptions (Table 3). While compound 1 has, in general, higher permeability coefficients that favors cytoplasm accumulation (and thus an increased probability of interacting allosterically with P-gp at the TMD/NBD interface), the lower permeabilities calculated for both (R)-and (S)-enantiomer of compound 2 suggest that they accumulate within the lipid bilayer, thus becoming more accessible for P-gp to extract them from the hydrophobic core of the membrane. Following, when already inside the polyspecific drug-binding pocket the more favorable binding energies and calculated CIc at the M-site enables the (S)-enantiomer of compound 2 as a more potent efflux inhibitor, in agreement with experimental data (Figure 3).

Discussion and conclusion
The utilization of dihydropyridine analogues as P-gp efflux modulators has been proposed since the early 1980s, being dexniguldipine one of the best known second-generation modulators (Hofmann et al., 1995;Reymann et al., 1993). Previous studies mostly focused on i) the identification of the binding site for dihydropyridine analogues, either experimentally (Borchers et al., 2002;Bruggemann et al., 1989;Ferry et al., 1992;Isenberg et al., 2001) or theoretically (Ranjbar et al., 2017;Shahraki et al., 2018) or ii) in the synthesis and evaluation of novel 1,4-dihydropyridine analogues Hilgeroth et al., 2013;Shekari et al., 2014;Szabon-Watola et al., 2014;Voigt et al., 2007), but quite surprisingly none of the above studies investigated the possible effect of stereochemistry in the activity of the synthesized compounds.
Herein, and in order to guide us to further develop novel DHPs as efflux inhibitors, we aimed at a thorough study on how the substitution pattern and stereochemistry may affect the MDR-reversal activity of dihydropyridine derivatives. The accumulation of R123 in the resistant cells by our DHP analogs clarified that both symmetric and asymmetric DHPs are able to inhibit R123 efflux by P-gp, with the asymmetric compound 2 having greater potency than the symmetric one (1) or verapamil.
While the initial molecular docking approach clearly identified that both symmetric and asymmetric compounds could bind P-gp either at its putative drug-binding sites or at the TMD/NBD interface (where the dihydropyridine binding site was experimentally demonstrated), only through molecular dynamics (MD) simulations we were able to fully characterize how both properties may influence DHPs interaction with Pgp. Accordingly, relative free-energies of binding suggested that, despite the higher affinities at the drug-binding M-site, the symmetric compound 1 can also favorably bind to the TMD/NBD interfaces with similar energies. Oppositelly, the (S)-enantiomer of compound 2 revealed a stronger capability of impairing conformational changes due to its higher crossinteraction capability. Several amino acids, identified as interacting with DHPs in all three sites, are described in literature as important for the drug recognition, folding and activity of P-gp, thus providing an additional support for our conclusions. Interestingly, our data additionally suggests that, while stereochemistry play a very important role on how exert DHPs its inhibitory activities on P-gp, the nature and type of substituent may affect the intracellular accumulation.
Taken together, it is our hypothesis that thiophenyl-containing symmetric DHPs have higher permeation rates through the cellular membrane, which in turn favors interactions at the cytoplasmic TMD/NBD interface, while asymmetric DHPs tend to accumulate within the membrane (due to higher energetic barriers while crossing the hydrophobic membrane core) and to diffuse into the internal P-gp cavity, where the (S)-enantiomer has more potent activities due to lower binding affinities and higher cross-interaction capabilities.
Future steps will include the synthesis and separate evaluation of both enantiomers of compound 2 to in order to corroborate the herein described results. Although chiral synthesis and separation are not trivial tasks, our results strongly show that such procedures are crucial to increase  the knowledge on DHPs as efflux inhibitors and that further studies must encompass one (or both) approaches. Of course, passive permeation studies would also benefit from evaluating stereoisomers separately. Regarding compound 1, changes on the substituent type (e.g. ring size, presence of heteroatoms or substituents within the aromatic ring or length) or in the substitution pattern of the aromatic ring at C-4 would allow a better evaluation of both passive permeation and ability to interact at the TMD/NBD interfaces.
In conclusion, we provide the first evidences of a dual-acting mode for dihydropyridines, depending on the substitution pattern and stereochemistry at position C-4. Unfortunately, our data does not allow us to fully exclude (i) if the symmetric compound can also have any meaningful activity when at the drug-binding site or, (ii) when at the cytoplasmic compartment, if asymmetric compounds could also bind at the TMD/NBD interfface and impair efflux by decoupling signal transmission between the TMD and NBD domains in a similar way as iodipine or azidopine (Isenberg et al., 2001). Furthermore, it is also conceivable that synergistic effects may occur when both enantiomers are simultaneously administered.