Computational methods to analyze and predict the binding mode of inhibitors targeting both human and mushroom tyrosinase

Tyrosinase, a copper-containing enzyme critical in melanin biosynthesis, is a key drug target for hyperpig-mentation and melanoma in humans. Testing the inhibitory effects of compounds using tyrosinase from Agaricus bisporus (AbTYR) has been a common practice to identify potential therapeutics from synthetic and natural sources. However, structural diversity among human tyrosinase (hTYR) and AbTYR presents a challenge in developing drugs that are therapeutically effective. In this study, we combined retrospective and computational analyses with experimental data to provide insights into the development of new inhibitors targeting both hTYR and AbTYR. We observed contrasting effects of Thiamidol ™ and our 4-(4-hydroxyphenyl)piperazin-1-yl-deriv-ative ( 6 ) on both enzymes; based on this finding, we aimed to investigate their binding modes in hTYR and AbTYR to identify residues that significantly improve affinity. All the information led to the discovery of compound [4-(4-hydroxyphenyl)piperazin-1-yl](2-methoxyphenyl)methanone (MehT-3, 7 ), which showed comparable activity on AbTYR (IC 50 = 3.52 μ M) and hTYR (IC 50 = 5.4 μ M). Based on these achievements we propose the exploitation of our computational results to provide relevant structural information for the development of newer dual-targeting molecules, which could be preliminarily tested on AbTYR as a rapid and inexpensive screening procedure before being tested on hTYR.


