Potential Stereoselective Binding of Trans-(±)-Kusunokinin and Cis-(±)-Kusunokinin Isomers to CSF1R

Breast cancer cell proliferation and migration are inhibited by naturally extracted trans-(−)-kusunokinin. However, three additional enantiomers of kusunokinin have yet to be investigated: trans-(+)-kusunokinin, cis-(−)-isomer and cis-(+)-isomer. According to the results of molecular docking studies of kusunokinin isomers on 60 breast cancer-related proteins, trans-(−)-kusunokinin was the most preferable and active component of the trans-racemic mixture. Trans-(−)-kusunokinin targeted proteins involved in cell growth and proliferation, whereas the cis-(+)-isomer targeted proteins involved in metastasis. Trans-(−)-kusunokinin targeted CSF1R specifically, whereas trans-(+)-kusunokinin and both cis-isomers may have bound AKR1B1. Interestingly, the compound’s stereoisomeric effect may influence protein selectivity. CSF1R preferred trans-(−)-kusunokinin over trans-(+)-kusunokinin because the binding pocket required a ligand planar arrangement to form a π-π interaction with a selective Trp550. Because of its large binding pocket, EGFR exhibited no stereoselectivity. MD simulation revealed that trans-(−)-kusunokinin, trans-(+)-kusunokinin and pexidartinib bound CSF1R differently. Pexidartinib had the highest binding affinity, followed by trans-(−)-kusunokinin and trans-(+)-kusunokinin, respectively. The trans-(−)-kusunokinin-CSF1R complex was found to be stable, whereas trans-(+)-kusunokinin was not. Trans-(±)-kusunokinin, a potential racemic compound, could be developed as a selective CSF1R inhibitor when combined.


