Breast cancer is a prevalent malignant illness that primarily affects women and eventually leads to death (Albrand and Terret, 2008). Breast cancer has the second highest fatality rate after lung cancer. Breast cancer is more common in North America and Northern Europe and less common in Asia and Africa (Richie and Swanson, 2003). The female hormone estrogen is essential for the development and progression of breast cancer (Labrie et al., 1999). The proportion of estrogen receptor-a (ERa) expressed in tumor cells is two-thirds of that expressed in normal breast tissues in breast cancer (Chen et al., 2008). ERa-positive breast cancer was frequent in metastatic illness among the breast cancer types. For ERa-positive breast tumors, the estrogen signaling pathway is a primary target (Stebbing et al., 2011). Tamoxifen (antiestrogen), aromatase inhibitors (estrogen withdrawal), and fulvestrant (direct targeting on the ERa receptor) are now used to treat ERa-positive breast cancer. The main challenge in treating breast cancer is endocrine hormone resistance (Ali and Coombes, 2002).

The phosphorylation of ERa is primarily engaged in endocrine resistance via regulating transcription activity and altering stability. Endocrine resistance is more common in breast cancer patients, and it results in higher levels of Era (Rani et al., 2019). Based on SiRNA screening of multiple genes and an in vivo mouse model, protein kinase enzymes were discovered to be critical targets for ERa-positive cancer. Protein kinase inhibitors might be effective treatments for ERa-positive breast cancer. Serine/threonine kinases and tyrosine kinases are the two categories of protein kinase enzymes based on sequence and structural homology (Torres-Ayuso and Brognard, 2019).

Protein kinase domains are divided into 11 conserved sub-domains that fold into a tiny N-terminal lobe and a larger C-terminal lobe. The glycine-rich loop (P-loop) in the N-terminal lobe favors ATP binding, and the activation loop (T-loop) in the C-terminal lobe regulates phosphorylation activity (Robinson et al., 2000). Lemur tyrosine kinase-3 (LMTK3) is a class of serine–threonine–tyrosine kinase (Shi et al., 2014). The LMTK3 protein may have a role in the b-catenin pathway (Shi et al., 2014) and leukemic cell survival. Exon changes in human LMTK3 are responsible for its specific vulnerability to ERa-positive breast cancer (Stebbing et al., 2012). There is no exact close LMTK3 homolog in humans; however, it shows some degree of homology to the epidermal growth factor receptor (EGFR). Nevertheless, both proteins differ in detail in terms of the mechanisms by which they achieve an inactive state and prevent ATP binding. Different variants of LMTK3 are reported (VAR_057116 and VAR_028943) (Ditsiou et al., 2020). Positive selection of genes identifies LMTK3 as a novel potential biomarker for ERa-positive breast cancer (Anbarasu and Jayanthi, 2018).

In breast cancer cells, LMTK3 regulates ERa via phosphorylation activity and is directly implicated in the modulation of endocrine resistance. Furthermore, LMTK3 protects ERa from proteasomal degradation. The amount of LMTK3 is overexpressed in aggressive breast cancer and closely correlates with survival and responsiveness to hormonal therapy. Ki6713, a tumor proliferation marker, correlates with LMTK3 expression (Zhao et al., 2013). ERa is highly expressed in breast cancer, and LMTK3 phosphorylation of ERa promotes breast cancer cell proliferation, angiogenesis, migration, and progression (Xu et al., 2014).

LMTK3 inhibition in breast tumors affected ERa at two levels: mRNA production and protein stability. Targeting LMTK3 in ERa-positive breast tumors is thought to be more successful than downregulating ERa mRNA expression in cancer cells (Johnson and O’Malley, 2011). Lemur tyrosine kinase-3 (LMTK3; also known as LMR3, TYKLM3, and KIAA 1883) is a predicted dual-specificity protein kinase whose expression levels have been implicated in cancer cell invasion, endocrine resistance, poor prognosis, and overall tumor progression in different types of malignancies (Giamas et al., 2011; Zhao et al., 2013).

Computational methods have paved a new way for different approaches to finding suitable anticancer targets (Katsila et al., 2016). These approaches provide the basis for computer-aided drug design and allow having an insight into the dynamics of a suitable target (Barabási et al., 2011). Results generated from binding underwent further screening in which only a set of non-anticancer drugs were chosen. This will also enable experimentalists to directly experimentally validate the resulting compounds without requiring further validation assays. Molecular dynamics simulation was performed to verify the time-dependent behavior of the finally selected complexes and their pattern matching with the control inhibitor. This has become a widely used discipline by computational scientists and provides new insights into the time-bound behavior of biomolecules and creates measurably different simulation times, allowing to differentiate signals from different parts of a molecule (Leach, 2007) and using the MM(GB/PB)SA method to calculate binding free energies of a complex (Miller et al., 2012). This will save the cost of applying all experimental assays in repetitions and will strengthen the fact that such approaches can be applied in the future for identifying other therapeutic agents. Thus, the finding observed from the current work is a way forward in deciding the rational development of adjuvant molecules suppressing tumor formation and, in turn, addressing the curse of breast cancer.

