Machine learning based QSAR and Molecular Dynamics simulations in the structural design and mechanism of action of imidazole derivatives with anti-melanoma activity

Abstract

Cancer is defined as a group of diseases in which abnormal cells multiply and can invade other organs, requiring continuous studies for new drugs. A series of 177 imidazo[1,2-a]pyridine and imidazo[1,2-a] pyrazine synthetic derivatives were previously obtained, and their anti-melanoma IC₅₀ values have been determined. Here, Artificial Intelligence algorithms were used to select molecular descriptors and build a QSAR model, highlighting structural characteristics related to enhanced molecular potency. Additionally, the imidazopyrazine nucleus was compared to a known inhibitor of the Aurora Kinase enzyme, an important target in cancer therapy. Thus, strategic imidazopyrazines were subjected to comparative molecular dynamics calculations, providing inferences about their possible mechanisms of action. The QSAR model allows for the design and prediction of nine new analogues with favourable predicted IC₅₀ values. Molecular dynamics simulations and the estimated binding energies are consistent with the ranking of activities presented by representatives of the series.

Keywords:
Imidazole derivatives; Anti-melanoma activity; Bio-inspired algorithms; Random Forest; Molecular Dynamics; QSAR

INTRODUCTION

The discovery and development of new drugs is a very complex, expensive, and lengthy process, typically costing an estimated 2.6 billion USD and taking 12 years on average (Chan et al., 2019). Today, due to advances in technology, in silico methodologies are valuable and well-established in modern research. The most important techniques in the field of computer-aided drug development (CADD) can be divided into Structure-Based and Ligand-Based Drug Design (SBDD and LBDD). The first is represented mainly by Molecular Docking and Molecular Dynamics simulations, while the latter is generally represented by Pharmacophore search and Quantitative Structure-Activity Relationship models (Shaker et al., 2021).

Considering the computational medicinal chemistry scenario, machine learning (ML) techniques have been used in recent decades for the drug discovery process, and studies in this area have increased substantially, due in part to a phenomenon known as the “Big Data Revolution” (Klanbauer, Hochreiter, Rarey, 2019; Talevi et al., 2020). With the increased availability and decreasing costs of technologies to generate, store, and analyse large and diverse datasets, the field of ML in the drug discovery process to facilitate decision-making has grown exponentially (Sellwood et al., 2018). Machine learning is a subarea of Artificial Intelligence (AI) that aims to develop and apply computer algorithms that learn from raw, unprocessed data to later perform a specific task (Carracedo-Reboredo et al., 2021). When compared to other approaches, ML-based techniques can be easily scaled to big datasets without the need for extensive computational resources, which makes this technique one of the most important and rapidly evolving topics in computer-aided drug discovery (Klanbauer, Hochreiter, Rarey, 2019; Lo, Torng, Altman, 2018). The integration of ML aims to improve the accuracy and efficiency of the methods in CADD (Vamathevan et al., 2019). Building machine learning models from high-dimensional data in small samples is a challenge. One common approach to overcome this issue is applying feature selection, an essential step in data preprocessing, which aims to find a feature subset for effective supervised learning (Alan, Karla, Nicola, 2015). Several approaches to feature selection are found in the literature, one of the most popular being the wrapper method, which uses an ML algorithm to evaluate each feature subset and a search algorithm to select new subsets (Bolón-Canedo, Sánchez- Maroño, Alonso-Betanzos, 2013a).

The search, development, and discovery of drugs against cancer are very important, as novel cancer therapeutics comprise a large proportion of failed drug trials, resulting in very limited new treatment options for large cohorts of patients, many of whom are refractory to current treatments (Adams, Rashidieh, 2020; Ogilvie et al., 2017; Mak, Pichika, 2019). Cancer is broadly defined as a group of diseases in which abnormal cells multiply and invade other organs (WHO, 2022). It is a complex disease, in which cells undergo a series of changes in their genes, developing characteristics such as the ability to maintain themselves in a chronic proliferative state, evade cell growth suppressors and the immune system, induce angiogenesis, alter cell metabolism, and resist cell death (Hanahan, 2022). According to the World Health Organization, 10 million lives were lost to cancer worldwide in 2020 (WHO, 2022). Therefore, it is of paramount importance that new drugs are developed for the treatment of this disease. The main goal of cancer research is to discover the most effective method of treatment for each cancer patient, as not everyone responds equally to a specific treatment due to external factors, such as lifestyle-related habits, and internal factors, like heterogeneity of cancer cells and the immune system (Hanahan, 2022).

Melanoma is a type of skin cancer that arises from the malignant transformation of melanocytes, the cells that produce melanin (Cramer, 1991). Although melanoma accounts for only about 5% of skin cancers, it is responsible for more than 75% of deaths from the illness (Rebecca, Somasundaram, Herlyn, 2020). Historically a disease of rare occurrence, the increase in cases in the last 50 years, especially in fair-skinned populations of European ancestry, is alarming (Erdmann et al., 2013; Arnold et al., 2014). The most recent statistics indicate that in 2020 alone, approximately 325,000 cases were diagnosed, and 57,000 casualties were recorded (GCO - WHO, 2023). If this pace continues, experts project more than half a million cases per year and 96,000 deaths by 2040 (Arnold et al., 2022). Therefore, the search for new anti-melanoma agents is imperative.

Cell cycle dysregulation is one of the main causes of melanoma. Genomic studies show that there is a significant increase in the expression of Aurora-A and B proteins in melanoma cells; however, we still have no approved drugs acting on Aurora kinases for the treatment of this type of cancer (Fatma et al., 2021). Aurora kinase A (AURKA) belongs to the serine/threonine kinases family, whose activation is associated with cell division processes. Research involving this protein has been in focus in recent years, and it is considered a promising target in cancer therapy as well as in preventing resistance phenomena of tumours (Fatma et al., 2021; Du et al., 2021; Punt et al., 2021).

Molecules containing imidazole scaffolds have been studied as potential drugs in cancer therapy. Considering the anti-melanoma action, we found a series of 177 imidazo[1,2-a]pyridine and imidazo[1,2-a]pyrazine synthetic derivatives, whose IC50 values had been determined and are available in the PubChem database (Garamvölgyi et al., 2016). In this regard, establishing structure-activity relationships as well as generating a quantitative model for predicting more potent analogues may provide a useful way forward for rational planning. The aforementioned article concludes that the in vitro data obtained do not confirm the expected mechanism of action of the anti-proliferative effect via the inhibition of B-Raf protein (Garamvölgyi et al., 2016). Surprisingly, we found in the Protein Data Bank a crystallographic structure of the Aurora Kinase enzyme with a ligand that exactly holds the structural core and substitution pattern of the most active imidazole series of analogues.

