Alzheimer’s Disease (AD) is a progressive neurodegenerative disorder, causing memory loss and in 60–70% of cases leading to dementia (Caltagirone et al., 1983). The AD pathology is widely believed to be associated with the production of β-amyloid peptide (Aβ), which is responsible for the plaque formations in the brain, disrupting normal neuronal functions (Takahashi et al., 2017). This pathological hallmark of AD might be considered in rational drug design and discovery as a promising therapeutic target to develop effective medication against this disorder. However, there is still a lack of efficient treatment for this disabling and ultimately fatal disease, e.g., donepezil and memantine usually provide at best only temporary and incomplete symptomatic relief (Nie et al., 2011). Moreover, various attempts had been made to use different drug binding sites (subregion-targets) in Aβ that could stop its aggregation and formation of a senile plague (Nie et al., 2011). In particular, tramiprosate (3-Aminopropanesulfonic acid, TRA), a mimic of glycosaminoglycans, targets the HHQK subregion at the N-terminus of Aβ (Nie et al., 2011). The TRA treatment of TgCRND8 mice resulted in a 30% reduction in the brain plaque load and the same decrease in the cerebral levels of soluble and insoluble Aβ (Gervais et al., 2007). Additionally, a dose-dependent 60% reduction of plasma Aβ levels was also observed, suggesting that this influences the Aβ central pool, changing either its efflux or its metabolism in the brain (Gervais et al., 2007). Despite the structural simplicity, high specificity, and excellent in vivo Aβ inhibition properties of TRA, it subsequently failed in the late stages of phase III clinical trial (Rauk, 2008). However, the data obtained from in vitro/vivo experiments and clinical trials could provide valuable evidence that Aβ inhibitors and their scaffolds represent a viable drug designing methodology for AD treatment (Blazer and Neubig, 2009). Indeed, some scaffold-based techniques, such as scaffold hopping, have become a powerful tool to determine the most promising drug-like candidates and already approved drugs in a drug repurposing protocol for AD (Kowal et al., 2019; Ballard et al., 2020). For example, some researchers had already performed in silico screening of a virtual library of a biaryl scaffold-containing compounds to inhibit the primary targets of AD therapeutics, such as acetylcholinesterase, β-secretase, monoamine oxidases, and N-methyl-D-aspartate receptor (Khalid et al., 2018). On the other hand, the scaffold searching of FDA and EMA-approved drug libraries was not previously performed to identify lead candidates for drug repurposing in AD. Moreover, the blood-brain Barrier (BBB) permeation of drug-like molecules is often considered in various pharmacokinetics (PK) studies as a pivotal PK-related descriptor (Shityakov et al., 2012; Carpenter et al., 2014). Therefore, in our study, we wanted to apply the molecular search based on the TRA scaffold to examine a database of FDA and EMA-approved drugs using in silico computational approaches, such as virtual screening, molecular dynamics, and BBB-based descriptor analyses.