Introduction
Tyrosinase (TYR, EC 1.14.18.1) is a copper-containing oxidase able to convert phenolic substrates into polyphenolic compounds. Different isoforms of TYR are expressed in various organisms, each with specific physiological roles. The catalytic activity of human TYR (hTYR) regulates melanin formation in melanocytes and is thus responsible for skin pigmentation. In particular, hTYR produces o-quinone compounds by hydroxylation and oxidation of the amino acid tyrosine; the o-quinone intermediates are precursors of melanin in a complex biosynthetic pathway for which hTYR is the rate-limiting enzyme [1]. In addition, both TYR-related proteins 1 (TYRP-1) and 2 (TYRP-2) contribute in melanogenesis through the final steps of eumelanin production [1,2].
Abnormal melanin production and corresponding hTYR overactivity have been shown to be responsible for various skin diseases such as post-inflammatory hyperpigmentation, age spots, and malignant melanoma [3]. Consequently, there is ongoing research aimed at identifying hTYR inhibitors for the development of therapeutics for these skin disorders [3][4][5][6]. Over the past decades, a large class of natural and synthetic compounds has been identified by using the enzyme from Agaricus bisporus (AbTYR) as a cheap surrogate of hTYR to test inhibitory effects. Among them, β-arbutin (1) [7] and kojic acid (2) [7] (Fig. 1) are approved for cosmetic applications, although they have recently been shown to have serious adverse effects [7]. The use of AbTYR to identify new therapeutics for skin diseases has recently faced sharp criticism. Roulier, in particular, emphasized the need for screening procedures that consider the structural diversity between hTYR and the TYR isoforms found in other organisms [5]. In particular, the AbTYR is a cytosolic oligomeric enzyme; whereas hTYR is a highly glycosylated monomeric protein anchored in the melanosome membrane. Even if the dicopper-coordinating hystidines are well conserved, the second sphere of coordination displays a certain variation so that hTYR and AbTYR share very low overall identity. Near the hTYR active site, several residues (H304, K306, R308, T343, T352, I368, S375 or S380) are crucial in establishing interactions with known potent hTYR inhibitors; interestingly, the above-mentioned residues are not present in AbTYR; especially, residue S380 seems to be involved in substrate activation as well as in the inhibitory activity of potent hTYRIs. This structural diversity might explain the reason behind the occurrence of failure for AbTYR inhibitors in exerting high inhibitory effects towards hTYR [5] thus limiting their potential therapeutic applications in hyperpigmentation therapy. Currently, different methods have been developed to test potential inhibitors against hTYR, which is recombinantly expressed in bacterial, insect and human cells [5].
In an extensive screening campaign with a library of 50,000 compounds, N- [4-(2,4-dihydroxyphenyl)-1,3-thiazol-2-yl]-2-methylpropanamide (3) (Thiamidol, Fig. 1) proved to be one of the most promising inhibitors against hTYR, which was expressed in human endothelial kidney cells and displayed an inhibition constant (K i ) value of 0.25 μM and a half maximal inhibitory concentration (IC 50 ) of 1.1 μM [8]. Thiamidol has been found to reduce melanin production in both in vitro and in vivo assays [9]. It is an active ingredient in cosmetic products for the topical treatment of facial melasma [10]. Furthermore, potent inhibitory effects on hTYR have been discovered for aurone derivatives and a class of resorcinol-based compounds (i.e. compounds 4 and 5, Fig. 1). Compounds 4 and 5 possessed similar K i values (0.35 μM and 0.25 μM, respectively) and were effective in hTYR-containing human melanoma MNT-1 cell lysates, showing IC 50 values of 15 μM and 1.6 μM, respectively [11].
The rational drug discovery of new potent hTYR inhibitors is hampered because there is not exact information about threedimensional structure of hTYR. However, crystal adducts of inhibitors bound to TYR from bacteria and fungi (e.g. Streptomyces castaneoglobisporus, Bacillus megaterium and Agaricus bisporus) have provided insight into the kinetic and structural properties of the catalytic cavity of TYR isozymes [12][13][14]. To decipher the three-dimensional structure of hTYR, several homology models have been constructed based on its similarity to other proteins, such as TYRP-1 and TYRP-2 [1,15]. These studies have provided suggestions for the structural arrangement of hTYR and have described potential binding sites for highly efficient inhibitors, including Thiamidol, aurone-derivatives and resorcinol-based compounds [11]. Our recent endeavors involved the development of TYR inhibitors revealed that a new class of AbTYR inhibitors can be attributed to 4-hydroxyphenylpiperazine derivatives [16,17]. Among them, compound 2-(4-benzyl-1-piperidyl)-1-[4-(4-hydroxyphenyl)piperazin-1-yl] ethanone (6) (Fig. 1) has proven to be effective AbTYR inhibitor at low micromolar concentration (IC 50 value of 3.80 μM) [17].
Considering that developing a successful hTYR inhibitor is challenging due to the difficulty of transferring the AbTYR inhibitory property to the human isozyme, herein we conducted biochemical screening and developed computational models to understand and compare the binding modes of selected compounds on hTYR and AbTYR; our strategy consisted in building an in silico model capable to predict the chemical features for dual inhibitors interacting with these two distinct enzymes with the same potency.

