Computational approaches to predict binding interactions between mammalian tyrosinases and (S)-(+)-decursin and its analogues as potent inhibitors

Abstract

Tyrosinase is a metalloprotein involved in the production of melanin, but overexpression of this pigment is problematic in dermatological disorders. (S)-(+)-Decursin and its analogues have been shown to exhibit inhibitory activities against melanin production in B16 murine melanoma cells and specifically inhibit tyrosinase. In order to elucidate their binding modes with mouse tyrosinase, we used computational methods to predict the structures of the complexes on the atomic level. In addition, the structure of human tyrosinase, which has high amino acid sequence identity with mouse tyrosinase, was studied to identify the possible binding mode and activity of (S)-(+)-decursin and its analogues for human treatment. The results showed that (S)-(+)-decursin and its analogues contain ester derivatives that interact with hydrophobic and aromatic residues, particularly Ile256 and Val265 in the binding site. Although human tyrosinase displayed similar interactions with the compounds, there was a key structural difference in residue 222, which caused some compounds to bind with lower binding affinities. Our results provide a foundation for the identification of new tyrosinase inhibitors for treating dermatological disorders.

---

1. Introduction

The melanosome is a melanin cell organelle responsible for synthesizing and storing melanin.¹ Melanin protects human skin against harmful rays or chemicals; however, hyperpigmentation, a characteristic of aging skin, can be problematic.² Tyrosinase plays a major role in the production of melanin. It is a catechol oxidase belonging to the type 3 copper protein group and is expressed in all living organisms. Tyrosinase converts monophenol into o-diphenols, from which o-quinones are produced by oxidization. The substrate of tyrosinase is L-tyrosine, which is converted to L-dopaquinone. Melanin is the final product of this metabolic process and plays a role in determining the color of animal hair and skin. One peroxidase molecule and two copper ions are present at the tyrosinase active site and each copper ion is coordinated with three histidine residues.^3–6

Tyrosinase inhibitors have been studied experimentally and theoretically for several decades to develop cosmetic, agricultural, and therapeutic products.^7–16 Several natural or synthesized tyrosinase inhibitors have been reported such as arbutin, kojic acid, L-ascorbic acid, and phenylthiourea. However, most have shown limitations because of insufficient inhibitory activity, instability, and cytotoxicity, necessitating the development of a new tyrosinase inhibitor that is more potent and highly specific.

(S)-(+)-Decursin has been reported as a new anti-cancer drug candidate and was shown to display various biological activities following extraction from Angelica gigas, Nakai. In our previous study, (S)-(+)-decursin and its derivatives were found to be effective inhibitors of melanin production in B16 melanoma cells for the treatment of dermatological disorders.⁸ Furthermore, through 3D-QSAR analysis, we determined which parts of these compounds are important for their effectiveness.¹⁷

In this study, we analyzed the detailed interactions of (S)-(+)-decursin and its derivatives with mouse tyrosinase, for which crystal structures are unavailable, using homology modeling, molecular docking, quantum mechanics/molecular mechanics (QM/MM), molecular dynamics simulations, and molecular mechanics with generalized Born surface-area solvation (MM-GBSA) calculations. Human tyrosinase, which has high amino acid sequence identity, was also studied to identify the possible binding mode and activity of (S)-(+)-decursin and its analogues. This information can ultimately be used to develop inhibitors of tyrosinase for human use.

To study binding interactions between both mammalian tyrosinase structures with (S)-(+)-decursin and its analogues, we first determined the parts of the inhibitors and tyrosinase homology model that interact and the relationships of these findings with our previous QSAR results through docking simulations. Arbutin was also studied and compared with (S)-(+)-decursin and its analogues as a reference ligand. Force field-based MM modeling of the tyrosinase–ligand complex to obtain energy-minimized structures revealed an unstable structure because of improperly parametrized copper ion–histidine residues. Thus, QM-level calculations were included to model the part of the complex containing metal ions. QM/MM calculations were employed so that costly QM calculations could be minimized. After optimizing the structures by QM/MM, additional molecular dynamics simulations were performed to identify protein–ligand interactions. Next, the binding free energies were calculated using MM-GBSA calculations based on molecular dynamics (MD) trajectories. Since the crystal structures of mouse tyrosinase (mTyr) and human tyrosinase (hTyr) are not available, it is necessary to predict protein–ligand interactions using computational methods to design new tyrosinase inhibitors.

2. Materials and methods

2.1. Ligand sets

Arbutin, (S)-(+)-decursin, and (S)-(+)-decursin analogues as potent inhibitors of melanin formation, determined in our previous study, were used in this study as shown in Fig. 1.⁸ The biological activities of these compounds were measured at an inhibition rate at 100 μM and we converted these values to IC₅₀ values assuming linear relationships around the data points. This is reasonable given that most of the measured inhibition rates were greater than 50%. These converted values, along with their negative logarithms, were used to determine the correlations between MM-GBSA binding free energy and experimental IC₅₀ values for mTyr.

