Computational insight into the mechanisms of action and selectivity of Afraxis PAK inhibitors
Di Han1, Huiqun Wang2, Wei Cui1, Beibei Zhang1 & Bo-Zhen Chen*,1
1 School of Chemical Sciences, University of Chinese Academy of Sciences, No. 19A, YuQuan Road, Beijing 100049, PR China
2 Department of Medicinal Chemistry, School of Pharmacy, Virginia Commonwealth University, 800 E Leigh Street, Richmond, VA 23298, USA
*Author for correspondence: [email protected]

Aim: The p21-activated kinases (PAKs) are involved in many important biological activity regulations. FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 were identified as PAKs inhibitors. Their detailed inhibitory mechanisms deserve further investigation. Results: Molecular dynamics simulations and further calculations for the PAK1/inhibitor and PAK4/inhibitor complexes indicate that their binding free ener- gies are basically consistent with the trend of experimental activity data. Conclusion: The anchoring of residues Leu347PAK1 and Leu398PAK4 is the structural basis for designing Afraxis PAK inhibitors. This study discloses the inhibitory mechanisms of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 toward PAK1 and PAK4 and some clues to enhance kinase activities and selectivities, which will provide valuable infor- mation to the development of more potent and selective PAK inhibitors.

First draft submitted: 24 September 2019; Accepted for publication: 5 December 2019; Published
online: 17 February 2020
Keywords: FRAX1036 • FRAX597 • G-5555 • molecular dynamics simulation • p21-activated kinase
The p21-activated kinases (PAKs) are a kind of serine/threonine kinases [1], which act as the effectors of Rac and Cdc42 belonging to the Rho family of small G-proteins [2–4]. PAKs have been found to participate in cancer cell signal transduction pathways and regulate the cell functions, which makes them become potential therapeutic targets in cancer [4–6]. Several studies have discovered that phosphorylating their substrates [7,8] can make PAKs be involved in many physiological processes, such as cell survival, cell proliferation, angiogenesis, and so on [9–12].
The PAK family is comprised of six isoforms (PAK1–PAK6), which can be divided into two groups: group I (PAK1–PAK3) and group II (PAK4–PAK6) according to their sequence and structural homology [13]. The sequence similarity in the kinase domain between group I and group II PAK isoforms is approximately 50%, while more than 90% among members for group I and about 80% for group II [13]. As the research on PAKs continues, increasing evidences suggest that although some effectors of the two groups are similar, the two groups of PAKs may recognize some distinct substrates and have unique biological functions [5].
Currently, owing to the flexibility of the ATP binding domain of PAKs, the number of potential selective inhibitors to PAKs that have been reported is limited [13]. Among the ATP competitive inhibitors reported to date [14–27], PF-3758309, a pan-PAK inhibitor [28,29], is the only one that proceeded to clinical testing [29,30]. However, this compound did not progress beyond Phase I owing to its low oral bioavailability [30]. In recent years, several Afraxis PAK inhibitors have been reported, such as FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 (Figure 1). These inhibitors, especially FRAX597, FRAX1036 and G-5555, have been widely used in recent studies involving PAKs inhibition [31–44]. FRAX414 was identified through a structure–activity relationship (SAR) study from the high-throughput screen (HTS) hit FRAX019. Further structural modification made to FRAX414 resulted in FRAX597 which acted as a group I PAK inhibitor with significant potential for the treatment of neurofibromatosis type 2 (NF2) and other cancers [22]. Compared with FRAX019 and FRAX414, FRAX597 shows higher inhibitory potency to group I PAKs and higher selectivity for the group I PAKs over the group II PAKs. Efforts to remove the aromaticity of FRAX597 resulted in FRAX1036 with improved drug properties and

Figure 1. Chemical structures of Afraxis PAK inhibitors (FRAX019, FRAX414, FRAX597, FRAX1036, and G-5555) and the
root-mean-square fluctuations of the inhibitors’ heavy atoms. According to the structural characteristics, the inhibitors are divided into five parts: Part A in red (head), Part B in blue, Part C in orange, Part D in purple, Part E in green (tail). The root-mean-square fluctuations (RMSFs) of the heavy atoms of each inhibitor binding with PAK1 or PAK4 are shown in parts accordingly. For example, in figure A’, the green bars indicate the RMSF values of the Part E heavy atoms of FRAX019 in the PAK1/FRAX019 system, the purple bars indicate the RMSF values of the Part D heavy atoms of FRAX019, and so on. The situations of figures A”, B’, B”, and so on are similar.
PAK: p21-activated kinase.

appropriate kinase selectivity [45]. To improve the limitations of FRAX1036, which is chiefly associated with its highly basic amine moiety, G-5555 was designed [31].
To date, cocrystal structures of PAK1 with certain inhibitors have been reported, while merely static interactions between them have been discussed [13,22,31]. To the best of our knowledge, the detailed inhibitory and selectivity mechanisms of this series of inhibitors are not clear. Hence, in this study, the interaction mechanisms of these five inhibitors with PAK1 (as the representative of group I PAKs) and PAK4 (as the representative of group II PAKs) were systematically investigated using molecular docking, molecular dynamics (MD) simulation, binding free energy calculation, binding energy decomposition analysis and other methods. Some suggestions for further modification of these inhibitors are put forward to improve the potency or selectivity. The computational results would help us understand the reasons leading to the different activities and selectivities of the five Afraxis PAK inhibitors toward the PAKs. Thus, this work would provide valuable information for the future design of specific and efficient PAK inhibitors.

Materials & methods
Preparation of structures
The cocrystal structures of the PAK1K299R/FRAX597 (Protein Data Bank [PDB] ID: 4EQC [22]), PAK1K299R/FRAX1036 (PDB ID: 5DFP [31]), PAK1D389N,T423E/G-5555 (PDB ID: 5DEY [31]) and
PAK4/FRAX486 (PDB ID: 5VEE [46]) complexes were retrieved from the RCSB Brookhaven PDB [47]. In Sybyl-X
2.0 software package [48], the missing residues of the PAKs were added, and FRAX597, FRAX1036 and G-5555 were sketched. According to the structure of FRAX597, FRAX019 and FRAX414 were constructed. The struc- tures of PAK1K299R/FRAX019, PAK4/FRAX019, PAK1K299R/FRAX414, PAK4/FRAX414, PAK4/FRAX597, PAK4/FRAX1036 and PAK4/G-5555 complexes were obtained by docking the corresponding inhibitors into the
receptors of PAK1K299R/FRAX597 (complexes PAK1K299R/FRAX1036 and PAK1D389N,T423E/G-5555 have the same binding pocket as it, see below) and PAK4/FRAX486 complexes, respectively, using Sybyl-X 2.0 software.
The electrostatic potential of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 were computed by Gaussian09 program [49] at the HF/6-31G* level, and the atomic partial charges of them were accomplished through the restrained electrostatic potential (RESP) fitting technique [50]. The partial charges and force field parameters for these inhibitors were generated using the antechamber module in AMBER16 [51]. The force field parameters for the phosphorylated threonine residue of PAK1 and the phosphorylated serine residue of PAK4 were
obtained from AMBER parameter database [52]. Since the above mutations of PAK1 were engineered mutations, the wild-type PAK1 (hereinafter referred to as PAK1) were re-mutated the targeting residues into the desired amino acids. Before the MD simulations were started, the missing atoms of the PAKs were added using the tleap program
in AMBER16. The systems were immersed in a rectangular TIP3P [53] water box of 10 A˚ from any solute atoms.
To neutralize the systems, suitable counterion Na+ (Cl-) was added to the PAK1/inhibitor (PAK4/inhibitor) complexes.

MD simulations
Preceding the MD simulations, each system was minimized via three steps using the pmemd.cuda program of AMBER16 running on a Tesla P100 card: first, the water molecules and ions in the systems were relaxed for 10,000 cycles (5000 cycles of steepest descent and 5000 cycles of conjugate gradient minimizations) with the receptor and ligand restrained. Second, the side chains of protein were optimized for 10,000 cycles (5000 cycles of steepest descent and 5000 cycles of conjugate gradient minimizations) with its backbone restrained. Finally, the entire system was minimized for 10,000 cycles (5000 cycles of steepest descent and 5000 cycles of conjugate gradient minimizations) without any restraint. In the molecular mechanics (MM) optimizations and the following MD simulations, the leaprc.protein.ff14SB [54] AMBER force field was used for the proteins and the general AMBER force field (GAFF) [55] for the ligands.
After minimization, 350 ps MD simulations were performed to the systems by gradually heating from 0 to 310 K. Then, 100 ns MD simulations were performed under the constant temperature of 310 K in the canonical ensemble using the pmemd.cuda program. In the MD simulations, the long-range coulombic electrostatics were calculated by using the particle mesh Ewald (PME) method [56], and all bonds connected to the hydrogen atom are fixed using
the SHAKE algorithm [57]. Periodic boundary conditions were applied to avoid edge effects in all calculations. The time step and the nonbonded cutoff were set to 2.0 fs and 10 A˚ , respectively. During the sampling process, the
coordinates were saved every 1.0 ps. After the MD simulations, 1000 snapshots of the simulated complex within the last 20 ns stable MD trajectory at 20 ps intervals were extracted to perform the further binding free energy calculations and decomposition analysis.