Results and discussion
Initially, we focused on our active AbTYR inhibitor 6 and decided to test its ability to inhibit hTYR using a procedure described by Winder and Harris [18]. The enzymatic activity was assessed using a spectrophotometric method, which employs 3,4-dihydroxy-L-phenylalanin (L-DOPA) as substrate. Upon oxidation of L-DOPA by hTYR, dopaquinone is formed, which subsequently undergoes autoxidation to dopachrome with an absorbance peak at 475 nm. The formation of dopachrome can be suppressed by the addition of inhibitory compounds, enabling the evaluation of the inhibitory efficacy as IC 50 in a dose-response curve. The IC 50 values were estimated by performing a non-linear regression using GraphPad Prism 9. Our screening revealed that compound 6 was a weak inhibitor of hTYR with an IC 50 value of 1.1 mM. It resulted less active in comparison with its ability to inhibit AbTYR (IC 50 of 3.80 μM), as previously reported in our paper [17]. On the other hand, the well-known inhibitor Thiamidol (3) exhibited good activity on hTYR with an IC 50 of 3.8 μM (Fig. 2), compared to the inhibition of AbTYR with an IC 50 of 108 μM, as reported in literature [8].
This opposite behavior of inhibitors 3 and 6 towards the two isozymes prompted us to investigate the key factors controlling the recognition process within the two catalytic sites.
Looking for the similarity between hTYR and AbTYR we developed a computational process consisting in different steps: a) to build a homology model for hTYR as there is currently a lack of three-dimensional structural data for the human isozyme on RCSB Protein Data Bank (PDB); b) to use and compare this homology model and the already available X-ray model of AbTYR for docking Thiamidol and compound 6, which represents our prototype of AbTYR inhibitor; c) to analyze the  different and similar binding mode of the inhibitors on hTYR and AbTYR. In order to fine tune a reliable homology model of hTYR for our docking studies, we started by creating a pool of six different hTYR homology models using the following software: Alphafold (V.2) [20], SwissModel [21], MODELLER [22] and TopModel [23] (see Supplementary Material S1).
Selection of the best model was performed by (i) initial visual evaluation and refinement using the Molecular Operating Environment (MOE, Chemical Computing Group, Montreal, Canada), (ii) geometric analysis using PROCHECK [24] and WHATCHECK [25], and (iii) molecular dynamics simulation of docking complexes and selection of plausible frames from the resulting trajectories. A detailed description of the workflow is provided in Supplementary Material S1. As result from this rigorous selection procedure, we selected the model originally generated by MODELLER (named Model 4) for subsequent docking studies.
Firstly, Thiamidol was rigidly docked on Model 4 of hTYR and on AbTYR (PDB code: 2Y9X [12]) and the obtained complexes have been minimized; its binding modes are shown in Fig. 3. Thiamidol contains a resorcinol moiety that is known to play a crucial role in providing interactions within the catalytic cavity of hTYR (Fig. 3A). In fact, it is able to coordinate one of the two copper ions (CuA) through a hydroxyl group located in the para position of the phenyl ring, which forms a hydrogen bond interaction with S380, a crucial residue for hTYR catalytic activity [5]. The ortho-hydroxyl group is also involved in a hydrogen bond interaction with the side chain of N364 via the carbonyl oxygen. The resorcinol fragment also interacts with the hydrophobic residue V377 and ππ stacks H367. Both the thiazole and isopropylamide moieties are stabilized by hydrophobic contacts within the hydrophobic active site entrance lined by I368, F347 and Ala357. Lastly, the nitrogen of the amide linking group establishes a hydrogen bond interaction with the hydroxyl group of S360.
Analysis of the docking of Thiamidol to AbTYR (Fig. 3B) revealed the coordination of the CuA ion through the para-hydroxyphenyl ring via interaction with the copper-chelating histidines. In contrast, the orthohydroxyl group of the resorcinol moiety formed a hydrogen-bond interaction with the carbonyl oxygen of the backbone of residue M280. It is evident that Thiamidol is engaged in several hydrophobic interactions within the catalytic cavity of AbTYR. Specifically, the resorcinol ring interacted with A286 and V283, while the thiazole ring established contact with F264. Interestingly, the isopropylamide group was found to lack favorable interactions with the hydrophobic pocket of this region of the cavity, likely due to a possible clash with F264.
Upon comparing the poses of Thiamidol in the two isozymes ( Fig. 3C), it was possible to re-highlight that the para-hydroxyl group of the resorcinol moiety formed a hydrogen bond with S380 in hTYR, that is a crucial residue for hTYR, thus confirming its potential role in inhibiting this enzyme; in contrast, the ortho-hydroxyl group of Thiamidol interacted with N364 in hTYR, while it bound to M280 in the AbTYR, and we found no evidence suggesting the significance of this interaction in inhibiting these two proteins. With regards to the amide nitrogen, it interacted with two different residues, G281in AbTYR and S360 in hTYR, located in completely different regions of the binding site, thus, we cannot make any assumptions about the hydrogen bonds formed by this group in both enzymes. However, this interaction may stabilize the additional lipophilic interactions formed by the isopropyl group in hTYR, which interestingly are lost in AbTYR. The binding stability of the Thiamidol inside hTYR and AbTYR binding site has been evaluated performing MD simulation on the docking complexes (see Supplementary Material).
Overall, our computational findings are in good agreement with experimental data that show Thiamidol to be a more potent single digit micromolar inhibitor of hTYR ( Figure) [8,19] compared to AbTYR (IC 50 = 108 μM) [8]. These results support the idea that the pattern of favorable interactions identified in our computational analysis play a critical role in stabilizing hTYR inhibitors compared to those targeting the AbTYR isoform.
Then we performed the docking analysis of compound 6, that displays a complete flip in inhibitory effects when compared to Thiamidol. In detail, inhibitor 6 showed high effects against AbTYR (IC 50 value of 3.80 μM) [17], but it failed to produce high inhibition of hTYR (IC 50 value of 1.1 mM) (see Fig. 2).
To perform the in silico study we employed our hTYR Model 4. Fig. 4A illustrates a plausible docking pose of the weaker inhibitor 6 and its interactions with hTYR. The para-hydroxyl substituent formed hydrogen bond interactions with S380, while the aromatic ring establishes hydrophobic contacts with V377 and forms ππ stack with H367.
Overall, the orientation of the hydroxyphenolic fragment of inhibitor 6 was similar to that of the resorcinol moiety of Thiamidol. However, unlike Thiamidol, no additional interactions stabilized the connecting scaffold or the benzylpiperidine tail of the inhibitor 6. Nonetheless, lipophilic contacts involving the piperidine ring and I368 may exist, as suggested by LigPlus [26] and Ligandscout [27,28] receptor-binding surfaces, even though they are not stabilized by any interaction of the piperazine scaffold or of the benzyl tail. At first glance, the lack of relevant contacts with our modeled hTYR may explain the lower ability of inhibitor 6 against this protein compared to AbTYR.
The docking analysis of inhibitor 6 with AbTYR ( Fig. 4B) revealed that the inhibitor can coordinate CuA as well as interacts with residues H85 and H296 through hydrogen bonds. The oxygen atom of the amide group also formed a hydrogen bond with the backbone nitrogen atom of V283. Additionally, several hydrophobic contacts were observed, including the piperidine ring with F264 and V248, the aromatic ring of benzylpiperidine tail with M259 and T261 and a π-π T-shaped interaction with F264. By analyzing the two distinct poses of inhibitor 6 on hTYR and AbTYR, it was found that in the case of AbTYR, the lipophilic tail provides additional interactions at the entrance of the binding site. However, the same hydrophobic tail induced the loss of stabilizing interactions during the recognition within hTYR. Moreover, the inhibitor 6 projected its benzyl tail toward the highly hydrophilic and solventexposed the region of hTYR. Information about the binding stability of this compound inside hTYR and AbTYR binding site were obtained through MD simulation as reported in Supplementary material. By comparing the Thiamidol binding-mode in hTYR (purple) with the binding mode of compound 6 in AbTYR, it might be speculated that the lipophilic tails of each inhibitor are targeted towards the same region of the catalytic pocket of hTYR (for Thiamidol) and AbTYR (for inhibitor 6), as shown in Fig. 5. This region in AbTYR is characterized by an α-helix with a hydrophobic surface, whereas in hTYR, it is an extended loop with an unpredictable conformation. This suggests that the lipophilic interactions at the entrance of hTYR may be restricted to V377, I368 and F347.
The hypothesis regarding the different inhibitory activities of the two studied compounds (Thiamidol and 6) against the two distinct isoforms AbTYR and hTYR was confirmed by the binding energy values obtained from the molecular mechanics generalized Born surface area (MM-GBSA) calculation (Table 1). The MM-GBSA calculation of the difference in binding energy of the enzyme-inhibitor complex revealed that Thiamidol exhibited a stronger binding affinity for hTYR, while compound 6 exhibited a stronger binding affinity for AbTYR.
Based on the above-mentioned suggestions that hydrophobic interactions on the rim of hTYR cavity might play a crucial role, we were interested in exploring the impact of the length and flexibility of the spacer connecting the aromatic tail to the piperazine central core of inhibitor 6. We thought as a potential positive modification the shortening of the length of this chemical spacer to afford an improvement of hTYR activity and we took in to account another previous series of AbTYR inhibitors synthesized in our laboratories, namely (4-(4hydroxyphenyl)piperazin-1-yl)arylmethanone derivatives [16]. Therefore, we focused on a prototype of this series, the [4-(4-hydroxyphenyl) piperazin-1-yl](2-methoxyphenyl)methanone (7, so-called MehT-3), which exhibited an IC 50 of 3.52 μM for diphenolase activity of AbTYR by competitively antagonizing the protein [16].
Before to test the in vitro inhibitor effects of compound 7 (MehT-3) against hTYR we tried to predict its binding mode through our docking protocol. So, MehT-3 was docked on our predicted structure of hTYR (Model 4). Specifically, MehT-3 established contacts with V377 through the 4-hydroxyphenyl moiety, while both the piperazine core and phenyl tail created interactions with I368 and F347 (Fig. 6A). Interestingly, the MehT-3 was able to form almost the same lipophilic interactions of Thiamidol (Fig. 6B) suggesting its plausible efficacy as hTYR inhibitor. We also analyzed the binding pose of MehT-3 on AbTYR to explore its ability to act as dual binding inhibitor.
The Fig. 6C displays the binding-mode of compound MehT-3, which exhibited almost the same hydrophobic interactions as Thiamidol and inhibitor 6 (Fig. 6D) through the hydroxyphenyl ring. However, the benzoyl tail of MehT-3 assumed a hybrid orientation into AbTYR unlike the tail of the two other inhibitors (Thiamidol and parent compound 6). The isopropyl group of Thiamidol was not involved in any interaction with AbTYR, whereas in the case of MehT-3 the 2-methoxybenzoyl group interacts with F264 and P277. Due to the shorter length of the tail of inhibitor MehT-3, some lipophilic interactions are lost (e.g., with M257 and V248). The methoxy group of MehT-3 forms a hydrogen bond with the V283 backbone nitrogen, like compound 6, suggesting that a hydrogen-bond acceptor close to the lipophilic tail might form additional interactions with AbTYR. We may speculate that this latter interaction is not formed by any of the three inhibitors in hTYR due to the different backbone conformation of the homologous V377 and their flanking residues.
Overall, inhibitor MehT-3 was predicted to be able in maintaining the crucial contact with S380 and create additional hydrophobic interactions to respect parent compound 6. To gain information about the binding stability of this compound inside hTYR and AbTYR binding site we carried out MD simulation studies (see Supplementary material). The hypothesis that the inhibitor MehT-3 was a nice double ligand was corroborated by our MM-GBSA calculations showing negative binding energy values in both enzymes (− 36.15 kcal/mol and − 10.42 kcal/mol for hTYR and AbTYR, respectively).
Encouraged by this promising binding pose of MehT-3 in the hTYR cavity (as shown in Fig. 6A), we conducted an inhibitory assay on hTYR. The biochemical screening revealed that MehT-3 possessed a good activity (IC 50 of 5.4 μM) (Fig. 7), thus demonstrating a significant improvement of inhibitory effect with respect to parent compound 6 (IC 50 of 1.1 mM).
To examine the mode of hTYR inhibition of the best compound MehT-3 (7), we carried out kinetic studies using L-DOPA as substrate. The compound exhibited significant inhibitory activity in the low    Fig. 8. Furthermore, the competitive inhibition pattern exhibited by MehT-3 was found to be analogous to that observed for the reference compound Thiamidol, as illustrated in the right panel of Fig. 8.
Notably, this result aligned well with our computational studies, which indicated that the presence of a 2-methoxyphenyl was advantageous for the occurrence of additional interactions with the hydrophobic area at the entrance of the cavity of hTYR. A preliminarily SAR consideration reveals that the differences in the structure of the two 4hydroxyphenylpiperazine-derivatives 6 and MehT-3 (7) seem to be not relevant for the inhibition of AbTYR; actually, that compounds 6 and 7 possessed very similar inhibitory effects towards AbTYR to respect their divergent efficacy in hTYR inhibition. As for compound 7, the shorter linking moiety seems to be preferred to afford dual AbTYR/hTYR inhibitor active in the low micromolar range (3.52 μM and 5.4 μM, respectively) through a competitive inhibition of enzyme activity. The molecular modeling analyses suggested that the structural properties of MehT-3 facilitated hydrophobic interactions with crucial residues at the rim of enzyme cavity of AbTYR (M280, V283, F264) and hTYR (F347, H367, I368 and V377); the above-mentioned network of interactions was proven by docking simulations for various inhibitors from synthetic and natural sources that have been reported in literature [29][30][31][32][33][34]. Among them, 2,4-diphenyl-2,3-dihydro-1,5-benzothiazepines and 4-thioflavonols resulted AbTYR inhibitors at low micromolar concentration; for these inhibitors the π-π interaction with F264 was predicted to be relevant for improve the interaction with AbTYR; moreover, the additional hydrophobic contacts were effective in stabilizing the inhibitor complex, especially the interactions with M280 [30,33,34] and V283 [32]. In the case of 2-phenylchromone derivatives, hydrophobic π-sigma interactions with the M283 often occurred [32]. In addition, docking simulations revealed H-bonding interactions engaged by AbTYR inhibitors and the pair of crucial residues M280 and V283 [31,34].
Notably, MehT-3 maintained the canonical polar interactions in the depth cavity near to CuA and S380 polar as observed for the potent hTYR inhibitors 4 and 5 on the computational studies on modeled hTYR reported in previous papers [11].

