DhkM amino-acid sequence obtained in DictyBase contains 2388 amino acids. It also contains seven transmembrane domains (residue number 550–570, 589–609, 645–665, 679–699, 896–916, 953–973, and 1025–1045), and two response regulatory domains (residue number 1541–1656 and 2262–2383) and the binding site was considered to locate between these two regions (residue number 1046–1540) from a biological point of view. The functional residues in this core part (residue number 1073–1498) which kept very high identities was just between the transmembrane domain and the response regulatory domain which was recognized as a histidine kinase domain.^11,15,16 Therefore, the region of residue number 1073–1498 was considered as the core part of DhkM which includes the binding site for ligand.

The similarity between DhkM and the template was confirmed by Clustal Omega.²² Histidine kinase 4JAU.1.A¹⁷ was used as the template for the structure prediction. Fifty seven identical positions and eighty six similar positions were defined among two hundred and thirty eight residues, which were applied to predict the structure of histidine kinase domain in the core part. In particular, the arrangement for the eighty six residues reached the characteristic arrangement of the histidine kinase domain by the template. These residues were supposed to keep the properties of histidine kinase. From this result, the histidine kinase domain contained in the core, could take a tertiary structure similar to the histidine kinase domain in the template with reliability, because of the conservation of the same domains. Therefore, the employed template in the present work should have reliability provided by the experimental measurement. Of course, the completely randomized structure could be started using the random seeds, and moreover the reasonable structure would be determined by replica-exchange approach. But some reliable data obtained by the experimental measurement has been already known, which could also provide a reasonable structure. The predicted structure by SWISS-MODEL¹¹ is shown in Figure 2. As seen in the figure, the predicted structure of DhkM contains a long stretching chain except the core. Then the determination of initial structure for docking simulations and energy analysis could be performed by the MD simulations.

Figure 2. The core part of the protein structure which includes binding sites depicted as several colored spheres generated by the protomol. Colored dummy atoms represent predicted binding pockets.

The determination of initial structure of the core part of DhkM and the structural refinement were carried out by the production MD-simulation of one hundred nanoseconds and the final structure is shown in Figure 3. The time-profile of root mean square deviation (RMSD) for the early stage of the MD simulation, which is considered to be the most important part in this case, has also been calculated to verify the structural convergence. RMSD was calculated for the backbone: all carbon, alpha carbon and nitrogen atoms, of the core part of DhkM. Best-fit RMSD was performed, in which each structure is rotated and translated so as to minimize the RMSD to the reference structure. In this case, the first frame was set to be the reference. The result is shown in Figure 4, and it can be seen that the structure of protein changed as the simulation underwent the long-time-step MD simulation, the folding structure appeared gradually and was used for docking simulations with DIF-1. Docking scores were obtained at several different time steps.

Figure 3. The refined protein structure of DhkM after the long-time-step (100 ns) simulation.

Figure 4. The time-profile of root mean square deviation (RMSD) for the DhkM structural changes.

Since DhkM is a kind of histidine kinase which contains an ATP-binding site and the function is to induce the cell death, the hypothesis is as follows;

i) From a biological point of view, the receptor in a stable state should be inactive because the activation costs cell death in nature.

ii) Once the cell-death process is required, ATP-binding process provides the necessary energy for receptor-activation while the receptor should be inactive without ATP-binding in a stable state.

Therefore, it is supposed that the structures which contain activated binding sites can only be observed at the very early stage of the MD simulations in a comparatively unstable state. Furthermore, the activation and inactivation processes would be controlled by the structural changes of the binding site, which could be expressed in characteristic changes of the RMSD value.

As predicted, unnatural protein structure with an incompact chain disappeared at around 3–4 ns (Figure 5a) the folding structure including six binding pockets emerged (Figure 5b), which involved high docking scores in binding simulation as listed in Table 1.

Figure 5. (a) Structures obtained from MD simulation. (b) The structure obtained at 3.2 ns containing six binding pockets, which are represented by colored dummy atom.

Table 1. Docking properties obtained by Surflex-Dock. The employed structure was produced by the MD simulation around 3,200 ps. Total score was expressed in terms of −log (Kd) units to evaluate the docking results.

In case of docking simulation for DIF-1, a higher total score, lower crash value, and appropriate pocket size would be desirable. CYAN, MAGENTA and YELLOW pockets clearly show the potential binding sites for DIF-1 as the result of the comprehensive evaluations.

Further analysis for binding models which showed high scores have been carried out in terms of structural arrangement and hydrophilic interaction or hydrophobic effect between DhkM and DIF-1. The stability of the binding model was confirmed by calculating the binding energy with the MD simulation.

6.3.1 Structural Arrangement Analysis:

