In contrast to the complete structure, the TMDs do not obtain stabilization according to the RMSD values, which increase and decrease over 1000 ns. Moreover, the RMSD values for G187W remain increasing from 750 to 1000 ns, and such behavior could be due to the mutation is in the cytosolic loop 1 that connects the transmembrane helix (TMH) 1 and TMH2. The NBD1 conformation in WT-MRP4 remained unstable and the RMSD values kept increasing throughout the 1000 ns of CG-MDS, while the NBD1 conformation in G187W and Y556C was stable with an RMSD value around 6.0 Å (Figure 5). The RMSD values for NBD2 of the WT-MRP4 and its variants were similar even though there was no In contrast to the complete structure, the TMDs do not obtain stabilization according to the RMSD values, which increase and decrease over 1000 ns. Moreover, the RMSD values for G187W remain increasing from 750 to 1000 ns, and such behavior could be due to the mutation is in the cytosolic loop 1 that connects the transmembrane helix (TMH) 1 and TMH2. The NBD1 conformation in WT-MRP4 remained unstable and the RMSD values kept increasing throughout the 1000 ns of CG-MDS, while the NBD1 conformation in G187W and Y556C was stable with an RMSD value around 6.0 Å (Figure 5). The RMSD values for NBD2 of the WT-MRP4 and its variants were similar even though there was no stabilization through the simulation. In addition, the RMSD values of G187W and Y556C tended to increase while the RMSD values of WT-MRP4 tended to decrease at the end of the simulation. ues for G187W remain increasing from 750 to 1000 ns, and such behavior could be due to the mutation is in the cytosolic loop 1 that connects the transmembrane helix (TMH) 1 and TMH2. The NBD1 conformation in WT-MRP4 remained unstable and the RMSD values kept increasing throughout the 1000 ns of CG-MDS, while the NBD1 conformation in G187W and Y556C was stable with an RMSD value around 6.0 Å (Figure 5). The RMSD values for NBD2 of the WT-MRP4 and its variants were similar even though there was no stabilization through the simulation. In addition, the RMSD values of G187W and Y556C tended to increase while the RMSD values of WT-MRP4 tended to decrease at the end of the simulation.

Figure 5. RMSD plot for NBD1 and NBD2 of WT-MRP4 and its variants throughout the 1000 ns of MDS. According to the RMSD plot in Figure 6, the conformation in residues (r) 85 to 236 (r85-236 site corresponding to TMH1-TMH4) of WT-MRP4 remained stable from 250 to 1000 ns while in G187W it stabilized at the last 250 ns. r85-236 site in Y556C did not stabilize throughout the 1000 ns. In the same Figure 6, the RMSD plot indicates that the conformation in the r715-866 site (TMH7-TMH10) in WT-MRP4 and Y556C kept constant, with a tendency toward increasing motion in Y556C and decreasing motion in WT-MRP4. Besides this, the conformation in r715-866 site in G187W did not stabilize and the RMSD

values suggest high motion throughout the 1000 ns of simulation. In the case of ATP sites 1 and 2 (Figure 7), the conformation remained unstable, suggesting that it is a site with high motion, with exception of ATP site 2 of WT-MRP4, which kept stable through the simulation. Regarding the three different MRP4 structures, alignments on the NBDs, TMDs, and TMH1-TMH4 and TMH7-TMH10 sites were performed to determine differences among the structures. In all the alignments (Supplementary Material S1), it was observed that all the sites studied presented significant structural differences, according to the RMSD values, comparing WT-MRP4 vs. its variants and G187W vs. Y566C. Moreover, the TMHs are responsible for the specificity for the substrate, and the r-85-236 and r715-866 sites in WT-MRP4 were significantly different with respect to G187W and Y556C, which could lead to differences in the ligand affinity, the ligand binding site, and the motion of the protein [23,24]. high motion, with exception of ATP site 2 of WT-MRP4, which kept stable through the simulation. Regarding the three different MRP4 structures, alignments on the NBDs, TMDs, and TMH1-TMH4 and TMH7-TMH10 sites were performed to determine differences among the structures. In all the alignments (Supplementary Material S1), it was observed that all the sites studied presented significant structural differences, according to the RMSD values, comparing WT-MRP4 vs. its variants and G187W vs. Y566C. Moreover, the TMHs are responsible for the specificity for the substrate, and the r-85-236 and r715- 866 sites in WT-MRP4 were significantly different with respect to G187W and Y556C, which could lead to differences in the ligand affinity, the ligand binding site, and the motion of the protein [23,24]. simulation. Regarding the three different MRP4 structures, alignments on the NBDs, TMDs, and TMH1-TMH4 and TMH7-TMH10 sites were performed to determine differences among the structures. In all the alignments (Supplementary Material S1), it was observed that all the sites studied presented significant structural differences, according to the RMSD values, comparing WT-MRP4 vs. its variants and G187W vs. Y566C. Moreover, the TMHs are responsible for the specificity for the substrate, and the r-85-236 and r715- 866 sites in WT-MRP4 were significantly different with respect to G187W and Y556C, which could lead to differences in the ligand affinity, the ligand binding site, and the motion of the protein [23,24].

Molecules 2021, 26, x FOR PEER REVIEW 6 of 22

Molecules 2021, 26, x FOR PEER REVIEW 6 of 22

According to the RMSD plot in Figure 6, the conformation in residues (r) 85 to 236

According to the RMSD plot in Figure 6, the conformation in residues (r) 85 to 236

(r85-236 site corresponding to TMH1-TMH4) of WT-MRP4 remained stable from 250 to 1000 ns while in G187W it stabilized at the last 250 ns. r85-236 site in Y556C did not stabilize throughout the 1000 ns. In the same Figure 6, the RMSD plot indicates that the conformation in the r715-866 site (TMH7-TMH10) in WT-MRP4 and Y556C kept constant, with a tendency toward increasing motion in Y556C and decreasing motion in WT-MRP4. Besides this, the conformation in r715-866 site in G187W did not stabilize and the RMSD values suggest high motion throughout the 1000 ns of simulation. In the case of ATP sites 1 and 2 (Figure 7), the conformation remained unstable, suggesting that it is a site with

(r85-236 site corresponding to TMH1-TMH4) of WT-MRP4 remained stable from 250 to 1000 ns while in G187W it stabilized at the last 250 ns. r85-236 site in Y556C did not stabilize throughout the 1000 ns. In the same Figure 6, the RMSD plot indicates that the conformation in the r715-866 site (TMH7-TMH10) in WT-MRP4 and Y556C kept constant, with a tendency toward increasing motion in Y556C and decreasing motion in WT-MRP4. Besides this, the conformation in r715-866 site in G187W did not stabilize and the RMSD values suggest high motion throughout the 1000 ns of simulation. In the case of ATP sites 1 and 2 (Figure 7), the conformation remained unstable, suggesting that it is a site with high motion, with exception of ATP site 2 of WT-MRP4, which kept stable through the

Figure 6. RMSD plot for r85-236 and r715-866 of WT-MRP4 and its variants. r86-236 represents TMH1-TMH4 and r715-866 represents TMH7-TMH10. **Figure 6.** RMSD plot for r85-236 and r715-866 of WT-MRP4 and its variants. r86-236 represents TMH1-TMH4 and r715-866 represents TMH7-TMH10. TMH1-TMH4 and r715-866 represents TMH7-TMH10.

Figure 7. RMSD plot for ATP binding site 1 (ATP site 1) and ATP binding site 2 (ATP site 2) of WT-MRP4 and its variants. Figure 7. RMSD plot for ATP binding site 1 (ATP site 1) and ATP binding site 2 (ATP site 2) of WT-MRP4 and its variants. **Figure 7.** RMSD plot for ATP binding site 1 (ATP site 1) and ATP binding site 2 (ATP site 2) of WT-MRP4 and its variants.

#### *2.2. Molecular Docking in MRP4 and Variants*

2.2. Molecular Docking in MRP4 and Variants By using cluster 1 from CG-MDS for each structure of MRP4, molecular docking was performed to explore the effect of the MRP4 variants on the affinity of three different groups of molecules, previously reported as substrates or inhibitors in vitro. Table 1 presents the docking score (DS), expressed as kcal/mol, related to the interaction between 2.2. Molecular Docking in MRP4 and Variants By using cluster 1 from CG-MDS for each structure of MRP4, molecular docking was performed to explore the effect of the MRP4 variants on the affinity of three different groups of molecules, previously reported as substrates or inhibitors in vitro. Table 1 presents the docking score (DS), expressed as kcal/mol, related to the interaction between By using cluster 1 from CG-MDS for each structure of MRP4, molecular docking was performed to explore the effect of the MRP4 variants on the affinity of three different groups of molecules, previously reported as substrates or inhibitors in vitro. Table 1 presents the docking score (DS), expressed as kcal/mol, related to the interaction between endogenous substrates and WT-MRP4, Y556C, and G187W. In this work, significant differences between docking poses were considered when a difference greater than 1 kcal/mol in DS was

present. According to the DS, most endogen substrates significantly changed their DS in G187W and Y556C with respect to WT-MRP4. The Y556C mutation is located at NBD1 and leads to a different conformation with respect to WT-MRP4, which causes the most substantial change in the ligand binding site of all the molecules studied, even more than those molecules with significant changes in DS compared to WT-MRP4.

**Table 1.** Endogen substrate DS (kcal/mol) in the WT-MRP4 and MRP4 variants.

cAMP: cyclic adenosine monophosphate. \* Represents more than 1 kcal/mol of difference in the DS compared with WT.

Since cAMP is considered to be the main molecule transported by MRP4 [25], it was used as a control to compare the effect of nsSNPs and ATP binding. The cAMP DS was not considered significantly different in WT-MRP4 with respect to Y556C and G187W. Cholic acid DS was significantly different, with more than 3 kcal/mol when comparing WT-MRP4 with respect to G187W, while cholic acid DS in WT-MRP4 with respect to Y556C was significantly different at over 2 kcal/mol. The taurocholic acid DS was significantly different, with more than 4 kcal/mol when comparing WT-MRP4 with respect to G187W. Folic acid (FA) was the molecule with the highest variation in DS when comparing WT-MRP4 to G187W and Y556C, with 7.05 and 6.17 kcal/mol differences.

Table 2 presents the DSs of drug substrates. Cefazoline and olmesartan DSs were significantly different in WT-MRP4 with respect to G187W, while cefazoline, furosemide, leucovorin, methotrexate, tenofovir, and topotecan DSs were significantly different in WT-MRP4 with respect to Y556C. According to DS, most drug substrates could present more affinity for Y556C than for WT-MRP4.

**Table 2.** Drug substrate DS (kcal/mol) in the WT-MRP4 and MRP4 variants.

\* Represents more than 1 kcal/mol difference in the DS compared against WT.

The DS of inhibitors is presented in Table 3. Glafenine DS was significantly different in WT-MRP4 with respect to both G187W and Y556C. Ceefourin-1, indomethacin, and sildenafil DSs were significantly different in WT-MRP4 with respect to G187W, while losartan DS was significantly different between WT-MRP4 and Y556C. Endogen substrates were the group of molecules with more variation in the DS, and thus such mutations on the ABCC4 gene could lead to changes in cell metabolism related to changes in the distribution of endogen substrates across cell membranes. Several mutational studies on MRP1 have demonstrated that amino acids in several TMDs are involved in substrate binding and nsSNPs can modify the ligand binding site or the affinity of ligands in a selective manner [11]. Yet, the DS values are not totally related to function. Regarding drug substrates and inhibitors, a small number of molecules changed their DS significantly to G187W and Y556C as related to WT-MRP4, which could modify the pattern of transport of cefazoline, ceftizoxime, olmesartan, topotecan, and ceefourin-1, which could, possibly, mean that they act as substrates or inhibitors depending on the MRP4 variant.

**Table 3.** Drug inhibitor DS (kcal/mol) in the WT-MRP4 and MRP4 variants.

\* Represents more than 1 kcal/mol difference in the DS compared against WT.

#### *2.3. Differences in the Interaction Pattern in MRP4 Structures*

Ligand interaction diagrams (LIDs) represent the pattern of intermolecular interactions of molecules with MRP4 amino acids. Those molecules with a >2 kcal/mol difference in the three different MRP4 structures appear in LIDs in Supplementary Material S2–S5 and the LIDs for FA appear in Figures 8–10. The interaction sites in the three different MRP4 structures were different, mainly in Y556C, for all the molecules exhibited in the LIDs. Cholic acid interacts mainly with hydrophobic residues on WT-MRP4 and Y556C, while on G187W it interacts with hydrophobic and polar residues. H-bonds are only exerted through hydrophobic residues, except on Y556C, where arginine exerts an H-bond with cholic acid. In addition, cholic acid interacts with positively charged amino acids in the three MRP4 structures, lysine on G187W and WT-MRP4, and arginine on Y556C (Supplementary Material S2). The taurocholic acid binding site was different in WT-MRP4 with respect to G187W and was totally different with respect to Y556C. The intermolecular interactions in Y556C and G187W were H-bonds, polar, hydrophobic, interactions with positively and negatively charged amino acids, while interactions in WT-MRP4 were hydrophobic and polar but did not interact with negatively charged amino acids and did not exert H-bonds (Supplementary Material S3). The cefazoline binding site was quite different on the three different MRP4 structures. Hydrophobic and polar residues on WT-MRP4, Y556C, and G187W interact with cefazoline, but only on G187W do H-bonds with arginine occur. WT-MRP4 and Y556C interact with cefazoline through negatively charged amino acids (glutamate), while the G187W interacts through positively charged amino acids (arginine) (Supplementary Material S4). Ceefourin-1 interacts with hydrophobic, polar, negative, and positively charged residues on WT-MRP4 and Y556C, even though the binding site is different in each MRP4 structure, which leads to only one difference; ceefourin-1 interacts via H-bonds in Y556C. Additionally, ceefourin-1 interacts in G187W

