An Integrated QSAR-MD-DCCM Pipeline: A Predictive Computational Platform for the Rational Design and Dynamic Functional Validation of Dual-Target Directed Ligands

1. Introduction

The pharmaceutical landscape is continually shifting from the “one disease, one target, one drug” paradigm toward sophisticated therapeutic strategies capable of addressing the complex, multifactorial nature of prevalent age-related disorders. This transition has been mandated by the clinical limitations of single-target agents, which often succumb to pathway redundancy, compensatory biological mechanisms, and the high failure rate associated with modulating a single node within a dense biological network. Consequently, the concept of Multi-Target-Directed Ligands (MTDLs) has emerged as a cornerstone of modern drug discovery. MTDLs, also known as polypharmacology agents, are single chemical entities designed to simultaneously and beneficially modulate two or more validated targets implicated in a disease cascade. This approach is particularly advantageous for neurodegenerative and oncological diseases, where complex, chronic pathology necessitates intervention at multiple points to achieve clinical efficacy, improve patient adherence, and potentially mitigate drug resistance [1,2,3,4,5].

The strategic importance of MTDLs is perhaps most evident in the context of global health challenges posed by cancer and Alzheimer’s disease (AD). Together, these two pathologies represent a massive public health burden, characterized by high morbidity, soaring socioeconomic costs, and a global patient population that rises exponentially with age [6,7]. While seemingly disparate, accumulating epidemiological evidence points to an intriguing and often paradoxical inverse relationship between cancer and AD incidence, suggesting shared underlying molecular foundations [8,9,10,11,12,13]. At the cellular level, both diseases converge on dysregulated processes, including chronic inflammation, oxidative stress, mitochondrial dysfunction, and, critically, abnormalities in microtubule (MT) dynamics [14,15,16,17]. This shared pathology establishes a robust rationale for designing dual-acting agents, which can exploit these common mechanisms to provide synergistic therapeutic benefits and address the complications arising from co-morbidity in elderly patients.

One of the most compelling common denominators is the microtubule system. Microtubules, dynamic polymers of α- and β-tubulin heterodimers, are essential cellular components performing distinct but vital roles in both cancer and neurodegeneration. In cancer cells, microtubules form the mitotic spindle, and inhibitors targeting tubulin polymerization are staples of chemotherapy (vinca alkaloids and taxanes). These agents destabilize or hyper-stabilize MTs, leading to mitotic arrest and apoptosis in rapidly dividing cells [18,19]. In Alzheimer’s disease, however, MTs are central to neuronal architecture, axonal transport, and synaptic integrity. Hyperphosphorylation of the microtubule-associated protein tau causes MT destabilization, leading to the collapse of the axonal network and the formation of neurofibrillary tangles, which is a hallmark of AD pathophysiology [20,21,22]. Therefore, a single MTDL capable of binding to the tubulin interface, specifically the colchicine binding site on the β-tubulin subunit offers a promising mechanism: preventing uncontrolled proliferation in cancer and protecting against cytoskeletal collapse and neural degeneration in AD [23,24]

Complementing this structural target, the acetylcholinesterase (AChE) enzyme remains the primary symptomatic target for AD. AChE inhibitors (Donepezil, Rivastigmine, etc.) function by blocking the hydrolysis of the neurotransmitter acetylcholine (ACh) in the synaptic cleft, thereby enhancing cholinergic transmission and temporarily improving cognitive function [25,26] However, these single-target agents fail to address the underlying disease progression or the MT pathology. The inclusion of an AChE inhibitory domain in a tubulin-targeting scaffold creates a synergistic MTDL: the tubulin-modulating component targets the neurodegenerative etiology, while the AChE-inhibitory component provides necessary symptomatic relief. Crucially, the structural requirements for compounds targeting both the β-tubulin colchicine site and the AChE catalytic/peripheral anionic sites often involve common aromatic and trimethoxyphenyl scaffolds [27,28]. This chemical compatibility reduces the synthetic complexity of creating a potent MTDL, making the dual-target approach both biologically and chemically sound.

