Computational Investigations on Inhibitors of Mycobacterium tuberculosis Shikimate Kinase: Machine Learning, Docking, Molecular Dynamics and Free Energy Calculations

Abstract

Shikimate kinase emerges as an intriguing macromolecular target for the development of novel pharmaceutical agents for the treatment of tuberculosis. This study aimed to develop a neural network (NN) for the discovery of potential inhibitors of Mycobacterium tuberculosis shikimate kinase and to conduct molecular docking and molecular dynamics (MD) simulations. The NN model pointed out to a set of 810 molecules with anti-tuberculosis activity, wherein 86% of this set also demonstrated positive outcomes according to docking calculations. Among these, 54 molecules exhibited a docking score ranging from -9 to -9.8 kcal mol^-1. Subsequently, a subset of molecules was selected for molecular dynamics studies and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) calculations. Furthermore, it was possible to assess that the dataset with higher affinity shared a similar electronic profile, as evidenced by the analysis of global descriptors (electronic chemical potential, hardness, and electrophilicity). The molecules displaying the lowest Gibbs free energy (∆G)_binding values, therefore the highest affinity, were identified as CHEMBL1229147, CHEMBL4095667, and CHEMBL120640.

Keywords:
shikimate kinase; machine learning; tuberculosis; docking; molecular dynamics

Introduction

Tuberculosis (TB), an infectious disease of a chronic and debilitating nature, holds a deeply rooted historical presence in human health. Commonly caused by Mycobacterium tuberculosis (MTB), the disease primarily affects the lungs, though manifestations can occur in other organs as well. Its global impact has been profound, both in terms of public health and socio-economic dimensions.¹,²

MTB is a microorganism possessing a complex cell wall composed of lipids and peptidoglycans, contributing to its resilience against adverse conditions and therapeutic interventions. Bacterial dissemination primarily occurs through respiratory droplets, with inhalation being the primary mode of infection.³

Despite advancements in diagnosis, treatment, and prevention, tuberculosis remains a formidable challenge on the global health stage. The emergence of drug-resistant strains, known as multidrug-resistant tuberculosis (TB MDR) and extensively drug-resistant tuberculosis (TB-XDR),⁴ further complicates disease control. Additionally, co-infection with human immunodeficiency virus (HIV) increases susceptibility and severity of tuberculosis, especially in vulnerable populations.

Traditionally, drug-susceptible TB was treated with a 6-month regimen, which has proven effective in curing most patients. However, the duration of treatment has always been a challenge, both for patients and healthcare systems. Research has sought to reduce treatment time without compromising efficacy. The addition of fluoroquinolones, such as moxifloxacin or gatifloxacin, to conventional regimens was one approach explored.⁵,⁶ Although initially promising, larger studies failed to demonstrate that these shorter regimens were as effective as the standard 6-month treatment.⁶

Another strategy has been to optimize the pharmacokinetics of existing drugs. Increasing doses of rifampicin, for example, have been shown to increase bactericidal activity, but the safety and effectiveness of this approach in different patient groups are still being investigated.⁷ A significant advance was the discovery that rifapentine, a rifamycin with a longer half-life, in combination with other drugs, can reduce treatment to 4 months. This finding led the World Health Organization (WHO) to conditionally recommend this shorter regimen for patients with pulmonary drug susceptible TB.⁶,⁸

The challenges encountered in treating pre-extensively drug-resistant tuberculosis (pre-XDR TB) are that current regimens have limited efficacy, with a high profile of adverse events.⁹ Treatment is expensive, making it inaccessible to most countries. Regimens such as BPaL (bedaquiline, pretomanid and linezolid) cause significant side effects, such as myelosuppression and neurotoxicity.¹⁰ Additionally, monitoring side effects in resource-limited settings is challenging. The BEAT-India (bedaquiline-linezolid-delamanid-clofazimine) clinical trial specifically recruited people with pre-XDR TB and used a four-drug regimen (bedaquiline, linezolid, clofazimine, and delamanid) for 6-9 months, with good results but a high rate of side effects.¹¹

In this context, the pursuit of new macromolecular targets for drug development capable of combating different MTB strains becomes imperative.¹,²,¹²

The Mycobacterium tuberculosis enzyme shikimate kinase (SK) is a key component of the shikimate pathway metabolism, playing a crucial role in the synthesis of essential precursors for the biosynthesis of aromatic acids and other compounds fundamental to the survival and pathogenicity of this bacterium. The also called shikimate pathway is a unique metabolic route found in mycobacteria and some other clinically significant bacteria, making it a potential target for the development of new antimicrobial agents.¹³,¹⁴

The protein is synthesized through the expression of the gene designated as aroK, consisting of a sequence of 176 amino acid residues. SK plays a crucial role by catalyzing the fifth step of the shikimate acid metabolic pathway. Its enzymatic activity involves the conversion of adenosine triphosphate (ATP) and shikimate into adenosine diphosphate and shikimate-3-phosphate through a phosphate group transfer reaction.¹³,¹⁵

In a study,¹⁶ researchers used in silico methods to identify 15 benzothiazole-derived molecules as potential inhibitors of Mycobacterium tuberculosis shikimate kinase. Two of these showed significant inhibitory activity with maximal inhibitory concentration (IC₅₀) values of 10.69 ± 0.9 and 46.22 ± 1.2 μM. Similarly, another study¹⁷ using docking techniques were able to perform a virtual screening and find two compounds reported to have IC₅₀ values of 1.6 and 2.8 μM. The discovery of these inhibitors highlights the potential for developing new drugs targeting shikimate kinase.¹⁸ In silico screening and subsequent in vitro studies¹⁸,¹⁹ have shown that these inhibitors bind effectively to the target enzyme, making them promising candidates for drug development.

Within this context, SK is an intriguing macromolecular target for the development of new drugs. However, the drug discovery process is complex and dependent on numerous factors; thus, the use of computational techniques provides rapid and cost-effective methods for identifying new drug prototypes or candidates.²⁰ Among these techniques, machine learning offers tools for classification or prediction of properties that are both concise and robust.²¹

In recent years, the field of machine learning (ML) has witnessed exponential growth, emerging as one of the most relevant disciplines at the intersection of computer science and statistics. ML aims to develop algorithms and techniques that enable computational systems to learn patterns from data without the need for explicit programming. This approach holds invaluable potential for a broad range of applications, from speech and image recognition to medical diagnostics, financial data analysis, and drug discovery.²²,²³

ML algorithms have been employed to develop various models for predicting chemical, biological, and physical characteristics of compounds in drug discovery. These methods can be integrated into all stages of the drug discovery process. For instance, ML algorithms have been utilized to discover new uses for drugs, predict drug-protein interactions, assess drug efficacy, identify safety biomarkers, and optimize the bioactivity of molecules.²²

Two main approaches characterize these methods: supervised learning (SL) and unsupervised learning (UL).²⁴,²⁵ SL relies on knowledge acquisition through pre-labeled training examples, aiming for the inference of labels for new samples. On the other hand, UL focuses on identifying structures and regularities underlying a set of samples, typically devoid of labels.²⁶ Within the SL category, there are several popular algorithms used to predict new drugs, among which we have the k-nearest neighbors (KNN) and support vector machine (SVM).²¹,²⁶

Deep learning is a subset of machine learning, characterized by a neural network with a large number of layers and parameters. Most deep learning methods employ neural network architectures and are therefore also referred to as deep neural networks.²⁷,²⁸

In this structure, each module at a level transforms its input into a more abstract and elevated representation. In this context, the term “deep” is associated with the number of layers in the network, more layers imply a deeper network. The entire neural network is composed of multiple processing units, referred to as artificial neurons (AN). As depicted in equations 1 and 2, the AN has synaptic inputs, which are weighted by their respective synaptic weights. These input signals are summed, and subsequently, the activation function is applied, responsible for constraining the amplitude of the output of the neuron. The bias serves to increase or decrease the net input to the activation function.²⁶,²⁷