Hydrophilic or hydrophobic properties were mapped by lipophilic potential with the MOLCAD option of SYBYL¹³ and are shown in Figure 6. In the bar of the figures, red tendency shows high lipophilic potential and blue tendency shows low lipophilic potential. The result for pocket CYAN is shown in Figure 6a. The ligand and receptor fit each other quite well in the structure. The side chain of DIF-1 is located just inside the pocket like a key and there are six hydrogen bonds formed between DIF-1 and the binding site. It is well known that there is a trend in which hydrogen bond decides the binding strength between ligand and receptor. The six hydrogen bonds appearing in CYAN pocket would ensure the stability of the system. An analysis was also performed on the molecular surface which was determined by the Connolly surface. The side chain of DIF-1 and binding site fit each other so well that pocket CYAN was filled up with the side chain of DIF-1. Thus, the pocket CYAN could be an appropriate binding site for DIF-1.

Figure 6. Lipophilic potentials and molecular surfaces (a) for pocket CYAN (b) for pocket MAGENTA and (c) for pocket YELLOW (Ligand drawn in sticks). CYAN showed the best result in docking simulations.

The result for pocket MAGENTA is shown in Figure 6b. The ligand and receptor did not fit each other quite as well in this structure. The side chain of DIF-1 seems to be pulled toward the pocket but is mostly located outside it. Since there is only a single hydrogen bond formed, the interaction is expected to be quite weak. A remarkable point is that the shape of ligand DIF-1 is quite warped which means the molecule had to take such an unstable structure in order to bind with the site. Subsequently as a result of analysis of the molecular surface, the side chain of DIF-1 did not fit the binding site as well as in the case of pocket CYAN. Thus, it is considered that pocket MAGENTA could not be an appropriate binding site for DIF-1.

Pocket YELLOW showed almost the same result as pocket MAGENTA and is shown in Figure 6c; it is considered to not be an appropriate binding site for DIF-1 either. Generally, pocket CYAN (LEU111 ASP115 ARG150 VAL151 TYR351 VAL354) shows the best result among these three potential binding sites and a further analysis from the point of view of hydrophilic interaction and hydrophobic effect has been carried out.

6.3.2 Analysis of Hydrophilic Interaction and Hydrophobic Effect:

It has been reported that the presence of halogen atoms in a molecule increases its lipophilicity and hydrophobicity.²³ The structure of ligand molecule DIF-1 was optimized by ab initio calculation and it contains two lipophilic areas (around chlorine, Cl), one hydrophilic area (around oxygen, O) and one side chain which appeared slightly lipophilic (Figure 7a). On the other hand, there were two lipophilic areas outside the pocket, one hydrophilic area at the entrance of the pocket and a slightly lipophilic region inside the pocket (Figure 7b.i). As shown in Figure 7b.ii, the lipophilic and hydrophilic areas of DIF-1 and pocket CYAN matched each other well in the structure; therefore the binding model of DIF-1 and pocket CYAN is expected to be stable because of the hydrophilic interaction and hydrophobic effect and it can be argued that the binding model on pocket CYAN is suitable.

Figure 7. (a) Two lipophilic areas, one lipophilic region and one hydrophilic area were found in DIF-1. (b) Lipophilic and hydrophilic areas of the binding pocket fitted with the ligand well.

Our goal is to verify whether the whole system becomes more stable or not after binding, and one of the approaches is to evaluate the stability of the binding model. The binding energy analysis is done based on free energy ΔG obtained by the molecular mechanics Poisson–Boltzmann surface area (MMPBSA) method.¹⁰ MMPBSA typically gives excessively large energies but is shown by Genheden and Ryde to be useful for post-processing of docked structures or to rationalize observed differences.²⁴ In the PB/SA part, ionic strength (in mM) for the PB equation was set to be 0.200, and the modern nonpolar solvent model was applied. The time interval of the frame was 2 ps. The binding energy was estimated to be −11.09 kcal/mol. As a result, the whole system is found to be quite stable after binding and pocket CYAN in DhkM could be an appropriate and suitable binding site for DIF-1.

Based on the RMSD curve, the characteristic structure changes were expected to occur at 21–23 ns during the MD simulation up to 100 ns. However, an interval of twenty nanoseconds is said to be long enough for the structure change. The spatial area and the configuration of the pocket constituted by almost the same residues with pocket CYAN at 3 ns can still be observed at 4–21 ns but closes gradually. The structures around 22 ns were also investigated in the present work (Figure 8a), and the pocket was actually closed and the ligand could not get in and bind to the site in the docking simulations (Figure 8b). At 62 ns, the pocket CYAN could not be observed anymore and the protein was expected to be completely inactive.

Figure 8. (a) Pocket constituted by almost the same residues with pocket CYAN, LEU 111 ASP113 PRO117 VAL151 TYR351 VAL354 LEY421 remaining at 6,000 ps. (b) The pocket was closed. The ligand cannot get in to bind with the site (Ligand drawn in sticks).

It is supposed that the structures for activated binding sites could be observed at an early stage of the MD simulations. Therefore, our results show that the pocket containing the binding site could only be activated; the emergence as the system is being stabilized is thought to be quite reasonable for the binding process.