Conclusion
This retrospective study revealed that certain residues (either conserved or homologous) should be targeted to develop novel dual inhibitors active towards hTYR and AbTYR. Specifically, these residues include H367 (263 in AbTYR), I368 (F264 in AbTYR), V377 (283 in AbTYR) and S380 (A286 in AbTYR). Based on the binding modes of the studied inhibitors 6,7 and thiamidol, some considerations can be made to guide the development of newer hits. One key difference is the nature of the copper chelator moiety of active hTYR inhibitor, which should preferably possess a hydrogen bond donor or acceptor function to establish polar contact with crucial residue S380 located near to dicopper center.Additionally, the length and nature of the lipophilic tail should also be considered, as well as the presence of small hydrophobic groups seem to be sufficient to strongly inhibit hTYR, providing a better shape complementarity with the entrance of active site lined by I368 and F347, while in the case of AbTYR the occurrence of aromatic ring is suggested for a stronger binding to F264 located at greater distance from catalytic cavity to respect F347.

Inhibitory assay of hTYR and kinetic studies
The intramelanosomal domain of hTYR was expressed and purified as previously reported in literature [35]. The inhibitory efficacy of compounds on human Tyrosinase (hTYR) was determined using a Fig. 7. Dose-response curve showing the inhibitory efficacy of compound MehT-3 on hTYR in a L-DOPA oxidase assay [18,19]. Thiamidol was used as positive control, with an observed IC 50 of 3.8 μM (literature 1.1 μM) [8]. The IC 50 of MehT-3 was determined to be at 5.4 μM.  L-DOPA oxidase assay as described in literature [18,19]. Thiamidol was purchased from BLDpharm (Shanghai, China). Compounds 6 and 7 (MehT-3) were prepared following a previously reported procedure as detailed in Supplementary Material. Thiamidol was used as a positive control [8]. A) The reaction mix (PBS, 2% DMSO, 1 mM L-DOPA, pH 7.4) was prepared in a 384-well plate by mixing 21 μl of 1.43 mM L-DOPA (in PBS, 2% DMSO, pH 7.4) and 4 μl of inhibitor (in PBS, 2% DMSO, pH 7.4).
Additionally, the reaction mix was prepared without inhibitor substituting by 4 μl PBS, 2% DMSO, pH7.4. The experiment was performed in triplicates. With a multi-channel pipette, 5 μl of 1.2 μM recombinant hTYR (in PBS, 2% DMSO, pH 7.4) were added to each reaction mix and the formation of dopachrome was immediately measured by spectrophotometrically monitoring the absorption at 475 nm every 10s for 10 min using a SpectraMax Paradigm microplate reader. The slopes for each inhibitor concentration were normalized against the slope of the reaction mix without inhibitor, and the obtained enzymatic activity (%) was plotted against the respective inhibitor concentration. IC 50 values were estimated by performing a non-linear regression analysis of the dose-response curve using GraphPad Prism 9. B) The in vitro kinetic experiment was performed in a buffer solution consisting of PBS, 2% DMSO (v/v) at pH 7.4. Six different concentrations of L-DOPA, covering a range from 0.06 to 1.9 mM, and three concentrations of inhibitors (0.52 μM, 1.04 μM, and 2.08 μM) were tested. For each condition, the experiment was performed in duplicate using 200 nM tyrosinase. The change in absorption was measured at 475 nm every 10 s over a period of 10 min using a Spectramax Paradigm spectrophotometer. The estimation of the inhibition constant K i was carried out through non-linear regression analysis utilizing Graphpad Prism 9 software, employing the best-fit value derived from the competitive inhibition model.