Fig. 1 Structure and in vitro activity of tyrosinase inhibitors. The atom numbering is shown in compound 1.

2.2. Sequence alignment and homology modeling

The amino acid sequences of tyrosinases isolated from Mus musculus (mouse, accession number: NP_035791.1), Homo sapiens (human, accession number: NP_000363.1), Bacillus megaterium (bacteria, accession number: ACC86108.1), Streptomyces castaneoglobisporus (bacteria, accession number: AAP33665.1), Agaricus bisporus (fungus, accession number: ADE67052.1), and Asperugillus oryzae (fungus, accession number: 3W6W_A) were retrieved from the NCBI protein sequence database.^18,19 Global and local pairwise sequence alignment of the six amino acid sequences was performed using the EMBOSS package to calculate sequence identity and similarity among the species.²⁰ To identify a suitable template, the amino acid sequences of mouse and human tyrosinases were searched against the Protein Data Bank (PDB) database using the PSI-BLAST algorithm.^21,22 The crystal structures of bacterial and fungal tyrosinase (PDB ID: 1WX2 and 3NM8) were used as templates to build a mouse and human tyrosinase structure.^5,23 Homology modeling of the mouse and human tyrosinase was performed using the Prime homology modeling program of Schrödinger.²⁴ The crystallographic positions of the backbone atoms and conserved side chains were mapped from the template, while the side chain coordinates of all non-identical residues were predicted. Copper ions were added to the model in the model building step. Next, the stereochemical quality of the homology models was evaluated using RAMPAGE.²⁵ This program provides an all-atom contact analysis that identifies any steric problems within the molecule and uses high accuracy Ramanchandran and rotamer distributions to check the main chain and side chains for conformational outliers.

2.3. Molecular docking

We used Glide by Schrödinger for docking.^26,27 Glide uses grids for fast scoring, and we used its grid-generation module to generate grids for the mTyr and hTyr structures produced through homology modeling. We set the scaling factor of the van der Waals radii to 0.8 and the partial charge cutoff to 0.15. The active site of tyrosinase, which contains copper atoms, was included in the grid generation. Six compounds were selected from arbutin, (S)-(+)-decursin, and its derivatives as ligands.^8,17 These selected compounds were drawn by 2D sketcher and optimized with MacroModel.²⁸ Next, all possible ionization states and stereoisomer structures of the ligands were generated using the ionizer option in LigPrep.²⁹ We used the SP mode of Glide to produce 10 poses per ligand while performing the docking of these six compounds to tyrosinase. The ligand interaction diagram module of Glide was used to analyze ligand–protein interactions.

2.4. QM/MM minimization

Geometry optimizations of tyrosinase active sites were performed using QSite, the Schrödinger's QM/MM module based on the frozen-orbital method.^30,31 The inclusion of copper ions necessitated the use of QM-level calculations to optimize the active site since there are no proper force fields with which to describe copper ions. For our QM/MM calculations, we included two copper ions, and the six histidine residues to which the copper ions in the QM region are coordinated. Density functional theory with a B3LYP/lacv3p** hybrid functional was used for the QM calculations, whereas the OPLS 2005 force field was employed for the MM calculations.^32–36 The QM and MM calculations can be set separately for QM/MM minimization in QSite. We set the maximum number of iterations to 200 for QM and 7000 for MM using the truncated Newton algorithm method. After QM/MM minimization calculations, we calculated the partial charges of the atoms in the QM regions using electrostatic potential fitting. These charge values were fed into the structure files for two copper ions and six histidine residues.

2.5. Molecular dynamics simulations

Molecular dynamics simulations of mammalian tyrosinase and inhibitor complexes obtained from QM/MM minimization were performed using the DESMOND module available from D. E. Shaw Research Group using the OPLS-2005 force field.^37,38 The complex structures were immersed in an orthorhombic box containing TIP3P water molecules and neutralized by adding salt counter ions.³⁹ The distance between the box wall and structures were set to greater than 10 Å to avoid direct interactions with its own periodic image. The structures were relaxed before simulation using a series of energy minimizations and short MD simulations. There were six steps in the relaxation process; the first two steps were minimization steps. The first was solute-restrained and the second had no restraints. Steps three through six were short MD simulations of 12, 12, 24, and 24 ps using the NPT ensemble at 10 K, 10 K, 300 K, and 300 K, respectively. The relaxed system was simulated for 30 ns with a time step of 2 fs, the NPT ensemble using a Nose–Hoover thermostat at 300 K, and Martyna–Tobias–Klein barostat at 1.01325 bar pressure.^40–42 Recent studies demonstrated the importance of electrostatic interactions between transition metal ions such as copper ions and biomolecules.^43–45 Wise and Coskuner developed new force field parameters for metalloproteins including copper ions, which is in excellent agreement with available experimental data.⁴⁵ Hence, a positional restraint of 30 kcal mol⁻¹ was applied to the copper ions and six coordinated nitrogen atoms of the histidine residues. The partial charges of two copper ions and six coordinated histidine residues calculated by QM/MM calculations were used. Atomic coordinate data and system energies were recorded every 4.8 ps. Root-mean-square deviation (RMSD) values were calculated using the Simulation Event Analysis module of DESMOND. All frames were aligned with the starting structure prior to the calculations. Protein–ligand contacts were analyzed using the Simulation Interactions Diagram module of DESMOND.