The 3D molecular structure (PBD ID: 2NAO) of a disease-relevant Aβ fibril (1–42), containing chains A, B, and C, was downloaded from the Protein Data Bank (Walti et al., 2016) to be used in the study. The Ramachandran plot server (https://zlab.umassmed.edu/bu/rama/index.pl) was implemented for the stereochemical validation of the receptor molecule to investigate the ϕ–ψ dihedral angles in a Ramachandran plot. Altogether, observed statistics showed that 86.54% (90 residues) and 11.54% (12 residues) of all observed residues were in the core and allowed regions. Additionally, no steric clashes were detected in the peptide structure. The TRA molecule and its propanesulfonic scaffold were built by using the MarvinSketch software (ChemAxon, Hungary). A database, containing the FDA and EMA-approved drugs (n = 4,642), was obtained from the DrugCentral 2021 online drug compendium. The Molsoft ICM 3.8-3 scaffold search algorithm was used to filter the database to identify scaffold-containing drugs. All ligands were protonated at pH = 7.4 and T = 310 K using the MOE software. Prior to molecular docking, the CASTp (Computed Atlas of Surface Topology of proteins) algorithm (Naghibzadeh, 2001) was implemented to detect the location of the peptide-ligand binding site with Cartesian coordinates located at the grid center: x = 11.87 Å; y = 18.88 Å; z = −27.75 Å. The AutoDock molecular docking algorithm to calculate binding affinity (ΔG _bind ) was implemented via the Raccoon v1.0 modeling suite to perform virtual screening. The receptor and ligand structure preparations for molecular docking included Gasteiger partial charges assignment and rotatable bonds definition according to the standard protocol published elsewhere (Shityakov et al., 2014). The inhibition constants (Ki) and Ligand Efficiencies (LE) were calculated from the binding energy values as follows: Ki = exp [ Δ G bind RT ] LE = Δ G bind N atm where R (gas constant) is 1.98 cal(mol*K) ⁻¹; T (room temperature) is 298.15 K; N _atm is the number of non-hydrogen atoms in a molecule. AutoDock v.4.2.5.1 was used in the study since its previous version incorrectly calculates part of the intermolecular desolvation energy term (Goodsell et al., 1996). The docking grid with dimension size of 60 × 60 × 60 Å and a grid spacing of 0.375 Å were used in the study. The Glide molecular docking algorithm with the Prime MM-GBSA (generalized Born solvent-accessible surface area) approach to calculate the free energy of binding (ΔG _PR ) was implemented using the default settings, such as the OPLS3e force field, 0.25 of charge cutoff, and 0.8 scaling of the vdW radii of non-polar receptor and ligand atoms. The G-score value, as an empirical scoring function, was used to approximate the ligand binding free energy. All molecular dynamics (MD) simulations were performed using the AMBER 16 package with the FF99SB and GAFF force fields for the Aβ peptide and its ligands (Case et al., 2005). The Antechamber module of AmberTools was employed to calculate the partial charges of the ligands using the semi-empirical AM1-BCC function according to the standard protocol (Marques et al., 2021). The systems were solvated with the TIP3P water models and neutralized by adding the Na + ions using the tLEap input script available from the AmberTools package. Long-range electrostatic interactions were modeled via the particle-mesh Ewald method (Essmann et al., 1995). The SHAKE algorithm (Miyamoto and Kollman, 1992) was applied to constrain the length of covalent bonds, including the hydrogen atoms. Langevin thermostat was implemented to equilibrate the temperature of the system at 310 K. A 2.0-fs time step was used in all of the MD setups. For the minimization and equilibration (NVT and NPT ensembles) phases, 100,000 steps and a 1-ns period were used, respectively. Finally, 100-ns classical MD simulations, with no constraints as NPT ensemble, were performed for each of the peptide-ligand complexes using the molecular mechanics combined with the Poisson–Boltzmann (MM-PBSA) or generalized Born (MM-GBSA) augmented with the hydrophobic solvent-accessible surface area term (Kollman et al., 2000; Shityakov et al., 2017). The MM-PBSA/GBSA solvation models were applied as a post-processing end-state method to calculate the free energies (ΔG _PB and ΔG _GB ) together with the entropies (TΔS) and enthalpies (ΔH) for the analyzed molecules, namely: Δ G PB   =  Δ H PB – TΔS ;  Δ G GB   =  Δ H GB – TΔS

To calculate the blood-brain barrier partitioning coefficients (logBB), the in-house python script, based on the Clark (logBB _cl ) and Rishton (logBB _ri ) linear regression models, was ran according to the following equations (Clark, 2003; Rishton et al., 2006): logB B cl = 0.152 AlogP − 0.0148 PSA + 0.139 logB B ri = 0.155 AlogP − 0.01 PSA + 0.164 where PSA and AlogP are the polar surface area and atom-based octanol-water partitioning coefficient. The in-house PyMol script was applied to calculate the Buried Surface Area (BSA) of the peptide-ligand complexes according to the equation: BSA = ( AS A pep + AS A lig ) − AS A comp 2 where ASA _pep and ASA _lig and ASA _comp are the accessible surface areas of the peptide, ligand, and complex components. Molecular descriptors, such as molecular weight (MW), AlogP, H-bond acceptors (HBA), H-bond donors (HBD), PSA, quantitative estimate of drug-likeness (QED), and Tanimoto molecular similarity (T) indexes were calculated with the Biscu-it™ tools and Rcpi and Rcdk libraries within the Python and R environments (Voicu et al., 2020; Cao et al., 2015). Buried surface areas (BSA) have been calculated for 2D Aβ-ligand interaction diagrams using the NACCESS program (Hubbard and Thornton, 1993; Wallace et al., 1995). Molecular graphics and visualization were performed with the LigPlot + program (EMBL-EBI, Wellcome Trust Genome Campus, Hinxton, United Kingdom) in order to build two-dimensional interaction diagrams from three-dimensional coordinates. Linear regression analysis, followed by graphical representation was performed by PyMol and GraphPad Prism v.8 (GraphPad Software, San Diego, CA, United States). The differences were considered statistically significant at a p-value of <0.05.

Before molecular docking and MD simulations, the DrugCentral database (n = 4,642) was filtered to find the molecules, containing the propanesulfonic scaffold of TRA. The scaffold search protocol is based on the chemical similarity searching that can be used to screen a database of compounds for structural similarity to a query chemical structure. The small portion (n = 13) of FDA and EMA-approved drugs (Table 1), containing propanesulfonic scaffold of TRA, was identified as several hit molecules suitable for further investigation using molecular docking to eventually locate some lead candidates. In fact, our study relies on the modified hit-to-lead protocol to identify promising lead compounds, inhibiting Aβ in order to highlight the possibility for their optimization.

The identified hit molecules belong to different chemical classes and can be prescribed to treat various pathological conditions, such as bacterial infections (leprosy and tuberculosis), lead or mercury poisoning, metabolic disorders, and neurological conditions. Subsequently, the rigid-flexible molecular docking was performed to calculate the affinity between hit substances and Aβ, targeting the HHQK subregion at the N-terminus of the peptide. The lowest-energy Aβ conformer of the NMR models (n = 10) was chosen as a receptor molecule for molecular docking with the potential energy of -1,106.70 kcal/mol, which makes it more stable in the solution. Additionally, the 3D alignment of the Aβ peptide and HHQK region reviled the 2.0 Å structural deviation indicating a relatively “conserved” protein-receptor binding site in solution NMR models (Supplementary Figures S1A,B).

To accurately assess any correlations between the number of atoms in the ligand molecule and the conformational effects happening within the peptide-ligand binding site, we estimated the root-mean-square (RMS) difference between the top two conformations and its average value calculated between all conformations and the lowest-energy conformation in the lowest-energy cluster (Supplementary Figure S2). A strong positive correlation (r ² = 0.78) with reliable statistics (p-value = 0.02) between the clRMS, as an RMS difference between top two conformations in the largest cluster, and the number of torsions (N _tor ) was observed by linear regression analysis (Figure 1A). Furthermore, a significant standard deviation of the clRMSa variable, as an average of the RMS difference between all conformations and the lowest energy conformation, was detected for some hit compounds (3, 4, 8, and 14), describing the elevated conformational change (Figure 1B), which greatly affects the correlation coefficient (r ² = 0.63).

Finally, menadione and camphotamide lead candidates were found among the top binders (Table 2 and Figures 2A,B) with the binding energy values (−7.11 and −6.82 kcal/mol) significantly lower than for the reference compound with a ΔG _bind value of −6.19 kcal/mol (Table 2 and Figure 2C). It was previously published that menadione sodium bisulfite as a remedy against hemorrhagic disease caused by vitamin K deficiency could inhibit Aβ toxic formation and aggregation in C. elegans, extending its life span and reducing disruption of cellular membranes (Zhang et al., 2018). The co-authors also hypothesized that menadione might inhibit amyloid formation due to its backbone similarity to 1,4-naphthoquinone, which shows strong anti-aggregation effects on amyloidogenic proteins, such as insulin and α-synuclein. Similarly, the 6-hydroxy-nicotine intermediate from P. nicotinovorans as a precursor for camphotamide has been shown to have neuroprotective effects with putative applications in AD treatment (Wang et al., 2015; Hritcu et al., 2017). As a positive inotropic agent (cardiotonic), this substance could be potentially prescribed to treat AD-associated amyloid cardiomyopathy, inhibiting, directly and indirectly, Aβ accumulation in the heart. Indeed, the inhibition of Aβ and the formation of neurofibrillary tangles along the heart-brain axis could be reached through targeted supplementation of neurotrophic factors to the brain as it was hypothesized by Shityakov and coauthors Shityakov et al. (2021).

To investigate how molecular similarity contributes to the Aβ binding, the Tanimoto coefficients as an appropriate choice for a fingerprint-based similarity and ligand efficiency as the affinity normalized by the number of non-hydrogen atoms were implemented. A moderate negative correlation (r ² ≈ 0.5) with reliable statistics (p-value = 0.008) between the T and LE parameters was observed by linear regression analysis (Figure 3).

In other words, if chemical compounds match more to the TRA scaffold, they also have very similar binding modes calculated per number of heavy atoms in the molecule. In addition, all hit molecules were subdivided into 4 big groups or clusters by the hierarchical clustering according to their E-state molecular fingerprints and pairwise similarity matrix, where the lead molecules (compounds 4 and 6) ended up in the same cluster (Figure 3B).

Next, the AlogP, QED, PSA, and logBB values for menadione, camphotamide, and TRA as a reference, and donepezil as standard control were calculated to evaluate the drug ability to possess the optimal BBB permeation properties and druglikness (Table 3).

The AlogP parameter for lead candidates suggested their mild lipophilic properties, which were significantly higher than for the reference substance but not as good as for the standard control. On the other hand, the QED parameters for lead candidates were very similar to those determined for donepezil, reflecting the distribution of molecular properties, such as the number of hydrogen bond donors and acceptors, the number of aromatic rings, and the presence of unwanted chemical functionalities. However, for a CNS-active compound to permeate the BBB, a surface area less than 60–70 Å² is usually required to achieve the desired brain bioavailability for the administered drug (Shityakov et al., 2014; Kelder et al., 1999). Judging by relatively low AlogP and high PSA, there is a risk that the BBB permeation of the lead candidates might not be entirely sufficient as confirmed by the logBB, which optimally should be logBB > 0 (Shityakov et al., 2015). Therefore, further lead optimization by in situ enumeration or fragment replacement might be a plausible choice to make chemical modifications, improving brain bioavailability (logBB) together with selectivity, pharmacokinetic and pharmacodynamic parameters, and decreasing toxicity. Reaction-based in situ enumeration can be achieved by proposing synthetically feasible candidates from a reagent library in particular, which reagents are most likely to produce potential active compounds for the binding site (Mok et al., 2017). On the contrary, the fragment replacement protocol performs scaffold hopping by replacing part of the scaffold structure while maintaining the favorable binding between the receptor and the ligand (Vainio et al., 2013). The elevated experimental logBB index of 0.89 for donepezil was mainly due to its carrier-mediated transport to the brain (Kim et al., 2010). It could probably be accomplished by organic cation or choline transporters (OCT1-3 and CHT1) and not just crossing the BBB via passive diffusion through endothelial cells (Kim et al., 2010).

To validate further molecular docking results, the free energy of binding based on implicit solvation models was calculated for amyloid-drug complexes. The MM-PBSA/GBSA calculations (Table 4), using 100 ns MD trajectories, confirmed the previous data completely (GBSA) and partially (PBSA), revealing much higher binding affinities of lead compounds to Aβ in comparison to the reference.

The former protocol provided the best affinity for menadione (ΔG _GB = −19.71 kcal/mol), probably because in some extensive studies the MM-GBSA approach was computationally more efficient, achieving better accuracy but being less rigorous (Hou et al., 2011). Therefore, the second validation was done to clarify the MM-PBSA/GBSA discrepancy by utilizing the Prime algorithm with the MM-GBSA-based scoring function. As a result, the G-score data and the entropy-enthalpy compensation analysis confirmed the previous findings, describing the exothermic nature (ΔH < 0) of the binding process with the decreased disorder (TΔS <0), which could occur spontaneously depending on temperature. Additionally, the BSA values (Table 4 and Figures 4A–C), which measures the size of the Aβ-drug interface, also confirmed the previous binding affinity pattern, where the BSA elevation leads to an increase in binding as it was already published for the peptide-ligand and host-guest systems (Shityakov et al., 2020; Esmaeilpour et al., 2021).

To analyse the movements of the studied complexes, the root-mean-square deviation (RMSD) and fluctuation (RMSF) values together with the radius of gyration (R _g ) with respect to the initial conformation were plotted versus time (Supplementary Figures S3A–E). The receptor RMSD values were relatively high (RMSD = 10 Å) stabilized after about 20 ns for the compound 4 and 6/Aβ complexes and after 40 ns for TRA/Aβ (Supplementary Figure S3A). The ligand RMSD remained low within 1.0 Å and stabilized almost instantly (Supplementary Figure S3B). The receptor RMSF values produced some picks associated with the high flexibility of the Aβ termini (Supplementary Figure S3C). The atomic fluctuations of ligands showed more stable profiles especially for compounds 4 and 6 (Supplementary Figure S3D). The decreased R _g values were associated with the receptor compactness, which was elevated during the simulation (Supplementary Figure S3E). Finally, the number of H-bonds between the receptor and its ligands and the fraction of residues involved in H-bonding were assessed to find the H-binding contribution to the affinity. These parameters were the highest for compound 4 with more H-bonds formed, (Supplementary Figure S3F) significantly increasing the residue fraction (Supplementary Figure S3G).

Finally, per-residue and pairwise energy decomposition analyses were employed to evaluate the energetic contribution of drugs and the binding site residues. Some unique residues, such as Tyr10 (B) for TRA, Lys16 (A, B) for menadione, and Lys16 (B, C) for camphotamide were identified at the energetic threshold (ΔG = −3.0 kcal/mol) and below by one or both implicit solvation protocols (Figures 5, 6). In particular, all amino acid residues involved in the interaction with the lead candidates had exceeded the energetic threshold (Figure 5A,B), which was not observed for the reference molecule (Figure 5C). Meanwhile, the Lys16 targeting in Aβ by the oxidation of a catechol structure to the o-quinone forming the o-quinone-Aβ adduct is believed to be responsible for its anti-aggregation activity (Liu et al., 2017).

In the pairwise interaction, only Tyr10 (B) was detected slightly above the adjusted energetic threshold (ΔG = −6.0 kcal/mol), and Lys16 (B) had the lowest energy during the MD simulation of the Aβ-menadione complex (Figure 6). Highlighting the contribution of these residues for the inhibitor of this site, curcumin and resveratrol were found to interact with Arg5, Ser8, Tyr10, Gln15, Lys16, Leu17, and Phe20 (Fu et al., 2014). Besides, solvent-accessible residues, such as Phe20 at the N-terminal.

KLVFF stretch was reported as crucial in the interactions of these two inhibitors (Fu et al., 2014). Some studies have also documented the pairwise interactions for some amyloid variants (Aβ_10-35) to illustrate residue energetic contributions in a process of Aβ reorganization driven basically by inter-chain hydrophobic and hydrophilic interactions and also solvation/desolvation effects (Campanera and Pouplana 2010).

In this study, the scaffold searching approach based on known Aβ inhibitor tramiprosate to screen the DrugCentral database (n = 4,642) was employed to identify hit compounds (n = 13) and lead candidates (n = 2) for AD drug repurposing. Two lead candidates, namely menadione bisulfite and camphotamide, out of 13 hit compounds were identified as promising Aβ inhibitors with the improved ΔG _bind and logBB parameters. The binding affinity modes were also confirmed by molecular dynamics simulations using implicit solvation models, in particular MM-GBSA. We assume that this scaffold searching methodology in conjunction with pharmaceutical profiles (logBB) can be applied as a starting routine for drug repurposing in AD. Overall, the proposed computational pipeline can be implemented through the early stage rational drug design to nominate drugs for AD that, after additional in vitro and in vivo validation, could be readily evaluated in a clinical trial.