Binding free energy calculations
The binding free energy of each system was calculated by the Poisson–Boltzmann surface area (MM/PBSA) procedure based on the following equation [58]:
∆Gbind = Gcomplex − Gprotein − Gligand = ∆EMM + ∆GPB + ∆GSA − T∆S= ∆Ebind − T∆S
where ∆EMM is the gas–phase interaction energy between the receptor and the ligand, including ∆Eele (elec- trostatic energy), ∆Evdw (van der Waals energy) and ∆Eint (internal energy); ∆GPB and ∆GSA are the polar and nonpolar contributions to desolvation upon ligand binding, respectively; -T∆S represents the contribution of the conformational entropy (translation, rotation and vibration) at temperature T. The ∆GPB (polar desolvation energy)

here was calculated through the MM/PBSA method based on the Poisson–Boltzmann (PB) equations [59]. In the PB calculations, the grid size which was used to solve the PB equation was set to 2 A˚ , and the dielectric constant values 80 and 4 were set for solvent and solute, respectively. The nonpolar contribution to desolvation term was estimated from the Laboratoire de Chimie des Polymeres Organiques (LCPO) method: ∆GSA = 0.0072 ∆SASA [60], where
∆SASA is the variation of the solvent accessible surface area acquired with a solvent-probe radius of 1.4 A˚ . The sum
of ∆EMM, ∆GPB and ∆GSA was named as the binding energy (∆Ebind). The binding energy of each system was estimated by using the mm pbsa program in AMBER16. The conformation entropy (-T∆S) was calculated using normal-mode analysis via the nmode program in AMBER16 [51]. Considering expensive computational cost [61,62], only 100 snapshots evenly extracted from the last 20 ns were calculated here to estimate the conformation entropy.

Binding energy decomposition analysis
In order to figure out the detailed interactions between PAKs and the inhibitors, decomposition analysis of the binding energy was performed by using the mm pbsa program in AMBER16. In addition, the binding energy decomposition analysis for each system was computed using the generalized Born surface area (MM/GBSA) decomposition process because of the high computational demand of the PB calculations. The basic idea of decomposition was to decompose the total energy between each residue and the inhibitor into four energy terms: van der Waals contribution (∆Evdw), electrostatic contribution (∆Eele), polar solvation contribution (∆GGB) and nonpolar solvation contribution (∆GSA). The ∆Evdw and ∆Eele terms were computed by the sander program in AMBER16. The polar contribution of desolvation (∆GGB) was calculated using the generalized Born (GB) model [63] while the nonpolar contribution of desolvation (∆GSA) was computed based on SASA determined by the ICOSA algorithm [64]. All energy components were calculated using 1000 snapshots extracted from the MD trajectory from the last 20 ns. The molecular surface and cartoon representation models were generated using the PyMOL program [65].

Results & discussion
As shown in Figure 1, all the five PAK inhibitors (FRAX019, FRAX414, FRAX597, FRAX1036 and G- 5555) contain a pyrido[2,3-d]pyrimidin-7-one moiety. For convenience of description, we divide them into five parts (Parts A–E) according to the structural characteristics. As mentioned above, the crystal structures of the PAK1K299R/FRAX597, PAK1K299R/FRAX1036, PAK1D389N,T423E/G-5555 and PAK4/FRAX486 complexes
have been reported in the RCSB Brookhaven Protein Data Bank. We first analyzed and compared the structural characteristics of PAK1 and PAK4 based on the above crystal structures.

Overview of crystal structures
Considering the high mobility of the N-terminal lobe (N-lobe) and the structural conservation of the C-terminal lobe (C-lobe) [13,66], the cocrystal structures of PAK1K299R/FRAX597, PAK1K299R/FRAX1036, PAK1D389N,T423E/G-5555 and PAK4/FRAX486 are superimposed on the C-lobe (the following structures are superimposed in the same way) and the results are shown in Figure 2. Clearly, the PAK1 folds consist of a larger C-lobe dominated by α-helices, a smaller N-lobe dominated by β-sheets and a hinge region. Helix-A is located at the C-terminus of the kinase, P-loop (phosphate-binding loop) is located on the N-lobe side of the cleft and Helix-C is located at the bottom of the active pocket. The alignment reveals that the regions with significant differences near the pocket in the three PAK1/inhibitor complexes are Helix-A at the C-terminus, and P-loop, which is similar to the results reported by Joachim Rudolph – among others [13]. We speculate that these regions may be flexible and prone to conformational changes, affecting the activities of inhibitors. Furthermore, FRAX597, FRAX1036 and G-5555 are all anchored to the hinge of PAK1 by forming two hydrogen bonds (see Figure 2). One is formed between the N3 nitrogen of the pyrido[2,3-d]pyrimidin-7-one (denoted as Part C) and the amide nitrogen of Leu347, and the other is formed between the NH group of Part C and the carbonyl oxygen of Leu347. It is noted that the anchored positions are all the same for these three inhibitors and the head (Part A) faces Helix-C located at the bottom of the active pocket of PAK1. Obviously, the tails (Part E) of FRAX597 and FRAX1036 are exposed to the solvent, which should be detrimental to the formation of specific interaction between this moiety and the kinase [13].
It can be found from Figure 2 that the active pocket of PAK4 is structurally similar to that of PAK1 and FRAX486 is also anchored at the hinge area of PAK4 by two hydrogen bonds. In the following sections, molecular dynamics simulation and subsequent computational analysis were reported and discussed to further disclose the

Figure 2. Superimposed cocrystal structures of PAK1/FRAX597 (magenta), PAK1/FRAX1036 (yellow), PAK1/G-5555 (salmon) and PAK4/FRAX486 (slate). The dotted lines represent the hydrogen bonds. In the inset, the red dashed lines denote the hydrogen bonds formed between the residue Leu347 and the respective inhibitor in the PAK1/FRAX597, PAK1/FRAX1036 and
PAK1/G-5555 complexes, and the blue dotted lines represent the hydrogen bonds formed by the residue Leu398PAK4 with FRAX486.

binding modes of PAK1/FRAX019, PAK4/FRAX019, PAK1/FRAX414, PAK4/FRAX414, PAK1/FRAX597, PAK4/FRAX597, PAK1/FRAX1036, PAK4/FRAX1036, PAK1/G-5555 and PAK4/G-5555 complexes.

The stability & flexibility of system simulations
We performed classic MD simulations on the above ten PAKs/inhibitor complexes to predict their binding modes and explore the key residues associated with the activity and selectivity. To ensure the rationality of the sampling strategy and clarify the dynamic stability of these complexes, 100 ns MD simulations were performed under the constant temperature of 310 K for each complex. Then the root-mean-square deviations (RMSD) of the protein backbone atoms (CA, C and N) and ligands in the ten complexes were calculated based on their initial structures, and the results are plotted in Figure 3 and Supplementary Figure 1, respectively. It shows that all the systems generally reach equilibrium after 50 ns, and then remain stable to the end of the simulations. Even though some of them have small fluctuations, they are still within the acceptable scale. Moreover, from Supplementary Figure 2 we can see that the potential energies of the systems are almost unchanged after 20 ns MD simulations. Therefore, MD trajectories between 80 and 100 ns being applied to conduct the following binding free energy and binding energy decomposition calculations should be reasonable and reliable.