Homology modelling: protein preparation
The six initial models were generated giving as input the FASTA file P14679 downloaded from UniProt database. The output results of the four programs were then aligned and superimposed using the sequence and structural alignment method on MOE (Version 2020.09.1) using the Model1 from SwissModel (PDB template: 5M8L) [2] as reference. The models containing Zn ions in the active site taken from their respective pdb templates, were modified converting the Zinc ions into Cu ions, while in the two models obtained from Alphafold, not containing any metal ion in active site, the coordinates were sampled from the model used as reference for the alignment (Model 1). Structure preparation was performed using MOE by adding missing side chain atoms, applying protonation states at physiological pH using Protonate3D [36], calculating charges using the OPLS-AA force field, capping C termini with NME and N termini by transferring them from complete models.
Protonation of the copper-chelating histidines was further adjusted setting them as HID. Clashes and psi-phi outliers reported by the structure preparation wizard concerning residues out of the active-site or distant from the di-copper center were tried to be fixed through the protein builder minimization protocol (minimize residue or minimize side-chain when regarding clashes) using OPLS-AA forcefields. The ones regarding residues close to the copper-chelating histidines were accounted for seeing whether they were confirmed also by WHAT-CHECK [25] and PROCHECK [24] and finally by visual inspection of the van der Waals radii. D-amino-acids as well as cis-amide bonds were accounted to see whether they were confirmed by WHATCHECK (D-amino-acids) and PROCHECK(cis-amide).