(1) $${\mathrm{u}}_{\mathrm{k}}=\sum _{\mathrm{j}=1}^{\mathrm{m}}{\mathrm{w}}_{\mathrm{kj}}{\mathrm{x}}_{\mathrm{j}}$$

(2) $${\mathrm{y}}_{\mathrm{k}}=\phi ({\mathrm{u}}_{\mathrm{k}}+{\mathrm{b}}_{\mathrm{k}})$$

where x_j represents the input signals, w_kj represents the respective synaptic weights, u_k is the output of the linear combiner, b_k is the bias, φ (.) is the activation function, and y_k is the output of the neuron.

A dense neural network (DNN), also known as a fully connected neural network, is a type of artificial neural network architecture in which each neuron in a layer is connected to all neurons in the next layer. This implies that each neuron receives information from all neurons in the previous layer and sends its outputs to all neurons in the next layer. In other words, all inputs are connected to all outputs between consecutive layers in the network.²⁹

The approach employed through machine learning for pattern recognition and value prediction proves to be highly efficient when applied in the context of rational drug development, especially when used in conjunction with virtual screening and molecular dynamics calculations.²⁰,²³ This integrated approach provides more substantial insights into the studied macromolecular target.

The general scheme of a virtual screening (VS) strategy using molecular docking as a tool begins with the processing of three-dimensional (3D) structural information of the target of interest. The target structure can be derived from experimental data (X-ray, nuclear magnetic resonance (NMR), or neutron scattering spectroscopy), homology modeling, or molecular dynamics (MD) simulations.³⁰,³¹ Molecular docking attempts to efficiently predict the non-covalent binding of macromolecules or a macromolecule (receptor) and a small molecule (ligand) starting from their unbound structures.³² With the aid of molecular docking calculations, compounds can be classified, and these are then post-processed through the analysis of calculated binding scores and validation of the generated pose.

This approach allows us to filter and identify potential compounds with affinity for the target macromolecule. The results, however, are based on an environment where the protein is static, neglecting the dynamic aspects of the binding process. To address this, more in-depth studies are necessary using molecular dynamics calculations.

Molecular dynamics simulation is grounded in the principles of classical mechanics and provides information on the time-dependent microscopic dynamic behavior of individual atoms constituting the system. Newton’s laws of motion are applied to predict the spatial position of each atom over time. Specifically, as time progresses, forces on each atom are repeatedly calculated and then used to update the position and velocity of each atom. The resulting trajectory is essentially a three-dimensional movie describing the atomic-level configuration of the system at each point during the simulated time interval.³³,³⁴

The simulations enable the capture of the positions and movements of all atoms in a biomolecular system over time and provide precisely controlled conditions, encompassing aspects such as the initial conformation of a protein, the presence of ligands, mutations or post-translational modifications, as well as other molecules in the surrounding environment, protonation state, temperature, tension across membranes, and other relevant parameters. Through the comparison of simulations conducted under different scenarios, it is possible to identify the effects resulting from various molecular perturbations.³³,³⁴

Another important aspect for understanding the interaction between a protein and a ligand is the binding free energy, a property that enables us to estimate the strength of this interaction. The molecular mechanics with Poisson-Boltzmann and surface area solvation (MM/PBSA) approach is widely employed due to its computational cost-effectiveness and reasonable estimations.³⁵,³⁶

In this method, the estimation of the free energy can be calculated using the equations 3 and 4:³⁵

(3) $$\mathrm{\Delta }{\mathrm{G}}_{\text{binding}}=⟨{\mathrm{G}}_{\mathrm{PL}}⟩-⟨{\mathrm{G}}_{\mathrm{P}}⟩-⟨{\mathrm{G}}_{\mathrm{L}}⟩$$

(4) $$\mathrm{G}={\mathrm{E}}_{\mathrm{bnd}}+{\mathrm{E}}_{\mathrm{el}}+{\mathrm{E}}_{\mathrm{vdw}}+{\mathrm{G}}_{\mathrm{pol}}+{\mathrm{G}}_{\mathrm{np}}-\mathrm{TS}$$

where PL, P, and L are, respectively, simulations of the complex, apo protein, and free ligand. Equation 4 of the state is composed of terms E_bnd (bond, angle, and dihedral), electrostatic, and van der Waals originating from molecular mechanics. The G_pol and G_np are the polar and nonpolar contributions, where G_pol is estimated from the Poisson-Boltzmann equation and G_np is estimated from the solvent-accessible surface area. The T is the absolute temperature, and S is the entropy.

To compute all averages, the MM/PBSA method typically commences with a molecular dynamics (MD) simulation of the complex using a single-trajectory approach, or three separate MD simulations of the complex, receptor, and ligand, respectively, in the multitrajectory approach.³⁵,³⁶ Within this context, assessments conducted by the MM/PBSA approach can be utilized to better rationalize structure-activity relationships and ligand selectivity profiles.³⁶

Density functional theory³⁷ (DFT) provides a powerful tool to understand and predict how molecules react. Descriptors such as chemical potential, Fukui function, global hardness, and global softness are widely used to characterize this reactivity and identify the most reactive sites. According to this theory, interactions between molecules can be classified into two types: mutual polarization followed by orbital interactions and Coulomb forces.³⁸ The intensity of these effects in molecules can be evaluated locally by the Fukui function and local hardness. And the use of these descriptors can aid in understanding protein-ligand interactions.³⁹

Hence, the objective of this study was to develop a neural network for the discovery of potential inhibitors of Mycobacterium tuberculosis shikimate kinase and perform molecular docking, molecular dynamics simulations, binding free energy estimation and calculation of descriptors.

Methodology

A comprehensive prediction pathway was devised, involving the application of a dense neural network. In this approach, the model forecasts anti-tuberculosis activity. Subsequently, molecules predicted to be active against tuberculosis undergo virtual screening with the shikimate kinase protein of M. tuberculosis. Following this step, six molecules were selected based on their calculated affinity through molecular docking for molecular dynamics simulations and estimation of binding free energy.

Database preparation

The ChEMBL database⁴⁰,⁴¹ was employed to acquire molecular structures, totaling 1,372,848 molecules. Initial data filtration involved the following constraints: (i) Lipinski’s rule: molecular weight (MW) ≤ 500 g mol^-1, number of hydrogen bond acceptor (NHBA) ≤ 10, number of hydrogen bond donor (NHBD) ≤ 5, partition coefficient between n-octanol: water (Log P) ≤ 5, resulting in 1,234,604 molecules. Further refinement was carried out with the second restriction: (ii) aromatic rings ≥ 1, NHBA ≥ 2, NHBD ≥ 2 based on the structure of the shikimate kinase substrate, and MW ≤ 300 g mol^-1 based on the rule of 3,⁴² resulting in a total of 83,271 molecules. These molecules were organized into a dataset containing the names and molecular structures in SMILES format.

Training set

The training set was obtained through the CDD platform,⁴³ we selected only compounds that underwent in vivo or in vitro testing, as we wanted to focus on molecules that either penetrated or failed to penetrate the mycobacterial cell wall. This resulted in a total of 961 molecules, with 20% set aside for testing. The cutoff for labeling compounds as active or inactive was based on a minimum inhibitory concentration of ten micromolar.

Development of the dense neural network

Descriptors

ECFP4 and ECFP6 (extended-connectivity fingerprints) descriptors are circular topological fingerprints, the atom environment for each atom becomes a feature component of an ECFP. Each feature generates a hash number, and a collection of the hash numbers forms an ECFP.⁴⁴,⁴⁵ In this study we work with fingerprints of diameter 4 and 6 with 2048 bins constructed using the RDKit library in Python.⁴⁶ These descriptors were calculated for the training set and organized into a data frame. The same procedure was applied to the molecules in the aforementioned database, creating a new data frame used for model predictions.

Optimization, training, and prediction

The model was evaluated using the following statistics: accuracy, receiver operating characteristic (ROC) curve, area under the curve (AUC), specificity, and recall. The metrics are defined using the abbreviations TP (true positive), TN (true negative), FP (false positive), and FN (false negative).