2.6. Binding energy calculation

The Prime MM-GBSA (molecular mechanics based generalized born/surface area) model of the Schrödinger suite was used to calculate the free energy of binding of the mammalian tyrosinase–inhibitor complex from the MD trajectory snapshots.⁴⁶ The binding free energy (ΔG_bind) was evaluated as:

ΔG_bind = ΔE_MM + ΔG_solv + ΔG_SA

where ΔE_MM is the difference in the minimized energies between the tyrosinase–inhibitor complex and the sum of the energies of the free tyrosinase and inhibitor. ΔG_solv is the difference between the GBSA solvation energy of the tyrosinase–inhibitor complex and the sum of the solvation energies of free tyrosinase and inhibitor. ΔG_SA is the difference between the surface area energies of the complex and the sum of the surface area energies of free tyrosinase and inhibitor. Based on the MM/GBSA binding free energy, we developed a correlation model between the calculated binding free energy and the experimental IC₅₀ values for mTyr.

3. Results and discussion

3.1. Computational prediction of mTyr and hTyr structures

The mTyr and hTyr structures, for which X-ray structures are not available, were built using the homology modeling method (Fig. 2). To date, X-ray structures have been determined for tyrosinase isolated from B. megaterium (bmTyr), S. castaneoglobisporus (scTyr), A. bisporus (abTyr), and A. oryzae (aoTyr). The bacterial and fungal tyrosinases are cytosolic proteins, whereas mammalian tyrosinases are integral membrane proteins anchored in the melanosomal membrane. This structural difference is caused by differently sized N-terminal signal peptides and C-terminal tails, whereas the catalytic region is relatively well conserved across species.^47,48

Fig. 2 Homology models of mammalian tyrosinase. (a) mTyr (green) and (b) hTyr (orange) structure is superimposed with scTyr (gray) and bmTyr (yellow). Binding site of tyrosinase is well conserved among the species, while the rest of the loop structure is altered. The binding sites of (c) mTyr and (d) hTyr are also displayed. Two copper ions (red) are coordinated by six histidine residues of mTyr and hTyr, respectively, and the coordinated positions are comparable to those of scTyr and bmTyr.

Six amino acid sequences were aligned using the pairwise sequence alignment tool of the European Bioinformatics Institute and sequence identities and similarities were evaluated. We performed global sequence alignments among the six sequences using the Needleman–Wunsch algorithm (Table 1). Amino acid sequence identities and similarities are more common within the same species and less common between different species, and mammalian species are more closely related than others are. Thus, we predict that mTyr and hTyr have similar structures and related functions. Next, local alignment was performed using the Smith–Waterman algorithm, and higher numbers than those for global alignment were obtained (Table 2). This result indicates that the function of the tyrosinase catalytic site is conserved. Notably, the six histidine residues coordinating the copper ion in the active site are highly conserved.

Table 1 The Needleman–Wunsch global alignment for six tyrosinases. The first number is % identity and the second number is % similarity at amino acid level

	mTyr	hTyr	bmTyr	scTyr	abTyr	aoTyr
a % Identity at amino acid level.b % Similarity at amino acid level.
mTyr		85.7^a/92.1^b	18.8/28.4	17.3/28.0	13.5/21.4	15.3/22.8
hTyr			19.3/29.0	18.4/28.3	13.6/20.9	14.5/23.1
bmTyr				41.4/54.6	17.1/24.7	14.3/23.2
scTyr					15.1/22.2	12.3/20.5
abTyr						23.8/38.2
aoTyr

---

Table 2 The Smith–Waterman local alignment for six tyrosinases. The first number is % identity and the second number is % similarity at amino acid level

	mTyr	hTyr	bmTyr	scTyr	abTyr	aoTyr
a % Identity at amino acid level.b % Similarity at amino acid level.
mTyr		85.7^a/92.1^b	29.4/43.8	29.0/45.2	22.9/34.3	20.1/30.3
hTyr			29.8/44.5	29.7/45.2	22.8/33.8	20.8/33.5
bmTyr				43.7/57.7	24.9/36.8	20.9/34.6
scTyr					25.3/35.8	20.8/31.9
abTyr						24.7/39.6
aoTyr

---