The complexity inherent in designing a compound to fit two distinct binding sites and navigate the necessary pharmacokinetic barriers (such as the Blood-Brain Barrier (BBB), crucial for central nervous system targets like AChE) makes traditional high-throughput screening prohibitively costly and time-consuming. This challenge necessitates the application of advanced Computer-Aided Drug Design (CADD) methodologies, an approach central to the field of Computational Chemistry and Chemical Informatics [29,30]. Quantitative Structure-Activity Relationship (QSAR) modeling stands out as an indispensable tool, enabling the rapid and cost-effective prediction of biological activity based on key molecular descriptors. QSAR transforms empirical activity data into mathematical models that can rationally guide chemical modifications and filter large virtual libraries before synthesis [31,32]. Following QSAR-guided optimization, molecular docking provides an essential structural validation, predicting the preferred binding pose and affinity (ΔG_binding) of the MTDL for each target [33,34].

However, static docking predictions often fail to account for the essential dynamic nature of protein-ligand interactions within a physiological environment. Proteins, particularly flexible targets like tubulin and multi-domain enzymes like AChE, undergo significant conformational changes upon ligand binding. This methodological gap is bridged by Molecular Dynamics (MD) simulation, which provides a high-resolution, time-dependent view of complex stability, pose persistence, and the nature of continuous residue contacts (hydrophobic, hydrogen bonding) over hundreds of nanoseconds [35,36]. Furthermore, advanced post-MD analyses, such as the Dynamic Cross-Correlation Map (DCCM), move beyond simple stability checks (Root Mean Square Deviation, RMSD) to uncover essential allosteric or dynamic coupling effects between distant protein domains [37]. Applying this integrated QSAR-MD-DCCM pipeline ensures that predicted MTDL candidates are not only computationally potent but also dynamically stable and structurally valid under simulated physiological conditions, a level of rigor demanded by contemporary chemical informatics research.

Driven by the need for structurally validated MTDLs against the converging pathologies of cancer and AD, and leveraging the immense power of integrated computational methodologies, the present study introduces a robust QSAR-MD-DCCM workflow for the rational design and virtual validation of novel dual-acting inhibitors. Starting from a focused library of trimethoxyphenyl-containing analogues, we established two predictive QSAR models to guide the rational design of sixteen novel compounds. These compounds were subsequently subjected to extensive in silico ADMET and toxicity profiling to prioritize candidates with optimal drug-like properties and BBB permeability. The most promising candidates were then evaluated against the colchicine binding site of β-tubulin (PDB ID: 4O2B) and the active site of AChE (PDB ID: 1EVE) via molecular docking and rigorous 100 ns MD simulations, and binding free energy calculations (MM-GBSA/MM-PBSA). We conclude by presenting a detailed analysis of ligand stability, binding differentiation, and complex dynamics, using DCCM to provide unprecedented mechanistic insights into the functional allosteric consequences of MTDL binding (Figure 1).

2. Results and Discussion

2.1. Predictive Validation: QSAR Model Performance

2.1.1. Statistical Validation and Predictive Power

The systematic application of the GA coupled with Multiple Linear Regression (MLR) yielded two statistically robust QSAR models, Model 1 and Model 2, aimed at predicting the inhibitory activity (pIC₅₀) of trimethoxyphenyl analogs against tubulin polymerization. The performance metrics (Table 1) confirmed that both models adhered strictly to the stringent OECD validation criteria, demonstrating exceptional internal consistency and external predictive capability.

The high R² values (Model 1: 0.8319; Model 2: 0.8384) demonstrated that over 83% of the variance in the experimental pIC₅₀ was explained by the selected descriptors. Critically, the strong cross-validation metrics, including Q² (LOO) (>0.79) and $R_{ext}^2$ (>0.85), confirmed the excellent internal predictability and external generalization of both models, respectively, ensuring they are not overfitted to the training data. Furthermore, the Y-randomization test results (Model 1: $R_{Yscr}^2$= 0.0866; Model 2: $R_{Yscr}^2$ = 0.0861) demonstrated a clear absence of chance correlation, firmly validating the mathematical significance of the derived equations and underpinning the reliability of the predictive models.

Model 1 Equation (Topological/Electrostatic Fit):

$${\mathrm{p}\mathrm{I}\mathrm{C}}_{50}=-11.8572+0.9217·{\mathrm{R}3\mathrm{s}}^{++}+0.9359·\mathrm{m}\mathrm{a}\mathrm{p}4\_26+34.0115·\mathrm{o}\mathrm{f}0\mathrm{u}\mathrm{g}+6.4191·\mathrm{E}\mathrm{T}\mathrm{A}\_\mathrm{s}\mathrm{h}\mathrm{a}\mathrm{p}\mathrm{e}\_\mathrm{p}$$