Accuracy is defined according to equation 5:

(5) $$\text{accuracy}=\frac{(\mathrm{TN}+\mathrm{TP})}{(\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN})}$$

Specificity is defined according to equation 6:

(6) $$\text{specificity}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}$$

Recall is defined according to equation 7:

(7) $$\text{recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

The ROC curve is constructed based on recall versus false positive rate (equation 8):

(8) $$\text{false positive rate}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TP}}$$

AUC represents the percentage of the ROC curve, serving as a metric that encapsulates the performance of the classifier. Hyperparameter optimization was conducted using a three-layer dense neural network developed through the Python TensorFlow⁴⁷ library and the grid search method of Scikit-learn. The network architecture was adjusted, as a deep network with many hidden layers can lead to overfitting when trained on a small dataset. The hyperparameter values considered included batch size (16, 32, 128), number of epochs (20, 50, 100), learning rate (0.001, 0.0001), L2 regularizer to mitigate model overfitting (0.00001, 0.00002, 0.0001), and neurons for hidden layers (16, 8, 5). The optimal values for batch size, epochs, learning rate, L2, and neurons were determined as 16, 50, 0.0001, 0.00002, and 8, respectively, and were selected for model training.

The neural network architecture employed for training consisted of three hidden layers, utilizing the rectified linear unit (ReLU) activation function for the hidden layers and a sigmoid activation function for the output layer. Binary cross-entropy was employed as the loss function, and the Adam optimization algorithm was utilized.⁴⁸ The model training involved a 5-fold cross-validation. Subsequently, the trained model was applied to the aforementioned dataset to predict molecules with anti-tuberculosis activity.

Virtual screening

By combining virtual screening with a deep learning model, we aim to predict molecules with potential anti-tuberculosis activity. The model generates a set of molecules, from which we select only those that interact with our target protein, thus guiding the rational design of potential inhibitors.

For the implementation of virtual screening, the AutoDock Vina software⁴⁹,⁵⁰ was employed. The structure of the shikimate kinase complexed with the substrate was 2IYQ,⁵¹ obtained from the Protein Data Bank (PDB). The 3D structures of the molecules predicted by the model were sourced from the ChEMBL platform.⁴⁰,⁴¹ The grid box was centered on the protein substrate binding site, with the following parameters: center_x = 64.741, center_y = 15.77, center_z = 12.789, size_x = 30, size_y = 30, size_z = 30. Molecular docking calculations were performed in triplicate. The cutoff point was defined as the value obtained from the docking calculation of shikimate with the substrate. Molecules with a more negative affinity than this cutoff point were labeled as active, while those with an equal or less negative affinity than the cutoff point were labeled as inactive, concerning the shikimate kinase.

Subsequently, three groups were created based on the affinity values obtained from the molecular docking calculations, namely: highly active (A), moderately active (B) and weakly active (C). Two molecules from each group were chosen for further molecular dynamics simulations. These selected sets were then analyzed for binding modes using the LigPlot software.⁵²

Analysis of molecular properties

To analyze the molecular properties, Log P, MW, NHBA NHBD, NRB, ring count (RC), and topological polar surface area (TPSA) were calculated using the Python library RDKit.⁴⁶

Analysis of its absorption, distribution, metabolism, and excretion (ADME) was performed using a web tool SwissADME⁵³ and for the in silico toxicity calculation of compounds ProTox-II,⁵⁴ this analysis was performed only for the top three compounds.

Molecular dynamics calculations

The molecules selected for molecular dynamics simulation originated from the outcomes predicted by the dense neural network and virtual screening, constituting a total of six molecules labeled as groups A, B, and C.

The initial configuration of the protein was obtained from the Protein Data Bank, and it was the same configuration used in the previously described virtual screening. The ligand topologies were constructed using the ACPYPE tool.⁵⁵ The simulations were conducted using GROMACS software, version 4.5.5.⁵⁶ The first step involved using the pdb2gmx program⁵⁷ to generate the topology file and the Cartesian coordinate file of the protein under analysis. Subsequently, the systems were centered in a cubic box, with faces at least 1.5 nm away from the protein. The box was filled with TIP3P water molecules⁵⁷ and the AMBER99SB force field.⁵⁸ Ions NaCl were then added at a concentration of 0.154 mol L^-1. Following this, system minimization was performed, followed by system equilibration in the NVT (N: number of atoms; V: volume, T: temperature) ensemble for 100 ps and subsequently in the NPT (N: number of atoms; P: pressure, T: temperature) ensemble for another 100 ps. Finally, the production dynamics were conducted at constant temperature (310 K) and pressure (1 bar) using the V-rescale⁵⁹ and Parrinello and Rahman⁶⁰ algorithms. A 1 nm cutoff radius was defined for Coulomb and van der Waals interactions. The equations were integrated using the leap-frog algorithm with a time step of 2 fs. The total simulation time was 100 ns (within this interval, structural stability of the system has already been established), performed in duplicates.

Analyses were carried out using GROMACS software⁵⁶ programs, including root mean square deviation (RMSD), radius of gyration (RG), root mean square fluctuation (RMSF), number of hydrogen bonds, and average distance.⁶¹

MM/PBSA

The estimation of binding energy for the protein-ligand binary complex and the assessment of residue contributions to binding energy were conducted through the MM/PBSA method utilizing the g_MMPBSA software.⁶²,⁶³ Trajectories derived from the simulations were compressed into a condensed trajectory for the purpose of calculations, performed at a temperature of 310 K.

Docking on relaxed structure

The docking of the previously selected six inhibitors and substrate was also performed with the shikimate kinase structure obtained from molecular dynamics, which was extracted from the last frame of the simulation trajectory. AutoDock Vina software⁵⁰ was also employed for these calculations to assess the maintenance of molecular ranking and their interactions with the protein.

Calculation of global descriptors and local descriptors

For a more comprehensive analysis of electronic properties, global descriptors were computed, including electron affinity (EA), ionization potential (IP), electronic chemical potential (ECP), electrophilicity, hardness, and softness.⁶⁴,⁶⁵ Initially, geometry optimization calculations were performed on the molecular structures, followed by single-point calculations of the six selected ligands and the substrate. Both calculations used the semi-empirical PM7 method in the MOPAC software.⁶⁶ The PRIMoRDIA software⁶⁴ was utilized for descriptor calculations, employing the frozen orbital approach (FOA) method.

Subsequently, the obtained data was organized into a data frame, and a correlation analysis between variables was conducted. Variables such as electrophilicity, hardness, and ECP were selected for further analysis. Following this, principal component analysis (PCA) was performed. Interactions were analyzed using the Discovery Studio software.⁶⁷

An evaluation of local hardness and Fukui function descriptors was conducted for protein-ligand complexes, the six selected ligands, and the substrate, all obtained from molecular dynamics simulations. For all semi-empirical calculations, including geometry optimization and single-point calculations, we employed the linear scaling algorithm MOZYME implemented in the MOPAC software⁶⁶ with the PM7 Hamiltonian. The descriptors were computed using the PRIMoRDIA software.⁶⁴

Pharmacophore pattern determination

First, we extracted pharmacophoric features by ligand using the Python library RDKit.⁴⁷ We calculated the frequency of these features, selecting the most frequent ones, and subsequently organized this data into a data frame. Next, we generated ensemble pharmacophores by employing k-means clustering to group features by type: donors, acceptors, hydrophobics, and aromatics, with their respective k-values set to 2, 3, 3, and 3.

Results and Discussion

The methodology described in sub section “Development of the dense neural network” was applied to the test set outlined in sub section “Database preparation”, and the results indicated that architecture employed in the DNN proved to be capable of achieving satisfactory performance. When evaluating the diameter size of the fingerprint to assess potential improvements in classification, the results depicted in Table 1 indicate the absence of significant differences. Furthermore, in Figure 1, it is evident that the ROC curves associated with ECFP4 and ECFP6 exhibit similar behavior.

Table 1
Assessment of the fingerprint diameter