Ligand preparation
All the ligands were sketched using ChemDraw and saved as mol files. They were then processed on MOE Database Viewer (version 2020.09.1) using the following protocol: (i) Washing using Dominant protonation protocol setting a pH of 7.4, generating 3-D coordinates (ii) Processing with QuickPrep using default settings (iii) Calculation of charges using the OPLS-AA forcefield with default values, (iv) minimization using the OPLS-AA forcefield with default values and (v) saved in mol2 format.

Docking protocol for the co-crystalized ligands
The models retained from the first prefiltering (Model 4, Model3 and Model6) were submitted to rigid docking procedure on the program GOLD (2021) [37] using the ligands L-Tyrosine(substrate), L-DOPA (substrate) Hydroquinone(inhibitor) and Kojic-acid(inhibitor) earlier co-crystalized in BmTYR and TYRP1 [38,39]. The complex of L-Tyrosine-BmTYR (PDB: 4P6R, chain A) [38] was superimposed to each model and the coordinates of L-Tyrosine were sampled for defining the center of the binding-site with a radius of 10 Å. The ligands were submitted to 100 genetic algorithm runs and solutions were clustered using a threshold of 0.75 Å. Scaffold constraints, considering the coordinates of the phenolic ring of L-Tyrosine pose in the X-Ray structure (PDB: 4P6R, chain A), were applied for the docking of L-Tyrosine and L-DOPA. Flip planar N and flip NRR and NHR function was set to true. All other parameters were left as default.
The clusters (top ranked solution per cluster) obtained from docking were merged in each respective model analyzed and the X-ray structures of the ligands used for validating the models were superimposed for choosing the poses from docking to submit to MD simulation. The poses with the highest superimposition with the X-ray conformers were selected for each of the three models.