Introduction
Breast cancer is one of the leading causes of mortality in women worldwide. In 2020, breast cancer was estimated at 2.26 million new cases and 0.68 million deaths [1]. The common therapeutic options for breast cancer patients are surgery, radiotherapy, chemotherapy, or hormonal therapy, in combination with targeted medication [2]. Despite the more accurate methods of treatment and diagnosis being developed, the incidence and death rate of breast cancer continues to rise [1][2][3][4]. Due to the drug resistance mechanism and the high heterogeneity of breast cancers, additional targets are needed to be explored and highly specific drugs are continuously in need of investigation [2,[4][5][6][7]. The development of more effective targeted therapeutic approaches could serve as a new paradigm for breast cancer treatment in the personalized medicine era.
Despite the fact that the cis-(±)-kusunokinin isomers have a significant effect on biological activity, these isomers remain an intriguing area for drug development.
In this study, we identified the target of four kusunokinin enantiomers (trans-(−), trans-(+), cis-(−) and cis-(+)-isomer) on 60 candidate proteins related to the CSF1R-breast cancer progression pathway and proposed a preferable form of kusunokinin isomers using molecular docking and our scoring system. The potential target protein was then simulated via molecular dynamics (MD) for trans-(−)-kusunokinin or trans-(+)-kusunokinin to determine the binding mode and stereoisomeric selectivity of the target to each isomer within the trans-racemic mixture.

Target Protein(s) Screening of Four Kusunokinin Isomers among CSF1R-Related Breast Cancer Progression Proteins
To investigate the possible target, binding ability and binding character of each kusunokinin stereoisomer, 60 candidate proteins related to CSF1R-associated breast cancer progression were selected for a molecular docking study. The proteins were grouped based on their primary mechanism: anti-apoptosis and survival (7 proteins), cell growth and proliferation (21 proteins), metastasis (11 proteins) and signal molecules (21 proteins). All PDB codes of the candidate proteins and PubChem CIDs for all given inhibitors were summarized in Table S1. The docking score in kcal/mol of each kusunokinin isomer and a known inhibitor of each protein are shown in Table 1. The target selection criteria were based on a protein(s) that (1) had a lower docking score than a known inhibitor and (2) had a ligand that was located at a similar site to the known inhibitor. The docked pose and the aforementioned criteria were used to determine each isomer's potential inhibition. Trans-(−)-kusunokinin displayed a lower docking score than known inhibitors in 37 out of 60 proteins. However, trans-(−)-kusunokinin showed similar binding sites to the known inhibitor of 12 proteins, including anti-apoptosis and survival (COX2, Hsp90a and Hsp90b), cell growth and proliferation (cdc25A, CDK2, RagC, CSF1R and EGFR) and metastasis (ALP, PALP, AKR1B1, JAK3, MNK2 and p38) (docking score of −11.84 to −6.92 kcal/mol). The binding ability of trans-(−)-kusunokinin with COX2, Hsp90b, RagC, CSF1R, EGFR and p38 was the highest among the four kusunokinin isomers. Trans-(−)kusunokinin had the greatest potential on CSF1R, with a docking score of −11.84 kcal/mol, whereas the CSF1R known inhibitor, pexidartinib, had a docking score of −11.59 kcal/mol.
Meanwhile, trans-(+)-kusunokinin showed a lower docking score than known inhibitors in 28 proteins. Seven proteins passed the potency inhibition criteria, namely COX2, Hsp90a (anti-apoptosis and survival), CDK2 and EGFR (cell growth and proliferation), ALP, PALP and AKR1B1 (metastasis), JAK3 and MNK2 (signal molecules). However, when compared to other kusunokinin stereoisomers, the trans-(+)-isoform did not show the best binding possibility in any of these proteins. Trans-(+)-kusunokinin showed the highest possibility on only AKR1B1, with a docking score of −11.15 kcal/mol, whereas epalrestat (known inhibitor) had a docking score of −10.14 kcal/mol.

Docking Score Interpretation of the Kusunokinin Isomers
The docking scores of each kusunokinin isomer were re-analyzed in rank scores ranging from 4 to 1 using our criteria. Score 4 represented the highest binding affinity of the four isomers, followed by 3 and 2 and 1 was the lowest. The kusunokinin isoform rank scores were then multiplied by the protein's weight score, which ranged from 3 to 1 based on drug-accessibility prediction. Number 3 stood for membrane-bound proteins, which have the highest chance of drug exposure. Number 2 represented cytoplasmic proteins, whereas number 1 represented proteins inside the nucleus, which have the lowest chance of drug exposure ( Figure 2) [36,[48][49][50][51]. The most preferred form was determined by adding the weighted rank scores within the protein group, which is represented as a sum score (weighted) in Table 2. Trans-(−)-kusunokinin presented the highest sum score and weighted sum score in all four protein groups including anti-apoptosis and survival, cell growth and proliferation, metastasis and signaling proteins with a weighted score of 54, 151, 52 and 148, respectively.
Meanwhile, trans-(+)-kusunokinin had the lowest sum score in all four protein groups and the lowest weighted sum score in all groups except signaling molecules.
The sum scores of cis-(−)-and cis-(+)-isomers were not the highest in any protein groups. However, the cis-(+)-kusunokinin isomer gave the highest weighted sum score in the metastasis proteins (weighted score of 56). No significant variation in the weighted sum scores of cis-(−)-and cis-(+)-kusunokinin isomers was observed in the anti-apoptosis and survival or the cell growth and proliferation proteins. Table 2. Candidate proteins with rank score and weight score.

Target Name
Weight 1 Kusunokinin Isomer Rank Score 2

Signaling molecules
1 Score for Ligand Accessibility (3 for membrane-bound proteins, which were most likely to contact with the drug, 2 for cytoplasm proteins and 1 for proteins inside nucleus which were the least exposed). 2 According to Table 1, the ranking score for kusunokinin stereoisomers is: 4 for the best binding ability prediction, 3 for the second, 2 for the third and 1 for the worst binding ability prediction. The blanket's number represented the rank score multiplied by the protein weight score. 3 The sum score represents the total sum of the rank scores of each isoform within the same protein group. The blanket number represented the sum of rank scores multiplied by the weight score of each protein.
Based on our criteria, we proposed that trans-(−)-kusunokinin was the best stereoisomer for anti-apoptosis and survival, cell growth and proliferation, metastasis and signaling proteins, whereas cis-(+)-isomer was a potential candidate for the inhibition of metastasis proteins.

Target Selection Based on Ligand-Protein Interaction Analysis
The binding energies of some proteins differed between four kusunokinin isomers, as shown in Table 1, whereas others did not. One-way ANOVA was used to test the variability of the docking score among isomers within the same protein (Table S2). Of the 12 possible targets, CSF1R had the most variety (SD = 1.089, p < 0.01), followed by Hsp90b and cdc25A (SD = 0.52 and 0.48, p < 0.01, respectively). Although the variances of EGFR, MNK2, JAK3 and AKR1B1 were all within the same range (SD = 0.032 to 0.036, p < 0.01), EGFR had the smallest differences between the binding energies of the four kusunokinin isomers. The 12 potential targets' binding pockets were visualized for further investigation. Figure 3 depicts the docking poses of four kusunokinin isomers to each protein binding site.
With the exception of CSF1R and EGFR, we discovered that all 12 target pockets were relatively round and had fewer aromatic residues than CSF1R, allowing kusunokinin stereoisomers to bend to their preferred poses and fit in more freely than CSF1R, which had a narrow binding site ( Figure 3). Although both CSF1R and EGFR have a "glove-like" pocket ( Figure 3M), CSF1R's narrow pocket forced kusunokinin stereoisomers into a planar arrangement in order to fit in the binding site ( Figure 3A), whereas EGFR's binding site was wider than CSF1R's, allowing for binding arrangements without stereoselectivity ( Figure 3A,L). The protein selectivity was explained by the narrow, glove-shaped CSF1R binding pocket. CSF1R and EGFR were chosen as models for studying kusunokinin stereoisomeric effects because they had the highest selectivity toward the trans-(−)-isoform, whereas EGFR had no selectivity to any kusunokinin stereoisomer. CSF1R and EGFR were tyrosine kinase receptors that regulated breast cancer cell growth and proliferation [52][53][54]. Both targets revealed that the best binding possibility was to the trans-(−)-isomer, followed by the cis-(−)-, cis-(+)-and trans-(+)-isomers. Figures 4 and 5 show the CSF1R and EGFR binding positions of four kusunokinin isoforms and a known inhibitor, respectively. Figure S1 depicts the alignment of all docked ligand poses at the binding pocket of CSF1R or EGFR. The docked poses of pexidartinib-CSF1R and erlotinib-EGFR were found to be in good agreement with the experimental crystal structure, as shown in Figure S2. In addition, the interaction of all kusunokinin isomers as well as the inhibitor with CSF1R and EGFR was summarized in Table 3. (E) pexidartinib (red). The green dash lines represented hydrogen bonding. Trp550, Tyr665 and Phe797 were responsible for π-π interaction (black dash lines). Docking poses of the ligands and amino acids interacting with the ligands through hydrophobic interaction were shown as sticks. Table 3. Type of interaction between kusunokinin isomers/inhibitor with CSF1R and EGFR.  CSF1R provided π-π stacking with trans-(−)-kusunokinin, cis-(−)-isomer and pexidartinib via Trp550 ( Figure 4A,C,E) and hydrogen bonds with Arg549 and at the same time forced trans-(+)-kusunokinin and cis-(+)-isomer to make contacts with outer residues ( Figure 4B,D). On the contrary, hydrogen bonds were the most common type of interaction found between EGFR residues and ligands. The four isoforms and a known inhibitor, erlotinib, bind to EGFR in the same pocket and share Phe856 as a mutual residue ( Figure 5). The EGFR binding site was larger and contained fewer aromatic residues than the CSF1R binding site, allowing all four kusunokinin structures to interact with no stereo restriction.
A comparison of two trans-kusunokinin isomers in EGFR binding was also described. Both ligands interacted with the same residues and had a strikingly similar binding position ( Figure 7A). Erlotinib was also found to be related to trans-(−)-kusunokinin. Both ligand binding pockets were discovered to be identical but in different positions. However, when it came to forming hydrogen bonds with the residues Lys754, Asp855 and Phe856, the two ligands interacted similarly ( Figure 7B). Interestingly, although erlotinib was found to form π-sulfur with Met766 and Met790 residues, the binding affinity was not as high as that of kusunokinin isomers.

Molecular Dynamic Simulations and Relative Binding Affinities of Trans-(±)-Kusunokinin and Pexidartinib to CSF1R
According to previous findings, CSF1R clearly prefers trans-(−)-kusunokinin over trans-(+)-kusunokinin. Therefore, the two trans isoforms were chosen for MD simulations on CSF1R in order to assess the stereoisomeric effect on protein selectivity. The structures of trans-(−)-kusunokinin and trans-(+)-kusunokinin are available in the PubChem database and have been tested for anticancer activity against breast cancer cell lines, whereas the cis-(±)-isomers are not yet available. MD simulations were run on trans-(−)-kusunokinin-CSF1R, trans-(+)-kusunokinin-CSF1R, pexidartinib-CSF1R and ligand-free CSF1R models in isotonic 0.15 M NaCl solution at pH 7 and 25 • C in our study. The residue number was rearranged using the AMBER20 procedure and the reference residue number is shown in Table S3.

Molecular Dynamic Simulations and Relative Binding Affinities of Trans-(±)-Kusunokinin and Pexidartinib to CSF1R
According to previous findings, CSF1R clearly prefers trans-(−)-kusunokinin over trans-(+)-kusunokinin. Therefore, the two trans isoforms were chosen for MD simulations on CSF1R in order to assess the stereoisomeric effect on protein selectivity. The structures of trans-(−)-kusunokinin and trans-(+)-kusunokinin are available in the PubChem database and have been tested for anticancer activity against breast cancer cell lines, whereas the cis-(±)-isomers are not yet available. MD simulations were run on trans-(−)-kusunokinin-CSF1R, trans-(+)-kusunokinin-CSF1R, pexidartinib-CSF1R and ligand-free CSF1R models in isotonic 0.15 M NaCl solution at pH 7 and 25 • C in our study. The residue number was rearranged using the AMBER20 procedure and the reference residue number is shown in Table S3.

Molecular Dynamic Simulations and Relative Binding Affinities of Trans-(±)-Kusunokinin and Pexidartinib to CSF1R
According to previous findings, CSF1R clearly prefers trans-(−)-kusunokinin over trans-(+)-kusunokinin. Therefore, the two trans isoforms were chosen for MD simulations on CSF1R in order to assess the stereoisomeric effect on protein selectivity. The structures of trans-(−)-kusunokinin and trans-(+)-kusunokinin are available in the PubChem database and have been tested for anticancer activity against breast cancer cell lines, whereas the cis-(±)-isomers are not yet available. MD simulations were run on trans-(−)-kusunokinin-CSF1R, trans-(+)-kusunokinin-CSF1R, pexidartinib-CSF1R and ligand-free CSF1R models in isotonic 0.15 M NaCl solution at pH 7 and 25 • C in our study. The residue number was rearranged using the AMBER20 procedure and the reference residue number is shown in Table S3.
The root-mean square displacement (RMSD) of the 300 ns-CSF1R MD simulations was calculated to assess a structure's conformational stability during the simulation. All MD simulations were stable, according to the RMSD plot ( Figure 8A). The structure flexibility of ligand-bound and ligand-free CSF1Rs was also interpreted using the root-mean square fluctuation (RMSF) ( Figure 8B). Each amino acid α carbon atom in the CSF1R structure was used to create the RMSF. To detect and compare the effect of each bound ligand on protein flexibility, the patterns of four RMSF models were plotted together. The flexibility patterns of trans-(−)-kusunokinin differed from those of ligand-free CSF1R from residues 10th to 50th and 130th to 150th. The overall flexibility patterns of trans-(−)-kusunokinin, trans-(+)-kusunokinin and pexidartinib were distinct, indicating that these three ligands have different effects on CSF1R conformational changes.
The root-mean square displacement (RMSD) of the 300 ns-CSF1R MD simulations calculated to assess a structure's conformational stability during the simulation. All simulations were stable, according to the RMSD plot ( Figure 8A). The structure flexib of ligand-bound and ligand-free CSF1Rs was also interpreted using the root-mean sq fluctuation (RMSF) ( Figure 8B). Each amino acid α carbon atom in the CSF1R struc was used to create the RMSF. To detect and compare the effect of each bound ligan protein flexibility, the patterns of four RMSF models were plotted together. The flexib patterns of trans-(−)-kusunokinin differed from those of ligand-free CSF1R from resi 10th to 50th and 130th to 150th. The overall flexibility patterns of trans-(−)-kusunok trans-(+)-kusunokinin and pexidartinib were distinct, indicating that these three lig have different effects on CSF1R conformational changes.  To investigate protein-ligand dissociation, the distance between the protein's center of mass and the ligand at each time point was presented in Angstrom (Å) and plotted against simulation time ( Figure 9A). The three ligand-CSF1R models were also visualized at 0, 100, 200 and 300 ns ( Figure 9B). The distance analysis of the MD trajectories revealed that as simulation time was increased, trans-(+)-kusunokinin moved away from the binding pocket, whereas trans-(−)-kusunokinin and pexidartinib stably bound CSF1R. These findings indicated that trans-(−)-kusunokinin was the only isoform in the racemic mixtures that bound to CSF1R, not trans-(+)-kusunokinin.
Molecules 2022, 27, 4194 13 of To investigate protein-ligand dissociation, the distance between the protein's center mass and the ligand at each time point was presented in Angstrom (Å) and plotted again simulation time ( Figure 9A). The three ligand-CSF1R models were also visualized at 0, 1 200 and 300 ns ( Figure 9B). The distance analysis of the MD trajectories revealed that simulation time was increased, trans-(+)-kusunokinin moved away from the binding pock whereas trans-(−)-kusunokinin and pexidartinib stably bound CSF1R. These findin indicated that trans-(−)-kusunokinin was the only isoform in the racemic mixtures th bound to CSF1R, not trans-(+)-kusunokinin.  The relative binding free energy was calculated using the MM/GBSA and MM/PBSA methods based on the MD trajectory ( Table 4). The relative binding ability of synthetic trans-(±)-kusunokinin to CSF1R was measured in comparison to the relative binding free energy of pexidartinib. Pexidartinib had the highest binding ability of the selected ligands, followed by trans-(−)-kusunokinin and trans-(+)-kusunokinin. When comparing the two isomers within the trans-racemic mixtures, the results showed that CSF1R preferred trans-(−)-kusunokinin to trans-(+)-kusunokinin.

Conformational Effects of Trans-(±)-Kusunokinin and Pexidartinib on CSF1R
The distance measured from the α carbon of each CSF1R residue to the protein center of mass was expressed as a d value in Angstrom (Å). Then, using the following equation, the difference in distance (∆d) between ligand-bound and ligand-free CSF1R residues was calculated: when d (n, ligand-bound) was a distance between the nth amino acid in ligand-bound CSF1R to the center of mass, and d (n, ligand-free) was a distance between the nth amino acid in ligand-free CSF1R to the center of mass. An interpretation of the ∆d value is as follows: 1.
The ∆d value > 0 implied the residue moved further away from the protein center.

2.
The ∆d value < 0 implied the residue moved closer to the protein center.

3.
The ∆d value = 0 implied the residue stayed in the same position.
The ∆d values that diverged more than 2 Å were considered significant. Figure 10 depicts the ∆d patterns of the three ligand-CSF1R models. Trans-(−)kusunokinin was found to have a contraction effect on the residues 1st and 5th of CSF1R juxtamembrane region (JM), and it simultaneously moved the residues 135th to 145th in the kinase domain away from the reference position. The remaining ∆d patterns were deemed insignificant ( Figure 10A).
Pexidartinib influenced the position of CSF1R residues 59th and 274th, as well as residues 135th to 145th. It also displayed similar distance patterns at the JM and AL as trans-(−)-kusunokinin. Nonetheless, the majority of the ∆d patterns were distinct, implying that both compounds affected CSF1R conformational changes differently, leading to the conclusion that trans-(−)-kusunokinin affected the CSF1R conformation differently than pexidartinib ( Figure 10C).
Furthermore, trans-(−)-kusunokinin binding energies were lower than trans-(+)-kusunokinin binding energies in every protein studied, indicating that the trans-(−)-isomer was a potent form in the racemic compounds [43]. The different actions of extracted trans-(−)kusunokinin and synthetic trans-(±)-kusunokinin were observed and hypothesized to be due to racemic compound nature and purity [41,42]. Because most biological target sites require stereospecific recognition of the binding counterpart, stereochemistry is one of the factors that affect pharmacokinetics and pharmacodynamics [55][56][57]. Different chiral carbon configurations in kusunokinin may result in a stereoisomeric preference for target protein selection and binding mode from all four kusunokinin isomers.
In this study, we first used molecular docking to investigate target proteins of four kusunokinin stereoisomers, followed by molecular dynamics (MD). COX2, Hsp90a, Hsp90b, cdc25A, CDK2, RagC, CSF1R, EGFR, AKR1B1, JAK3, MNK2 and p38 were bound with all four kusunokinin stereoisomers. Trans-(−)-kusunokinin bound CSF1R with the highest specificity, whereas the other three forms bound AKR1B1. These findings suggested that when cells were treated with trans-(+)-kusunokinin, the trans-(−)-isoform had a higher chance of binding CSF1R than other candidates, whereas some trans-(−)-isoform and trans-(+)-kusunokinin may have crossed the cell membrane and inhibited other proteins such as AKR1B1. Figure 11 depicts a pathway containing the 12 proteins with the highest docking score and related molecules.
Considering the docking score alone was insufficient for target protein selection (Table 1) due to that some proteins did not differ in binding energies among the four kusunokinin isoforms. Therefore, we presented a novel sorting method for drug-target screening. The docking score of kusunokinin isomers was recalculated as a rank score from 4 to 1 based on the binding possibility of the isomer. The rank score was then given a weighted score of 3-2-1 based on the protein's likelihood of drug exposure. The weighted score was determined by target subcellular localization, which was an important factor in drug discovery [48][49][50][51]65]. Thus, we divided proteins into three groups and gave a score according to their role: membrane-bound proteins (3), cytosol proteins (2) and nuclear proteins (1). The membrane protein is the most suitable for drug targets that do not require membrane permeation [36,[66][67][68]. Our scoring method facilitated protein selection for in vitro research.
The weighted sum score of trans-(−)-kusunokinin indicated that it was the most potent compound in the racemic mixture, exhibiting high selectivity in binding to proteins involved in cell growth, signaling molecules and anti-apoptosis ( Table 2). The sum and weighted scores of trans-(−)-kusunokinin on CSF1R indicated that the trans-(−)-isoform was the active component in the racemic mixture, which accounts for the decision on CSF1R as the target of interest. Unlike trans-(±)-kusunokinin, the two cis-isomers exhibited a similar score on the weighted score on cell growth and anti-apoptotic proteins. Moreover, cis-(+)-kusunokinin isomers tend to be specific with metastasis proteins more than others. The action of cis-(+)-kusunokinin isomers should be verified by conducting additional research on the metastasis pathway in the near future. Differences in configuration of kusunokinin stereoisomers affected binding poses and distinguishable binding energies. The variability of four kusunokinin binding energies within the same protein was investigated via one-way ANOVA (Table S2). CSF1R had the highest selectivity toward trans-(−)-kusunokinin, whereas EGFR had the smallest difference between the docking score of the four isomers. Although the variances of EGFR, MNK2, JAK3 and AKR1B1 were all within the same range, EGFR was selected for binding pocket comparison to CSF1R because it is a membrane-bound protein from the proliferation and cell growth group, like CSF1R. Both CSF1R and EGFR have "glove-shaped" pockets; however, free binding of stereoisomers was allowed due to wide EGFR's allosteric site and lack of aromatic residues. The stereoselectivity of CSF1R resulted from the pocket shape, which required a planar arrangement of the ligand to fit the pocket. The proteins of all isomers with good binding energies had spherical binding pockets, allowing isomers to bend into their preferred pose with less restriction than the narrow pocket in CSF1R.
The number and position of aromatic residues within the binding site also contribute to kusunokinin isoform stereoselectivity. The bond rotation of the kusunokinin aromatic rings had little effect on the binding energies due to the small number of aromatic residues. The proteins lacking isomeric selectivity relied primarily on hydrogen bonds and hydrophobic interactions. The differences in docking score among isomers with spherical protein pockets were caused by the different binding poses, which formed different hydrogen bonds between the ligand's methoxy groups and appropriate residues within the binding site.
All four isomers bound CSF1R at the allosteric pocket, but they were unable to stretch their formaldehyde ethylene acetal rings into the ATP binding area like pexidartinib's pyrimidine rings. The ligands did, however, bind CSF1R at the comparatively binding site to the crystal structure of CSF1R-pexidartinib, hinting at the possible inhibitory effects. Similar to pexidartinib, trans-(−)-and cis-(−)-isomers formed π-π stacking with Trp550 in JM [42,69]. Trans-(+)-and cis-(+)-isomers, on the other hand, could not fit into the JM region and were thus forced out to form the π-π interaction with the outer residues (Phe797 or Tyr665, respectively). Phe797 at the start of the activation loop (AL) in DFG motifs indicated whether CSF1R was activated (DFG-out) or inactivated (DFG-in) [70]. The π-π interaction of trans-(+)-kusunokinin's benzene rings with Phe797 could lead to an inhibition of phosphorylation of the DFG motif and thus prevent CSF1R activation. However, the interaction of trans-(+)-kusunokinin with CSF1R was found to be less stable than that of trans-(−)-kusunokinin. MD simulation also illustrated that, in the JM region, the π-π interaction between the benzene rings of trans-(−)-kusunokinin and Trp550 stabilized the CSF1R conformation, which was a critical structure for improving potency and/or selectivity towards CSF1R. However, trans-(−)-kusunokinin, trans-(+)-kusunokinin and pexidartinib caused distinct conformational changes in CSF1R. An in vitro study is required to confirm trans-(±)-kusunokinin-CSF1R binding.
Lower capital costs and/or difficulty in chiral synthesis rationalized the usable racemic drug form over a single active enantiomer [56,57,71,72]. Pure enantiomers are usually expensive to chemically synthesize on a commercial scale [72,73]. As a result, the majority of available drugs are marketed as racemates, which are equimolar mixtures of two enantiomers, such as ibuprofen [74]. The main target protein, CSF1R, demonstrated strong selectivity toward the trans-(−) isoform and at the same time forced the trans-(+) isoform out of the pocket, as the MD simulation over time suggested. Therefore, trans-(±)-kusunokinin appeared to be a promising selective CSF1R racemic inhibitor. These statements validated the use of racemic trans-(±)-kusunokinin as a potential candidate or developed scaffold for targeting CSF1R and other proteins in the CSF1R breast cancer progression pathway.

Target Protein Selection and Preparation
The candidate protein was selected based on breast cancer-CSF1R associated pathways, then grouped as anti-apoptosis and survival, cell growth and proliferation, metastasis, or signal molecules. The protein crystal structure, in a PDB format, was acquired from the Protein Data Bank (https://www.rcsb.org, accessed on 19 August 2021). The monomeric form was chosen with co-crystallized ligand(s), solvent molecules and/or other nonprotein molecules removed. Only cofactor molecules such as Zn 2+ and Mg 2+ were retained. The structure preparation was performed using Visual Molecular Dynamics (VMD) [76]. All polar hydrogens were added corresponding to the protonation state at pH 7 using Autodock Auxiliary Tool (ADT) version 1.5.6.

Ligand Preparation
The 3D structure of trans-(−)-kusunokinin (Pubchem CID: 384876), trans-(+)-kusunokinin (Pubchem CID: 92154432), as well as each target protein inhibitor, were retrieved from Pub-Chem (https://pubchem.ncbi.nlm.nih.gov, accessed on 19 August 2021). All retrieved files were in SDF format, and the files were then converted into PDB format via the Online SMILE Translator and Structure File Generator (https://cactus.nci.nih.gov/translate/index.html, accessed on 19 August 2021). For the cis-isomers, both structures were optimized in Open Babel software using a steepest descent algorithm with a convergence of 10 −10 kcal/mol. The ligand structure was finally changed into a PDBQT file format using ADT prior to molecular docking.

Docking Procedure and Binding Energy Calculation
The x-y-z grid point numbers of 120-120-120 with the center at the macromolecule were exploited in the docking study. Since we performed binding docking assumption, these chosen x-y-z grid point numbers had covered the protein receptor, ensuring that randomization was done without bias. The grid spacing was set as 0.375 Å. A genetic algorithm (GA) with 50 runs and a population size of 200 was operated with default parameters in AutoDock4. A Lamarckian genetic algorithm 4.2 was applied to predict bound conformations. The conformation with the lowest docking score in kcal/mol was considered as the best pose. The docking score of kusunokinin stereoisomers and its inhibitor to each target protein were reported and then assessed by the scoring system.

Target Protein Selection for Ligand-Protein Interaction Analysis
The docking score of the target protein known inhibitor, as well as the character of protein binding pocket and the interacted amino acid residues, were considered in this sorting step in order to reduce the number of candidate proteins to the most promising kusunokinin-derived drug targets for future study. A target(s)-choosing criteria for postdocking analysis required that the protein had a lower docking score of kusunokinin isomer than those of known inhibitor and the binding pocket of the said kusunokinin isomer and the known were similar. These criteria are based on a hypothesis that the same binding site and a similar docked pose between the ligand of interest and the known inhibitor could imply the inhibitory effect of the ligand to that protein.

Statistical Analysis
The relative binding energies from the four kusunokinin isoforms within the same protein were assessed by the one-way ANOVA using the Calc package in Ubuntu. p < 0.01 was considered to indicate a statistically significant difference among isomer docking scores.

Scoring Procedure
The scoring procedure was carried out to investigate the preferred form of kusunokinin and its potential target proteins, which could be the cause of the contradictory effects observed in previous studies. The docking score of each kusunokinin stereoisomeric forms and the individual proteins was re-categorized as a rank score in 4-3-2-1 manner from the highest binding possibility (4) to the lowest binding possibility (1).
Weight scores ranging from 3 for transmembrane proteins, 2 for cytoplasmic protein and 1 for nuclear protein, based on drug-accessibility, were applied to the rank score as a multiplier to emphasize the differences in each kusunokinin isoform's binding ability. This criterion was earlier illustrated in Figure 2.
The summations of the rank score and weighted rank score within each protein group were computed and presented as the sum score and weighted sum score, respectively. The weighted sum score indicated the preferred-kusunokinin form and the protein group which most likely to be the target of the kusunokinin most-preferred isoform.

System Preparation of the Compound-Protein Complex
AMBER20 was used for molecular dynamics (MD) simulation to simulate the dynamic aqueous condition as a function of temperature and pressure. The docked PDB structure of the kusunokinin stereoisomer-protein complex was used as a starting point. All of the previous polar hydrogen from AutoDock4 was removed and full hydrogen atoms were added using the AMBER20 package's LeAP module [77]. At pH 7, the corresponding state is all ionizable amino acid side chains. If necessary, the disulfide bond(s) were introduced. TIP3P water solvated the complex at a distance of 14 from the protein surface, resulting in approximately 17,400 water molecules. Na + or Cl − were used to neutralize the solvated system. To simulate a 0.15 M NaCl solution, 32 NaCl pairs were used.

Simulation Protocol of the Compound-Protein Complex
The AMBER20 force field was used to parameterize the protein and solvent. The previous study's kusunokinin RESP charge was used [42]. The protein solution system was energy-minimized for 1000 steps in a periodic boundary condition using steepest descent and conjugate gradient methods. The energy-minimized structure was then equilibrated at 298 K (25 • C) using Langevin Dynamics, also known as the NVT ensemble. Each NVT ensemble lasted 200 ps and had a time step of 1 fs. The system was then placed into an isothermal and isobaric (NPT) ensemble using the Berendsen algorithm [78] for 300 ns at 298 K and 1.013 bar (1 atm) with a time step of 2 fs. The cut-off value of 12 was used to handle non-bonded interactions. The electrostatic interaction was calculated using the PME summation. All H-X bonds were restrained using the SHAKE technique.
Finally, a snapshot MD trajectory of 15,000 equidistant points was used. The first 10,000 snapshots from the 0 to 200 ns MD trajectory were omitted and the last 100 ns of the NPT trajectory was used to calculate the conformational parameters and relative binding energy. VMD was used to visualize the structure and analyze the binding pocket. The RMSD plot acquisition was calculated from VMD and visualized using the Grace package in Ubuntu.

Relative Binding Free Energy Evaluation of the Compound-Protein Complex
The average relative binding free energies (∆G) in kcal/mol were calculated from MD trajectories using molecular mechanics/generalized Born-surface area (MM/GBSA) and molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) [78]. The computation protocol was summarized in the previous study [79].

Conclusions
The molecular simulation and our newly proposed ranking system predicted that trans-(−)-kusunokinin targeted CSF1R, whereas trans-(+)-kusunokinin and cis-(±)-isomer inhibited AKR1B1. In the cell growth and proliferation group, trans-(−)-kusunokinin was the most potent form targeting proteins. The possibility of cis-(±)-isomer targeting proteins in the metastasis group was to be further investigated. Because of the narrow binding pocket and the aromatic Trp550, the stereoisomeric preference of trans-(−)-kusunokinin was observed in CSF1R via selective π-π interaction with a planar arranged trans-(−)kusunokinin. Other targets, such as EGFR binding pockets, were larger than those of CSF1R and lacked aromatic residues, allowing kusunokinin isomers to bend and bind freely. Our study concluded that the racemic trans-(±)-kusunokinin may have inhibited breast cancer cells primarily through the binding and suppression of CSF1R by trans-(−)-kusunokinin, in addition to the binding of trans-(−) or trans-(+) isomers to other target proteins.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/molecules27134194/s1, Table S1: PDB identification codes of CSF1Rderived cancer progression associated proteins and Pubchem CIDs of selected known inhibitors, Table S2: One-way ANOVA of four kusunokinin isomers docking score to the same protein, Table S3: Reference-and MD rearranged-residue number of CSF1R kinase domain, Figure S1: Selected candidate proteins binding pocket alignment, Figure S2: Alignment of docking structure to the reported crystal structure.