Exploring aspartic protease inhibitor binding to design selective antimalarials

Selectivity is a major issue in the development of drugs targeting pathogen aspartic proteases. Here we explore the selectivity determining factors by studying specifically designed malaria aspartic protease (plasmepsin) open-flap inhibitors. Metadynamics simulations are used to uncover the complex binding/unbinding pathways of these inhibitors, and describe the critical transition states in atomistic resolution. The simulation results are compared to experimentally determined enzymatic activities. Our findings demonstrate that plasmepsin inhibitor selectivity can be achieved by targeting the flap loop with hydrophobic substituents that enable ligand binding under the flap loop, as such behaviour is not observed for several other aspartic proteases. The ability to estimate compound selectivity before they are synthesized is of great importance in drug design, therefore, we expect that our approach will be useful in selective inhibitor design not only against aspartic proteases, but other enzyme classes as well.


Introduction
Aspartic proteases are wide spread in eukaryotes, and they are attracting attention as drug targets [1,2] . This family of enzymes uses an activated water molecule bound to the catalytic dyad -two aspartic acid residues -to attack the scissile peptide bond in a substrate. The aspartic protease catalysed peptide bond cleavage can be prevented by inhibitors that replace this catalytic water molecule or in which the scissile peptide bond is replaced by a non-cleavable moiety. There are numerous inhibitors reported in literature targeting aspartic proteases involved in HIV/AIDS, malaria, Alzheimer's disease, etc [3,4] . Over the recent years, aspartic proteases from the malaria parasite -plasmepsins (plms) -have gained particular interest from medicinal chemists [5][6][7][8][9][10][11][12][13] due to increasing resistance of known antimalarials [14,15] . An important issue when developing pathogen aspartic protease inhibitors is their selectivity over human enzymes. Plms share high sequence similarity with several human aspartic proteases, the most similar being the lysosomal enzyme cathepsin D (catD). Therefore, it is commonly used as a marker for cross-inhibition. Plm II and catD share 35% sequence identity, and even higher similarity is observed in the active site. Due to the high protein structural similarity, design of potent plasmepsin inhibitors that do not inhibit catD remains a challenging task. Recent studies [16] have shown that plm inhibitor selectivity over human aspartic proteases can be achieved by targeting the flap loop region of plms. Flap loop is a single long β-hairpin structure that lies perpendicularly over the aspartic dyad (see Figure 1.A), and is believed to be involved in the substrate recognition [17][18][19][20] . There have been numerous, mostly computational, studies focused on the characterisation of aspartic protease flap loop dynamics [18,19,[21][22][23][24][25][26][27][28][29][30][31][32] . It is commonly acknowledged [25] , that the flap-loop of aspartic proteases occupies a semi-open conformation in apo state. However, the hinge like motion of the flap loop is also capable of exposing a hydrophobic pocket (see Figure 1.B) under the flap loop producing so called open-flap conformation [33][34][35][36][37] . Inhibitors that bind partly under the flap loop (also known as open-flap inhibitors) have improved selectivity over human aspartic proteases, and the selectivity gain is attributed to differences in enzyme flexibility -plms easily reach openflap conformation where an additional hydrophobic pocket is formed, whereas cathepsins either do not have or do not easily reach an open-flap conformation. The precise mechanism of enzyme reorganisation for the inhibitor binding to open-flap conformation, however, is not yet known. It is unclear what are the protein structural features and key protein-inhibitor interactions that enable flap loop opening, and lead to inhibitor binding under the flap loop. Such information would allow us to design inhibitors with improved selectivity not only against plms, but other aspartic proteases, and, possibly, other enzymes with highly mobile loops near the binding site. An accurate description of the flap loop motion upon ligand binding at a molecular scale is essential to understand selectivity determining factors of open-flap inhibitors. Currently, however, there are no experimental methods that would enable such dynamic process characterisation at a molecular scale. Lately, in-silico techniques have emerged as a viable alternative to experimental methods and are becoming increasingly effective in providing atomic-level description of protein-ligand interactions [38][39][40][41][42][43][44][45][46][47] . Dynamic processes like these are typically studied using molecular dynamics (MD) simulations, however, ligand binding/unbinding is a rare event that can be hard to sample within a reasonable time [48,49] . Therefore, enhanced sampling methods like metadynamics (MetaD) are utilised to accelerate the occurrence of rare events by depositing a history dependent biasing potential to selected degrees of freedom (also known as collective variables (CVs)) [39,40,[50][51][52][53][54][55][56] . The bias potential within CV sub-space is built up by adding Gaussians along the sampled trajectory, and bias potential applied encourages the system to visit higher energy states and discourage the system from revisiting configurations already sampled. This approach leads to a significantly improved exploration of the simulation space, and, once the exploration is complete, an efficient reconstruction of the free energy surface (FES) as a function of the CVs used [57,58] . MetaD in its various modifications has been applied in numerous ligand-protein binding studies [38][39][40][41][42][43][44][45][46][47] . One of these MetaD variations is funnel metadynamics (FMetaD), where a funnel shaped restraint potential is applied to the system to reduce the space explored by the ligand once it is unbound, thus reducing computational cost even more [46,[59][60][61][62][63] . In FMetaD, the bias potential is typically applied to the distance between the ligand and the binding site, thus, sampling multiple binding/unbinding events within simulation. FMetaD allows identification of the ligand binding mode, absolute protein-ligand binding free energy, insight into the dynamics of the ligand binding, including the presence of alternative ligand binding modes and other energetically relevant states (e.g., transition states). Furthermore, it has been demonstrated that FMetaD can be used to study ligand binding process in challenging systems where binding site is buried or protein undergoes significant conformational changes upon binding [60,64] . FMetaD has demonstrated to be able to predict ligand binding free energy with <1.5 kcal/mol accuracy [60,61] . It is important to note, however, that FMetaD is an MDbased protocol and the correctness of the results depend on the accuracy of the protein and ligand force-field parameters, as well as system amino acid and ligand protonation state selection. By examining known non-peptidomimetic plm inhibitors, we previously hypothesized [16] , that one of the key features enabling the ligand binding to the enzyme in open flap conformation might be hydrophobic substituent capable of aromatic interactions with Trp41, Tyr77 and/or Phe111 (plm II numbering) under the flap loop. To test this hypothesis, we designed and synthesized a series of 2-aminoquinazolin-4(3H)-one based [65][66][67] inhibitors (see Figure 2) with hydrophobic substituents bearing benzene group at various locations, determined their activity against plm II, IV, V and cat D, and explored their binding pathway using funnel metadynamics simulations. This showed that the inhibitor binding mode can be estimated using metadynamics simulations, and information obtained can be used in selective inhibitor design.  IDs 2BJU and 4Z22, respectively). Residues W41 and Y77, believed to play crucial role in flap loop dynamics, and catalytic dyad (D34 and D214) are shown as sticks (plm II numbering). Note that W41 changes its conformation upon flap opening/closing. B Hydrophobic flap pocket exposed upon flap loop opening (PDB ID 4Z22). Ligand is shown as green sticks, hydrophobic flap pocket residues as white sticks. Flap loop is in salmon. Hydrogens are omitted for clarity.