The methodology depicted in Figure 1 provides a concise summary of the process for prioritizing potential drug targets. These steps enabled us to select the current target for computer-aided drug design (CADD) analysis.

Figure 1. Flow chart of the current study.

Herein, we use an anticancer library. The library was imported to PyRx 0.8 software (Dallakyan and Olson, 2015), where the library was filtered using the Lipinski rule of five. Afterward, the library was energy-minimized using an MM2 force field using ChemDraw 3D (Cousins, 2011). The energy-minimized structures were then converted into the pdbqt format to make them ready for docking.

The anticancer library was explored to identify potential target inhibitors. In the process of structure-based virtual screening, the library, which included inhibitors derived from natural sources, was utilized https://www.selleckchem.com/screening/fda-approved-anticancer-drug-library.html. Natural compounds, renowned for their diverse chemical scaffolds, offer a promising avenue for the discovery of novel anticancer agents with potentially reduced adverse effects. The choice to concentrate on natural inhibitors was made to harness the inherent diversity of natural products, aligning with the specific objectives of our investigation. Approximately 5,000 compounds were retrieved and were filtered against drug- and lead-likeness. Prior to that, the library was screened for drug-like inhibitors first and then lead-like inhibitors. This was accomplished using Ligand Scout 4.1’s (Wolber and Langer, 2005) drug-like and lead-like rules, which included the following drug-like rules: the five rules of Lipinski (MW ≤ 500, HBA ≤10, MLogP ≤4.15, HBD ≤5, and TPSA 40–130 Å2) (Lipinski, 2004), Veber filter (rotatable bonds ≤10, TPSA ≤140) (Veber et al., 2002), Ghose filter (LogP ≥0.4–≤ 5.6, MW ≥ 160–≤ 480, Atoms ≥20–≤ 70, and MR ≥ 40–≤ 130) (Ghose et al., 1999), Egan rule (WLogP ≤5.88 and TPSA ≤131.6) (Egan et al., 2000), and Muegge rule (TPSA ≤150, numbers of ring ≤7, numbers of carbon >4, number of heteratoms >1, HBA ≤10, MW ≥ 200–≤ 600, number of rotatable bonds ≤15, HBD ≤5, and XLogP ≥ −2–≤ 5). Filters for lead-likeness (250 MW 350, XLOGP 3.5, and rotatable bonds 7) were used to further examine drug-like compounds, as shown in Tables 1, 2. The protocol outlined in this section adheres to a widely accepted and standardized approach utilized in the screening and filtering of lead-like inhibitors against the specific protein target. This approach follows established practices in the field, ensuring consistency and comparability with other studies in the domain of lead compound discovery (Oprea, 2002).

Table 1. Physicochemical properties of the top complexes along with the control inhibitor.

Table 2. ADMET properties of the selected small drug-like molecules by SwissADME.

Drug targets were selected based on their mechanism of function and involvement in cancer, especially in breast cancer. The best targets from the Protein Data Bank were retrieved with PDB ID: 6seq. In this study, only the available LMTK3 kinase domain was used in the structure-based virtual screening and biophysics analysis. The chosen target protein models were subjected to energy reduction to improve their quality. Herein, we did not explicitly consider post-translational modifications or delve into the broader aspects of LMTK3’s interactions with other proteins. Our study is driven by the specific research goals of identifying and characterizing potent LMTK3 inhibitors. The crystallographic experiment showed that the LMTK3 structure derived from the database exists in a monomeric form under certain circumstances. Our investigation focused on the monomeric form of LMTK3 as shown in the crystal structure, despite previous literature indicating possible dimerization under certain circumstances. Furthermore, the visualization tool UCSF Chimera (Chen et al., 2014) was used to carry out this procedure, and the structures were minimized using Gasteiger charges. Structural restrictions were removed using 1,500 rounds of minimization runs, with 750 steps using the steepest descent approach and 750 steps using the conjugate gradient method (Ahmad and Azam, 2020b). The ff03.rl force field was applied with a step size of 0.02. The minimized proteins were then submitted to a validation process to determine their quality before being used in docking investigations. The adjustments made to the cut-off values were not arbitrary or based on ad hoc thumb rules; rather, they are all adjusted by a comprehensive review of literature rules and further refined through optimization (Bashir et al., 2024).