In our previous study, (S)-(+)-decursin and its analogues were found to effectively inhibit melanin formation in mouse melanoma cells.⁸ A mTyr structure was built to identify interactions in the protein–inhibitor complex using the Prime homology modeling program. We also built an hTyr structure, which had high sequence similarity with mTyr to predict binding interactions. The mTyr and hTyr amino acid sequences were obtained from the NCBI database and run through PSI-BLAST searching the PDB. 1WX2 from scTyr and 3NM8 from bmTyr in the PDB were determined to be appropriate template structures for homology modeling. Homology with mammalian tyrosinase around the active site is higher for bacterial than for fungal tyrosinase. Therefore, the structural analysis concentrating on the helix bundle and adjacent loop contributing the copper-binding residues was based on bacterial tyrosinase, while the remaining structure was generated based on fungal tyrosinase.

The predicted structures of mTyr and hTyr consist of 330 amino acids of the full sequence of 533 and 529 amino acids, respectively. The aligned region was from 113 to 433. Cordes et al. constructed a number of differently truncated soluble hTyr variants and found that enzyme forms lacking the C-terminal domain and transmembrane region are readily secreted by HEK cells and that K_m-values of the variants were not significantly different from those of the wild-type enzyme.⁴⁹ We removed the residues beyond 443, which make up the transmembrane domain, and the N-terminal signal peptides before 133 because they do not influence catalytic activity. After building the homology models, RAMPAGE was used to evaluate their validity with a Ramachandran plot. We found that 94.5% of mTyr residues and 96.6% of hTyr residues were within either favored or allowed regions (Fig. S1†). Since the few remaining residues were located far from the active site, the homology models appeared to be suitable for studying ligand–protein interactions.

The generated mTyr and hTyr models were compared to identify differences of binding site residues between the mammalian species (Fig. S2†). In the binding site, two residues differed between species. Instead of Thr75 and Arg222 of mTyr, there are Ala75 and Lys222 in the binding site of hTyr. Residue 75 is distant from the binding site and has little effect on ligand binding compared with residue 222. The residues 222 in mTyr and hTyr are both located near the edge of the binding sites and form intra-hydrogen bonds with Glu91. The orientation of the side chain of Glu91 was altered according to residue 222. The side chain of Arg222 in mTyr was oriented towards the binding site and followed the side chains of Glu91 to the binding site. In contrast, in hTyr, the side chain of Lys222 faced towards the side chain of Glu91 and formed a hydrogen bond. This structural difference influenced to a significant extent the ligand binding. Thus, we found that residue 222, which is located near the binding site, could influence the interaction with a ligand between mammalian species.

3.2. Binding mode prediction of tyrosinase inhibitors with mTyr and hTyr

Molecular docking was used to predict the binding mode of S-(+)-decursin and its analogues, which are potent inhibitors of mouse melanoma cells, to the mTyr homology structure to understand the activities of the compounds (Fig. 3). In addition, we identified the binding mode of the ligands to the hTyr homology structure to predict the activity for human treatment. Arbutin, which showed low activity in our experimental study, was docked to both mammalian tyrosinase structures as a reference compound. Docking predicted the binding positions of (S)-(+)-decursin and its analogues in two different ways: one with the two ester groups at positions 1 and 9 pointing towards Arg222 of mTyr and Lys222 of hTyr, and the other with the ester groups pointing in the exactly the opposite direction (Fig. 3a). This difference in binding poses is significant, as it leads to different inhibitor contacts with the residues around the binding site; thus, the binding mode must be determined. Our previous experimental and computational results^8,17 showed that substitution of the hydroxyl group at position 7 of (+)-decursinol by ester derivatives plays a key role in the interaction with tyrosinase. Therefore, the second binding mode enabling the formation of hydrogen bonds with the ligands was selected for further analysis.

The model structures of mTyr and hTyr were constructed from the template structures 1WX2 and 3NM8. 1WX2 was co-crystallized with a caddie protein, which carries the copper ions for tyrosinase. The caddie protein places a tyrosine residue at the tyrosinase active site. Thus, the conformational states of mTyr and hTyr would differ when bound to small molecules. We performed minimization of the protein–ligand complex after docking. However, the active site including the copper ions and histidine residues were destabilized by MM minimization with Schrödinger's IMPACT using the OPLS_2005 force field. It is well known that force field-based MM minimization is unable to model metal ion-including structures. Thus, it was necessary to use QM level theory to correctly predict the tyrosinase active site structure. We used the QM/MM method because modeling the entire protein–ligand complex with QM would be too costly. Treating only the active site with QM and the remaining structure with MM allowed us to model the tyrosinase active site in which metal ions are present with sufficient accuracy in a reasonable amount of time. We used QSite, Schrödinger's QM/MM calculation module, designating the active site as the QM region.