In this sense, the present work aims to use Artificial Intelligence techniques in QSAR modelling of the imidazole derivatives reported by Garamvölgyi and collaborators (Garamvölgyi et al., 2016), allowing structural rationalization and proposals for potentially more active architectures. Furthermore, we present Molecular Dynamics calculations for the elucidation of the mechanism of action of the series via Aurora Kinase inhibition, comparing the results with the co- crystallographic imidazo-pyrazine analogue. We expect to open a path for the future exploration of the structural pattern presented here, with both in silico and in vitro approaches, in order to generate and test promising compounds for application in the treatment of human melanoma.

METHODS

Initial computational procedures

The compounds with anti-melanoma activity investigated in this article were obtained from the PubChem database (Kim et al., 2023), which includes a set of 117 imidazole derivatives reported by Garamvölgyi and collaborators (Tables 1-9 of the original paper) (Garamvölgyi et al., 2016), along with their respective 3D structures and IC₅₀ values. These 3D representations were computed by PubChem using a specific algorithm that optimizes the geometries at the MMFF94s force field (Kim, Bolton, Bryant, 2013; Bolton et al., 2011). Other designed molecules were built using the ACD/ Chemsketch program and optimized by the Universal Force Field (ACD/Labs, 2022).

AI algorithms and QSAR modelling

The matrix of molecular descriptors was obtained using Dragon 7 software, which can generate 5,270 molecular descriptors, including the simplest atom types, functional groups and fragment counts, topological and geometrical descriptors, three-dimensional descriptors, among others (Kode, 2017). After generation, descriptors with constant or near-constant values, as well as those with missing values, were discarded. The biological activity of the compounds, expressed by IC₅₀, was converted to -Log(IC₅₀) or “pIC₅₀”, reducing the standard deviation and making higher values correspond to better activities. A previous external test set was not defined in this study, as this procedure could compromise the final model, considering the low number of samples. The literature indicates that the validation of models using powerful Artificial Intelligence algorithms can be satisfactory by procedures such as leave-one-out and K-fold techniques (Majumdar, Basak, 2018), especially when the satisfactory model and metrics also lead to the interpretation of physical sense and practical utility in the design of molecules.

After completing the above procedures, a 117 x 3,232 matrix was created. The first column contains the pIC₅₀ dependent term (Y). This matrix was then used in the Wrapper variable selection approach, which employs a machine learning algorithm to evaluate each feature subset and a search algorithm to select new subsets. We utilized a two-step Wrapper approach in this study (Figure 1).

FIGURE 1
Flowchart of variable selection procedure used in this paper.

The first step used the Artificial Bee Colony (ABC) algorithm as a search method. The ABC algorithm is a bio-inspired global optimization algorithm proposed by Karaboga, Basturk (2007). This algorithm works with a population of solutions simultaneously that evolve over multiple cycles. The two main elements of this algorithm are: 1) the food source corresponding to a possible solution to the optimization problem and 2) the bees that guide the search. This algorithm has three types of bees: the employed, the onlooker, and the scout. The search in ABC consists of three steps: 1) Sending the employed bees to a food source and evaluating the quality; 2) Sending onlooker bees to choose food sources after obtaining information from employed bees and updating their quality; 3) Sending the scouts to the search area to discover new food sources.

In the feature (variable) selection, each food source represents a feature subset. A detailed description of the ABC search procedure can be found in Bolaji et al. (2014). In this study, the ABC configuration included twenty cycles, fifteen employed bees, fifteen onlooker bees, and one scout bee.

Step 02 (Figure 1) was based on the Best First Algorithm (BFA), which is a graph-based local search algorithm (Hart, Nilsson, Raphael, 1968). When used for feature selection in a dataset with N features, the process starts by building N models, one for each feature. The feature associated with the highest-performing model is selected to be at the graph root. In the next iteration, N-1 models will be created with two input features: the one selected in the previous iteration and another of the remaining N-1 features. Again, the feature combination that provides the best performance is selected. The procedure is repeated until a stopping criterion is met. For this study, the criterion was five expansions without improvement.

In both Wrapper approaches, the machine learning algorithm used was Random Forest, an ensemble-type machine learning algorithm proposed by Breiman (2001). The Random Forest is a set of decision trees that make up a forest. This algorithm has two parameters: N, which represents the number of trees that make up the forest, and M, which indicates the number of variables that will be used to build the trees. The Random Forest setup used in both steps was N = 50 and M = Log(2) of the number of features. The performance evaluation metric for each feature subset was the correlation coefficient.

The two-step feature selection method was used to address three issues that would arise if each approach were used independently: 1) Using Random Forest without feature selection for this dataset generates a solution with low performance because the number of features is much greater than the number of samples, which makes convergence difficult. This problem is known in the literature as the curse of dimensionality (Brown et al., 1995). 2) Using only the first step, which involves a global optimization algorithm like ABC, on a vast search space with a small sample and correlated features, makes it difficult to fine-tune the solutions due to the stochastic nature of the algorithm. Therefore, we used the first step as a pre-filter for feature selection. 3) Using only the second step, which involves a local search algorithm like BFA, results in poor-performing solutions because the algorithm tends to get stuck in local optima early in the search. Thus, we used the second step after an initial selection with a global search algorithm.

The final model was built using the twelve features selected after the two steps and a Random Forest with one hundred trees and four features per tree. This model was superior to the three isolated approaches.

Validation parameters calculated for the final model were based on: Leave-One-Out cross-validation (coefficient of determination “Q² _LOO”); K-Fold cross- validation, considering K=5 and 50 runs with the respective coefficient of determination “Q² _KFOLD”. These metrics are very relevant considering non-linear machine learning techniques applied to QSAR (Majumdar, Basak, 2018).

Finally, it is important to point out that five outliers were detected in an initial inspection, which correspond, in a more accurate analysis, to compounds structurally similar to others but with very different activity values (14c, 15f, 15i, 15d, 6i in the article by Garamvölgyi et al., 2016). In this scenario, these samples can be considered true outliers, whose explanation may even be due to inaccuracies in the test. Therefore, it was decided to remove them so as not to damage the model.

Molecular Dynamics simulations

