A total of 3,620 compounds approved by the FDA were obtained from the DrugBank database (https://go.drugbank.com/) for the docking study (Knox et al., 2024). DrugBank was prioritized for its rigorously curated collection of FDA-approved drugs, which includes detailed pharmacological, chemical, and clinical data. The compounds were preprocessed for their structural minimization, protonation, and appropriate atom assignment types. We utilized Paroxetine as a reference molecule because it is a well-studied SERT inhibitor (Davis et al., 2016). We obtained the SERT crystal structure (PDB ID: 5I6X) from the Protein Data Bank (http://www.rscb.org) at a resolution of 3.14 Å (Coleman et al., 2016). The docking accuracy required standard protein preparation procedures, which started with solvent molecule removal followed by hydrogen atom addition and stereochemical optimization of the structure. MGL AutoDock Tools (Goodsell et al., 2021) were used to optimize the protein structure for additional optimization, culminating in energy minimization to achieve structural stability. The docking procedure was carried out using InstaDock v1.2 (Mohammad et al., 2021) in a blind search space of 86 × 67 × 83 Å grid system that concentrated at coordinates −36.109, −22.026, and 2.803 along the X, Y, and Z-axes.

The docking protocol was validated through a retrospective redocking study, in which co-crystallized Paroxetine was redocked into the SERT binding site. The redocked Paroxetine aligned closely with its original co-crystallized pose (Supplementary Figure S1), demonstrating the accuracy and reliability of the docking protocol in reproducing ligand binding conformations within the SERT pocket. The redocked Paroxetine aligned closely with its original co-crystallized pose, demonstrating an RMSD of 1.286 Å, which confirms the reproducibility and precision of the docking protocol. After the docking procedure, the binding poses were analyzed to identify compounds that showed favorable binding affinities and proper orientation and alignment with the native ligand for further evaluation.

Bioinformatics technologies enable scientists to study drug-target interactions at a detailed level, thus enabling the exploration of new therapeutic ligands. The selected compounds from the molecular docking screening underwent biological evaluation through PASS (Prediction of Activity Spectra for Substances) analysis and drug profiling (Lagunin et al., 2000). The PASS web server (https://www.way2drug.com/passonline/) compares molecular structures with millions of experimentally validated bioactive molecules to determine compound biological activities. The analysis generates two probability values, Pa and Pi, where higher Pa scores indicate greater chances of biological effects occurring. A comprehensive examination of the selected compounds’ binding conformations and molecular interactions was carried out after PASS analysis. The visualization of polar interactions, hydrogen bonding, hydrophobic contacts, and π-π stacking interactions between selected molecules and SERT was performed using PyMOL (DeLano, 2002), followed by Discovery Studio Visualizer (2005) for analyzing detailed binding site molecular interactions.

MD simulation is a leading computational approach for researching biological macromolecules through atomic-scale modeling and conformational movement analysis (Ali et al., 2024). The technique delivers essential details about protein-ligand bindings, protein conformations, and molecular motion patterns, making it vital for contemporary drug discovery operations and biomolecular research (Mohammad et al., 2020). We utilized GROMACS 2022.3β software (Van Der Spoel et al., 2005), operated from an Ubuntu 2024 LTS system. We developed a simulation environment with precision to model molecular interactions correctly. The simulation system was hydrated using TIP3P explicit water molecules (Mark and Nilsson, 2001), without detergents or lipids. The systems were neutralized by adding suitable Na⁺ and Cl⁻ counterions. Both protein and ligand structures utilized the CHARMM36-July2022 force field (Huang et al., 2017) for calculations, while the CGenFF server (Zhu, 2019) was used to generate topology parameters for the ligands to ensure accurate force field results.

The execution of MD followed a sequence of three distinct protocols, which began with energy minimization, followed by equilibration, and finished with production simulation. The steepest descent algorithm ran for 5,000 steps to optimize the initial molecular configuration by removing steric clashes. The equilibration process occurred in two stages through NVT before transitioning to the NPT ensemble simulation. A gradual heating process from 0 K to 300 K was performed while applying weak positional restraints on the protein backbone to stop sudden structural distortions. The systems were equilibrated for 1 ns for a stable production run. The 500 ns production MD simulations enabled a complete structural dynamics analysis between free protein and protein-ligand complexes. The simulation trajectories were recorded every 10 ps to observe essential molecular interactions and conformational changes throughout the simulation period.

Principal component analysis (PCA) is a statistical tool that discovers major motion patterns in MD trajectory data (Moradi et al., 2024). PCA of protein-ligand complexes helps to understand and visualize the system’s collective motions by reducing the data’s dimensionality (Yang et al., 2009). It helps identify key conformational changes and visualize the dynamics in a lower-dimensional space. PCA was used to investigate the atomic motion and global conformational dynamics of SERT and its drug-bound complexes. The methodology enabled the detection of necessary protein movements, alongside detecting meaningful functional protein conformations (Yang et al., 2009). The eigenvector coordinates were determined using a positional covariance matrix calculated from MD simulation data. The principal motions emerged when the covariance matrix underwent an orthogonal transformation for diagonalization to create an eigenvalue diagonal matrix.

The principal components (PCs) obtained from eigenvalues and eigenvectors represented the protein-ligand system’s most significant large-scale atomic motions. PC1 and PC2 emerged as the leading factors that explained atomic displacement. The system’s conformational transitions received deeper analysis through free energy landscape (FEL) generation. We can study protein stability and conformational states through FEL analysis by exploring free energy changes across principal components (Abdelsattar et al., 2021). This method demonstrated the capability to show stable conformational states as energy minima while identifying possible transition states involved in ligand-induced conformational changes.

Molecular docking is a standard computational approach that predicts how ligands orient and fit inside the target protein’s active sites (Naqvi et al., 2018). This methodology reveals the assessment of binding strength, interaction patterns, and modulating potential, enhancing its utility for pharmaceutical research. We used a set of 3620 FDA-approved drug molecules from the DrugBank database for virtual screening against SERT. We utilized Paroxetine, a well-known SERT inhibitor and FDA-approved SSRI, to provide reference standards for binding affinity and interaction profile assessments. The docking analysis identified the top 10 compounds by evaluating their binding energy scores, which indicate how strongly the compounds bind to SERT (Table 1). The selected molecules showed appreciable binding affinities through docking scores ranging from −10.5 to −11.2 kcal/mol, while Paroxetine exhibited a binding affinity of −8.0 kcal/mol. The strength of ligand-protein interactions during molecular docking increases as the binding energy values decrease, indicating better potential for effective inhibition (Sousa et al., 2006). These compounds show significant potential to function as competitive SERT inhibitors because they demonstrate stronger binding affinity than Paroxetine.

The PASS analysis represents a commonly utilized computational method that uses structural properties to forecast biological activity profiles of molecules (Lagunin et al., 2000). The server evaluates molecular structures against numerous bioactive compounds in databases to predict possible pharmacological results. We utilized PASS analysis to determine biological activity properties of the top 10 molecules selected from molecular docking screening against SERT. The study of ten molecules made Flunarizine stand out as the best candidate based on its drug profile associated with depression, migraine, and neuroprotective activities. The PASS results showed that Flunarizine possesses a strong potential for treating neurodegenerative and psychiatric disorders, which supports its development as an SERT inhibitor (Table 2). The PASS predictive tool generates two probability values known as Pa and Pi to assess activity potential. A compound demonstrates promising biological effectiveness if its Pa score exceeds Pi. The antidepressant and neurological disorder treatment capability of Flunarizine was confirmed through Pa values spanning from 0.618 to 0.927. However, further experimental validation of Flunarizine as a repurposed SERT inhibitor should be carried out for MDD and related neuropsychiatric conditions. The reference drug Paroxetine also showed antidepressant, nootropic, and mood disorder treatment with considerably high potential, which validated the results from the PASS predictions.

The interaction analysis of Flunarizine and the reference inhibitor was carried out to understand their binding mechanism with SERT (Figure 1A). The study demonstrated that Flunarizine binds in the SERT active site, showing binding characteristics similar to Paroxetine (Figure 1A). Flunarizine forms various stabilizing chemical bonds with critical SERT binding pocket residues, demonstrating its strong potential to act as a promising inhibitor (Figure 1B). Similarly, Paroxetine also interacts with this binding site of SERT, which showed various common interactions (Figure 1C). The spatial orientation of Flunarizine alongside its binding conformation showed that it binds with critical function-related residues, including Thr439, a serotonin binding site (Yang and Gouaux, 2021), like Paroxetine (Figure 1D). The Flunarizine binding pattern indicates that this compound would function as an effective serotonin-competitive inhibitor of SERT by analogously inhibiting its binding site to current prescribed SSRIs.

Flunarizine establishes a promising complex with SERT through multiple hydrogen bonds combined with π-stacking and hydrophobic interactions, thus potentially increasing its inhibitory effects. Flunarizine forms various other interactions with critical residues more than Paroxetine does, demonstrating a potentially superior binding affinity and more stable interaction with SERT. The residues such as Tyr175, Ser438, Thr439, and Phe341 play essential roles in the inhibitory process by forming crucial bonding interactions that sustain SERT inhibition (Yang and Gouaux, 2021). Flunarizine’s binding interactions surpass Paroxetine’s because it establishes multiple additional interactions. The additional molecular interactions increase the binding energy, indicating that Flunarizine demonstrates a better potential for inhibiting SERT. The wide usage of Paroxetine as an SSRI reveals that Flunarizine reflects the potential to become a repurposed antidepressant drug with new applications.

Further, a detailed analysis showed that the binding mechanism of Flunarizine with critical residues in the SERT binding cavity occurs through multiple interaction types, including hydrogen bonds, halogen bonds, π-interactions, and van der Waals contacts (Figure 2). Flunarizine establishes close interactions with the SERT binding site (Figure 2A). The ligand stability was supported by two interactions: hydrogen bonding with Tyr175 and halogen bonding with Ser438. The π-sigma interaction occurred between Ile179 and Flunarizine, while π-π stacking interactions formed with Tyr176, Trp103, and Phe341. The SERT-Flunarizine complex received stabilization from multiple hydrophobic interactions between Leu99, Trp103, Ile172, Pro403, and Val501 as well as van der Waals interactions that included Asp98, Arg104, Trp182, Phe334, Phe335, Gly338, Phe407, Thr439, Gly442, Val489, Glu493, and Thr497 (Table 3). At the same time, the reference compound Paroxetine also showed a similar interaction profile as Flunarizine (Figure 2B).

It forms hydrogen bonds with Ala96, Ala169, and Ser336 and a single π-sigma interaction with Ile172. Additionally, Phe341 and Ser438 formed π-π stacking interactions, and Tyr95, Ala169, Ile172, and Ala173 formed alkyl/Pi-alkyl interactions. The binding of Paroxetine with SERT formed various van der Waals contacts to Asp98, Tyr176, Phe335, Gly338, Ser439, Gly442, Thr497, and Val501. Flunarizine forms a hydrogen bond with Tyr175 (2.8 Å) and a halogen bond with Ser438 (3.1 Å), residues critical for SERT’s conformational transitions. In contrast, Paroxetine interacts weakly with Ala96 (3.2 Å) and lacks halogen bonding, explaining its lower binding affinity. The analysis revealed that Flunarizine formed various common and more contact points with the serotonin binding site residues of SERT than Paroxetine, thus indicating a potentially more substantial inhibitory potential (Yang and Gouaux, 2021).

The MD simulation method delivers essential information about protein and protein-ligand complexes while demonstrating structural shifts and helping to evaluate binding persistence throughout time (Xue et al., 2022). We performed an all-atom 500 ns MD simulation through GROMACS to analyze SERT and its complexes with SERT-Flunarizine and SERT-Paroxetine. The results of energy calculations showed that SERT potential energy decreased from −1,717,590 kJ/mol for SERT to −1,190,080 kJ/mol for SERT-Flunarizine and −1,190,630 kJ/mol for SERT-Paroxetine. The protein-ligand complex stability improved compared to unbound SERT, as shown by total energy values of −1,380,460, −948,548 kJ/mol, and −949,102 for SERT, SERT-Flunarizine, and SERT-Paroxetine, respectively. The binding stability of Flunarizine indicates it to be higher than Paroxetine’s, thus indicating its potential to act as a competitive SERT inhibitor. The analysis demonstrates how Flunarizine forms a solid complex with SERT, which supports its potential application as a treatment drug for MDD. Various structural deviations and flexibility parameters were analyzed from the MD simulations trajectories, as discussed in the ensuing sections.

PCA was employed to analyze SERT’s molecular motion, conformation flexibility, and ligand complexes (Figure 7). PCA can identify dominant motion patterns by reducing the complexity of molecular dynamics trajectories into principal components (PCs). During the simulation, the first principal components (PC1 and PC2) were employed to plot the vibration space traversed by SERT, SERT-Flunarizine, and SERT-Paroxetine complexes. The 2D subspace PCA plot illustrated the unbound SERT protein traversing a more expansive conformation space, significantly along PC2, reflecting enhanced flexibility and dynamic motion (Figure 7A). The SERT-Flunarizine and SERT-Paroxetine complexes demonstrated more constrained motion distributions along PC1, reflecting reduced conformation flexibility upon ligand binding. The motion ranges observed were PC1: −6.6–2.5 nm and PC2: −2.9–3.8 nm for SERT, PC1: −2.8–4.1 nm and PC2: −4.2–3.9 nm for SERT-Flunarizine, and PC1: −4.8–4.9 nm and PC2: −4.0–3.4 nm for SERT-Paroxetine. To further analyze the conformational state stability, eigenvector plots against time were made for the three systems (Figure 7B). The observation demonstrated SERT with prominent structural fluctuations, while SERT-Flunarizine and SERT-Paroxetine complexes assumed different and stable conformations with time. The stabilization of PC1 and PC2 fluctuations with time reflects the inhibition of excessive protein motion upon ligand binding and the stabilization of SERT into a structured conformation. The results confirm the orientation of a stable conformational state by Flunarizine comparable with Paroxetine, further establishing its potential as a SERT inhibitor.

FEL analysis was performed to investigate further the thermodynamic stability and conformation transitions of SERT and its ligand-bound complexes. FEL maps were obtained from PC1 and PC2 from PCA trajectory data and yielded insights into transition states from metastable towards the native conformation. The 3D Gibbs FEL maps depict changes in energy across conformation states. The areas represented in red correspond to high-energy unstable states, and the regions described in dark blue correspond to low-energy conformations (Figure 8). The unbound SERT protein displayed a broad dark blue basin corresponding to a dynamic and flexible native conformation (Figure 8A). The SERT-Flunarizine complex displayed tighter single dark blue basins corresponding to the constraining of the conformation space by the ligand and stabilization of SERT into singular low-energy conformations (Figure 8B). The SERT-Paroxetine displayed a wider low-energy state confined within multiple basins (Figure 8C). Here, metastable states were also observed in both ligand-bound complexes, corresponding to minor conformation adjustments upon the binding of the ligand. Taken together, the ligand-bound SERT in the presence of Flunarizine exhibited stability in MD simulations and essential dynamics. These findings suggest that Flunarizine can be further explored experimentally as a potential repurposed drug targeting SERT inhibitors for therapeutic development against MDD.

Overall, the results of this study demonstrate that a calcium channel blocker, Flunarizine, showed promise as an antidepressant drug that can be repurposed to target SERT. SERT regulates synaptic serotonin levels, and its dysregulation is implicated in MDD’s emotional and cognitive deficits. Flunarizine is a selective channel entry blocker with calmodulin binding properties and histamine H1 blocking activity (Stubberud et al., 2019). It is well-established that calcium channel blockers have been explored for their potential antidepressant effects (Harrison et al., 2020; Bozorgi et al., 2020). By competitively inhibiting SERT, Flunarizine prolongs serotonin availability, mirroring SSRIs’ mechanism but with enhanced binding stability. The research uses computational screening technology to show how existing FDA-approved drugs can be evaluated for new pharmacological purposes, shortening the drug discovery timeline. While docking identifies initial hits, MD simulations validate their dynamic stability, a critical step absent in standalone docking studies. This synergy mitigates false positives and confirms Flunarizine’s binding persistence. Recent studies corroborate the utility of integrated computational approaches in small-molecule inhibitor discovery (Vikhar Danish Ahmad et al., 2024). Further experimental evaluations on Flunarizine must be performed to verify computational results and assess the drug’s effectiveness for MDD treatment. Future work should prioritize in vitro assays (e.g., serotonin uptake inhibition in neuronal cell lines) and in vivo behavioral studies in rodent models of depression to validate Flunarizine’s antidepressant effects. Additionally, pharmacokinetic studies and safety profiling are essential to assess its suitability for repurposing in MDD. These steps will bridge the gap between computational insights and clinical application.