Fingerprint	Accuracy	AUC	Recall	Specificity
ECFP4	0.80	0.88	0.77	0.84
ECFP6	0.80	0.87	0.76	0.84

Figure 1
5-Fold cross-validation ROC plots of ECFP6 (a) and ECFP4 (b).

Both achieved interesting results, despite the generated models having reasonable recall; specificity, AUC, and precision yielded acceptable outcomes. Then diameter 4 was selected. Consequently, it was feasible to obtain a suitable model for predicting molecules with anti-tuberculosis activity.

Subsequently, we applied the obtained model to the dataset acquired from the ChEMBL platform,⁴⁰,⁴¹ which had been pre-selected, comprising a total of 83,271 molecules. The outcome revealed a prediction of 810 molecules with anti-tuberculosis activity.

Following this, the 810 molecules underwent virtual screening against the shikimate kinase of M. tuberculosis. The results revealed that approximately 86% (699 molecules) of the total predicted molecules exhibited significant outcomes concerning the target enzyme. Within this molecular subset, approximately 7% (54 molecules) demonstrated an affinity ranging from -9 to -9.8 kcal mol ¹, approximately 46% (323 molecules) achieved a score between -8 and -8.9 kcal mol^-1, and another ca. 46% (322 molecules) yielded a score within the range of -7 to -7.9 kcal mol^-1. The obtained result for the substrate was -6.9 kcal mol^-1, a value defined as the cutoff point.

When analyzing the molecular properties of these 810 molecules, we observed significant differences in the descriptors Log P, MW, NHBA, RC, and TPSA, as shown in Table 2, between those labeled as active and inactive. A TPSA value < 140 Ų indicates good intestinal absorption, while TPSA values < 60 Ų suggest good blood-brain barrier penetration.⁶⁸ Our findings indicate that the set predicted by the model, on average, exhibits favorable values for this descriptor. Moreover, we successfully obtained a group of molecules with interesting structural diversity and distinct molecular properties in each active and inactive group.

Table 2
Difference between the means of the calculated molecular properties

Log Pa	MWb	NHBAb	NHBD	NRB	RCb	TPSAb	Group
2.59	281.97	4.33	2.08	4.15	2.89	73.84	active
2.34	262.37	3.88	2.03	4.63	2.14	62.99	inactive

• ^a
  Tukey’s test, alpha 0.05, p < 0.0048;
• ^b
  Tukey’s test, alpha 0.05, p < 0,001; Log P: partition coefficient between n-octanol: water; MW: molecular mass; NHBA: number of hydrogen bond acceptor atoms; NHBD: number of hydrogen bond donor atoms; NRB: number of rotatable bonds; RC: ring count; TPSA: topological polar surface area.

And subsequently, choose two molecules each from three groups based on the affinity values obtained from the molecular docking calculation (Figure 2). Group A, whose values ranged from -9 to -9.8 kcal mol^-1, included the molecules CHEMBL1229147 and CHEMBL4095667. For group B, the molecules CHEMBL239030 and CHEMBL120640 were selected, with values in the range of -8 to -8.9 kcal mol^-1. Group C comprised molecules CHEMBL273205 and CHEMBL3237485, with values ranging from -7 to -7.9 kcal mol^-1 (molecular structures are provided in the Supplementary Information (SI) section, Figure S1).

Figure 2
Selected molecules, (group C) CHEMBL3237485 = -7.7 kcal mol^-1 and CHEMBL273205 = -7.7 kcal mol^-1; (group B) CHEMBL120640 = -8.2 kcal mol ¹ and CHEMBL239030 = -8.6 kcal mol^-1; (group A) CHEMBL4095667 = -9.5 kcal mol^-1 and CHEMBL1229147 = -9.8 kcal mol^-1.

Based on these results, CHEMBL1229147 exhibited the highest affinity for the target protein. As depicted in Figure 3, it is evident that all molecules interacted in the same region of the protein, which is the region where the substrate of shikimate kinase also interacts.

Figure 3
Visualization of ligands interacting with shikimate kinase, with molecules from group A depicted in blue, molecules from group B in magenta, and molecules from group C in orange. Figure generated using PyMOL software.⁶⁹

The clustering results of the pharmacophores of the selected ligands revealed the presence of four predominant features: donors with a total of 12, acceptors with a total of 18, hydrophobics with a total of 17, and aromatics with a total of 16 (Figure S14, SI section). This suggests the potential types of interactions between the ligand and the receptor. Moreover, based on the data obtained from the molecular property calculations, we can confirm the significance of hydrogen bond acceptor groups.

When analyzing the binding mode of molecules with SK, it was possible to identify the amino acid residues with which the ligands could be interacting. For group C, we identified the residues that may be involved in hydrogen bonding, namely Arg58 and Gly81. It is noteworthy that only the compound CHEMBL273205 interacted with both residues, while the ligand CHEMBL3237485 interacted solely with Arg58. For group B, we identified the residues potentially engaged in hydrogen bonding as Gly80, Ser16, Lys15, Arg136, and Arg58. Compound CHEMBL239030 interacted with all these residues, whereas no hydrogen bonding was observed with the ligand CHEMBL120640. Concerning group A, the residues capable of forming hydrogen bonds were identified as Gly80, Gly12, Ser16, Arg136, and Arg117. Compound CHEMBL4095667 interacted with most of these residues, with the exception of Gly80, while ligand CHEMBL1229147 interacted only with Gly12, Ser16, and Gly80. It is noteworthy that the most frequent amino acid residues (AA) were Arg58, Gly80, Ser16, Arg136. When comparing these residues with those interacting via hydrogen bonding with shikimate (Arg58, Gly80, Arg136), in addition to those involved in van der Waals interactions, it suggests that these compounds function as competitive inhibitors with the substrate. Similar results were observed by Aoki and co workers,⁷⁰ where they found the same pattern for the inhibitor studied in their work.

The SK comprises three main domains: a flexible L-dopa-induced dyskinesia (LID) domain (amino acid residues 112 124), which displays an opening and closing characteristic; an SB (substrate binding) binding domain, responsible for substrate binding in SK; and a more rigid CORE domain composed of α-β-α motifs. The binding pocket of SK contains several conserved residues, such as Asp33, Arg57, Arg116, and Arg132, which have been identified as crucial for catalysis, as mutations in these amino acids lead to inactive enzymes.¹³,¹⁴

When assessing the overall behavior of the shikimate kinase structure with its substrate throughout molecular dynamics, as depicted by the RMSD in Figure S2a (SI section), the system exhibits pronounced fluctuations at the onset of the dynamics, followed by stabilization. This fluctuation arises from the high flexibility of the LID domain. The data obtained from the RMSF calculations highlight the mobility in this region as well, as illustrated in Figure S2b (SI section). During the process of substrate binding and unbinding, the substrate dissociates from the protein and subsequently rebinds, ultimately remaining in the binding site. This behavior is elucidated by the data obtained from the calculation of the average distance between the center of mass of the protein and the ligand, as presented in Figure S2d (SI section).

The RMSD results of compound CHEMBL1229147, the molecule with the strongest interaction as predicted by the docking results as illustrated in Figure 4a, exhibit structural fluctuations attributed to the flexibility of the LID domain. Despite these fluctuations, the system remained stable, as further confirmed by the analysis of the RG (see SI section), which indicated that the system maintained stability with small amplitude oscillations. The RMSF calculations also revealed significant fluctuations in the LID domain region (Figure 4b).

Figure 4
(a) RMSD from the MD simulation of the SK complexed with the ligand CHEMBL1229147, with the moving average represented in red. (b) RMSF from the MD simulation of the SK complexed with the ligand CHEMBL1229147. Noteworthy is the fluctuation in the LID domain (112-124). (c) Number of hydrogen bonds between SK and compound CHEMBL1229147 during MD. (d) Average distance between the center of mass of the protein and the ligand CHEMBL1229147, with the moving average represented in red.