The binding mode of tyrosinase inhibitors with mTyr and hTyr was analyzed after QM/MM minimization (Fig. 3b–e). The compounds, which all contain an ester group at the 7-position, formed a hydrogen bond with Arg222 in mTyr, except compound 1, which contains a hydroxyl group at the 7-position. In addition, all compounds, including compound 1, formed a hydrogen bond with Lys222 of hTyr. The arginine residue was conformationally restricted compared to the lysine residue, and this mutation between mammalian tyrosinases influenced hydrogen bond formation with compound 1. The substitutions at the 7-position of the compounds added additional hydrophobic or aromatic groups, which interact with hydrophobic residues, particularly Ile256 and Val265, in the tyrosinase binding site. The carbonyl group at 2-position of the compounds interacts with copper ions.

Fig. 3 Binding modes of inhibitors in mammalian tyrosinase structures. (a) Docking predicted the binding poses in two different ways: one with the two ester groups at positions 1 and 9 pointing towards Arg222 of mTyr and Lys222 of hTyr; and the other with the ester groups pointing exactly the opposite direction. (b) Compound 1 and (c) compound 2–5 bound to mTyr after QM/MM minimization. The compound 1 does not form hydrogen bond with Arg222, while the compound 2–5 does. (d) Compound 1 and (e) compound 2–5 bound to hTyr after QM/MM minimization. Mutations from Arg222 in mTyr to Lys222 in hTyr influence the activity of compound 1. All the compounds form hydrogen bonds with Lys222. The compounds are colored as follows: compound 1: gray, compound 2: turquoise, compound 3: dark green, compound 4: yellow and compound 5: purple. Residues which are important for binding the compounds are indicated in bold and hydrogen bonds are colored in purple dotted line.

Arbutin also docked to both tyrosinase structures and QM/MM was minimized as shown in Fig. S3.† Arbutin formed a hydrogen bond with Arg222 in mTyr and Lys222 in hTyr. In addition, the sugar part of arbutin formed a few hydrogen bonds with Glu91 and Asn252 in mTyr and Ser263 in hTyr.

Based on these results, we found that substituting chemical groups at the 7-position created highly potent tyrosinase inhibitors that formed hydrogen bonds with Arg222 in mTyr and Lys222 in hTyr and hydrophobic interactions with the hydrophobic residues Ile256 and Val265 in the binding site. Although mTyr and hTyr show high sequence identities, the inhibitors including arbutin bound differently to both structures because of their structural differences. Thus, we expect that potent inhibitors of mTyr can be altered for hTyr.

3.3. Molecular dynamics simulations of mTyr and hTyr with inhibitors

We performed MD simulations using mammalian tyrosinase and inhibitor complexes to analyze binding interactions and identify different structural dynamics between mTyr and hTyr. To clarify the dynamic stability of the system, RMSD values were calculated using the starting structures as a reference (Fig. 4). The RMSD plot indicated that the RMSD of the backbone C_α atoms in the complexes increased for 10 ns and then remained stable until the end of the simulation. The results showed that the trajectories of the MD simulations of the complexes after equilibrium were reliable.

Fig. 4 The RMSD of the backbone C_α atoms of (a) mTyr and (b) hTyr in inhibitor complex structures against their initial minimized complex structures.

The binding interactions of the mTyr–inhibitor complexes were identified by trajectory analysis (Fig. 5 and S4†). Compound 1 continuously bound to the active site of mTyr. The compound showed not only a hydrophobic interaction with His90 and His255, but also an electrostatic interaction with copper ions. Thus, the scaffold and carbonyl group at the 2-position of the compound is important for the binding of mTyr. Compound 2, which contains an ester and isobutylene group at the 7-positon, has additional interactions with Arg222, Leu239, Asn252, and Ile256. Hydrogen bond and water bridge interactions were observed between the compound and Arg222. The isobutylene group was found to have a hydrophobic interaction with Leu239 and Ile256. Compound 3, which contains the but-1-ene group, interacts with Ile256, Ser263, and Val265 as compared with compound 1. Compound 4 contains an acrylic acid derivative group and frequently interacts with Ile256 and Val265. The potency of compound 4 can be explained by the interaction of its acrylic acid derivative group, which fully covers the Val265 residue. Compound 5 contained a phenyl derivative group and showed a pi–cation interaction with Lys196. Arbutin mainly interact with Glu233 by forming hydrogen bonds in the binding site. Based on these results, we determined that arbutin and compound 1 sufficiently bound to the binding site of mTyr and that the substitution of the hydroxyl group at the 7-position of (+)-decursinol confers potency to the inhibitors, which interact with the residues around the binding site of mTyr.

Fig. 5 Binding poses of tyrosinase inhibitors in mTyr. Each complex structure is last snapshot from the molecular dynamics simulations. The compounds are colored as follows: compound 1: gray, compound 2: turquoise, compound 3: dark green, compound 4: yellow, compound 5: purple and arbutin: aquamarine. Residues which are important for binding of the compounds during molecular dynamics simulations are indicated in bold and shown in CPK model.