Model 2 Equation (Polarizability/Flexibility):

$${\mathrm{p}\mathrm{I}\mathrm{C}}_{50}=-12.1211+0.2414·\mathrm{M}\mathrm{o}\mathrm{r}15\mathrm{i}+0.9746·\mathrm{m}\mathrm{a}\mathrm{p}4\_26+33.5921·\mathrm{o}\mathrm{f}0\mathrm{u}\mathrm{g}+6.8459·\mathrm{E}\mathrm{T}\mathrm{A}\_\mathrm{s}\mathrm{h}\mathrm{a}\mathrm{p}\mathrm{e}\_\mathrm{p}$$

2.1.2. Mechanistic Interpretation of Key Descriptors

The efficacy of the QSAR approach, a pillar of this integrated methodology, is strongly validated by the chemical relevance of the selected descriptors, which collectively characterize the optimal molecular features for inhibiting tubulin polymerization at the colchicine binding site (Table S3, SI).

of0ug (Molecular Polarizability) emerged as the most dominant determinant, featuring the highest coefficient in both models (Model 1: 34.0115; Model 2: 33.5921). This descriptor quantifies a molecule’s ability to redistribute its electron density, indicating that high electronic polarizability is essential for strong activity. This property facilitates strong van der Waals and induced dipole interactions, which are characteristic of non-covalent binding within the highly aromatic environment of the tubulin colchicine site.

map4_26 (Topological Complexity), present in both models, represents the topological complexity and molecular architecture. Its positive contribution suggests that a precise degree of molecular complexity and specific branching pattern is required for optimal recognition and fit within the constricted colchicine binding pocket.

ETA_shape_p (3D Conformational Flexibility), also highly weighted in both models (coefficient around 6.5), quantifies the responsiveness of the molecule’s three-dimensional shape to conformational changes. The strong positive correlation confirms that a certain level of molecular flexibility and shape adaptability is advantageous, enabling the ligand to optimally adjust its conformation to the dynamic structural changes that occur in β-tubulin during polymerization inhibition.

The two models are further differentiated by unique descriptors: R3s++ (Steric Accessibility) in Model 1 emphasizes that reduced steric hindrance is a prerequisite for effective binding within the spatially confined pocket. Conversely, Mor15i (Electronic Density) in Model 2, a 3D-MoRSE descriptor, reinforces the need for electronic complementarity and specific mass distribution to fine-tune the inhibitory potential.

Given its slightly superior R² and the inclusion of descriptors emphasizing both electronic density (Mor15i) and flexibility, Model 2 was prioritized as the primary design tool. Model 1 served as an essential cross-validation tool focused on steric fit, ensuring that the subsequent designed analogs were simultaneously robust across key electronic and spatial parameters, a dual requirement that critically informed the rational structural modifications detailed next.

2.2. Design and Preclinical Filtering: Novel Analogs and ADMET Analysis

2.2.1. QSAR-Guided Rational Design

Leveraging the mechanistic insights generated by the predictive QSAR models, particularly the high importance placed on molecular polarizability (of0ug), topological complexity (map4_26), and conformational flexibility (ETA_shape_p), a panel of 16 novel compounds was rationally designed around the trimethoxyphenyl core. The primary objective was to systematically perturb these key descriptors to optimize the predicted pIC₅₀ and enhance the compound’s dual-targeting profile. Structural modifications included the addition of polar groups (amino, benzoic acid), substitution of heteroatoms (chlorine, sulfur), and the incorporation of fused aromatic rings (naphthalene, pyridone-fused systems). These changes were targeted to refine lipophilicity, rigidity, and the capacity for π-stacking interactions essential for both tubulin and AChE binding sites. The specific design rationale for each compound (Table 2), systematically links the chemical modification to the targeted descriptor perturbation. All newly conceived structures were optimized using the MMFF94 force field in MarvinSketch to obtain energetically favorable 3D coordinates, which served as the consistent input for all subsequent in silico validation steps.

Figure 2. Designed compounds using generated QSAR model.

Figure 2. Designed compounds using generated QSAR model.

2.2.2. ADMET Profile and Lead Candidate Selection