During the molecular dynamics simulations, the maximum number of hydrogen bonds formed by the SK-ligand binary complex was three. The interaction remained consistently stable, maintaining a consistent number of one bond throughout the simulation period, as depicted in Figure 4c. While the ligand remained complexed with the protein during the dynamics, in the last ca. 25 ns, the complex exhibited periodic transitions between the bound and unbound states. This behavior is further elucidated by analyzing the data obtained from calculating the average distance between the center of mass of the protein and the ligand was 1.7 ± 0.2 nm, as depicted in Figure 4d.

The results of the molecular dynamics of the compound CHEMBL4095667, another molecule belonging to the highest affinity group, derived from the calculation of RMSD (Figure 5a) and RG (see Figure S13, SI section), indicate that the system remained stable during the simulation. Similar to other studied systems, this compound also exhibited fluctuations in the LID domain region, as evidenced by the data obtained from the RMSF calculation (Figure 5b).

Figure 5
(a) RMSD from the MD simulation of the SK complexed with the ligand CHEMBL4095667, with the moving average represented in red. (b) RMSF from the MD simulation of the SK complexed with the ligand CHEMBL4095667. (c) Number of hydrogen bonds between SK and compound CHEMBL4095667 during MD. (d) Average distance between the center of mass of the protein and the ligand CHEMBL4095667, with the moving average represented in red.

Throughout the molecular dynamics, the SK-ligand binary complex formed a maximum of eight hydrogen bonds, and this change remained consistently stable, maintaining an average of three bonds (Figure 5c). The ligand remained bound to the protein during the simulation. Analysis of the data obtained from the calculation of the average distance between the center of mass of the protein and the ligand was 0.75 ± 0.06 nm, which highlights this behavior (Figure 5d).

For compounds belonging to group B, ligand CHEMBL120640 (see SI section) remained bound to the protein for a longer period during the initial 50 ns, subsequently exhibiting similar dissociation and association behavior. Throughout the molecular dynamics simulation, it achieved a maximum of seven hydrogen bonds, but maintained an average of two bonds along the trajectory. In contrast, ligand CHEMBL239030 remained bound during the initial 25 ns and between 70 to 80 ns, displaying coupled/uncoupled behavior outside these intervals. It formed a maximum of six hydrogen bonds and maintained an average of two bonds throughout the trajectory. Analysis of the data obtained from the calculation of the average distance between the center of mass of the protein and the ligand was 0.9 ± 0.06 nm (see SI section for RMSD, RG, number of hydrogen bonds, RMSF, and average distance figures).

For group C ligands, there was greater oscillation in coupled/uncoupled states, and the protein retained its opening and closing behavior due to the movement of the LID domain. Both molecules, over the course of the trajectory, formed an average of one hydrogen bond and had a maximum of four bonds (see SI section for RMSD, RG, number of hydrogen bonds, RMSF, and average distance figures).

The estimate binding free energy using the MM/PBSA method for compound CHEMBL1229147 was -67.3 kJ mol^-1, whereas for compound CHEMBL4095667, it was -82.2 kJ mol^-1, and for ligand CHEMBL120640, it was -84.9 kJ mol^-1. The compounds of group C exhibited values ranging from -42 to -62 kJ mol^-1 (Table 3). Additionally, an analysis of residue contributions to the binding energy was conducted for the three compounds exhibiting the highest affinity with SK. Despite not aligning with the order derived from docking calculations, a consistent pattern of activity persisted, with compounds from group A, B, and C exhibiting, towards the receptor, a superior affinity than the substrate, which registered a value of -25 kJ mol^-1.

Table 3
MM/PBSA analysis of protein-ligand complexes

Molecule	VDW / (kJ mol-1)	EEL / (kJ mol-1)	PSE / (kJ mol-1)	SASA / (kJ mol-1)	∆Gbinding / (kJ mol-1)
CHEMBL1229147	-98.4	-23.7	65.1	-10.4	-67.3
CHEMBL4095667	-136.6	-65	133.1	-13.8	-82.2
CHEMBL120640	-129.4	-72.7	131.2	-14	-84.9
CHEMBL239030	-112.1	-33.4	95.6	-12.3	-62.2
CHEMBL273205	-80.5	-25.1	72.3	-9.1	-42.4
CHEMBL3237485	-97.5	-21.6	74.3	-11.7	-56.5
Substrate	-64	-88.4	134.9	-8.4	-25.8

According to a MM/PBSA residue decomposition analysis for compound CHEMBL1229147, residues Ile48, Arg117, Phe49, Ser44, and Thr115 were identified as the top five most negatively contributing residues. For ligand CHEMBL4095667, these residues were Val116, Arg117, Leu10, Leu119, and Pro118. Similarly, for compound CHEMBL120640, the identified residues were Val116, Arg117, Leu10, Leu119, and Pro118. When assessing the binding modes for these compounds, it was observed that certain residues engage in hydrophobic interactions. This implies that not only hydrogen bonds play a crucial role but also weaker interactions such as van der Waals (VDW) interactions. Notably, VDW interactions made the most substantial contributions to the stability of the binary complex, as evidenced by the contributions presented in Table 3. In terms of magnitude, VDW interactions exerted the most significant influence on lowering the calculated binding energy variation in the MM/PBSA analysis.

Moreover, it is noteworthy that ligand CHEMBL120640 exhibited higher electrostatic contributions when compared to both CHEMBL1229147 and CHEMBL4095667, but the values found for VWD were also high in terms of magnitude. This observation underscores the importance of considering not only hydrogen bonding but also weaker interactions, particularly VDW interactions, in understanding the overall stability of the ligand-receptor complexes. During the molecular dynamics simulations, the CHEMBL1229147 ligand presented an average distance value of 1.7 nm, a higher value when compared to the others, that is, this complex was not as stable when compared to the other two ligands, CHEMBL120640 with average distance value of 0.9 nm and CHEMBL4095667 with average distance value of 0.75 nm.

We used the final structure of the SK obtained in the molecular dynamics simulations as receptor for a docking using the same methodology as before. Through the results of docking calculations on the relaxed structure, molecules were reclassified based on the new affinity values presented in Table 4. The pattern concerning molecules with higher affinity towards the substrate remained consistent. Ligand CHEMBL1229147 exhibited the most negative affinity, as observed in the previous analysis, while compound CHEMBL4095667 ranked third. Interaction residue results for CHEMBL1229147 included Ser16, Gly14, Gly9, and Lys15, involving hydrogen bonding for the first three interactions, and for the last interaction which is of the π-alkyl type. For CHEMBL4095667, Tyr140 participated in hydrogen bonding. In the case of CHEMBL239030, the interacting residues were Gly112, Gly12, Ser77, Asp32, and Ser16.

Table 4
Results of global descriptors and docking calculations on the relaxed structure

Molecule	ECP / eV	Hardness / eV	Electrophilicity / eV	Docking score / (kcal mol-1)	Group
CHEMBL1229147	-4.43	3.74	2.63	-8.2	A
CHEMBL239030	-4.83	3.90	2.98	-7.7	A
CHEMBL4095667	-4.69	3.79	2.90	-7.5	A
CHEMBL120640	-4.24	4.06	2.21	-6.7	B
CHEMBL273205	-4.70	3.90	2.83	-6.4	B
CHEMBL3237485	-4.01	3.83	2.10	-6.7	B
SKH	-5.88	5.01	3.44	-5.5	S

The PCA analysis illustrated how these groups of molecules cluster based on their electronic profile (Figure 6), with an explained variance of 99.9%. The collective description of this profile, encompassing electronic properties such as ECP, hardness, and electrophilicity (Table 4), indicates not only electronic similarities but also structural resemblance. This justifies why compound CHEMBL273205 is positioned closer to those classified as group A.

Figure 6
Results of the analysis of the main components of the descriptors: electronic chemical potential, hardness, and electrophilicity.

When we analyze the global hardness as shown in Table 4, it is possible to observe that those who had a greater contribution in the electrostatic term of the estimation of the variation of the binding free energy, also had a greater value of hardness, and the substrate that had the highest value of these parameters, indicates how Coulomb forces are important in the binding process with shikimate kinase.