Binding energies predicted by MM/PBSA & decomposition analysis by MM/GBSA
The predicted binding free energies, the energy components of all systems, and IC50/Ki derived from the experi- ments [22,31,45] are summarized in Table 1. The calculated results indicate that the binding free energies (∆Gbind in Table 1) of FRAX597, FRAX1036 and G-5555 toward PAK1 (-61.30, -56.55 and -55.88 kcal/mol) are all lower than those of FRAX019 and FRAX414 (-37.81 and -53.24 kcal/mol), which is in good agreement with the trend of experimental activity data. Regrettably, the predicted binding free energies of the PAK4/inhibitor systems are not consistent with the experimental trend, which may be because that the drug–receptor interactions are controlled both thermodynamically and kinetically. Since the calculation of kinetic processes is very time- consuming, here we only focus on the binding mode of each inhibitor from the perspective of thermodynamics. It is obvious that the van der Waals (∆Evdw) component is the major favorable contributor in the binding energy and plays a vital role for the inhibitor binding in all the ten complexes (-53.11, -40.40, -58.55, -54.37, -73.47,
-57.88, -66.80, -63.47, -61.90 and -50.56 kcal/mol for PAK1/FRAX019, PAK4/FRAX019, PAK1/FRAX414, PAK4/FRAX414, PAK1/FRAX597, PAK4/FRAX597, PAK1/FRAX1036, PAK4/FRAX1036, PAK1/G-5555
and PAK4/G-5555, respectively), followed by the electrostatic (∆Eele) contribution. As shown in Supplemen- tary Table 1, the net electrostatic contributions (∆Eele+∆GPB) are also favorable for the binding of receptors with ligands (-4.19, -6.21, -8.54, -10.02, -5.57, -12.78, -10.16, -16.56, -15.37 and -17.87 kcal/mol for these ten systems, respectively) although the polar contribution of desolvation (∆GPB) could partly offset the electrostatic interaction (∆Eele). The percentage ratios of the nonpolar interaction part (∆Evdw+∆GSA) in the binding energy
(∆Ebind in Table 1) are more than 76% (93.46, 88.20, 88.61, 86.05, 93.67, 83.70, 88.12, 81.27, 81.97 and
76.35 for the ten systems, respectively), indicating that the binding of these complexes are mainly dependent on

3 PAK1/FRAX019 2
1
3 PAK4/FRAX019 2
1

3 PAK1/FRAX414 2
1
3 PAK4/FRAX414 2
1

3 PAK1/FRAX597 2
1

4 PAK4/FRAX597 3
2

3 PAK1/FRAX1036 2
1

3 PAK4/FRAX1036 2
1

3 PAK1/G-5555
2
1

3 PAK4/G-5555
2
1
0 20

40 60 80
Time (ns)

100

Figure 3. The root-mean-square deviations of the protein backbone atoms (CA, C and N) in the ten complexes relative to the respective first snapshot.
RMSD: Root-mean-square deviation.

nonpolar interactions. Furthermore, the influences of nonpolar interactions on the inhibitor’s kinase selectivity ([∆(∆Evdw+∆GSA)/∆(∆Ebind)] PAK1/inhibitor-PAK4/inhibitor) are analyzed. As shown in Supplementary Table 1, the percentage values of [∆(∆Evdw+∆GSA)/∆(∆Ebind)] PAK1/inhibitor-PAK4/inhibitor for all the inhibitors are much larger
than 100% except for FRAX1036, indicating that the non-polar interactions should be responsible for the gen- eration of kinase selectivities toward PAK1, while the polar interactions are detrimental from the viewpoint of

Table 1. The predicted binding free energies (kcal/mol) and individual energy terms and biological activities of the PAK1/FRAX019, PAK4/FRAX019, PAK1/FRAX414, PAK4/FRAX414, PAK1/FRAX597, PAK4/FRAX597, PAK1/FRAX1036,
PAK4/FRAX1036, PAK1/G-5555, PAK4/G-5555 complexes, together with the biological activity reported in the
literature.
Complex Δ Eele † Δ Evdw ‡ Δ GPB § Δ GSA ¶ Δ Ebind# -TΔ S Δ Gbind Biological activity [μM]
PAK1/FRAX019 -14.30 ± 4.22 -53.11 ± 3.24 10.11 ± 1.41 -6.93 ± 0.23 -64.24 ± 4.84 26.43 ± 4.81 -37.81 0.483††
PAK4/FRAX019 -14.44 ± 3.90 -40.40 ± 3.34 8.23 ± 1.58 -5.93 ± 0.29 -52.53 ± 4.13 19.62 ± 6.89 -32.91 0.580††
PAK1/FRAX414 -20.19 ± 5.15 -58.55 ± 3.43 11.65 ± 1.38 -7.87 ± 0.20 -74.96 ± 5.20 21.72 ± 11.32 -53.24 0.020††
PAK4/FRAX414 -19.93 ± 4.00 -54.37 ± 3.72 9.91 ± 1.25 -7.38 ± 0.27 -71.76 ± 4.80 21.63 ± 10.49 -50.13 0.559††
PAK1/FRAX597 -20.15 ± 4.09 -73.47 ± 3.58 14.58 ± 1.54 -8.97 ± 0.19 -88.01 ± 4.98 26.71 ± 6.01 -61.30 0.008††
PAK4/FRAX597 -24.14 ± 4.91 -57.88 ± 3.95 11.36 ± 1.26 -7.74 ± 0.34 -78.40 ± 6.08 23.70 ± 6.82 -54.70 >10††
PAK1/FRAX1036 -24.11 ± 5.40 -66.80 ± 3.54 13.95 ± 1.64 -8.64 ± 0.24 -85.61 ± 5.80 29.06 ± 7.22 -56.55 0.023‡‡
PAK4/FRAX1036 -29.13 ± 5.94 -63.47 ± 3.54 12.57 ± 1.32 -8.33 ± 0.32 -88.35 ± 6.22 23.50 ± 5.03 -64.85 2.4‡‡
PAK1/G-5555 -32.05 ± 7.50 -61.90 ± 3.16 16.68 ± 1.96 -7.95 ± 0.16 -85.21 ± 6.80 29.33 ± 4.45 -55.88 0.004‡‡
PAK4/G-5555 -30.27 ± 5.39 -50.56 ± 3.25 12.40 ± 1.46 -7.14 ± 0.27 -75.57 ± 4.76 23.99 ± 5.67 -51.58 unknown
†Electrostatic contribution.
‡van der Waals contribution.
§The polar contribution of desolvation.
¶The nonpolar contribution of desolvation.
# The binding energy.
††IC50 [22].
‡‡Ki derived from the experiment [31,45].
IC50 : Half maximal inhibitory concentration; -TΔ S: Conformation entropy.

thermodynamics. The opposite is true for FRAX1036. In addition, from Supplementary Table 1, it can be found that solvation is unfavorable for the binding of these inhibitors with the kinases.
In order to further identify which residues are responsible for the binding of these protein/inhibitor complexes,
MM/GBSA binding energy decomposition analysis was performed to decompose the binding energy into the contribution of each residue. For each complex, the key residues with great contributions (∆Ebind <-1.00 kcal/mol) to the binding energy are shown in Figure 4. Furthermore, intermolecular hydrogen bonds and key intramolecular hydrogen bonds have been calculated during the last 20 ns MD simulations because hydrogen bonds play an important role in protein–drug interactions, and the results are summarized in Table 2. Detailed analysis of the above data will be carried out in the following sections.

The binding modes of FRAX019, FRAX414, FRAX 597, FRAX 1036, & G-5555 with PAK1
Figure 5 shows the binding modes of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 with PAK1. As shown in Figure 5A, the binding modes of the five inhibitors with PAK1 are similar to each other. The calculations indicate that the binding pocket of PAK1 consists of three parts: Gly350, Leu396 and Lys538 in C-lobe; Ile276, Val284, Ala297, Lys299, Glu315 (belonging to helix-C) and Met344 in N-lobe; Tyr346 and Leu347 in the hinge region, which is consistent with the crystal structures. Obviously, this pocket is big enough to accommodate Parts A, B and C of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555. It gives a possible explanation of the fact that these five inhibitors have inhibitory activities to PAK1. Residues Leu396 from C-lobe, Ile276 and Val284 from N-lobe and Tyr346 and Leu347 from the hinge region are surrounding the pyrido-[2,3-d]pyrimidin-7-one cores (Part C) of the inhibitors. They establish extensive interactions with the inhibitors, which is confirmed by the results of binding energy decomposition (Figures 4 & 5A). Ala297, Lys299, Glu315 and Met344 deep in the active pocket of PAK1 can interact effectively with FRAX414, FRAX597, FRAX1036 and G-5555 containing a large head (Figure 1, Part A). In addition, our calculations reproduce well that FRAX597, FRAX1036 and G-5555 can form two anchoring hydrogen bonds with PAK1, respectively (Table 2). Furthermore, the simulation results show that FRAX019 and FRAX414 can also form these two anchoring hydrogen bonds with PAK1 at the same position. This anchoring effect may be the structural basis of substitution modification of the inhibitors to effectively change their PAK1 binding affinities. Since the binding energy decomposition does not include the conformational entropy, the values of binding energy (∆Ebind in Table 1) are used in the following discussion.