We performed virtual screening, a computer method used to identify possible lead compounds at an early stage. This technique included swiftly evaluating a large collection of anticancer drugs via computer-based algorithms and simulations. Virtual screening of the ligand library against LMTK3 was performed using PyRx (Dallakyan and Olson, 2015). To account for protein flexibility, especially in regions with disorder, we leveraged various sources of structural information. Experimental crystal structures, obtained from reliable databases https://www.rcsb.org/structure/6SEQ, played a pivotal role in providing insights into the potential conformations of the protein, encompassing both ordered and disordered regions (Stebbing et al., 2018). This approach allowed us to consider the inherent flexibility of the protein and enhance the realism of our docking simulations. The docking coordinates set depend on the active site information obtained from the active site prediction phase. The coordinates of an active site residue Asn313 x = −19.236 y = 23.590 and z = −56.880F were applied during the docking study for each compound. The number of iterations carried out is 100. Each docking conformation was assigned with a binding energy score in kcal/mol. The one with the lowest energy score was considered the best binder and complexed with the enzyme for further investigation. The visualization of the best docked complexes was carried out using UCSF Chimera v1.15 and Discovery Studio v.2020 (Biovia, 2017).

The starting conformation for the molecular dynamics (MD) simulations was based on the coordinates from the Protein Data Bank (PDB) and was prepared for simulation using the AMBER16 simulation program (Case et al., 2016). The initial parameterization was carried out using the antechamber program of AMBER. The force field that was used is FF14SB (Case et al., 2014) for the enzyme, while for ligands, GAFF was used (Sprenger et al., 2015). The simulation was performed in the following order: energy minimization of the complexes, heating, equilibrium, and production run of at least 300 ns. Temperature control was achieved via Langevin dynamics. SHAKE algorithm was applied for constraining hydrogen bonds (Kräutler et al., 2001). Simulation trajectories were evaluated through the CPPTRAJ program (Roe et al., 2013). The coordinates of alpha carbon (C) are commonly thought to indicate an amino acid’s location in three-dimensional space. RMSD is a metric that allows comparing the relative locations of protein C atoms by computing their averaged distance over a period. It is written mathematically as

where N is the number of compared atoms and di is the distance between the ith pair of atoms.

Another key metric of structural changes is the RMSF. It is used to determine the backbone atoms of the docked target (N, C, and C). It is the root mean square of the averaged distance between an atom and its average geometric location in a particular set of dynamics, and it may be read as the set of atom positions recorded over a specific time scale. The RMSF is calculated using the following equation:

where T represents the time interval, xi represents the position of an atom at a particular time, and x represents the averaged position of the atom.

The radius of gyration is used to assess the overall packing quality and density of a structure. It is a physical characteristic that may be estimated experimentally, most commonly via small-angle X-ray scattering (SAXA). The following equation was used to quantify the compactness of a macromolecular system:

where N is the total number of atoms, mi denotes the mass of atom I, ri denotes the position vector of atom I, and rcm denotes the molecule’s center of mass.

The amount of solvent-accessible surface area (SASA) is critical in determining the conformation and functionalities of biological macromolecules. Typically, the amino acid residues on a protein’s surface operate as active sites and/or interact with other molecules and ligands. This allows researchers to obtain a better understanding of the molecule’s behavior in a solvent environment, distinguishing between its hydrophilic and hydrophobic properties and revealing the components involved in protein–ligand interactions (Gromiha and Ahmad, 2005).

WaterSwap was carried out to calculate the binding energies of the complex along with solvated systems inside the cavity. WaterSwap samples conformational space using MC simulations [using the Sire program (Woods et al., 2014)]. To begin the MC computation, a starting configuration is necessary, which is commonly obtained via MD simulations of the protein–ligand complex (Ahmad and Azam, 2020a). WaterSwap computations were used to construct the most representative structures from the ensemble obtained by the MD simulations. The free energy was then calculated for all three complex systems along with the control system using the free energy perturbation (FEP) (Zwanzig, 2004), thermodynamic integration (TI) (Oostenbrink et al., 2000), and Bennett acceptance ratio (BAR) (Bennett, 1976) methods. Three approaches were chosen to calculate free energy because the similarity of the energy values obtained from these methods indicates the dependability and convergence of the results.