When analyzing local hardness descriptors, we were able to identify residues with the most significant values: Lys15, Arg58, Arg136, and Tyr140, interestingly, some of these had been mapped as possible participants in hydrogen bonds. For the Fukui function, we identified residues Arg136, Tyr140, and Lys15.

Among the three primary compounds examined, CHEMBL4095667 exhibited a π-alkyl interaction with residue Arg136. This interaction was the most pronounced as indicated by the Fukui function, supporting the likelihood of such an interaction. Residues Arg117 and Leu10 were identified as major contributors to the binding energy of the complex and also displayed high intensity according to the Fukui function. Furthermore, Arg136, Arg117, and Tyr140 demonstrated elevated local hardness. The latter was previously identified as a potential hydrogen bond site.

In CHEMBL1229147, a π-alkyl interaction was observed between Lys 15 and the ligand. Based on local hardness and Fukui function reactivity descriptors, this interaction was determined to be one of the strongest. Residues Arg 117 and Leu 10 were identified as major contributors to the binding energy of the complex. Furthermore, these residues also exhibited high intensity according to the Fukui function and local hardness descriptors.

In CHEMBL120640, a π-alkyl interaction was observed between the ligand and Arg136. This interaction was the strongest as indicated by the Fukui function. Arg117 was identified as a major contributor to the binding energy of the complex and also exhibited high intensity according to the local hardness and Fukui function descriptors. Furthermore, a π-cation interaction was observed with Lys15, which also displayed high intensity based on the Fukui function and local hardness.

During the absorption, distribution, metabolism, excretion and toxicity (ADMET) analysis, we evaluated the potential for blood-brain barrier penetration, P-glycoprotein interaction, gastrointestinal absorption, water solubility, toxicity, and CYP inhibition. Interestingly, all molecules showed good gastrointestinal absorption. However, among the three compounds, CHEMBL4095667 exhibited the best results (Table 5), being the most soluble and having the highest median lethal dose (LD₅₀), thus being the least toxic compared to the others.

Table 5
ADMET analysis results for the three main ligands

Molecule	WS	GI	BBB	Pgp S	CI	Toxicity (LD50) / (mg kg-1)
CHEMBL1229147	MS	high	yes	yes	yes	100
CHEMBL4095667	S	high	yes	yes	yes	1500
CHEMBL120640	MS	high	yes	no	yes	250

Conclusions

In conclusion, the process of developing or discovering a new drug involves numerous variables, and this procedure is inherently slow, not always culminating in success. The technological advancements of recent decades have propelled scientific research towards approaches enabling rapid analyses and the handling of vast datasets, among other capabilities. Within this context, machine learning has demonstrated significant potential as a tool for analyzing medicinal chemistry data. In the current study, the identified model predicted 810 molecules with anti-tuberculosis activity. One of the challenges encountered in developing the deep learning model was minimizing overfitting. Both the model implementation and hyperparameter choices were carefully designed to mitigate this effect. Additionally, we had to adapt the network architecture. Given the relatively small dataset, a deep network with many hidden layers would have exacerbated overfitting. Therefore, we opted for a neural network with a shallower architecture to minimize overfitting in this study. When combined with virtual screening aimed at identifying molecules with an affinity for MTB shikimate kinase, it was observed that 86% of these model-predicted molecules exhibited an affinity for the target protein. Among this subset, 54 molecules demonstrated a docking score ranging from -9 to -9.8 kcal mol^-1. The group of molecules selected for molecular dynamics calculations, MM/PBSA, and global descriptors exhibited a similar binding mode, occupying the substrate-binding region, suggesting a mechanism of competitive inhibition. However, the need for future studies to conduct experimental tests to validate the findings is evident. Furthermore, it was possible to ascertain that the dataset with higher affinity shared a similar electronic profile. The molecules with the lowest ∆G_binding values were identified as CHEMBL1229147, CHEMBL4095667, and CHEMBL120640, which can subsequently be selected as prototypes for the development of a new shikimate kinase inhibitor.

Supplementary Information

Supplementary information is available free of charge at http://jbcs.sbq.org.br as PDF file.

Acknowledgments

We thank Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPQ) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) for their support.

References

• ¹ Basso, L. A.; da Silva, L. H. P.; Fett-Neto, A. G.; de Azevedo, W. F.; Moreira, I. D. S.; Palma, M. S.; Calixto, J. B.; Astolfi Filho, S.; dos Santos, R. R.; Soares, M. B. P.; Santos, D. S.; Mem. Inst. Oswaldo Cruz 2005, 100, 475. [Crossref]
  » Crossref
• ² Andrade, C. H.; Pasqualoto, K. F. M.; Zaim, M. H.; Ferreira, E. I.; Rev. Bras. Cienc. Farm. 2008, 44, 167. [Crossref]
  » Crossref
• ³ Seung, K. J.; Keshavjee, S.; Rich, M. L.; Cold Spring Harbor Perspect. Med. 2015, 5, a017863. [Crossref]
  » Crossref
• ⁴ Pieters, J.; Cell Host Microbe 2008, 3, 399. [Crossref]
  » Crossref
• ⁵ Merle, C. S.; Fielding, K.; Sow, O. B.; Gninafon, M.; Lo, M. B.; Mthiyane, T.; Odhiambo, J.; Amukoye, E.; Bah, B.; Kassa, F.; N’Diaye, A.; Rustomjee, R.; de Jong, B. C.; Horton, J.; Perronne, C.; Sismanidis, C.; Lapujade, O.; Olliaro, P. L.; Lienhardt, C.; N. Engl. J. Med. 2014, 371, 1588. [Crossref]
  » Crossref
• ⁶ Motta, I.; Boeree, M.; Chesov, D.; Dheda, K.; Günther, G.; Horsburgh, C. R.; Kherabi, Y.; Lange, C.; Lienhardt, C.; McIlleron, H. M.; Paton, N. I.; Stagg, H. R.; Thwaites, G.; Udwadia, Z.; Van Crevel, R.; Velásquez, G. E.; Wilkinson, R. J.; Guglielmetti, L.; Motta, I.; Kherabi, Y.; Van Crevel, R.; Guglielmetti, L.; Clin. Microbiol. Infect. 2024, 30, 1107. [Crossref]
  » Crossref
• ⁷ Brake, L. H. M.; de Jager, V.; Narunsky, K.; Vanker, N.; Svensson, E. M.; Phillips, P. P. J.; Gillespie, S. H.; Heinrich, N.; Hoelscher, M.; Dawson, R.; Diacon, A. H.; Aarnoutse, R. E.; Boeree, M. J.; Eur. Respir. J. 2021, 58, 2000955. [Crossref]
  » Crossref
• ⁸ World Health Organization (WHO), WHO Consolidated Guidelines on Tuberculosis. Module 4: Treatment-Drug-Resistant Tuberculosis Treatment, 2022 update; WHO: Geneva, 2022. [Link] accessed in January 2025
  » Link
• ⁹ Dheda, K.; Gumbo, T.; Gandhi, N. R.; Murray, M.; Theron, G.; Udwadia, Z.; Migliori, G. B.; Warren, R.; Lancet Respir. Med. 2014, 2, 321. [Crossref]
  » Crossref
