Introducing a two‐dimensional graph of docking score difference vs. similarity of ligand‐receptor interactions

Observation of molecular docking results was generally performed by analyzing the docking score and the interacting amino acid residues separately either in tables or graphs. Sometimes it was not easy to rank the tested ligands’ docking results, especially if there were many ligands. This study aims to introduce a new way to analyze docking results with a two‐dimensional graph between the difference in docking score and the similarity of ligand‐receptor interactions. Molecular docking was performed with one reference ligand and several test ligands. The docking score difference was obtained between the test and the reference ligands as the graph’s x‐axis. Meanwhile, the y‐axis contains the similarity of ligand‐receptor interactions, obtained from the ratio of amino acid residues and the types of interactions between the test and reference ligands. Docking result analysis was more straightforward because two critical parameters were presented in one graph. This graph could be used to support the analysis of the docking results.


Introduction
Molecular docking (or simply docking) is an in silico method used to analyze the interactions between two molecules. Of the two molecules, one will act as a test compound or ligand, while the other will act as a target or known as a receptor. In its use, the docking method is widely used in various purposes in the field of drug design and discovery, especially for screening in the discovery of potential compounds with certain potential activities, as well as to explain the mechanism of action of the inter actions that occur between drug compounds with known activity against the target protein (Lin et al. 2020; Meng et al. 2012. Compared to several other in silico methods, molec ular docking is one of the most popular and widely used, both as the primary method and for confirming other meth ods. From the beginning of 2020 to October 2020 alone, there have been more than 10,000 articles from the Sco pus database published with titles, keywords, or abstracts containing the keyword "molecular docking." This num ber is more than 1,000 articles compared to the previous year, confirming molecular docking's popularity. This is mainly due to molecular docking's various advantages, in cluding a brief analysis process, low cost, and guaranteed safety because it does not use dangerous chemical reagents (Deshpande et al. 2020; Pinzi andRastelli 2019).
Aside from its popularity, molecular docking also presents various challenges in its analysis. These obstacles are generally related to the type of software used, consid ering that much software can be used to perform molec ular docking, both free and paid. Apart from the techni cal problems associated with the software used, one of the biggest challenges in analyzing the docking results is the ligand ranking based on the docking result (Pagadala et al. 2017).
Observations of molecular docking results are gen erally carried out by analyzing two main parameters of the docking results separately: the docking score and the interacting amino acid residues. These two parame ters are equally important in the analysis of docking re sults. If the docking score is often associated with lig and affinity for the receptor, the ligandreceptor interac tion becomes an indicator of whether the resulting interac tion can cause activity or not, compared to reference/co crystal ligands (Pantsar andPoso 2018; Ferreira et al. 2015). However, analyzing the two parameters simultane ously requires more effort, considering that one is quanti tative while the other is qualitative (Vieira andSousa 2019; Ramírez andCaballero 2016). This may not be a concern if only a few ligands are tested, but it will cause problems if the number of test ligands is large. The problem will arise when presenting the data, where there will be a very long table to present all the data.
An approach that can be taken to facilitate the anal ysis of the docking results' two parameters is to change the ligandreceptor interaction parameters into quantita tive parameters so that the two parameters can be plotted as a graph with two parameters. The rational way is to com pare the interacting amino acids and the types of interac tions in the test ligands with those shown by the reference ligand, then expressed in terms of the percentage of simi larity. This approach assumes that the more similar amino acid interactions will provide a greater chance of a similar mechanism of action, which will result in the same activity (Li et al. 2019; Ramírez andCaballero 2016). To the best of our knowledge, this approach has not been previously reported.
Based on this background, this study aims to intro duce a new way to analyze docking results with a two dimensional graph between the difference in docking score and the similarity of ligandreceptor interactions.

Materials and Methods
For demonstration purposes, the docking process was per formed using proteins that bind to the highly selective reference ligands, and the test ligands were compounds known to target these proteins. These reference ligands include a cocrystal ligand that binds to the receptor or compounds known to interact with the receptor's binding site (Kolb and Irwin 2009). The docking score and the ligandreceptor interactions will be used to compile a two dimensional graph between the difference in docking score and the similarity of ligandreceptor interactions.
The protein chosen was prostaglandin G/H synthase 2 or cyclooxygenase 2 (COX2), an enzyme as a target re ceptor with PDB ID 3LN1, which binds celecoxib, a selec tive COX2 inhibitor. The receptor was chosen because of the availability of cocrystal ligand that were selective to the receptor, a high enough crystal resolution (2.40 Å) with 0 Ramachandran outliers (Wang et al. 2010), the results of redocking with low RMSD and ΔG suggesting good bind ing affinity (Mandour et al. 2016), and has been reported to be used in docking studies more than 10 times in the past three years in the Scopus database. However, this ap proach was expected to be applied to other receptors as well.

Ligands preparation
As a representative of the test ligands, 15 anti inflammatory compounds were used either known to act on the COX2 pathway (e.g., etoricoxib) or not (ac etaminophen, aspirin). The test ligand structure was ob tained from PubChem (https://pubchem.ncbi.nlm.nih.gov /) and then downloaded in SDF format. All test ligands were then prepared according to the method reported by Pratama et al. (2020) then saved in .pdbqt format.

Receptor preparation
The receptor used in the demonstration is COX2, which binds to celecoxib (PDB ID 3LN1) and was downloaded from the Protein Data Bank website (https://www.rcsb.org /). The receptor consists of four chains (A, B, C, and D), with the chains used for the docking process was chain A (Wang et al. 2010). The parts of the receptors that were not used (e.g., water, solvent, unused chains) were then removed and given polar hydrogen as well as chargesand finally saved in .pdbqt format using AutoDockTools 1.5.6.

Hardware and software
The hardware and software used in this study were the same as the research reported by Pratama et al. (2020), with AutoDock Vina for docking and Discovery Studio Vi sualizer for visualization. Twodimensional graphical cre ation of ligandreceptor interactions was carried out using Microsoft Excel 2019.

Validation of docking protocol
The docking protocol validation was carried out by the redocking method reported by Morris et al. (2009). The observed parameter was a rootmeansquare deviation (RMSD). The RMSD value less than 2 Å indicating a valid docking protocol and can be used for the docking process.

Molecular docking
Docking for all test ligands performed in the same way as the validation process with similar sizes and positions of the grid box. The results were grouped under two param eters: free energy of binding (ΔG; kcal/mol) and ligand receptor interactions. The ligandreceptor interactions are recorded based on two parameters: the amino acids that interact and the types of interactions that occur. The dock ing process was repeated five times, and the average ΔG value and the deviation were determined. The maximum allowable deviation value was ±0.05 kcal/mol to avoid high variation. The two parameters were then compared their similarity to celecoxib as a reference ligand then the average was calculated and expressed as a percentage.

Twodimensional graph of ligandreceptor in teractions
The difference in ΔG values and the ligandreceptor interactions obtained earlier was then used to create a twodimensional graph. The xaxis was filled with the reduction in the ΔG value of each test ligand against the ΔG value of the reference ligand (celecoxib). The difference from the ΔG value of each test ligand against ΔG celecoxib was calculated based on the following Equation 1: The higher the difference in value, the lower the poten tial for affinity than the reference ligand, and vice versa. If the difference value is negative, it indicates that the test ligand has a higher potential for affinity than the reference ligand.
Meanwhile, the yaxis is filled with the average per centage of the similarity of the amino acids that interact and the types of interactions. The similarity of the interact ing amino acids was calculated by dividing the number of amino acids interacting with both the test and the reference ligand, divided by the total amino acids interacting with the reference ligand. Meanwhile, the similarity of interac tion types is calculated by dividing the number of amino acids with the same types of interactions for both the test ligand and the reference ligand, divided by the number of amino acid interactions on the reference ligand. The two parameters were then averaged and expressed as a percent age as shown in the following Equation 2: ×100% (2) %similarity = similarity of ligandreceptor interaction nAAtest = the number of amino acids of the test ligand that also interacts with the reference ligand nAAref = the number of amino acids of the reference ligand intAAtest = the number of amino acid interactions of the test ligand that is also present in the reference ligand intAAref = the number of amino acid interactions of the reference ligand The number 0.5 in the ligandreceptor interaction sim ilarity equation shows the impact of each parameter on the similarity of the ligandreceptor interaction. At the time of writing, it was still unknown how the influence of the amino acids that interact and the types of interactions that occur on ligandreceptor interactions. The initial assump tion used was that the two parameters had the same effect on the ligandreceptor interaction, so the impact weight given was the same. This opens up opportunities for fur ther research regarding the comparison of the effects of the two on ligandreceptor interactions. The two parame ters were then used to quantitatively express the similarity of the ligandreceptor interactions on the yaxis, as shown in Figure 1.

Validation of docking protocol
The RMSD value of the redocking process obtained was 0.857 Å, indicating that the docking protocol was valid. The visualization of ligand overlays from redocking with reference ligand from crystallographic results was pre sented in Figure 2. The redocking ligands' appearance shows the same orientation as the crystallographic ligands, apart from a slight shifted in position. The results of the FIGURE 2 Overlays of redocking ligands (blue) with reference ligands from crystallography data (green) at receptors 3LN1 with RMSD 0.857 Å.  Table 1. The dimensions of the grid box used were relatively small with a size of 32 x 20 x 26 Å, ad justing to the size of the reference ligands, which are also small. The redocking results showed that 24 amino acids interacted with celecoxib. The interactions were dom inated by weak interactions such as Van der Waals (11 amino acids) but had considerably strong interactions such as hydrogen bonds (six amino acids).

Molecular docking
Docking of all test ligands shows the results as predicted, in which the ligands that were known to target COX2, such as etoricoxib, had the lowest AG value of 11.4 ± 0 kcal/mol. This value is only a 1.1 kcal/mol difference from that shown by the celecoxib from the redocking process which is 12.5 ± 0 kcal/mol, while the other test ligands show a difference between 2.8 to 6.24 kcal/mol. Etori coxib also had a high similarity of interacting amino acid residues with 79.17%, but interaction type was lowered with 25%. The average of both was 52.08%, still low ered than some other test ligands such as indomethacin and nimesulide, with 68.75% and 56.25%, respectively. Apart from having the highest similarity, indomethacin also had the secondlowest docking value after etoricoxib with 9.7 ± 0 kcal/mol, making both etoricoxib and indomethacin the two test ligands with the closest parameter of celecoxib as reference ligands. In other words, these two ligands were the strongest candidates compared to other test lig ands. The docking of all test ligands at the binding site of the 3LN1 receptor was presented in Table 2 and then used to filled in the twodimensional graph plotted previously, as presented in Figure 3.

Discussion
The docking results were consistent with other studies re ported previously by Sadasivam et al. (2020), in which a ligand such as etoricoxib, designed as a selective COX2 inhibitor, show the highest potency compares to other lig ands. The docking results also show that indomethacin has potential, although it is not a selective COX2 inhibitor. These results are consistent with research from Abuelizz et al. (2017) and Oniga et al. (2017), who reported in domethacin's potency against COX2. The 3LN1 recep tor selection itself was based on the consistency of results from several previous studies (Coybarrera 2020; Molinari et al. 2019; Shrivastava et al. 2017. In other words, the docking results obtained for the demonstration are appro priate and in line with other studies. However, there are more interesting points to discuss.
As mentioned in the previous section, the more the number of test ligands will cause more data to be pre sented. However, by presenting the docking result data in a twodimensional graph, the data presented is much more concise while still presenting essential data for anal ysis. Readers do not have to read the entire docking re sults table to find the ligands with the best docking results. Otherwise, they only need to look for the test ligand lo cated on the upper left of the graph to determine which test ligand has the smallest difference in docking values and the highest ligandreceptor interaction similarity with the reference ligand. This also makes the presentation of the results shorter even though the number of ligands tested is immense, considering that the graph space remains the same and does not depend on the number of test ligands.
There are times when the ligands with the smallest docking score difference and the highest ligandreceptor interaction similarity compared to reference ligands are two different test ligands, as presented in the results of this study (etoricoxib and indomethacin). In that case, the recommendations that can be given for both, consid ering these two parameters are equally crucial in molec ular docking. It has been occasionally reported that the test ligand with the best docking score has an unstable in teraction with the receptor when tested by molecular dy namics, compared to the test ligand with a ligandreceptor interaction closer to the ligand reference (Lam et al. 2018; Salmaso andMoro 2018).
One of the essential points in presenting the docking results in this twodimensional graph is the obligation to provide reference ligands for comparison, considering that some of the receptors available in the protein data bank are aporeceptors (Forli et al. 2016). For new target receptors for which no proven reference ligands have been found, it is not possible to present the docking results in a two dimensional graph. Under these conditions, it is highly recommended to conduct in vitro and in vivo studies first to ensure that the reference ligand to be used is proven.
This twodimensional graph has been used before to present the docking results of several types of receptors (Pratama et al. 2020(Pratama et al. , 2021. Research by Pratama et al. (2020) even presents the docking results of two receptors ΔG±SD (kcal/mol) -6.26 ± 0.05 -6.4 ± 0 -8.5 ± 0 -8.72 ± 0.04 -11.4 ± 0 -7.7 ± 0 -9.7 ± 0 -9.44 ± 0.05 -9.08 ± 0.04 -8.92 ± 0.04 -8.8 ± 0 -8.3 ± 0 -9.1 ± 0 -8.1 ± 0 -9.  in one graph simultaneously by providing different color dots, showing the advantage that this twodimensional graph can be used to present the docking results of multi ple receptors at once. This is advantageous because it can streamline research reports, considering that some studies with molecular docking are carried out at more than one receptor, while some scientific journals impose limits on the number of supporting illustrations such as tables and figures in the manuscript. Thus, it is expected that the pre sentation of docking results in this twodimensional graph can become a new trend in the presentation of docking re sults. Apart from its simplicity, two things can still be im proved in the presentation with these twodimensional graphics. First, given that some of the test ligands can have very similar docking scores and ligandreceptor in teractions often due to similar structures (Gimeno et al. 2019), it is possible that some ligands are very close to gether and even overlap on the graph. An enlarged version of the graph is required to overcome the above mentioned in the area where the ligands overlap. Second, as previ ously described, the calculation between the similarity of the interacting amino acids and the types of interactions averaged is assumed to have the same impact. Meanwhile, no previous research has proven and calculated these two parameters' impact on the docking results. Hence, there is an opportunity for further research on how to compare the impact of the two parameters on docking results, whether they have the same impact or whether one has more impact than the other.

Conclusions
The twodimensional graph between the difference in docking score and the similarity of ligandreceptor inter actions can support the analysis of the docking results by presenting the docking results briefly.