The binding interaction of hTyr–inhibitor complexes were also identified by trajectory analysis (Fig. 6 and S5†). Compound 1, which continuously bound to the active site of hTyr, showed similar binding interactions to those with mTyr. A hydrophobic interaction was observed with His90 and Val265, and an electrostatic interaction with the copper ions. Compound 2 is pulled out of the binding site; this binding difference is caused by a mutation from Arg in mTyr to Lys in hTyr. The lysine residue is more flexible than the arginine residue, and thus it can interact with other residues. The potency of compound 2 may be reduced by loose hydrogen bonding. Arbutin is also pulled out of the binding site because of the structural differences between mammalian species. The phenol group of arbutin forms a hydrogen bond with Lys222, while the rest of the sugar group hydrogen bonds with Ser263. The binding interactions with hTyr are relatively reduced compared with mTyr, which forms more hydrogen bonds with arbutin. Compound 3 interacts with Ile256 and Val265, compound 4 has more frequent interactions with Val265, and compound 5 extensively interacts with the surface of the structure.

Fig. 6 Binding poses of tyrosinase inhibitors in hTyr. Each complex structure is last snapshot from the molecular dynamics simulations. The compounds are colored as follows: compound 1: gray, compound 3: dark green, compound 4: yellow and compound 5: purple. Residues which are important for binding of the compounds during molecular dynamics simulation are indicated in bold and shown in CPK model.

Based on the results, we found that S-(+)-decursin and its analogues are bound to mTyr and hTyr via electrostatic and hydrophobic interactions with copper ions and residues around the binding site. In both mammalian tyrosinase structures, Ile256 and Val265 are important for ligand binding. Although the hydrogen bond with Arg222 in mTyr or Lys222 in hTyr is important for the potency of compounds in docking, the interaction only affects compound 2 in MD simulations. hTyr amino acid sequences showed high identity to mTyr. We expect that potent inhibitors for mTyr have a similar inhibitory effect on hTyr. All compounds except compound 2 showed similar binding interactions with hTyr. Compound 2 is pulled out of the binding site of hTyr because a mutation occurred from Arg in mTyr to Lys in hTyr. Arbutin is also pulled out because of the structural difference between mammalian species. The arginine residue is more restricted than the lysine residue, enabling compound 2 to form a stable hydrogen bond with mTyr. This structural difference may explain the different inhibitory activities of the species. Overall, further studies are required to identify new inhibitors of hTyr to treat dermatological disorders.

3.4. Binding free energies calculated by MM-GBSA

Binding free energy was calculated for the mTyr–inhibitor complexes using MM/GBSA methods (Table 3). The 625 snapshots were extracted at a time interval of 4.8 ps from the last 3 ns of MD trajectories to analyze the binding free energy. As shown in Table 3, the binding free energies of the compounds were −30.65, −44.96, −58.91, −50.64, −50.81, and −16.67 kcal mol⁻¹, respectively. Furthermore, the ranking of the experimental binding free energies is consistent with our predictions, showing that the current analysis by MM-GBSA is reliable (Fig. 7). According to the energy components of the binding free energies, the major favorable contributors to inhibitor binding were van der Waals energies. That van der Waals energies mostly drive the binding of the four inhibitors to mTyr agrees well with previous experimental analyses. It also suggests that optimization of the van der Waals interactions between inhibitors and mTyr may lead to the development of potent small molecule inhibitors of mTyr.

Table 3 Binding free energy components for the mouse tyrosinase–inhibitor complexes determined by using the MM-GBSA method. The binding free energy (ΔG_bind) includes four energy terms: electrostatic contribution (ΔE_ele), van der Waals contribution (ΔE_vdW), polar desolvation term (ΔG_GB) and nonpolar desolvation term (ΔG_SA)

Compound	ΔGbind	ΔEele	ΔEvdW	ΔGGB	ΔGSA	pIC50
Compound 1	−30.65	−10.08	−32.26	33.12	−21.43	3.39
Compound 2	−44.96	−12.06	−38.48	35.39	−29.81	4.15
Compound 3	−58.91	−2.12	−40.83	14.01	−29.97	4.2
Compound 4	−50.64	0.04	−41.16	17.71	−27.23	4.21
Compound 5	−50.81	−6.38	−41.63	28.39	−31.19	4.1
Arbutin	−16.67	−19.38	−31.72	51.56	−17.13	3.21

---

Fig. 7 Correlation between the experimental pIC₅₀ and the MM-GBSA binding energies of mTyr with arbutin, (S)-(+)-decursin and its analouges.

4. Conclusions