Compound design
The 2-aminoquinazolin-4(3H)-one derivatives were designed to have linear, hydrophobic flap pocket substituents of various alkyl chain length, benzene ring position and benzene ring count (see Figure 2). The initial compound subset varied in length of the linear alkyl linker between the ligand core and terminal benzene ring (2)(3)(4)(5)(6)(7)(8)(9). Once the optimal flap pocket substituent length was identified, another ligand subset, where the substituent length was kept constant while the position of the benzene ring was varied, was synthesized (11)(12)(13)(14)(15). To evaluate the contribution of the terminal benzene ring on the inhibitor activity, compound 10, bearing linear alkyl chain of the same length as the most active compound was synthesized. Moreover, several compounds bearing two benzene rings at various positions were synthesized to evaluate the impact of additional benzene ring in the flap pocket substituent (16)(17)(18)(19). For comparison purposes, compound without flap pocket substituent (1) -the original fragment hit identified in our screening campaign [67] -was prepared and characterised as well.

Chemistry
The synthesis of compounds 1, 2, 10, 11, 16 and 21 has been previously described [65,67] . Compounds 3-9, 12-15 and 17-20 were synthesized as shown in Scheme 1. Building block 21 was subjected to Sonogashira reaction with alkynes 24a-e, followed by the saturation of alkyne and the cleavage of benzoyl group to obtain products 4, 6, 7, 13 and 18. Suzuki-Miyaura reaction of 21 with the corresponding vinyl dioxaborolanes 25ac, or aryl or alkyl dioxaborolanes 26a-c with subsequent saturation of unsaturated bonds and the cleavage of benzoyl group gave products 5, 14, 19 and 3, 12, 20. The target compound 17 was synthesized from pinacol borane derivative 22 and benzylbromide 27 using Suzuki-Miyaura type reaction with subsequent hydrogenation of alkyne moiety and benzoyl group cleavage. Coupling of 21 with propargyl alcohol (28) followed by oxidation of alcohol to aldehyde and reduction of alkyne to alkene provided 23. Building block 23 was further converted into products 8, 9 and 15 using the Wittig reaction, the saturation of double bonds and the cleavage of benzoyl group (see Supporting Information for details).