Residue
0

-1

-2

-3

-4

Residue
0

-1

-2

-3

Residue
0

Residue
0

-1
-1
-2
-2
-3

-4 -3

Residue
0

PAK1/FRAX597
0

Residue

-1 -1

-2 -2

-3 -3

Residue
0

PAK1/FRAX1036
0

Residue

-1
-1
-2
-2
-3

-3 -4

Residue
0

Residue
0

-1 -1

-2 -2

-3 -3

Figure 4. The binding energy spectra of key residues in complexes according to the generalized Born surface area decomposition analysis. (A) PAK1/FRAX019, (B) PAK4/FRAX019, (C) PAK1/FRAX414, (D) PAK4/FRAX414, (E) PAK1/FRAX597, (F) PAK4/FRAX597, (G) PAK1/FRAX1036, (H) PAK4/FRAX1036, (I) PAK1/G-5555 and (J) PAK4/G-5555.

FRAX019
FRAX414 FRAX019

FRAX597 FRAX1036 G-5555

FRAX414

FRAX414 FRAX597

FRAX597 FRAX1036

FRAX1036 G-5555

Figure 5. Binding modes of the inhibitors with key residues in the active pocket of PAK1. (A) FRAX019 (green), FRAX414 (cyan), FRAX597 (magenta), FRAX1036 (yellow) and G-5555 (salmon) in the active pocket of PAK1; (B) PAK1/FRAX019 (green) versus PAK1/FRAX414 (cyan); (C) PAK1/FRAX414 (cyan) versus PAK1/FRAX597 (magenta); (D) PAK1/FRAX597 (magenta) versus PAK1/FRAX1036 (yellow); (E) PAK1/FRAX1036 (yellow) versus PAK1/G-5555 (salmon).

Table 2. Visible percentage of hydrogen bonds during the last 20 ns molecular dynamics simulations between the Afraxis p21-activated kinase inhibitors and the p21-activated kinases.†
Complex Donor Acceptor Occupied (%) Distance (A˚ ) Angle (deg)
PAK1/FRAX019 :FRAX019@N4:FRAX019@H18 :L347@O 99.12 2.927 157.60
:L347@N:L347@H :FRAX019@N2 60.48 3.432 151.00
PAK4/FRAX019 :FRAX019@N4:FRAX019@H18 :L398@O 89.82 2.978 146.38
:L398@N:L398@H :FRAX019@N2 97.91 3.093 158.47
PAK1/FRAX414 :FRAX414@N4:FRAX414@H26 :L347@O 98.47 2.874 150.40
:L347@N:L347@H :FRAX414@N2 95.38 3.171 158.26
PAK4/FRAX414 :FRAX414@N4:FRAX414@H26 :L398@O 94.16 3.052 153.89
:L398@N:L398@H :FRAX414@N2 95.39 3.110 152.74
PAK1/FRAX597 :FRAX597@N3:FRAX597@H16 :L347@O 98.58 2.915 152.62
:L347@N:L347@H :FRAX597@N5 98.63 3.086 160.32
:H545@N:H545@H :I276@O 0.63 2.87 160.08
:N302@N:N302@H :S281@O 0.57 2.87 162.75
PAK4/FRAX597 :FRAX597@N7:FRAX597@H28 :L398@O 98.66 2.970 155.93
:L398@N:L398@H :FRAX597@N3 96.99 3.133 160.07
PAK1/FRAX1036 :FRAX1036@N4::FRAX1036@H8 :L347@O 95.40 2.925 152.91
:L347@N:L347@H :FRAX1036@N1 98.40 3.080 158.23
:N544@ND2:N544@HD22 :K275@O 95.14 2.951 157.55
PAK4/FRAX1036 :FRAX1036@N5:FRAX1036@H30 :L398@O 69.50 2.850 156.86
:L398@N:L398@H :FRAX1036@N7 20.09 2.933 159.86
PAK1/G-5555 :G-5555@N2:G-5555@H2 :L347@O 98.32 2.965 150.34
:L347@N:L347@H :G-5555@N3 94.88 3.145 159.59
PAK4/G-5555 :L398@N:L398@H :G-5555@N1 47.79 2.916 160.66
:G-5555@N2:G-5555@H6 :L398@N 47.28 2.870 154.28
†A hydrogen bond was defined if the donor–acceptor (H) distance was less than 3.5 A˚ and the donor–acceptor (H) angle was larger than 120◦.

PAK1/FRAX019
In the PAK1/FRAX019 complex, the cyclopentyl group (Figure 1, Part B of FRAX019) is directed toward the crack between the C-lobe and N-lobe and the methylpiperazine moiety (Part E) is oriented toward the outlet of the active pocket (Figure 5B). Both Part B and Part E point to the same side of the pyrido-[2,3-d]pyrimidin-7-one core (Part C), and the entire inhibitor has a C-shaped profile. By forming two hydrogen bonds with Leu347 of the hinge region (Figure 5B & Table 2), FRAX019 is anchored near the outlet of the active pocket and establishes extensive interactions with the surrounding residues Ile276 (-4.02 kcal/mol), Val284 (-1.64 kcal/mol), Tyr346 (-2.00 kcal/mol), Leu347 (-1.63 kcal/mol), Gly349 (-1.23 kcal/mol), Gly350 (-1.07 kcal/mol), Leu396 (-2.72 kcal/mol), Thr406 (-1.67 kcal/mol) and Lys538 (-1.1 kcal/mol). The van der Waals interactions between these residues and the inhibitor are dominant (Supplementary Table 2). Among these key residues, Ile276, Val284, Tyr346, Leu347, Leu396 and Thr406 well surrounding Part B and Part C of FRAX019 (Figures 1 & 5B) are significantly favorable for FRAX019 binding, which can be confirmed by the residue energy contributions (Figure 4A). However, although the anchoring hydrogen bonds can fix the inhibitor in the kinase active pocket, they may prevent the inhibitor from entering the deep part of the active pocket. As can be seen in Figure 5B, FRAX019 can only occupy nearly half of the pocket cavity, leaving its tail (Part E) exposed outside the active pocket. We also calculated the root-mean-square fluctuation (RMSF) of the inhibitor’s heavy atoms in the PAK1/FRAX019 complex, and the results are summarized in Figure 1. It is clear that the RMSF value of Part E of FRAX019 is much larger than those of the other parts, indicating that the methylpiperazine moiety is relatively free. The results of the binding energy decomposition calculation indeed show that FRAX019 cannot form effective interactions with Ile298, Lys299, Val342, Glu315 and Met319 deep in the active pocket and that Part E lacks specific interactions with the kinase. These binding characteristics should be the reason why FRAX019 has a low PAK1 activity.

PAK1/FRAX414
The results of binding energy decomposition show that FRAX414 can bind to ten surrounding key residues of PAK1 (Ile276, Val284, Ala297, Lys299, Met344, Tyr346, Leu347, Gly350, Leu396 and Lys538, and their energy contributions are -3.36, -1.69, -1.52, -2.00, -2.02, -2.26, -1.62, -1.44, -2.40 and -1.53 kcal/mol, respectively). FRAX414 is obtained by replacing the hydrogen atom with a 2,4-dichlorophenyl group (Part A) and the cyclopentyl group with an ethyl group (Part B) (see Figure 1, FRAX414 and FRAX019). Apparently, FRAX414 increases the inhibitor’s head (Part A). Compared with FRAX019, the substituent 2,4-dichlorophenyl (Part A) of FRAX414 buries more deeply into the active pocket of PAK1 (Figure 5B). This moiety establishes extensive interaction with the surrounding hydrophobic residues Ala297, Ile298, Lys299, Val342 and Met344. The energy contribution of these residues was 0.71, 0.53, 1.76, 0.89 and 1.18 kcal/mol larger than those in PAK1/FRAX019, respectively. However, the energy contributions of Ile276, Gly349 and Thr406 are 0.66, 0.70 and 1.05 kcal/mol smaller than those of PAK1/FRAX019, respectively. This is attributed to the weakening of the van der Waals interactions between these residues and FRAX414 (Supplementary Table 2) due to the substitution of the high molecular weight cyclopentyl (Part B) of FRAX019 by the low molecular weight ethyl (Part B) of FRAX414. Therefore, FRAX414 showing the higher experimental activity than FRAX019 is primarily ascribed to the introduction of the 2,4-dichlorophenyl group (Part A). However, the FRAX414 with an increased head still does not interact well with the residues located at the bottom of the pocket, such as Glu315 (-0.18 kcal/mol) and Met319 (-0.36 kcal/mol).