MD simulation of the co-crystallized ligands on hTYR
Three replicas of 50 ns were performed for each docking complex selected using Desmond software (Version 2021, D.E. Shaw Research, 120 W. 45th St, NY, USA). Each docking complex was prepared prior the simulation as follows: (i) Removed N-terminus (signal peptide, residues: 1-18) and C-terminus (cytosolic domain, residues: 498-529). For Model 4, which lacked of the transmembrane domain (residues: 476-497), MOE (version 2020.09.01) has been used to align Model 3 on Model 4 and for sampling atom coordinates of the target region from the former, starting from the flanking region with highest superimposition in the two models, and (ii) Preprocessed the docking complexes on Maestro (Schrödinger, NY, USA) using the Structure Preparation Wizard, recapping both N and C-termini and leaving the other settings as default.
Each system was created through the Desmond System builder tool using the following settings. The simulation was performed at 310 K, using NPγT ensemble and the Noose-Hoover chain thermostat. The atom coordinates for trajectory analysis were recorded each 50ps (1000 frames).

Docking protocol for thiamidol, compounds 6 and 7 on hTYR
The same settings employed for the docking of the ligands used to validate the homology models of hTYR were applied for the docking Thiamidol and inhibitor 6. The top ranked pose of each ligand was selected for MM-GBSA binding energy calculation on Prime (Schrodinger) prior binding-mode comparison.