Following the rational design, all 16 novel compounds underwent comprehensive ADMET and toxicity screening using ADMETlab 3.0 and ProTox-II platforms to ensure their preclinical viability (Table 3). This pre-screening filter was critical for identifying candidates with favorable pharmacokinetic profiles before proceeding to resource-intensive structural simulations. Given the dual requirement for systemic efficacy (cancer) and CNS penetration (AChE targeting), selection criteria prioritized candidates with high BBB permeability, favorable lipophilicity (logP ideally 2-5) for good oral bioavailability, and minimal predicted safety liabilities.

This systematic, multi-parameter evaluation led to the identification of Compounds 15 and 16 as the optimal lead candidates. Compound 15 displayed favorable lipophilicity (log P of 2.417, log D of 2.718) and a predicted BBB value of 0.564, indicating high bioavailability and the potential to cross the BBB for efficient AChE targeting. Its ProTox analysis suggested a low probability of key toxicities. Compound 16 showed optimized lipophilicity (log P of 3.229, log D of 3.073) and an exceptionally high BBB penetration value (0.998), suggesting highly effective CNS targeting. While it exhibited a safe profile for neurotoxicity and cardiotoxicity, its elevated hepatotoxicity score and CYP2E1 interaction were noted as areas requiring future experimental investigation. The strong, yet complementary, dual-target pharmacokinetic profiles of Compounds 15 and 16 justified their advancement for the next rigorous stage: molecular docking and dynamic simulations to validate their structural binding hypotheses.

2.3. Structural Validation: Dual-Target Molecular Docking and Affinity Assessment

The initial structural validation of the rigorous computational workflow was performed by assessing the static binding affinity (ΔG) and detailed interaction profiles of the QSAR-filtered lead candidates, compounds 15 and 16, against their respective targets. For the anticancer target, the colchicine-binding site of Tubulin (PDB ID: 4O2B) was used, and for the anti-Alzheimer’s target, the active site of AChE (PDB ID: 1EVE) was employed. The accuracy of the docking protocol was rigorously validated by successfully reproducing the crystallographic poses of the reference ligands, Colchicine (RMSD: 0.172 Å) and Donepezil (RMSD: 0.532 Å), confirming the geometric accuracy of the method (Figure 3).

The calculated binding free energies (Table 4) revealed a critical target-differentiated potency, strongly supporting the feasibility of the dual-targeting strategy. Compound 16 showed a marginally superior predicted affinity for the Tubulin binding site (ΔG: -10.0 kcal/mol), making it the prime anticancer candidate, while compound 15 exhibited superior predicted affinity for the AChE active site (ΔG: -10.6 kcal/mol), confirming its strong potential as an anti-Alzheimer’s agent.

2.3.1. Tubulin Binding Profile (Anticancer Target)

The affinity of compound 16 (ΔG: -10.0 kcal/mol) was nearly identical to the positive control, Colchicine (ΔG: -10.1 kcal/mol), suggesting highly effective molecular recognition at the interface of the α- and β-tubulin subunits. Compound 16 achieved this potent binding via extensive interaction profile: Hydrogen Bonding with key residues Gln11, Asn101, and Ala180 (Chain A) and Asn258 (Chain B); and π-Alkyl Contacts with hydrophobic residues Leu248, Ala250, Lys254, and Leu255 (Chain B). This widespread engagement, spanning both chains, explains its superior affinity and strong propensity to destabilize the microtubule structure. Compound 15 (ΔG: -9.7 kcal/mol) displayed a similar but less extensive profile, relying primarily on π-Alkyl interactions (Chain B) and a single hydrogen bond to Asn258.

2.3.2. Acetylcholinesterase Binding Profile (Anti-Alzheimer’s Target).

Compound 15 demonstrated a strong affinity (ΔG of -10.6 kcal/mol), approaching the clinical standard Donepezil (ΔG: -11.8 kcal/mol) and confirming its potential for cholinesterase inhibition. Its interaction profile was diverse and indicative of a dual-binding site mechanism: Hydrophobic Interactions, including π-Alkyl contacts at Phe290 and Phe331; π-π Stacking interactions with aromatic residues, notably Tyr121, Trp279, Phe330, and Tyr334; and Hydrogen Bonding with Asp72 and Gly118. The trimethoxyphenyl scaffold effectively engaged the peripheral anionic site (PAS) residues, while the rest of the molecule penetrated deeper into the active site Gorge. Compound 16 (ΔG: -9.7 kcal/mol) also formed favorable contacts, but lacked the extensive π-π stacking seen in compound 15, resulting in a diminished affinity. The complementary target affinity profiles observed in the docking study were used to guide the subsequent rigorous dynamic validation step.