with hydrophobic, polar, and positively charged residues, and via H-bonds and π-π stacking (Supplementary Material S5). The Y556C mutation led to a substantial change in the binding site of most of the molecules analyzed, while G187W mutation did not substantially change the binding site with respect to WT-MRP4. Although the molecule's binding site was different among the MRP4, the pattern of intermolecular interactions could be similar to that observed in the LIDs. The molecules presented in the LIDs had a substantial modification to their binding site and intermolecular interactions in Y556C, which could be related to changes in the transport-rate, IC50, entrance to the inner cavity, and effect on the Y556C conformational movement. did not change its binding site throughout all the AA-MD simulation, but it kept moving throughout 25 ns to obtain a DS greater than that of the most stable conformation at 8 ns (Supplementary Material S12–S17). The G187W mutation is located at cytosolic loop 1, close to the entrance of the inner cavity [26], and leads to changes in the MRP4 conformation which could block or interfere with ligand binding or the entrance to its binding site. Notwithstanding this, it is not possible to determine with this study whether there is a relevant effect on the FA transport by G187W. Figure 9 shows the FA LID in G187W at T0 and 8.0 ns in AA-MDS. The FA binding site was slightly different, with intermolecular interactions by an H-bond and π-cation at T0, while at 8.0 no intermolecular interactions were observed, but the pocket binding was composed of hydrophobic amino acids in both time steps. did not change its binding site throughout all the AA-MD simulation, but it kept moving throughout 25 ns to obtain a DS greater than that of the most stable conformation at 8 ns (Supplementary Material S12–S17). The G187W mutation is located at cytosolic loop 1, close to the entrance of the inner cavity [26], and leads to changes in the MRP4 conformation which could block or interfere with ligand binding or the entrance to its binding site. Notwithstanding this, it is not possible to determine with this study whether there is a relevant effect on the FA transport by G187W. Figure 9 shows the FA LID in G187W at T0 and 8.0 ns in AA-MDS. The FA binding site was slightly different, with intermolecular interactions by an H-bond and π-cation at T0, while at 8.0 no intermolecular interactions were observed, but the pocket binding was composed of hydrophobic amino acids in both time steps.

S6–S23). As mentioned above, FA did not change its conformation significantly in WT-MRP4 from 0 to 25 ns, while in G187W the conformation was significantly different at 20 and 25 ns, as the carboxyl group was responsible for the FA conformational changes and multiple intermolecular interactions, such as H-bonds, π-π, and π-cation. FA in G187W

S6–S23). As mentioned above, FA did not change its conformation significantly in WT-MRP4 from 0 to 25 ns, while in G187W the conformation was significantly different at 20 and 25 ns, as the carboxyl group was responsible for the FA conformational changes and multiple intermolecular interactions, such as H-bonds, π-π, and π-cation. FA in G187W

Molecules 2021, 26, x FOR PEER REVIEW 10 of 22

Molecules 2021, 26, x FOR PEER REVIEW 10 of 22

Figure 8. FA binding interaction diagram in WT-MRP4 at T0 (A) and 16.20 ns (B) in AA-MDS. The nomenclature for the intermolecular interactions is shown in the bottom of the figure. **Figure 8.** FA binding interaction diagram in WT-MRP4 at T0 (**A**) and 16.20 ns (**B**) in AA-MDS. The nomenclature for the intermolecular interactions is shown in the bottom of the figure. Figure 8. FA binding interaction diagram in WT-MRP4 at T0 (A) and 16.20 ns (B) in AA-MDS. The nomenclature for the intermolecular interactions is shown in the bottom of the figure.

Figure 9. FA binding diagram interactions on G187W at T0 (A) and at 8.0 ns (B) in AA-MDS. **Figure 9.** FA binding diagram interactions on G187W at T0 (**A**) and at 8.0 ns (**B**) in AA-MDS. different across the AA-MDS (Supplementary Material S18–S23).

Figure 9. FA binding diagram interactions on G187W at T0 (A) and at 8.0 ns (B) in AA-MDS.

Figure 10 shows the FA LID in Y556C at T0 and 15.00 ns in AA-MDS. In this case, the

Figure 10. FA binding interaction diagrams in Y556C at T0 (A) and 15.0 ns (B) in AA-MDS. **Figure 10.** FA binding interaction diagrams in Y556C at T0 (**A**) and 15.0 ns (**B**) in AA-MDS.

The AA-MDS studies were conducted in cAMP as a substrate control to determine if

First, AA-MDS were performed on the wild type and variants for MRP4-cAMP complex to obtain the C1 conformations. Further, ATP and Mg2+ were added into the MRP4 structures to obtain a more complete system, and further AA-MDS were performed to determine a new C1. This procedure was applied in the MRP4-FA complexes as well. Both C1 conformations before ATP and after ATP were used for umbrella sampling studies to obtain the free binding energy (ΔG) of each ligand. Since two ATP molecules bind to MRP4, we referred to them as ATP1 and ATP2 for the NBD1 and NBD2, respectively. In all the MRP4-FA and MRP4-cAMP complexes with ATP bonds, the binding site and the ligand conformation were slightly different concerning the MRP4-FA and MRP4-cAMP complexes without ATP. The presence or the absence of ATP modifies the pattern of intermolecular interactions as well. The H-bonds, π-π stacking, and π-cation interactions are present together or individually in all the complexes, except in G187W-FA without ATP. Given the fact that the ATP influences the ligand conformation, it modifies the intermolecular interactions as well. Figures 11 and 12 show the ATP binding site in NBD1 and

Figure 11. Representation of the ATP binding site 1 in the complete MRP4 structure and an amplification of such site. Color pink in the complete MRP4 structure represents the ATP site 1. Color red

NBD2 from the WT-MRP4, G187W, and Y556C structures.

Table 4 shows the most important residues required for the interactions of endogen substrates, drug substrates, and inhibitors in each MRP4 structure. It is worth noting that in the inhibitor group, only Lys 329 and Arg 951 appear in both endogen and drug substrates, suggesting that such residues may play an important role in the binding site. In the endogen and drug substrate groups, at least two residues are repeated—Arg 946 and Lys 702.

**Table 4.** Residues considered as important for ligand–MRP4 interactions.

#### *2.4. All-Atom Molecular Dinaymics (AA-MD) Simulations and Umbrella Sampling Studies*

FA was the molecule with the highest DS variation in G187W and Y556C related to WT-MRP4 in molecular docking studies, considering that ATP was not bound; hence, to study the effect of the presence of ATP and the mutations, changes in the pattern of intermolecular interactions and affinity to MRP4 structures, 25 ns AA-MDS, and 10 ns of Umbrella sampling simulations were carried out on the MRP4-FA complexes and compared with cAMP as a control molecule. The C1 in AA-MDS was at 16.2 ns in WT-MRP4, 8.0 ns in G187W, and 15.0 ns in Y556C. Figure 8 shows the FA LID in WT-MRP4 at T0 and 16.2 ns. In this simulation, the differences in the FA binding site and the pattern of intermolecular interaction can be observed according to the WT-MRP4 conformation at a given time. The FA binding site was the same at T0 and 16.2 ns, suggesting that the protein can be in an inward-facing conformation for 25 ns or even more. Moreover, the pattern of intermolecular interactions between FA and WT-MRP4 is quite similar, consisting of H-bonds; π-π stacking; and interactions with polar, hydrophobic, and negatively and positively charged residues. The differences in the pattern of intermolecular interactions rely on the H-bonds, with 4 H-bonds at T0 with Phe698 and Leu835, Glu1002 and Thr839. Meanwhile, at 16.2 ns the interactions were mainly hydrophobic (π-π stacking) and there was one H-bond with Glu1002. Therefore, it seems that FA from 0 to 16.2 ns interacts with MRP4 to achieve its optimal bonding only. A longer simulation will help to determine all the FA binding sites across the WT-MRP4. Besides this, the AA-MDS studies were performed on each MRP4-FA complex to determine differences in the patterns of interactions and conformations in frames from 0 to 25 ns every 5 ns referred to T0 (Supplementary Material S6–S23). As mentioned above, FA did not change its conformation significantly in WT-MRP4 from 0 to 25 ns, while in G187W the conformation was significantly different at 20 and 25 ns, as the carboxyl group was responsible for the FA conformational changes and multiple intermolecular interactions, such as H-bonds, π-π, and π-cation. FA in G187W did not change its binding site throughout all the AA-MD simulation, but it kept moving throughout 25 ns to obtain a DS greater than that of the most stable conformation at 8 ns (Supplementary Material S12–S17). The G187W mutation is located at cytosolic loop 1, close to the entrance of the inner cavity [26], and leads to changes in the MRP4 conformation which could block or interfere with ligand binding or the entrance to its binding site. Notwithstanding this, it is not possible to determine with this study whether there is a relevant effect on the FA transport by G187W. Figure 9 shows the FA LID in G187W at T0 and 8.0 ns in AA-MDS. The FA binding site was slightly different, with intermolecular

interactions by an H-bond and π-cation at T0, while at 8.0 no intermolecular interactions were observed, but the pocket binding was composed of hydrophobic amino acids in both time steps. Figure 10. FA binding interaction diagrams in Y556C at T0 (A) and 15.0 ns (B) in AA-MDS.

Molecules 2021, 26, x FOR PEER REVIEW 11 of 22

different across the AA-MDS (Supplementary Material S18–S23).

contrary, the pattern of intermolecular interactions as well as the FA conformation was

Figure 10 shows the FA LID in Y556C at T0 and 15.00 ns in AA-MDS. In this case, the FA binding site was the same, but the intermolecular interactions by H-bonds at 15.00 ns were higher than at T0, which seems to confer to FA a stronger binding over time. Y556C promoted the greatest change in the FA binding site with respect to WT-MRP4 and G187W, but throughout the 25 ns in AA-MDS the binding site did not change. On the contrary, the pattern of intermolecular interactions as well as the FA conformation was different across the AA-MDS (Supplementary Material S18–S23). The AA-MDS studies were conducted in cAMP as a substrate control to determine if the nsSNPs affect the pattern of intermolecular interactions in the same manner as FA. First, AA-MDS were performed on the wild type and variants for MRP4-cAMP complex to obtain the C1 conformations. Further, ATP and Mg2+ were added into the MRP4 structures to obtain a more complete system, and further AA-MDS were performed to deter-

The AA-MDS studies were conducted in cAMP as a substrate control to determine if the nsSNPs affect the pattern of intermolecular interactions in the same manner as FA. First, AA-MDS were performed on the wild type and variants for MRP4-cAMP complex to obtain the C1 conformations. Further, ATP and Mg2+ were added into the MRP4 structures to obtain a more complete system, and further AA-MDS were performed to determine a new C1. This procedure was applied in the MRP4-FA complexes as well. Both C1 conformations before ATP and after ATP were used for umbrella sampling studies to obtain the free binding energy (∆G) of each ligand. Since two ATP molecules bind to MRP4, we referred to them as ATP1 and ATP2 for the NBD1 and NBD2, respectively. In all the MRP4-FA and MRP4-cAMP complexes with ATP bonds, the binding site and the ligand conformation were slightly different concerning the MRP4-FA and MRP4-cAMP complexes without ATP. The presence or the absence of ATP modifies the pattern of intermolecular interactions as well. The H-bonds, π-π stacking, and π-cation interactions are present together or individually in all the complexes, except in G187W-FA without ATP. Given the fact that the ATP influences the ligand conformation, it modifies the intermolecular interactions as well. Figures 11 and 12 show the ATP binding site in NBD1 and NBD2 from the WT-MRP4, G187W, and Y556C structures. mine a new C1. This procedure was applied in the MRP4-FA complexes as well. Both C1 conformations before ATP and after ATP were used for umbrella sampling studies to obtain the free binding energy (ΔG) of each ligand. Since two ATP molecules bind to MRP4, we referred to them as ATP1 and ATP2 for the NBD1 and NBD2, respectively. In all the MRP4-FA and MRP4-cAMP complexes with ATP bonds, the binding site and the ligand conformation were slightly different concerning the MRP4-FA and MRP4-cAMP complexes without ATP. The presence or the absence of ATP modifies the pattern of intermolecular interactions as well. The H-bonds, π-π stacking, and π-cation interactions are present together or individually in all the complexes, except in G187W-FA without ATP. Given the fact that the ATP influences the ligand conformation, it modifies the intermolecular interactions as well. Figures 11 and 12 show the ATP binding site in NBD1 and NBD2 from the WT-MRP4, G187W, and Y556C structures.

Figure 11. Representation of the ATP binding site 1 in the complete MRP4 structure and an amplification of such site. Color pink in the complete MRP4 structure represents the ATP site 1. Color red **Figure 11.** Representation of the ATP binding site 1 in the complete MRP4 structure and an amplification of such site. Color pink in the complete MRP4 structure represents the ATP site 1. Color red and represent the residues interacting with ATP, and rainbow colors represent ATP. The yellow dotted lines indicate polar interactions.

Molecules 2021, 26, x FOR PEER REVIEW 12 of 22

dotted lines indicate polar interactions.

and represent the residues interacting with ATP, and rainbow colors represent ATP. The yellow

Figure 12. Representation of the ATP binding site 2 in the complete MRP4 structure and an amplification of such site. Color cyan in the complete MRP4 structure represents the ATP site 2. Color red represent the residues interacting with ATP, and rainbow colors represent ATP. The yellow dotted lines indicate polar interactions. **Figure 12.** Representation of the ATP binding site 2 in the complete MRP4 structure and an amplification of such site. Color cyan in the complete MRP4 structure represents the ATP site 2. Color red represent the residues interacting with ATP, and rainbow colors represent ATP. The yellow dotted lines indicate polar interactions.

Figures 13–15 show the FA binding sites in WT-MRP4, G187W, and Y556C, respectively. All the afore mentioned figures contain the full MRP4 structure and a close-up of the ligand binding site. As observed in Figures 13 and 14, the FA binding site is almost the same and it is surrounded by the TMH1-TMH3; meanwhile, Figure 15 shows that the FA binding site in Y556C is surrounded by TMH2-TMH7 deep in the inner cavity. This finding is interesting, due to the fact that the Y556C mutation could promote the greater affinity of FA for the Y556C variant by switching into a "high affinity inward-facing orienta-Figures 13–15 show the FA binding sites in WT-MRP4, G187W, and Y556C, respectively. All the afore mentioned figures contain the full MRP4 structure and a close-up of the ligand binding site. As observed in Figures 13 and 14, the FA binding site is almost the same and it is surrounded by the TMH1-TMH3; meanwhile, Figure 15 shows that the FA binding site in Y556C is surrounded by TMH2-TMH7 deep in the inner cavity. This finding is interesting, due to the fact that the Y556C mutation could promote the greater affinity of FA for the Y556C variant by switching into a "high affinity inward-facing orientation", making the binding of FA easier [27]. Nevertheless, this does not mean that the FA can be transported properly, since MRP4 has an altered conformation within the NBD2, which blocks the ATP hydrolysis but not ATP binding. However, MRP4 can function with the energy provided by one ATP hydrolysis, but the Y556C mutation could interfere with such ATP hydrolysis [28].