Enzymatic activity
Compounds synthesized were tested for their inhibitory activity against digestive vacuole plms II and IV, endoplasmic reticulum protease plm V, and selectivity cross-marker -human catD (see SI for experimental details). The plm II inhibition potency (IC50) of the compounds tested varied in a range from 3.5 µM down to 64 nM depending on the length of hydrophobic substituent and position of aromatic ring (see Figure 2 and Table SI1 for full inhibition data). There were minor variations between plm II and plm IV inhibition potencies, whereas the potency against plm V and human catD was considerably lower. The catD inhibition potency followed plms II and IV inhibition pattern with 10 to 70 fold lower activities (55 to 2 µM), whereas potency against plm V varied in a narrow range (90 to 25 µM) and size of the flap pocket substituent seems to play no role in inhibitor activity against plm V. More detailed analysis of structure activity relationships (SAR) of the inhibitor series indicates that installation of the flap pocket substituent greatly increases the inhibitor potency (15-fold activity increase from compound 1 to 2). By increasing the distance between the ligand core (2-aminoquinazolin-4(3H)-one group) and the terminal aromatic ring, the inhibitory activity initially decreases, reaching lowest activity for inhibitors with methylene and ethylene linker (3,4), whereas further increase of linker length led to increase in inhibition potency. The highest potency was achieved for the compound bearing phenylpentyl group (7). Linker length increase beyond pentylene group led to reduced activity (8,9). The inhibition potency of the compound 10 suggests that presence of the terminal aromatic group in the flap pocket substituent is not necessary for high activity, as inhibition potency of both compounds (7 and 10) was comparable. A ligand subset (11)(12)(13)(14)(15) bearing the linkers with the same overall length as compound 7, but varied in positions of aromatic ring was tested in order to understand, whether the activity of the inhibitor is affected solely by the length of the hydrophobic flap pocket substituent or position of the aromatic ring. The general inhibitory activity trend was the same as in initial inhibitor set, however, these inhibitors were 5 to 20 times more potent. The last ligand subset had terminal benzene ring fixed at a specific distance, and an additional aromatic ring was installed in the flap pocket substituent to introduce stacking interactions with Trp41, Tyr77, Phe111, Tyr115 and/or Phe120. Inhibitors bearing two benzene rings in the flap pocket substituent showed low micromolar/nanomolar potency, however, no clear SAR was observed for these compounds.

