Virtual Screening and Molecular Simulation of Drug Candidates Targeting the Human SLC4A4 Protein for Colorectal Cancer Therapy

2.1. Data Acquisition

To investigate the expression pattern of the gene SLC4A4 in colorectal cancer (CRC), we accessed a public resource, Expression Atlas, established by EMBL-EBI (https://www.ebi.ac.uk/gxa/home) [15]. The gene-specific data were retrieved by searching for the gene symbol (SLC4A4) or the Ensembl gene ID (ENSG00000080493). From the gene overview page, we selected the “Differential expression” tab, which reports changes in gene expression across different experimental conditions and tissues. This tabular result contains major features such as experiment accession E-GEODseries, comparisons of disease versus control, log2FC, and p-value (adjusted p-value) [15]. Then, the findings were reanalyzed to identify fold changes in gene expression. It was performed by evaluating SLC4A4 expression in CRC tissues compared with normal tissues. The analysis was performed to confirm differential expression reliably enabling to analyze the role of SLC4A4 in CRC.

2.2. Protein Preparation

NBCe1, the human SLC4A4, a sodium-coupled acid-base transporter, was prioritized in this study from a pool of several colorectal cancer-associated gene expression data obtained from the Expression Atlas database. Homology modeling was performed to refine the CryoEM structure of the human SLC4A4 protein (PDB ID: 6CAA) [30] and to identify potential structural gaps/ missing amino acid residues for subsequent investigation using the SWISS-MODEL server [31, 32]. In this procedure, the SWISS-MODEL template library (SMTL 2024-03-27 version) was searched for evolutionarily related protein structures as templates. The best template is then selected based on a heuristic selection strategy.

2.3. Screening (Receptor-based)

For receptor-based screening, the 3D structure of the target protein (PDB ID: 6CAA) retrieved from SWISS-MODEL was used to find potential binding pockets on the protein surface. The DrugRep platform, an automated, parameter-free virtual screening server, was used to facilitate this process. DrugRep employs CurPocket, an in-house curvature-based cavity-detection method, to automatically predict potential pockets on the protein surface. This method is efficient for extracting docking boxes for molecular docking. Next, a high-throughput virtual screening was performed on a curated set of 2470 FDA-approved drugs from the DrugBank (version 5.1.7) collection. Docking procedure. The batch docking over the drug library was performed using DrugRep, which provides a procedure for batch docking using the widely used molecular docking tool AutoDock Vina (version 1.1.2) [33]. The platform sets the docking parameters, including the grid center coordinates and the dimensions of the search region, based on the predicted pockets.

The exhaustiveness parameter was set to 8, the default for DrugRep, which balances computational speed and accuracy. Docking results were ranked by Vina binding affinity scores in kcal/mol. Compounds with docking scores lower than –8.0 kcal/mol were considered promising hits for further analysis. This threshold was selected as a stringent cutoff to identify compounds with strong predicted binding affinities, and it is a common practice in virtual screening studies to filter a large library down to a manageable number of lead candidates [34]. The top hits were then passed to a more refined docking step with CBDock2 [35].

2.4. Curation and Preparation of FDA-Approved Ligands for Virtual Screening

The DrugRep cultivated approved drug library refers to a collection of FDA-approved drugs that have been modified for curation and made available in virtual screening campaigns or repurposing studies. This library includes current, up-to-date medications that are already FDA-approved for human use, along with detailed pharmacological profiles and safety data. As part of the high-throughput docking workflow, the online tool DrugRep was used with default parameters to prepare ligands from a curated library of FDA-approved drugs obtained from the DrugBank database (version 5.1.7).This intentional emphasis on a known library of FDA-approved drugs is one of the important advantages of our drug repurposing approach, as it obviates the need for the comprehensive safety and pharmacokinetic studies required for new chemical entities. The compounds were ranked by docking score, with higher scores indicating stronger predicted interactions with the protein's binding site. The three best-scoring candidates were taken forward for further studies to identify the best drug [34].

2.5. Molecular Docking Using CB-DOCK2

To enhance the docking analysis and verify the binding poses of the top 3 drugs (DB01396 [Digitoxigenin], DB00762 [Irinotecan], and DB04868 [nilotinib]), we used the CB-Dock2 web server. CB-Dock2, a two-tier blind docking algorithm with template-based and structure-based blind docking approaches. The server searches for homologous templates with similar protein and ligand structures to facilitate docking if the above first steps are successful. Yet, in this study, with the protein and ligands of interest, the CB-Dock2 server reported that no suitable templates were found in its database [35].

Therefore, all three ligand docking was performed via the ligand-protein blind docking protocol. This procedure, inherited from the original CB-Dock server, is based on the AutoDock Vina engine and guided by the surface curvature of a protein to find potential binding cavities. It further docks molecular structures within the identified pockets for detection to optimize ligand and ligand conformations and evaluate binding affinities. The server generated multiple binding poses for each ligand in some of the predicted cavities. Vina binding affinity scores (in kcal/mol) were used for the final ranking of the complexes, and the best-scored pose for each ligand was considered for subsequent analysis. This protocol successfully identified both the most likely binding sites and the types of interactions involved between ligands and an SLC4A4 protein [35].

2.6. Molecular Dynamics Simulation

The MD simulations were performed with GROMACS 2024.1- dev version (www.gromacs.org) in combination with the CHARMM36-jul2022.ff version of the force field under CHARMM36-jul2022.ff [36]. The initial protein-ligand complex conformation was achieved with CB-DOCK2. Protein and ligand structures were prepared using the pdb2gmx utility, and missing hydrogens were added during the process. Ligand topology and parameters files were created with the latest version of CGenFF through CHARMM-GUI [37]. Using the TIP3P water model, the simulation system was solvated in a cubic water box, with a 10 Å buffer between the protein and the box edges. Charge neutrality was achieved by adding sodium ions. Both systems were treated with periodic boundary conditions for a 50 ns run using the Leapfrog integrator, along with equilibration at NVT and NPT conditions. Electrostatics were calculated by ParticleMesh Ewald (PME). Invite a force-based switching function to capture non-bonded interactions over 10 and 12 Å. To remove the steric clashes and bad contacts, 50,000 steps of steepest descent were run for energy minimization. Structural analyses, such as RMSD and RMSF, were performed using gmxrms and gmxrmsf, respectively. Hydrogen bond analysis was performed with gmxhbond. Additionally, Rg and SASA were calculated using gmx gyrate and gmx sasa, respectively. Simulation data were visualized in PyMOL and plotted using Grace software [38, 39]. Ligand-receptor dynamics were simulated with 50 ns on the GridMarkets platform, which took ~22 hours.

2.7. Study Design

An exploratory quantitative and computational analysis is conducted that combines publicly available transcriptomic data with in silico drug screening tools. The main goal was to investigate the importance of the SLC4A4 gene in colorectal cancer (CRC) and to identify potential FDA-approved drugs targeting its protein product. The study was a structured multi-stage computational process:

1. Gene Expression Analysis – SLC4A4 expression in CRC vs. normal tissue was assessed for differential expression using Expression Atlas public data sets, with log2 fold change and adjusted p-values.
2. Protein Modeling – The Cryo-EM structure of SLC4A4 protein (pdb id: 6caa) was modeled by SWISSMODEL.
3. Virtual Screening and Docking – The modeled protein was subjected to docking with a library of FDA-approved drugs obtained from DrugBank (via DrugRep) using AutoDock Vina and redocked with CB-DOCK2.
4. Molecular Dynamics Simulation – The lead drug candidate (Nilotinib) was further validated for its binding stability using a 50-ns MD simulation with GROMACS, with analysis of RMSD, RMSF, hydrogen bonding, radius of gyration, and SASA metrics.

As this study is observational and computational in nature, it did not involve any live or human subjects. All datasets were retrieved from publicly accessible databases, and the analytical tools employed are open-source or web-based. Such an approach promotes transparency and ensures the reproducibility of computational drug repurposing research.

3.1. Expression of SLC4A4 from Expression Atlas Datasets

Analysis of SLC4A4 Expression from Expression Atlas Datasets. To investigate SLC4A4 expression in colorectal cancer (CRC), we performed a comprehensive analysis using RNA-sequencing (RNA-seq) and array-based data from the Expression Atlas database. When the analysis was restricted to datasets with CRC, significant differences in SLC4A4 were observed. The SLC4A4 expression was downregulated consistently in these datasets. It was identified with the thresholds of Log2 Fold Change (log2FC)<-1 and an adjusted p-value (adj. p-val) ≤0.05. This downregulation of SLC4A4 was observed in many CRC datasets and is shown in Table 1, suggesting a possible role for SLC4A4 in CRC development.

3.2. Molecular Modelling

A 3D model of SLC4A4 with template Q9Y6R1 was built as provided in Table 2. The model was analysed using the Protein Structure and Model Quality Assessment Tool of Swiss-Model Server, which resulted in favorable scores: GME (Global Model Quality Estimate) = 0.71; Ramachandran Favoured (%) = 90.69. Furthermore, the MolProbity Score was 1.64 (Table 3) [32]. These evaluations highlighted the feasibility of employing this model for further structure-based 6CAA analysis and design.

3.3. Molecular Docking

Drug repurposing for approved drugs was done with a DrugRep tool running with default parameters, considering compounds in DrugBank (version 5.1.7). The findings from virtual screening and docking were the identification of a key binding pocket in the protein SLC4A4, which was the major ligand-binding site. The features of this pocket are summarized in Table 4. The pocket (Pocket 1) is a substantial 3676 Å3 at (−1.4, −18.4, −0.4) in the xyz space. It contains critical amino acid residues, including hydrophilic residues such as Aspartic Acid (ASP), Glutamic Acid (GLU) , and Lysine (LYS), and hydrophobic residues such as Phenylalanine (PHE) and Leucine (LEU). This mixed composition means that the pocket is diverse enough to bind multiple ligands via various types of interactions, which is imperative to understand, as many of the best ligands for our four target drugs are characterized by highly positive binding (which is in our favor).

Redocking was carried out with selected compounds using CBDOCK2 to obtain the highest-ranked docking pose for each compound, based on the scoring functions: DB01396 (Digitoxin), DB00762 (Irinotecan), and DB04868 (Nilotinib). Interestingly, Nilotinib had the best docking score, suggesting a higher binding affinity (Table 6). The top three hits, identified based on cavity-based docking, are listed in Table 5.

The docking complex between the modeled 6CAA protein and Nilotinib (Fig. 1) was then used for subsequent molecular dynamics simulation was also performed to assess its stability and binding properties.

3.4. Molecular Dynamics Simulation

For the protein–nilotinib complex, the system showed stable dynamics, and the protein reached a maximum RMSD of ~1.5 nm over a 50-ns simulation (black line in Fig. 2A). The ligand also maintained fluctuations at ~0.30 nm RMSD over the same period (red line in Fig. 2A), indicating a stable positioning within the binding pocket. In general, the RMSD results suggest that the ligand remained stably bound to the protein's active site.

The RMSF analysis showed that the N-terminal part of the protein fluctuated slightly at about 1.0 nm and that the C-terminal was more flexible with its fluctuations up to ~2.0 nm as represented by the black line in the Fig. (2B). In contrast, the ligand was much less flexible, staying within ~0.5 nm during the whole simulation (the red line).

The radius of gyration(Rg) results showed ligand (red line) fluctuated up to 0.5 nm while protein(black line) showed a slow decreasing from 4.0 nm to about 3.8 nm at the end of the simulation. This indicates that the ligand is stably located in the binding site of the protein (Fig. 2C).

Hydrogen bonding analysis showed that the ligand established up to seven hydrogen bonds intermittently at 4, 8, 9, 11, 12, 14, and 16ns during the simulation. At the end of the 50-ns run, three hydrogen bonds survived within the protein-ligand complex, as shown in Fig. (2D).

The Solvent Accessible Surface Area (SASA) for the protein only (black line) and the protein-ligand complex (red line) during the entire MD simulation. The two graphs superimpose closely, indicating little fluctuation, and demonstrating that the system was stable over the entire simulation period (Fig. 2E).

Upon energy minimization, the Nilotinib–protein complex showed a number of significant bonds. At the initial conformation (0 ns), the compound makes an H-bond with Ser38, Arg40, Pro105, Ile132, Gln107, Ser349, and Glu351. Pi-interactions were observed with residues such as Arg40, Ala104, Lys277, Pro279, His351, Thr861, Ser862, Ala863, and Lys869. Additionally, van der Waals interactions involved Val35, Lys37, Tyr39, Arg41, Pro106, Leu108, Lys133, Met348, Thr859, and Glu861 (Fig. 3A).

At the 25 ns mark, Nilotinib still had four hydrogen bonds with Cys38 and Glu351. Pi-interactions were observed with Lys37, Ser38, and His354, and van der Waals interactions with Pro36, Tyr39, Arg40, Ala104, Pro105, Pro106, Gln107, Leu108, Ser349, and Asp350 (Fig. 3B).

At 50 ns (end of simulation), the ligand was hydrogen-bonded to Ser38 and Ser349. Pi-interactions with Pro36, Ala104, Pro105 and His354 and van der Waals interactions with Tyr39, Arg40, Arg41, Pro106, Gln107, Leu108, Glu351, and Lys814 (Fig. 3C).

From 0 ns to 50 ns, Nilotinib remained buried in the protein's binding pocket, suggesting stable binding with diverse favorable interactions throughout the simulation.