Allosteric modulation of laeviganoid-based clerodane diterpenes derivatives in muscarinic acetylcholine M1 receptor against tinnitus: a structure-based virtual screening approach

Abstract

Background

Chronic tinnitus is a complication that affects the central nervous system, specifically the auditory cortex, causing a phantom perception of sounds and noises without any external acoustic stimulus. It is more frequent in men than in women and can be caused by excessive exposure to auditory stimuli. The main modulator of auditory functions, particularly in terms of neuroplasticity in the auditory system, is the M1 muscarinic acetylcholine receptor (mAChR M1). In the literature, natural oxygenated heterocyclic compounds have been used to develop drugs that act on the central nervous system (CNS), including clerodane diterpenes. The aim of this study was to evaluate the modulatory action of a series of naturally occurring clerodane diterpenes against chronic tinnitus.

Results

The structure-based virtual screening revealed that Laeviganoid derivatives L1-8 share structural similarities with other oxygenated heterocyclic compounds that modulate mAChR M1. The prediction of pharmacokinetic properties highlighted the L4 derivative as a potential candidate for distribution in the CNS due to its high cell permeability (P_app,A→B = 1.9 × 10⁻⁵ cm/s) and metabolic stability. Molecular docking simulations indicate that the ligand interacts with the active site of mAChR M1 through hydrophobic interactions with residues Tyr106, Trp378, Tyr381 and Tyr404, with an affinity energy of approximately − 8.7 kcal/mol. Molecular dynamics simulations have shown that the L4/M1 complex is stable as a function of time (200 ns).

Conclusion

The in silico results suggest that the L4 can perform allosteric modulation of mAChR M1 in the treatment of tinnitus, as it can bind to the same interaction site as the tiotropium.

Background

Chronic tinnitus is a condition of the human auditory system that causes the perception of sounds and noises (also known as phantom perception) without any external acoustic stimuli. It affects approximately 7 million people worldwide and can be debilitating. It is more common in men [1,2,3,4]. The development of other underlying disorders, including anxiety, stress, and depression, is a worrying aspect of this disorder. These disorders are reported in at least 60% of affected individuals [5, 6]. Additionally, around 20% of the population with tinnitus exhibit suicidal tendencies due to the coexistence of these disorders with this condition [7, 8].

Chronic tinnitus is caused by neuroplasticity in the auditory cortex. This suggests that excessive exposure to acoustic stimuli can cause these changes. Cholinergic receptors are associated with neuroplasticity and are located in both the central nervous system (CNS) and the peripheral nervous system. Recent studies have shown that muscarinic acetylcholine receptors of the M1 type (mAChR M1) are crucial in developing drugs that modulate the activity of the auditory cortex, including those that cause chronic tinnitus [8,9,10].

It is in this context that several classes of oxygenated heterocyclic organic compounds derived from natural products are favoured pharmacophores in the development of new CNS-active compounds [11], including flavonoids [12], chalcones [13] and clerodane diterpenes [14]. Recent studies indicate the efficacy of clerodane diterpenes from Casearia corymbosa in the allosteric modulation of GABA_A receptors, showing anxiolytic activity similar to benzodiazepines [15]. Other studies have shown that this class of compounds inhibits acetylcholinesterase in treating Alzheimer’s disease [16]. Although the effect of clerodane diterpenes on tinnitus has not been elucidated, the neuromodulatory effect of this class of compounds is characteristic.

Clerodane diterpenes, in particular, have a chemical structure rich in sp3 hybridisation atoms (single bonds) distributed between fused cyclic structures, which can be hydroxylated or have oxygenated heterocyclic groups. In the literature, it can be seen that this class of compounds has a rich therapeutic potential against disorders that affect the central nervous system, such as anxiety and depression [15, 17].

In 2023, Wang and colleagues isolated several clerodane diterpenes from the roots of Croton laevigatus [18]. The plant exhibits various biological activities, such as anti-inflammatory [19], cytotoxic [20] and tyrosine phosphatase inhibitory activity [21]. Its heterocyclic phytochemicals consist of a furanyl portion linked to hydroxylated fused bicyclic structures (L1-4 and L8) and a furanyl portion linked to carbonylated fused tricyclic structures (L5-7) (Fig. 1).

Fig. 1

Adapted from Wang et al., (2023) [18]

Two-dimensional representation of the chemical structures of selected Laeviganoid derivatives (L1-8).

In this context, the aim of this study was to characterise the pharmacokinetics and pharmacodynamics of a series of laeviganoid derivatives (L1-8) and to estimate their bioactivity against chronic tinnitus via allosteric modulation of mAChR M1.

Methods

This study employed a virtual target selection screening approach and drug-target interaction analyses, as outlined in Fig. 2. Initially, we selected the eight least polar derivatives of Laeviganoids from Wang et al., (2023) study [18]. These compounds were then plotted on a 2D model and converted into linear SMILES notation for virtual therapeutic target selection screening. A similarity test was conducted with bioactive compounds from the ChEMBL database (more than 300,000 compounds). The compounds were filtered based on their predicted activity in the CNS, determined by lipophilicity and polarity descriptors, as well as their susceptibility to being P-glycoprotein substrates, based on the medicinal chemistry trends of the Pfizer rule [22].

Fig. 2

Methodological design of virtual screening and drug-interaction implemented to discover the anti-Tinnitus potential of a series of derivatives of Laeviganoids 1–8 (L1-8)