Metadynamics simulations
Schematic representation of the funnel shaped restraint potential used in our FMetaD simulations is given in Figure 3.A. The restraint potential was placed in a way that does not interfere with ligand binding and unbinding events, and its dimensions were chosen as small as possible to accelerate convergence, but large enough to avoid artefacts in the sampling. The funnel defining parameters are summarised in Figure 3.A. Here, the CV controlling the ligand binding/unbinding was the distance (d1) between the centre of mass of ligand transition state mimetic atoms and the centre of mass of the aspartic dyad Asp34 and Asp214 Cγ atoms (see Figure 3.B,C). Distance CV easily discriminates between bound and unbound states, however, trajectory obtained can be reweighed to remap the free energy as a function of alternative CVs and gain more detailed insights, revealing bound minima, which are overlying in the maps corresponding to the initial CV choice [58,68] . Several unbiased and exploratory MetaD simulations were first performed to get an idea of the preferred path of unbinding with statistical significance. During our preliminary simulations we observed that flap loop does not close completely upon ligand unbinding, and rotation of the Trp41, which is known to change conformation when bound to open-flap inhibitors [37,67] (see Figure 1.A), was not observed. To overcome this obstacle, we introduced PathCV [69] that allows to study transitions between two or more states. Here we constructed the flap loop opening/closing pathway (including Trp41 rotation) by running trial Path MetaD simulations using X-ray crystal structures of closed and open flap conformations [67,70] as end states (PDB IDs 1LF4 and 4Z22, respectively; see Figure 3.D,E). Equally spaced loop opening/closing pathway was constructed from the preliminary simulations and was further used in ligand binding/unbinding studies. Additionally, soft harmonic restraining walls were applied at both ends of the funnel axis to limit the exploration of the ligand inside the region enclosed by the funnel potential. The overall structure of the protein was preserved throughout simulations, during which the loops and the binding site residues showed the necessary flexibility to allow binding and unbinding of the open-flap inhibitors. Binding to plm II. The performed metadynamics simulations directly provide an estimate of where the system spends most of its time [52,55] . The FES computed as a function of distance d1 and torsion φ (ligand alignment with respect to flap pocket axis, see Figure SI1) indicate that there are several stabilised transition states along the ligand binding pathway (see Figure 4). Most of these states are observed for all of the plm inhibitors synthesized and tested. The deepest basin A in all simulations corresponds to the state where inhibitor is bound in a pose analogue to that found in crystal structure (PDB ID: 4Z22). The 2-aminoquinazolin-4(3H)one core that mimics the transition state is involved in hydrogen bonding with Asp34 and Asp214, adjacent ligand core aromatic ring is interacting with the flap loop Tyr77 sidechain, whereas flap pocket substituent is buried under the flap loop. The dominant interactions under the flap loop depend on the flap pocket substituent length, with longest substituents interacting with Val105, Phe111 and/or Tyr115 sidechains deeper in the flap pocket, while shorter inhibitor substituents are interacting with Trp41, Phe111 and/or Phe120 sidechains. In all cases, ligand binding have disrupted the hydrogen bonding between Trp41 and Tyr77, which is common for aspartic proteases in closed flap conformation [16,71] . Once this interaction is not present, enzyme flap loop is able to open-up and Trp41 is able to rotate to expose the hydrophobic depth of the hydrophobic pocket. This behaviour is observed upon binding of all inhibitors, even those with short flap pocket substituents. It is worth noting that Tyr77 stacking interactions with inhibitor core in basin A were observed on both sides of the inhibitor (e.g. inhibitor locked in-between Tyr77 and Phe111, or inhibitor interacting with Tyr77 on one side and solvent exposed on the other). Basin B corresponds to the binding mode where ligand core (2-aminoquinazolin-4(3H)-one group) is occupying pose similar to that of crystal structure, while the flap pocket substituent is not buried under the flap loop. Flap pocket substituent here is interacting with Tyr77 and Trp41, however, it is situated on the other side of these residues and is partly solvent exposed. Since this binding mode is observed when the flap loop is in semi-open conformation, we believe that this might be one of the transition states in the ligand binding and flap loop opening pathway. More specifically, the flap loop opening is initiated by the inhibitor binding in the substrate binding groove, and is followed by the inhibitor rearrangement were hydrophobic flap pocket substituent is guided underneath the flap loop while maintaining interactions with aspartic dyad. During this flap opening motion the common aspartic proteases hydrogen bonding network is disrupted and Trp41 rotates to expose the terminal end of the hydrophobic flap pocket and engage in aromatic stacking interactions with flap pocket substituent. Besides previously described poses, where ligand core is interacting with catalytic dyad, two additional binding modes were identified for majority of the inhibitors. One, where aromatic group of the hydrophobic flap pocket substituent is wedged in-between Tyr77 and Phe111/Tyr115 while ligand core is solvent exposed (basin C in Figure 4), and another one, where the inhibitor is interacting with aromatic residues Phe11, Tyr17 in the binding grove (basin D in Figure 4). Ligand-protein interactions observed in basin D take place ~2 nm away from the aspartic dyad, therefore, it is unlikely that they play crucial role in inhibitor recognition and inhibitor binding, however, these intermediate states might play role in ligand guiding to the binding site. Basin C corresponds to the state where inhibitor is situated closer to the binding site, however, does not interact directly with the aspartic dyad. Interactions here are dominated by the hydrophobic stacking between the inhibitor terminal benzene ring and Tyr77, Phe111 and/or Tyr115 sidechains in flap loop and α helix right next to the loop. These interactions are mostly observed when the flap loop is in closed conformation, and just like complex in basin B, could be a transition state on the flap loop opening pathway. Here flap loop opening is initiated by the hydrophobic interactions that direct the flap pocket substituent under the flap loop by successively interacting with aromatic residues Phe120, Tyr115, Phe111, Tyr77 and finally Trp41, while the transition state mimetic remains solvent exposed. The presence of alternative binding modes observed in FMetaD simulations is supported by the X-ray diffraction data. Attempts were made to co-crystallise plm II and compounds synthesized (see SI for details), however, resolution for the majority of the crystals obtained were suboptimal. The only crystal that was solved was with compound 7 (PDB ID: 7QYH, to be released upon publication). Due to poorly defined electron density in the plm II binding site, it was not possible to solve the position of compound 7, however, unmodelled electron density in the binding site indicates ambiguous flap pocket substituent binding (see Figure  5). Electron density map suggest that most probable phenlypentyl substituent placement is under the flap loop on either side of the Trp41, which matches FES basins A and B observed in simulations. The states described above are observed for all inhibitors, however, there are some notable differences between inhibitors. The FES of the 2-aminoquinazolin-4(3H)-one without flap pocket substituent (1) is the smoothest one, indicating that the inhibitor binding energy is low and energetic penalty associated with the transition from one binding mode to another one is also low (~2 kcal/mol). Several energetically similar binding modes with negligible transition energies mean that the inhibitor is highly mobile and fluctuates between multiple competing binding modes, and, as a result, its inhibition potency is low. Introduction of more extensive flap pocket substituent reduces the accessible binding space and shifts the binding preference to the binding mode observed in the crystal structure. This results in improved inhibitor potency, since role of the alternative binding modes is diminished. This is in line with experimentally determined activity data, where inhibition potency increases as flap pocket substituent length increases, however, the impact of the flap pocket length increase is not linear. We believe that lower activity for compounds bearing phenylmethyl-and phenylethylsubstituents (2, 3) is a result of sub-optimal terminal benzene ring placement in the flap pocket due to competitive stacking interactions. More specifically, shorter linker length prevents simultaneous ligand core stacking interactions with Tyr77 and/or Phe111, and terminal benzene ring interactions with Trp41 and/or Tyr115. Thus, compounds bearing methylene and ethylene linkers are able to form optimal interactions either with catalytic dyad and Tyr77/Phe111, or with Tpr41/Tyr115, but not both. In addition, inhibitors with shorter flap pocket substituents bind to the open flap conformation, however, they are not able to fill the pocket completely, thus, leading to partly solvated hydrophobic pocket and/or introducing thermodynamic penalty due to partly filled cavity. The increase of the linker length enables simultaneous interactions for both groups, reaching optimal positioning with pentylene linker (comp. 7). Moreover, the phenylpentyl substituent size is optimal to fill the flap pocket completely. Flap pocket substituents bearing longer alkyl chain linker can fit in the cavity (8,9), however, substituents that are too long in their fully extended conformation are expected to adopt energetically unfavourable conformations with internal gauche strains, leading to reduced potency [72] . The potency of the inhibitors with fixed flap pocket length and various benzene ring position (11)(12)(13)(14)(15)7) indicate that the flap pocket depth filling improves the inhibitor activity to sub-micromolar level. The fact that SAR with respect to benzene ring position for this inhibitor subset was identical to initial inhibitor subset indicate that aromatic stacking and/or steric effects provided by the benzene ring play the same, decisive role in inhibitor activity. The additional terminal n-alkyl substituent for this inhibitor series improves the overall activity with respect to corresponding shorter substituents by filling the depth of the flap pocket, thus, reducing the thermodynamic penalty due to unfilled and/or partly solvated hydrophobic pocket [72] .  [58] as a function of the ligand core-catalytic dyad distance d1 (see Figure 3.B) and a torsion φ (ligand alignment with respect to flap pocket axis, see Figure SI1). Isosurfaces are shown every 1 kcal/mol. Deepest FES basins are indicated as A, B, C and D, and corresponding binding modes are shown on the right. Binding mode in basin A matches the pose observed in the crystal structure 4Z22, whereas binding modes B, C and D are transition states on the ligand binding pathway. The binding modes estimated from the un-modelled electron density in the plm II-compound 7 co-crystal are indicated with black dots. Ligand is shown as green sticks, selected binding pocket residues and catalytic dyad residues as white sticks. Flap loop is in salmon. Hydrogens are omitted for clarity. Poor activity against catD and plm V. The FES of 2-aminoquinazolin-4(3H)-ones binding to plm V ( Figure  6 top) and catD (Figure 6, bottom) are slightly different from those of plm II binding. The majority of the binding modes observed during the inhibitor binding to plm II are present also on the catD and plm V binding pathway, however, there is one noticeable difference -binding under the flap loop is negligible. For both enzymes, the energetically favoured binding mode is similar to that observed in basin B for plm II, where transition state mimetic group is interacting with aspartic dyad while the flap pocket substituent is located between the flap loop and unstructured region right next to it (see compound 7 binding modes B in Figure 6). Here hydrophobic flap pocket substituent is interacting with Phe74, Ile76 and/or Val144 of catD, and Tyr286, Tyr288 and/or Leu311 of plm V. Another, similarly favoured binding mode corresponds to the inhibitor bound in a binding groove in the opposite direction, where hydrophobic part of the inhibitor interacts with Phe126 and/or Phe131 of catD, and Leu179 and/or Tyr182 of plm V (see Figure 6, binding modes A). Despite plm V highly flexible flap loop tip (since Tyr139 is not involved in direct hydrogen bonding with protein core β sheet motif amino acids like plm II and catD, see Figure 7) inhibitor binding in a pose where flap pocket substituent is fully buried under the flap loop was not observed. Minor interactions under the plm V flap loop tip (basin C) were similar to those observed in basin C for plm II, where hydrophobic flap pocket substituent was stacked in-between Tyr139 and Phe180 while the transition state mimetic was solvent exposed. The rigidity at the plm V flap loop hinge region does not allow further flap opening, and as such inhibitor remains in non-optimal conformation. Inhibitor interactions under the catD flap loop tip were similar -terminal benzene ring interacted with Tyr78 and Phe126, however it did not lead to inhibitor binding under the flap loop. Similar binding behaviour was observed for all inhibitors. Taking FMetaD simulation data into account, it is safe to assume, that poor plm V and catD inhibition is a result of inhibitor inability to bind under the flap loop, and enzyme inhibition is achieved by inhibitor binding in an alternative pose. These alternative binding modes, however, lack optimal (coplanar) 2-aminoquinazolin-4(3H)-one core placement (due to flap pocket substituent clashing with the plm V and catD flap loops) and, thus, ligand activity is reduced. We believe that structural feature that plays crucial role in inhibitor binding under the flap loop is hydrophobic amino acid region at the flap pocket "entrance" between the flap loop and adjacent α helix. All enzymes have phenylalanine situated at the same place right next to the flap loop tip (Phe111 for plm II, Phe180 for plm V and Phe126 for catD, see Figure 8), which in unbound state can interact with flap loop tip Tyr and approaching inhibitor. Such transition states were observed in our simulations, and correspond to previously described states in basin C. Once initial interactions are established, the hydrophobic inhibitor substituent is then guided deeper in the flap pocket to interact with Tyr115 and Trp41 of plm II. In plm II, the Tyr115 sidechain is rigidly locked in the position by hydrogen bonding with Asp121 main chain, whereas Trp41 is known to change the conformation upon non-peptidomimetic inhibitor binding/flap opening. The flap pocket residue rearrangement results in a conformation where hydrophobic flap pocket substituent is locked between Tyr115 and Trp41. The hydrophobic flap pocket remains un-solvated during the flap pocket filling step, and required flap opening is aided by gradual flap pocket filling with hydrophobic inhibitor tail. In plm V and catD, the position of the Tyr115 sidechain is occupied by Glu176 and Gln121, respectively, and, therefore, hydrophobic interactions with the inhibitor are discouraged. This prevents ligand flap pocket substituent guiding deeper under the flap loop and efficient inhibitor binding. Figure 6. The FES of compound 7 binding to plm V (top) and catD (bottom) obtained using FMetaD simulations and reweighed [58] as a function of the ligand core-catalytic dyad distance d1 and a torsion φ. Isosurfaces are shown every 1 kcal/mol. Deepest FES basins for each system are indicated as A and B, and corresponding binding modes are shown on the right. Ligand is shown as green sticks, selected binding pocket residues and catalytic dyad residues as white sticks. Flap loop is in salmon. Hydrogens are omitted for clarity.  . Superimposed plm II (white, PDB ID 1SME), plm V (pink, PDB ID 4ZL4) and catD (purple, PDB ID 1LYA) binding sites. Protein backbone is depicted in cartoon representation, selected binding pocket residues are shown as sticks. Hydrogens are omitted for clarity.