PAK1/FRAX597
Compared with FRAX414, FRAX597 has an additional thiazolyl group in Part A (Figure 1). Since FRAX597 is anchored at the same position as FRAX019 and FRAX414, the longer head (Part A) of the FRAX597 has significant interactions with Glu315 (-0.80 kcal/mol) and Met319 (-0.94 kcal/mol) on the helix-C located at the bottom of the active pocket, as well as with Val342 (-1.44 kcal/mol) deep in the active pocket (Figure 5C). The energy contributions of the above three residues, compared with that in the PAK1/FRAX414 system, are increased by 0.62, 0.58 and 0.52 kcal/mol, respectively. The above facts indicate that the longer Part A of FRAX597 effectively enhances the affinities of the inhibitor with the residues deep in the pocket as well as the residues at the bottom of the pocket. It is noteworthy that FRAX597 can stably capture the C-terminus of PAK1 and promote the formation of an effective hydrogen bond between the terminal residues His545 and Ile276 on the N-lobe β-sheet (Table 2), leading to closer distances between FRAX597 and the two residues Lys542 and His545. This conformation changes lead to that the binding affinities of FRAX597 with Lys542 and His545 are increased by 0.78 and 2.45 kcal/mol, respectively, compared with FRAX414. Whereas, the contribution of Ile276 to the binding energy in PAK1/FRAX597 is 0.92 kcal/mol smaller than that in PAK1/FRAX414. Through analysis, an effective hydrogen bond is formed between Ser281 and Asn302 in the P-loop (Table 2). This hydrogen bond may pull the P-loop (connected with the β-sheet) toward the bottom of the active pocket, resulting in a change in the position of Ile276 (Supplementary Figure 3A). As mentioned above, Part A in FRAX597 is long enough to establish a stable interaction with the bottom of the binding pocket. Therefore, it is no need to further increase its length.

PAK1/FRAX1036
As shown in Figure 5D, the binding position of FRAX1036 in PAK1 is highly coincident with that of FRAX597. Compared with FRAX597, FRAX1036 has a slightly bigger head group (see Part A of FRAX1036 and FRAX597 in Figure 1), resulting in an enhanced van der Waals interaction with Glu315 at the bottom of the pocket, while also enhancing the electrostatic and van der Waals interactions with Lys299 (Supplementary Table 2). Although the aromaticity-removing substitution (Part D of FRAX1036) enhances its conformational stability in the active pocket (as can be seen from the RMSF values in Figure 1), it also leads to a weakening of its affinities with the surrounding residues, such as Ile276 (-2.03 kcal/mol) and Gly350 (-1.16 kcal/mol). In addition, FRAX1036 can help the N-lobe β-sheet of PAK1 capture the C-terminus. Nonetheless, hydrogen bonds between Lys275 and Asn544 are formed (Table 2), which makes it impossible for FRAX1036 to establish effective interactions with Lys542, Asn544 and His545. These binding characteristics determine that Part A of FRAX1036 has a slightly stronger interaction with PAK1, compared with FRAX597, while Part D have slightly weaker interactions with PAK1, maintaining a considerable affinity overall (-85.61 and -88.01 kcal/mol). However, the drug properties of FRAX1036 have been improved, so this effort to remove aromaticity is relatively successful.

PAK1/G-5555
Our calculation results show that the binding affinities of PAK1 with G-5555 is very similar to that with FRAX1036 (-85.21 and -85.61 kcal/mol). Compared with the other four inhibitors, G-5555 removes Part E at the tail and its Part D (Part B) is much smaller (larger) than the others (see Figure 1). Since G-5555 is anchored at the same position as the other inhibitors, removing Part E should not affect its binding affinity with PAK1 (see above). However, the smaller Part D may reduce its binding affinities with Ile276, Tyr346, Ala348, Gly350 and Lys538 located at the exit of the pocket (Figure 5E). This is indeed confirmed by that their energy contributions are 0.84, 0.34, 0.80, 0.86 and 1.51 kcal/mol smaller than that in PAK1/FRAX1036, respectively (Supplementary Table 2). Interestingly, without the spatial barrier of the big tail, PAK1’s Helix-A tends to fill the outlet of the active pocket, making the binding affinity of G-5555 with Thr541 (-1.23 kcal/mol) on this helix significantly increase compared with the other four inhibitors. At the same time, the interactions of the large Part B in G-5555 with Asp393 (-0.48 kcal/mol) at the C-lobe crack and Leu396 (-2.68 kcal/mol) are significantly enhanced (the other four inhibitors cannot interact with Asp393). Thus, the modification from FRAX1036 to G-5555 is advisable.
In addition, it can be observed from Figure 5A that Part Bs of FRAX414, FRAX597 and FRAX1036 are almost free to rotate, which should be attributed to that this portion is a small ethyl group (see Figure 1). The binding energy decomposition results indicate that the van der Waals interactions of these three inhibitors with the Part B adjacent residues Asp393 (-0.06PAK1/FRAX414, -0.05PAK1/FRAX597 and -0.05PAK1/FRAX1036 kcal/mol) and Leu396 (-2.07PAK1/FRAX414, -2.22PAK1/FRAX597 and -1.91PAK1/FRAX1036 kcal/mol) are indeed significantly weaker than that of FRAX019 and G-5555 with Asp393 (-0.13PAK1/FRAX019 and -0.46PAK1/G-5555) and Leu396 (-2.47PAK1/FRAX019
and -2.41PAK1/G-5555). Therefore, an appropriate increase in the volume of Part B should help to increase the affinity
of an inhibitor with PAK1.
As a result, for the binding of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 with PAK1, Leu347 plays a crucial role. This residue can form two hydrogen bonds with Part Cs of these five inhibitors and anchor them at the same position in the active pocket of PAK1. The anchoring effect is very important for qualitative modification of these inhibitors. It is due to the anchoring effect of hydrogen bonds that increasing the length of the inhibitor’s head (Part A) and removing its tail (Part E) are good choices for designing efficient inhibitors on the basis of FRAX019. In addition, adding appropriate volume and positive charge (see below) to Part B should increase the affinity of this part.

The selectivity mechanisms of FRAX019, FRAX414, FRAX597, FRAX1036, & G-5555 toward PAK1 & PAK4
As shown in Table 1, the calculated binding free energies of FRAX414, FRAX597 and G-5555 with PAK4 are 3.11,
6.60 and 4.30 kcal/mol larger than those with PAK1, respectively, showing that they all have certain selectivity to PAK1 and the latter two are better than the former, which is in agreement with the experiment. However, the calculation results for FRAX019 and FRAX1036 are inconsistent with the experiment. This situation may be caused by a kinetical effect on the process of binding between the inhibitor and PAKs being neglected computationally. Considering that the kinetic calculations for the ten complexes are time consuming, here we discuss the selectivities only from a thermodynamic point of view.
Figure 6 shows the binding modes of the five inhibitors with the key residues in PAK4. For convenience of comparison, the binding modes of the five PAK1/inhibitor complexes are also shown in Figure 6. Obviously, for each inhibitor, the binding mode with PAK4 is very similar to that with PAK1. Furthermore, our calculations show that Part C in each of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 is anchored in the hinge region located at the exit of the binding pocket of PAK4 by forming two hydrogen bonds with Leu398 (Table 2), similar to the case of PAK1. In the following sections, we will compare the binding modes of each of the five inhibitors to
PAK4 and PAK1 in detail, respectively, and try to find clues to enhance selectivities.

Comparison of the binding modes of PAK4/FRAX019 & PAK1/FRAX019 complexes
From Figure 4 and Supplementary Figure 4, it can be found that the key residues in PAK4/FRAX019 are identical (or conservatively substituted residues) to the corresponding key residues in PAK1/FRAX019. Differently, the hinge region with the anchoring residue (Leu398) of PAK4 is closer to the interior of the cavity compared with the situation of PAK1, resulting in FRAX019 closer to the crack in PAK4 than in PAK1 (Figure 6A). The hinge region sequence residues of PAK4 and PAK1 are ‘M395-EFLE-G400’ and ‘M344-EYLA-G349,’ respectively. From the superposition of the simulated structures (Figure 6A), it could be found that, compared with Glu345PAK1, the

PAK1/FRAX019 PAK4/FRAX019

PAK1/FRAX597 PAK4/FRAX597

PAK1/FRAX414 PAK4/FRAX414

