Intrinsic Deformability of Human Acetylcholinesterase in Presence of Diverse Compounds: A Biophysical-Computational Study Based on Classical Dynamics, Elastic Network Models, Statistical Potentials, Energy Frustration, and Cavity Volumetry

1. Introduction

The enzyme known as acetylcholinesterase (AChE) is responsible for hydrolyzing acetylcholine and has an esterase activity. It plays a crucial role in the proper functioning of neurons through cholinergic pathways. AChE is primarily found in the synaptic spaces of both the central and peripheral nervous systems, as well as in the membranes of red blood cells. Its main function is to halt nerve impulses by catalyzing the hydrolysis of acetylcholine (ACh), which is a neurotransmitter involved in various cognitive processes within the cholinergic system. Consequently, the inhibition of the AChE enzyme using specific inhibitors has emerged as a promising therapeutic approach for managing a range of neurological disorders, including Lewy body dementia and Alzheimer’s disease [1,2]. Moreover, the accumulation of ACh in synapses can lead to pathological conditions. This occurs when ACh receptors, namely nicotinic and muscarinic receptors, are excessively stimulated, resulting in neurological ailments such as depression, dizziness, headaches, nausea, breathing difficulties, and even sudden death [3].

AChE inhibitors (IECA) are employed in the treatment of these disorders to enhance cholinergic activity by increasing the levels of ACh in cholinergic synapses [3]. Nonetheless, apprehensions regarding the toxicity and adverse effects associated with synthetic substances have spurred the quest for novel, safe, and effective natural compounds. In this regard, plant phenolic compounds have gained significant attention as extensively studied natural products due to their diverse biological applications, including the inhibition of AChE [4].

Indeed, numerous studies have indicated that huperzine A, a derivative of the Chinese herb Huperzia serrata, exhibits selectivity, potency, and reversibility as an inhibitor of the AChE receptor antagonist. Both preclinical research and clinical trials have demonstrated the potential therapeutic effect of huperzine A in Alzheimer’s disease (AD) treatment, although there is currently insufficient evidence to recommend its clinical use. Presently, huperzine A is utilized as a dietary supplement in certain countries [1]. Moreover, our research group has recently investigated other natural compounds, such as cyanidin and resveratrol, for their diverse therapeutic properties, particularly their potential as anti-Alzheimer’s alternatives, given their ability to interact with and inhibit AChE [5,6,7].

Cyanidin, known scientifically as 3,5,7,3’,4’-Pentahydroxyflavilium, 2-(3,4-Dihydroxyphenyl)-3,5,7-trihydroxy-1-Benzopyrylium, is a naturally occurring organic compound that falls under the categories of both flavonoids and anthocyanins. It serves as a pigment and is commonly found in a variety of fruits including blackberries, raspberries, grapes, cherries, blueberries, and even purple corn. Additionally, cyanidin can be found in fruits like apples and plums. The highest concentrations of this compound are typically found in the fruit peel [8].

Resveratrol, chemically known as 3,5,4’-trihydroxystilbene, belongs to the class of stilbenoids, which are natural phenolic compounds. It is classified as a phytoalexin, naturally synthesized by different plants, particularly in the skin of grapes, blueberries, raspberries, and blackberries [9].

The significance of assessing the inherent flexibility or deformability of protein complexes has been well-established. Understanding this intrinsic deformability, which is determined by their three-dimensional structures or the inter-residue contact topology, is crucial for comprehending the flexibility and dynamics of cooperative movements that are closely associated with the structural alterations observed during their biological functions [10,11,12,13,14].

Intrinsic conformational movements play a crucial role in protein function and are typically energetically advantageous. These movements are conserved throughout evolution among orthologs. As a result, it has been observed that the sequences and structures of proteins, shaped by evolutionary processes, rely on their intrinsic dynamics to carry out their specific functions. These intrinsic dynamics can be influenced by various physicochemical perturbations, including ligand-induced changes. Such changes have the potential to alter the preferred modes or conformations that the protein adopts during its functional processes. Particularly, the binding of ligands can regulate the intrinsic dynamics by modifying the energy landscape of the protein, often by influencing its curvature [14].

The intrinsic deformability of AChE has been extensively studied, revealing that the magnitude of its structural fluctuations can significantly increase when specific inhibitors bind to the primary entry channel leading to the active site. This binding effectively abolishes the enzyme’s activity towards natural substrates like acetylcholine. Exploring the intrinsic flexibility of AChE has also been studied by other authors in relation to synthetic compounds obtained through “click” chemistry, resulting in the development of exceptionally potent binding compounds [15,16,17]. This approach enables the design of inhibitors that elicit a pronounced “induced fit” response in AChE, leading to the exposure of previously solvent-hidden regions. Consequently, ligand-induced conformational rearrangements can facilitate the formation of enzyme complexes that optimize their inhibitory activity on AChE by modulating the transient movements associated with both the active and peripheral sites. These interactions can even occur in time-dependent conformations [15,16].