Herein, we predicted the binding interactions between mammalian tyrosinase and its inhibitors using computational methods. Amino acid sequence analysis showed that the catalytic region of tyrosinase is highly conserved across species and that mammalian species showed higher amino acid sequence similarity than other species. In our previous study, (S)-(+)-decursin and its analogues effectively inhibited melanin formation in mouse melanoma cells. We generated a three-dimensional homology model of mouse tyrosinase to analyze the binding interactions with its inhibitors. Since human tyrosinase has high sequence similarity, its structure was used to predict binding interactions and identify structural differences between mammalian species, which influences the activity of inhibitors. Molecular docking simulations were used to select possible binding modes, and QM/MM minimization was performed on mammalian tyrosinase and inhibitor complexes. We found that substitution of the hydroxyl group at the 7-position of (+)-decursinol by ester derivatives plays key role in the interaction. The inhibitors can interact with the hydrophobic and aromatic residues around the binding site of tyrosinase structures. MD simulations were performed to study the dynamic behaviors of the complexes and to identify different structural dynamics between the different tyrosinases and inhibitors. The results showed that residual differences might explain the differences in inhibitory activities. In both mammalian tyrosinase structures, Ile256 and Val265 were important for ligand binding. According to residue 222 in mTyr or hTyr, ligands can selectively bind to the structure. Thus, this difference is significant for human treatment. Glu91, which forms a hydrogen bond with residue 222, may form a hydrogen bond with the ligand. Binding free energies were calculated for the mouse tyrosinase–inhibitor complexes using the MM/GBSA method. Biological activities were related to the calculated binding free energies of the inhibitors. Our study provides a foundation for the development of tyrosinase inhibitors for treating dermatological disorders.

Acknowledgements

This work was supported by NRF Grants No. 2013R1A2A2A01067638 and 2012M3C1A6035 362. We also thank Schrödinger, LLC for providing software for our research.

References