Conclusions
The characterisation of the malaria aspartic protease flap loop dynamics upon open-flap inhibitor binding provided several conclusions. First, several alternative binding modes are observed for majority of 2aminoquinazolin-4(3H)-one based plm inhibitors. These binding modes are believed to be key transition states on the ligand binding pathway, and plays role in the ligand flap pocket substituent guiding under the flap loop and flap loop opening. Second, introduction of larger flap pocket substituent reduces the conformational space available by the ligand and shifts the binding preference to the binding mode observed in the crystal structure. This results in improved inhibitor activity, since role of the alternative binding modes is diminished. Third, poor 2-aminoquinazolin-4(3H)-one based inhibitor potency against plm V and catD is a result of ligands inability to bind under the flap loop. Binding modes identified suggest that the ligand binding under the flap loop in plm II is initiated by the hydrophobic interactions in-between the flap loop and adjacent α-helix; such interactions in plm V and catD are unfavourable. Taken together, our findings confirm that plm inhibitor selectivity can be achieved by targeting the flap loop with hydrophobic substituents that enables ligand binding under the flap loop, as such behaviour is not observed for several other aspartic proteases. The insight on the molecular determinants of ligand selectivity between plm II, plm V and catD provide a clear indication for the structure-based rational design of selective inhibitors. The ability to estimate compound selectivity before they are synthesized is of great importance in drug design, therefore, we expect that approach used here will be a valid support in selective inhibitor design not only against aspartic proteases, but other enzymes as well.