2.4. Dynamic Validation: 100 ns MD Simulation

2.4.1. MD Simulation and Stability Analysis

To move beyond static binding predictions and perform a rigorous validation of the integrated workflow, 100 ns all-atom MD simulations were conducted for the four final complexes, demonstrating the dynamic behavior of the lead compounds under simulated physiological conditions. This process served as a high-resolution filter, confirming the stability and permanence of the strong binding poses observed in the docking study.

The RMSD relative to the protein served as the primary metric for tracking conformational stability (Figure 4). The compound 16-Tubulin complex exhibited superior stability, maintaining RMSD values largely within the 1.0 - 4.0 Å range, with the main body of the compound remaining fixed inside the colchicine binding site for the entire simulation. The Compound 15-Tubulin complex also proved highly rigid, staying well below 2.0 Å RMSD for most of the trajectory. Conversely, the AChE complexes showed slightly higher, yet constrained, fluctuations: Compound 15-AChE fluctuated between 1.8 - 4.8 Å, settling into a stable secondary conformation after an initial ring rotation, while compound 16-AChE remained rigid, oscillating between 1.6 - 3.2 Å.

The RMSF of the ligands further differentiated the binding dynamics, indicating that the compounds utilize complementary stabilization mechanisms. Compound 15 displayed high flexibility in the AChE site (RMSF: 3.0 Å), suggesting adaptation to engage multiple hydrophobic residues and continuous H-bonding, but low flexibility in the Tubulin site (RMSF: 1.5 Å). Compound 16 showed the reverse trend: low flexibility in the AChE pocket (RMSF: 1.5 Å) but higher flexibility in the Tubulin pocket (RMSF: 3.0 Å), indicating the motion required to react to the fluctuating tubulin structure and interact with multiple residues.

2.5. Energetic Validation: MM-GBSA/MM-PBSA Analysis

The MM-GBSA/MM-PBSA decomposition of the MD trajectories provided valuable insights into the energetic basis of binding for compounds 15 and 16 with both tubulin and AChE, supporting their potential as dual-target anticancer and anti-Alzheimer’s agents. In all four complex systems, the gas-phase contributions were markedly favorable and were dominated by van der Waals and electrostatic factors. This highlighted the importance of tight packing and aromatic stacking with residues embedded within both the colchicine site of tubulin and the aromatic active site region of AChE. The consistently stabilizing van der Waals profiles reflected the contribution of the tri-methoxy substituted phenyl core enabling the extensive hydrophobic complementarity. The polar solvation penalties observed are typical for the positioning of ligand within deep hydrophobic binding sites. Both compounds display comparably favorable total binding energies across tubulin and AChE, supporting their suitability as dual-acting candidates. Compound 15, with an extended biphenyl substituent, achieved stable binding through enhanced π-stacking and hydrophobic enclosure. This is particularly advantageous for tubulin interactions but also remains beneficial in AChE. Compound 16, with a heterocyclic imide-linked fused ring system, showed a good balance between electrostatic interactions and solvation effects. This behavior works together with the hydrophobic contribution of its trimethoxyphenyl group, allowing it to interact favorably within both binding pockets. This difference in structural features between the two ligands helps explain why both maintained favorable binding energies across the two targets despite of their varying energetic contributions. Overall, the MM-GBSA/MM-PBSA results supported the selection of 15 and 16 as strong candidates for dual inhibition. Both compounds demonstrated the ability to bind effectively to tubulin as well as AChE. They used slightly different but complementary interaction patterns. Their binding energies remained stable over the full 100 ns MD simulations. Their structural flexibility allowed them to adapt well within two very different binding environments. Together, these findings highlighted their promise as multifunctional therapeutic leads with potential relevance in both cancer treatment and neurodegenerative disease management.

Figure 5. MM-GBSA/MM-PBSA analysis of binding free energy (ΔG_bind) for compounds 15 and 16. The analysis compares the energetic contributions (van der Waals, Electrostatic, Solvation) to the stable binding of the lead candidates against the Anti-Alzheimer (AChE) and Anticancer (Tubulin) targets. Favorable ΔG_bind values across both MM-GBSA (gb) and MM-PBSA (pb) methods support their dual-inhibitory potential.