PAK1/FRAX1036 PAK4/FRAX1036

PAK1/G-5555 PAK4/G-5555

Figure 6. Binding modes of the five inhibitors with key residues in PAK4, together with those in PAK1. (A) PAK1/FRAX019 (marine) versus PAK4/FRAX019 (hotpink); (B) PAK1/FRAX414 (marine) versus PAK4/FRAX414 (hotpink); (C) PAK1/FRAX597 (marine) versus PAK4/FRAX597 (hotpink); (D) PAK1/FRAX1036 (marine) versus PAK4/FRAX1036 (hotpink); (E) PAK1/G-5555 (marine) versus PAK4/G-5555 (hotpink).

side chain orientation of Glu396PAK4 occurs at nearly 90◦, and the Glu396PAK4 also shifts toward the inside of the cavity. In addition, the side chains of both Phe397PAK4 and Tyr346PAK1 are almost parallel to Part C of FRAX019,
and the position of the Phe397PAK4 also shifts toward the interior of the cavity. Beyond that, there are no other significant differences in the position and side chain orientation between both the hinge residues of PAK4 and PAK1. The difference in the side chain orientation of Glu396PAK4 and Glu345PAK1 causes that the distances of FRAX019 with the residues Ile327PAK4, Gly400PAK4, Leu447PAK4 and Ser457PAK4 be longer than that with the corresponding residues Ile276PAK1, Gly349PAK1, Leu396PAK1 and Thr406PAK1 (Supplementary Figure 3B, C, D & E). This is also confirmed by the difference in the van der Waals interactions of these residues with FRAX019 in both complexes (Supplementary Table 2). It is noted that, different from the case of PAK1, FRAX019 does not effectively capture the C-terminus of PAK4 and thus cannot establish effective interactions with Arg589PAK4 (-0.01 kcal/mol), but Arg586PAK4 (-1.68 kcal/mol) at the C-terminus also has a strong affinity with FRAX019 due to its long side chain. Therefore, capturing the C-terminus or not should have nonsignificant effects on the kinase selectivity of FRAX019. Overall, FRAX019 shows certain PAK1 selectivity.

Comparison of the binding modes of PAK4/FRAX414 & PAK1/FRAX414 complexes
It can be found from Figure 6B that in PAK4/FRAX414, the opening degree of the β-sheet on the N-lobe side of the pocket crack is slightly larger than that in PAK1/FRAX414, resulting in the distance of the inhibitor with Ile327PAK4 being slightly longer than that with Ile276PAK1 (Supplementary Figure 3F). This difference in distance should be beneficial for the selectivity to PAK1, which is confirmed by the binding affinity between the residues and the inhibitor in PAK4/FRAX414 being smaller than that in PAK1/FRAX414 (Figure 4C & D). The side chain of Ile337PAK4 is significantly larger than that of the corresponding residue Thr286PAK1, leading to the distance of the inhibitor with the former being slightly shorter than that with the latter (Supplementary Figure 3G). This directly results in the contribution of Ile337PAK4 to the binding energy being 0.68 kcal/mol larger than that of Thr286PAK1, reducing the selectivity of FRAX414 to PAK1. Similar to the case of FRAX019, FRAX414 also fails to capture the C-terminus of PAK4, which should not have any effects on FRAX414’s kinase selectivity.

Comparison of the binding modes of PAK4/FRAX597 & PAK1/FRAX597 complexes
The simulation results (Figure 6C) show that FRAX597 is not as deep in the active pocket of PAK4 as in that of PAK1; because of this, the binding affinity of FRAX597 with Lys350PAK4 is 0.35 kcal/mol smaller than that with Lys299PAK1. As shown in Figures 4E, F & 6C, the key residues that contribute differently to the receptor– ligand binding energy in the two systems concentrate on the C-terminus of the kinases. Similar to FRAX019 and FRAX414, FRAX597 also fails to capture the C-terminus of PAK4 but can form a strong interaction with Arg586PAK4 (-1.28 kcal/mol). However, the interaction of FRAX597 with the C-terminus of PAK4 is apparently weaker than that of PAK1, different from the cases of FRAX019 and FRAX414. This should be attributed to the formation of the hydrogen bond between His545 and Ile276 in PAK1/FRAX597 complex which makes the inhibitor to be hugged by the nearby C-terminal residues Lys538 (-1.21 kcal/mol), Lys542 (-1.22 kcal/mol), Asn544 (-0.83 kcal/mol), and His545 (-2.46 kcal/mol) (see above). Consistent with the work of Zhang et al. [46], our calculations could not prove the hypothesis that the selectivity comes from the conformational preferences of Met395PAK4, Met344PAK1, Lys350PAK4 and Lys299PAK1 [13].

Comparison of the binding modes of PAK4/FRAX1036 & PAK1/FRAX1036 complexes
As shown in Figure 6D, the binding mode of PAK4/FRAX1036 is very similar to that of PAK1/FRAX1036. Since the side chain of Ile337PAK4 is slightly larger than that of Thr286PAK1, the van der Waals interaction of FRAX1036 with the former is 0.55 kcal/mol larger than that with the latter. Compared with PAK1, the head (Part A) of FRAX1036 is clearly biased toward the N-lobe of the kinase in PAK4. This difference results in significantly stronger affinities of FRAX1036 with the nearby Lys350PAK4 (-3.19 kcal/mol) and Val393PAK4 (-2.17 kcal/mol) than with the corresponding residues Lys299PAK1 (-2.32 kcal/mol) and Val342PAK1 (-1.70 kcal/mol), and weaker interactions of FRAX1036 with Glu366PAK4 (-0.63 kcal/mol) on helix-C than with Glu315PAK1 (-1.08 kcal/mol). In addition, FRAX1036 cannot capture the C-terminus of PAK4, which should have no effects on its kinase selectivity, consistent with the cases of FRAX019 and FRAX414. It seems that FRAX1036 has certain selectivity toward PAK4 thermodynamically, which is not consistent with the experiment.

Figure 7. FRAX597 and surface charge distributions generated using PyMOL. (A) PAK1 and (B) PAK4.

Comparison of the binding modes of PAK4/G-5555 & PAK1/ G-5555 complexes
Comparing the simulated structures of PAK4/G-5555 and PAK1/G-5555 (Figure 6E), we can find that the head (Part A) of the G-5555 is more biased toward the C-lobe side in PAK4, which is exactly opposite of the case of FRAX1036. This binding characteristic results in that the affinities of G-5555 with Val393PAK4 (-1.25 kcal/mol) is significantly smaller than that of the corresponding residue Val342PAK1 (-1.79 kcal/mol), while Phe459PAK4 (-0.52 kcal/mol) is slightly larger than Phe408PAK1 (-0.17 kcal/mol). At the same time, Part B of G-5555 is a little far from the C-lobe side in PAK4, resulting in the affinity of Asp444PAK4 here being 0.40 kcal/mol smaller than that of Asp393PAK1, also accompanied by the affinity of Leu447PAK4 (-2.16 kcal/mol) being smaller than that of Leu396PAK1 (-2.68 kcal/mol). When G-5555 is combined with PAK4, the β-sheet at the crack of N-lobe slightly moves outward compared with PAK1, leading to a slightly lower affinity of Ile327PAK4 (-0.80 kcal/mol) than that of Ile276PAK1 (-1.19 kcal/mol). Similar to FRAX597, G-5555 fails to capture the C-terminus of PAK4 and has weaker interactions with the C-terminus of PAK4 than PAK1. Therefore, G-5555 has a much stronger binding affinity with PAK1 than PAK4, showing a good selectivity, which is in agreement with the experimental result.
Furthermore, taking PAK1/FRAX597 and PAK4/FRAX597 complexes as examples, we analyzed the surface charge distribution maps of PAK1 and PAK4 generated using PyMOL and the results were shown in Figure 7. It can be seen that the charge distributions of both active pockets are very similar. The cracks of the kinase active pockets near Part B of FRAX597 are negative charge enriched mainly due to the existence of Asp393PAK1 and Asp444PAK4. Therefore, adding appropriate positive charge to Part B should increase the affinity of this part and thus probably enhance the activity. It can only be seen that PAK1 has more negative charge on the C-lobe near Part D of the inhibitor compared with PAK4. Therefore, we speculate that switching Part D with a positively charged substituent may help to improve the PAK1 selectivity of the inhibitor.
Based on the above analyses, we can obtain some conclusions. Leu398 of PAK4 had similar anchoring effects on these five Afraxis PAK inhibitors by forming two hydrogen bonds with each of them. The five inhibitors all cannot capture the C-terminus of PAK4, but can capture the C-terminus of PAK1. Capturing the C-terminus or not should have nonsignificant effects on the kinase selectivities of FRAX019, FRAX414 and FRAX1036. Whereas, both FRAX597 and G-5555 have weaker interactions with the C-terminus of PAK4 than PAK1, beneficial for their PAK1 selectivities. Additionally, we speculate that switching Part D with a positively charged substituent may help to improve the PAK1 selectivity of the inhibitor.