The top compounds from the library were screened and further tested for ADMET characteristics. All the compounds in each library had already been subjected to Lipinski’s rule of five and ADMET descriptions prior to docking. If the medicine meets the physiochemical property requirements of Lipinski’s rule of five, it will be declared optimum. This test assesses if a chemical compound is drug-like, defined as having biological action that may be administered orally. A drug-like chemical molecule should contain hydrogen bond donors (HBDs) > 5, molecular weight (MW) 500 g/mol, and hydrogen bond acceptors (HBAs) sites >10, and log p-value 5 reflects hydrophobicity of a compound, according to the rule of thumb (RO5). Furthermore, two new factors were added: rotatable bonds 10 and a polar surface area (PSA) 140 A, both of which are concurrent with drug flexibility and permeability. SwissADME (Daina et al., 2017) was used to determine ADMET descriptors to compute absorption, distribution, metabolism, elimination, and toxicity (ADMET). The ADMET blood–brain barrier was used in this study. The ADMET plasma protein binding descriptor was used to determine if a medication molecule will form a strong link with the blood carrier protein. The ADMET CYP2D6 binding model was used to evaluate cytochrome P450 2D6 enzyme inhibition using a compound’s 2D chemical structure. ADMET hepatotoxicity assesses the likelihood of hepatotoxicity proliferation in humans caused by a wide range of structurally varied substances. Compounds should correspond with the needs of AlogP98 5 and PSA 140 2 for good cell permeability. Tables 1, 2 provide the Lipinski and ADMET descriptions of the top eight complexes along with the control molecule.

The selected ligand molecules were docked into the active site of the target lemur tyrosine kinase-3 (LMTK3) using PyRx. Around 4,000 compounds were utilized in the docking study. The active site of the target was supported by the literature where the coordinates x = −19.236, y = 23.590, and y = −56.880 of the Asp313 were applied. After docking, the top inhibitor molecules were selected based on their binding affinities (Table 3). The top1 compound (R)-5,7-dihydroxy-2-(4-hydroxy-3-methoxyphenyl)-6,8-bis(3-methylbut-2-en-1-yl)chroman-4-one results in a high binding affinity of −9.8 kcal/mol. The top2 compound 1,3,6,7-tetrahydroxy-2,5,8-tris(3-methylbut-2-en-1-yl)-9H-xanthen-9-one has a binding affinity of −9.5 kcal/mol, and the top3 compound (R)-5,7-dihydroxy-2-(4-hydroxyphenyl)-6-(2-methylbut-3-en-2-yl)chroman-4-one has a binding affinity of −8.8 kcal/mol, whereas the control inhibitor molecule 6-(4-pyridinyl)-3-[3-(trifluoromethoxy)phenyl]imidazo [1,2-b]pyridazine results in a −8.23 kcal/mol binding affinity. Detailed visualization analysis carried out through UCSF Chimera and the preferred orientation of the ligand binding, as shown in Figure 2. The 2D interactions of the complexes were found via the Discovery Studio visualizer, which results in hydrogen bonding and other interactive linkages, as shown in Figure 3. The docked poses show the active residues of the protein-binding site. These hotspot residues, which are the active site domain residues, include Tyr163, Asn138, Asp133, Tyr56, Glu52, Ser132, Asp313, and Asp151, making strong hydrogen bonds during the docking simulation. Furthermore, validation of the hotspot residues has been investigated through molecular dynamics simulation.

Table 3. Docked score of the top complexes along with the control inhibitor against the target protein.

Figure 2. Docked poses of the top complexes along with the control inhibitor against the target protein.

Figure 3. Depicting the top complexes along with the control inhibitor against the target protein.