Efforts to investigate the intrinsic deformability of AChE have primarily utilized classical dynamics tools [15,16], and elastic network models [17,18]. These approaches have been employed to analyze the inherent conformational movements of this protein system in the presence and absence of ligands [11,19,20,21]. This type of comparative analysis is particularly valuable due to two main reasons. Firstly, the conformational flexibility of target proteins poses a significant challenge in accurately modeling and describing protein-inhibitor interactions [22,23,24]. Secondly, in theoretical inhibition assays, it is essential to determine whether the observed conformational changes, beyond those naturally controlled by the protein for its functional interactions, could be induced by the inhibitor [25].

In order to tackle the challenges and considerations mentioned above, it is important to highlight that a crucial element in the biophysical exploration of intrinsic conformational movements lies in the utilization of complementary analyses based on statistical energy functions or potentials. These functions describe local interactions along the protein chain, hydrophobic forces, and tertiary interactions. This approach has demonstrated its effectiveness in predicting and interpreting critical regions within proteins with respect to conformational changes, protein folding, and perturbations at complex interfaces [23,26,27].

Likewise, the utilization of algorithms inspired by energy landscape theory provides a means to quantify the localization and distribution of energy within protein structures. This concept proves to be valuable in gaining a deeper comprehension of the biological behavior of proteins and the impact of perturbations on their intrinsic energy. This intrinsic energy, which is conserved, is closely linked to conformational changes and plays a significant role in understanding protein dynamics [28,29,30].

In this study, our objective was to conduct a comprehensive biophysical-computational comparative analysis to investigate the potential impact of natural compounds of interest on the intrinsic deformability of AChE human. To achieve this, we employed models based on classical dynamics, elastic networks, statistical potentials, energy frustration, and volumetric analysis. This computational approach provides an integrative perspective by utilizing various tools to assess the extent of conformational changes induced by the ligands of interest on AChE. This protein system holds significant biomedical relevance, particularly in the context of neurodegenerative diseases.

2. Results and Discussion

2.1. Docking, Inhibitory Potency and Molecular Dynamics (MD) Analysis

Table 1 shows the results of molecular docking obtained by machine learning using the DockThor algorithm for the compounds considered in this study against AChE (Figure 1), and also shows the values of the post-classification performed to validate the viability of the links of the complexes formed using MM/PBSA before and after MD analysis at 200 ns.

All docking results were thermodynamically favorable regardless of the conditions considered in this study, with an overall mean binding energy of ΔG ≈ -12.2 kcal.mol^–1. According to DockThor calculations, the thermodynamic mean was ΔG ≈ -9.7 kcal.mol^–1, with a minimum of ΔG ≈ -8.9 kcal.mol^–1 and a maximum of ΔG ≈ -12.0 kcal.mol^–1. TZ2PA6 was the compound with the most favorable binding free energy (ΔG = -12.03 kcal.mol^–1) followed by Cyanidin (ΔG = -9.50 kcal.mol^–1), while the rest of the compounds presented similar energetics (ΔG ≈ -8.9 kcal.mol^–1) as predicted by DockThor. It is important to highlight that both the TZ2PA6 control and Cyanidin bound to residues adjacent to the active site reported for AChE, even highlighting the interaction with the aromatic residue Tyr332, which is also present in the active site of the highly related enzyme BuChE [31,32]. These results are significant because most of the binding affinities predicted by DockThor have been described to be moderate in the range of high to low micromolar affinity units (scores < -8 kcal.mol^–1 are associated with submicromolar affinities) [33].

On the other hand, MM/PBSA calculations showed an increase in the binding energies of all compounds after MD with a global mean that increased from ΔG_MM/PBSA ≈ -5.8 kcal.mol^–1 to ΔG_MM/PBSA ≈ -21.1 kcal.mol^–1. The thermodynamic mean of MM/PBSA was ΔG_MM/PBSA ≈ -13.5 kcal.mol^–1, with a minimum of ΔG_MM/PBSA ≈ -1.4 kcal.mol^–1 and a maximum of ΔG_MM/PBSA ≈ -38.8 kcal.mol^–1. As predicted by DockThor, TZ2PA6 was the compound with the most favorable relative binding energy calculated by MM/PBSA (ΔG_MM/PBSA =-38.78 kcal.mol^–1), followed by Huperzine A (ΔG_MM/PBSA =-22.74 kcal.mol^–1) and Cyanidin (ΔG_MM/PBSA =-17.21 kcal.mol–1) after 200 ns MD (Table 1; Figure 2). Resveratrol was one of the compounds with the least favorable relative binding energy in both DockThor predictions and after 200 ns MD according to MM/PBSA calculations (Table 1). These results are important because MM/PBSA calculations have been described to have an overall success rate of ≈74%, even improving the predictions made by efficient scoring functions such as Glide SP (≈59% success rate) [34].