Molecular Dynamics (MD) calculations were performed using the Gromacs 2019.2 package (Abraham et al., 2015) and the GROMOS 54a7 force field. The protein used in simulations was downloaded from the RCSB Protein Data Bank (www.rcsb.org) (Berman et al., 2000) under the code PDB ID: 3NRM (Belanger et al., 2010). This structure corresponds to the Aurora Kinase enzyme complexed with an imidazo[1,2-a]pyrazine, the same structural profile of some of the compounds related in this study. Five ligands were used in MD calculations (Figure 2 and Table I): the co-crystallographic (D1); the less active (D2); the second most active (D3); and the most active (D4) from the initial set of 117 imidazole derivatives (respectively, compounds 6o, 17e, and 18h from the article by Garamvölgyi et al., 2016); and one ligand designed after QSAR modelling with the best values of predicted pIC50 (Q3 - Table I). The last four compounds had the imidazo-pyrazine/pyridine nucleus and chains aligned with the co-crystallographic ligand at the binding site before the analyses. All the molecules were parameterized using the ATB platform (Automated Topology Builder) (Malde et al., 2011).

FIGURE 2
Four molecules used in Molecular Dynamics simulations to study the mechanism of action of imidazole derivatives in melanoma cell lines (D1 is the co-crystallographic ligand in Aurora Kinase (PDB ID: 3NRM); D2, D3 and D4 are the less, second and most active compounds in melanoma activity). The compound Q3 is present in Table II, after QSAR results.

TABLE I
QSAR based designed compounds and predicted pIC₅₀ values

Designed compound	Predicted pIC50
	7.14
	6.64
	7.42
	6.65
	6.76
	6.62
	6.76
	6.73
	7.01

The initial geometries for the dynamics of the tested compounds were generated based on the co- crystallographic ligand structure. For each molecule, the imidazopyrazine system was manually adjusted and superimposed on the same location as the original ligand. The complex was minimized by steepest descent, as described below, before the equilibration and production stages. For all the systems, a cubic box involving the entire complex was traced and filled with SPC water molecules, following periodic boundary conditions. Na+ and Cl- were used to neutralize the total charges, replacing solvent molecules with these ions, considering a final physiological concentration of 0.1 M. The long-range electrostatic interactions were modelled using the particle mesh Ewald method, with a cutoff of 1.2 nm (Darden, York, Pedersen, 1993). The same cutoff was used for the calculation of the van der Waals interactions. Bond lengths involving hydrogens and water molecules were initially restrained using the P-LINCS (Hess et al., 2008) and SETTLE (Miyamoto, Kollman, 1992) algorithms, respectively.

The leap-frog algorithm (Van Gunsteren, Berendsen, 1998) was applied to integrate the motion equation. The systems were submitted to an initial energy minimization using the steepest descent algorithm, and then an NVT ensemble was performed at 310 K, with a time step of 2 fs for 100 ps, using the modified Berendsen thermostat. Subsequently, we performed an NPT ensemble under the same conditions, with a temperature of 310 K and a pressure of 1.0 bar, using a Parrinello-Rahman barostat (Parrinello, Rahman, 1991) and the modified Berendsen thermostat. Finally, a full MD simulation of 200 ns was performed at 310 K and 1 bar of pressure, with an integration time step of 2 fs, collecting data every 10 ps. This total time interval was chosen because it is appropriate for specific ligand-protein interactions and due to the equilibrium observed throughout the simulation in RMSD plots and ligand-protein distances.