Structural convergence words include RMSD, RMSF, and Rg. The examination of RMSD with backbone atoms demonstrated early equilibration up to 5 ns, and after 6 ns, the structure began to converge and achieved the stable state, generating configuration with an RMSD of 4.5 MD simulation. The RMSD plot demonstrated that LMTK3 was in a stable conformation with minimal deviation from its original structure. The Cα atom changes were discovered by examining the RMSF data in connection to the LMTK3 residues. The RMSF proved that certain variations occurred in the protein’s loop regions with one exception which is flexible throughout the simulation time intervals. The root mean square deviation (RMSD) analysis was used to evaluate the docked protein’s backbone dynamics and Cα atom motions throughout a 200 ns simulation time interval. Some variations were noticed within the first 15 ns of the simulation, but the subsequent portion of the simulation indicated general stability. As illustrated in Figure 4A, the average RMSD value for the Top1 docked protein was estimated to be 4.4 Å, with a maximum peak at 5.85 Å. According to the full RMSD graph, no substantial domain alterations occurred within the physical setting of the protein–ligand combination during the top1 complex simulated system. We also investigated the top2 complex system, which remains more stable than top1 at the start of the simulation but remains a bit fluctuated during the mid-time intervals, as shown in the black line in Figure 4A. The average RMSD was recorded as 3.8 Å with a maximum value of 4.3 Å at 145 ns but remained stable till the end of the simulation. The third complex shows stability in protein at 4.5 Å, and the ligand became stable at 48 ns. The ligand shows different contacts with different residues at the start of the simulation and then it remains static in its binding cavity with some angular displacement at the end of the simulation. Rg values, which are inversely linked to protein compactness, were critical in determining ligand stability. Stable protein folding is characterized by little variation in Rg values, indicating static ligand behavior within the cavity. The rGyr values of compound top1 within the binding pocket varied from 36.8 Å to 4.0 from 0 to 200 ns and then fluctuated between 40.1 Å and 38.5 Å for all the top complexes throughout 200 ns (Figure 4B). Compound 2 had rGyr values ranging from 35.6 Å to 39Å and stabilized at 39.2 Å (Figure 4B). Compound 3, on the other hand, exhibited considerable oscillations, with rGyr values ranging from 37.7 Å to 38.6 Å throughout the simulation, indicating stable behavior (Figure 4B). On other hand, the control inhibitor shows considerable fluctuation throughout the simulation time intervals, inferring an average rGyr value of 45.6 Å with a maximum value of 49.3 Å. Considering these complexes in comparison with the control inhibitor, it has been inferred that our selected top complexes are more stable during the simulation time intervals. As shown in the figures, the control inhibitor shows more fluctuation, and the ligand inside the cavity moves inside the binding cavity with more structural changes, as seen in the RMSF graph.

Figure 4. Depiction of the stability of the bound complexes, i.e., top1, top2, and top3 LMTK3 proteins in terms of (A) root mean square deviation (RMSD), (B) radius of gyration (C), and root mean square fluctuation (RMSF). (D) SASA analysis. The structural analysis of MD simulations was performed using trajectory files, including time series values.

The proteins’ local changes were analyzed using RMSF, which shows an average RMSF of 1.81 Å, 2.38 Å, and 2.35 Å for complex-1, complex-2, and complex-3 for 200 ns, respectively (Figure 4C), whereas high variations were found in the N-terminal and C-terminal areas, as well as in the loop regions. The interacting binding site residues at the loop region at the C-terminus fluctuate, resulting in an open binding pocket conformation. The assessment of all the complexes at different ns revealed that complex-1 deviates by 1.4 Å from the docked complex while preserving the ligand in its binding pocket. Complexes 2 and 3 have deviations of 2.1 Å and 2.2 Å, respectively; however, this deviation and opening of the binding pocket did not allow the ligand to be retained outside the surface. However, this change has been observed in the control inhibitor, which deviates more as compared to the top complexes; herein, the control shows more deviation with a maximum RMSF value of 5 Å, and the ligand moves around the surface of protein-binding site residues, as shown in Figure 4C. Furthermore, as shown in Figure 4, the ligand remained well positioned inside the binding site throughout the simulation and did not cause any protein destabilization among all the top complexes.

A SASA analysis was performed to investigate the behavior of LMTK3’s hydrophilic and hydrophobic SASA. SASA findings of 5,000–5,500 nm indicated that accessibility was maintained with few modifications during the simulation and validated that the residues had been thoroughly exposed to the solvent (see Figure 4D). The findings revealed that in loop1, during the simulations, the area (182–202 residues) demonstrated less flexibility, with modest differences in a-helix and turn occurrence. The loop2 area (212–232 residues) showed more adaptability during all the complexes, but these structural modifications have been observed more fluctuated during the control complex system with a SASA value of more than 12,500 nm.

Binding energy analysis provides insight into individual residue contributions, which aids in the exploration of the protein–ligand complex’s composition. Binding energies for top1, top2, top3, and control complexes estimated by MM-GBSA are −33.74 kcal/mol, −44.07 kcal/mol, - 42.35 kcal/mol and - 26.95 kcal/mol, respectively, whereas the MM-PBSA approach projected binding energies of −33.47 kcal/mol for the top1 complex, −44.98 kcal/mol for the top2 complex, −43.09 kcal/mol for the top3 complex, and −27.41 for the control inhibitor complex, respectively. According to MM-GBSA, the van der Waals score observed (−33.63 kcal/mol) and electrostatic interactions (−12.14 kcal/mol) contribute the most to the binding energy in the top1 complex. High binding scores for the other top complexes were also investigated along the control inhibitor complex. Furthermore, gas phase and solvation phase interactions have contributed to the binding energy of all the complexes, according to MM-GB/PBSA (Table 4).

Table 4. Binding energy calculation for the top complexes along the control inhibitor after simulation time intervals of 200 ns.