To better understand the biological impact of pharmacokinetic properties, the parallel artificial membrane permeability Assay (PAMPA) was used to estimate properties that are usually only tangible in a laboratory environment and with specific cellular languages, including passive cell permeability, clearance rate, and systemic distribution volume [23]. The most favourable ligands were then subjected to molecular docking simulations, where the structure-based virtual screening approach is employed to assess the potential of L1-8 derivatives to modulate M1 muscarinic acetylcholine receptors (mAChR M1), given the neuromodulatory potential elucidated in the literature.

Ligand selection and optimisation

The ligands were chosen from Laeviganoids, the chemical constituents of Croton laevigatus, as identified in Wang et al., (2023) study [18]. These polar furanyl-based clerodane diterpene derivatives (L1-8) were extracted with EtOAc after being soaked in 95% EtOH. The structures were visualised using MarvinSketch® version 23.17.0, an academic licence program by Chemaxon© (https://chemaxon). The ligands’ pre-hydrogenised chemical structure was optimised for structure using the classical molecular mechanics formalism Merk Molecular Force Field 94 (MMFF94) for structural optimisation.

The optimisation process was very strict and returned only the molecule structure with the lowest potential energy. The most stable form of the optimised ligands was then subjected to molecular dynamics simulations with the therapeutic target indicated by the structure-based virtual screening (SBVS) stage. The methodology of da Rocha et al., (2022) [24] was followed for target prediction and prediction of absorption, distribution, metabolism and excretion (ADME) attributes using the Simplified Molecular Input Line Entry System (SMILES) form for structure-based virtual screening (SBVS).

Virtual screening of target selection

The chemical structures of the L1-8 derivatives were submitted to the SwissTargetPrediction online server for SBVS prediction using linear SMILES notation (http://www.swisstargetprediction.ch/). The server is set up to calculate the 3D similarity between two molecules, i and j (as shown in Eq. 1), using x and y vectors for the 20 most stable conformations of bioactive compounds deposited in the ChEMBL database. The shortest Manhattan distance, d, characterises high structural similarity between two compounds (as shown in Eq. 2). The compounds used have known activity from the organism Rattus norvegicus.

$$ \frac{1}{{\left( {1 + \frac{1}{18}d_{ij} } \right)}} $$

(1)

$$ d = \mathop \sum \limits_{s = 1}^{18} \left| {x_{s} - y_{s} } \right| $$

(2)

In silico ADME study

Topological analysis and druglikeness

The chemical structures optimised by the MMFF94 formalism were taken to the molecular lipophilicity potential (MLP) plot in the Python Molecular Graphics (PyMOL) program, as shown in Eq. 3:

$$ MLP = \mathop \sum \limits_{i = 1}^{N} F_{i} \cdot f\left( {d_{ik} } \right) $$

(3)

where F is the lipophilicity contribution of a fragment i, in a molecule with N fragments, and f(d_ik) is the distance function between two atoms i and k in a chemical bond, and the sum includes all structural contributions. The surface map displays a spectrum ranging from polar fragments (red colour spectra) to hydrophobic or apolar fragments (blue colour spectra). The plots were analysed in relation to the descriptors of lipophilicity (AlogP), topological polar surface area (TPSA), and molecular weight (MW), and were evaluated according to the biopharmaceutical classification criteria of the Pfizer, Inc. rule and the GlaxoSmithKline (GSK) rule, constituting the druglikeness analysis stage [25, 26].

Predicting PAMPA descriptors

To corroborate the empirical druglikeness estimates, the molecules in Simplified Molecular Input Line Entry System (SMILES) notations were submitted to the prediction of Parallel Artificial Membrane Permeability Assay (PAMPA) descriptors using the ADMETlab 2.0 (https://admetmesh.scbdd.com/), AI Drug Lab (https://ai-druglab.smu.edu/admet) and SwissADME (http://www.swissadme.ch/) web servers, through the methodology of da Rocha et al., (2022) [24].

The predicted descriptors include in vitro absorption, distribution, metabolism and excretion (ADME) properties such as apparent permeability (P_app,A→B) in Madin-Darby Canine Kidney (MDCK) cell line, P-glycoprotein (P-gp) substrate, intrinsic clearance in the human liver microsome (CL_int,u), volume of distribution of steady-state (Vd_ss), Plasma Protein Binding (PPB), Ames mutagenicity and drug-induced liver injury (DILI).

Site of metabolism prediction

The site of metabolism prediction protocol followed a similarity test of the molecular fragments read by SMILES with biomolecules with known biotransformation sites for different isoforms of cytochrome P450 (CYP450) in the human liver microsome (phase I) [27]. The test considers the alignment between structural specificity and the sensitivity of the molecular fragment being a site of metabolism. To do this, the SMILES code was read from the XenoSite online server (https://xenosite.org/), and the results were related to the metabolism descriptors, which include substrates and inhibitors of CYP450 types, and toxicity from the ADME predictive test.

Molecular docking simulations

To understand the biochemical interactions of the L1-8 derivatives with the M1 muscarinic acetylcholine receptor (mAChR M1) biological receptor identified by virtual screening in the target prediction phase, a series of molecular docking simulations was performed to verify the structural contributions that potentiate the target’s modulatory effect.

Thus, the ‘Structure of the human M1 muscarinic acetylcholine receptor bound to antagonist Tiotropium’ was obtained from the RCSB Protein Data Bank server (https://www.rcsb.org/) based on the methodology of da Rocha et al. (2023) [28]. The membrane receptor in the Homo sapiens organism and Spodoptera frugiperda expression system has been classified and deposited under PDB code ID: 5CXV. The receptor’s allosteric modulation control ligand is the benzodiazepine bromazepam (BZP), as determined by the molecular docking protocol of Bojić et al. (2017) [10]. The protein preparation step involved removing water molecules (H₂O) and co-crystallised molecules, such as cholesterol hemisuccinate, triethylene glycol, glycerol, 1,2-ethanediol, and the inhibitor tiotropium (TTP), using the UCSF Chimera® program [29].

The AutoDockTools® program (https://autodocksuite.scripps.edu/adt/) was used to add Gasteiger charges and configure the grid box to cover the entire conformational space of the protein. The dimensions of the grid box were x = 65 Å, y = 82 Å and z = 100 Å, adjusted to the coordinates x = -15.140 Å, y =  − 16.115 Å and z = 61.145 Å. The AutoDockVina® code [30] was configured to perform 50 independent simulations, each with 25 poses. The best pose selection criterion was based on the formation of the ligand–receptor complex with an affinity energy of less than −6.0 kcal/mol, within a statistical validation interval of root-mean-square deviation (RMSD) of less than 2.0 Å [31].

Molecular dynamics simulations

Molecular dynamics (MD) simulations were conducted to evaluate the stability of the ligand–protein complex formed by the ligand exhibiting the highest degree of alignment with the receptor in the molecular docking simulations about the mAChR M1 receptor. The simulations entail predicting the relative dynamics of each atom in the system comprising the ligand–protein complex over a specified period by convergence criteria based on variations in RMSD, affinity energy, and the formation of favourable interactions at the site of interest [32].

In order to conduct the MD simulations, the complexes formed by the best poses were employed in the molecular docking simulations (ligand–protein complex). For each system, water molecules were incorporated by the TIP3P format within a cubic box that encompasses the entire complex [33, 34]. After adding water molecules, the system must undergo equilibration by introducing Na⁺ and Cl⁻ ions [35]. Following the equilibration and solvation stages, the systems were configured to act under the influence of the CHARMM36 force field [36, 37], utilising GROMACS 2020.4 software for the preparation stages and prediction calculations (https://www.gromacs.org/).

For more refined criteria for preparing MD systems, two parameters have been standardised, which are closely related to the movement of the atoms present in the system, while the other parameter is related to the MD simulation time; the first parameter works on the temperature at which the system will be under the influence, with a value set equal to 310 k and consequently kept constant until the end of the simulation, for this criterion the V-rescale integrator was used [38], the second parameter used is the adjustment of the internal pressure that the system will be under the influence. The second parameter used is to set the internal pressure equal to 1 bar, as indicated by the Parrinello-Rahman method [39]. About the time that the simulations will be carried out, a time of 200 ns was set, with the appropriate use of the leapfrog integrator [40, 41], in order to present and highlight the similarity resulting from the simulations and subsequently deduce the coherence of the data predicted in the simulation.

Once the MD simulations have been carried out, the root-mean-square deviation for molecular dynamics (RMSD_MD) values are presented and inserted in Eq. 5. For this study; this equation highlights the use of the x, y, and z coordinate axes associated with a pair of molecules inserted in the simulated system. The equation presents the values of the axes (x, y, and z), deriving the use of all mathematical operations to determine the RMSD values, referring to each estimated system [42].

$$ RMSD_{MD} = \sqrt {\frac{1}{n}\mathop \sum \limits_{i = 1}^{N} \left( {x_{i1} - x_{i2} } \right)^{2} + \left( {y_{i1} - y_{i2} } \right)^{2} + \left( {z_{i1} - z_{i2} } \right)^{2} } $$

(5)

Still referring to the MD simulations, the molecular mechanics of generalised born surface area (MM/GBSA) values were simulated, highlighted in Eq. 6, where its values consist of the use of the initial files of the MD systems in comparison with the simulation trajectory file, the use of this equation is intended to point out the variations in free binding energy (ΔG_bind), of the systems formed with complex (ligand and protein), the energy value shows an association with the affinity resulting from the interaction between ligand and protein, and thus pointing out which system shows a more excellent favour about the interaction potential [43, 44].

$$ \Delta G_{{{\text{bind}}}} = E_{{{\text{vdW}}}} + E_{{{\text{ele}}}} + G_{{{\text{GB}}}} + G_{{{\text{SA}}}} - T\Delta S $$

(6)

The estimation of the energy values has been influenced by van der Waals interactions (E_vdW), which are associated with electrostatic energy contributions (E_ele), taking into account the polar effects in the complexes formed between the ligand and the protein (G_GB). The term G_SA represents the non-polar components of these complexes, while TΔS corresponds to the product of the absolute temperature and the entropy variation of the complex formed.

Results

Virtual screening of target selection

The possible therapeutic action route of the L1-8 derivatives based on their chemical structure was identified through the application of the structure-based virtual screening (SBVS) protocol (Fig. 3). The ligands were observed to exhibit between 40 and 50% of their bioactivity as ligands for GPCRs, also known as membrane receptors, and proteins of the enzyme class, with the exception of the L5-6 derivatives, where this range can expand to 60–70% (Fig. 3A). The ligands can bind to drug transporters such as P-gp (in cell membrane) and plasma proteins (in systemic circulation), as shown by the blue colour bar in the graph in Fig. 3A. It is important to maintain objectivity and avoid subjective evaluations.

Fig. 3

Virtual structure-based screening results expressed as A percentage bioactivity of each L1-8 derivative by biological target class and B 3D similarity with known bioactive compounds from the ChEMBL database for identification of the specific biological target, where the colour spectrum ranges from red (few or no similar compounds) to blue (up to 240 similar bioactive molecules)

The ligands exhibited a high structural similarity with other mAChR M1 ligands in the ChEMBL database, particularly with oxygenated heterocyclic compounds. This class includes the L1-8 Laeviganoids (Fig. 3B), with the L1, L4 and L8 derivatives standing out due to their structural similarity with at least 114 compounds in the database that are mAChR M1 ligands. These derivatives are strong candidates against Tinnitus (Table 1). Furthermore, the compounds exhibited structural similarity with at least 32 bioactive compounds found in the ChEMBL database that serve as substrates for the CYP450 isoform type 17A1 (CYP17A1) during phase I metabolism (refer to Table 1). The heatmap in Fig. 3B shows that the majority of the compounds (represented by white to blue spectra) undergo biotransformation in the human liver microsome system, except for derivatives L2 and L7 (represented by red spectra), indicating that these substances are not biotransformed in pre-systemic metabolism.

Table 1 Number of similar compounds (known 3D structures) to the L1-8 derivatives for different specific targets with a considerable number of similar compounds

Topological analysis and druglikeness

In the topological analysis of the MLP shown in Fig. 4, it is evident that the polar oxygenated H-bond donor and acceptor groups (OH and O groups) have an impact on the molecule’s surface accessibility to both aqueous and organic environments of the L1-8 derivatives. In this analysis, it was observed that the L4 ligand had the highest AlogP value of 4.62 (Table 2). This was due to its monohydroxylated fused bicyclic structure (blue to green colour spectra) and polar H-bond acceptor contribution in the portion of the carbonyl linked to the furan (yellow to red colour spectra), compared to the dihydroxylated L2 derivative.

Fig. 4

Ligands L1-8 with their molecular lipophilicity potentials (MLP) plotted and analysed in terms of the polarity of the molecular fragments, where the colour spectrum ranges from red (polar groups) to blue (hydrophobic or apolar groups), resulting in the calculated values of AlogP

Table 2 Physicochemical properties of lipophilicity (AlogP), molecular weight (MW in g/mol) and polarity (TPSA in Ų) calculated and applied to the GSK and Pfizer biopharmaceutical classification systems

On the other hand, the derivatives L5-6 exhibited an AlogP value of 2.91 (Table 2), indicating moderate lipophilicity. The difference in the position of the substituted (R) and (S)-OH groups (Fig. 1) did not affect the distribution of the hydrophilic surface. The hydrophilic surface showed a strong contribution from the carbonyl groups of the fused tricyclic portion of these derivatives (red colour spectra), as shown in Fig. 4.

When considering the polarity values (TPSA), it is evident that the lipophilicity of the ligands has a significant impact on the distribution of the compounds in the CNS [45]. The L4 derivative, with a TPSA of approximately 50, is the least polar. The value of 44 Ų reflects the contribution of the polar surfaces of the carbonyl portion conjugated to the furan ring. The H-bond acceptor groups of the C = O and R-O-R type are less polar than OH-type donor groups [46]. This results in a low susceptibility to being a P-gp substrate, as indicated by the red colour in Fig. 5. Therefore, this compound is more effective at the CNS than the other analogues [47]. Compound L4 scored positively ( +) for the Pfizer rule, indicating that it resides in a physicochemical space formed by compounds of high lipophilicity (AlogP > 3) and low polar surface area (TPSA < 75 Ų), suggesting distribution to the CNS [25]. However, the GSK rule raises a structural concern regarding its lipophilicity, which may limit its pharmacokinetics (Table 2). Furthermore, all derivatives possess an optimal polarity for effective gastrointestinal absorption due to the presence of oxygenated groups such as R-OH and R₂O [48]. The estimated bioavailability (%F) is unlikely to be less than 30% as shown in Table 3.

Fig. 5

Alignment graph between lipophilicity (AlogP) and polar surface (TPSA) for the estimation of the BBB activity, where the red region is the physicochemical space of access to the CNS, formed by compounds with high AlogP and low TPSA, and the region outlined up to TPSA = 75 Ų is the region CNS safety. The blue indicators are compounds whose presence of OH groups results in a polarity that makes the molecules susceptible to efflux by P-gp, while the red indicators highlight passively permeable molecules

Table 3 Pharmacokinetic descriptors predicted by the consensus test of ADME properties between the SwissADME, ADMETlab 2.0 and AI Drug Lab databases

Although this step constitutes an empirical prediction of druglikeness, pharmacokinetic descriptors were provided to corroborate the results, thereby ensuring greater alignment between the predictive techniques. The results of the pharmacokinetic prediction demonstrate a strong correlation with the physicochemical attributes evaluated in this section, which can be seen in Table 3.

Predicting ADME properties

The predictive analysis of ADME and the topological analyses of druglikeness based on the chemical structure of the L1-8 derivatives are in agreement. It was observed that compounds located in a physicochemical space defined by AlogP > 3 and TPSA ≤ 75 Ų may exhibit activity in the CNS. Here, compound L4 is notable for its low polarity, which results in greater cellular permeability compared to efflux.

In this way, the parallel artificial membrane permeability assay (PAMPA) descriptors were employed to estimate the cell permeability potential of the L1-8 derivatives [23, 49]. The results demonstrated that the pharmacokinetic properties align with the physicochemical and structural trends observed in the compounds. This supports the idea that the compound is unlikely to be a P-gp substrate, as well as the P_app,A→B descriptor of 1.9 × 10^–5 cm/s, this suggests that the compound can penetrate more selective cell membranes, such as the MDCK cell model, commonly employed to estimate the permeability of active drug candidates in the central nervous system (CNS). The L2, L6, and L8 compounds exhibited P_app,A→B MDCK values exceeding 2.0 × 10⁻⁵ cm/s, yet they may function as P-gp substrates. This suggests they may undergo passive efflux to the extracellular environment, potentially reducing the viability of small molecules acting in the CNS. Although derivatives L1 and L7 exhibited a high predicted apparent permeability coefficient (P_app,A→B MDCK in the order of 5.6 and 2.1 × 10⁻⁵ cm/s, respectively), the compounds may be highly polar as a consequence of the substituted carboxyl-based groups (R-COO-R) in the Laeviganoid substructure (Fig. 5).

Furthermore, it was observed that the compounds exhibited a lipophilicity range indicative of high hydrophobicity, which may facilitate the formation of chemical entities widely distributed in the organic phase of physiological environments. It was observed that compound L4, likely to be active in the CNS without being a P-gp substrate, exhibited a predicted relative Vdss value of 4.91 L/kg. This indicates that the high lipophilicity of the compound facilitates a broader distribution of the bioavailable molecular fraction in biological tissues than in blood plasma, thereby corroborating the likelihood of the substance’s access to the CNS, where the blood–brain barrier (BBB) is a highly selective biological membrane. Moreover, all the compounds exhibited Vdss values between 3.0 and 6.7 L/kg, indicating their high lipophilicity.

These findings are by the PPB prediction, whereby the estimated values below 55% indicate that a substantial proportion of the molecular fraction is not bound to the serum proteins of the blood plasma, thereby allowing for a more widespread distribution and the potential for a more significant biological effect. Furthermore, the predicted clearance rate (CL_int,u) of 7.70 mL/min/kg is a significant factor influencing the oral bioavailability of the L4 derivative. This suggests that the compound may exhibit metabolic stability at the level of hepatic microsomes, resulting in a high concentration of the non-biotransformed molecular fraction in the systemic circulation (Table 3). The compound’s metabolic stability is directly related to the formation of low-reactive metabolites in the human liver microsome system. This can be observed in the section on the site of metabolism (Fig. 6).

Fig. 6

Prediction of the site of metabolism of the derivatives A L1, B L2, C L3 and D L4, considering the probability that they are substrates of the CYP450 isoforms in phase I metabolism

Site of metabolism prediction

The test for structure-dependent metabolism site prediction allowed for the estimation of potential secondary metabolites resulting from the oral administration of L1-8 derivatives. This test aligns the descriptors of structural specificity and sensitivity of the fragments with CYP450 substrate isoforms [27].

The observation made here is that the furan ring has unsaturated sites that are susceptible to aromatic hydroxylation (Fig. 6). This process can lead to the formation of reactive epoxide intermediates that are capable of intercalating into protein and DNA structures, as shown in Fig. 6. Notably, the ligands L3 (Fig. 6C) and L4 (Fig. 6D) feature an isolated alkene within the bicyclic ring of the Laeviganoids substructure. The presence of an aliphatic hydroxyl group at this position results in the formation of hydroxylated metabolites that exhibit reduced reactivity compared to the epoxides derived from the hydroxylation of the furan ring. This finding supports the hypothesis that the compounds in question have a low probability of causing liver damage or the formation of reactive metabolites that could interact with proteins and DNA (Table 3).

These fragments were identified through a test for detecting substrates of CYP450 isoforms in pre-systemic metabolism (Phase I). This test is based on structure–activity relationship (SAR) descriptors, which relate the structural specificity and sensitivity of the molecular fragment to its capacity to act as a CYP450 substrate. These are derivatives that can be derived from specific isoforms, such as CYP17A1, identified in the virtual screening of the target prediction stage. Consequently, they are processes that can reduce the systemic bioavailability of these substances. However, they do not negatively affect the order of hepatic clearance, leading to the formation of chemical species of low toxicity by metabolic activation (Table 3).

Molecular docking against mAChR M1

At the end of the cycle of 50 independent molecular docking simulations, each consisting of 20 poses for each of the ligands, it was observed that all the ligands achieved a statistical adjustment of RMSD of less than 2.0 Å. This indicates optimal specificity for the amino acid residues of their binding sites on the mAChR M1 (Fig. 7A). The redocking carried out with the inhibitor TTP achieved an RMSD of around 0.99 Å, which was used as the statistical standard for the other simulations (orange bar). The L4 ligand formed a ligand–receptor complex with mAChR M1 with an affinity energy of approximately −8.7 kcal/mol, which meets the ideality standard (affinity energy < −6.0 kcal/mol) and provides more favourable energy conditions than the TTP inhibitor (−7.157 kcal/mol) (Fig. 7B).

Fig. 7

Data from the molecular docking simulations expressed as A RMSD in Å and B affinity energy in kcal/mol, for the L1-8 derivatives and the TTP inhibitor against the mAChR M1

This affinity energy analysis is of great biological importance in understanding the biochemical interactions of small ligands with receptors expressed in the CNS. According to Shityakov et al., (2014) [50], the formation of ligand–protein complexes with an affinity energy lower than −6.0 kcal/mol is associated with CNS-active compounds that readily penetrate the BBB. This analysis suggests that compound L4 is a lead compound due to its CNS activity characteristics and affinity energy against the mAChR M1 receptor (Fig. 7B).

When analysing the ligand–receptor interactions, it was possible to observe that the L4 ligand complexed to the mAChR M1 interacting with amino acid residues in common with the TTP inhibitor, revealing that the compound binds to the same interaction site as the co-crystallised inhibitor (Fig. 8A). Here, it was possible to observe that the L4 ligand showed hydrophobic interactions in common with the TTP inhibitor, including interactions with the aromatic side chains of the Tyr106, Trp378, Tyr381 and Tyr404 residues (Fig. 8B), where the calculated distances between the ligand and these amino acid residues are around 3.1–3.4 Å (Table 4), characterising interactions of moderate strength (blue to white spectra in Fig. 8C) [51]. It is interesting to note that the inhibitor interacts with the aromatic centre of the Tyr106 residue by hydrophobic interaction with a contribution from the π electrons of its thiophenyl portion (Fig. 8D), while the L4 derivative acts as an H-bond acceptor for the hydroxyl of the residue through its furanyl portion (Fig. 8E).

Fig. 8

A Three-dimensional representation of mAChR M1 with the ligands TTP and L4 in their binding sites, B interactions between the ligands and the residues of the receptor’s active site, C heatmap relating the distances between the ligands and the amino acid residues to the interaction forces, where the colour spectrum ranges from blue (strong interactions) to red (weak interactions) and two-dimensional representations of the structural contributions of the ligands D TTP and E L4 in the interactions with the residues of the receptor’s active site

Table 4 General data on the interactions between the L1-8 derivatives and the mAChR M1, compared to the TTP and BZP controls, expressed as types of interaction, interaction residues and ligand–receptor distance involved in each type of interaction characterised in the molecular docking simulations

On the other hand, the other ligands (L1-3, and L5-8) bound to mAChR M1 at a site believed to be the allosteric modulation site of the benzodiazepine BZP. They interact mainly with the aromatic residues Tyr82, Thr189, Tyr381, Trp400 and Tyr404 (Fig. 9A). Here, interactions of similar strength between ligands L1-3, L5-6, and L8 with the aromatic side chain of the Tyr404 residue are highlighted. The distances calculated between 3.0 and 4.0 Å reveal interactions of moderate to weak strength (clear blue to white spectra in Fig. 9B).

Fig. 9

A Three-dimensional model of the coupling of the L1-3 and L5-8 derivatives to the BZP allosteric modulation site in relation to the TTP inhibitor B heatmap relating the distances between the ligands and the amino acid residues to the interaction forces, where the colour spectrum varies from blue (strong interactions) to red (weak interactions)

Molecular dynamics simulations

RMSD analysis

Following an analysis of the Molecular Dynamics (MD) simulations, Fig. 10 presents the results of the root-mean-square deviation (RMSD) variations of the simulated systems, with duplicates classified as Run 1 (black) and Run 2 (red).

Fig. 10

RMSD_MD variations of the simulated systems A M1 receptor B System formed in the presence of the inhibitor TTP C System composed of the ligand L4

Figure 10A shows the variations of the biological receptor without ligands in its cavities. The RMSD data show that the two simulations are similar. In the first simulation (black), between 0 and 17 ns, the RMSD variations were around 2.5 Å. From 18 ns onwards, the RMSD values of both simulations remained between 2.8 Å and 3.1 Å until the end of the simulation. This behaviour indicates that the receiver showed minor structural deformation. In contrast, in the second run associated with the M1 receptor, there was a more pronounced structural variation with values of 4.6 Å at 51 ns and a further deformation of 4.3 Å observed at approximately 73 ns. After this point, the system showed minor structural variations and remained cohesive until the end of the simulation, with values of around 3.0 Å.

The variations of the system with the TTP inhibitor associated with the M1 receptor are shown in Fig. 10B. Both runs showed values of around 2.5 Å up to 20 ns. After this period, up to 38 ns, the first run increased structural deformation, reaching 4.2 Å, while the second run maintained values close to 2.6 Å. After 40 ns, the RMSD values of the two runs remained relatively constant and similar, ending the simulation with values of 2.6 Å and 3.0 Å for Run 1 and Run 2, respectively.

In Fig. 10C, the results of the simulations with the L4 ligand inserted show that both runs showed similar behaviour in terms of structural deformation of the system. In the first run, the RMSD variations remained around 2.8 Å up to 40 ns and remained constant until the end of the simulation, reaching around 4.0 Å at 200 ns. The second run showed pronounced oscillations, reaching values of up to 4.2 Å at certain times and remaining consistent until 88 ns, with deformations of around 3.4 Å. Between 90 and 100 ns, the second run showed pronounced deformation, reaching 4.3 Å. After 100 ns, the complex with L4 showed less conformational variation, remaining relatively stable until the end of the simulation, with values around 3.8 Å.

These results show the consistency and relative stability of the MD simulations for the different systems analysed, indicating the conformational variations with and without the presence of the ligands.

According to the data from the MD simulations, the system with the M1 receptor alone showed minor conformational variations, reaching RMSD values of up to 4.6 Å, indicating low structural variation. When the inhibitor TTP was added to the M1 receptor, the results showed a progressive increase and more coherent values over the 200 ns of the simulation. After approximately 140 ns, the system showed stability with less structural deformation, maintaining values close to 4.2 Å.

In the simulation with the L4 ligand, the structural variations were similar to those observed with the TTP inhibitor, reaching values of up to 4.3 Å. These data suggest that the L4 ligand can induce a structural deformation in the M1 receptor similar in magnitude to that caused by the TTP inhibitor. Thus, the presence of the L4 ligand appears to play a comparable role in the structural stability of the receptor, as indicated by the analysis of the RMSD deviations.

Hydrogen bond analysis

The results of the MD prediction simulations were used to evaluate the percentage of hydrogen bond formation between the initial (0 ns) and final (200 ns) states of the simulation. Interactions with amino acid residues in the inhibition site (Asp105, Tyr106, Trp157, Ala196, Val113, Phe197, Asn382, and Trp378) were considered as criteria for analysing favourable interactions to determine an interaction profile and estimate the viability of the compounds analysed. Only contributions greater than or equal to 5% of the total simulation were included for the analysis of hydrogen bond formation.

Figure 11A shows the percentage of hydrogen bonds formed between the TTP ligand and the M1 receptor. The results show that throughout the simulation, TTP interacted with 13 amino acid residues in the receptor, the most important being Tyr404 (55.34%), Asn382 (42.61%), Cys407 (29.72%), Tyr106 (17.23%), and Ala196 (13.24%). Other residues were Phe197 (8.49%), Tyr408 (8.54%), Tyr381 (12.79%), Asp105 (5.04%), Gln110 (5.45%), Trp378 (8.19%), Trp157 (6.54%) and Thr192 (6.99%).

Fig. 11

Percentage of hydrogen bond formation A system composed of the TTP/M1 complex B system formed in the presence of the L4/M1 complex

In Fig. 11B, the data for the system with the ligand L4 associated with the M1 receptor show that L4 formed hydrogen bonds with ten amino acid residues, most frequently with residues Tyr106 (71.03%), Trp378 (35.57%), Ala196 (19.03%) and Tyr381 (16.99%). Other residues involved were Ser109 (15.93%), Asn382 (13.64%), Phe197 (11.29%), Thr192 (6.95%), Leu183 (6.75%) and Trp157 (5.85%).

These results provide an overview of the hydrogen interactions between the ligands and the M1 receptor and highlight the amino acid residues that play a crucial role in stabilising the TTP/M1 and L4/M1 complexes throughout the simulation.

The molecular dynamics simulations indicated that both ligands formed hydrogen bonds with residues in the inhibition site of the target receptor. The TTP ligand established hydrogen bonds with seven amino acid residues in this site (Asp105, Tyr106, Trp157, Ala196, Phe197, Asn382 and Trp378), with binding percentages ranging from 5.04 to 42.61%. These findings indicate that TTP exhibits high specificity for the receptor cavity and maintains its presence in the interaction site as predicted by molecular docking simulations, thereby substantiating its potential as an inhibitor.

The L4 ligand demonstrated the formation of hydrogen bonds with five amino acid residues within the inhibition site (Trp157, Ala196, Asn382, Tyr106 and Phe197), with frequencies ranging from 5.85 to 71.03%. These values indicate that L4 also exhibits high specificity for the predicted cavity, thereby substantiating its viability as an M1 receptor inhibitor, in a manner analogous to TTP.

These findings underscore the potential of both ligands as M1 receptor inhibitors, exhibiting structural specificities that align with the predictions derived from molecular docking.

MM/GBSA analysis

The results of the molecular dynamics simulations indicated which complex systems (ligand and receptor) exhibited the most favourable free energy indices, as determined by ΔG_bind (kcal/mol) through the MM/GBSA calculations. These findings are detailed in Table 5. The free energy values for the TTP/M1 and L4/M1 systems were -10.72 ± 2.78 kcal/mol and −26.93 ± 2.96 kcal/mol, respectively. In the TTP/M1 system, the terms with the greatest contributions to the free energy of binding were E_vdW and Eele, with values of −53.90 kcal/mol and 41.79 kcal/mol, respectively. In the L4/M1 system, the most significant contributions came from EvdW and GGB, with values of −38.95 kcal/mol and 21.26 kcal/mol, respectively.

Table 5 Free energy values associated with the TTP/M1 and L4/M1 complexes

These values highlight the most relevant interactions and suggest greater stability for the L4/M1 system in comparison to TTP/M1, as evidenced by the differences in free energy of binding between the complexes.

The TTP/M1 complex exhibited three favourable contributions to the free energy of binding, derived from the E_vdW, G_GB and G_SA terms. Conversely, the E_ele and -TΔS terms were identified as unfavourable contributions. For the L4/M1 complex, the E_vdW, E_ele and G_SA terms enhanced the affinity energy, whereas the G_GB and -TΔS terms exerted unfavourable effects.

The results of the molecular dynamics simulations indicate that both systems are viable for interaction, with binding free energy values of −10.72 ± 2.78 kcal/mol for TTP/M1 and −26.93 ± 2.96 kcal/mol for L4/M1. A more detailed analysis revealed that the L4/M1 complex exhibited greater viability of interaction due to its lower binding free energy value, indicating enhanced stability for this system.

Discussion

The SBVS is a prediction technique that aids in the early stages of drug discovery. It estimates a biological mechanism of action based on the molecular structure of small bioactive molecules [52, 53]. The SBVS algorithm employed in this study estimated the bioactivity of a substance based on the structure of the input compound. The final 3D similarity with over 300,000 compounds with known biological activity deposited in the ChEMBL database was used to select the biological target for the pharmacodynamics simulations [54, 55].

In this study, it was possible to observe that the L1-8 derivatives present a specific molecular structure similar to a series of mAChR M1 allosteric modulators, including heterocyclic compounds based on oxygenated organic functions. This finding suggests a mechanism of action susceptible to pharmacodynamic evaluation by molecular docking simulations, as well as an evaluation of ADME properties convenient to the CNS [28].

Characterising the pharmacokinetics and bioactivity of these drug candidates is crucial to the pharmaceutical industry. This is because around 50% of drug failures in clinical trials are due to pharmacokinetics and toxic effects resulting from undesirable interactions [56]. Systematic analysis of physicochemical descriptors, aided by computers, can predict the in vitro pharmacokinetic spectrum and bioactivity of new drug candidates, such as enzyme inhibitors and G-protein coupled receptor (GPCR) ligands [57]. Wager and coworkers (2016) [22] observed that among a group of drugs, drug candidates, and drugs in experimental stages, substances residing in a physicochemical space formed by high lipophilicity and low polarity (AlogP > 3 and TPSA ≤ 75 Ų) showed CNS activity with high apparent cell permeability (P_app,A→B > 10 × 10⁻⁶ cm/s), low passive efflux by P-gp, and low hepatic clearance rate (CL_int,u < 100 mL/min/kg). The alignment of these attributes ensures high oral absorption and good metabolic stability. This enables the selection of new drug candidates with optimal cell viability [58, 59]. The study shows that the derivatives meet Pfizer, Inc. biopharmaceutical classification criteria. However, due to the low polarity of the L4 derivative, it is not a substrate of P-gp. This suggests that the bioavailable molecular fraction of this substance is more widely distributed in the CNS compared to the other derivatives.

The results of the topological analyses and the analysis of the ADME descriptors corroborate the analyses of the molecular docking simulations, where the L4 ligand stands out for its favourable alignment between the statistical and energetic parameters when interacting with the mAChR M1.

Recent studies on molecular recognition (structure-based virtual screening) have revealed the relationship between heterocyclic organic compounds with organic functions and the modulation of cholinergic receptors. These compounds act as agonists in different types of receptors in the CNS [60,61,62]. It can be observed that Laeviganoids derivatives are based on oxygenated clerodane diterpenes that contain a furanyl portion. This reveals the high similarity that these compounds share with other oxygenated organic compounds that modulate mAChR M1. Molecular docking simulations confirmed that the L4 derivative interacts with the TTP inhibitor through the aromatic side chains of the Tyr106 and Tyr404 residues [63]. Additionally, the L4 derivative forms a hydrophobic interaction with Ser109 and an H-bond interaction with the hydroxyl group of the side chain of the Tyr381 residue, with a contribution from the oxygenated group of the furanyl portion. The finding suggests that the L4 derivative may be involved in the allosteric modulation of mAChR M1 in isolation. This makes it a promising active principle for the treatment of Tinnitus.

The binding free energy data (ΔG_bind) obtained from the molecular dynamics simulations indicates a correlation between the energy parameters and the conformational stability of the complexes formed. It can be observed that the L4/M1 complex exhibits a greater capacity for interaction in comparison to TTP/M1, as evidenced by the more negative ΔG_bind value. This behaviour may be related to the structural profile of the L4 ligand, which appears to have a higher affinity for the key residues of the M1 receptor inhibition site. This is evidenced by the more frequent interactions with specific residues, such as Tyr106 and Phe197 [63]. Therefore, the presence of L4 in the complex appears to enhance the stability of the system, as evidenced by both the E_vdW values and the impact of E_ele, which may indicate a heightened inhibitory potential for the M1 receptor. This distinction in energetic interactions lends support to the proposition that structural and electronic variations between ligands exert a considerable influence on binding specificity and affinity, thereby underscoring the significance of adapting ligand properties to enhance efficacy and selectivity in receptor-ligand complexes.

Conclusion

In conclusion, the structure-based virtual screening identified mAChR M1 as the target based on the unique structural features of laeviganoid derivatives L1-8. Among these, predictive ADME testing indicated L4 as the most favourable for CNS-targeted action due to its alignment of high lipophilicity and low polarity, which enhanced cellular permeability (high P_app,A→B) and facilitated BBB penetration. Molecular docking simulations further revealed that only L4 effectively interacted with active site residues, suggesting its potential as a promising anti-tinnitus agent.

Molecular dynamics simulations demonstrated that both TTP and L4 ligands caused minimal structural deformation to the receptor, with consistent results between the two. In terms of hydrogen bonding, both ligands formed interactions with residues within the receptor’s inhibition site, thereby supporting the viability of both compounds as potential mAChR M1 binders. However, MM/GBSA energy evaluations showed that L4 had a lower binding free energy, indicating a more favourable interaction compared to TTP.