tion", making the binding of FA easier [27]. Nevertheless, this does not mean that the FA can be transported properly, since MRP4 has an altered conformation within the NBD2, which blocks the ATP hydrolysis but not ATP binding. However, MRP4 can function with the energy provided by one ATP hydrolysis, but the Y556C mutation could interfere with such ATP hydrolysis [28]. Figure 16 shows that the cAMP binding site in Y556C is quite different compared to WT-MRP4 and G187W. The cAMP is located out of the normal binding site out of the TMDs and it may not be transported at a proper rate or efficacy. It seems that Y556C is in the "low-affinity outward-facing orientation", where the main characteristic in this stage is that the affinity of the transported entity switches from high affinity (low chemical potential of the substrate) to low affinity (high chemical potential). Such change in affinity is due to the gate to the inside is closed, and the gate to the outside of the membrane is opened [23]. The global effect will depend on the cell type. Focusing on leukemia cells, could be harmful to them, because the intracellular cAMP levels would be increased and thus, leading to apoptosis [25]. On the other hand, the cAMP binding site in Y556C could represent the initial interaction only and further re-location such as the WT-MTP4. Longer AA-MDS are required to test this hypothesis.

To contrast the effect of nsSNPs in MRP4, AA-MDS and umbrella sampling studies were performed in ceefourin-1, which is possesses a high selectivity over MRP4 inhibition. It has been suggested that ceefourin-1 may act as a competitive inhibitor in MRP4 [7], which is consistent to the AA-MDS where the binding site for cAMP and ceefourin-1 is similar

56

Figure 13. FA binding site in WT-MRP4 and amplification of such site. Rainbow colors in both the complete MRP4 structure and amplification of the binding site represent the FA structure while the

in WT-MRP4, G187W and Y556C with ATP in the complex (Figures 17–19). Interestingly, ceefourin-1 binds at the same abnormal binding site as cAMP in Y556C-ATP. Moreover, Y556C seems to be in the same "low-affinity outward-facing orientation" where the low affinity (high chemical potential) of the ligand leads to switch in the binding site [23]. As mentioned before, such abnormal binding site could be the initial binding site but longer AA-MDS is required to demonstrate it. tion", making the binding of FA easier [27]. Nevertheless, this does not mean that the FA can be transported properly, since MRP4 has an altered conformation within the NBD2, which blocks the ATP hydrolysis but not ATP binding. However, MRP4 can function with the energy provided by one ATP hydrolysis, but the Y556C mutation could interfere with such ATP hydrolysis [28].

Molecules 2021, 26, x FOR PEER REVIEW 12 of 22

dotted lines indicate polar interactions.

lines indicate polar interactions.

and represent the residues interacting with ATP, and rainbow colors represent ATP. The yellow

Figure 12. Representation of the ATP binding site 2 in the complete MRP4 structure and an amplification of such site. Color cyan in the complete MRP4 structure represents the ATP site 2. Color red represent the residues interacting with ATP, and rainbow colors represent ATP. The yellow dotted

Figures 13–15 show the FA binding sites in WT-MRP4, G187W, and Y556C, respectively. All the afore mentioned figures contain the full MRP4 structure and a close-up of the ligand binding site. As observed in Figures 13 and 14, the FA binding site is almost the same and it is surrounded by the TMH1-TMH3; meanwhile, Figure 15 shows that the FA binding site in Y556C is surrounded by TMH2-TMH7 deep in the inner cavity. This finding is interesting, due to the fact that the Y556C mutation could promote the greater affinity of FA for the Y556C variant by switching into a "high affinity inward-facing orienta-

Figure 13. FA binding site in WT-MRP4 and amplification of such site. Rainbow colors in both the complete MRP4 structure and amplification of the binding site represent the FA structure while the **Figure 13.** FA binding site in WT-MRP4 and amplification of such site. Rainbow colors in both the complete MRP4 structure and amplification of the binding site represent the FA structure while the yellow spheres and stick represent the ATP. The yellow dotted lines indicate polar interactions. The LID of WT-MRP4-FA is presented at the bottom. yellow spheres and stick represent the ATP. The yellow dotted lines indicate polar interactions. The LID of WT-MRP4-FA is presented at the bottom.

Figure 14. FA binding site in G187W and amplification of such site. Rainbow colors in both the complete WT-MRP4 and amplification of the binding site represent the FA structure while the yellow spheres and stick represent the ATP. The yellow dotted lines indicate polar interactions. The LID of G187W-FA is presented at the bottom. **Figure 14.** FA binding site in G187W and amplification of such site. Rainbow colors in both the complete WT-MRP4 and amplification of the binding site represent the FA structure while the yellow spheres and stick represent the ATP. The yellow dotted lines indicate polar interactions. The LID of G187W-FA is presented at the bottom.

Figure 15. FA binding site in Y556C and amplification of such a site. Rainbow colors in both the complete G187W and amplification of the binding site represent the FA structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar interactions. The LID of Y556C-

Figure 16 shows that the cAMP binding site in Y556C is quite different compared to WT-MRP4 and G187W. The cAMP is located out of the normal binding site out of the TMDs and it may not be transported at a proper rate or efficacy. It seems that Y556C is in the "low-affinity outward-facing orientation", where the main characteristic in this stage is that the affinity of the transported entity switches from high affinity (low chemical potential of the substrate) to low affinity (high chemical potential). Such change in affinity is due to the gate to the inside is closed, and the gate to the outside of the membrane is opened [23]. The global effect will depend on the cell type. Focusing on leukemia cells, could be harmful to them, because the intracellular cAMP levels would be increased and thus, leading to apoptosis [25]. On the other hand, the cAMP binding site in Y556C could

FA is presented at the bottom.

Figure 14. FA binding site in G187W and amplification of such site. Rainbow colors in both the complete WT-MRP4 and amplification of the binding site represent the FA structure while the yellow spheres and stick represent the ATP. The yellow dotted lines indicate polar interactions. The

yellow spheres and stick represent the ATP. The yellow dotted lines indicate polar interactions. The

Molecules 2021, 26, x FOR PEER REVIEW 13 of 22

LID of WT-MRP4-FA is presented at the bottom.

LID of G187W-FA is presented at the bottom.

Figure 15. FA binding site in Y556C and amplification of such a site. Rainbow colors in both the complete G187W and amplification of the binding site represent the FA structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar interactions. The LID of Y556C-FA is presented at the bottom. **Figure 15.** FA binding site in Y556C and amplification of such a site. Rainbow colors in both the complete G187W and amplification of the binding site represent the FA structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar interactions. The LID of Y556C-FA is presented at the bottom. represent the initial interaction only and further re-location such as the WT-MTP4. Longer AA-MDS are required to test this hypothesis.

Figure 16. cAMP binding site in Y556C and amplification of such site. Rainbow colors in both the complete Y556C and amplification of the binding site represent the cAMP structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar interactions. The LID of Y556C-**Figure 16.** cAMP binding site in Y556C and amplification of such site. Rainbow colors in both the complete Y556C and amplification of the binding site represent the cAMP structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar interactions. The LID of Y556CcAMP is presented at the bottom.

site but longer AA-MDS is required to demonstrate it.

To contrast the effect of nsSNPs in MRP4, AA-MDS and umbrella sampling studies

The cAMP and ceefourin-1 abnormal binding site in Y556C was one of the most re-

markable findings in this work and showed the importance of studying the global effect

were performed in ceefourin-1, which is possesses a high selectivity over MRP4 inhibition. It has been suggested that ceefourin-1 may act as a competitive inhibitor in MRP4 [7], which is consistent to the AA-MDS where the binding site for cAMP and ceefourin-1 is similar in WT-MRP4, G187W and Y556C with ATP in the complex (Figures 17–19). Interestingly, ceefourin-1 binds at the same abnormal binding site as cAMP in Y556C-ATP. Moreover, Y556C seems to be in the same "low-affinity outward-facing orientation" where the low affinity (high chemical potential) of the ligand leads to switch in the binding site [23]. As mentioned before, such abnormal binding site could be the initial binding

of Y556C in different cell lines.

site but longer AA-MDS is required to demonstrate it.

Molecules 2021, 26, x FOR PEER REVIEW 14 of 22

AA-MDS are required to test this hypothesis.

cAMP is presented at the bottom.

of Y556C in different cell lines.

represent the initial interaction only and further re-location such as the WT-MTP4. Longer

Figure 16. cAMP binding site in Y556C and amplification of such site. Rainbow colors in both the complete Y556C and amplification of the binding site represent the cAMP structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar interactions. The LID of Y556C-

To contrast the effect of nsSNPs in MRP4, AA-MDS and umbrella sampling studies

The cAMP and ceefourin-1 abnormal binding site in Y556C was one of the most re-

markable findings in this work and showed the importance of studying the global effect

were performed in ceefourin-1, which is possesses a high selectivity over MRP4 inhibition. It has been suggested that ceefourin-1 may act as a competitive inhibitor in MRP4 [7], which is consistent to the AA-MDS where the binding site for cAMP and ceefourin-1 is similar in WT-MRP4, G187W and Y556C with ATP in the complex (Figures 17–19). Interestingly, ceefourin-1 binds at the same abnormal binding site as cAMP in Y556C-ATP. Moreover, Y556C seems to be in the same "low-affinity outward-facing orientation" where the low affinity (high chemical potential) of the ligand leads to switch in the binding site [23]. As mentioned before, such abnormal binding site could be the initial binding

**Figure 17.** Ceefourin-1 binding site in WT-MRP4 and amplification of such site. Rainbow colors in both the complete WT-MRP4 and amplification of the binding site represent the ceefourin-1 structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar interactions. The LID of WT-MRP4-ceefourin-1 is presented at the bottom. Figure 17. Ceefourin-1 binding site in WT-MRP4 and amplification of such site. Rainbow colors in both the complete WT-MRP4 and amplification of the binding site represent the ceefourin-1 structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar interactions. The LID of WT-MRP4-ceefourin-1 is presented at the bottom.

Figure 18. Ceefourin-1 binding site in G187W and the amplification of such a site. Rainbow colors in both the complete G187W and amplification of the binding site represent the ceefourin-1 structure, while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar interactions. The LID of G187W-ceefourin-1 is presented at the bottom. **Figure 18.** Ceefourin-1 binding site in G187W and the amplification of such a site. Rainbow colors in both the complete G187W and amplification of the binding site represent the ceefourin-1 structure, while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar interactions. The LID of G187W-ceefourin-1 is presented at the bottom.

teractions. The LID of Y556C-ceefourin-1 is presented at the bottom.

Figure 19. Ceefourin-1 binding site in Y556C and the amplification of such a site. Rainbow colors in both the complete Y556C and amplification of the binding site represent the ceefourin-1 structure, while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar in-

Umbrella sampling simulations were performed to determine the theoretical affinity,

expressed as the free energy of binding (ΔG), of an event related to protein–ligand interactions along a reaction coordinate. Table 5 shows the ΔG obtained by Umbrella sampling through 10 ns of simulation on each complex. According to Umbrella sampling results, the Y556C-FA complex is the most stable with or without ATP bound to MRP4. Interestingly, ΔG for WT-FA and Y556C-FA complexes were lower in those MRP4 without ATP. In addition, ΔGs for G187W-cAMP and Y556C-cAMP complexes were higher without ATP, suggesting that the ATP bound is required for both FA and cAMP binding

Figure 18. Ceefourin-1 binding site in G187W and the amplification of such a site. Rainbow colors in both the complete G187W and amplification of the binding site represent the ceefourin-1 structure, while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar

Figure 17. Ceefourin-1 binding site in WT-MRP4 and amplification of such site. Rainbow colors in both the complete WT-MRP4 and amplification of the binding site represent the ceefourin-1 structure while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar

interactions. The LID of G187W-ceefourin-1 is presented at the bottom.

Molecules 2021, 26, x FOR PEER REVIEW 15 of 22

interactions. The LID of WT-MRP4-ceefourin-1 is presented at the bottom.

Figure 19. Ceefourin-1 binding site in Y556C and the amplification of such a site. Rainbow colors in both the complete Y556C and amplification of the binding site represent the ceefourin-1 structure, while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar interactions. The LID of Y556C-ceefourin-1 is presented at the bottom. **Figure 19.** Ceefourin-1 binding site in Y556C and the amplification of such a site. Rainbow colors in both the complete Y556C and amplification of the binding site represent the ceefourin-1 structure, while yellow spheres represent the ATP. The yellow dotted lines indicate polar and non-polar interactions. The LID of Y556C-ceefourin-1 is presented at the bottom.

Umbrella sampling simulations were performed to determine the theoretical affinity, The cAMP and ceefourin-1 abnormal binding site in Y556C was one of the most remarkable findings in this work and showed the importance of studying the global effect of Y556C in different cell lines.

expressed as the free energy of binding (ΔG), of an event related to protein–ligand interactions along a reaction coordinate. Table 5 shows the ΔG obtained by Umbrella sampling through 10 ns of simulation on each complex. According to Umbrella sampling results, the Y556C-FA complex is the most stable with or without ATP bound to MRP4. Interestingly, ΔG for WT-FA and Y556C-FA complexes were lower in those MRP4 without ATP. In addition, ΔGs for G187W-cAMP and Y556C-cAMP complexes were higher without ATP, suggesting that the ATP bound is required for both FA and cAMP binding Umbrella sampling simulations were performed to determine the theoretical affinity, expressed as the free energy of binding (∆G), of an event related to protein–ligand interactions along a reaction coordinate. Table 5 shows the ∆G obtained by Umbrella sampling through 10 ns of simulation on each complex. According to Umbrella sampling results, the Y556C-FA complex is the most stable with or without ATP bound to MRP4. Interestingly, ∆G for WT-FA and Y556C-FA complexes were lower in those MRP4 without ATP. In addition, ∆Gs for G187W-cAMP and Y556C-cAMP complexes were higher without ATP, suggesting that the ATP bound is required for both FA and cAMP binding stabilization. ∆G for WT-MRP4-cAMP with or without ATP was the same. This represents an interesting finding because in WT-MRP4, the binding of ATP may promote the stabilization and conformational changes on MRP4 protein to promote an adequate interaction with FA to be properly transported out of the cell. Regarding the information provided by the LIDs, it is not possible to link the intermolecular interactions to ∆G because of the lack of a pattern in all the complexes. The intermolecular interactions do not define the increase or decrease of ∆G in this case. To confirm these observations, experimentally testing the FA efflux on cells expressing equal levels of WT-MRP4, G187W, and Y556C, as well as MRP4 inhibition by ceefourin-1 to measure the transport rate would afford further information of the differences among mutants and wild-type MRP4, experiments that are currently being carried out and will be reported in future articles.

