Theoretical Studies on Selectivity of HPK1/JAK1 Inhibitors by Molecular Dynamics Simulations and Free Energy Calculations

Hematopoietic progenitor kinase 1 (HPK1) is a negative regulator of T cell receptor, which has been regarded as a potential target for immunotherapy. Yu et al. observed the off-target effect of the high-throughput screening HPK1 kinase inhibitor hits on JAK1 kinase. The off-target effect is usually due to the lack of specificity of the drug, resulting in toxic side effects. Therefore, exploring the mechanisms to selectively inhibit HPK1 is critical for developing effective and safe inhibitors. In this study, two indazole compounds as HPK1 inhibitors with different selectivity towards JAK1 were used to investigate the selectivity mechanism using multiple computational methods, including conventional molecular dynamics simulations, binding free energy calculations and umbrella sampling simulations. The results indicate that the salt bridge between the inhibitor and residue Asp101 of HPK1 favors their selectivity towards HPK1 over JAK1. Information obtained from this study can be used to discover and design more potent and selective HPK1 inhibitors for immunotherapy.


Introduction
Hematopoietic progenitor kinase (HPK1), also known as MAP4K1, is a serine/threonine kinase that is a member of the mitogen-activated kinase kinase kinase kinase (MAP4K) family [1]. HPK1 is predomainantly expressed in hematopoietic cells and mainly regulates T cell activation through the TCR signaling pathway [2,3]. After TCR activation, HPK1 is phosphorylated on tryrosine 379 by kinases such as lymphocyte-specific protein tyrosine kinase (LCK), the activated HPK1 feedback phosphorylates Ser376 of Src homology 2 domain-containing leukocyte protein of 76 kDa (SLP76), which binds to negative regulator complex 14-3-3 [4][5][6][7]. This results in ubiquitination of lysine 30 on SLP76 and degration of SLP76 and subsequently limits strength and duration of TCR signals [8]. HPK1 kinase dead mice and HPK1 knockout mice models showed that loss of HPK1 kinase function enhanced T cell proliferation in response to TCR stimulation and reduced tumor growth and did not show fatal inflammation [3,[9][10][11][12]. These studies have demonstrated that HPK1 inhibitors have potential applications in improving T cell function and combating immunosuppressive tumor microenvironments. Furthermore, the combination of HPK1 inhibitor and anti-programmed death-ligand-1 (anti-PD-L1) performed better anti-tumor immune effects than blocking individual targets [9,13,14]. Therefore, HPK1 inhibitors have potential roles in cancer immunotherapy.
Due to the structural similarities of the ATP binding site, kinase inhibitors typically display high promiscuity. Small inhibitors targeting specific kinase can be off-target in kinases of different families. Yu et al. observed off-target effects of the highpthroughput screening HPK1 kinase inhibitor hits on JAK1 kinase [15]. JAK1 is a member of the Janus kinase (JAK) family and plays a key role in cytokine-mediated inflammatory and autoimmune responses through the JAK/STAT signal pathway [16,17]. Off-target effect is usually due to the lack of specificity of the drug, resulting in toxic side effects. Therefore, autoimmune responses through the JAK/STAT signal pathway [16,17]. Off-target effe usually due to the lack of specificity of the drug, resulting in toxic side effects. There designing potent HPK1 inhibitors with selectivity is very important. Yu et al. report series of indazole inhibitors as HPK1 inhibitors with different degrees of selectivity wards JAK1, such as XHS (751 fold selectivity over JAK1) and XHV (112 fold) [15]. selective inhibitors may provide us with information to reveal the mechanism of th lectivity of ligand for HPK1 over JAK1.
Molecular dynamics simulations have emerged as a powerful approach to exp mechanisms of inhibitor selectivity [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. In the present study, two indazole inhibi XHS and XHV, were selected to study their selective mechanisms towards HPK1 JAK1. XHS and XHV are highly (751-fold; HPK1 IC50 = 2.6 nM, JAK1 IC50 = 1952.6 nM) moderately (112-fold; HPK1 IC50 = 89 nM, JAK1 IC50 = 9968 nM) selective inhibitors, res tively, for HPK1 over JAK1. The chemical structures and IC50 of XHS and XHV are sh in Figure 1. We studied the selectivity of XHS and XHV by multiple molecular comp tional methods. MM-GBSA was used to determine the binding free energy. The umb sampling simulations were used to study the dissociation process of the two inhib from the HPK1/JAK1 binding pocket. The results demonstrated that the salt bridge in action between the residue Asp101 of HPK1 and the inhibitor seems to be importan the selective inhibition of HPK1. Furthermore, the umbrella sampling results suggest the dissociation process is important in kinase selectivity. Consequently, we suggest in addition to focusing on binding affinity, the dissociation process should also be ta into consideration.

Sequence Alignment of HPK1 and JAK1
To compare the difference between HPK1 and JAK1 from sequence, the Clusta pairwise sequence alignment program tools for aligning two sequences in a pairwise m ner, and the resulting sequence alignment was displayed by the online tool ESPript3.0 in Figure 2. In it, the white characters with red background indicated that the amino sequences are completely identical, and the red characters with white background i cated that the amino acids are different but similar in properties. HPK1 shared 23.6% quence identity and 53.5% sequence similarity with JAK1 in kinases domain and sh 46.5% sequence identity and 76.7% sequence similarity with JAK1 at ATP binding po

Sequence Alignment of HPK1 and JAK1
To compare the difference between HPK1 and JAK1 from sequence, the ClustalW2 pairwise sequence alignment program tools for aligning two sequences in a pairwise manner, and the resulting sequence alignment was displayed by the online tool ESPript3.0 [32] in Figure 2. In it, the white characters with red background indicated that the amino acid sequences are completely identical, and the red characters with white background indicated that the amino acids are different but similar in properties. HPK1 shared 23.6% sequence identity and 53.5% sequence similarity with JAK1 in kinases domain and shared 46.5% sequence identity and 76.7% sequence similarity with JAK1 at ATP binding pocket.

Figure 2.
Sequence alignment of HPK1 and JAK1. The white characters with red background indicated that the amino acid sequences are completely identical, and the red characters with white background indicated that the amino acids are different but similar in properties.

The Conformational Stabilities of the Four Systems
A detailed analysis of root mean square deviation (RMSD) showed that the backbone atoms, active pocket atoms and ligand heavy atoms of the three systems HPK1-XHS, HPK1-XHV, JAK1-XHV tend to converge after 50 ns, indicating that the systems became sufficiently stable through 200 ns of simulations ( Figure 3A,C,D). The protein backbone of the JAK1-XHS system is not as stable as the other three systems, but the active pocket atoms and ligand heavy atoms are stable throughout the simulations ( Figure 3B). In general, through 200 ns conventional molecular dynamics simulations, we obtained more reasonable binding conformation for the two aligned systems JAK1-XHS and JAK1-XHV.

The Conformational Stabilities of the Four Systems
A detailed analysis of root mean square deviation (RMSD) showed that the backbone atoms, active pocket atoms and ligand heavy atoms of the three systems HPK1-XHS, HPK1-XHV, JAK1-XHV tend to converge after 50 ns, indicating that the systems became sufficiently stable through 200 ns of simulations ( Figure 3A,C,D). The protein backbone of the JAK1-XHS system is not as stable as the other three systems, but the active pocket atoms and ligand heavy atoms are stable throughout the simulations ( Figure 3B). In general, through 200 ns conventional molecular dynamics simulations, we obtained more reasonable binding conformation for the two aligned systems JAK1-XHS and JAK1-XHV.

Binding Free Energy Calculations
To compare the binding affinities of two inhibitors towards HPK1 and JAK1, the binding free energy (∆G bind ), which includes gas-phase energy, solvation energy, and entropy term, was computed using MM-GBSA method. Table 1 summarized ∆G bind for four systems and experimental binding free energies calculated from IC 50 values. Results showed that the calculated ∆G bind of XHS (−26.60 ± 7.08 kcal/mol) and XHV (−18.68 ± 6.92 kcal/mol) binding to HPK1, which corresponds to the order of experimental binding affinities (HPK1: XHS, IC 50 = 2.6 nM; XHV, IC 50 = 89 nM). The calculated ∆G bind of inhibitors XHS (−18.09 ± 8.87 kcal/mol) and XHV (−10.48 ± 4.98 kcal/mol) binding to JAK1 are also in agreement with the experimental binding affinity order (JAK1: XHS, IC 50 = 1952.6 nM; XHV, IC 50 = 9968 nM). In addition, HPK1-selective inhibitor XHS has a lower calculated binding energy for HPK1 (∆G bind = −26.60 ± 7.08 kcal/mol) than JAK1 (∆G bind = −18.09 ± 8.87 kcal/mol) and shows a more significant difference in ∆G bind values for HPK1 over JAK1 than XHV. In general, the difference in ∆G bind reflects that XHS is a more potent selective HPK1 inhibitor over JAK1 than XHV.
To further study the driving force for selective binding of two inhibitors on HPK1 over JAK1, the binding enthalpy change (ΔH) was decomposed into independent binding free energy (Table 1) by MM-GBSA method. The results revealed that the primary driving force between inhibitors and kinases (HPK1 and JAK1) was the non-polar interaction energy, which is the sum of ΔEvdW and ΔGnonpl,sol, while the polar interaction (the sum of ΔGele and ΔGGB) had a negative impact on binding. The ΔGbind difference between inhibitor XHS and two kinases is polar interaction (9.52 ± 1.43 kcal/mol for HPK1 kinase, 17.47 ± 0.71 kcal/mol for JAK1 kinase). The ΔGbind difference between inhibitor XHV and two kinases is non-polar interaction (−56.48 ± 2.90 kcal/mol for HPK1 kinase, −49.16 ± 3.99 kcal/mol for  To further study the driving force for selective binding of two inhibitors on HPK1 over JAK1, the binding enthalpy change (∆H) was decomposed into independent binding free energy (Table 1) by MM-GBSA method. The results revealed that the primary driving force between inhibitors and kinases (HPK1 and JAK1) was the non-polar interaction energy, which is the sum of ∆E vdW and ∆G nonpl,sol , while the polar interaction (the sum of ∆G ele and ∆G GB ) had a negative impact on binding. The ∆G bind difference between inhibitor XHS and two kinases is polar interaction (9.52 ± 1.43 kcal/mol for HPK1 kinase, 17.47 ± 0.71 kcal/mol for JAK1 kinase). The ∆G bind difference between inhibitor XHV and two kinases is non-polar interaction (−56.48 ± 2.90 kcal/mol for HPK1 kinase, −49.16 ± 3.99 kcal/mol for JAK1 kinase).

Hydrogen Bond and Salt Bridge Analysis
Hydrogen bond (H-bond) occupancy was calculated to track the formation of hydrogen bonds between the inhibitors and their corresponding targets during the MD simulations. According to Table 2, in HPK1-XHS system, the stable H-bond is formed between backbone NH of Cys94 and the N atom of indazole scaffold (90.80% occupancy). In JAK1-XHS system, the stable H-bond is formed between backbone NH of Leu959 and the N atom of indazole scaffold (86.64% occupancy). In HPK1-XHV system, the stable H-bond is formed between backbone NH of Cys94 and the N atom of indazole scaffold (82.36% occupancy). In JAK1-XHV system, the stable H-bond is formed between backbone NH of Leu959 and the N atom of indazole scaffold (65.76% occupancy). The H-bond occupancy of HPK1 systems is higher than that of JAK1 systems, which can explain the fact that XHS and XHV have better inhibition activity for HPK1 over JAK1. The salt bridge is double-charged hydrogen bond based on the classification of hydrogen bond proposed by Gilli et al. [33][34][35][36]. Next the salt bridge occupancy was calculated to monitor the formation of salt bridge between the inhibitors and their corresponding targets. As shown in Table 3, salt bridges were detected in only two systems, and the salt bridge between the protonated NH+ of piperazine and the side chain carboxyl of Asp101 (OD with 62.03% occupancy) in the HPK1-XHS system is more stable than the salt bridge between protonated NH+ of piperazine with the side chain carboxyl of Glu966 (OE with 7.41% occupancy) in the JAK1-XHS system. In general, the salt bridge between inhibitor and Asp101 is critical for enhanced potency and selectivity.

The Binding Mode between Kinases (HPK1/JAK1) and Inhibitors (XHS/XHV)
In order to study the detailed binding mode between the two inhibitors (XHS/XHV) and the two kinases (HPK1/JAK1) and explain the difference in binding affinity, the last frame of the four system trajectories were extracted. There are similar binding pockets between HPK1 and JAK1. The binding modes of inhibitors with kinases (HPK1 and JAK1) are depicted in Figure 5. According to Figure 5A, the structure of XHS bound to HPK1 demonstrated (1) that the residues Val31, Lys46, Ala141, Asn142 and Ala154 surround 2-fluoro-6-methoxyphenyl group and more lipophilic fluorine group point to P-loop of protein; (2) that the residues Ala44, Val75, Met91, Glu92, Phe93, Cys94 and Leu144 surround indazole core, and the backbone NH of Cys94 forms a H-bond with the N of the indazole; (3) that the residues Leu23, Gly24, Gly97, Ser98 and Asp101 surround meta piperazine substituted pyridine group, and salt bridge is formed between Asp101 and protonated NH+ on piperazine. According to Figure 5B, the structure of XHS bound to JAK1 demonstrated (1) that the residues Val859, Lys908, Arg1007, Asn1008, Gly1020 and Asp1021 surround 2-fluoro-6-methoxyphenyl group and the more lipophilic fluorine group point to P-loop of protein; (2) that the residues Ala906, Val938, Met956, Glu957, Phe958, Leu959 and Leu1010 surround indazole core and the backbone NH of Leu959 forms a H-bond with the N of the indazole core; and (3) that the residues Leu881, Gly882, Gly962, Ser963 and Glu966 surround meta piperazine substituted pyridine group and an unstable salt bridge is formed between Glu966 and protonated NH+ on piperazine. As depicted in Figure 5C, the structure of XHV bound to HPK1 demonstrated that (1) the residues Val31, Lys46, Ala141, Asn142, Ala154 and Asp155 surround 2-fluoro-6-methoxyphenyl group, the more lipophilic fluorine group point to P-loop of protein; (2) that the residues Ala44, Val75, Met91, Glu92, Phe93, Cys94 and Leu144 surround indazole core and the backbone NH of Cys94 forms a H-bond with the N of the indazole; and (3) that the residues Leu23, Ala96, Gly97, Ser98 and Asp101 surround para piperazine substituted phenyl group, which is located at the entrance of the binding site. As depicted in Figure 5D, the structure of XHV bound to JAK1 demonstrated (1) that the residues Val889, Lys908, Arg1007, Asn1008, Gly1020 and Asp1021 surround 2-fluoro-6-methoxyphenyl group and the more lipophilic fluorine group point to P-loop of protein; (2) that the residues Ala906, Val938, Met956, Glu957, Phe958, Leu959 and Leu1010 surround indazole core and the backbone NH of Leu959 forms a H-bond with the N of the indazole core; and (3) that the residues Leu881, Ser961, Gly962, Ser963 and Glu966 surround para piperazine substituted phenyl group, which is located at the entrance of the binding site. According to the results of our simulation, there are similar binding sites between HPK1 and JAK1. The comparative results of the binding sites suggest that the salt-bridge interaction between Asp101 of HPK1 kinase and the inhibitor plays a dominant role in the binding of HPK1 selective inhibitors. The inhibitor XHS is shown as green sticks, the inhibitor XHV is shown as yellow sticks and the key residues of HPK1 and JAK1 kinase binding to inhibitors are shown as cyan sticks. Hydrogen bonds and salt bridges kinases and inhibitors are indicated by red and black dashes, respectively.

Unbinding Pathways of Inhibitors (XHS and XHV) Dissociating from Kinases (HPK1 and JAK1)
In this study, the last frame of the four system trajectories were used as the initial structures for the umbrella sampling simulations. As depicted in Figure 6, each system reaches convergence of PMFs after performing 7 ns umbrella sampling simulations. As depicted in Figures 7 and 8, the PMF depth derived from the umbrella sampling simulations of the HPK1 system is higher than that of the JAK1 system. To further determine the detailed process of dissociation, the last 1 ns in each system was used for analysis. As indicated in Figure 7F, minimum value of PMF (0 kcal/mol) was calibrated to 0 Å on RC to indicate the binding state of XHS in HPK1/JAK1. The RCs were extended from 0 to 20 Å represent XHS 20 Å away from the RCs. XHS dissociated from RC with 22.43 ± 0.07 kcal/mol for HPK1 and 15.09 ± 0.86 kcal/mol for JAK1 (Table 4); the difference is 7.34 The inhibitor XHS is shown as green sticks, the inhibitor XHV is shown as yellow sticks and the key residues of HPK1 and JAK1 kinase binding to inhibitors are shown as cyan sticks. Hydrogen bonds and salt bridges kinases and inhibitors are indicated by red and black dashes, respectively.

Unbinding Pathways of Inhibitors (XHS and XHV) Dissociating from Kinases (HPK1 and JAK1)
In this study, the last frame of the four system trajectories were used as the initial structures for the umbrella sampling simulations. As depicted in Figure 6, each system reaches convergence of PMFs after performing 7 ns umbrella sampling simulations. As depicted in Figures 7 and 8, the PMF depth derived from the umbrella sampling simulations of the HPK1 system is higher than that of the JAK1 system. To further determine the detailed process of dissociation, the last 1 ns in each system was used for analysis. As indicated in Figure 7F, minimum value of PMF (0 kcal/mol) was calibrated to 0 Å on RC to indicate the binding state of XHS in HPK1/JAK1. The RCs were extended from 0 to 20 Å represent XHS 20 Å away from the RCs. XHS dissociated from RC with 22.43 ± 0.07 kcal/mol for HPK1 and 15.09 ± 0.86 kcal/mol for JAK1 (Table 4); the difference is 7.34 kcal/mol. This finding suggested that the inhibitor XHS dissociated from JAK1 much easier than that of XHS from HPK1 during the dissociation process. When XHS escapes from the binding pocket of HPK1 at point B of the PMF curve (~1.6 Å, Figure 7F), an obvious energy barrier is first displayed, where the 2-fluoro-6-methoxyphenyl group moves outward along the RC direction and the H-bond between the N of the indazole core and the NH of Cys94 is broken. Afterward, the piperazine moiety of XHS flips outward and the salt bridge with Asp101 is broken (ascending region B-C in Figure 7F). Next, the meta piperazine substituted pyridine group is exposed to the solvent (point D in Figure 7F). Finally, XHS completely dissociated from HPK1 (point E in Figure 7F). For the JAK1-XHS system, the PMF curve rises significantly in region A'-B' (Figure 7F'), where the 2-fluoro-6-methoxyphenyl group moves outward along the RC direction the H-bond between the N of indazole core and the NH of Leu959 is broken. Then the meta piperazine substituted pyridine group flipped and exposed to solvent (point C', Figure 7F Figure 7F), an obvious energy barrier is first displayed, where the 2-fluoro-6-methoxyphenyl group moves outward along the RC direction and the H-bond between the N of the indazole core and the NH of Cys94 is broken. Afterward, the piperazine moiety of XHS flips outward and the salt bridge with Asp101 is broken (ascending region B-C in Figure 7F). Next, the meta piperazine substituted pyridine group is exposed to the solvent (point D in Figure 7F). Finally, XHS completely dissociated from HPK1 (point E in Figure 7F). For the JAK1-XHS system, the PMF curve rises significantly in region A'-B' (Figure 7F'), where the 2-fluoro-6-methoxyphenyl group moves outward along the RC direction the H-bond between the N of indazole core and the NH of Leu959 is broken. Then the meta piperazine substituted pyridine group flipped and exposed to solvent (point C', Figure 7F     The PMF curve of inhibitor XHV exhibits a completely different behavior compared to the PMF curve of inhibitor XHS. According to Table 4, XHV dissociated from RC with 17.68 ± 0.16 kcal/mol for HPK1 and 10.02 ± 0.29 kcal/mol for JAK1; the difference is 7.66 kcal/mol. For the HPK1-XHV system, the PMF curve first rises significantly (regions A-B-C, Figure 8F); the piperazine moiety of the inhibitor XHV is flipped, and the H-bond between the N of indazole core and the NH of Cys94 is broken. Then the piperazine moiety is released from the binding pocket, and the inhibitor XHV eventually dissociated from HPK1 (regions D-E in Figure 8F). Distinct from that of HPK1, two stages of sharp rise in PMF (regions A'-B'-C' and D'-E' in Figure 8F') are observed before the release of the inhibitor XHV from the binding pocket of JAK1. With the increase in PMF, the 2-fluoro-6methoxyphenyl group of inhibitor XHV moves to entrance of the binding pocket, and then the H-bond between indazole core and Leu959 is broken.(regions A'-B'-C' in Figure 8F'). The indazole core is then flipped and exposed to the solvent, and then the inhibitor XHV is completely dissociated from the binding pocket of JAK1 (region D'-E' in Figure 8F'). The difference in PMF during dissociation of inhibitor XHV from HPK1 and JAK1 is due to the difference in the hydrophobic interaction of the 2-fluoro-6-methoxyphenyl group and piperazine moiety of XHV with kinases. In summary, the dynamic and energetic information obtained in this research is helpful for understanding binding/unbinding mechanism of the inhibitors XHS/XHV to and from HPK1/JAK1.

Sequence Alignment and Structure Superposition
The crystal structure of human HPK1 kinase with each of ligands XHS and XHV (PDB ID: 7L25 and 7L24) [15] and JAK1 kinase bound with inhibitor (PDB ID: 6SMB) [37]   The PMF curve of inhibitor XHV exhibits a completely different behavior compared to the PMF curve of inhibitor XHS. According to Table 4, XHV dissociated from RC with 17.68 ± 0.16 kcal/mol for HPK1 and 10.02 ± 0.29 kcal/mol for JAK1; the difference is 7.66 kcal/mol. For the HPK1-XHV system, the PMF curve first rises significantly (regions A-B-C, Figure 8F); the piperazine moiety of the inhibitor XHV is flipped, and the H-bond between the N of indazole core and the NH of Cys94 is broken. Then the piperazine moiety is released from the binding pocket, and the inhibitor XHV eventually dissociated from HPK1 (regions D-E in Figure 8F). Distinct from that of HPK1, two stages of sharp rise in PMF (regions A'-B'-C' and D'-E' in Figure 8F') are observed before the release of the inhibitor XHV from the binding pocket of JAK1. With the increase in PMF, the 2-fluoro-6methoxyphenyl group of inhibitor XHV moves to entrance of the binding pocket, and then the H-bond between indazole core and Leu959 is broken.(regions A'-B'-C' in Figure 8F'). The indazole core is then flipped and exposed to the solvent, and then the inhibitor XHV is completely dissociated from the binding pocket of JAK1 (region D'-E' in Figure 8F'). The difference in PMF during dissociation of inhibitor XHV from HPK1 and JAK1 is due to the difference in the hydrophobic interaction of the 2-fluoro-6-methoxyphenyl group and piperazine moiety of XHV with kinases. In summary, the dynamic and energetic information obtained in this research is helpful for understanding binding/unbinding mechanism of the inhibitors XHS/XHV to and from HPK1/JAK1.

Sequence Alignment and Structure Superposition
The crystal structure of human HPK1 kinase with each of ligands XHS and XHV (PDB ID: 7L25 and 7L24) [15] and JAK1 kinase bound with inhibitor (PDB ID: 6SMB) [37] were downloaded from the Protein Data Bank (https://www.rcsb.org/, accessed on 12 April 2021). The ClustalW2 [38] pairwise sequence alignment tools were used to align the human HPK1 (PDB: 7L25) and JAK1 (PDB: 6SMB) in FASTA format to determine the proportion of conserved and unique residues. To construct complex structures of 6SMB and inhibitors XHS and XHV, models 7L25, 7L24 and 6SMB were structurally aligned using PyMOL Molecular Graphics System, and the inhibitors XHS and XHV were merged into the binding pocket of JAK1. The crystal structure of HPK1-XHS, HPK1-XHV and aligned structure JAK1-XHS, JAK1-XHV were used as the initial structure for the following calculations.

Molecular Dynamics Simulation
The complexes structure of HPK1-XHS (crystal structure), JAK1-XHS, HPK1-XHV (crystal structure) and JAK1-XHV were used as starting coordinates for the MD simulations. The Amber18 package [39] was used to perform all of the MD simulations. The restrained electrostatic potential (RESP), calculated by Hartree-Fork (HF) method with 6-31G* basis set in the Gaussian09 [40], was used to fit the partial charges for the inhibitor XHS and XHV. The kinase parameters were described using the ff14SB force field, and the inhibitors were parameterized using the general AMBER force field (gaff). TIP3PBOX water molecules were used to solvate the complex 12 Å away to the boundary. The systems were minimized, heated and equilibrated after neutralizing the systems with chloride ions and sodium ions. At 300 K and 1.0 atmosphere pressure, 200 ns MD simulations were performed in NPT ensemble. The trajectory analysis was carried out using AMBER18 s cpptraj module.

Binding Free Energy Calculations
A total of 1000 snapshots extracted from each system's last 50 ns trajectories were used to calculate binding free energy using the MM-GBSA [41][42][43] method. The calculation for the binding free energy is as follows.
In Equation (1), G complex , G receptor , G ligand represent the free energy of complex, receptor and ligand molecules, respectively. In Equation (2), E gas and G sol represent the gas-phase energy and the solvation free energy. T∆S represents the change of conformational entropy when the ligand binds at temperature T. In Equation (3), E gas consists of internal energy (E int ), electrostatic energy (E ele ) and van der Waals energy (E vdW ). In Equation (4), G sol consists of polar solvation energy (G GB ) and nonpolar solvation energy (G nonpl , sol ). In Equation (5), G nonpl , sol was calculated by the SASA using a water probe radius of 1.4 Å. The surface tension constant γ was determined to be 0.0072 kcal/(mol·Å2) [44]. Due to the long computational time, the entropy term (−T∆S) was calculated using only 50 of the 1000 snapshots.

Free Energy Decomposition Analysis
In order to gain detailed information of the protein-ligand binding, MM-GBSA was utilized to decompose the interaction energies to each residue without consideration the contribution of entropy term. The essential residues of the indazole ihibitors for HPK1 over JAK1 binding are identified by the individual residue contributions derived by MM-GBSA energy decomposition, which also help us design more selective HPK1 inhibitors.

Umbrella Sampling Simulations
To characterize the unbinding process of each inhibitor from its corresponding targets, umbrella sampling simulations [45] were carried out to elucidate the drug binding selectivity using the pmemd program in AMBER18. The last frame of the four system trajectories were used as the initial structures for the umbrella sampling simulations. The distance between the carbon atom of inhibitors (red atom in Figure 1) and CA atom of Phe at hinge (HPK1 is Phe93, JAK1 is Phe958) was selected as reaction coordinates (RCs) using the tunnel prediction tool CAVER3.0 [46][47][48][49]. The RC was divided into 41 windows with a step of 0.5 Å. To guarantee convergence for each system, 7 ns MD simulations were run for each window. The residues at the end of each system are subjected to a force constant of 5 kcal mol −1 Å −2 to avoid the drift of receptor-ligand complexes. The umbrella sampling force constant was set to 5 kcal mol −1 Å −2 . The weighted histogram analysis method (WHAM) [50,51] was used to determine the potential of mean force (PMF) values.

Conclusions
In order to investigate the selective inhibition mechanism for HPK1 inhibitors over JAK1, two indazole inhibitors' binding selectivity was clarified using a computational approach that include all-atom conventional MD simulation, MM-GBSA free energy calculation and umbrella sampling simulations. The difference in HPK1 inhibitor activity between the two inhibitors is mainly due to the fact that the para-piperazine in ligand XHV extends to the solvent front, while the piperazine in ligand XHS changes direction to interact with Asp101, and the binding preference of the inhibitor XHS to HPK1 over JAK1 is controlled by the salt-bridge interaction of the piperazine ring with the carboxyl group on Asp101. The results of this study will offer useful recommendations for the rational development of HPK1 inhibitors with better selectivity over JAK1.
Author Contributions: H.G. and X.Y. designed the research. H.G., C.T. and Y.P. performed the research and analyzed the data. H.G. and X.Y. prepared the manuscript. All authors have read and agreed to the published version of the manuscript.