Docking protocol for thiamidol, compounds 6 and 7 on AbTYR
The three compounds object of the study were rigidly docked in AbTYR (PDB code: 2Y9X). The protein structure was prepared using the Protein Preparation Wizard using the default setting. Again, 100 GA runs per compound were performed on GOLD software (2021), using the default settings and clustering docking solution using a threshold of 0.75 Å. A 10 Å radius from the following coordinates was set to define the binding-site: 10.021; − 28.823; − 43.596. A corrected version of ASP scoring function which has an additional metal coordination term, was employed for ranking the docking poses, since a better superposition during re-docking of the native co-crystallized ligand (Tropolone) was observed. The top ranked pose of each ligand was selected for MM-GBSA binding energy calculation on Prime (Schrodinger) prior binding-mode comparison.

MM-GBSA calculation
The MM-GBSA calculation was performed using the Prime MM-GBSA [40,41] tool of the Schrödinger Suite [release 2020-3] [42]. First, the Refine Protein-Ligand Complex function of Prime was used to optimize the complex from the docking calculation, considering the refinement of the residues within 5 Å of the ligands. VSGB and OPL3e were employed as the solvation model and force field, respectively. In addition, constraints were applied to the two copper ions and, in the case of hTYR, the membrane was set up at the transmembrane residues 476-497 level. Then, for each complex the binding energy was predicted by means the MM-GBSA method. Again, VSGB was used as the solvation model, OPLS3e was set as the force field, and the implicit membrane model was employed only for the hTYR complexes. The other parameters were kept as default values.

Interaction analysis
The minimized complexes of Thiamidol, 6 and 7 with hTYR and AbTYR obtained from MM-GBSA calculation were analyzed on three different platforms, whose results were merged in order to develop our SAR hypothesis.
-Discovery Studio Ligand-receptor interaction tool to account lipophilic (described as Alkyl, phi -Alkyl, phi -sigma or phisulfur), phi-phi (stacked or T-shaped) and ionic interactions, as well as solvent exposure; -LigPlus to account for metal coordination and van der Waals interactions -Ligandscout (v.4.4.8): pharmacophores were employed to account hydrogen-bonds, lipophilic, phi-phi and ionic interactions involving specific moieties/features/chemical-physics properties of the ligand and their environmental patterns.
Receptor binding surfaces (data not shown) of each complex were used to rationalize vdW interactions reported by LigPlus according to the nature of the environment/binding-site surfaces (hydrophobic or hydrophilic); The results obtained by the two software were merged to trace the SAR of the inhibitors for two the enzymes examined.
All the settings were left as default.