WaterSwap analysis was used to determine the absolute binding free energy of each complex, including the control inhibitor molecule. The computed binding free energy values were shown in the WaterSwap analysis (in kcal/mol) for top inhibitors in relation to critical hotspot residues in the complex system. The computations were carried out for 1000 iterations with the cluster cut-off size set to 1 and a maximum of 25 members allowed in each cluster. A total of 1.6 109 MC motions were conducted in order to calculate the Bennett value, thermodynamic integration (TI), free energy perturbation (FEP), and TI quadrature, which are some of the approaches used. Figure 5 depicts the calculated free energy values for each complex. When the preceding procedures converge with a difference of 1 kcal/mol, the final value obtained was regarded as excellent. The drug complexes, including top1, top2, top3, and control, have an FEP value of −47.1 kcal/mol, −48.53 kcal/mol, −51.8 kcal/mol, and −27.14 kcal/mol, respectively, whereas a Ti value of −47.59 kcal/mol for top1, −48.53 kcal/mol for top2, −51.36 for top3, and −26.85 for the control system was recorded. Another approach used to find out the binding energies calculation was a BAR value for all the complexes; herein, −46.87 kcal/mol for top1, −49.68 kcal/mol for top2, −51 kcal/mol for top3, and −27.55 kcal/mol for the control system were calculated. Hence, these investigations suggest that all the top inhibitor molecules show the best binding free energies in terms of the control inhibitor system. These results further evidence and predict the scale of the potency of inhibitors for the experimental studies in a real-time system.

Figure 5. WaterSwap analysis of the simulated complex systems along with the control system.

Human lemur tyrosine kinase-3 (LMTK3) is a dual-specificity kinase with a well-established oncogenic involvement in several tumor types, indicating its potential as a therapeutic target. In several breast cancer cohorts, LMTK3 expression is markedly higher in high-grade breast tumors and is linked with poor survival rates (Giamas et al., 2011; Stebbing et al., 2012). Given its function in endocrine resistance, more research into its potential usefulness as a therapeutic target is necessary. In breast cancer, neoadjuvant and adjuvant chemotherapy has enhanced overall survival (Gianni et al., 2014; Hatzis et al., 2016). Resistance to cytotoxic chemotherapy, on the other hand, is the leading cause of treatment failure and mortality in women with breast cancer (Yu et al., 2013; von Minckwitz et al., 2016). An abundance of LMTK3 in human cancer was related to disease-free survival and suggested a good response to endocrine treatment. The findings also suggested that targeting LMTK3 might be useful in the treatment of drug-resistant malignancies. Kinases are phosphorylation-catalyzing enzymes that play an important part in cell activity. Kinomes are now the focus of large-scale genomics and RNA interference (RNAi) screening in drug discovery efforts, particularly in the quest for anticancer medicines (Collins and Workman, 2006).

LMTK3 deletion resulted in an 80% decrease in ERa protein in MCF-7 cells. ERa breakdown is proteasome-dependent, and inhibiting its degradation raises ERa protein levels. LMTK3 deletion decreases ERa half-life and promotes ubiquitination. LMTK overexpression, on the other hand, phosphorylates ERa, protecting it from proteasomal destruction and, therefore, stabilizing the ERa protein. LMTK3 knockdown decreased estrogen receptor 1 (ESR1) mRNA expression (Cilibrasi et al., 2021). We have used in silico drug design strategies to discover new inhibitors against the LMTK3 target to stop the mechanistic behavior in the cancer pathway. This is carried out via a molecular docking strategy, where we obtained top complexes with high binding affinities with strong hydrogen bond interactions. In addition, a control inhibitor was also docked against the protein LMTK3, showing strong interactions. These results were then further validated via molecular dynamics simulation studies which inferred the complexity of the top complexes, with mechanistic time interval behavior.

The comprehensive investigation of complex-1, complex-2, and complex-3 found that complex-2 is more stable because it mimics the natural secondary structural elements (SSEs). The SSE composition for each trajectory frame across the simulation duration demonstrates that all three top complexes are structurally stable.