The complexes were used in MD simulations to evaluate the stability of the compounds bound to AChE and select low-energy conformations to study the deformability and conformational and energy-structural changes of AChE in the presence of each ligand [24,35]. From the minimum energy structures generated by MD, empirical data were obtained for 1) MM/PBSA (Table 1), 2) RMSD (Figure 3a), 3) RMSF (Figure 3b), 4) Rg (Figure 3c), and 5) interaction depletion (Figure 3d) as suggested [24,36,37]. The means of the dynamic parameters calculated after 200 ns MD showed stable binding of all the compounds considered in this study to AChE regardless of the nature of the compounds (RMSD ≤ 0.55 Å, RMSF ≤ 0.32 Å and Rg ≤ 0.02 Å) (Figure 3). Cyanidin was the compound that presented the slightly more stable binding oriented towards structural folding (RMSD = 2.33 Å and RMSF = 1.50 Å) compared to the free protein (RMSD = 2.88 Å and RMSF = 1.83 Å) and the rest of the compounds (Figure 3) in terms of both RMSD and RMSF. In all cases, an increase in the type and number of interactions was predicted after 200 ns of MD, with the largest number of predicted interactions being of the hydrophobic type (Figure 3).

By calculating theoretical kinetic parameters of the compounds considered for which inhibition constants (K_i) or IC₅₀ values against AChE have not been described, such as Cyanidin and Resveratrol, a higher inhibitory potency was predicted for Cyanidin (pIC₅₀ = 12.2) compared to Resveratrol (pIC₅₀ = 7.4) (Table 2). When comparing the experimental kinetic calculations described for the rest of the compounds, it was observed that the theoretical kinetic calculations for Resveratrol are similar to those reported experimentally for Huperzine A after an average of 18 assays [38], but less favorable than those reported for TZ2PA6 [17]. The compound Bis-DMM was the compound with the least favorable experimentally reported kinetic parameters [39] (Table 2). These inhibition potency results are in line with most of the binding affinities predicted from the docking algorithm used before MD and based on machine learning, in which the compounds that presented the most favorable energies such as Cyanidin and TZ2PA6 tend to also present submicromolar affinities [33].

2.2. Analysis of Intrinsic Deformability Associated with AChE Conformational Changes

In relation to the structural deformability of the complexes formed between the compounds and AChE after 200 ns of MD, empirical data were obtained using coarse-grained computational strategies recommended to determine possible conformational changes in these dynamic systems (Table 3) [23]. In this sense, and in terms of statistical potentials, all the compounds induced local changes in the protein residues in terms of distance and solvent accessibility. Cyanidin was the compound that mediated the greatest perturbation of the intrinsic accessibility of the residues to the solvent compared to the rest of the complexes, and it was also able to alter the psi (Ψ) angles of the AChE system (Table 3).

On the other hand, it was additionally predicted that the greatest change induced by TZ2PA6 was in terms of the torsion of both phi (φ) and psi (Ψ) angles, while Huperzine A and Resveratrol affected the psi (Ψ) angles to a greater extent (Table 3). These results show how the compounds can perturb the stability of protein biomolecules in relation to the structural fluctuations related to their energy landscape in both free and bound states. This approach allows us to analyze the contribution of residues to the global folding of a conformation, as well as the relationship between the free energy of that conformation and the probability of observing that conformation, as has been reported [23,26].

The predictions of the intrinsic quasi-rigid regions or domains (Q) showed that Cyanidin, in addition to perturbing both specific angles and ASA, was the only compound that increased the number of Q regions, while the rest of the compounds were able to alter these regions but promoted a decrease in them. The increase in the predicted intrinsic Q domains in AChE in the presence of Cyanidin could be related to a redistribution of the intrinsic quasi-rigid regions of AChE, indicating a probable gain in rigidity of the protein system, as opposed to the probable flexibilization of the protein by the rest of the compounds. All complexes presented valid rigidity quality (QS) values, being similar or higher than that calculated for the free protein (Table 3). These perturbations in the Q regions are determinants because, thanks to the specific properties of each complex and the alterations of the free energy landscapes, they can reliably reproduce the structural fluctuations induced by the compounds. The Q region prediction method is based on the notion that, in genuinely rigid regions (validated by the internal quality parameter QS), the distances between two points remain strictly unchanged during MD [23,40].