Conclusion
In the present study, to clarify the inhibition mechanisms and find clues to enhance kinase activity and selec- tivity of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555, MD simulations, MM/PBSA free energy calculations, and MM/GBSA binding energy decomposition analysis and other methods were conducted for the PAK1/FRAX019, PAK4/FRAX019, PAK1/FRAX414, PAK4/FRAX414, PAK1/FRAX597, PAK4/FRAX597,
PAK1/FRAX1036, PAK4/FRAX1036, PAK1/G-5555 and PAK4/G-5555 complexes. The calculated results in-

dicate that the binding free energies of FRAX597, FRAX1036 and G-5555 toward PAK1 are lower than those of FRAX019 and FRAX414, which is in good agreement with the trend of experimental activity data. Moreover, the calculated binding free energies of FRAX414, FRAX597 and G-5555 with PAK4 are 3.11, 6.60 and 4.30 kcal/mol larger than those with PAK1, respectively, showing they all have certain selectivity to PAK1 and the latter two are better than the former, which is also consistent with the experiment. Van der Waals interaction is the major favorable contributor in the binding energy and plays a crucial role for the activity and selectivity.
FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 are all anchored at the hinge area by forming two hydrogen bonds between their pyrido-[2,3-d]pyrimidin-7-one moieties and Leu347PAK1. It is due to this anchoring effect that increasing the length of the inhibitor’s head (Part A) and decreasing its tail (Part E) are good choices for designing efficient inhibitors on the basis of FRAX019. In addition, our calculation predicted that adding appropriate volume and positive charge to Part B should increase the affinity of this part.
Similar to the case with PAK1, these five Afraxis PAK inhibitors are also anchored at the hinge area by forming two hydrogen bonds between their pyrido-[2,3-d]pyrimidin-7-one moieties and Leu398PAK4. Differently, they can capture the C-terminus of PAK1 but cannot capture the C-terminus of PAK4. Capturing the C-terminus or not should have nonsignificant effects on the kinase selectivities of FRAX019, FRAX414 and FRAX1036. Whereas, it has a significant effect on the selectivies of FRAX597 and G-5555. Additionally, we speculate that switching Part D with a positively charged substituent may help to improve the PAK1 selectivity of the inhibitor.
Our work revealed the inhibitory mechanisms of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 toward PAK1 and PAK4. Thus, this work would give valuable information for the future design of specific and efficient PAK inhibitors.

Future perspective
PAKs, involved in many diseases including cancer, are considered as potential therapeutic targets. Numerous molecular and structural biology studies of PAK inhibitors have been carried out in the past decade. Currently, the discovery and development of PAK inhibitors, especially inhibitors with satisfactory kinase selectivity, is at a rapid development stage and it is expected to develop a new generation of antitumor drugs. The inhibitory mechanisms of FRAX019, FRAX414, FRAX597, FRAX1036 and G-5555 toward PAK1 and PAK4 revealed in our work would give valuable information for the future design of specific and efficient PAK inhibitors.

Supplementary data
To view the supplementary data that accompany this paper please visit the journal website at: www.future-science.com/doi/suppl
/10.4155/fmc-2019-0273

Financial & competing interests disclosure
The project was supported by the National Natural Science Foundation of China (grant no. 41530315). The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.
No writing assistance was utilized in the production of this manuscript.

References
Papers of special note have been highlighted as: • of interest; •• of considerable interest
1. Ottilie S, Miller P, Johnson D et al. Fission yeast pak1+ encodes a protein kinase that interacts with Cdc42p and is involved in the control of cell polarity and mating. EMBO J. 14(23), 5908–5919 (1995).
2. Zhao Z-S, Manser E. PAK family kinases: physiological roles and regulation. Cell. Logist. 2(2), 59–68 (2012).

3. Rudolph J, Crawford JJ, Hoeflich KP, Chernoff J. p21-activated kinase inhibitors. Enzymes 34(Pt B), 157–180 (2013).
4. Radu M, Semenova G, Kosoff R, Chernoff J. PAK signalling during the development and progression of cancer. Nat. Rev. Cancer 14(1), 13–25 (2014).
• Mechanisms of p21-activated kinase (PAK) activation in cancer.
5. Ye DZ, Field J. PAK signaling in cancer. Cell. Logist. 2(2), 105–116 (2012).
6. Crawford JJ, Hoeflich KP, Rudolph J. p21-activated kinase inhibitors: a patent review. Expert Opin. Ther. Pat. 22(3), 293–310 (2012).
7. Li X, Wen WH, Liu KD et al. Phosphorylation of Caspase-7 by p21-activated protein kinase (PAK) 2 inhibits chemotherapeutic drug-induced apoptosis of breast cancer cell lines. J. Biol. Chem. 286(25), 22291–22299 (2011).
8. Berger A, Hoelbl-Kovacic A, Bourgeais J et al. PAK-dependent STAT5 serine phosphorylation is required for BCR-ABL-induced leukemogenesis. Leukemia 28(3), 629–641 (2014).
9. Menges CW, Sementino E, Talarchek J et al. Group I p21-activated kinases (PAKs) promote tumor cell proliferation and survival through the AKT1 and Raf–MAPK pathways. Mol. Cancer Res. 10(9), 1178–1188 (2012).
10. Chan PM, Manser E. PAKs in human disease. Prog. Mol. Biol. Transl. Sci. 106, 171–187 (2011).
11. Semenova G, Chernoff J. Role of group I Paks in MPNST cell proliferation, migration and invasion. Presented at: AACR 106th Annual Meeting. PA, USA, 18–22 April 2015.
12. Fryer BH, Field J. Rho, Rac, Pak and angiogenesis: old roles and newly identified responsibilities in endothelial cells. Cancer Lett. 229(1), 13–23 (2005).
13. Rudolph J, Crawford JJ, Hoeflich KP, Wang W. Inhibitors of p21-activated kinases (PAKs) miniperspective. J. Med. Chem. 58(1), 111–129 (2015).
•• A comprehensive review on PAKs inhibitors.
14. Maksimoska J, Feng L, Harms K et al. Targeting large kinase active site with rigid, bulky octahedral ruthenium complexes. J. Am. Chem. Soc. 130(47), 15764–15765 (2008).
15. Blanck S, Maksimoska J, Baumeister J, Harms K, Marmorstein R, Meggers E. The art of filling protein pockets efficiently with octahedral metal complexes. Angew. Chem. Int. Ed. 51(21), 5244–5246 (2012).
16. Tonani R, Bindi S, Fancelli D, Pittala V. Danello MWO2004007504 A1 (2003).
17. Guo C, Mcalpine I, Zhang J et al. Discovery of pyrroloaminopyrazoles as novel PAK inhibitors. J. Med. Chem. 55(10), 4728–4739 (2012).
18. Zhang J, Yang A, Kephart SE et al. US7884117B2 (2011).
19. Guo C, Johnson MC, Li H et al. WO2007/023382A2 (2006).
20. Abdel-Magid AF. PAK1: a therapeutic target for cancer treatment. ACS Med. Chem. Lett. 4(5), 431–432 (2013).
21. Abad-Zapatero C, Metz JT. Ligand efficiency indices as guideposts for drug discovery. Drug Discov. Today 10(7), 464–469 (2005). 22 . Licciulli S, Maksimoska J, Zhou C et al. FRAX597, a small molecule inhibitor of the p21-activated kinases, inhibits tumorigenesis of
neurofibromatosis type 2 (NF2)-associated Schwannomas. J. Biol. Chem. 288(40), 29105–29114 (2013).
•• Report of the x-ray structure of PAK1K299R/FRAX597 complex.
23. Koval AB, Wuest WM. An optimized synthesis of the potent and selective Pak1 inhibitor FRAX-1036. Tetrahedron Lett. 57(3), 449–451 (2016).
24. Mccoull W, Hennessy EJ, Blades K et al. Identification and optimisation of 7-azaindole PAK1 inhibitors with improved potency and kinase selectivity. MedChemComm 5(10), 1533–1539 (2014).
25. Staben ST, Feng JA, Lyle K et al. Back pocket flexibility provides group II p21-activated kinase (PAK) selectivity for type I 1/2 kinase inhibitors. J. Med. Chem. 57(3), 1033–1045 (2014).
26. Aliagas-Martin I, Crawford J, Lee W, Mathieu S, Rudolph J. WO2013026914A1 (2013).
27. Karaman MW, Herrgard S, Treiber DK et al. A quantitative analysis of kinase inhibitor selectivity. Nat. Biotechnol. 26(1), 127–132 (2008).
28. Zhang J, Guo C, Bouzida D et al. US8067591B2 (2011).
29. Murray BW, Guo C, Piraino J et al. Small-molecule p21-activated kinase inhibitor PF-3758309 is a potent inhibitor of oncogenic signaling and tumor growth. Proc. Natl Acad. Sci. USA 107(20), 9446–9451 (2010).
• The first PAK inhibitor to proceed into clinical trials.
30. Rosen L, Blumenkopf T, Breazna A et al. Phase 1, dose-escalateion, safety, pharmacokinetic and pharmacodynamic study of single agent PF-03758309, an oral PAK inhibitor, in patients with advanced solid tumors. Presented at: 22nd EORTC-NCI-AACR Symposium on Molecular Targets and Cancer Therapeutics. PA, USA, 16–19 November 2010.
31 . Ndubaku CO, Crawford JJ, Drobnick J et al. Design of selective PAK1 inhibitor G-5555: improving properties by employing an unorthodox low-p K a polar moiety. ACS Med. Chem. Lett. 6(12), 1241–1246 (2015).