Figure 5. MM-GBSA/MM-PBSA analysis of binding free energy (ΔG_bind) for compounds 15 and 16. The analysis compares the energetic contributions (van der Waals, Electrostatic, Solvation) to the stable binding of the lead candidates against the Anti-Alzheimer (AChE) and Anticancer (Tubulin) targets. Favorable ΔG_bind values across both MM-GBSA (gb) and MM-PBSA (pb) methods support their dual-inhibitory potential.

2.6. Functional Validation: Post-MD Mechanistic Insights

The final methodological steps involved performing functional analysis on the 100 ns trajectories, serving as the ultimate verification point of the predictive workflow.

2.6.1. The Ramachandran Plot Analysis

This analysis was performed on the final (100 ns) snapshot of each molecular dynamics trajectory using the MolProbity server to rigorously validate the stereochemical quality and conformational integrity of the target proteins post-ligand binding (Figure 6). The results consistently demonstrated that the ligand-induced conformations for all four complexes were structurally sound and highly reliable. Specifically, the complexes showed high percentages of residues in the favored regions (ranging from 86.5% to 91.2%) and allowed regions (ranging from 96.8% to 99.1%). A quality threshold of >85% of residues in the favored regions is considered excellent for post-MD analyses, and the low number of outliers, confined primarily to flexible loop regions, unequivocally confirmed the stereochemical validity of the models.

2.6.2. Dynamic Cross-Correlation Matrix (DCCM) Analysis

The final, high-value component of the integrated pipeline was the DCCM Analysis. This analysis quantifies the concerted motions between residue pairs, providing a functional, allosteric context for the observed stability. The DCCM maps (Figure 7) confirmed the compounds’ complementary dynamic signatures.

For Tubulin, compound 15 exhibited a complex pattern of both positive and negative correlations, suggesting it engages in broader, more widespread structural adjustments across the α/β interface (greater structural flexibility). Conversely, compound 16 showed a more simplified and localized correlation profile, aligning with its rigid, focused binding mode.

For AChE, both complexes maintained positive correlations in critical pairs (Tyr70-Trp84, Tyr121-Glu199, and Tyr334-His440). Crucially, compound 15 showed a distinct negative correlation between Gly441-Tyr70, a key pair regulating the gorge entrance. The presence of this anti-correlation in 15 suggests it is more effective at modulating the dynamic opening/closing of the AChE gorge, reinforcing its superior predicted inhibitory affinity and deeper mechanistic engagement observed in the docking study.

The DCCM analysis serves as the functional validation of the entire predictive workflow, reinforcing that compound 15 employs a broader, more flexible engagement, while compound 16 uses a more focused, rigid strategy, confirming the versatility and therapeutic promise of the leads.

3. Materials and Methods

3.1. Dataset Preparation and QSAR Modeling Rationale

The development of robust QSAR models commenced with an exhaustive literature search utilizing the ChEMBL and BindingDB databases to create a high-quality dataset of 57 molecules (Table S3, SI). The key inclusion criterion was the presence of the trimethoxyphenyl scaffold, a pharmacophore known for its affinity to the colchicine binding site on tubulin and documented tubulin polymerization inhibitory activity (IC₅₀). Biological activity values (IC₅₀) were uniformly converted to the negative logarithmic scale (pIC₅₀) to ensure a normal distribution and reduce data skewness, which is essential for stable Multiple Linear Regression (MLR) modeling. The final dataset of 57 compounds was systematically divided to ensure rigorous validation: 10 compounds (≈ 17.5%) were designated as the external validation set to assess predictive capability on unseen data, while the remaining 47 compounds formed the training set for model generation and internal cross-validation. This division strategy maximized the training set size, ensuring model stability and adequate descriptor space coverage while maintaining substantial external validation power, adhering strictly to OECD principles for QSAR development [31,38].

3.2. Descriptor Generation & Pre-Filtering

Molecular structures for all 57 compounds were rigorously optimized using the MMFF94 force field within MarvinSketch to ensure consistent, low-energy three-dimensional conformations. Subsequently, a massive initial pool of 11,829 molecular descriptors was generated using specialized software packages, including PaDEL-Descriptor, alvaDesc, MORDRED, and MERA via the OCHEM server [39]. This extensive descriptor space, encompassing constitutional, topological, geometrical, electronic, and hybrid properties, was intentionally broad to provide the Genetic Algorithm (GA) with the widest possible feature set for robust structure-activity relationship identification. To mitigate the inherent risks of data overfitting and multicollinearity associated with such large descriptor pools, a systematic two-stage pre-filtering procedure was implemented. First, descriptors exhibiting the same value across more than 80% of the dataset were removed; second, a threshold of 95% correlation was applied to eliminate redundant descriptor pairs. This methodical pre-filtering eliminated 7,943 descriptors, resulting in a refined and highly informative set of 3,886 descriptors for the subsequent QSAR model construction phase.