The ∆G for the WT-MRP4-ATP1 complex was lower than that for FA and cAMP, which is reasonable considering its role in MRP4 functioning and needs to be bound to NBDs for longer time than ligands. In the case of ATP1/2 complexes with MRP4 variants, the most remarkable finding was that the Y556C-ATP1 complex was the least stable, theoretically demonstrated through the Umbrella sampling results, and it may be related to the mutation in the NBD1. The mechanism for the MRP4 function requires the

transmission of the molecular motion from the NBDs to the TMDs. At this point, the ATP-binding can be considered as the power stroke in which the chemical potential of the transported entity changed, and ATP hydrolysis leads to the formation of an extra negative charge, thus opening the closed nucleotide sandwich structure and the opening of the nucleotide sandwich structure facilitates Pi release and ADP dissociation, which in turn allows the TMDs and access gates to reset to the high-affinity orientation on the original side of the membrane to continue the transport cycle [23,29]. If the Y556C mutation promotes a decreased affinity of ATP1 for its binding site and, in turn, a blockade of ATP1 hydrolysis, the Y556C variant activity could be diminished or truncated, also considering that WT-MTP4 and its variants lack the ability to hydrolyze ATP2 in NBD2; thus, the lack of Y556C efficacy in the transport substrates. Additionally, ATP2 had the highest affinity for Y556C, and it would be interesting to test if high affinity could reestablish the ATP2 hydrolysis in NBD2.

**Table 5.** Binding free energies (∆Gs) obtained by umbrella sampling.

FA-ATP, cAMP-ATP, and ceefourin-1-ATP: the ∆G is referred to FA, cAMP or ceefourin-1 in complex with ATP. FA-ATP1, FA-ATP2, cAMP-ATP1, cAMP-ATP2, ceefourin-1-ATP1, and ceefourin-1-ATP2: ∆G calculation is focused on ATP1 or ATP2

The abnormal cAMP and ceefourin-1 binding in Y556C (Figures 16 and 19) is consistent with the ∆G for Y556C-cAMP and Y556-ceefourin-1, which had the lowest affinity with respect to the complexes with WT-MRP4 and G187W, and it indicates a low affinity at that binding site, suggesting two possibilities: (a) it is an initial cAMP or ceefourin-1 binding site to further binding at the inner cavity, or (b) possible deficiency in cAMP and ceefourin-1 transport. The latter is the most feasible due to cAMP interacts to G187W and WT-MRP4 in its normal binding site, but Y556C seems to be in the "low affinity outwardfacing orientation" as mentioned before. To predict the effect of nsSNPs on the efficacy of chemotherapeutics, it is important to determine in a further study, the relation between ∆G in silico and the transport rate in vitro of several substrates. This suggests that ∆Gs higher than that for the WT-cAMP-ATP complex may be related to a low rate of cAMP transport. The ∆Gs for ceefourin-1 in all the complexes with MRP4 and its variants had the best affinity compared to those complexes of cAMP, supporting the idea that ceefourin-1 may act as a competitive inhibitor at least with cAMP. However, the presence of ATP in WT-MRP4 and its variants promotes a better affinity to cAMP than ceefourin-1 in the complexes with G187W and Y556C. A competitive assay in vitro is required to determine if ceefourin-1 had the best affinity to WT-MRP4 and its variants. The Figures 20 and 21 represent the LIDs for FA and cAMP, respectively. In both figures, the intermolecular interactions of ligands with the absence of ATP are presented in the upper panel while the intermolecular interactions of ligands with the presence of ATP are shown in the bottom panel. The FA binding site is almost the same with or without ATP, but the number and the type of intermolecular interactions is different. It seems that ATP promotes a greater number of H-bonds between FA and all WT-MRP4 and its variants. Regarding the cAMP

binding site, the ATP binding did not modify it; thus, the effect of ATP seems to be related with changes in the intermolecular interactions and affinity that, according to Umbrella sampling results, the ATP decreases the ∆G in the G187W-cAMP-ATP and Y556C-cAMP-ATP complexes. Figure 22 presents the LIDs for ceefourin-1. The binding site is different in all the complexes; interestingly, such a binding site does not change in the presence of ATP. The main difference induced for ATP is the intermolecular interaction pattern, which led to great amount of π-π stacking and H-bonds. ATP complexes. Figure 22 presents the LIDs for ceefourin-1. The binding site is different in all the complexes; interestingly, such a binding site does not change in the presence of ATP. The main difference induced for ATP is the intermolecular interaction pattern, which led to great amount of π-π stacking and H-bonds. The MRP4 variants may predispose the population to a given disease regarding the site of the MRP4 that was affected by the mutation and the change in the affinity of a given substrate. The clinical implications of MRP4 variants have been studied over the past years and it is crucial to describe the relation of the MRP4 variants with diseases [30].

Molecules 2021, 26, x FOR PEER REVIEW 17 of 22

binding site to further binding at the inner cavity, or (b) possible deficiency in cAMP and ceefourin-1 transport. The latter is the most feasible due to cAMP interacts to G187W and WT-MRP4 in its normal binding site, but Y556C seems to be in the "low affinity outwardfacing orientation" as mentioned before. To predict the effect of nsSNPs on the efficacy of chemotherapeutics, it is important to determine in a further study, the relation between ΔG in silico and the transport rate in vitro of several substrates. This suggests that ΔGs higher than that for the WT-cAMP-ATP complex may be related to a low rate of cAMP transport. The ΔGs for ceefourin-1 in all the complexes with MRP4 and its variants had the best affinity compared to those complexes of cAMP, supporting the idea that ceefourin-1 may act as a competitive inhibitor at least with cAMP. However, the presence of ATP in WT-MRP4 and its variants promotes a better affinity to cAMP than ceefourin-1 in the complexes with G187W and Y556C. A competitive assay in vitro is required to determine if ceefourin-1 had the best affinity to WT-MRP4 and its variants. The Figures 20 and 21 represent the LIDs for FA and cAMP, respectively. In both figures, the intermolecular interactions of ligands with the absence of ATP are presented in the upper panel while the intermolecular interactions of ligands with the presence of ATP are shown in the bottom panel. The FA binding site is almost the same with or without ATP, but the number and the type of intermolecular interactions is different. It seems that ATP promotes a greater number of H-bonds between FA and all WT-MRP4 and its variants. Regarding the cAMP binding site, the ATP binding did not modify it; thus, the effect of ATP seems to be related with changes in the intermolecular interactions and affinity that, according to Umbrella sampling results, the ATP decreases the ΔG in the G187W-cAMP-ATP and Y556C-cAMP-

Figure 20. LIDs for FA with the different MRP4 structures with or without ATP. (A) WT-MRP4-FA complex. (B) G187W-FA complex. (C) Y556C-FA complex. (D) WT-MRP4-FA-ATP complex. (E) G187W-FA-ATP complex. (F) Y556C-FA-ATP complex. **Figure 20.** LIDs for FA with the different MRP4 structures with or without ATP. (**A**) WT-MRP4-FA complex. (**B**) G187W-FA complex. (**C**) Y556C-FA complex. (**D**) WT-MRP4-FA-ATP complex. (**E**) G187W-FA-ATP complex. (**F**) Y556C-FA-ATP complex. Molecules 2021, 26, x FOR PEER REVIEW 18 of 22

Commented [M61]: Please replace with a sharper

image.

Figure 21. LIDs for cAMP with the different MRP4 structures with or without ATP. (A) WT-MRP4 cAMP complex. (B) G187W-cAMP complex. (C) Y556C-cAMP complex. (D) WT-MRP4-cAMP-ATP complex. (E) G187W-cAMP-ATP complex. (F) Y556C-cAMP-ATP complex. **Figure 21.** LIDs for cAMP with the different MRP4 structures with or without ATP. (**A**) WT-MRP4 cAMP complex. (**B**) G187W-cAMP complex. (**C**) Y556C-cAMP complex. (**D**) WT-MRP4-cAMP-ATP complex. (**E**) G187W-cAMP-ATP complex. (**F**) Y556C-cAMP-ATP complex.

Figure 22. LIDs for ceefourin-1 with the different MRP4 structures with or without ATP. (A) WT-MRP4-ceefourin-1 complex. (B) G187W-ceefourin-1 complex. (C) Y556C-ceefourin-1 complex. (D) WT-MRP4-ceefourin-1-ATP complex. (E) G187W-ceefourin-1-ATP complex. (F) Y556C-ceefourin-1-

Structure prediction by protein threading for MRP4 was performed using its primary sequence (code: O15439) from the UniProt database [31,32]. MRP4 mutant models were

3.1. Protein Threading for WT-MRP4 and Its Variants G187W and Y556C

ATP complex.

3. Materials and Methods

Molecules 2021, 26, x FOR PEER REVIEW 18 of 22

Figure 22. LIDs for ceefourin-1 with the different MRP4 structures with or without ATP. (A) WT-MRP4-ceefourin-1 complex. (B) G187W-ceefourin-1 complex. (C) Y556C-ceefourin-1 complex. (D) WT-MRP4-ceefourin-1-ATP complex. (E) G187W-ceefourin-1-ATP complex. (F) Y556C-ceefourin-1- ATP complex. **Figure 22.** LIDs for ceefourin-1 with the different MRP4 structures with or without ATP. (**A**) WT-MRP4-ceefourin-1 complex. (**B**) G187W-ceefourin-1 complex. (**C**) Y556C-ceefourin-1 complex. (**D**) WT-MRP4-ceefourin-1-ATP complex. (**E**) G187W-ceefourin-1-ATP complex. (**F**) Y556C-ceefourin-1- ATP complex.

3. Materials and Methods 3.1. Protein Threading for WT-MRP4 and Its Variants G187W and Y556C Structure prediction by protein threading for MRP4 was performed using its primary sequence (code: O15439) from the UniProt database [31,32]. MRP4 mutant models were The MRP4 variants may predispose the population to a given disease regarding the site of the MRP4 that was affected by the mutation and the change in the affinity of a given substrate. The clinical implications of MRP4 variants have been studied over the past years and it is crucial to describe the relation of the MRP4 variants with diseases [30].

#### **3. Materials and Methods**

#### *3.1. Protein Threading for WT-MRP4 and Its Variants G187W and Y556C*

Structure prediction by protein threading for MRP4 was performed using its primary sequence (code: O15439) from the UniProt database [31,32]. MRP4 mutant models were made by the substitution of amino acids, G187W, or Y556C into the WT-MRP4 primary sequence. Each primary sequence was uploaded into the I-TASSER [17,18] server for the calculations of the models. The crystallography of MRP1 (PDBD code: 5UJ9) from *Bos taurus* was used as a template in all cases [33,34].

#### *3.2. Coarse-Grained Molecular Dynamics (CG-MD) Simulations*

All the systems (WT, G187W, and Y556C) for the simulations were built in the Martini Maker module [35] from CHARMM-GUI [36,37] using the Martini2.2p force field and adding a phosphatidylcholine (POPC) lipid bilayer membrane in an isothermal-isobaric ensemble (NPT) at 310.15 K. The simulations of 1 µs CG-MD were carried out in the Gromacs 2018.7 program [38,39], after of the minimizations and equilibrium protocols suggested by CHARM-GUI server. The module "cluster" of Gromacs 2018.7 was used to find the relevant conformations of the simulation using the "gromos" algorithm and the backbone for alignment. The most representative structure of the largest conformation cluster of the three simulations was converted into all atoms in the Martini Maker/All-Atom converter from CHARMM-GUI for the following calculations.

All the CG-MD simulations were performed in the ADA cluster of the National Laboratory of Advanced Scientific Visualization at Campus Juriquilla of the National Autonomous University of Mexico (LAVIS-UNAM).

#### *3.3. Molecular Docking on WT-MRP4, G187W and Y556C*

The 3D structures of the selected ligand groups, substrate drugs, endogenous substrates, and inhibitors were obtained from the PubChem public database [40]. Molecular docking was performed with AutoDock 4.2.6 optimized for graphical-processing units using a total of 50 runs and 25,000,000 evaluations; a grid box of 22.5 × 22.5 × 22.5 Å<sup>3</sup> centering on the relevant amino acids reported by Ravna 2008 and 2009, El-Sheink 2008, and Chen 2018 [9,41–43]; in a Lamarckian genetic algorithm and Solis-Wets local search [44].

#### *3.4. All-Atom Molecular Dynamics (AA-MD) Simulations*

For molecular dynamic studies with FA as a ligand, the protein–ligand complexes WT-MRP4-FA, G187W-FA, or Y556C-FA, obtained from the molecular docking, were used, and AA-MD simulations were performed in Desmond 3.6 as an application of the Maestro software [45,46] as graphical interface. The AA-MDS systems were built using the System Setup module with an OPLS force-field (Optimized Potentials for Liquid Simulations), adding a POPC lipid membrane and simple point charge (SPC) water model in an NPT assembly at 310.15 K. Once the system was built, the standard relaxation protocol for system relaxation with increasing temperatures and decreasing restraints was used and an MD production simulation of 25 ns (100 frames) was performed in the Molecular Dynamics module. Clustering was performed in the Desmond Trajectory Clustering module to obtain the most representative conformation of the largest cluster (cluster1, C1). Trajectory analyses were performed in the Simulation Event Analysis module in Maestro.

#### *3.5. Umbrella Sampling*

Using the C1 of the 25 ns AA-MD trajectory, the system was built under the previously mentioned conditions. Once the system was obtained and relaxed using the AA-MD relaxation protocol, a 10 ns (100 frames) MD simulation was performed in the Metadynamics module of Desmond using the protein and ligand center-of-mass distance as the collective variable, with 0.3 kcal/mol height and 0.1 kcal/mol width as the Gaussian parameters for the umbrella protocol, on an NPT ensemble at 310.15 K and 1.01325 bar. Finally, the analysis was performed in the Metadynamics Analysis module of Desmond.

#### *3.6. Manipulation of the Complexes and Figures*

All the protein figures and alignments presented in this work were made in PyMOL software (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger) [47]. The manipulation of the WT-MRP4, G187W, and Y556C structures was performed in Maestro (Schrödinger Release 2019-3: Maestro, Schrödinger, LLC, New York, NY, USA) [48].

#### **4. Conclusions**