• ¹⁰ Padmapriyadarsini, C.; Vohra, V.; Bhatnagar, A.; Solanki, R.; Sridhar, R.; Anande, L.; Muthuvijaylakshmi, M.; Rana, M. B.; Jeyadeepa, B.; Taneja, G.; Balaji, S.; Shah, P.; Saravanan, N.; Chavan, V.; Kumar, H.; Ponnuraja, C.; Livchits, V.; Bahl, M.; Alavadi, U.; Sachdeva, K. S.; Swaminathan, S.; Padmapriyadarsini, C.; Jeyadeepa, B.; Lakshana; Akbar, N.; Arulraj, E.; Karthikeyan; Muthukumar; Tamizharasan; Balaji, S.; Shivakumar, S.; Muthuvijayalakshmi, M.; Gayathri; Ponnuraja, C.; Kumar, H.; Saravanan, N.; Sridhar, R.; Kumar, R.; Ramesh; Vohra, V.; Rana, M. B.; Singla, N.; Myneedu, V. P.; Lawrence, A.; Kushwaha, D.; Shivam, D. K.; Sarin, R.; Bhatnagar, A. K.; Taneja, G.; Rawat, A.; Haniff, M.; Rahul; Rai, P.; Saini, S.; Mathur, K. K.; Solanki, R. N.; Patel, P. G.; Prajapati, V.; Parmar, B.; Wadkar, K.; Shah, P. L.; Parmar, S.; Vyas, P.; Mistri, K.; Anade, L.; Chavan, V.; Bhui, N. K.; Tipre, P.; Shah, D.; Patwa, S. K.; Nhavakar, A.; Brito, A.; Keny, K.; Karanjkar, V.; Pal, K.; Godam, K.; Huje, M.; Ghadge, S.; Udmalle, M.; Posture, V. V.; Bansode, J.; Bhal, M.; Ranjan; Pillai, D.; Semwal, S.; Livchits, S. L. V.; Alavadi, U.; Swamikan, R.; Nasubo, D. N.; Parmar, M.; Sahu, S.; Mukadi, Y.; Swaminathan, S.; Clin. Infect. Dis. 2023, 76, e938. [Crossref]
  » Crossref
• ¹¹ Mok, J.; Lee, M.; Kim, D. K.; Kim, J. S.; Jhun, B. W.; Jo, K. W.; Jeon, D.; Lee, T.; Lee, J. Y.; Park, J. S.; Lee, S. H.; Kang, Y. A.; Lee, J. K.; Kwak, N.; Ahn, J. H.; Shim, T. S.; Kim, S. Y.; Kim, S.; Kim, K.; Seok, K. H.; Yoon, S.; Kim, Y. R.; Kim, J.; Yim, D.; Hahn, S.; Cho, S. N.; Yim, J. J.; Lancet 2022, 400, 1522. [Crossref]
  » Crossref
• ¹² Rosado, L. A.; Vasconcelos, I. B.; Palma, M. S.; Frappier, V.; Najmanovich, R. J.; Santos, D. S.; Basso, L. A.; PLoS One 2013, 8, e61918. [Crossref]
  » Crossref
• ¹³ Dias, M. V. B.; Faím, L. M.; Vasconcelos, I. B.; de Oliveira, J. S.; Basso, L. A.; Santos, D. S.; de Azevedo, W. F.; Acta Crystallogr., Sect. F: Struct. Biol. Cryst. Commun. 2007, 63, 1. [Crossref]
  » Crossref
• ¹⁴ Ojeda-May, P.; Biophysica 2022, 2, 194. [Crossref]
  » Crossref
• ¹⁵ Lipinski, C. F.; Maltarollo, V. G.; Oliveira, P. R.; da Silva, A. B. F.; Honorio, K. M.; Front. Robot AI 2019, 6, 108. [Crossref]
  » Crossref
• ¹⁶ Mehra, R.; Rajput, V. S.; Gupta, M.; Chib, R.; Kumar, A.; Wazir, P.; Khan, I. A.; Nargotra, A.; J. Chem. Inf. Model. 2016, 56, 930. [Crossref]
  » Crossref
• ¹⁷ Gordon, S.; Simithy, J.; Goodwin, D. C.; Calderón, A. I.; Perspect. Med. Chem. 2015, 7, S13212. [Crossref]
  » Crossref
• ¹⁸ Chagaleti, B. K.; Reddy, M. B. R.; Saravanan, V.; Shanthakumar, B.; Pryia, D.; Kumar, S. K.; Kathiravan, M. K.; J. Biomol. Struct. Dyn. 2023, 41, 14582. [Crossref]
  » Crossref
• ¹⁹ Almeida, A. M.; Marchiosi, R.; Abrahão, J.; Constantin, R. P.; dos Santos, W. D.; Ferrarese-Filho, O.; Phytochem. Rev. 2024, 23, 421. [Crossref]
  » Crossref
• ²⁰ Vamathevan, J.; Clark, D.; Czodrowski, P.; Dunham, I.; Ferran, E.; Lee, G.; Li, B.; Madabhushi, A.; Shah, P.; Spitzer, M.; Zhao, S.; Nat. Rev. Drug Discovery 2019, 18, 463. [Crossref]
  » Crossref
• ²¹ Lo, Y. C.; Rensi, S. E.; Torng, W.; Altman, R. B.; Drug Discovery Today 2018, 23, 1538. [Crossref]
  » Crossref
• ²² Gupta, R.; Srivastava, D.; Sahu, M.; Tiwari, S.; Ambasta, R. K.; Kumar, P.; Mol. Diversity 2021, 25, 1315. [Crossref]
  » Crossref
• ²³ Jiang, T.; Gradus, J. L.; Rosellini, A. J.; Behav. Ther. 2020, 51, 675. [Crossref]
  » Crossref
• ²⁴ Haq, A. U.; Li, J. P.; Saboor, A.; Khan, J.; Wali, S.; Ahmad, S.; Ali, A.; Khan, G. A.; Zhou, W.; IEEE Access 2021, 9, 22090. [Crossref]
  » Crossref
• ²⁵ Patel, L.; Shukla, T.; Huang, X.; Ussery, D. W.; Wang, S.; Molecules 2020, 25, 5277. [Crossref]
  » Crossref
• ²⁶ Haykin, S. S.; Haykin, S. S.; Neural Networks and Learning Machines, 3^rd ed.; Prentice Hall: New York, 2009.
• ²⁷ Xu, Y.; Yao, H.; Lin, K.; Expert Opin. Drug Discovery 2018, 13, 1091. [Crossref]
  » Crossref
• ²⁸ Bre, F.; Gimenez, J. M.; Fachinotti, V. D.; Energy Build. 2018, 158, 1429. [Crossref]
  » Crossref
• ²⁹ Lionta, E.; Spyrou, G.; Vassilatis, D.; Cournia, Z.; Curr. Top. Med. Chem. 2014, 14, 1923. [Crossref]
  » Crossref
• ³⁰ Cheng, T.; Li, Q.; Zhou, Z.; Wang, Y.; Bryant, S. H.; AAPS J. 2012, 14, 133. [Crossref]
  » Crossref
• ³¹ Kitchen, D. B.; Decornez, H.; Furr, J. R.; Bajorath, J.; Nat. Rev. Drug Discovery 2004, 3, 935. [Crossref]
  » Crossref
• ³² Namba, A. M.; da Silva, V. B.; da Silva, C. H. T. P.; Ecletica Quim. 2008, 33, 13. [Crossref]
  » Crossref
• ³³ Rapaport, D. C.; The Art of Molecular Dynamics Simulation, 2^nd ed.; Cambridge University Press: Cambridge, 2004. [Crossref]
  » Crossref
• ³⁴ Wang, C.; Greene, D.; Xiao, L.; Qi, R.; Luo, R.; Front. Mol. Biosci. 2018, 4, 87. [Crossref]
  » Crossref
• ³⁵ Genheden, S.; Ryde, U.; Expert Opin. Drug Discovery 2015, 10, 449. [Crossref]
  » Crossref
• ³⁶ Poli, G.; Granchi, C.; Rizzolio, F.; Tuccinardi, T.; Molecules 2020, 25, 1971. [Crossref]
  » Crossref
• ³⁷ Geerlings, P.; De Proft, F.; Langenaeker, W.; Chem. Rev. 2003, 103, 1793. [Crossref]
  » Crossref
• ³⁸ Torrent-Sucarrat, M.; De Proft, F.; Geerlings, P.; Ayers, P. W.; Chem. - Eur. J. 2008, 14, 8652. [Crossref]
  » Crossref
• ³⁹ Rocha-Santos, A.; Chaves, E. J. F.; Grillo, I. B.; de Freitas, A. S.; Araújo, D. A. M.; Rocha, G. B.; ACS Omega 2021, 6, 8764. [Crossref]
  » Crossref