3.3. Molecular Sketching and Geometry Optimization

Leveraging the mechanistic insights derived from the QSAR models, a series of 16 novel compounds were rationally designed based on the original trimethoxyphenyl scaffold. The design strategy was guided by the positive contribution of the identified key descriptors. For instance, structural modifications such as the addition of chlorine or ethoxy groups were introduced to modulate Mor15i (electronic density) and MAP4_26 (topological complexity), while the introduction of fused rings or changes in substituent positions were aimed at optimizing ETA_shape_p (conformational flexibility) and steric fit (R3s++). Specific modifications included altering the aromaticity, introducing flexible linkers, and changing the polarity of key functional groups to maximize the predicted pIC₅₀ for tubulin inhibition. All newly sketched structures were subsequently subjected to geometry optimization using the MMFF94 force field in MarvinSketch to obtain energetically favorable 3D coordinates, which served as the input for all subsequent ADMET and docking calculations [40,41,42,43].

3.4. ADMET Prediction and Selection of Lead Candidates

Following the rational QSAR-guided design, all 16 novel compounds were subjected to comprehensive ADMET profiling to rigorously filter candidates before proceeding to resource-intensive structural simulations. This critical in silico screen employed the ADMETlab 3.0 and ProTox-II platforms to predict essential pharmacokinetic parameters, including Caco-2 permeability, BBB penetration, aqueous solubility (logS), lipophilicity (logP and logD), metabolism via CYP isoforms, and general toxicity risks such as hepatotoxicity and mutagenicity. Given the dual-target nature of the study, requiring systemic efficacy for cancer and CNS penetration for Alzheimer’s Disease, the selection criteria were stringent, prioritizing candidates with high BBB permeability (essential for AChE targeting), favorable logP (ideally 2-5) for oral bioavailability, and minimal predicted safety liabilities. This systematic, multi-parameter evaluation led to the identification of Compounds 15 and 16 as the optimal lead candidates. Compound 15 was selected for its balanced physicochemical profile and predicted BBB permeability, while Compound 16 was advanced due to its exceptional BBB penetration prediction, justifying their prioritization for the subsequent molecular docking and dynamics simulations to validate their structural binding hypotheses [44].

3.5. Molecular Docking & Validation

The dual-target potential of the two selected leads, Compounds 15 and 16, was structurally evaluated using molecular docking simulations against their respective targets. For the anticancer target, the crystal structure of β-tubulin in complex with colchicine (PDB ID: 4O2B) was used. For the anti-Alzheimer’s target, the crystal structure of Acetylcholinesterase (AChE) (PDB ID: 1EVE) was utilized. The proteins were pre-processed (removing water molecules, standardizing protonation states) using the Discovery Studio Visualizer and SwissPDB Viewer. Docking calculations were performed using AutoDock Vina integrated within the CB-Dock2 server, which automatically defines the optimal binding site cavity, ensuring an unbiased “blind docking” approach [45,46]. Validation of the docking methodology was performed by redocking the co-crystallized ligands (Colchicine for 4O2B and Donepezil for 1EVE). The method was accepted only after the Root Mean Square Deviation (RMSD) between the docked and crystal poses was less than 2.0 Å [47]. The final predicted binding affinities (ΔG in kcal/mol) and the nature of the specific protein-ligand interactions (H-bonds, π-stacking, hydrophobic contacts) were meticulously analyzed using Discovery Studio Visualizer and UCSF Chimera for both compounds against both targets [48].

3.6. Molecular Dynamics (MD) Simulation and Post-MD Analysis