The mean RMSD for the top1 was 4.4 Å, with a maximum value of 5.85 Å. We also looked at top2 complex systems, which are more stable than top1 at the start of the simulation but are more fluctuated over the mid-time intervals, as shown by the black line in Figure 4A. The average RMSD was 3.8, with a maximum of 4.3 at 145 ns, although it remained constant until the end of the experiment. The third compound demonstrates protein stability at 4.5 and ligand stability at 48 ns. The ligand makes distinct interactions with different residues at the start of the simulation and then remains static in its binding cavity with some angular movement at the conclusion. Rg values were crucial in assessing ligand stability since they are inversely related to protein compactness. The RMSF value also determined the thermal flexibility of the system during the simulation time intervals, which inferred that the top1 complex was more stable as compared to the control system. The other systems were also stable during simulations with a bit of fluctuation at some time intervals but remained static at the end of the simulation. Following the trajectory analysis, the top two compounds were subjected to atom contact parallel to SASA to determine non-native atom contacts as well as ligand–solvent interactions. Furthermore, a unique and more advanced WaterSwap approach was employed to cross-verify the MMGB/PBSA results. WaterSwap agrees on the inhibitor affinity for the receiver domain, as evidenced by the minimal convergence of the Bennett, TI, and FEP algorithms.

These investigations further open a way toward pilot-scale experimental analysis for further validation of compounds against cancer therapies.

Ahmad, F., and Azam, S. S. (2020a). From pan-genome to protein dynamics: a computational hierarchical quest to identify drug target in multi-drug resistant Burkholderia cepacia. J. Mol. Liq. 317, 113904. doi:10.1016/j.molliq.2020.113904

CrossRef Full Text | Google Scholar

Ahmad, F., and Azam, S. S. (2020b). Role of ring positioning and preferential occupation of ligand obtained through molecular dynamics simulation of peptidoglycan associated lipoprotein (Pal). J. Mol. Graph. Model. 98, 107585. doi:10.1016/j.jmgm.2020.107585

PubMed Abstract | CrossRef Full Text | Google Scholar

Anbarasu, K., and Jayanthi, S. (2018). Identification of curcumin derivatives as human LMTK3 inhibitors for breast cancer: a docking, dynamics, and MM/PBSA approach. 3 Biotech. 8, 228. doi:10.1007/s13205-018-1239-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bashir, Y., Noor, F., Ahmad, S., Tariq, M. H., Qasim, M., Tahir ul Qamar, M., et al. (2024). Integrated virtual screening and molecular dynamics simulation approaches revealed potential natural inhibitors for DNMT1 as therapeutic solution for triple negative breast cancer. J. Biomol. Struct. Dyn. 42, 1099–1109. doi:10.1080/07391102.2023.2198017

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennett, C. H. (1976). Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys. 22, 245–268. doi:10.1016/0021-9991(76)90078-4

CrossRef Full Text | Google Scholar

Biovia, D. S. (2017). Discovery studio visualizer. San Diego, CA: USA, 936.

Google Scholar

Case, D. A., Babin, V., Berryman, J. T., Betz, R. M., Cai, Q., Cerutti, D. S., et al. (2014). The FF14SB force field. Amber 14, 29–31.

Google Scholar

Case, D. A., Betz, R. M., Cerutti, D. S., Cheatham, T., Darden, T., Duke, R. E., et al. (2016). Amber 16. San Francisco: University of California. doi:10.13140/rg.2.2.27958.70729

CrossRef Full Text | Google Scholar

Chen, J. E., Huang, C. C., and Ferrin, T. E. (2014). RRDistMaps: a UCSF Chimera tool for viewing and comparing protein distance maps. Bioinformatics 31, 1484–1486. doi:10.1093/bioinformatics/btu841

PubMed Abstract | CrossRef Full Text | Google Scholar

Cilibrasi, C., Ditsiou, A., Papakyriakou, A., Mavridis, G., Eravci, M., Stebbing, J., et al. (2021). LMTK3 inhibition affects microtubule stability. Mol. Cancer 20, 53. doi:10.1186/s12943-021-01345-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Cousins, K. R. (2011). Computer review of ChemDraw ultra 12.0. ACS Publications.

Google Scholar

Daina, A., Michielin, O., and Zoete, V. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 7, 42717. doi:10.1038/srep42717

PubMed Abstract | CrossRef Full Text | Google Scholar

Ditsiou, A., Cilibrasi, C., Simigdala, N., Papakyriakou, A., Milton-Harris, L., Vella, V., et al. (2020). The structure-function relationship of oncogenic LMTK3. Sci. Adv. 6 (46), eabc3099. doi:10.1126/sciadv.abc3099

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghose, A. K., Viswanadhan, V. N., and Wendoloski, J. J. (1999). A knowledge-based approach in designing combinatorial or medicinal chemistry libraries for drug discovery. 1. A qualitative and quantitative characterization of known drug databases. J. Comb. Chem. 1, 55–68. doi:10.1021/cc9800071

PubMed Abstract | CrossRef Full Text | Google Scholar