In addition to the predicted stability at the MD level, statistical potentials, and Q regions, it was observed that all the compounds altered the number of intrinsic nodes (#Nodes) and links (#Links) associated with the alpha carbons and their interactions in AChE, with Cyanidin being the compound that mediated the greatest increase in both the number of nodes and their links, compared to the rest of the complexes and free AChE (Table 3). This greater number of predicted amino acid interaction networks in the AChE + Cyanidin complex have been shown to be valuable for both protein folding and establishing residue contributions to complex binding. In this sense, the stability of Cyanidin binding is due to its having a greater number of interactions of the residues that act as nodes (intersecting, connecting, or binding points) [22,41].

To further investigate and improve the comprehensiveness of the search for possible conformational changes in the complexes, the DynOmics-ENM tool was used, which integrates the two most widely used elastic network models (GNM and ANM). After applying these elastic network models to the complexes considered, it was observed in the same way that all the compounds altered the intrinsic conformational displacement of AChE at the level of the “collectivity” data of the modes, which is indicative of the proportion of conformations that move collectively together in the same mode. Cyanidin was the compound that caused the greatest collective displacement compared to the rest of the complexes and free AChE (Table 3), promoting increased flexibility at the ends of the chain (Figure 4).

The data suggest that the increase in the collectivity of the AChE + Cyanidin complex could be associated with the stable binding of Cyanidin, which promotes an increase in the local rigidity of the complex by displacing the intrinsic flexibility of the system to the ends of the chain, being able to drive the conformations of the native structure towards a folded state around the ligand as predicted in the MD analyses, and as predicted by this ENM method, which considers unique characteristics, relating the environment, the prediction of the functional regions of the protein and the modeling of dynamic structures using the coarse-grained approach [42].

ENM also revealed that Cyanidin was the only molecule considered in this study that slightly altered the intrinsic sensitivity (or sensor-like behavior) of the AChE protein system differently from the rest of the compounds. Cyanidin caused the AChE + Cyanidin complex to have a greater propensity to behave as a “receptor” molecule (the greater the magnitude, the more likely a receptor-type response), that is, a greater propensity for residues to receive signals due to the effect of the ligand affecting the native contact time between the residues (Figure 4). The rest of the compounds, although they altered the native response of AChE to perturbations, the measurements were oriented towards the sensor-like behavior of the system, as in the case of the model drug TZ2PA6 (Figure 4).

This behavior can in turn alter the system’s allosteric communication, as shown by the increased number of steps (Communication) required for residues/nodes to transfer interaction-mediated perturbation, as indicated by the stochastic Markov model of information diffusion. This model has been developed to explore communication between residues in proteins and allows the generation of response maps that indicate the regions where residues potentially act as sensors (capable of detecting and responding to changes in the environment) and effectors (participating in the catalysis of chemical reactions, among others) [31,42].

Theoretical methods such as the energy frustration model allow us to locate and quantify the energy frustration present in native protein structures by comparing the contribution to the additional stabilization energy attributed to a given pair of amino acids in the native protein with the energies that would be found through a perturbation (in this case, a ligand) that is capable of creating a different environment for the pair of interacting amino acids in the protein system [30,43]. Based on the energy frustration approaches, a trend similar to that observed in previous deformability analyses was predicted, where Cyanidin was the only one of the compounds whose contribution to folding led to sufficient additional stabilization for normalized native pairs of residues in terms of typical energy fluctuation (according to the global scoring criterion established by the model for minimum frustration) (Figure 5). Cyanidin managed to decrease local interaction and generate an increase in the protein regions called “minimally frustrated” or stabilized [30], as observed in terms of MD and ENM.

2.3. Volumetric Analysis of Internal Cavities

All compounds altered the distribution of AChE’s internal cavities with an average of ≈ 132 ų, which represented a difference of ≈ 5 ų compared to the average volume of AChE’s intrinsic internal cavities after 200 ns MD, with an average minimum of ≈ 40 ų and a maximum of ≈ 528 ų. Cyanidin was the compound that produced the greatest variability in the distribution of internal cavity volume (Figure 6), generating the largest cavities compared to the other compounds, and even exceeding up to 10 times the volume of the largest internal cavity predicted for AChE (Table 4).

These results were further supported by the analysis of the fluctuation of the volume of the internal cavities. ${\langle \delta V_V^2\rangle }^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$, which showed that Cyanidin is capable of promoting an increase of almost double the fluctuation (${\langle \delta V_V^2\rangle }^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$≈ 222 ų) compared to the native volume (${\langle \delta V_V^2\rangle }^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$≈ 117 Å3), followed by Bis-DMM and unlike the rest of the compounds that decreased this fluctuation (Table 4).

In this sense, when analyzing the probability that the formed internal cavities can be modulated by a drug-like molecule (Druggability Probability), it was found that 95% of the internal cavities generated in AChE by the effect of interaction with Cyanidin, followed by Huperzine A and above the probability predicted for the rest of the compounds considered in this study. In fact, it is important to note that the internal cavities with druggability probability by the effect of Cyanidin were not only those with the highest average number of residues in the site, but 75% corresponded to internal cavities with a Druggability Probability = 1.0, followed by 40% for TZ2PA6 (Figure 7; Table 4).

These druggability predictions are important because this parameter is a key factor in determining whether a drug candidate progresses from an early stage to a more advanced stage. Therefore, druggability allows us to represent the ability of a drug to bind to cavities with high affinity [44].

Therefore, since Cyanidin can generate guided pockets by its proximity, as well as non-guided pockets by said proximity according to the considered druggability prediction model [45], it could serve to mediate synergistic inhibitions as reported [46], or to consider transport assays for other ligands by co-pigmentation using compounds even like those considered in this study [5]. These results are relevant because predicting the ability of protein cavities to bind to drugs or drug-like molecules is of great interest in the target identification phase in drug discovery or repurposing [47].

As shown by the dynamic and energy models, as well as the various volumetric analysis strategies considered in this study to predict the intrinsic deformability of this type of protein system, it is evident that studying the intrinsic flexibility of AChE can be used both for system characterization and for the selection of potential inhibitors of interest such as Cyanidin, as has been previously done by applying similar approaches with the compound TZ2PA6 [15,16]. In addition, the volumetric results offer a mechanistic insight as has been suggested [48], and in turn support the observations of experimental reports demonstrating the inhibitory activity of the compounds tested in this study [17,38,39].

3. Materials and Methods

3.1. Selected Molecules and Preparation from Databases for Docking, Theoretical Inhibition, and Molecular Dynamics (MD)

The crystal structure of human acetylcholinesterase (AChE) (EC3.1.17) (PDB: 4EY5) was obtained in PDB format from the RCSB Protein Data Bank (https://www.rcsb.org/). For this study, we considered protein structure chain A and two drugs: the polyphenol Resveratrol (CID: 445154) and the flavonoid Cyanidin (CID: 128861). These compounds have been computationally investigated as potential alternatives for Alzheimer’s treatment due to their ability to interact with and inhibit AChE [5,6,7]. Additionally, we included the alkaloid as bis-(-)-8-demethylmaritidine (CID: 169450457) [39], and Huperzine A (CID: 854026) [38], and the click chemistry compound TZ2PA6 (CID: 5289508) [38] as experimental controls capable of inhibiting AChE in vitro. All compounds were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) in SDF format and converted to PDB format using OpenBabel-3.0 [48]. The Molegro Molecular Viewer (MMV) [49] and BIOVIA Discovery Studio Visualizer [50] packages were utilized for visualization and file preparation.

To explore the interaction between AChE and the compounds, we constructed complexes using directed molecular docking based on the coordinates of the previously reported binding pocket [39]. To carry out this task, we employed the DockThor tool (https://www.dockthor.lncc.br/), a machine learning-based molecular docking algorithm [51]. In order to enhance accuracy, the program was executed with approximately 1 × 10⁶ evaluations per run, while default parameters were utilized for the remaining configuration settings. One of the advantages of using DockThor is that it provides additional steps for the automatic preparation of ligands, including the addition of partial charges, hydrogen atoms, adjustment of rotatable bonds, among other parameters. During the preparation stage, multiple PDB files were generated from the most favorable poses predicted by DockThor, representing AChE with and without ligands bound. During the sampling and ligand pose generation process, the selection criterion was limited to the three most thermodynamically favorable binding poses (with the highest negative magnitudes). This approach allowed us to focus on the most feasible and thermodynamically favorable positions within the complexes.

Based on this criterion, the selected complexes underwent further analysis, including calculations of inhibitory potency and molecular dynamics (MD). The inhibition constant for binding of ligand to proteins ($K_i$) (in units of M) was obtained as:

$$K_i=e^{\left({\displaystyle \frac{{-\Delta G}_{dock}}{RT}}\right)}$$

(1)

The equation provided defines the binding affinity ($\Delta G_{dock}$) in kcal mol^–1, where R represents the universal gas constant, and T represents the absolute temperature. According to this equation, a higher $K_i$ value indicates weaker binding between the inhibitor and the protein, suggesting a higher likelihood of dissociation of the protein-inhibitor complex. To determine the IC₅₀ values (concentration of the inhibitor required to achieve 50% inhibition), we utilized the IC50-to-Ki web tool (https://bioinfo-abcc.ncifcrf.gov/IC50_Ki_Converter/index.php) assuming competitive inhibition as recommended [24]. In this study, we assumed a hypothetical quantitative or stoichiometric substrate-inhibitor ratio of 1:1 in order to simulate the same amount of protein and inhibitors and thus avoid concentration-driven interactions, based on the assumption of structural similarities between the ligand and the substrate as previously suggested [24,52]. In contrast, the experimental inhibitory concentration values were obtained in molar expression, as described later for each of the compounds considered. MD simulations were conducted on the docking hits with two primary objectives: 1) to assess the relative stability of the ligand within the binding pocket, and 2) to explore the conformations with minimum energy and analyze their impact on the thermodynamic and structural stability of the complexes.

For each protein-ligand complex, a series of procedures were carried out to relax the MD system. This process consisted of three phases: relaxation, equilibrium, and sampling, following recommended protocols [24]. The MD simulation took place in a water system with explicit representation, where the system was immersed in a cubic water box with dimensions of 70 Å for each axis (X, Y, Z) and a minimum distance from the protein to the edge of the box at the beginning of the simulation of 8 Å. The system comprised of 8016 water molecules, 32 Na+ ions, and 22 Cl- ions, resulting in a molar density of Na/Cl = 2.76 × 10^–3. Each MD system contained one instance of the protein system and one instance of the docking ligand. The simulation used the Amber99SB-ILDN force field and the TIP3P water model, as widely suggested for the study of AChE [53,54,55,56]. The ligands were automatically described using GAFF from the myPresto package, and both the topology files for the ligands and the proteins were generated automatically using the myPresto package [57]. The system was neutralized, and water molecules were treated as rigid bodies, with periodic boundary conditions applied. The simulation was initiated using the steepest descent algorithm and the conjugate gradient energy minimization method, with position restrictions on the atoms of the protein-ligand complex. An initial 100 ps simulation was conducted, allowing water molecules to diffuse around the protein and reach equilibrium with the protein-ligand system. The PME method was used to calculate the electrostatic contribution to non-bonded interactions, with a cutoff of 14 Å. The NVT ensemble was applied for 100 ps to balance all systems and thermalize them at 300 K, followed by a second run of 100 ps using the NPT ensemble to equilibrate the system at 1 atm and 300 K. The SHAKE algorithm was applied to satisfy link geometry constraints. The production step was performed using the output from the NPT ensemble as the initial configuration for an MD production series at a constant temperature of 300 K, with a total simulation time of 200 ns. Minimum energy structures were obtained every 20 ns as target structures for subsequent analyses. The system was neutralized as suggested to study the ligand-protein systems [24]. In addition, several MD simulation analyses were included, such as RMSD (root mean square deviation), RMSF (root mean square fluctuation), Rg (radius of gyration), and interaction depletion, alongside comparisons to positive controls as recommended [24,58,59]. All MD simulations and additional adjustments were carried out using myPresto package, which can be accessed at https://www.mypresto5.jp/en/.

To evaluate the complexes and calculate the binding free energy ($\Delta G_{bind}$) from the MD trajectories, we employed the widely used molecular dynamics/Poisson-Boltzmann surface area (MM/PBSA) method, which utilizes a thermodynamic integration approach. The final snapshots of the minimized system were subjected to MM/PBSA rescoring, implemented in AMBER22. In line with the aforementioned suggestion, we selected 2500 frames from the last 20 ns of the MD simulation to compute the solvation free energy and molecular mechanics potential energy, which are used to estimate the overall binding affinity of the complexes of interest. The equations to calculate the MM/PBSA binding free energy, as well as the detailed parameters for the MM/PBSA calculation procedures [60], can be found in Supplementary Material. The ΔGMM/PBSA of the complexes was determined based on the frames extracted after the MD simulation. For visualization and analysis, we employed MMV and BIOVIA. For more details on the method, we recommend referring to [24].

3.2. Determination of Conformational Changes of Dynamized Complexes Using Statistical Potentials, Elastic Network Models, and Energy Frustration

3.2.1. Statistical Potentials of Complexes

To determine conformational changes in dynamized complexes, we utilized the CUPSAT (http://cupsat.tu-bs.de/index.jsp) [61] and SWOTein (http://babylone.ulb.ac.be/SWOTein/) packages [62]. These algorithms are inspired by potential statistics and are specifically designed to study protein stability in relation to fluctuations in the energy landscape. CUPSAT (Cologne University Protein Stability Analysis Tool) is a powerful tool used for analyzing and predicting protein stability changes. It employs structural environment-specific atom potentials and torsion angle potentials to predict the unfolding free energy (ΔG) in kcal mol–1 of the complexes. The prediction model of CUPSAT is built upon a dataset of over 4000 non-redundant protein structures with 50% sequence identity. The atomic potentials are derived from a radial pair distribution function and an atomic classification system. Boltzmann energy values are then calculated based on the distribution of radial pairs of amino acid atoms [61]. The details of the equations and theoretical foundations for calculating energy functions of CUPSAT model, which are derived from mean force potentials and are based on the inverse Boltzmann principle (establishing a close relationship between probability densities and energies), are shown in the Supplementary Material.

SWOTein (Strengths and Weaknesses of prOTEINs) is a research tool that investigates how individual residues contribute to the overall folding free energy of proteins. It accomplishes this by utilizing statistical potentials derived from databases that incorporate structural elements, including inter-residue distances. Positive contributions to the free energy are interpreted as stability weaknesses, while negative contributions represent stability strengths. Statistical potentials have proven useful in studying various biophysical properties of proteins, such as protein-protein interactions. The SWOTein model employs knowledge-based mean force potentials derived from datasets of experimentally solved 3D protein structures, employing the principles of the inverse Boltzmann law. The central objective of this model is to establish a relationship between the free energy of a given state and the probability of observing the system in that state. These statistical potentials are based on a coarse-grained representation of protein structures, focusing on main-chain heavy atoms and the amino acid-dependent geometric average center of side chains [62]. The details of the equations and theoretical foundations for calculating energy functions of SWOTein model, are shown in the Supplementary Material.

By combining atomistic MD trajectories and minimum energy structures, we gain a comprehensive understanding of the dynamics and atomic fluctuations within a complex [41]. The low-frequency modes of the Elastic Network Model (ENM) align closely with the dynamic modes or minimum energy structures predicted by MD simulations, exhibiting consistent directions, relative amplitudes of motion, and distances. As a result, several ENM-based approaches enable various applications, including direct MD simulations, and the incorporation of receptor flexibility in docking approaches [63].

3.2.2. Analysis of Intrinsic Deformability Associated with AChE Conformational Changes

In this study, tools based on Elastic Network models (ENM), a technique known as coarse-grained (CG) normal mode analysis (NMA) [64,65], were considered. These tools offer a valuable approach for studying the vibrational dynamics of protein systems near their energy minimum state. In the ENM model, the protein structure is represented using a reduced set of atoms, specifically the Cα atoms, which are treated as nodes. The interactions between pairs of nodes are described by a single-term Hooke harmonic potential. ENM approaches utilize two methods for computing the cross-correlation of atomic motion. The first method, known as the linear cutoff ENM, assumes a force constant for pairwise interactions between Cα atoms within a specified cutoff distance. The second method, referred to as Kovacs-ENM, adjusts the force constant based on the distance between the interacting particles [35,41,66,67,68,69]. The details of the equations and theoretical foundations of ENM, are shown in the Supplementary Material.

Specifically, in this study, we utilized the SPECTRUS (SPECTral-based Rigid Units Subdivision) server (http://spectrus.sissa.it/#home) [40], which performs a decomposition of proteins or protein complexes into quasi-rigid domains. This decomposition is based on analyzing the fluctuations in distances between pairs of amino acids. To compare the functional dynamics of protein complexes with varying degrees of structural similarity, we employed the MD trajectories obtained under the aforementioned conditions. SPECTRUS utilizes an ENM approach, which takes into account the specific properties of each protein complex and its free energy landscape. This allows for reliable reproduction of structural fluctuations [40]. Quasi-rigid domain decomposition methods are based on the concept that, for truly rigid bodies, the distances between any two constituent points remain constant during motion in space [40]. The details of the equations and theoretical foundations of NMA and SPECTRUS, are shown in the Supplementary Material.

In order to further explore and improve the comprehensiveness of our search for potential conformational changes in the complexes, we utilized the DynOmics-ENM server (http://gnm.bahargroup.org/index.php). This server offers a dynamic analysis of biomolecular systems using two popular elastic network models (ENM): Gaussian network model (GNM) and anisotropic network model (ANM). By integrating these models, the server calculates the dynamics of structural coordinates. The DynOmics-ENM method takes into account various unique characteristics, such as the environment, functional site prediction, and the reconstruction of coarse-grained deformed structures. For more detailed information about the method, we recommend referring to [42]. GNM and its extension ANM are elastic network models that operate at a coarse-grained residue level. GNM focuses on predicting the relative magnitudes of fluctuations, while ANM goes a step further by predicting both the magnitudes and directionalities of collective motions. GNM results are generally more robust and are therefore preferred for evaluating square displacements in low-frequency modes. GNM is specifically used to calculate mean-square fluctuations and correlations between residue fluctuations. On the other hand, ANM is utilized to generate conformations that describe residue fluctuations based on the average X-ray structure, specifically in the principal directions of motion [69].

3.2.3. Local Energy Frustration of Complexes

Proteins undergo biologically optimized folding processes to achieve stability and carry out their functions effectively. Consequently, it is expected to observe a certain level of energetic conflict within their local structures. These conflicts enable proteins to explore various conformations within their native set, ultimately contributing to the emergence of “biological function.” In order to gain insights into the configurational patterns, dynamics, and transitions of complex systems, an examination of frustration and its role in the functional dynamics of these complexes was conducted [30]. To quantify the degree of local protein frustration, we utilized the Frustratometer software (http://frustratometro.qb.fcen.uba.ar/) [28,30]. This software is built on the energy landscape theory and provides a model for assessing protein frustration. The Frustratometer employs the Associative Memory, Water Mediated, Structure and Energy Model (AWSEM), a non-additive and coarse-grained force field that can accurately predict protein and complex structures. AWSEM incorporates various factors, including hydrogen bonding and a bioinformatics-based local structure bias term, which considers the effects of neighboring bodies influenced by the local sequence. The local frustration index provided by the Frustratometer measures the contribution of a residue or pair of residues to the overall energy of each structure, particularly when they deviate from their native positions. This index effectively captures the energy variations associated with the molten globule conformations of the polypeptide chain. The local frustration index is particularly valuable for studying tertiary structures and allows for the analysis of energy variations by manipulating parameters such as the distances and densities of the interacting amino acids. For more detailed information about the methodology, we recommend referring to the provided references [28,30]. The details of the equations and theoretical foundations of Frustratometer, are shown in the Supplementary Material.

3.2.4. Volumetric Analysis of Internal Cavities

To determine the changes in the internal cavities of the AChE protein in each complex, the MOLE 2.5 tool (https://mole.upol.cz/) was used. This tool is specialized in the rapid and fully automated localization and characterization of internal cavities in (bio)macromolecular structures [70]. Additionally, the PockDrug tool (https://pockdrug.rpbs.univ-paris-diderot.fr/cgi-bin/index.py?page=Home) was employed to predict pocket druggability or the probability of a drug-like molecule interacting with a specific cavity. This tool provides consistent pharmacological results using different pocket estimation methods. It is robust with respect to pocket boundaries and estimation uncertainties, making it efficient even for pockets that are difficult to estimate. It clearly distinguishes druggable pockets from non-druggable ones using different estimation methods [45]. Besides, the web tool CavityPlus (http://www.pkumdl.cn:8000/cavityplus/index.php#/) was also used for the visualization of the internal cavities. It features a CAVITY module that generates all potential ligand binding sites on protein surfaces. This tool allows assigning different colors to the cavities for easier visual perception. The entire graphical visualization is implemented using the JSmol subprogram on the website [71,72].

On the other hand, knowing that the intrinsic protein flexibility is relevant to the interaction with the ligand and can contribute to the binding process and that the internal cavities play an important role in this structural flexibility and stability [48]; the fluctuation of the cavity volume ${\langle \delta V_V^2\rangle }^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$ is estimated using the following equation:

$${\langle \delta V_V^2\rangle }^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{n=1}^n{\left(V_{V\left(n\right)}-\langle V_V\rangle \right)}^2}}$$

(32)

where, $V_{V\left(n\right)}$ is the cavity volume of each true cavity estimated through MOLE 2.5 tool, $\langle V_V\rangle$ is the mean value of $V_{V\left(n\right)}$ and N the total number of cavities [48].

4. Conclusions

All binding affinity, theoretical inhibitory potency, and dynamics analyses showed thermodynamically favorable and stable binding of the tested compounds against AChE during the considered molecular dynamic simulation time in this study. TZ2PA6 showed the most favorable binding to AChE, followed by the alkaloid Huperzine A and the flavonoid Cyanidin. On the other hand, intrinsic deformability and frustration analyses revealed that all compounds also induced changes in the intrinsic flexibility and rigidity of AChE, with Cyanidin causing the greatest perturbation in the structure of AChE compared to the compounds used as controls due to its in vitro inhibitory activity. Additionally, volumetric analysis of the internal cavities revealed that although all compounds altered the distribution and size of these cavities within AChE, Cyanidin generated the greatest variability in the distribution of the internal cavity volume. Cyanidin significantly increased the fluctuation of the internal cavity volume. These findings suggest that variations in the binding affinity of the compounds can impact AChE differently. On the other hand, Cyanidin could play an important role as a potential inhibitor of AChE, showing a different biophysical and molecular mechanism than other compounds demonstrated as inhibitors, primarily targeting the alteration of the internal cavities of the AChE system. In this sense, Cyanidin could offer a new avenue for future research and therapeutic applications, especially due to its ability to induce conformational changes, particularly at the level of the internal cavities. This behavior is essential, as the ability of protein cavities to bind to drugs is of great interest in the target identification phase of drug discovery.