To obtain the 3D structure of MRP4 and its variants, which are not resolved by NMR or X-ray, protein threading was performed in this study and relaxing the structures was performed by CG-MDS to carry out the molecular docking, while MDS and umbrella sampling studies were performed to yield relevant information regarding the residues involved in the binding of the studied molecules groups and changes in the ∆G of FA and cAMP in the presence or the absence of ATP, which also allowed us to observe the relevance of the mutations in the binding and movement of MRP4 and its variants. According to our results, the nsSNPs G187W and Y556C led to changes in the ligand binding site, DS, and binding energy (∆G). In addition, the ATP binding to MRP4 significantly modifies the intermolecular interactions (at least, for FA and cAMP) and the binding energy compared to the complexes where ATP was not bound to MRP4. The effect of the abnormal binding site of cAMP in Y556C is consistent with the highly selective MRP4 inhibitor ceefourin-1, which makes it interesting to study such mutations in vitro. The affinity of ceefourin-1 for WT-MRP4 and its variants is higher than the affinity of cAMP. Cofactors such as ATP and Mg2+ should be included in the in silico analyses related to MRP4. It is well known that nonsynonymous mutations usually affect the protein function or activity and its conformation, but this is the first report that suggests that most endogen substrates change their affinity and binding site in G187W and Y556C, which could modify the cell metabolism. We will report in further works the measure of substrate efflux and the relation between G187W or Y556C expression, location, and cell viability to determine the overall effect of the nsSNPs G187W and Y556C in vitro.

**Supplementary Materials:** The following are available online: Table S1: RMSD values obtained for the alignment of C1 cluster from AA-MDS for different sites of WT-MRP4, G187W and Y556C, Figure S2: Binding site and intermolecular interaction of cholic acid in WT-MRP4 (A), G187W (B) and Y556C (C). D Nomenclature of LIDs, Figure S3: Binding site and intermolecular interaction of taurocholic acid in WT-MRP4 (A), G187W (B) and Y556C (C), Figure S4: Binding site and intermolecular interactions of cefazoline in WT-MRP4 (A), G187W (B) and Y556C (C), Figure S5: Binding site and intermolecular interaction of ceefourin-1 in WT-MRP4 (A), G187W (B) and Y556C (C), Figure S6: WT-MRP4-FA complex at T0 in AA-MDS, Figure S7: WT-MRP4-FA complex at 5 ns in AA-MDS, Figure S8: WT-MRP4-FA complex at 10 ns in AA-MDS, Figure S9: WT-MRP4-FA complex at 15 ns in AA-MDS, Figure S10: WT-MRP4-FA complex at 20 ns in AA-MDS, Figure S11: WT-MRP4-FA complex at 25 ns in AA-MDS, Figure S12: G187W-FA complex at T0 in AA-MDS, Figure S13: G187W-FA complex at 5 ns in AA-MDS, Figure S14: G187W-FA complex at 10 ns in AA-MDS, Figure S15: G187W-FA complex at 15 ns in AA-MDS, Figure S16: G187W-FA complex at 20 ns in AA-MDS, Figure S17: G187W-FA complex at 25 ns in AA-MDS, Figure S18: Y556C-FA complex at T0 in AA-MDS, Figure S19: Y556C-FA complex at 5 ns in AA-MDS, Figure S20: Y556C-FA complex at 10 ns in AA-MDS, Figure S21: Y556C-FA complex at 15 ns in AA-MDS, Figure S22: Y556C-FA complex at 20 ns in AA-MDS, Figure S23: Y556C-FA complex at 25 ns in AA-MDS.

**Author Contributions:** Conceptualization: E.B., L.B., and G.G.-A.; calculations: E.B., G.A.-D., and A.R.-M.; data analysis: E.B., G.A.-D., and A.R.-M.; writing—original draft preparation: E.B., G.A.-D.; writing—review and editing: E.B., L.B., G.G.-A., G.A.-D., L.B., and A.R.-M.; project administration and funding acquisition: L.B. and G.G.-A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the University of Queretaro (UAQ) through the FOFI-UAQ research grant FCQ2018-39 and the National Council of Research and Technology (CONACYT) doctoral scholarship no. 335389 for Edgardo Becerra and no. 335496 for Giovanny Aguilera-Durán.

**Institutional Review Board Statement:** Not Applicable.

**Informed Consent Statement:** Not Applicable.

**Acknowledgments:** The authors gratefully acknowledge the technical support from Luis Aguilar, Alejandro de León, Carlos Flores, and Jair García of the National Laboratory of Advanced Scientific Visualization (LAVIS-UNAM).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Biocomputational Prediction Approach Targeting FimH by Natural SGLT2 Inhibitors: A Possible Way to Overcome the Uropathogenic Effect of SGLT2 Inhibitor Drugs**

**Mutaib M. Mashraqi 1,†, Navaneet Chaturvedi 2,3,† , Qamre Alam <sup>4</sup> , Saleh Alshamrani <sup>1</sup> , Mosa M. Bahnass 1,5 , Khurshid Ahmad <sup>6</sup> , Amany I. Alqosaibi <sup>7</sup> , Mashael M. Alnamshan <sup>7</sup> , Syed Sayeed Ahmad <sup>6</sup> , Mirza Masroor Ali Beg <sup>6</sup> , Abha Mishra <sup>2</sup> , Sibhghatulla Shaikh 6,\* and Syed Mohd Danish Rizvi 8,\***

**Abstract:** The Food and Drug Administration (FDA) approved a new class of anti-diabetic medication (a sodium–glucose co-transporter 2 (SGLT2) inhibitor) in 2013. However, SGLT2 inhibitor drugs are under evaluation due to their associative side effects, such as urinary tract and genital infection, urinary discomfort, diabetic ketosis, and kidney problems. Even clinicians have difficulty in recommending it to diabetic patients due to the increased probability of urinary tract infection. In our study, we selected natural SGLT2 inhibitors, namely acerogenin B, formononetin, (−)-kurarinone, (+)-pteryxin, and quinidine, to explore their potential against an emerging uropathogenic bacterial therapeutic target, i.e., FimH. FimH plays a critical role in the colonization of uropathogenic bacteria on the urinary tract surface. Thus, FimH antagonists show promising effects against uropathogenic bacterial strains via their targeting of FimH's adherence mechanism with less chance of resistance. The molecular docking results showed that, among natural SGLT2 inhibitors, formononetin, (+)-pteryxin, and quinidine have a strong interaction with FimH proteins, with binding energy (∆G) and inhibition constant (ki) values of −5.65 kcal/mol and 71.95 µM, −5.50 kcal/mol and 92.97 µM, and −5.70 kcal/mol and 66.40 µM, respectively. These interactions were better than those of the positive control heptyl α-Dmannopyranoside and far better than those of the SGLT2 inhibitor drug canagliflozin. Furthermore, a 50 ns molecular dynamics simulation was conducted to optimize the interaction, and the resulting complexes were found to be stable. Physicochemical property assessments predicted little toxicity and good drug-likeness properties for these three compounds. Therefore, formononetin, (+)-pteryxin, and quinidine can be proposed as promising SGLT2 inhibitors drugs, with add-on FimH inhibition potential that might reduce the probability of uropathogenic side effects.

**Citation:** Mashraqi, M.M.; Chaturvedi, N.; Alam, Q.; Alshamrani, S.; Bahnass, M.M.; Ahmad, K.; Alqosaibi, A.I.; Alnamshan, M.M.; Ahmad, S.S.; Beg, M.M.A.; et al. Biocomputational Prediction Approach Targeting FimH by Natural SGLT2 Inhibitors: A Possible Way to Overcome the Uropathogenic Effect of SGLT2 Inhibitor Drugs. *Molecules* **2021**, *26*, 582. https://doi.org/10.3390/ molecules26030582

Academic Editors: Marco Tutone and Anna Maria Almerico Received: 18 December 2020 Accepted: 19 January 2021 Published: 22 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** sodium–glucose co-transporters 2; FimH; uropathogenic bacteria; urinary tract infections; diabetes

#### **1. Introduction**

Globally, diabetes mellitus is one of the most prevalent metabolic diseases and is estimated to increase to 552 million cases by 2030 [1]. Diabetes has been considered to enhance vulnerability to infectious diseases and has also been linked with an enhanced risk of death from infectious illness in some [2,3], but not all [4], investigations.

Sodium–glucose cotransporter-2 (SGLT2) inhibitors are a novel group of drugs used to treat patients with type 2 diabetes mellitus (T2DM). These drugs exert their effects by preventing glucose reabsorption at the proximal renal tubule and enhancing the excretion of urinary glucose [5]. Owing to the elevated urine glucose levels, SGLT2 inhibitors enhance the risk of urinary tract infections (UTIs) [6]. In addition, pharmacologicallyinduced urine glucose levels with SGLT2 inhibitors in diabetic patients might further sustain bacterial growth [7]. By themselves, SGLT2 inhibitors can possibly enhance the risk of UTIs and susceptibility to genital infections when used to manage patients [8]. SGLT2 inhibitors might cause serious UTIs, as the FDA warned in December 2015 [9]. Empagliflozin and canagliflozin are the preferred drugs suggested for T2DM patients with established cardiovascular disease [10,11]. The United States-based public safety advisory reported 19 cases of fatal renal or blood infection from March 2013 to October 2014. The origin of these cases was traced to a UTI induced by SGLT2 inhibitor intake [9].

Bacterial pili are proteinaceous projections extending from the bacterial cell surface and are used for attachment and cell motility [12]. Gram-negative bacteria use Type 1 pili to adhere to the host tissue and, therefore, Type 1 pili have been established as an important virulence factor in UTIs. Type 1 pili are made up of repeated subunits of FimA protein. These subunits are arranged to form a helical wound cylinder that composes a thick pilus rod. The distal flexible tip of the pilus rod is comprised of a single copy of proteins—FimF and FimG—while the tip adhesin bears FimH. The distinct binding of FimH (terminal adhesin) to mannosylated host glycoproteins mediates the adhesion of bacterial pathogens to the host tissue. UTI pathogenesis is caused by the FimH and, therefore, is a promising curative target [13].

Acerogenin B is a cyclic diarylheptanoid obtained from the bark of *Acer nikoense*, and it has been found to inhibit both SGLT1 and SGLT2 [14]. Kurarinone and formononetin are flavonoids isolated from the dried root of *Sophora flavescens.* Kurarinone has demonstrated inhibitory activity against both SGLT1 and SGLT2; however, formononetin was reported to inhibit only SGLT2 and not SGLT1 [15]. Quinidine is isolated from the bark of the cinchona tree and (+)-pteryxin is extracted from the plant *Peucedanum* spp. Recently, both of these natural compounds have been found to inhibit both SGLT1 and SGLT2 [16]. Heptyl α-D-mannopyranoside is a FimH antagonist [17] and, in this study, it was used as a reference compound.

In the present study, we explored natural SGLT2 inhibitors (acerogenin B, kurarinone, formononetin, quinidine, and (+)-pteryxin) that might have less severe uropathogenic side effects than does the approved SGLT2 inhibitor canagliflozin, also known as gliflozins. We speculated that formononetin, (+)-pteryxin, and quinidine would be promising SGLT2 inhibitors with less severe uropathogenic side effects.

#### **2. Methodology**

#### *2.1. SGLT2 Inhibitors and Target Protein Structure Retrieval*

The 3-dimensional structure of FimH was taken from the protein data bank (PDB ID: 4AV5), while the 3-dimensional structure of SGLT2 was made by employing the SWISS-MODEL Workspace after retrieving its amino acid sequence from Uniprot (P31639). The PDB structure 2XQ2 was used as a template and the generated model was validated

using various in silico tools, viz., SAVES v6.0 and VADAR (Volume, Area, Dihedral Angle Reporter). The PDB structure of the compounds canagliflozin (CID: 24812758), acerogenin B (CID: 10913542), formononetin (CID: 5280378), (−)-kurarinone (CID: 10812923), (+) pteryxin (CID: 511787), quinidine (CID: 441074), and the FimH antagonist heptyl α-Dmannopyranoside (CID: 11300413) were retrieved from the PubChem database.

#### *2.2. Physicochemical Properties and Toxicity Potential Prediction*