The binding energies for complexes were calculated based on the MM/PBSA (Molecular Mechanics/Poisson- Boltzmann Surface Area) method (Kollman et al., 2000), using the g_mmpbsa package (Kumari, Kumar, 2014) (https://rashmikumari.github.io/g_mmpbsa/), which neglects the entropic part of Gibbs energy equations but is a user-friendly and open-source tool implemented to interface with GROMACS. The section of the simulation used to calculate the average binding energy and deviation was the last 10 ns of production, with a 200 ps time interval. Hydrogen bonds were identified throughout the simulation based on geometric criteria as described previously by our group (Silva et al., 2017).

RESULTS AND DISCUSSION

QSAR modeling

The molecular descriptors calculated by Dragon software and selected after the Wrapper methodology were: NNRS, GATS5v, GATS7i, P_VSA_LogP_4, P_ VSA_s_5, Eta_betaP, Chi1_EA(dm), SM11_AEA(ri), Mor13v, B07[N-F], F05[N-F], F05[O-O], LLS_01. The relative importance of the variables for the Random Forest analysis can be visualized in Figure 3. Feature importance is computed as the mean of the accumulation of impurity decreases within each tree, also called “Gini Importance.” The 2D structural representations and IC₅₀ values for all 117 imidazo-pyrazines and imidazo-pyridines can be found in the original paper by Garamvölgyi et al. (2016), specifically in Tables 1-9 of that paper.

FIGURE 3
Importance of each variable selected by the Wrapper method.

The validation parameters expressed by Leave-One- Out (LOO) cross-validation (Figure 4) as well as K-Fold analysis (Figure 5) showed good results, indicating a robust and predictive model (coefficient of determination Q² _LOO = 0.641 and Q² _KFOLD with a maximum deviation of 0.1 in relation to LOO). The Y-scrambling analysis (Figure 6) also shows the absence of chance correlation with the selected descriptors, where all the values for LOO in scrambled Y are much smaller than 0.4. The values of Q² for true models in the graph (blue points) show a small variation due to the random nature of the Random Forest algorithm, which may vary slightly during weighting.

FIGURE 4
Leave-One-Out cross-validation results.

FIGURE 5
K-fold results after 50 runs.

FIGURE 6
Y-scrambling results.

Considering the structure-activity relationships based on the most relevant variables, we can initially observe that the most active molecules tend to have higher values of the P_VSA_LogP_4. P_VSA descriptors are based on the approximation of the molecular van der Waals surface area (VSAi) at the atomic level, combined with another molecular property (Pi) within a specific bin (in our case, bin 4) (Todeschini, Consonni, 2009; Labute, 2000). Compounds with more halogen atoms and CH3 groups in the benzene rings tend to exhibit higher values of the descriptor and show better activities (pIC₅₀ > 6.5).

The SM11_AEA(ri) descriptor corresponds to the “spectral moment of order 11 from augmented edge adjacency matrix weighted by resonance integral,” which integrates the “edge adjacency indices” type. This descriptor stores only two types of values for the data series, around 4.609 and 4.616, respectively. The highest values are only present in compounds with pIC₅₀ higher than 6.5, although both values are present in all activity ranges. This fact may reveal some interesting characteristics captured by this variable. The “edge adjacency matrix” (E) is derived from the “molecular graph” (G), which in turn corresponds to a molecular representation where vertices and edges are chemically interpreted as atoms and covalent bonds (Todeschini, Consonni, 2009). Then, E contains information about the connectivity between atoms and is also called the bond matrix. The descriptor considers order 11 for the spectral moment, which is related to the number of paths definable from the molecular edge and is weighted by the resonance integral from the Hückel matrix (Todeschini, Consonni, 2009). We observed that smaller molecules, without substituents in the six-membered ring of the benzimidazole system, have lower values of the descriptor. This must be linked to the smaller number of possible paths as well as the smaller conjugation of these molecules.

The GATS5v corresponds to the “Geary autocorrelation - lag 5 weighted by atomic van der Waals volumes” descriptor (Todeschini, Consonni, 2009). It is part of the set of Geary Autocorrelation Descriptors, defined by:

$GATSdw=\left(\frac{n-1}{2}\right)\left(A\right)/\left(B\right)$

$A=\left(\frac{1}{\mathrm{\Delta }}\right)(\sum _{i=1}^n\sum _{j=1}^n{\delta }_{ij}{\left(w_i-w_j\right)}^2$

$B=\sum _{i=1}^n{\left(w_i-W\right)}^2$

where w _i and w _j are the weights of atoms “i” and “j”, “w” can be represented by mass, volume, electronegativity, among others, and “W” is the average over the entire molecule, being δ_ij a Kronecker delta (Todeschini, Consonni, 2009). In the case of our specific descriptor, this represents the correlation between van der Waals volumes of two atoms separated by a topological distance of 5 bonds. The variable, in general, presents a range of values between 0.85 and 1.2, showing higher values in the most active molecules. Highlighting the three compounds with pIC₅₀ above 7, their values for the descriptor are all above 0.98, while the worst compounds in the series have most values between 0.88 and 0.95. From a structure- activity relationship perspective, and considering this descriptor alone, the presence of halogens, particularly chlorine, promotes a greater difference between atomic van der Waals radii, with that atom being present in multiple forms in the most active compounds. It is also important to highlight that the presence of extra nitrogen in the imidazopyrazine pattern brings greater numerical differences in relation to the carbons separated by 5 bonds, and this descriptor is also related to this characteristic. This is interesting, considering that a structural qualitative analysis of the dataset points to the additional nitrogen in that system as essential for better activity.

GATS7i follows the same concept, involving the correlation between the ionization potentials (IP) of two atoms at a distance of 7 bonds. The three most active molecules contain fluorine and generally have lower values for this descriptor due to the correlation of the fluorines in positions 3 or 4 of the ring with the nitrogen of the urea linker (both atoms with high IP), which does not happen with the less active compounds in the series. The presence of additional nitrogen in the imidazo-pyrazine pattern, together with additional halogens, raises the mean IP, increasing the term in the denominator of the descriptor. Therefore, lower values tend to be associated with more active compounds.

The seventh most important descriptor is the variable F05[N-F], which represents the frequency of fluorine atoms separated from nitrogen atoms by a distance of 5 bonds. This variable is easy to interpret and shows high values in compounds with a pIC₅₀ above 7.0. Observing the structures of these molecules, we found that the presence of urea is again important for the biological effect, especially when associated with fluorine atoms in the immediately neighbouring ring. Additionally, using an amide linker instead of urea does not provide the optimal distance indicated by the descriptor. This is another feature that can be exploited in the design of new compounds.

The other less important descriptors that allowed the final calibration of the model were: Eta_BetaP - “eta pi and lone pair VEM count”; P_VSA_s_5 - “P_VSA-like on I-state, bin 5”; B07[N-F], “Presence/absence of N - F at topological distance 7”; NNRS, “ normalized number of ring systems”; F05[O-O], “Frequency of O - O at topological distance 5”; LLS_01, “modified lead-like score from Congreve et al. (2003) (6 rules)” (Todeschini, Consonni, 2009).

QSAR based design of new compounds

When considering potentially more active structures while preserving the characteristics highlighted by QSAR modelling, some interesting possibilities include obtaining halogenated substitutions in the two phenyl systems at the ends and testing the duplication of the urea linker at both ends of the compounds. We sought to investigate possible analogues for synthesis by considering the synthesis routes from the original article and verifying the commercial availability of the reagents: 3-Chloro-5-(trifluoromethyl)benzoic acid; 3-Chloro-4- (trifluoromethyl)benzoic acid; 3,5-Bis(trifluoromethyl) benzoic acid; 3,5-Bis(chloro)benzoic acid; and 3-Isocyanate-2-methyl-6-(trifluoromethyl)pyridine. We considered the combination of amide/urea linkers as well as urea/urea analogues. These substances are available on the commercial website https://www.fishersci.se/shop. Table I presents nine designed compounds and their respective pIC₅₀ values predicted by the model. The original activities of 117 molecules comprise a pIC₅₀ range between 4.5 and 8.0. All the designed new compounds have pIC₅₀ values higher than 6.5, denoting good activities, around 0.3 μM. The compound Q3 presented a value of 7.42 (~0.04 μM), being a strong candidate for future tests. The use of two urea groups as linkers does not improve activity. However, the presence of 2-methylpyridine results in a compound with a pIC₅₀ of 7.01 (Q9), making it a scaffold that can be used to explore new synthetic analogues.

It is important to emphasize that generating analogous compounds with strong indications of activity is crucial for the success of a new drug. This is because, in drug development, factors such as solubility, pharmacotechnical requirements, and choice of route of administration, among others, require a range of possibilities. The molecule with the highest affinity or activity in vitro does not always become the final drug. Therefore, the model enabled the design of several structures with predicted high activity, providing options for future tests.

Considering the possible mechanism of action for the imidazole derivatives, molecular dynamics simulations bring interesting elucidations to explain the potency of the training molecules, as well as the best designed by QSAR.

Molecular Dynamic simulations

The original article that promotes the synthesis and biological evaluations of the imidazole series studied here was inspired by the structural profile of a B-Raf kinase inhibitor. However, tests of enzymatic inhibition led to the conclusion that, despite the phenotypic cytotoxic activity observed in a melanoma cell line, this proposed mechanism of action was not confirmed. In this sense, and after a search in the Protein Data Bank, we found that a compound with an imidazo[1,2-a]pyrazine nucleus and the same substitution profile at the rings was obtained with the Aurora Kinase A enzyme. This is an important target in studies involving cancer treatment (Tanaka et al., 1999; Cheung et al., 2009). Therefore, we decided to investigate this route as a possible mechanism of action through comparative molecular dynamics simulations.

Figure 7 presents comparative graphs obtained after simulations. Except for RMSF, all the plots include lines representing the 5 ns-period running average. The comparative RMSD plot for C-alphas and the ligands (Figure 7 - A and B, respectively) shows that all the complexed ligands/proteins achieve relative stabilization after 170 ns, with the co-crystallographic (D1) ligand being the least distant from the initial structure. The minimum distance graph (Figure 7-C) shows that the co-crystallographic ligand D1 has the shortest ligand- protein distances along the trajectory, followed by the most active of the imidazole series (D4). However, no significant difference was observed between all of them (maximum range of 0.05 nm considering all the lines).

FIGURE 7
Comparative graphs from Molecular Dynamics simulations (lines in black corresponds to the co-crystallographic ligand (D1); in red, imidazole derivative with low activity (D2); in blue, the second most active compound (D3); in green, the derivative with highest activity (D4); in yellow, QSAR designed best compound (Q3). A - RMSD for protein C-alpha of five complexes; B - RMSD for ligands of five complexes; C - Ligand-protein minimal distances for the five complexes; D - RMSF for residues along the trajectory; E - Ligand-protein hydrogen bonds formed during simulations.

The comparative RMSF graphs do not show significant differences, with attention warranted only by an exclusive peak between 150-200 for the most active compound D4 (Figure 7-D in green), which corresponds to a residue in an alpha-helix near the protein’s binding site. The designed compound (Q3 - yellow) also has a particular fluctuation (residue numbers 200-250), which, however, corresponds to a terminal glycine with no influence on the binding site. Thus, the perturbation in the alpha-helix region of the binding site can be pointed to as important in the affinity and efficacy (D4 - IC₅₀ = 0.01 nM and best binding energy), considering that it does not occur for the other molecules with lower affinity.

Regarding the number of hydrogen bonds, the less active molecule (D2) has a running average of no more than 1 interaction (Figure 7-E). The crystallographic ligand shows a moving average between 2 and 3 interactions. The most active and second most active molecules show a range of 1-3 interactions. The designed QSAR molecule fluctuates around no more than 2 interactions.

Table II shows the estimated binding energy for all the compounds and the specific contributions. Although the binding energy values obtained by the MM/PBSA method are often overestimated, previous studies state that they can be used in a relative way when comparing the stability of complexes (Kumari et al., 2014). As noted, surprisingly, the co-crystallographic ligand presented the highest binding energy (least negative). This can be explained by its smaller size and fewer possible interactions, considering that the compounds in the imidazole series still fill a pocket of more hydrophobic nature, despite showing fewer hydrogen bonds. Electrostatic terms were the most favourable for D1 compared to the other compounds, while van der Waals interactions favoured the investigated molecules. Its larger surface area and extensive contact formation ultimately outweighed the advantage of the rich presence of heteroatoms in the co-crystallographic ligand. Additionally, even the least active molecule in the series tested by Garamvölgyi and co-workers (D2) shows a relatively higher affinity value. The other energies align with the structure-activity relationships in Aurora Kinase, with the least active, the second most active, and the most active compounds showing increasing stabilization at the binding site. These findings support the structural mechanism of action hypothesis.

TABLE II
Binding energies for all the complexes, calculated as the average from frames between 190-200 ns, with an interval dt = 0.2 ns

Terms of Energy*	D1	D2	D3	D4	Q3
Van der Waals	-175.938	-224.848	-265.238	-318.953	-269.889
	+/-9.567	+/-10.180	+/-13.351	+/-15.213	+/-16.619
Electrostatic	-80.919	-36.222	-39.394	-48.017	-43.647
	+/-15.354	+/-12.685	+/- 7.471	+/-8.579	+/-11.465
Polar solvation	177.139	149.791	184.089	203.261	183.177
	+/-24.029	+/-23.119	+/-20.656	+/-15.421	+/-32.250
SASA	-16.224	-23.434	-26.179	-31.435	-29.407
	+/-1.007	+/- 1.349	+/- 1.623	+/- 1.436	+/-1.758
Binding energy	-95.942	-134.713	-146.722	-195.145	-159.766
	+/-13.413	+/-11.822	+/-18.690	+/-22.821	+/-20.526

• *in kJ/mol

Graphical representations of D1-D4 and Q3 within the active site of the protein can be visualized in Figure 8. The pocket has an entry region with two positively charged residues, Arg220 and Arg137, which interact with the ligand region that is more external to the pocket. In the central region, there are residues that have hydrophobic side chains (represented in yellow), such as Leu194, Ala213, and Leu263. Going deeper into the pocket, an Asp274 residue (depicted in red) and a Lys162 residue (depicted in green) make the region more charged and prone to hydrogen bonding.

FIGURE 8
Last frames of simulations showing some interactions in the Aurora Kinase binding site. A: Generic binding pocket of the Aurora Kinase enzyme with crystallographic ligand; B-F: Molecular interactions for D1-D4 and Q3, respectively.

The analysis of the interactions between the ligands and the pocket amino acids can elucidate the inhibition pathway and explain the energy data presented in Table II. The interactions observed with the co-crystallographic ligand D1 were similar to those described by Belanger et al. (2010), with hydrogen bonds between the Ala213 main chain and the imidazo[1,2-a]pyrazine core, and at the innermost end of the pocket between the 3-(4-pyrazolo) group and the side chain of Asp274. In addition to those already described, we also observed a hydrogen bond between the 8-aminoisothiazole group and the Arg220 side chain at the entrance to the pocket. Leu194 and Leu263 (not shown in the image) were closer to the region between the imidazo[1,2-a]pyrazine and 3-(4-pyrazolo) groups with poor stabilization.

As for D2, which presented the second worst binding energy, the interactions were mostly hydrophobic in nature, with van der Waals type interactions prevailing. This finding corroborates the data that make up the total energy value (Table II), as it presented lower solvation energy. The greatest contribution to the result comes from hydrophobic interactions due to its size. It only maintained a stable hydrogen bond with the Arg137 residue at the entrance of the pocket and does not present the characteristic interaction of imidazo[1,2-a]pyrazine ligands. At the end of the pocket, it features a hydrophobic six-membered ring. The hydrophobic residues Ala213, Leu263, and Leu194 (not shown in the image) were not close to D2, resulting in weak and unstable interactions.

D3 has a similar size to D2, with an aromatic ring linked to the imidazo[1,2-a]pyrazine core. However, it has halogens at its ends, which confer less hydrophobicity to the molecule. Additionally, it has a carbonyl diamide group instead of a simple amide as in D2. As a result, the interactions extend the pocket and allow the involvement of other residues such as Glu260, enabling the formation of two hydrogen bonds between its main chain carbonyl and the nitrogens of the diamide carbonyl group of D3. Other donor-acceptor hydrogen bonds are formed with the N-H groups of the Arg220 residue and the amide- like nitrogen of D3. The hydrophobic interactions with Ala213, Leu263 (not shown in the image), and Leu194 were stabilized as they are in a hydrophobic region of the ligand with the aromatic rings.

The D4 ligand was the most active among those observed. It presents a very stabilizing interaction between its carbonyl diamide moiety and the carboxylate group of the side chain of residue Asp274. This residue has previously been reported as important for anchoring imidazo[1,2-a]pyrazine ligands (Belanger et al., 2010). Unlike D3, in which this group interacts with the carbonyl oxygen of the main chain, in D4 it forms a double hydrogen bond between the N-Hs and the oxygens of the aspartate carboxylate. Furthermore, the neighbouring residue, Phe275, also interacts with the ligand, forming a hydrogen bond between its main chain nitrogen and the carbonyl of the carbonyl diamide portion, making this anchoring very efficient. At the entry of the pocket, there is also some stabilization with the interaction with the Pro214 main chain. The other interactions with hydrophobic residues were not highlighted but make up a large part of the binding energy of D4 and Aurora kinase.

The ligand generated from machine learning data, Q3, has hydrogen bonds in its two amide-like groups: one with the N-H and the oxygen of the main chain of Glu260, exactly like D3, and another with the N-H group of Arg137 and its oxygen, as seen in D2. Furthermore, it also replicates an interaction seen in D4, between the carbonyl oxygen of the main chain of Pro214 and the N-H of the amide-like moiety. The structure enabled the formation of a hydrogen bond between the oxygen of the main chain of Leu139 and the N-H of the carbonyl diamide group of Q3, which had not been observed in the other ligands. The hydrophobic rings were stabilized with interactions between amino acids in the centre of the pocket, as highlighted with Leu263, also present in D2 and D3.

It is important to note that the 3,5-trifluoro pattern on the phenyl urea ring, with the same substitution on the opposite ring of the structure, was not reported in the original article. Thus, the Q3 molecule designed from the QSAR, which showed greater affinity than all the original series except for the most active D4, can be exploited in cytotoxicity tests and for the generation of potentially promising future analogues.

By analysing other structural observations, we can observe that the less active one (D2) positions the imidazopyrazine ring opposite to the D1 ligand along the equilibrated trajectory, remaining in an internal occupation region of the protein similar to this co-crystallographic ligand by taking the dimethyl group out of the general pocket (Figure 8). The second (D3) and the most active of the series (D4) maintain the orientation of the rings according to the original ligand. These compounds fill an entire additional pocket at the bottom of Figure 7 that D1 and D2 do not occupy. This may explain the activity and the order of the activity values against melanoma cells (IC50).

The QSAR-projected molecule (Q3), on the other hand, achieves greater stabilization compared to the second most active compound (D3), surpassed only by D4. We observe in Figure 8 that the 3,5-trifluoro groups occupy the associated pocket region. However, a perspective of better accommodation in this region was not obtained compared to D4 (Table II). Comparing the images in Figure 8, the new groups seem to displace the loop formed between GLY142 and GLY145, probably due to a steric disturbance from the bis-trifluoro. However, the RMSD chart along the trajectory does not indicate that this influence was significant in the overall average. Thus, in general, even though it does not have greater affinity than D4, the QSAR-designed Q3 deserves attention and serves to guide new structural proposals. This detailed work can be carried out in continuity with the data presented here.

Analysis of ADME profile

To preview the pharmacokinetic profile of the compounds tested in vitro, comparing it with the co- crystallographic ligand used in the MD calculations, as well as the best compound estimated in our QSAR analyses, some ADME parameters were calculated with the help of the Swiss-ADME platform (Daina, Michielin, Zoete, 2017). The Table III below presents a few theoretical filters that summarize different properties within a drug- likeness analysis. In Lipinski’s rule of five (Lipinski et al., 2001), compounds with a molar weight (MW) < 500 Da, Log(P) < 5, number of hydrogen bond acceptors (NHBA) < 10, and number of hydrogen bond donors (NHBD) < 5 have more desirable characteristics within pharmacokinetic criteria for drugs. The crystallographic compound D1 and the first imidazole studied, D2, do not violate this rule. The compounds with the best in vitro activity, D3, D4, and the QSAR-designed compound Q3, presented some violations. In fact, the large number of atoms and the presence of various halogens led to MW and Log(P) values above those recommended.

TABLE III
ADME parameters calculated for some compounds of this paper

Parameter	D1	D2	D3	D4	Q3
Lipinski	YES	YES	No; MW>500, MLOGP>4.1	No: MW>500, MLOGP>4.1	No; MW>500, MLOGP>4.1
Ghose	YES	No; WLOGP>5.6, MR>130	No; MW>480, WLOGP>5.6, MR>130	No; MW>480, WLOGP>5.6, MR>130	No; MW>480, WLOGP>5.6, MR>130
Veber	YES	YES	YES	No; Rotors>10	No; Rotors>10
Egan	YES	YES	No; WLOGP>5.8	No; WLOGP>5.8	No; WLOGP>5.8
Muegge	YES	No; XLOGP3>5	No; XLOGP3>5	No; MW>600, XLOGP3>5, H-acc>10	No; MW>600, XLOGP3>5, H-acc>10
PAINS	0	0	0	0	0
Brenk	0	0	0	0	0
Synthetic acessibility	3.18	3.14	3.5	3.87	3.96

However, the good pharmacodynamic profile of these compounds opens up the possibility of using carriers, liposomal forms, and other methods to overcome these pharmacokinetic drawbacks. The other parameters include terms such as molar refractivity, number of rotatable bonds, total polar surface area, among others that can be consulted in the references (Ghose, Viswanadhan, Wendoloski, 1999; Veber et al., 2002; Egan, Merz, Baldwin, 2000; Muegge, Heald, Brittelli, 2001). Objectively, analysing the table, both the compounds in the article and the QSAR-designed compound tend to have similar issues with Log(P) and MW. They also show problematic values for molar refractivity (MR), number of rotatable bonds (NRB), and NHBA. The MR parameter can be implicated in a good balance between hydro/lipophilic character, but it can also represent an additional capacity to fill a binding site, as it is related to volume, which correlates positively with the results of the molecular dynamics simulation. The high NRB is due to the high number of single bonds. However, considering that they always represent connections between unsaturated rings and other double bonds or pairs of electrons, this parameter ends up being masked, not revealing the resonance capacity and double character of these single bonds.

The high NHBA can be due to the large number of halogen atoms, which contribute to affinity but can impact solubility. This should be improved by nanoemulsion-type formulations or specific carriers, opening opportunities for studies in this direction. Finally, the PAINS (Baell, Holloway, 2010) and Brenk et al. (2008) alerts show no present violations, meaning that there is maximum selectivity for the target and that the compounds do not present structural characteristics related to drug development problems, respectively. Moreover, the synthetic accessibility of the compounds (Ertl, Schuffenhauer, 2009) was considered at a similar level, with the respective routes well described in the original article.

CONCLUSIONS

Artificial intelligence algorithms, using bioinspired techniques and Random Forest approaches, were able to construct a robust and predictive QSAR model for a series of 117 imidazole derivatives with anti-melanoma in vitro activity, as previously reported in the literature. These results and the structure-activity relationships obtained allowed the design of a series of potentially active analogues, whose IC₅₀ predictions placed them in a good efficiency profile. The generated model can still be used for future designs and tests, in continuation of what is presented here. Additionally, molecular dynamics simulations comparing a known inhibitor of the Aurora Kinase enzyme with the compounds under study revealed a possible anticancer mechanism.

The calculations were able to demonstrate that some representatives of the imidazole series studied in this article showed higher stabilization at the binding site of the protein, providing an insight into this mechanism of action at the molecular level, which had not yet been revealed. This opens up several opportunities to explore the structural pattern presented here, using both in silico and in vitro approaches, to generate promising compounds for application in the treatment of human melanoma.

ACKNOWLEDGEMENTS

The authors are grateful to the National High- Performance Processing Centre of the Federal University of Ceará - CENAPAD-UFC and the High-Performance Computing Centre at UFRN (NPAD/UFRN) for the computer cluster facility. Molecular graphics and analyses were performed with UCSF Chimera, developed by the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco, with support from NIH P41-GM103311.

REFERENCES

• Abraham MJ, Murtola T, Schulz R, Páll, S, Smith JC, Hess B, et al. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. Software X. 2015;1-2:19-25.
• ACD/Labs. ChemSketch, 2022;version 2022.1.2.
• Adams RC, Rashidieh B. Can computers conceive the complexity of cancer to cure it? Using artificial intelligence technology in cancer modelling and drug discovery. Math Biosci Eng. 2020;17(6):6515-6530.
• Alan J, Karla B, Nikola B. A review of feature selection methods with applications. 38th International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO). 2015; 10.1109/MIPRO.2015.7160458.
  » https://doi.org/10.1109/MIPRO.2015.7160458
• Arnold M, Holterhues C, Hollestein LM, Coebergh JWW, Nijsten T, Pukkala E, et al. Trends in incidence and predictions of cutaneous melanoma across Europe up to 2015. J Eur Acad Dermatol Venereol. 2014;33(7):1214-1223.
• Arnold M, Singh D, Laversanne M, Vignat J, Vaccarella S, Metheus F, et al. Global Burden of Cutaneous Melanoma in 2020 and Projections to 2040. JAMA Dermatol. 2022;158(5):495-503.
• Baell JB, Holloway GA, New substructure filters for removal of pan assay interference compounds (PAINS) from screening libraries and for their exclusion in bioassays. J Med Chem. 2010;53:2719-2740.
• Belanger DB, Curran PJ, Hruza A, Voigt J, Meng Z, Mandal AK, Discovery of imidazo[1,2-a]pyrazine- based Aurora kinase inhibitors. Bioorg Med Chem Lett. 2010;20(17):5170-4.
• Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, et al. The Protein Data Bank. Nucleic Acids Res. 2000;28(1):235-42.
• Bolaji ALA, Khader AT, Al-Betar MA, Awadallah MA. University course timetabling using hybridized artificial bee colony with hill climbing optimizer. J Comput Sci. 2014;5(5):809-818.
• Bolón-Canedo V, Sánchez-Maroño N, Alonso-Betanzos A. A review of feature selection methods on synthetic data. Knowl Inf Syst. 2013a;34:483-519.
• Bolton EE, Chen J, Kim S, Han L, He S, Shi W, et al. PubChem3D: a new resource for scientists. J Cheminform. 2011;3(1):32.
• Breiman L. Random Forests. Mach Learn. 2001;45:5-32.
• Brenk R, Schipani A, James D, Krasowski A, Gilbert IH, Frearson J, et al. Lessons learnt from assembling screening libraries for drug discovery for neglected diseases. ChemMedChem 2008;3:435-444.
• Brown M, Bossley KM, Mills DJ, Harris CJ. High dimensional neurofuzzy systems: overcoming the curse of dimensionality” Proceedings of 1995 IEEE International Conference on Fuzzy Systems. 1995;10.1109/ FUZZY.1995.409976.
  » https://doi.org/10.1109/ FUZZY.1995.409976
• Carracedo-Reboredo P, Liñares-Blanco J, Rodríguez- Fernández N, Cedrón F, Novoa FJ, Carballal A, et al. A review on machine learning approaches and trends in drug discovery. Comput Struct Biotechnol J. 2021;12(19):4538- 4558.
• Chan HCS, Shan H, Dahoun T, Vogel H, Yuan S. Advancing Drug Discovery via Artificial Intelligence. Trends Pharmacol Sci. 2019;40(8):592-604.
• Cheung CHA, Coumar MS, Hsieh HP, Chang JY, Coumar MS, Cheung CH, et al. Aurora kinase inhibitors in preclinical and clinical testing. Expert Opin Investig Drugs. 2009;18(4):379-98.
• Congreve M, Carr R, Murray C, Jhoti H. A ‘rule of three’ for fragment-based lead discovery? Drug Discov Today. 2003;8(19):876-7.
• Cramer SF. The origin of epidermal melanocytes. Implications for the histogenesis of nevi and melanomas. Arch Pathol Lab Med. 1991;115:115-119
• Daina A, Michielin O, Zoete V. SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci Rep. 2017;7:42717.
• Darden T, York D, Pedersen L. Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. J Chem Phys. 1993;98:10089-10092.
• Du R, Huang C, Liu K, Dong Z. Targeting AURKA in Cancer: molecular mechanisms and opportunities for Cancer therapy. Mol Cancer. 2021;20(1):15.
• Egan WJ, Merz KM, Baldwin JJ. Prediction of Drug Absorption Using Multivariate Statistics. J Med Chem . 2000;43:3867-3877.
• Erdmann F, Lortet-Tieulent J, Schüz J, Zeeb H, Greinert R, Breitbart EW, et al. International trends in the incidence of malignant melanoma 1953-2008-are recent generations at higher or lower risk? Int J Cancer. 2013;132(2):385-400.
• Ertl P, Schuffenhauer A. Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions. J Cheminform . 2009;1:8.
• Fatma S, Cagla K, Besra YO, Aycan A, Roya G, Sezgi K, et al. The Evaluation of Effect of Aurora Kinase Inhibitor CCT137690 in Melanoma and Melanoma Cancer Stem Cell. Anticancer Agents Med Chem. 2021;21(12):1564- 1574.
• Garamvölgyi R, Dobos J, Sipos A, Boros S, Illyés E, Baska F, et al. Design and synthesis of new imidazo[1,2-a] pyridine and imidazo[1,2-a]pyrazine derivatives with antiproliferative activity against melanoma cells. Eur J Med Chem . 2016;108:623-643.
• GCO - Global Cancer Observatory, WHO - World Health Organization. 2023. https://gco.iarc.fr Accessed 29 Jun 2023
  » https://gco.iarc.fr
• Ghose, AK, Viswanadhan, VN, Wendoloski, JJ. 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. 1999;1:55-68.
• Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov. 2022;12(1):31-46.
• Hart PE, Nilsson NJ, Raphael, B. A formal basis for the heuristic determination of minimum cost paths. Systems Science and Cybernetics, IEEE Transactions. 1968;10.1109/TSSC.1968.300136.
  » https://doi.org/10.1109/TSSC.1968.300136
• Hess B, Kutzner C, Van Der Spoel D, Lindahl E. GROMACS 4: algorithms for highly efficient, load- balanced, and scalable molecular simulation. J Chem Theory Comput. 2008;4(3):435-47.
• Karaboga D, Basturk B. A powerful and efficient algorithm for numerical function optimization: artificial bee colony (ABC) algorithm. J Glob Optim. 2007;39:459-471.
• Kim S, Bolton EE, Bryant SH. PubChem3D: conformer ensemble accuracy. J Cheminform . 2013;5(1):1.
• Kim S, Chen J, Cheng T, Gindulyte A, He J, He S, et al. PubChem 2023 update. Nucleic Acids Res . 2023;51(D1):D1373-D1380.
• Klanbauer G, Hochreiter S, Rarey M. Machine Learning in Drug Discovery. J Chem Inf Model. 2019; 59(3):945-946.
• Kode Cheminformatics. Dragon version 7.0.10, 2017.
• Kollman PA, Massova I, Reyes C, Kuhn B, Huo S, Chong L, et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc Chem Res. 2000;33(12):889-97.
• Kumari R, Kumar R. Open-Source Drug Discovery Consortium, Lynn A. g_mmpbsa - A GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inf. 2014;54(7)1951-1962.
• Labute P. A widely applicable set of descriptors. J Mol Graph Model. 2000;18(4-5):464-77.
• Lipinski, CA, Lombardo F, Dominy, BW, Feeney, PJ. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev. 2001;46: 3-26.
• Lo YC, Rensi SE, Torng W, Altman RB. Machine learning in chemoinformatics and drug discovery. Drug Discov Today . 2018;23(8):1538-1546.
• Majumdar S, Basak SC. Beware of External Validation! - A Comparative Study of Several Validation Techniques used in QSAR Modelling. Curr Comput Aided Drug Des. 2018;14(4):284-291.
• Mak KK, Pichika MR. Artificial intelligence in drug development: present status and future prospects. Drug Discov Today . 2019;24(3):773-780.
• Malde AK, Zuo L, Breeze M, Stroet M, Poger D, Nair PC, et al. An Automated Force Field Topology Builder (ATB) and Repository: Version 1.0. J Chem Theory Comput 2011;7(12):4026-37.
• Miyamoto S, Kollman PA. Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J Comput Chem. 1992; 13(8)952-962.
• Muegge I, Heald SL, Brittelli D. Simple selection criteria for drug-like chemical matter. J Med Chem . 2001;44;1841-1846.
• Ogilvie LA, Kovachev A, Wierling C, Lange BMH, Lehrach H. Models of models: A translational route for cancer treatment and drug development. Front Oncol. 2017;19:7:219.
• Parrinello M, Rahman A. Polymorphic transitions in single crystals: A new molecular dynamics method. J Appl Phys. 1981;52:7182-7190.
• Punt S, Malu S, McKenzie JA, Manrique SZ, Doorduijn EM, Mbofung RM, et al. Aurora kinase inhibition sensitizes melanoma cells to T-cell-mediated cytotoxicity. Cancer Immunol Immunother. 2021;70(4):1101-1113.
• Rebecca VW, Somasundaram R, Herlyn, M. Pre-clinical modeling of cutaneous melanoma. Nat Commun. 2020;11(1):2858.
• Sellwood MA, Ahmed M, Segler MHS, Brown N. Artificial intelligence in drug discovery. Future Med Chem. 2018;10(17)2025-2028.
• Shaker B, Ahmad S, Lee J, Jung C, Na D. In silico methods and tools for drug discovery. Comput Biol Med. 2021;137:104851.
• Silva SB, Pinheiro MP, Fuzo CA, Silva SR, Ferreira TL, Lourenzoni MR, et al. The role of local residue environmental changes in thermostable mutants of the GH11 xylanase from Bacillus subtilis. Int J Biol Macromol. 2017;97:574-584.
• Talevi A, Morales JF, Hather G, Podichetty JT, Kim S, Bloomingdale PC, et al. Machine Learning in Drug Discovery and Development Part 1: A Primer. CPT Pharmacometrics Syst Pharmacol. 2020;9(3):129-142.
• Tanaka T, Kimura M, Matsunaga K, Fukada D, Mori H, Okano, Y. Cancer Res. 1999;59:2041-2044.
• Todeschini R, Consonni V. Molecular descriptors for chemoinformatics. Wiley Online Library. 2009;10.1002/9783527628766.
  » https://doi.org/10.1002/9783527628766
• Vamathevan J, Clark D, Czodrowski P, Dunham I, Ferran E, Lee G, et al. Applications of machine learning in drug discovery and development. Nat Rev Drug Disco. 2019;18(6):463-477.
• Van Gunsteren WF, Berendsen HJC. A leap-frog algorithm for stochastic dynamics. Mol. Simul. 1988;1(3):173-185.
• Veber DF, Johnson SR, Cheng HY, Smith BR, Ward KW, Kopple KD. Molecular properties that influence the oral bioavailability of drug candidates. J Med Chem . 2002;45:2615-2623.
• WHO - World Health Organization, Cancer. 2022 2022 https://www.who.int/news-room/fact-sheets/detail/cancer Accessed 30 Jun 2023.
  » https://www.who.int/news-room/fact-sheets/detail/cancer

Edited by

• Associated Editor:
  Silvya Stuchi Maria-Engler

Publication Dates

• Publication in this collection
  20 Jan 2025
• Date of issue
  2025

History

• Received
  20 June 2024
• Accepted
  20 June 2024