•• Report of the cocrystal structures of PAK1K299R/FRAX1036 and PAK1D389N,T423E/G-5555.
32. Mortazavi F, Lu J, Phan R et al. Significance of KRAS/PAK1/Crk pathway in non-small-cell lung cancer oncogenesis. BMC Cancer
15(1), 381 (2015).
33. Chow H-Y, Dong B, Duron SG et al. Group I Paks as therapeutic targets in NF2-deficient meningioma. Oncotarget 6(4), 1981–1994 (2015).
34. Beauchamp RL, James MF, Desouza PA et al. A high-throughput kinome screen reveals serum/glucocorticoid-regulated kinase 1 as a therapeutic target for NF2-deficient meningiomas. Oncotarget 6(19), 16981–16997 (2015).
35. Nuche-Berenguer B, Ramos-A´lvarez I, Jensen R. The p21-activated kinase, PAK2, is important in the activation of numerous pancreatic acinar cell signaling cascades and in the onset of early pancreatitis events. Biochim. Biophys. Acta Mol. Basis. Dis. 1862(6), 1122–1136 (2016).
36 . Yeo D, He H, Patel O, Lowy AM, Baldwin GS, Nikfarjam M. FRAX597, a PAK1 inhibitor, synergistically reduces pancreatic cancer growth when combined with gemcitabine. BMC Cancer 16(1), 24 (2016).
• Application of FRAX597 in reducing the growth of pancreatic cancer.
37. Ruiz R, Jahid S, Harris M et al. The RhoJ-BAD signaling network: an Achilles’ heel for BRAF mutant melanomas. PLoS Genet. 13(7), e1006913 (2017).
38. Tan Y, Sementino E, Chernoff J, Testa JR. Targeting MYC sensitizes malignant mesothelioma cells to PAK blockage-induced cytotoxicity. Am. J. Cancer Res. 7(8), 1724–1737 (2017).
39. Chang JK, Ni Y, Han L et al. Protein kinase D1 (PKD1) phosphorylation on Ser203 by type I p21-activated kinase (PAK) regulates PKD1 localization. J. Biol. Chem. 292(23), 9523–9539 (2017).
40. Yeo D, Phillips P, Baldwin GS, He H, Nikfarjam M. Inhibition of group 1p21-activated kinases suppresses pancreatic stellate cell activation and increases survival of mice with pancreatic cancer. Int. J. Cancer 140(9), (2017).
41. Lyons J, Brubaker DK, Ghazi PC et al. Integrated in vivo multiomics analysis identifies p21-activated kinase signaling as a driver of colitis. Sci. Signal. 11(519), eaan3580 (2018).
42. Prudnikova TY, Chernoff J. The group I Pak inhibitor Frax-1036 sensitizes 11q13-amplified ovarian cancer cells to the cytotoxic effects of Rottlerin. Small GTPases 8(4), 193–198 (2017).
43. Araiza-Olivera D, Feng Y, Semenova G, Prudnikova TY, Rhodes J, Chernoff J. Suppression of RAC1-driven malignant melanoma by group A PAK inhibitors. Oncogene 37(7), 944–952 (2017).
44. Rudolph J, Murray LJ, Ndubaku CO et al. Chemically diverse group I p21-activated kinase (PAK) inhibitors impart acute cardiovascular toxicity with a narrow therapeutic window. J. Med. Chem. 59(11), 5520–5541 (2016).
45 . Ong CC, Gierke S, Pitt C et al. Small molecule inhibition of group Ip21-activated kinases in breast cancer induces apoptosis and potentiates the activity of microtubule stabilizing agents. Breast Cancer Res. 17(1), 59 (2015).
• Application of FRAX1036 in the detection of PAK1 loss of function.
46 . Zhang EY, Ha BH, Boggon TJ. PAK4 crystal structures suggest unusual kinase conformational movements. Biochim. Biophys. Acta Proteins Proteom. 1866(2), 356–365 (2018).
•• First report of the cocrystal structures of Afraxis inhibitor (FRAX486) with PAK4.
47. Berman HM, Battistuz T, Bhat T et al. The protein data bank. Acta Crystallogr. Sect. D. Biol. Crystallogr. 58(6), 899–907 (2002).
48. Sybyl-X 2.0. Tripos International, MO, USA (2012).
49. Frisch MJ, Trucks GW, Schlegel HB et al. Gaussian 09. Revision A.01; Gaussian, Inc., CT, USA (2009).
50. Bayly CI, Cieplak P, Cornell W, Kollman PA. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 97(40), 10269–10280 (1993).
51. Case DA, Betz RM, Cerutti DS et al. AMBER 2016. University of California, CA, USA (2016).
52. Homeyer N, Horn AHC, Lanig H, Sticht H. AMBER force-field parameters for phosphorylated amino acids in different protonation states: phosphoserine, phosphothreonine, phosphotyrosine, and phosphohistidine. J. Mol. Model. 12(3), 281–289 (2006).
53. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79(2), 926–935 (1983).
54. Maier JA, Martinez C, Kasavajhala K, Wickstrom L, Hauser KE, Simmerling C. ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11(8), 3696–3713 (2015).
55. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. Development and testing of a general amber force field. J. Comput. Chem.
25(9), 1157–1174 (2004).
56. Darden T, York D, Pedersen L. Particle mesh Ewald: an N· log (N) method for Ewald sums in large systems. J. Chem. Phys. 98(12), 10089–10092 (1993).
57. Ryckaert J-P, Ciccotti G, Berendsen HJ. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23(3), 327–341 (1977).

58. Kollman PA, Massova I, Reyes C et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33(12), 889–897 (2000).
59 . Rocchia W, Alexov E, Honig B. Extending the applicability of the nonlinear Poisson−Boltzmann equation: multiple dielectric constants and multivalent ions. J. Phys. Chem. B 105(28), 6507–6514 (2001).
60. Weiser J, Shenkin PS, Still WC. Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem.
20(2), 217–230 (1999).
61. Hou T, Yu R. Molecular dynamics and free energy studies on the wild-type and double mutant HIV-1 protease complexed with amprenavir and two amprenavir-related inhibitors: mechanism for binding and drug resistance. J. Med. Chem. 50(6), 1177–1188 (2007).
62. Wang J, Hou T, Xu X. Recent advances in free energy calculations with a combination of molecular mechanics and continuum models.
Curr. Comput. Aid Drug 2(3), 287–306 (2006).
63. Onufriev A, Bashford D, Case DA. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 55(2), 383–394 (2004).
64. Gohlke H, Kiel C, Case DA. Insights into protein–protein binding by binding free energy calculation and free energy decomposition for the Ras–Raf and Ras–RalGDS complexes. J. Mol. Biol. 330(4), 891–913 (2003).
65. The PyMOL Molecular Graphics System. Version 2.0 Schro¨dinger, LLC (2002). http://pymol.org
66. Lei M, Lu WG, Meng WY et al. Structure of PAK1 in an autoinhibited conformation reveals a multistage activation switch. Cell 102(3), 387–397 (2000).