The physicochemical properties and toxicity potential of SGLT2 inhibitors and FimH antagonists were calculated by applying the Orisis Datawarrior property explorer tool (http://www.openmolecules.org/datawarrior/download.html). At first, molecular weight, the number of rotatable bonds, the number of hydrogen bond acceptors and donors, cLogP value, topological polar surface area (TPSA), and the Lipinski's rule violation [18] were estimated to check physicochemical parameters. Later, the method outlined by Zhao et al. [19] was applied to calculate percentage of absorption; here, the following formula was used: absorption % = 109 − (0.345 × TPSA). Toxicity was also predicted for SGLT2 inhibitors and FimH antagonists by using toxicity features such as irritability, reproductive effects, tumorigenicity, and mutagenicity. Orisis Datawarrior property explorer tool toxicity predictions are based on comparative analysis of our tested compounds with the preestimated set of known structural molecules.

#### *2.3. Molecular Docking*

The ligands were docked to the SGLT2 and FimH proteins using "AutoDock 4.2" by following the protocol described by Rizvi et al. [20]. To minimize the energy usage of the ligand molecules, a Merck molecular force field (MMFF94) was employed. The ligand atoms were added with Gasteiger partial charges. Docking calculations were done on the target proteins. Essential hydrogen atoms, Kollman united atom type charges and solvation parameters were added by using AutoDock tools. Consequently, the binding pocket was added with conserved water molecules to mimic the *in vivo* environment. An auto grid program was used to generate the affinity (grid) maps sized at 60 × 60 × 60◦ , the aim of which was to target the grid coordinates in the catalytic site of the target protein (SGLT2 and FimH). The x, y, and z coordinate values for the FimH protein targeting the catalytic site were taken as 3.305, −16.68, and −13.57, respectively. Docking simulations were done using the Lamarckian genetic algorithm and the Solis and Wets local search method. The initial position, orientation, and torsions of the ligands were set randomly. Each docking experiment was derived from 100 different runs that were set to terminate after a maximum of 2,500,000 energy evaluations. The population size was set to 150. The Discovery Studio Visualizer 2.5 (Accelrys, San Diego, CA, USA) was used to generate the final figures.

#### *2.4. LIGPLOT+ ANALYSIS*

After the completion of all docking experiments, the best combination of ligand–FimH and ligand–SGLT2 was selected and analyzed using LIGPLOT+ Version v.2.1 (EMBL-EBI, Cambridgeshire, UK). The hydrogen and hydrophobic interactions between the important amino acid residues of FimH/SGLT2 with each ligand were analyzed by LIGPLOT. For each interaction, the 3-D structures generated were converted into 2-D figures by applying the LIGPLOT algorithm.

#### *2.5. Molecular Dynamics Simulation*

System building: GROMACS 4.6.7 packages [21,22] were used to prepare the system and perform molecular dynamics (MD) simulations using the gromos53a6 force field [23]. The protein solute was solvated by explicit SPC216 water [24] in a dodecahedral box with a margin of 10 Å between the solute and the box walls. Systems were brought to neutrality by the addition of sodium counter ions.

Simulation detail: A 10 Å cut-off distance was taken under the particle-mesh Ewald method [25] to calculate long-range electrostatic interactions, and a 10 Å cut-off distance was also considered to evaluate van der Waals interactions. The LINCS algorithm of fourth-

order expansion was used to constrain bond lengths [26]. The steepest descent algorithm was applied to optimize the removal of steric clashes between atoms for 10,000 steps. The system was equilibrated for 1 ns with position restraints on all heavy atoms. Berendsen weak coupling schemes were used to maintain the system at 300 K and 1 atom, using separate baths for the system. Initial velocities were generated randomly using a Maxwell– Boltzmann distribution corresponding to 300 K. Finally, the production run was performed for 20 ns. Furthermore, xmgrace (http://plasma-gate.weizmann.ac.il) was used for preparing graphs. Ligand topology preparation was implemented by using the PRODRG server, using the option specifying no chirality, full charge, and no energy minimization [27].

Trajectory analysis: The g\_rms tool of the GROMACS package was used to calculate the root-mean-square deviation (RMSD) of each trajectory. The covariance matrix and eigenvectors of the trajectories were calculated through g\_covar and g\_anaeig programs.

#### **3. Results and Discussion**

#### *3.1. Physicochemical Properties and Molecular Docking*

Many varieties of plant-based natural compounds have been reported, which have significant anti-diabetic effects. In streptozotocin-stimulated diabetic mice, bakuchiol (a polyphenol compound) decreases glucose levels and enhances serum insulin levels [28]. Caffeine (an alkaloid) lowers glucose, induces insulin secretion in vitro, and improves glucose absorption in skeletal muscle [29]. Vanillic acid isolated from *Fagara tessmannii Engl.* (Family: Rutaceae) exhibits α-glucosidase inhibitory actions in vitro [30]. Christinin-A is a triterpenoidal saponin glycoside that exhibits an anti-hyperglycemic effect in both type 1 and type 2 diabetic rats [31].

Molecular docking has an important role in drug discovery, assisting in digging out the active or lead compounds from a library of natural compounds [32]. It is one of the most widely used virtual screening tools, particularly when the three-dimensional structure of the target protein is available. Docking enables the prediction of both ligand–target binding affinity and the structure of the protein–ligand complex, which are useful for optimizing the lead [33,34].

Prior to the molecular docking study, we checked the physicochemical properties and toxicity potential of SGLT2 inhibitors and FimH antagonists (Tables 1 and 2). The selected SGLT2 inhibitors were known drug canagliflozin and natural SGLT2 inhibitors, namely, acerogenin B, formononetin, (−)-kurarinone, (+)-pteryxin, and quinidine, while heptyl α-Dmannopyranoside was selected as an FimH antagonist. During physicochemical property assessment, we found that out of all the compounds tested, (−)-kurarinone showed one violation of the Lipinski rule [18], i.e., a cLogP value (Logarithm of compound partition coefficient between n-octanol) higher than 5 (Table 1). On the other hand, all the tested compounds showed no toxicity except (+)-pteryxin. (+)-pteryxin was predicted to have a high irritant effect with no mutagenic, tumorigenic, or reproductive toxicity (Table 2).

**Table 1.** Physicochemical properties of natural sodium–glucose co-transporter 2 (SGLT2) inhibitors and control compounds.

\* Control drugs/compounds; \*\* Percentage of Absorption (% of Absorption) was calculated by the following equation: % of Absorption = 109 − (0.345 × Topological Polar Surface Area); \*\*\* Logarithm of the compound partition coefficient between n-octanol and water.

**Table 2.** Toxicity potential of natural SGLT2 inhibitors and control compounds.

\* Control drugs/compounds.

The predicted model of SGLT2 revealed that 81.25 percent of the residues had an average 3D-1D score of ≥0.2, resulting in a "pass" in SAVES v6.0.0 by the VERIFY3D tool. In addition, the Ramachandran plot (showing 93% of the residues in the allowed region), fractional accessible surface area, stereo/packing quality index, fractional residue volume, and 3-D profile quality index (produced by the VADAR 1.8 server) showed that the predicted 3-D model was within an appropriate range (Figure 1).

**Figure 1.** *Cont*.

**Figure 1.** *Cont*.

**Figure 1.** Validation of the predicted 3-D structure of SGLT2.

The inhibition of SGLT2 has been considered a novel pharmacotherapy for T2DM treatment [35]. Accordingly, molecular docking studies have revealed that the natural SGLT2 inhibitors formononetin, (+)-pteryxin, and quinidine were efficiently bound with SGLT2. Formononetin was found to interact with the F98, E99, A102, L149, K152, T153, V286, S287, Y290, W291, I456, Q457, and S460 binding pocket residues of SGLT2 (Figure 2a); while the S74, G79, H80, K154, I155, V157, D158, S161, S393, I397, and I456 residues of SGLT2 were observed to bind with (+)-pteryxin (Figure 2b). Furthermore, quinidine was found to interact with the S74, N75, H80, L84, F98, E99, V286, S287, Y290, W291, F453, I456, and Q457 residues of SGLT2 (Figure 2c). Consistent with this, the amino acid residues H80, F98, T153, K154, V157, D158, V290, I397, F453, I456, Q457, and S460 were shown to be important for the inhibition of SGLT2 [36]. Amino acid residues H80, V286, Y290, W291, and I456 are the main hydrophobic residues of SGLT2, interacting with the dock inhibitors formononetin, (+)-pteryxin and quinidine (Figure 3a–c). This concurs with a previous report wherein the amino acid residues H80, Y290, and I456 of SGLT2 have also been reported have a hydrophobic interaction with the inhibitor [36,37]. The binding energy (BE) for the catalytic domain interactions of formononetin–SGLT2, (+)-pteryxin–SGLT2, and quinidine–SGLT2 were found to be −7.63 kcal/mol, −9.01 kcal/mol, and −8.77 kcal/mol,

respectively, while their inhibition constants (Ki) were 2.57 µM, 0.245 µM, and 0.371 µM, respectively (Table 3).

**Figure 2.** The amino acid residue of SGLT2 interacting with formononetin (**a**), (+)-pteryxin (**b**), quinidine (**c**), and canagliflozin (**d**). The ligands (formononetin, (+)-pteryxin, quinidine, and canagliflozin) are represented as green stick forms and hydrogen bonds are indicated as green dashed lines.

**Table 3.** The docking results of the molecular interactions of SGLT2 and FimH with natural SGLT2 inhibitors and control compounds.

\* Control compound for SGLT2; \*\* Control compound for FimH.

Heptyl α-D-mannopyranoside \*\* - - −4.46 109.49

(+)-pteryxin −9.01 0.248 −5.50 92.97 Quinidine −8.77 0.371 −5.70 66.40

> The FDA approved canagliflozin, an SGLT2 inhibitor, for use in T2DM treatment in 2013 [38]. In the present study, canagliflozin was used as a positive control against SGLT2. Canagliflozin was observed to bound with G79, H80, Y150, K154, D158, W289, Y290, D294, S393, I397, S460, and I456 residues of SGLT2 (Figures 2d and 3d). Interestingly, these amino acids of SGLT2 have also been found to interact with SGLT2 inhibitors (formononetin,

(+)-pteryxin, and quinidine). I456 was one of the most reactive common hydrophobic residues of SGLT2, interacting with formononetin, (+)-pteryxin, and quinidine, as well as canagliflozin (Figure 3a–d). The Y290 residue of SGLT2 was found to make H-bonds with quinidine, while the same residue was observed to hydrophobically interact with formononetin and canagliflozin (Figure 3a,c,d). In addition, G79, H80, and I397 were the common hydrophobic interacting residues of SGLT2 with (+)-pteryxin and canagliflozin (Figure 3b,d).

**Figure 3.** Ligplot analysis of formononetin (**a**), (+)-pteryxin (**b**), quinidine (**c**), and canagliflozin (**d**) in terms of their interaction with SGLT2. The amino acid residues forming hydrophobic interactions are shown as red arcs while the hydrogen bonds are shown as green dashed lines with indicated bond lengths.

The function of human SGLT1 protein is dramatically affected by amino acid at position 457. It has been shown that residue 457 (i.e., Q457) in human SGLT1 directly interacts with sugar for its reabsorption [39,40]. The amino acid sequences of SGLT1 and SGLT2 revealed that both of these proteins have glutamine at the residue 457 position, and glucose–galactose malabsorption occurs due to mutation in the glutamine residue (Q457) [41]. Consistent with this, in the present study, formononetin and quinidine were observed to interact with the Q457 residue that is suggested to impair SGLT2 function.

Since the discovery of the first natural SGLT2 inhibitor (i.e., phlorizin), several other synthetic glucoside analogs have been developed [42]. Tofogliflozin is a selective SGLT2 inhibitor that enhances urinary glucose excretion in a dose-dependent manner [43]. Luseogliflozin is an orally active second-generation SGLT2 inhibitor that has protective effects on pancreatic beta-cell mass and function [44]. Furthermore, efforts have been made to find new active natural ingredients from *S. Flavescents* [15], alkaloids from *A. macrophylla* [45], and *Schisandrae Chinensis Fructus* for the development of SGLT2 inhibition [46]. A 4-O-methyl derivative of sergliflozin-A has been reported to exhibit SGLT2 inhibition activity [47].

SGLT2 inhibitors are some of the most promising anti-diabetic agents introduced into clinical practice over the last decade. The therapeutic benefits of SGLT2 inhibitors include weight loss, a reduction in blood pressure, and an enhancement in high-density lipoprotein level [48,49]. However, the SGLT2 inhibitors dapagliflozin, canagliflozin, and empagliflozin have been found to cause UTIs and genital infections in diabetic patients [50–52]. Additionally, the FDA has issued warnings about the possible UTI-inducing side effects of SGLT2 inhibitors in December 2015 [9].

UTIs pose a severe medical issue worldwide, with more than 85% of UTIs caused by uropathogenic *Escherichia coli* [53]. FimH is a bacterial adhesin that facilitates the colonization of uropathogenic *E. coli* on the cell surface of the human and murine bladder and leads to the formation of biofilm [54]. Therefore, this adhesin has been considered a virulence factor and a promising therapeutic target for the treatment of UTIs [55]. Interestingly, the docking results indicate that the SGLT2 inhibitors formononetin, (+)-pteryxin, and quinidine show strong binding with FimH. FimH was found to interact with formononetin through 11 amino acid residues, namely F1, N46, D47, Y48, I52, D54, Q133, N135, Y137, N138, and D140 (Figure 4a), while (+)-pteryxin was found to interact with 6 amino acid residues, namely Y48, T51, I52, Y137, N138, and D140 (Figure 4b). Similarly, 10 amino acid residues, namely the F1, I13, Y48, I52, D54, Q133, N135, Y137, N138, and D140 residues of FimH, were found to interact with quinidine (Figure 4c). The BE for formononetin–FimH, (+)-pteryxin–FimH, and quinidine–FimH interactions were found to be −5.65 kcal/mol, −5.50 kcal/mol, and −5.70 kcal/mol, respectively. The corresponding estimated inhibition constants (Ki) for the above-mentioned complexes were determined to be 71.95 µM, 92.97 µM, and 66.40 µM, respectively (Table 3).

**Figure 4.** The amino acid residue of FimH's interaction with formononetin (**a**), (+)-pteryxin (**b**), quinidine (**c**), and heptyl α-D-mannopyranoside (**d**). The ligands (formononetin, (+)-pteryxin, quinidine, and heptyl α-D-mannopyranoside) are represented as green stick forms and hydrogen bonds are represented as green dashed lines.

The FimH protein has two domains, the C-terminal pilin domain and the N-terminal lectin domain. The FimH lectin domain has a mannose-binding pocket (N46, D47, D54, Q133, N135, N138, and D140), in which sugar contributes to the formation of various hydrogen bonds with FimH. Hydrophobic regions are present in the surrounding region of the mannose-binding pocket and consist of hydrophobic support (F1, I13, and F142), the tyrosine gate (Y48, I52, and Y137), and the residue T51 [56]. Consistent with this, in the present study, the D47, D54, and N138 of the FimH protein were involved in hydrogen binding with formononetin (Figure 5a), while Y48, I52, Y137, and D140 were the common hydrophobic amino acid residues interacting with formononetin, (+)-pteryxin and quinidine (Figure 5a–c).

**Figure 5.** Ligplot analysis of the interactions of formononetin (**a**), (+)-pteryxin (**b**), quinidine (**c**), and heptyl α-Dmannopyranoside (**d**) with SGLT2. The amino acid residues forming hydrophobic interactions are shown as red arcs while the hydrogen bonds are shown as green dashed lines with indicated bond lengths.

In the present study, heptyl α-D-mannopyranoside (a FimH antagonist) was used as a positive control against the FimH protein. Molecular docking analysis revealed that amino acid residues, namely the F1, I13, Y48, T51, Q133, N135, N138, D140, and I52 residues of the FimH protein, have a vital role in binding with heptyl α-D-mannopyranoside (Figures 4d and 5d). The main FimH binding pocket amino acid residues are F1, N46, D47, D54, E133, N135, D140, and F142, and mutation in these individual main residues leads to them affecting FimH function and reducing its virulence [57,58]. Interestingly, we found that these amino acid residues of FimH were also determined to interact with natural SGLT2 inhibitors (formononetin, (+)-pteryxin, and quinidine).

O- and C-linked alpha-D-mannosides with hydrophobic and aryl substituents are effective FimH antagonists. The substitution of biphenyl derivatives may provide additional advantages for the FimH antagonist molecule. The addition of 1,10-biphenyl pharmacophore and various aglycone atoms enhanced the alpha-D-mannose derivatives' suitability as FimH inhibitors [59,60]. Further, glycomimetics based on mannose scaffolding

has also been synthesized and widely studied for their aptitude as FimH antagonists [61,62]. Besides synthetic compounds, natural substances, like cranberry and its derivatives (such as myricetin, proanthocyanidin (PAC)-standardized cranberry extract, and polyphenol metabolites extracted from PAC), have anti-adhesive effects on uropathogenic *E. coli* [63].

In docking experiments, it is usually crucial to look for a ligand that can bind efficiently with the receptor, using Gibbs free energy as a parameter of better binding [64,65]. The strength of an interaction between a ligand and a receptor is measured in terms of the free energy of binding. The lowest BE is the output of the efficient binding of the drug/ligand to the active site of the receptor [66]. A higher (negative) BE is a sign of efficient interaction between the ligand and the receptor [67]. Accordingly, in the present study, formononetin, (+)-pteryxin, and quinidine exhibited strong interaction with the FimH protein, with a high BE relative to that of the positive control heptyl α-D-mannopyranoside, suggesting that these compounds could be promising SGLT2 inhibitors with less severe uropathogenic side effects.

Diabetic ketoacidosis and its associated events have been reported at a low frequency in T2DM patients treated with canagliflozin [68]. Interestingly, in the present study, the BE of canagliflozin with FimH relative to the other SGLT2 inhibitor indicates that canagliflozin has less potency in inhibiting the FimH protein, thereby having a high susceptibility to causing diabetic ketosis and UTIs.

#### *3.2. Root Mean Square Deviation*

Root mean square deviation (RMSD) is the most significant dynamic to explore in terms of conformational changes by means of stability in the structure and dynamic behavior of the receptor [69]. The complexes of formononetin, quinidine, and (+)-pteryxin with FimH obtained in molecular docking that interacted the best were further probed by MD simulation. The binding of ligands to their receptor protein can result in large conformational deviations in the resulting complex, specifically at the binding site. In this study, values of RMSD were increased at the beginning, with respect to the native structure of the RMSD. Figure 6A reports the RMSD of the backbone atoms of FimH protein as they interacted with each ligand molecule, where formononetin–FimH, quinidine–FimH and (+)-pteryxin–FimH interactions are represented with a black, red, and green color code, respectively. RMSD changes in all systems were initiated at the same point (~0.125 nm); however, in the formononetin–FimH complex, the conformational changes took place throughout the stabilization process. The change in the RMSD values in the formononetin– FimH complex indicated a small fluctuation until ~30 ns, and, afterwards, RMSD was found to be stable, remaining at around 0.35 nm in value. Due to the initial ups and downs in the RMSD values of the formononetin–FimH complex, the topologies of the structures were observed to significantly change during simulation. The observed initial RMSD change predicted the large conformational changes in regions near the binding pocket. Although the simulation lasted for 50 ns, the discrepancy from the initial structure within the first 30 ns was sufficient to point out the protein structures that were substantially denatured at binding sites and adjacent areas. While the RMSD values of the quinidine– FimH complex showed little turbulence until ~20 ns, the values reached a plateau later on. The value observed was greater than 0.3 nm initially; however, after 35 ns it went steady at 0.37 nm. In general, overall RMSD values fluctuated within 0.5 nm ± 0.05 nm throughout the simulation.

Like the above complexes, RMSD values showed frequent conformational changes and binding residue incorporation during the simulation analysis of the pteryxin–FimH complex. Conformational changes were observed at a functional domain near the binding site of the complex, and a considerable movement of the pteryxin ligand was also predicted at that site. Moreover, the RMSDs of all three ligand candidates were calculated to evaluate the behavior of the ligands within the binding site of the protein. Figure 6B shows the plot representation of RMSD for each ligand. Ligands (formononetin and quinidine) have generally shown high stability within the close contact residues of protein, whereas, in the

case of pteryxin, this ligand demonstrated several fluctuations within the binding site. Due to the movement of the pteryxin ligand, an initial sudden elevation (~0.13 nm) followed by fluctuations until 0.15 ns was observed during the simulation. Although the RMSDs maintained a plateau later on, an average value of 0.21 nm was calculated throughout the simulation. Moreover, the difference in RMSDs due to the movement of pteryxin from the initial to the final frame at the binding site was observed as 0.15 nm, which revealed that pteryxin shows considerably weaker binding with interacting amino acid residues during the simulation. On the other hand, it was observed that formononetin and quinidine stably interacted during simulation, and the pattern of contact between the binding residues and both ligands (formononetin and quinidine) predicted regulation of the activity of the FimH protein.

**Figure 6.** Representation of an RMSD plot of each complex in terms of their relationship between RMSD values and simulation time, shown in three different color codes. Panel (**A**) shows an RMSD plot of the backbone atoms of FimH protein receptor interacting with three ligand candidates, while panel (**B**) represents the RMSD plots of ligand candidates during simulation, which are color coded black (Formononetin), red (Quinidine), and green (Pteryxin). FimH protein complexes with Formononetin and Quinidine possessed high stability during simulation, while the third complex with Pteryxin demonstrated comparatively less stability.

Consequently, the RMSDs of the formononetin–FimH and quinidine–FimH complexes were found to have comparatively more stable trajectories than the (+)–pteryxin–FimH complex during simulation analysis (Figure 6).

#### *3.3. Principle Component Analysis*

To identify the overall patterns of motion in the complexes during simulation, we utilized principal components analysis (PCA) for the three complexes. The covariance matrices of the simulations were calculated to produce the eigenvectors, and then a screening of the trajectories was performed around each of the diverse eigenvectors. Furthermore, the dominant motions during a four-vector RMSD for each complex were observed (Figure 7). A large portion of the overall variations in the receptor protein can be explained by a few reduced amplitude eigenvectors with large eigenvalues. The analysis illustrated that the fourth eigenvector accounts for the depletion of the RMSD value, although the rest of the three vectors showed elevation at the initial time of the simulation. It was shown that after ~8 ns the RMSD values reached a plateau and maintained their value until the end of the simulation (Figure 7).

To know the movement of the backbone atoms of the receptor protein, 100 frames of the first principal component were collectively loaded into VMD, which show the motility of the residues (Figure 8), although there might be differences between each complex with respect to the motions sampled. The width of the band is proportional to amplitude; a narrow band signifies those segments that hardly moved, while wide bands signify the sections most affected by the transitions. The middle domain of the protein showed less movement than the adjacent domains that predicted the stability of the middle domain.

The rest of the protein exhibited comparatively wider bands, which might be due to an increase in functionality.

**Figure 7.** Each panel (**A**–**C**) reports the calculation of projection of particular trajectory on eigenvectors for each three complexes. The RMSDs for every four eigenvectors are depicted as four subplots in each panel. Plot shows the projections of each vector are fitted to the structure in the eigenvector.

**Figure 8.** The motion of the backbone atoms of the protein receptor during simulation of the trajectory as calculated by principal component analysis. The width of the bands is proportional to amplitude.

It is worth mentioning that the binding energy values and MD simulations obtained can only demonstrate the binding efficacy and stability of inhibitors with the target protein. Further benchwork experiments are required to optimize these compounds (formononetin, (+)-pteryxin, and quinidine) as promising SGLT2 inhibitors with add-on FimH inhibition potential that might reduce the probability of uropathogenic side effects.

#### **4. Conclusions**

SGLT2 inhibitors are a newer class of drugs that have enormous anti-diabetic potential. Unfortunately, the side effects, specifically those that contribute to the development of UTIs, have impeded their success rate. Several clinicians, as well as diabetic patients, are still in a dilemma over the use of this class of anti-diabetic medication. In contrast, FimH is blooming as a potential alternative target for UTI treatment, and FimH inhibitors are currently under development. These inhibitors act against the uropathogenic bacterial strains adherence to the mucosal surface of the urinary tract, resulting in a lower chance of encountering resistance. In the present study, we docked natural SGLT2 inhibitors with FimH target proteins to predict their FimH interaction potential. To the best of our knowledge, this is the first time that SGLT2 inhibitors have been explored in terms of their use against FimH. Our findings suggest that among all SGLT2 inhibitors examined, formononetin,

(+)-pteryxin, and quinidine exhibited strong interactions with the FimH protein with high binding energy, in comparison to the positive control heptyl α-D-mannopyranoside. On the other hand, the FDA-approved SGLT2 inhibitor canagliflozin indicated lower interaction with FimH. Thus, we hypothesize that if explored further, natural SGLT2 inhibitors might reduce the probability of UTIs when used against FimH. The authors anticipate researchers to duly use the findings of this study to design more versatile and potent dual inhibitors against SGLT2 and FimH to cope with the uropathogenic side effects of SGLT2 inhibitor class anti-diabetic medications.

**Author Contributions:** Conceptualization, S.S. and S.M.D.R.; methodology, S.S., N.C., and S.M.D.R.; writing—Original draft preparation, S.S. and S.M.D.R.; writing—Review and editing, M.M.M., Q.A., S.A., M.M.B., K.A., S.S.A., M.M.A.B., A.M., A.I.A., and M.M.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Deanship of Scientific Research, Najran University, Najran, Saudi Arabia grant number NU/MID/17/091.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Doxorubicin Encapsulation in Carbon Nanotubes Having Haeckelite or Stone–Wales Defects as Drug Carriers: A Molecular Dynamics Approach**

**Leonor Contreras 1,\* , Ignacio Villarroel <sup>2</sup> , Camila Torres <sup>2</sup> and Roberto Rozas <sup>1</sup>**

**Abstract:** Doxorubicin (DOX), a recognized anticancer drug, forms stable associations with carbon nanotubes (CNTs). CNTs when properly functionalized have the ability to anchor directly in cancerous tumors where the release of the drug occurs thanks to the tumor slightly acidic pH. Herein, we study the armchair and zigzag CNTs with Stone–Wales (SW) defects to rank their ability to encapsulate DOX by determining the DOX-CNT binding free energies using the MM/PBSA and MM/GBSA methods implemented in AMBER16. We investigate also the chiral CNTs with haeckelite defects. Each haeckelite defect consists of a pair of square and octagonal rings. The armchair and zigzag CNT with SW defects and chiral nanotubes with haeckelite defects predict DOX-CNT interactions that depend on the length of the nanotube, the number of present defects and nitrogen doping. Chiral nanotubes having two haeckelite defects reveal a clear dependence on the nitrogen content with DOX-CNT interaction forces decreasing in the order 0N > 4N > 8N. These results contribute to a further understanding of drug-nanotube interactions and to the design of new drug delivery systems based on CNTs.

**Keywords:** carbon nanotubes; Stone–Wales defects; haeckelite defects; doxorubicin encapsulation; drug delivery system; binding free energies; noncovalent interactions; molecular dynamics

#### **1. Introduction**

Doxorubicin (DOX), an antineoplasic drug, approved for medical use by the FDA [1,2], has been used for more than 40 years to combat various types of cancers despite the cardiological risks associated with its use. Researchers at the Mayo Clinic [3] have reviewed its mechanism of action and, mainly thanks to the work of Denard et al. and Zhang et al. [4,5], have proposed two alternatives that explain its function. In one of them, DOX would stabilize a complex formed by double-stranded DNA and topoisomerase, which later it would cut both strands of DNA. The alternative is the production of a larger quantity of ceramides which would produce the translocation of a CREB3L-1 protein from the endoplasmic reticulum to the Golgii apparatus. There, some proteases would break the CREB3L-1 protein in such a way that its N-terminal fragment would be translocated to the nucleus where it would direct the DNA transduction, to finally express p21 proteins, which would be those that inhibit tumor growth. Other mechanisms of action of DOX reviewed by Ferreira et al. consider the intercalation of DOX in nuclear DNA and mitochondrial DNA, inhibition of topoisomerase-IIß, and epigenetic factors that involve methylation and deacetylation reactions [6].

In order to increase drug bioavailability and avoid its adverse effects, several types of drug carriers have been used, among which carbon nanotubes (CNTs) have shown

**Citation:** Contreras, L.; Villarroel, I.; Torres, C.; Rozas, R. Doxorubicin Encapsulation in Carbon Nanotubes Having Haeckelite or Stone–Wales Defects as Drug Carriers: A Molecular Dynamics Approach. *Molecules* **2021**, *26*, 1586. https:// doi.org/10.3390/molecules26061586

Academic Editors: Marco Tutone and Anna Maria Almerico

Received: 2 March 2021 Accepted: 11 March 2021 Published: 13 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to form stable associations with DOX [7]. However, CNTs also exhibit some toxicity problems [8–11]. Fortunately, CNTs can be functionalized with fragments to increase their water solubility, which prevents them from being deposited as agglomerates in the body. Under certain conditions of concentration and purity CNTs are non-toxic [12]. Additionally, the functionalization of the CNTs facilitates the anchoring of the nanotubes right in the tumor to be attacked. The physicochemical and conductive properties of CNTs give them a versatility of applications in various fields such as electronics, photonics, catalysts, drug carriers, biotechnology, bone tissue engineering and others [13–15].

In the current work we are interested in the ability of CNTs to adsorb drugs and transport them to the target site. It is important to determine the structural parameters that facilitate the DOX-CNT association and to allow the development of strong DOX-CNT intermolecular interaction forces, which help to inhibit the drug from being released before reaching its target. Once there, the acidic pH of the tumor environment causes the release of the drug. Indeed, several studies show experimentally that DOX release is favored at pH of 5 or less [16–19], which is also demonstrated at a theoretical level [20].

Another technologically important characteristic of CNTs is their chirality which significantly modify their conductive properties. For example, armchair (n, n) nanotubes are conductive while zigzag (n, 0) and chiral (n, m) nanotubes are semiconductors except when the value of the difference (n − m) is a multiple of 3, since nanotubes become conductive [21–24]. The diameter of the nanotube has also been shown to be an important structural point since, depending on the diameter, the degree of curvature of the nanotube can be controlled, which seems to be a decisive factor in the stability of the DOX-CNT associations.

As can be deduced, the chirality and the diameter of the nanotubes are properties that determine their behavior and also the presence of structural defects which can change the chirality [25]. The carbon rings of the nanotube that differ from the hexagons are called defects, and depending on their ordering and distribution in the nanotube they confer different properties. For example, the five- and seven-membered rings, when distributed around the perimeter of the nanotube, constitute defects called bumpy. If they are distributed axially in the nanotube they are called zipper defects [26]. Other fourmembered rings, along with eight-membered rings, are called haeckelite defects [22,27]. These defects are formed by the addition of two carbon atoms or ad-dimers. However, there are other defects that are formed by rearrangement of their bonds, such as the Stone–Wales defects formed by a pair of rings of five and eight members [28].

The presence of some of these defects has significant technological importance. For example, zigzag nanotubes in the presence of ad-dimers can induce plastic transformations in a material that would otherwise be brittle [29]. Chiral nanotubes stand out, which in the presence of bumpy defects considerably increase their conductivity over armchair and zigzag nanotubes according to DFT studies considering dispersion-corrected B3LYP-D3 functional [30]. Zigzag nanotubes that contain bumpy defects show greater conductive ability and capacity to be reduced [30]; nitrogen doping increases the conductive ability of armchair nanotubes [30]. In addition, armchair nanotubes with bumpy defects notably increase their ability to adsorb hydrogen with very convenient hydrogen adsorption energy values, for their use in the management of clean energy [30].

CNTs form stable associations with doxorubicin [7]. Various theoretical and molecular dynamics studies predict a high capacity of CNTs as carriers [20,31–36] and several molecular dynamics studies of functionalized DOX-CNT systems with various organic groups have contributed to the study of DOX loading and release [37–43]. Although there are several works on the adsorption of DOX in CNT, there are no studies that report on the optimal structure that a nanotube should have to behave as a DOX nanocarrier. The situation is complicated because there are also no experimental data available on the formation and characterization of DOX-CNT complexes and the determination of their DOX-CNT binding energies. For non-functionalized nanotubes, Wang and Xu [20], systematically studied the adsorption and encapsulation of DOX in armchair nanotubes of different diameters,

using the theoretical methods PM6-DH2 and M06-2X in the ONIOM scheme and found that the diameter of the nanotube at which the best DOX encapsulation occurred was 14 Å and corresponded to (10,10) armchair nanotubes. This same behavior was confirmed through a study of molecular dynamics for armchair, zigzag and chiral nanotubes, finding that the strongest DOX-CNT interactions were produced for 14 Å in diameter nanotubes, regardless of chirality [32]. A different situation occurs in the presence of defects in the nanotube. In the case of bumpy defects, a dependence on chirality is observed, since armchair nanotubes with bumpy defects present weaker DOX-CNT interactions than armchair nanotubes without defects. In contrast, bumpy defects in chiral nanotubes favor the DOX-CNT interaction [32].

Several methods of synthesis of CNTs have been reviewed [44]. However, there is a lack of comparative systematic experimental antecedents on this issue, which makes it possible to pose as a valid hypothesis for a theoretical study that the presence of defects in the nanotube, the type and number of defects and their position, also modify the DOX-CNT association properties, along with the chirality and size of the nanotube.

The previous antecedents also lead us to investigate if there is a general trend for some type of nanotube, for example, with chirality or type of defect that accounts for the degree of DOX-CNT association. Our research questions include: (i) how does structural or nitrogen doping defects affect the ability of CNTs as drug carriers, in this case, DOX. (ii) Is the effect produced by the defects the same, regardless of the chirality and the size of the nanotube? (iii) Is there any type of defect that has better characteristics than the others? (iv) How does it affect the number of defects present?

In this work DOX-CNT binding energies are determined for chiral nanotubes with haeckelite defects (with rings of 4 and 8 carbon atoms) and for armchair and zigzag nanotubes with Stone–Wales defects (SW), by means of the MM/PBSA and MM/GBSA methods implemented in the AMBER program of molecular dynamics.

#### **2. Results**

Below are the results obtained by molecular dynamics (MD) simulation for DOX encapsulation systems in chiral CNTs with haeckelite (Hk) defects and also in armchair and zigzag CNTs with Stone–Wales (SW) defects. The Hk defects consist of a pair of rings of 4 and 8 members each, while the SW defects are made up of a pair of rings of 5 and 7 members each, as shown in Figure 1.

**Figure 1.** Representation of carbon nanotube structural defects and their C–C bond distances in Å. (**a**) haeckelite defect; (**b**) Stone–Wales defect.

#### *2.1. Chiral Nanotubes with Hk Defects*

Chiral nanotubes Ch(13,08) with one Hk defect (named Hk1) and two Hk defects (named Hk2) having 0N, 4N and 8N were studied considering different initial positions of the DOX: in the region of the defect (D), in the regular region of the nanotube (R) (there

Å

are no defects in that area) with the DOX NH<sup>2</sup> group pointing towards the center of the nanotube (v1 orientation) or to the inverse direction (v2 orientation) as shown in Figure 2 for Hk2 chiral nanotubes. Other additional DOX orientations refer to Hk1 nanotubes: when the DOX NH<sup>2</sup> group is oriented in a direction proximal to the defect (p) or is in the direction opposite to the defect (o) as shown in Figure 3.

Å

**Figure 2.** DOX orientation when encapsulated in a chiral nanotube having two haeckelite defects. (**a**) v1 orientation; (**b**) v2 orientation. Lateral views.

**Figure 3.** The encapsulated DOX orientations into a chiral nanotube having one haeckelite defect. (**a**) Hk1-DoxDIn.v1p; the DOX NH<sup>2</sup> group is oriented in a direction proximal to the defect; (**b**) Hk1-DoxDIn.v1o; the DOX NH<sup>2</sup> group is oriented in an opposite direction. Lateral and frontal views shown.

2.1.1. Chiral Nanotubes with HK1 Defects

For chiral nanotubes with Hk1 defects (Hk1 chiral nanotubes), the results predict a similar behavior for both undoped, 0N, and nitrogen doped nanotubes having 4N and 8N. In all three cases, the DOX-CNT interaction is favored when the DOX is located in the defect area with the NH<sup>2</sup> group pointing towards the center of the nanotube (v1 orientation)

π

Å

in the proximal direction close to the defect as shown in Table 1 (runs 2, 10 and 18 for 0N, 4N and 8N, respectively) with Poisson–Boltzman (PB) binding energies of −102, −99 and −102 kcal/mol, respectively. Coherently, most of these systems exhibit equilibrium distances with values between 3.2 and 3.6 Å evidencing stronger DOX-CNT interactions which are favored by the orientation of the DOX that facilitates the NH-π interaction. In Figure 4 the initial conformation of the Ch(13,08)8N-Hk1-DoxDIn.v1p complex is shown together with the final conformations after 2 ns of MD simulation and after 100 ns. It is observed that DOX does not move towards the regular part of the nanotube but interacts with the defect and as a result, in that area, there is a significant deformation of the nanotube. These results obtained for chiral nanotubes show the same behavior as was reported for armchair nanotubes with more favorable DOX-CNT interactions for systems in which DOX is located in the defect region and when it is oriented with its nitrogen atom directed towards the center of the nanotube. However, armchair nanotubes having one haeckelite defect exhibit DOX-CNT binding energies that are more exothermic suggesting stronger DOX-CNT interactions [33].

**Table 1.** DOX-CNT Poisson–Boltzman (PB) and generalized bond (GB) binding energies (in kcal/mol) for the nitrogen doped and undoped chiral nanotubes (34 Å length) having one haeckelite defect, Hk1, considering encapsulated system. dp-NT are the distances from the DOX anthraquinonic plane to the nanotube wall; dN-NT is the distance from the DOX nitrogen atom to the nanotube wall. All distances are expressed in Å.

<sup>1</sup> DoxD means DOX position is in the defect zone; DoxR is for the DOX in the regular part of the nanotube; v1 means the nitrogen atom of the DOX is oriented towards the center of the tube; v2 indicates the inverse orientation; v1p indicates that the DOX nitrogen atom is located in a proximal space regarding the defect meanwhile v1o is used to indicate that the DOX nitrogen atom is located in an opposite space regarding the defect as shown in Figure 3.

DOX-CNT systems, doped with 4N and containing Hk1 defect exhibit quite similar DOX-CNT PB binding energy values between −79 and −77 kcal/mol, probably accounting for an electronic distribution that interacts with the drug in a similar way regardless of DOX position and orientation. This could be due to the arrangement of nitrogen atoms which are part of two pyrimidine rings placed opposite each other on the walls of the nanotube.

**Figure 4.** Representation of encapsulated DOX-CNT complex for nitrogen-doped chiral nanotube with one haeckelite defect, Ch(13,08)8N-Hk1-DoxDIn.v1p, at different simulation times. (**a**) 0 ns; (**b**) 2 ns; (**c**) 100 ns. Two laterals and one frontal view are shown.

#### 2.1.2. Chiral Nanotubes with Hk2 Defects

Carbon nanotube diameter effect. The nitrogen doped and undoped Hk2 chiral (13,08) CNTs of 14 Å diameter and 19 Å length showed better DOX-CNT PB and GB binding energies than the corresponding Hk2 chiral (13,10) CNTs of 16 Å diameter (calculated with RESP charges for DOX) as shown in Figure 5. This was an expected result considering PM6-DH2 and M06-2X theoretical calculations in the scheme of ONIOM for the DOX encapsulation in armchair CNTs without defects [20], and also molecular dynamics studies on armchair, zigzag and chiral nanotubes with reported values of 14 Å as an optimal value of the nanotube diameter for encapsulating the DOX [32]. A diameter of 14 Å allows the proper curvature of the nanotube for the formation of different attractive and complementary non-covalent interactions between the nanotube and the DOX that stabilize the entire system which is also fulfilled in this case of nanotubes containing two haeckelite defects in their structure.

Hk2 chiral nanotubes of diameter 14 Å exhibit less favorable DOX-CNT PB binding energies for nitrogen doped nanotubes in the order 0N > 4N > 8N with values of −101, −97 and −74 kcal/mol, respectively (see Figure 5), which predicts stronger DOX-CNT interactions for undoped Hk2 chiral (13,08) nanotubes. These values were calculated using RESP charges for DOX. RESP (restrained electrostatic potential) approach to derive partial charges has been reported as having a lower average error than MM3 and CHARMm in a study considering 55 molecules [45]. Mean SD for the PB binding energy values between 2.6 and 3.2 kcal/mol and between 2.6 and 3.1 kcal/mol for the GB binding energy values were observed, being in all cases less than 4.2%.

**Figure 5.** DOX-CNT binding energies for nitrogen-doped and undoped chiral nanotubes of 19 Å length with different diameter, having Hk2 defects and calculated with RESP charges for DOX. (**a**) PB binding energies; (**b**) GB binding energies. Blue is for 14 Å diameter and red is for 16 Å diameter.

Carbon nanotube length and DOX pose effects. Longer Hk2 chiral nanotubes (34 Å length) exhibit more exothermic DOX-CNT PB binding energy values than Hk2 chiral nanotubes of 19 Å length, for both nitrogen doped and undoped nanotubes. For the undoped nanotubes (0N) and those doped with 4N, a clear preference of the DOX for the nanotube defect zone is shown with DOX-CNT PB binding energy values of −109 and −104 kcal/mol, respectively (runs 1 and 5, Table 2). When the DOX is in the regular zone of the nanotube, no significant differences are observed between the v1 or v2 orientations of the DOX. However, in cases where the DOX is initially located in the defect zone, stronger interactions are predicted for v1 DOX orientation for both undoped and 4N doped systems as shown in Table 2 (runs 1 vs. 2 and 5 vs. 6).

**Table 2.** Poisson–Boltzman (PB) and generalized bond (GB) binding energies (in kcal/mol) for the nitrogen doped and undoped chiral nanotubes (34 Å length) having two haeckelite defects, Hk2, considering encapsulated system. dp-NT are the distances from the DOX anthraquinonic plane to the nanotube wall; dN-NT is the distance from the DOX nitrogen atom to the nanotube wall. All distances are expressed in Å.

<sup>1</sup> DoxD means DOX position is in the defect zone; DoxR is for the DOX in the regular part of the nanotube; v1 means the nitrogen atom of the DOX is pointing towards the center of the tube; v2 indicates the inverse orientation.

The best system for DOX encapsulation in Hk2 chiral nanotubes is therefore Ch(13,08)0N-HK2-DoxDIn.v1 (run 1, Table 2) with the DOX in the defect region and v1 orientation (with the nitrogen pointing towards the center of the nanotube). In this conformation, the formation of the non-covalent DOX-CNT interactions is facilitated, which are mainly constituted by π–π stacking interactions complemented by NH-π, CH-π, C=O-π and van der Waals interactions [20,31,46]. Figure 6 shows the non-covalent interactions for the most favorable case with the DOX in the area of the nanotube defect (a large green surface is observed) and as a comparison, the same nanotube with the DOX encapsulated in the regular area (less green regions and more red regions are observed), with PB DOX-CNT binding energy values of −109 and −80 kcal/mol, respectively (runs 1 and 3, Table 2). A program specially developed for the visualization of non-covalent interactions (NCI) was used [47].

Under similar MD simulation conditions but considering RESP charges for DOX, Hk2 chiral (13,08) nanotubes (−101 kcal/mol, Figure 5) predict stronger DOX-CNT interactions than reported Hk2 armchair (10,10) nanotube (−83.4 kcal/mol) with the encapsulated DOX located in the defect zone in v1 orientation [33]. Hk2 chiral nanotubes predict stronger DOX-CNT interactions than reported Hk2 zigzag (18,0) nanotubes which exhibit values of DOX-CNT PB binding energies of −78.7 kcal/mol for undoped nanotubes [33]. The three types of CNTs in comparison have diameters of 14 Å. In this way, in terms of chirality and according to the indicated results, Hk2 nanotubes exhibit the following order of ability to encapsulate the DOX: chiral > armchair > zigzag, despite chiral nanotubes are shorter than zigzag and armchair nanotubes (19 vs. 34 Å length). The enhanced ability of chiral nanotubes with respect to other nanotubes to encapsulate DOX was reported also for perfect CNTs through MD simulation studies considering RESP charges for DOX [32].

**Figure 6.** Representation of non-covalent interactions (NCI) for DOX encapsulation inside chiral carbon nanotubes having two haeckelite defects. (**a**) DOX D position; (**b**) DOX R position. Cutplot 0.01 0.1. Blue surfaces indicate strong interactions; green means weak interactions; red means repulsion [47]. Different comparative lateral and frontal views shown.

#### *2.2. Nanotubes with Stone–Wales Defects*

The encapsulation of the DOX was studied in armchair and zigzag nanotubes that have one and two Stone–Wales defects (SW1 and SW2, respectively) as shown in Figure 7.

**Figure 7.** Representation of carbon nanotubes having one and two Stone–Wales defects. (**a**) SW1 armchair; (**b**) SW2 armchair; (**c**) SW1 zigzag; (**d**) SW2 zigzag. Side and front views shown.

#### 2.2.1. SW1 and SW2 Armchair Nanotubes

Armchair (10,10) nanotubes of 20 Å and 34 Å in length having Stone–Wales defects were studied, which showed different behaviors in DOX encapsulation. The shorter nanotubes (20 Å length) exhibit a significant stronger interaction with the DOX in two situations: (i) when they are of the SW2 type (with two defects, doped and undoped) in comparison with SW1 nanotubes and (ii) when they are doped with 4N (SW1 and SW2) as is clearly depicted in Figure 8a. Nitrogen doped SW1 and SW2 armchair (10,10) nanotubes of 20 Å length predict stronger interactions with the DOX than corresponding longer SW1 and SW2 nanotubes of 34 Å in length as shown in Figure 8 particularly for 4N doped nanotubes.

**Figure 8.** DOX-CNT PB binding energies for nitrogen-doped and undoped armchair (10,10) nanotubes of different length having one and two Stone–Wales defects. (**a**) 20 Å length; (**b**) 34 Å length.

In contrast, longer SW1 armchair nanotubes predict somewhat stronger interactions than longer SW2 armchair nanotubes the difference being more significant for the undoped SW1 and SW2 nanotubes with PB DOX-CNT binding energy values of −105 and −80 kcal/mol for, respectively, with DOX v1 orientation, and for DOX v2 orientation −92 vs. −81 kcal/mol, respectively, as shown in Figure 8b.