• ⁴⁰ Zdrazil, B.; Felix, E.; Hunter, F.; Manners, E. J.; Blackshaw, J.; Corbett, S.; de Veij, M.; Ioannidis, H.; Lopez, D. M.; Mosquera, J. F.; Magarinos, M. P.; Bosc, N.; Arcila, R.; Kizilören, T.; Gaulton, A.; Bento, A. P.; Adasme, M. F.; Monecke, P.; Landrum, G. A.; Leach, A. R.; Nucleic Acids Res. 2024, 52, D1180. [Crossref]
  » Crossref
• ⁴¹ Davies, M.; Nowotka, M.; Papadatos, G.; Dedman, N.; Gaulton, A.; Atkinson, F.; Bellis, L.; Overington, J. P.; Nucleic Acids Res. 2015, 43, W612. [Crossref]
  » Crossref
• ⁴² Jhoti, H.; Williams, G.; Rees, D. C.; Murray, C. W.; Nat. Rev. Drug Discovery 2013, 12, 644. [Crossref]
  » Crossref
• ⁴³ Ekins, S.; Bradford, J.; Dole, K.; Spektor, A.; Gregory, K.; Blondeau, D.; Hohman, M.; Bunin, B. A.; Mol. BioSyst. 2010, 6, 840. [Crossref]
  » Crossref
• ⁴⁴ Rogers, D.; Hahn, M.; J. Chem. Inf. Model. 2010, 50, 742. [Crossref]
  » Crossref
• ⁴⁵ Morgan, H. L.; J. Chem. Doc. 1965, 5, 107. [Crossref]
  » Crossref
• ⁴⁶ Landrum, G.; Tosco, P.; Kelley, B.; Rodriguez, R.; Cosgrove, D.; Vianello, R.; Gedeck, P.; Gareth, J.; Schneider, N.; Kawashima, E.; Nealschneider, D.; Dalke, A.; Swain, M.; Cole, B.; Turk, S.; Savelev, A.; Vaucher, A.; Wójcikowski, M.; Take, I.; Scalfani, V. F.; Walker, R.; Probst, D.; Ujihara, K.; Godin, G.; Lehtivarjo, J.; Bérenger, F.; Bisson, J.; Zenodo 2024 [Crossref]
  » Crossref
• ⁴⁷ Abadi, M.; Agarwal, A.; Barham, P.; Brevdo, E.; Chen, Z.; Citro, C.; Corrado, G. S.;Davis, A.; Dean, J.; Devin, M.; Ghemawat, S.;Goodfellow, I.; Harp, A.; Geoffrey Irving, G.; Isard, M.; Jozefowicz, R.; Jia, Y.; Kaiser, L.; Kudlur, M.; Levenberg, J.; Mané, D.; Schuster, M.; Monga, R.; Moore, S.; Murray, D.; Olah, C.; Shlens, J.; Steiner, B.; Sutskever, I.; Talwar, K.; Tucker, P.; Vanhoucke, V.; Vasudevan, V.; Viégas, F.; Vinyals, O.; Warden, P.; Wattenberg, M.; Wicke, M.; Yu, Y.; Zheng, X.; arXiv 2015 [Link] accessed in January 2025
  » Link
• ⁴⁸ Kingma, D. P.; arXiv 2014 [Link] accessed in January 2025
  » Link
• ⁴⁹ Trott, O.; Olson, A. J.; J. Comput. Chem. 2010, 31, 455. [Crossref]
  » Crossref
• ⁵⁰ Trott, O.; AutoDock Vina, version 1.2; Molecular Graphics Lab, The Scripps Research Institute, San Diego, USA, 2021.
• ⁵¹ Hartmann, M. D.; Bourenkov, G. P.; Oberschall, A.; Strizhov, N.; Bartunik, H. D.; J. Mol. Biol. 2006, 364, 411. [Crossref]
  » Crossref
• ⁵² Wallace, A. C.; Laskowski, R. A.; Thornton, J. M.; Protein Eng., Des. Sel. 1995, 8, 127. [Crossref]
  » Crossref
• ⁵³ SwissADME; Molecular Modeling Group, Swiss Institute of Bioinformatics, Switzerland, 2017; Daina, A.; Michielin, O.; Zoete, V.; Sci. Rep. 2017, 7, 42717. [Crossref]
  » Crossref
• ⁵⁴ Banerjee, P.; Eckert, A. O.; Schrey, A. K.; Preissner, R.; Nucleic Acids Res. 2018, 46, W257. [Crossref]
  » Crossref
• ⁵⁵ da Silva, A. W. S.; Vranken, W. F.; BMC Res. Notes 2012, 5, 367. [Crossref]
  » Crossref
• ⁵⁶ Abraham, M. J.; Van Der Spoel, D.; Lindahl, E.; Hess, B.; GROMACS, version 4.5; Stockholm University, Stockholm, Sweden, 2010; Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E.; Bioinformatics 2013, 29, 845. [Crossref]
  » Crossref
• ⁵⁷ Jorgensen, W. L.; Madura, J. D.; J. Am. Chem. Soc. 1983, 105, 1407. [Crossref]
  » Crossref
• ⁵⁸ Showalter, S. A.; Brüschweiler, R.; J. Chem. Theory Comput. 2007, 3, 961. [Crossref]
  » Crossref
• ⁵⁹ Bussi, G.; Donadio, D.; Parrinello, M.; J. Chem. Phys. 2007, 126, 014101. [Crossref]
  » Crossref
• ⁶⁰ Parrinello, M.; Rahman, A.; J. Appl. Phys. 1981, 52, 7182. [Crossref]
  » Crossref
• ⁶¹ Turner, P. J.; XMGRACE, version 5.1.19; Center for Coastal and Land-Margin Research, Oregon Graduate Institute of Science and Technology, Beaverton, 2005.
• ⁶² Kumari, R.; Kumar, R.; Lynn, A.; J. Chem. Inf. Model. 2014, 54, 1951. [Crossref]
  » Crossref
• ⁶³ Kumari, R.; g_mmpbsa MM/PBSA method for GROMACS, version 1.6,2014; Open Source Drug Discovery Consortium (OSDD), Nova Delhi, India, 2024.
• ⁶⁴ Grillo, I. B.; Urquiza-Carvalho, G. A.; Rocha, G. B.; J. Chem. Inf. Model. 2020, 60, 5885. [Crossref]
  » Crossref
• ⁶⁵ Grillo, I. B.; Urquiza-Carvalho, G. A.; Bachega, J. F. R.; Rocha, G. B.; J. Chem. Inf. Model. 2020, 60, 578. [Crossref]
  » Crossref
• ⁶⁶ Stewart, J. J. P.; J. Comput. -Aided Mol. Des. 1990, 4, 1 [Crossref]; Stewart, J. J. P.; MOPAC2016, version 23.0.3; Stewart Computational Chemistry, Colorado Springs, CO, USA, 2016.
  » Crossref
• ⁶⁷ BIOVIA, version 4.5; Discovery Studio Visualizer, Dassault Systèmes, San Diego, USA, 2022. [Link] accessed in January 2025
  » Link
• ⁶⁸ Amin, N.; Sharma, R. K.; Katiyar, D.; Kannaujiya, V. K.; Sinha, R. P.; J. Proteins Proteomics 2024, 15, 135. [Crossref]
  » Crossref
• ⁶⁹ The PyMOL Molecular Graphics System; Pymol, version 3; Schrödinger, LLC, USA, 2024
• ⁷⁰ Kawamoto, S.; Hori, C.; Taniguchi, H.; Okubo, S.; Aoki, S.; Tuberculosis 2023, 141, 102362. [Crossref]
  » Crossref

Edited by

• Editor handled this article:
  Paula Homem de Mello (Associate)

Publication Dates

• Publication in this collection
  17 Feb 2025
• Date of issue
  2025

History

• Received
  12 Aug 2024
• Accepted
  31 Jan 2025