The prediction of the mechanism of the binding process is based on the result of docking simulations. Two lipophilic areas outside the pocket, one hydrophilic area at the entrance of the pocket and a lipophilic region inside the pocket are observed in our simulation. It is considered that the lipophilic areas play the important role of ligand-recognition and translocation which agrees with the results in drug design research.^25,26 On the other hand the hydrophilic site at the entrance of the pocket and a lipophilic region inside the pocket make the system stable and have relevance to the function of the receptor. The verification is performed using the molecule DIF-2 and DIF-3, shown in Figure 9, which have a similar molecular structure to DIF-1.

Figure 9. Molecular structures of DIF-2 (left) and DIF-3 (right).

DIF-2 contains a shorter alkyl chain structure with a butyl group than DIF-1 having a pentyl group, and DIF-3 takes the structure of monochloro-DIF-1. In our prediction, since DIF-2 contains a shorter alkyl chain structure than DIF-1, DIF-2 could be recognized and bind to DhkM with a lower stability than DIF-1 which may lead to a lower bioactivity. On the other hand, in contrast to the dichloride structure, lipophilic area is lost in DIF-3 decreasing the lipophilicity of the molecule. As a result, DIF-3 could not be normally recognized and the binding process would occur with a low probability.

Docking simulations of DIF-2 and DIF-3 have also been carried out under the same conditions as DIF-1. Results for the three MD simulation trajectories of DIF-1 and DIF-2 are summarized in Table 2. DIF-2 showed lower score than DIF-1 such as 3.9492 in trajectory 1, but it was still able to constitute a well binding model. The ligand and receptor also match each other quite well in molecular structure (Figure 10). DIF-3 could bind to the protein in the docking simulation and the score was 3.2542 in trajectory 1. But as shown in Figure 11, DIF-3 is unrealistically warped in order to bind with the site, which is not observed in the case of DIF-1 and DIF-2. The respective bond-length of the C1-C2 and C1-O bonds were 0.35135 and 4.27622 Å; the ring structure of benzene was extremely distorted with the dissociation of C-OH bond. In general, the biosystem should keep the ring structure of benzene, which means DIF-3 is quite unrealistic as a result of docking simulation. Note that this is the result given by docking simulation by SYBYL, which means no molecular dynamics. The protein-structure keeps the shape in docking simulation, while the ligand only can change the shape. However, we could analyze some properties using the Sybyl package. Each Cl atom of DIF-1 and DIf-2 was contained in the receptor-binding site and affects the binding process as a necessary factor. Because of the lack of hydrophobic interaction of chlorine (when performing the docking simulation of DIF-3), the obtained docking model seems to form hydrogen bonds in order to maximize the score. As a result of steric effects, and the bonds of the DIF-3 molecule extending or shrinking to bind with the binding-pocket of the protein, the structure is unrealistically warped.

Figure 10. DIF-2 could bind to DhkM as well but showed a lower bioactivity than DIF-1 because of the weaker interaction and compatibility (Ligand drawn in sticks).

Figure 11. The molecule of DIF-3 is unrealistically warped in order to bind with the site. The structure would not be stable enough to compose a good binding model because of the lack of hydrophobic interaction (Ligand drawn in sticks).

Table 2. Summary of the results for the docking simulations of DIF-1 and DIF-2

On the other hand, the binding energies of DIF-2 and DIF-3 have been also calculated with MMPBSA as DIF-1. The energy was estimated as 0.5856 kcal/mol and 2.035 kcal/mol. As a result, we could remarkably see that DIF-1&DhkM (−11.09 kcal/mol) is a more favorable complex.

The reason for lower scores of DIF-2 and DIF-3 than DIF-1 is thought to be:

Furthermore, it has been reported that DIF-2 possessed 40% of the activity of DIF-1 and the monochloro-DIF-1 is much less active than DIF-1.²⁷ And our speculation and the result of docking simulations agree very well with the experimental data²⁷ as summarized in Table 3, which shows our binding model and proposed binding mechanism could be correct.

Table 3. Summary of the results for prediction, docking simulation and experimental data

The folding mechanism of DhkM has also been studied based on the binding results. It is very interesting that whether the differences between the structures with and without ligand binding would occur or not. In our work, two MD simulations were performed (Figure 12). One is the simulation for DhkM only without DIF-1, another one is the simulation for DhkM bound with DIF-1. These two results were compared, and they folded in different tertiary structures at 90 ns. It is considered that ligand binding could obviously affect protein folding process, with a totally different tertiary structure appearing, and finally producing an unknown function. By referring to Figure 5a, the tertiary structure seems to keep the shape up to 90 ns if docking with the ligand, DIF-1.

Figure 12. The result of the continued MD simulations for DhkM without DIF-1 (left) and DhkM bound with DIF-1 (right) folded in different tertiary structures at 90 ns.