1. M. S. Marks and M. C. Seabra, The melanosome: membrane dynamics in black and white, Nat. Rev. Mol. Cell Biol., 2001, 2(10), 738–748.
2. A. Slominski, et al., Melanin pigmentation in mammalian skin and its hormonal regulation, Physiol. Rev., 2004, 84(4), 1155–1228.
3. H. Decker, T. Schweikardt and F. Tuczek, The first crystal structure of tyrosinase: all questions answered?, Angew. Chem., Int. Ed. Engl., 2006, 45(28), 4546–4550.
4. M. Rolff, et al., Copper-O₂ reactivity of tyrosinase models towards external monophenolic substrates: molecular mechanism and comparison with the enzyme, Chem. Soc. Rev., 2011, 40(7), 4077–4098.
5. M. Sendovski, et al., First structures of an active bacterial tyrosinase reveal copper plasticity, J. Mol. Biol., 2011, 405(1), 227–237.
6. C. Olivares and F. Solano, New insights into the active site structure and catalytic mechanism of tyrosinase and its related proteins, Pigm. Cell Melanoma Res., 2009, 22(6), 750–760.
7. M. Criton and V. Le Mellay-Hamon, Analogues of N-hydroxy-N′-phenylthiourea and N-hydroxy-N′-phenylurea as inhibitors of tyrosinase and melanin formation, Bioorg. Med. Chem. Lett., 2008, 18(12), 3607–3610.
8. K. Lee, et al., Synthesis of (S)-(+)-decursin and its analogues as potent inhibitors of melanin formation in B16 murine melanoma cells, Eur. J. Med. Chem., 2010, 45(12), 5567–5575.
9. S. Takahashi, et al., Structural insights into the hot spot amino acid residues of mushroom tyrosinase for the bindings of thujaplicins, Bioorg. Med. Chem., 2010, 18(22), 8112–8118.
10. M. T. H. Khan, Molecular design of tyrosinase inhibitors: A critical review of promising novel inhibitors from synthetic origins, Pure Appl. Chem., 2007, 79(12), 2277–2295.
11. C.-J. Pei, et al., Inhibition of tyrosinase by gastrodin: An integrated kinetic-computational simulation analysis, Process Biochem., 2013, 48(1), 162–168.
12. R. J. Deeth and C. Diedrich, Structural and mechanistic insights into the oxy form of tyrosinase from molecular dynamics simulations, JBIC, J. Biol. Inorg. Chem., 2010, 15(2), 117–129.
13. Y. X. Si, et al., An integrated study of tyrosinase inhibition by rutin: progress using a computational simulation, J. Biomol. Struct. Dyn., 2012, 29(5), 999–1012.
14. W. J. Hu, et al., Kinetic, structural and molecular docking studies on the inhibition of tyrosinase induced by arabinose, Int. J. Biol. Macromol., 2012, 50(3), 694–700.
15. Z. J. Wang, et al., The effect of fucoidan on tyrosinase: computational molecular dynamics integrating inhibition kinetics, J. Biomol. Struct. Dyn., 2012, 30(4), 460–473.
16. K. Bagherzadeh, et al., A new insight into mushroom tyrosinase inhibitors: docking, pharmacophore-based virtual screening, and molecular modeling studies, J. Biomol. Struct. Dyn., 2015, 33(3), 487–501.
17. K. Lee, et al., 3D-QSAR Study of Melanin Inhibiting (S)-(+)-Decursin and its Analogues by Pharmacophore Mapping, Bull. Korean Chem. Soc., 2012, 33(1), 149–152.
18. D. A. Benson, et al., GenBank, Nucleic Acids Res., 2013, 41, D36–42.
19. E. W. Sayers, et al., Database resources of the National Center for Biotechnology Information, Nucleic Acids Res., 2012, 40, D13–D25.
20. P. Rice, I. Longden and A. Bleasby, EMBOSS: the European Molecular Biology Open Software Suite, Trends Genet., 2000, 16(6), 276–277.
21. H. M. Berman, et al., The Protein Data Bank, Nucleic Acids Res., 2000, 28(1), 235–242.
22. S. F. Altschul, et al., Gapped BLAST and PSI-BLAST: a new generation of protein database search programs, Nucleic Acids Res., 1997, 25(17), 3389–3402.
23. Y. Matoba, et al., Crystallographic evidence that the dinuclear copper center of tyrosinase is flexible during catalysis, J. Biol. Chem., 2006, 281(13), 8981–8990.
24. Schrödinger, Prime, version 3.5, LLC, New York, NY, 2014.
25. S. C. Lovell, et al., Structure validation by Calpha geometry: phi, psi and Cbeta deviation, Proteins, 2003, 50(3), 437–450.
26. R. A. Friesner, et al., Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy, J. Med. Chem., 2004, 47(7), 1739–1749.
27. T. A. Halgren, et al., Glide: a new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening, J. Med. Chem., 2004, 47(7), 1750–1759.
28. Schrödinger, MacroModel, version 9.9, LLC, New York, NY, 2012.
29. Schrödinger, LigPrep, version 2.5, LLC, New York, NY, 2012.
30. R. B. Murphy, D. M. Philipp and R. A. Friesner, A mixed quantum mechanics/molecular mechanics (QM/MM) method for large-scale modeling of chemistry in protein environments, J. Comput. Chem., 2000, 21(16), 1442–1457.
31. D. M. Philipp and R. A. Friesner, Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide, J. Comput. Chem., 1999, 20(14), 1468–1494.
32. A. D. Becke, Density-functional thermochemistry. III. The role of exact exchange, J. Chem. Phys., 1993, 98(7), 5648–5652.
33. C. Lee, W. Yang and R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B, 1988, 37(2), 785–789.
34. P. J. Stephens, et al., Ab initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields, J. Chem. Phys., 1994, 98(45), 11623–11627.
35. P. J. Hay and W. R. Wadt, Ab initio effective core potentials for molecular calculations. Potentials for K to Au including the outermost core orbitals, J. Chem. Phys., 1985, 82(1), 299–310.
36. J. L. Banks, et al., Integrated Modeling Program, Applied Chemical Theory (IMPACT), J. Comput. Chem., 2005, 26(16), 1752–1780.
37. Desmond Molecular Dynamics System, version 3.7, D. E. Shaw Research, New York, NY, 2014; Schrödinger, Maestro-Desmond Interoperability Tools, version 3.7, New York, NY, 2014.
38. K. J. Bowers, et al. Scalable algorithms for molecular dynamics simulations on commodity clusters, in SC 2006 Conference, Proceedings of the ACM/IEEE, IEEE, 2006.
39. W. L. Jorgensen, et al., Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79(2), 926–935.
40. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81(1), 511–519.
41. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A, 1985, 31(3), 1695–1697.
42. G. J. Martyna, D. J. Tobias and M. L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys., 1994, 101(5), 4177–4189.
43. B. Graham, et al., An examination of the binding behavior of histidine-containing peptides with immobilized metal complexes derived from the macrocyclic ligand, 1,4,7-triazacyclononane, JBIC, J. Biol. Inorg. Chem., 2007, 12(1), 11–21.
44. D. F. Hansen and J. J. Led, Determination of the geometric structure of the metal site in a blue copper protein by paramagnetic NMR, Proc. Natl. Acad. Sci. U. S. A., 2006, 103(6), 1738–1743.
45. O. Wise and O. Coskuner, New force field parameters for metalloproteins I: Divalent copper ion centers including three histidine residues and an oxygen-ligated amino acid residue, J. Comput. Chem., 2014, 35(17), 1278–1289.
46. Z. Yu, M. P. Jacobson and R. A. Friesner, What role do surfaces play in GB models? a new-generation of surface-generalized born model based on a novel gaussian surface for biomolecules, J. Comput. Chem., 2006, 27(1), 72–89.
47. J. C. Garcia-Borron and F. Solano, Molecular anatomy of tyrosinase and its related proteins: beyond the histidine-bound metal catalytic center, Pigm. Cell Res., 2002, 15(3), 162–173.
48. N. Wang and D. N. Hebert, Tyrosinase maturation through the mammalian secretory pathway: bringing color to life, Pigm. Cell Res., 2006, 19(1), 3–18.
49. P. Cordes, et al., Expression in non-melanogenic systems and purification of soluble variants of human tyrosinase, Biol. Chem., 2013, 394(5), 685–693.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra09365e

---