To comprehensively assess the stability, flexibility, and permanence of the protein-ligand interactions observed in the static docking poses, 100 ns all-atom MD simulations were conducted for the four final complexes: 15-Tubulin, 16-Tubulin, 15-AChE, and 16-AChE. The simulations were performed using the Desmond module within the Schrödinger Suite. Each system was prepared by placing the docked complex in a periodic boundary condition box solvated with TIP3P water molecules and neutralized with Na⁺ and Cl^- ions. The OPLS3e force field was used to describe all atomic interactions. After a thorough minimization and pre-equilibration phase, the production run was executed for 100 ns under constant temperature (300 K) and pressure (1 bar) conditions, with trajectories saved every 50 ps [49].

Post-MD Analysis was critical for interpreting the dynamic stability and functional integrity: trajectory analysis (Protein/Ligand RMSD and residue RMSF) monitored stability and flexibility; interaction timelines tracked the persistence of specific residue contacts; the Ramachandran Plot validation ensured the stereochemical quality and conformational suitability of the final complex structures using the MolProbity server [50]; and finally, the advanced Dynamic Cross-Correlation Matrix (DCCM) analysis, performed with the MD-TASK suite, was used to quantify the functional connectivity and allosteric communication within the protein induced by ligand binding [51].

To provide the energetic basis for the observed dynamic stability, the binding free energy (ΔG_bind) for the four protein-ligand complexes (15-Tubulin, 16-Tubulin, 15-AChE, and 16-AChE) was calculated using the MM-GBSA and MM-PBSA methodologies. This was performed on the molecular dynamics trajectories using the Prime module of the Schrödinger Suite. The analysis was conducted on the final 50 ns (500 frames) of the 100 ns production run to ensure the inclusion of only equilibrated conformations. The free energy of binding was calculated according to the general equation:

$$∆G_{bind}=∆E_{MM}+∆G_{solv}-T∆S$$

where ΔEMM represents the change in molecular mechanics energy (sum of ΔE_vdw and ΔE_elec), ΔG_solv represents the change in solvation energy (polar and non-polar components), and TΔS is the change in conformational entropy. The non-polar contribution to solvation energy was calculated using the Solvent Accessible Surface Area (SASA) model. The entropy term (TΔS) was neglected due to its computational intensity, which is common practice for comparative studies. The focus was on the most favorable non-covalent terms: van der Waals and electrostatic energies, as they dominate the gas-phase energy in these systems.

4. Conclusions

This study successfully established an integrated QSAR-MD-DCCM pipeline for the rational discovery and dynamic validation of multi-target-directed ligands (MTDLs), satisfying all stringent OECD criteria. Two highly accurate, statistically robust QSAR models (R² > 0.83) were developed from a focused trimethoxyphenyl scaffold library, where mechanistic interpretation demonstrated that molecular polarizability (of0ug), topological complexity (map4_26), and conformational flexibility (ETA_shape_p) are the indispensable key determinants for effective tubulin polymerization inhibition. These quantitative insights directly guided the rational design of 16 novel analogues, demonstrating the pipeline’s utility in early-stage Machine Learning application. The leads identified through this initial filtering, compounds 15 and 16, were selected based on optimized drug-likeness and their high predicted capacity for BBB penetration, an essential feature for effective AChE targeting. Molecular docking and rigorous 100 ns MD simulations confirmed dynamically stable and persistent binding for both leads against both Tubulin (PDB ID: 4O2B) and AChE (PDB ID: 1EVE). Furthermore, MM-GBSA/MM-PBSA binding free energy calculations provided essential energetic validation, showing highly favorable total binding energies (ΔG_bind) that primarily stemmed from stabilizing van der Waals interactions. This analysis revealed a complementary, target-differentiated potency: Compound 16 showed marginally superior predicted affinity for the Tubulin binding site (ΔG: -10.0 kcal/mol), while compound 15 exhibited superior predicted affinity for the AChE active site (ΔG: -10.6 kcal/mol). Post-MD Ramachandran analysis confirmed the high stereochemical integrity of all four protein-ligand complexes. Furthermore, the DCCM analysis provided the final and highest level of functional validation, revealing complementary dynamic signatures: Compound 15 was found to induce a broader range of dynamic interactions and allosteric coupling, suggesting a more flexible engagement, while compound 16 exhibited a more focused and rigid inhibition strategy. In conclusion, this rigorous, integrated QSAR-MD-DCCM pipeline represents a highly predictive and robust platform for the de novo design and dynamic validation of multi-target therapeutics. The successful virtual validation of two distinct dual-acting candidates (15 and 16) serves as compelling proof-of-concept for this advanced computational workflow, ready for immediate experimental synthesis and in vitro validation.