Giamas, G., Filipović, A., Jacob, J., Messier, W., Zhang, H., Yang, D., et al. (2011). Kinome screening for regulators of the estrogen receptor identifies LMTK3 as a new therapeutic target in breast cancer. Nat. Med. 17, 715–719. doi:10.1038/nm.2351

PubMed Abstract | CrossRef Full Text | Google Scholar

Gianni, L., Eiermann, W., Semiglazov, V., Lluch, A., Tjulandin, S., Zambetti, M., et al. (2014). Neoadjuvant and adjuvant trastuzumab in patients with HER2-positive locally advanced breast cancer (NOAH): follow-up of a randomised controlled superiority trial with a parallel HER2-negative cohort. Lancet Oncol. 15, 640–647. doi:10.1016/S1470-2045(14)70080-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Gromiha, M., and Ahmad, S. (2005). Role of solvent accessibility in structure based drug design. Curr. Comput. - Aided Drug Des. 1, 223–235. doi:10.2174/1573409054367664

CrossRef Full Text | Google Scholar

Hatzis, C., Symmans, W. F., Zhang, Y., Gould, R. E., Moulder, S. L., Hunt, K. K., et al. (2016). Relationship between complete pathologic response to neoadjuvant chemotherapy and survival in triple-negative breast cancer. Clin. Cancer Res. 22, 26–33. doi:10.1158/1078-0432.CCR-14-3304

PubMed Abstract | CrossRef Full Text | Google Scholar

Katsila, T., Spyroulias, G. A., Patrinos, G. P., and Matsoukas, M. T. (2016). Computational approaches in target identification and drug discovery. Comput. Struct. Biotechnol. J. 14, 177–184. doi:10.1016/j.csbj.2016.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kräutler, V., Gunsteren, W. F. V., and Hünenberger, P. H. (2001). A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 22, 501–508. doi:10.1002/1096-987X(20010415)22:5<501::AID-JCC1021>3.0.CO;2-V

CrossRef Full Text | Google Scholar

Labrie, F., Labrie, C., Bélanger, A., Simard, J., Gauthier, S., Luu-The, V., et al. (1999). EM-652 (SCH 57068), a third generation SERM acting as pure antiestrogen in the mammary gland and endometrium. J. Steroid Biochem. Mol. Biol. 69, 51–84. doi:10.1016/s0960-0760(99)00065-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Leach, A. (2007). Ligand-based approaches: core molecular modeling.

Google Scholar

Miller, B. R., McGee, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012). MMPBSA.py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput. 8, 3314–3321. doi:10.1021/ct300418h

PubMed Abstract | CrossRef Full Text | Google Scholar

Oostenbrink, B. C., Pitera, J. W., van Lipzig, M. M., Meerman, J. H., and van Gunsteren, W. F. (2000). Simulations of the estrogen receptor ligand-binding domain: affinity of natural ligands and xenoestrogens. J. Med. Chem. 43, 4594–4605. doi:10.1021/jm001045d

PubMed Abstract | CrossRef Full Text | Google Scholar

Rani, A., Stebbing, J., Giamas, G., and Murphy, J. (2019). Endocrine resistance in hormone receptor positive breast cancer–from mechanism to therapy. Front. Endocrinol. 10, 245. doi:10.3389/fendo.2019.00245

CrossRef Full Text | Google Scholar

Shi, H., Li, Q., Ji, M., Wu, J., Li, Z., Zheng, X., et al. (2014). Lemur tyrosine kinase-3 is a significant prognostic marker for patients with colorectal cancer. Int. J. Clin. Exp. Pathol. 7, 1101–1107.

PubMed Abstract | Google Scholar

Sprenger, K. G., Jaeger, V. W., and Pfaendtner, J. (2015). The general AMBER force field (GAFF) can accurately predict thermodynamic and transport properties of many ionic liquids. J. Phys. Chem. B 119, 5882–5895. doi:10.1021/acs.jpcb.5b00689

PubMed Abstract | CrossRef Full Text | Google Scholar

Stebbing, J., Filipovic, A., Ellis, I. O., Green, A. R., D’Silva, T. R., Lenz, H.-J., et al. (2012). LMTK3 expression in breast cancer: association with tumor phenotype and clinical outcome. Breast Cancer Res. Treat. 132, 537–544. doi:10.1007/s10549-011-1622-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Stebbing, J., Shah, K., Lit, L. C., Gagliano, T., Ditsiou, A., Wang, T., et al. (2018). LMTK3 confers chemo-resistance in breast cancer. Oncogene 37, 3113–3130. doi:10.1038/s41388-018-0197-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Veber, D. F., Johnson, S. R., Cheng, H. Y., Smith, B. R., Ward, K. W., and Kopple, K. D. (2002). Molecular propert