Applying multi-state modeling using AlphaFold2 for kinases and its application for ensemble screening

1 1 Applying multi-state modeling using AlphaFold2 for 2 kinases and its application for ensemble screening 3 4 Jinung Song1¶, Junsu Ha2¶, Juyong Lee1,2,3, Junsu Ko2, Woong-Hee Shin2,4* 5 1 College of Pharmacy, Seoul National University, Seoul, Republic of Korea 6 2 Arontier Co., Seoul, Republic of Korea 7 3 Department of Molecular Medicine and Biopharmaceutical Sciences, Graduate School of 8 Convergence Science and Technology, Seoul, Republic of Korea 9 4 Department of Biomedical Informatics, Korea University College of Medicine, Seoul, Republic 10 of Korea 11 12 * Corresponding author 13 E-mail: [email protected] (WHS) 14 15 ¶: These authors contributed equally to this work. 16 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 2 17 Abstract 18 Structure-based virtual screening (SBVS) is a pivotal computational approach in drug 19 discovery, enabling the identification of potential drug candidates within vast chemical libraries 20 by predicting their interactions with target proteins. The SBVS relies on the receptor protein 21 structures, making it sensitive to structural variations. Kinase, one of the major drug targets, is 22 known as one of the typical examples of an active site conformation change caused by the type 23 of binding inhibitors. Examination of human kinase structures shows that the majority of 24 conformations have the DFGin state. Thus, SBVS using the structures might cause a favor of 25 type of ligand type I inhibitors, bind to the DFGin state, rather than finding the diverse scaffolds. 26 Recent advances in protein structure prediction, such as AlphaFold2 (AF2), offer promising 27 solutions but may still be possibly influenced by the structural bias in existing templates. 28 To address these challenges, we introduce a multi-state modeling (MSM) protocol for kinase 29 structures. We apply MSM to AF2 by providing state-specific templates, allowing us to 30 overcome structural biases and thus apply them to kinase SBVS. We benchmarked our MSM 31 models in three categories: quality of predicted models, reproducibility of ligand binding poses, 32 and identification of hit compounds by ensemble SBVS. The results demonstrate that MSM- 33 generated models exhibit comparable or improved structural accuracy compared to standard AF2 34 models. We also show that MSM models enhance the accuracy of cognate docking, effectively 35 capturing the interactions between kinases and their ligands. 36 In virtual screening experiments using DUD-E compound libraries, our MSM approach 37 consistently outperforms standard AF2 modeling. Notably, MSM-based ensemble screening 38 excels in identifying diverse hit compounds for kinases with structurally diverse active sites, .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 3 39 surpassing standard AF2 models. We highlight the potential of MSM in broadening the scope of 40 kinase inhibitor discovery by facilitating the identification of chemically diverse inhibitors. 41 42 Author Summary 43 One of the main problems with structure-based virtual screening is structural flexibility. 44 Ensemble screening is one of the conventional approaches to solving the issue. Gathering 45 experimental structures or molecular simulations could be used to compile the receptor 46 structures. Recent developments in algorithms for predicting protein structures, like AlphaFold2, 47 suggest that different receptor conformations could be produced. However, the prediction 48 approaches produce biased structures because of the bias in the structure database. In order to 49 solve the problem, we developed a protocol called multi-state modeling for kinases. Rather than 50 supplying multiple sequence alignments as an input, we gave the AlphaFold2 a specific template 51 structure and the sequence alignment between the template and query. 52 Our findings imply that our technique can yield a particular structural state of interest with an 53 enhanced or comparable structural quality to AlphaFold2 and predict highly accurate protein- 54 ligand complex structures. Lastly, compared to the typical AlphaFold2 models, ensemble 55 screening using the multi-state modeling approach improves the structure-based virtual screening 56 performance, particularly for diverse active molecular scaffolds. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 4 57 Introduction 58 Structure-based virtual screening (SBVS) is one of the most widely used computational 59 drug discovery approaches to identify novel active compounds against the target from a virtual 60 molecule library. By predicting the interaction between the ligand and the target protein, the 61 technique ranks the ligands based on their scores that mimic the binding affinity between the 62 molecules. The SBVS is a cost- and time-efficient way to explore the vast chemical space, by 63 significantly narrowing down the number of drug candidates that need to be synthesized and 64 tested in experiments. Calculating the interaction between the receptor and ligand can be done in 65 various ways: molecular docking, molecular dynamics, fingerprint, pharmacophore matching, 66 and so forth. As its name implies, the method requires the target protein structure, and the 67 performance depends on the protein conformation. For targets whose experimental structures are 68 not available, researchers should predict the three-dimensional structures using methods such as 69 homology modeling. 70 One of the major obstacles in the SBVS method, especially molecular docking, is caused by its 71 static picture of a receptor structure. Proteins are flexible molecules, so they can change their 72 shape depending on the binding partners. The structural change for receptor protein might lead to 73 a failure of molecular docking [1-3]. One of the techniques for treating receptor flexibility is 74 ensemble screening. The method uses pre-generated diverse receptor structures gathered from 75 experimental structure databases or simulation trajectories such as molecular dynamics or normal 76 mode analysis. Ligands in the screening library are docked to individual structures and ranked by 77 their representative scores. There are various methods to get the representative score: arithmetic 78 mean, harmonic mean, etc. Since the ensemble method uses multiple target structures, it is 79 important to reflect the structural diversity when selecting the receptor ensemble. However, the .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 5 80 crystal structures of the target would be biased to thermodynamically stable states or the major 81 type of inhibitor-bound form. This might make the screening for obtaining diverse scaffolds or 82 even lead to a failure in the SBVS even structural ensembles are used. 83 Kinases are one of the typical examples having structural diversity. They play a key role 84 in biological processes in phosphorylation transferring a phosphate group from ATP to proteins. 85 Since the process modulates key cellular operations like cell cycle regulation, metabolism, and 86 apoptosis, kinases have become attractive targets in drug discovery. According to Santos et al. 87 [4], kinases belong to one of the four main target families, which are G protein-coupled receptors 88 (GPCRs), ion channels, nuclear receptors, and kinases. In the ChEMBL database [5], the kinases 89 are targeted by 14.2% of compounds. 90 The kinase domain, a structural domain with catalytic function, has a highly conserved structure. 91 It is composed of an N-lobe and a C-lobe, linked by a hinge region [6]. The N-lobe is structured 92 with five β-strands and a single alpha helix, known as the C-helix, while the C-lobe is comprised 93 of multiple alpha helices. The catalytic activity of kinase domains is regulated by two loops, 94 namely the activation loop and the catalytic loop. The His-Arg-Asp (HRD) motif of the catalytic 95 loop directly interacts with the hydroxyl group of the target protein residue (serine, threonine, or 96 tyrosine) set for phosphorylation. The Asp-Phe-Gly (DFG) motif found in the N-terminal of the 97 activation loop plays a crucial role in anchoring ATP to the active site. The conformational states 98 of the active site of kinase domains are classified as DFGin, DFGinter, or DFGout, based on the 99 orientation of the aspartic acid of the DFG motif relative to the site. Fig 1 illustrates two distinct 100 structural states, DFGin (Fig 1A) and DFGout (Fig 1B), of BRAF. The DFGin state locates Phe 101 of the DFG motif into the ATP binding pocket, allowing kinase to hold ATP, thus it is called the .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 6 102 active state. On the other hand, DFGout conformation directs the Phe out of the ATP binding 103 pocket, so it is classified as an inactive state. 104 105 Fig 1. The structure of kinase domain and key structural elements. The active state (DFGin- 106 BLAminus) of BRAF (PDB ID: 6UAN, A) and the inactive state (DFGout-BBAminus) of BRAF 107 (PDB ID: 8C7X, B). The active site is focused in the red box on the right of each panel. The 108 activation loop and DFG motif changes conformation depending on the state. 109 110 The kinase inhibitors can be categorized depending on the binding site conformation of kinases 111 and the place the compounds bind. Type I inhibitors are the ones that bind to the ATP-binding 112 site and thus compete with ATP. Most type I inhibitors bind to the DFGin state. We observed 113 that the majority of experimentally determined human kinase structures form DFGin state (87%, 114 as of May 2023). On the other hand, type II inhibitors often associate with the kinase to the 115 DFGout state. The type II compounds tend to occupy ATP-binding site partially and a 116 hydrophobic pocket close to the ATP-binding site, which is opened when the activation loop 117 forms the DFGout conformation. In general, type II inhibitors have more selectivity than type I 118 inhibitors. According to Hari et al. [7], the type II inhibitors' selectivity is influenced by the 119 inherent variations in kinase's capacity to adopt the DFGout conformation. Lastly, type III 120 inhibitors, sometimes referred to as allosteric inhibitors, bind to a kinase site that is not directly 121 connected to the ATP-binding site. Therefore, it is crucial to take into account as many different 122 structural states of kinases as possible to find diverse hit molecules for SBVS. 123 Recent advances in protein structure prediction using deep learning techniques, such as 124 AlphaFold2 (AF2) [8] and RoseTTAFold [9], allow for accurate modeling of protein structures. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 7 125 However, these methods rely on pre-trained models from the PDB database, so the generated 126 models might be affected by the conformational state distribution in the PDB if the protein can 127 form diverse structures. Therefore, the methods could produce kinase structures similar to the 128 DFGin state, and thus SBVS targeting kinases using the predicted models could potentially yield 129 outcomes favoring type I inhibitors even using the modeled structures. 130 To solve the issue, Heo and Feig [10] suggest a method called multi-state modeling (MSM) to 131 predict GPCRs with high accuracy of the desired structural state. The authors found that 132 predicted models by AF2 tend to have either active or inactive states depending on GPCR classes 133 due to the small number of experimentally determined structures of GPCRs. Instead of providing 134 multiple sequence alignment (MSA) as an input of AF2, the authors align a query sequence to a 135 template sequence with the structural state of interest. The MSM technique showed an improved 136 performance of modeling for GPCRs. Cognate docking using MSM models also showed 137 enhanced accuracy with smaller root-mean-square-distance (RMSD) from the crystal structure. 138 Inspired by the work, we established an MSM protocol for modeling the kinase structures 139 to overcome the structural bias of kinases for SBVS by giving state-specific templates to AF2. 140 All human kinase experimental structures were classified by the active site conformation using 141 KinCoRe [11] to construct a state-specific template database. Our protocol was able to predict 142 kinase conformations with the desired structural state with a high accuracy. Then we 143 benchmarked cognate docking to observe the MSM models produced accurate binding modes 144 with lower RMSD than standard AF2 models. Finally, we performed ensemble SBVS with 145 generated models by MSM. The ensemble method showed higher performance compared to 146 screening with modeled structures using AF2, especially when the active molecules are diverse. 147 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 8 148 Results and Discussion 149 Conformation Distributions of Experimental Kinase Structures and 150 AlphaFold2 Predicted Models 151 Proteins are flexible molecules, so they may change their structures when bind to their 152 partners. This conformational change could act as an obstacle to SBVS. The active site of 153 kinases also forms diverse conformations depending on the type of the binding compound. For 154 example, the DFGin state binds to the type I inhibitor, while the DFGout state binds to the type II 155 inhibitor. Dunbrack and his colleagues introduced a rule, called KinCoRe [11], to classify the 156 kinase structure. The rule classifies structural states of kinase structures into 12 categories based 157 on the spatial state of the activation loop and the dihedral angle of DFG motif. Details of the 158 criteria can be found in the Materials and Methods section and Modi et al [11]. 159 Fig 2 shows a distribution of kinase conformational states annotated by KinCoRe scheme [11] in 160 the PDB database (blue bar). Based on the KinCoRe notation, more than half (53.6%) of 161 experimentally determined kinase structures have DFGin-BLAminus conformation. Details of 162 the structure distribution are shown in Supplementary Material (S1 Table). Other than the major 163 state, the other conformational states occupy less than 10% of PDB structures. 164 165 Fig 2. Distribution of crystal structures and standard AF2 structures for each kinase 166 conformation. The X-axis is conformational states annotated by KinCoRe and the y-axis is the 167 percentage of each state. All human crystal kinase structures deposited in RCSB-PDB, DUD-E 168 target protein crystal structures, and predicted AF2 models with default parameters (standard 169 AF2) are colored as blue, orange, and green, respectively. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 9 170 171 The preference for the DFGin state in experimentally determined structures might be influenced 172 by the thermodynamic stabilities of the conformational states. Meng et al [12], studied the 173 transition between the structural states of c-Abl and c-Src kinases using umbrella sampling and 174 the potential of mean force. With a difference of 1.4 kcal/mol, the DFGin state of c-Abl has a 175 lower Gibbs free energy than the DFGout state. The calculation shows the active conformation 176 (DFGin) is the dominant population while the DFGout state occupied only 9% in their 177 simulation. Similarly, for c-Src, the DFGin state is more favored than the DFGout state. The free 178 energy difference between the states is calculated as 5.4 kcal/mol. The authors also studied a 179 thermodynamic barrier of the transition between the states [13]. It is estimated as the order of 2-4 180 kcal/mol making the transition from highly populated the DFGin state to the DFGout state hard. 181 Levy and his colleagues [14] collected 2,896 kinase structures and multiple sequence alignment 182 and applied Potts model to predict structural propensities from sequences. With the statistical 183 potential, most of the kinases are predicted to have preferences for the DFGin state. The penalty 184 for forming the DFGout state reaches 2-3 kcal/mol for the extreme case. 185 Not only for the thermodynamic preference, but the type of kinase inhibitor distribution might 186 also affect the skewness of kinase structures. We counted the number of compounds in kinase 187 inhibitor types from PKIDB (assessed August 2023) [15]. PKIDB is a curated database of kinase 188 inhibitors in clinical trials. Out of 369 compounds in the database, only 84 molecules have their 189 inhibition type annotation, because annotating the ligand type needs the complex structure. Type 190 I inhibitor (DFGin bound) is the dominant form, occupying 66% (55 compounds) of the 191 annotated inhibitors. In contrast, 17 molecules are labeled as type II (DFGout bound). Thus, 192 DFGin conformation (active state) might have more chance to be crystallized than the DFGout .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 10 193 state. The biased conformational states of kinases would make the discovery of chemically 194 diverse kinase inhibitors harder. 195 We also examined the conformational state distribution of 25 DUD-E targets [16], that are 196 benchmarked throughout this study, in the PDB database (orange bar). DUD-E is a widely used 197 benchmark dataset for evaluating virtual screening methods, composed of pharmaceutically 198 important targets such as GPCRs, kinases, and nuclear receptors. Among the set, we extracted 25 199 kinase targets for benchmarking. The experimental structures are less biased than all PDB 200 structures, about 30% of DUD-E proteins have DFGin-BLAminus conformation. The other 201 states tend to have a higher proportion than the PDB structure distribution of the states. The 202 difference in the state distribution between all kinases and DUD-E targets potentially means that 203 the highly biased nature of kinase structures might not be suitable for finding diverse types of 204 kinase inhibitor discovery. 205 From an MSA of a given sequence, AF2 [8] extracts coevolution information by the MSA 206 Transformer and predicts the three-dimensional structure based on the information and deep- 207 learning models trained on existing protein structures. We modeled 25 kinase catalytic domains 208 provided in the DUD-E kinase subset with the default parameters of AF2 resulting in 125 209 structures (five models per target), then assigned the conformational states of the models using 210 KinCoRe. Throughout this paper, AF2 with default parameters is called standard AF2. The 211 average plDDT and MolProbity [17] score of the predicted structures is 89.38 and 1.04, 212 respectively. This implies that they are properly modeled. Out of 125 predicted models, 91 213 structures (72%) are annotated to have DFGin-BLAminus conformation (green bar). The 214 distribution of predicted models is more skewed than all human kinase structures and the DUD-E 215 set proteins. However, the other states have a lower proportion than the experimental structures. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 11 216 The standard AF2 did not produce DFGin-ABAminus, DFGin-BLBtrans, DFGinter-BABtrans, 217 and DFGout-Unassigned conformations, which occupy a small portion of the human kinase PDB 218 structures (7.3%, 3.0%, 0.3%, and 3.2%, respectively, S1 Table). Since the experimentally 219 determined kinase structures are biased to the DFGin-BLAminus, predicted kinase structures by 220 AF2 might have a chance to be biased to the conformation, which can be observed in previous 221 research [10]. In addition to the biased trained models, the template selection process in AF2 222 modeling does not take the structural state of the kinase into account. 223 224 Predicting State-specific Kinase Structures using Multi-state 225 Modeling Protocol 226 AF2 with MSM protocol provides conformational state-specific structures as templates 227 for AF2 models. All human kinase structures from the KLIFS database [18] and further 228 classified by their structural state following KinCoRe rules [11]. For each state, the five highest 229 sequence similar structures were selected as templates. Then AF2 was executed with the 230 template information to predict five models for each template. The models with undesired 231 conformational state were removed, and then the model with the highest plDDT was selected for 232 our benchmark. The overall workflow is illustrated in Fig 3. Details of the kinase MSM protocol 233 are elaborated in the Materials and Methods section. 234 235 Fig 3. Workflow of multi-state modeling of kinases 236 237 We prepared two template sets, a trivial template (TT) set containing 100% sequence identical 238 structural template and a nontrivial template (NT) set not having identical protein. With our .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 12 239 MSM protocol, AF2 was able to model structures in a specific state, producing 8.5 and 8.2 240 models with TT set and NT set on average, respectively. The number of predicted models for 241 each target ranges from three (WEE1) to 11 (ABL1 and LCK, S2 Table). The average plDDT 242 scores of structures are 89.63 and 88.08 for models from the TT set and NT set, respectively, 243 showing comparable values to the standard AF2 predictions (89.38). As expected, TT set models 244 have higher accuracy on average than NT set, but their difference is marginal. The distribution of 245 the MolProbity score is given in S1 Fig. We also examined the quality of models for each 246 structural state. The average plDDT values range from 86.62 (DFGin-BLAplus) to 94.48 247 (DFGin-BLAminus), meaning that the quality of the models does not depend much on the 248 structural state of kinases. The highest plDDT score of DFGin-BLAminus might be caused by 249 the highest frequency of the state in the PDB database (S1 Table). The Pearson’s correlation 250 coefficient between the average plDDT and percentage in the PDB structure of each structural 251 state is 0.739. Even though the MSM protocol provides a structural template for AF2, the model 252 quality might be influenced by the pre-trained models of AF2, so the average plDDT of each 253 model follows the distribution of the crystal structure. 254 The TT set model and the NT set model have average MolProbity scores of 1.21 and 1.24, 255 respectively. These figures represent a modest decline from the baseline AF2 models (1.04). We 256 also examined each structural state's average MolProbity score. 1.06 (DFGin-ABAminus) is the 257 lowest value, and 1.33 (DFGinter-Unassigned) is the highest. For DFGin-BLAminus, the highest 258 populated states for both PDB and standard AF2, MSM models show the average MolProbity 259 score as 1.08. The average MolProbity score and the distribution in the PDB structure have a - 260 0.59 Pearson's correlation coefficient, indicating a weak correlation, and the high score of 261 MolProbity might be caused by the states with the small PDB populated states. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 13 262 To investigate the accuracy of the models, we measured the TM-Score [19] of predicted models 263 to crystal structures given in the DUD-E set, called ‘reference structure’ throughout this paper. 264 TM-Score assesses a structural similarity between two given proteins, ranging from zero (not 265 similar) to one (identical). For MSM models, we used the predicted structures in the same 266 structural state as the reference. The average TM-Scores are 0.92 and 0.90 for the models 267 predicted using TT and NT sets, respectively. Compared with standard AF2 models (average 268 TM-Score: 0.87), the MSM technique provided more similar models than the standard AF2 269 protocol, which is expected since the MSM provides structural templates. Models with TT sets 270 generally have more accurate structures than those with NT sets as also expected. 271 Consequently, by utilizing structures that represent a variety of states for the target kinase, the 272 MSM is able to generate diverse structures as desired with high accuracy. Thus, it could provide 273 a proper structure set for kinase ensemble SBVS. 274 275 Cognate Docking Accuracy of a Compound to the Multi-state 276 Modelled Structures 277 To examine whether modeled structures are suitable for molecular docking and thus 278 structure-based virtual screening or not, we first conducted a cognate docking experiment on the 279 predicted structures, both standard AF2 and MSM models. Ligands from complex crystal 280 structures of the DUD-E kinase subset were used for this benchmark. For the standard AF2 281 model, the highest plDDT model for each protein, which is also used for performing virtual 282 screening in the next section, was selected for evaluation. For evaluating the MSM protocol, we 283 used the modeled structure with the same KinCoRe annotation as the reference PDB structure 284 provided by DUD-E. Out of 25 kinases in DUD-E, IGF1R was removed, since the crystal .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 14 285 structure was not assigned any of the structural states by KinCoRe, Therefore, the number of 286 target proteins becomes 24. AutoDock-GPU [20] was employed to predict the 50 binding poses. 287 The RMSDs of all predicted poses to the crystal binding pose were calculated. To analyze, we 288 took three values: the RMSD of the best AutoDock score pose, that of the closest pose to the 289 crystal binding mode, and the average RMSD of all 50 poses. 290 Table 1 summarizes the docking accuracy evaluation results. Taking the best AutoDock score 291 models, the standard AF2, MSM structures modeled with TT, and NT have average RMSD of 292 2.74 Å, 2.15 Å, and 3.49 Å, respectively. In addition, the success cases with an RMSD cutoff of 293 2.0 Å, a standard criterion for judging docking success [21-23], are 11 (standard AF2), 16 (with 294 TT), and 9 (with NT) out of 24 receptors. Individual RMSD values are given in S3 Table. 295 296 Table 1. Docking accuracy benchmark result. The values are the average values of 24 proteins 297 and the numbers in the parentheses are the number of success cases with RMSD < 2 Å. Structure Standard AF2 MSM with TT MSM with NT Best AutoDock Score 2.74 (11) 2.15 (16) 3.49 (9) Lowest RMSD 1.50 (20) 1.40 (19) 2.14 (14) Average RMSD 2.93 2.49 3.34 298 299 Comparing the MSM models and AF2 models, the multi-state models with TT sets have the most 300 accurate docking poses. As observed in the previous section, template information influenced the 301 quality of the docking poses, i.e., the predicted docking poses of the models with TT sets have 302 smaller RMSD and more successful cases than those with NT sets. One successful example is .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 15 303 AKT2 (Fig 4). Although the TM-scores of MSM with TT set (0.98) and standard AF2 (0.96) are 304 similar, the best AutoDock score models docked to the MSM using TT sets has more accurate 305 binding pose than standard AF2, RMSD of 0.86 Å (magenta, Fig 4B) and 3.24 Å (cyan, Fig 4C), 306 respectively. 307 308 Fig 4. Predicted cognate docking structures of AKT2. A. Superimposed structures of crystal 309 structure (PDB ID: 3D0E, grey), MSM model (DFGin-BLAminus, magenta; TM-Score: 0.98), 310 and standard AF2 (cyan; TM-Score: 0.96). B. Predicted docking poses of the cognate ligand on 311 the AF2 with MSM. The RMSD of the predicted pose is 0.86 Å. C. Predicted docking poses of 312 the cognate ligand on the standard AF2. The RMSD of the predicted pose is 3.24 Å. 313 314 Although the docking is not successful (RMSD > 2 Å), MSM structure was able to retrieve the 315 interaction between protein and ligand found in the crystal structure for some cases. Out of eight 316 cases, the interacting residues in the reference structures were successfully captured more than 317 50% in five proteins analyzed by PLIP [24] (S4 Table). For example, the predicted binding pose 318 of the best AutoDock score conformation of PRKCB cognate docking was 2.79 Å when the 319 MSM with TT structure was used as the receptor structure. However, out of eleven interacting 320 residues identified in the reference structure, ten residues were retrieved in the MSM-ligand 321 complex model. By analyzing the interaction pattern, the predicted docking pose to the multi- 322 state modeled structure has hydrophobic interaction with L348, F353, V356, A369, K371, A483, 323 and D484 and hydrogen bonding with T404, E421, and V423. These interactions are also 324 observed in the crystal structure (S2 Fig). The cognate docking benchmark result would imply .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 16 325 that MSM models could be used to predict binding poses of kinase-ligand complexes, thus they 326 are suitable for virtual ensemble screening. 327 Similar trends are observed for the lowest RMSD conformations and average RMSD of 50 328 conformations. Among the AF2 predicted structures, MSM with TT shows the most accurate 329 models (average of the lowest RMSD pose: 1.40 Å, success case: 19) and the standard AF2 330 performs slightly lower (average of the lowest RMSD pose: 1.50 Å, success case: 20). The 331 docking accuracy of MSM with NT sets is the worst among the receptor structure sets (average 332 of the lowest RMSD pose: 2.14 Å, success case: 14). In most cases, the best AutoDock score 333 conformations do not match to the lowest RMSD conformations, which means that the 334 AutoDock score could not be able to find the optimal docking poses. The average RMSD of 50 335 conformations follows the same order as the other two metrics: MSM with TT is the smallest, 336 and MSM with NT is the highest. 337 338 Virtual Screening Performance with Multi-state Models 339 To investigate the advantage of using the MSM technique for SBVS, the compound 340 library for each kinase protein from DUD-E docked to structures with diverse states generated by 341 using our method and compared to the structures modeled with standard AF2. AutoDock-GPU 342 was used for the benchmark. For ensemble docking using generated models by MSM, since a 343 molecule was docked to a couple of receptor structures and each compound-structure pair had a 344 docking score, we needed to decide representative scores of the compounds to rank the 345 molecules. We employed two types of representative scores: the AutoDock best (ADB) and the 346 Boltzmann-weighted (BW) scores. The ADB score picks the lowest AutoDock score among the .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 17 347 docked results, while the BW score is a weighted average of the scores across all structures. The 348 details of the BW score are given in the Materials and Methods section. 349 Table 2 shows the performance of SBVS using the various receptor structure sets. To evaluate 350 the performance, five metrics are used: enrichment factors at 1%, 5%, and 10% (EFX%), an area 351 under ROC curve (AUC), and Boltzmann-enhanced discrimination of receiver operating 352 characteristic (BEDROC). EF at X% indicates a capability of finding active molecules within the 353 top X%, and AUC represents the discrimination power of a screening method between active and 354 decoy molecules. BEDROC puts exponential weights to the early rank of molecules, thus it is 355 able to solve the ‘early recognition problem’ caused in AUC [25]. Details of the metrics are 356 illustrated in the Materials and Methods section. As observed in the cognate docking benchmark, 357 MSM structures are better or equal to standard AF2 results, regardless of the scoring method or 358 the template set for modeling. Details of individual results are provided in S5 Table. Comparing 359 EF1% values target-by-target, the MSM performed better than or equal to standard AF2 in 18 360 proteins out of 25 targets (72%) using the ADB score screened to the structures modeled with the 361 NT set. Even with the lowest EF1% combination, TT models with BW scoring, the MSM 362 performed better than or equal to standard AF2 in more than half of proteins (13 proteins). 363 Interestingly, although the TT set models with the same KinCoRe classification as the reference 364 have more accurate structure and docking poses than the NT set models, the virtual screening 365 performance is slightly worse in both scoring schemes. 366 367 Table 2. Performance of structure-based virtual screening on various receptor models. MSM with TT MSM with NT ADB BW ADB BW Standard AF2 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 18 DUD-E Kinase Subset (25 proteins) EF1% 7.2 7.2 8.2 8.0 6.6 EF5% 3.5 3.7 3,7 3.7 3.5 EF10% 2.7 2.7 2.6 2.7 2.6 AUC 0.660 0.639 0.661 0.664 0.659 BEDROC 0.193 0.196 0.183 0.189 0.183 Average Dissimilarity ≥ 0.7 (14 proteins) EF1% 5.3 5.6 6.4 6.9 4.6 EF5% 3.0 3.1 3.2 3.2 2.8 EF10% 2.4 2.4 2.2 2.3 2.1 AUC 0.646 0.652 0.648 0.650 0.641 BEDROC 0.166 0.171 0.154 0.161 0.146 Average Dissimilarity < 0.7 (11 proteins) EF1% 9.6 9.3 10.6 9.3 9.1 EF5% 4.2 4.4 4.4 4.3 4.3 EF10% 3.0 3.1 3.1 3.1 3.3 AUC 0.678 0.681 0.678 0.681 0.681 BEDROC 0.228 0.228 0.220 0.226 0.229 368 369 To elucidate the variance in performance between MSMs using TT and NT models, we 370 identified kinases that exhibited notable differences in EF1% between the two sets. Of the kinases 371 studied, both ABL1 and KDR demonstrated superior performance using the MSM with the NT 372 set compared to the TT set, across both ensemble scoring methods (S5 Table). We measured .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 19 373 TM-Scores of the individual models with the same structural state, those of templates used to 374 predict the structures, and EF1% of the models (S6 and S7 Tables). Among the kinase states 375 analyzed, the DFGout-BBAminus for both proteins produced a significant difference in both 376 templates and models. When the kinase forms an active state, the activation loop forms an 377 extended conformation to facilitate the catalytic function of the protein, thus the DFGin 378 conformation has a conserved structure. On the other hand, for the DFGout conformation, the 379 activation loop collapsed on the protein surface with high flexibility [11]. Therefore, although 380 the protein structures have the same KinCoRe notation with DFGout, the activation loop can 381 have different structures. The TM-Scores of the templates to model the DFGout-BBAminus state 382 are 0.88 and 0.96 for ABL1 and KDR, respectively. The structural difference of templates 383 influenced the predicted models, resulting in TM-Scores 0.89 for both proteins. S3 Fig shows a 384 structural difference of ABL1 predicted models with DFGout-BBAminus conformation. The 385 difference of activation loop location also might affect to the virtual screening performance, 386 leading TT set (3.28 and 8.78 for ABL1 and KDR, respectively) has lower EF1% value than NT 387 set (15.83 and 15.85 for ABL1 and KDR, respectively). Significantly, this difference of 388 performance in the DFGout-BBAminus state impacted the overall ensemble docking results. 389 One benefit of using the MSM models is that the predicted models are diverse, so it is potentially 390 useful for discovering various scaffolds of hit chemicals. To examine whether the hypothesis is 391 true or not, we divided the DUD-E kinase subset into two based on the diversity of active 392 compounds in the screening library. The pairwise Tanimoto coefficients (Tc) between the active 393 compounds were calculated by using RDKitFP fingerprint in RDKit [26]. Then the pairwise 394 distances between the compounds are calculated as (1 – Tc). The diversity of active compounds 395 is defined as the average distances of the active compounds (S5 Table). With a threshold of 0.7, .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 20 396 the kinases were classified into two groups, 14 proteins with higher than or equal to the cutoff 397 and the remaining 11 targets. 398 For the kinases with more diverse active compounds, the ensemble screening with MSM models 399 showed higher performance than standard AF2 models in all metrics. For EF1%, the MSM 400 performed much better than standard AF2 (Table 2). Out of 14 proteins, the MSM with NT set 401 and BW scoring has higher or equal EF1% values than standard AF2 models in 11 targets. On the 402 other hand, when the active compounds become less diverse, MSM still performed better than 403 the standard AF2, but the gap between them declined for EF1%. This implies that our approach 404 would be powerful for discovering diverse molecular scaffolds. 405 One of the diverse active hit examples is ABL1. The average dissimilarity of the active 406 compounds is 0.73. Regardless of the template set and the scoring scheme, EF1% of MSM 407 ensemble screening showed higher performance (TT models with ADB: 9.3, TT models with 408 BW: 10.4, NT models with ADB: 15.3, and NT models with BW: 17.5) than standard AF2 409 model (6.0). We also examined the diversity of active compounds ranked within the top 1% 410 ranked molecules. Our ensemble protocol tends to find diverse hits: 0.51, 0.53, 0.63, and 0.62 for 411 TT models with ADB, TT models with BW, NT models with ADB, and NT models with BW, 412 respectively. In contrast, the diversity of active compounds using the standard AF2 is 0.27. 413 CSF1R is another example showing MSM ensemble screening is able to find diverse hit 414 compounds. The average dissimilarity of hit molecules is 0.73. The gap of EF1% between MSM 415 models (5.4) and the standard AF2 model (3.6) is smaller than the case of ABL1. The diversities 416 of discovered hit compounds within the top 1% are 0.70, 0.65, 0.73, and 0.73 in the order of TT 417 models with ADB, TT models with BW, NT models with ADB, and NT models with BW, 418 similar to the average dissimilarity of all active compounds. However, the hit compounds within .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 21 419 the top 1 % identified using the standard AF2 model is 0.23. Fig 5 shows the docking pose of 420 one of the active compounds, ChEMBL245377. The MSM ranked the active compound within 421 the top 1% (10th using TT models and BW scoring, out of 12316 compounds) which is not 422 ranked by the standard AF2 structure (167th). The compound is not successfully docked using the 423 the standard AF2 model. 424 425 Fig 5. Predicted binding poses of ChEMBL245377 and CSF1R structures. A. The 426 superimposed complex structures of ChEMBL245377 and CSF1R. The crystal structure of 427 CSF1R (PDB ID:3KRJ), standard AF2 predicted structure, and MSM model (DFGout- 428 BBAminus, with TT sets) are represented as ribbon diagrams colored as gold, pink, and orange, 429 respectively, while the predicted docking poses of the compound are represented as sticks. B,C. 430 Focused binding sites and docking poses of ChEMBL245377. The compound is ranked 10th by 431 ensemble docking (B) and 167th in docking for the standard AF2 predicted structure (C), 432 433 Although the ensemble screening with MSM structures performed generally better than standard 434 AF2, there is an issue with selecting the representative score of a compound. For instance, for 435 FGFR1, the MSM model with TT set and BW score showed EF1% as 2.1. When we observed 436 EF1% values of individual models, the DFGin-BLBplus state outperformed than any other 437 structures including standard AF2 (Table 3). Thus, a proper method for selecting or calculating 438 representative scores for a compound should be designed to achieve high performance for MSM 439 ensemble screening. 440 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 22 441 Table 3. EF1% for FGFR1 structures of specific states. The MSM structures are built from TT 442 set. Structural State EF1% DFGin-BLAminus 2.9 DFGin-BLAplus 6.4 DFGin-BLBminus 0.7 DFGin-BLBplus 8.6 DFGin-BLBtrans 2.1 DFGout-BBAminus 3.6 DFGout-Unassigned 1.4 MSM with BW score 2.1 Standard AF2 1.4 443 444 Conclusion 445 The receptor conformation affects SBVS performance. Like other proteins, kinase adjusts 446 the conformation of its binding site in response to the binding ligand. Therefore, it is crucial to 447 have adequate kinase structure to obtain inhibitors with the required mode of action or diversity. 448 For human kinase structures that were identified through experiments, however, there is a clear 449 bias toward the active state. The prediction of the AF2 structure could be influenced by the bias 450 in the PDB database. We noticed that there is a bias toward the active state in the predicted 451 structures with standard AF2 protocol. The results of the SBVS using the predicted structure 452 would be compromised by this bias in receptor structure. To overcome the bias, we applied the 453 MSM technique by providing a structural template to AF2 to generate structures with diverse .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 23 454 states and using the models for ensemble docking. Compared with standard AF2 models, MSM 455 protocol produced more accurate or comparable models although it did not give MSA as an 456 input. Also, in cognate docking study, MSM models provided ligand docking poses close to the 457 crystal structure. With the diverse predicted kinase structure, we performed ensemble screening. 458 The ensemble method showed enhanced or comparable SBVS performance to the standard AF2 459 modeled structure result. We also observed that our method would be more suitable when the 460 ligands are diverse, leading to the identification of a diverse range of kinase inhibitors. Even for 461 the targets that ensemble docking method could not find active compounds, we found that some 462 of the models outperformed the standard AF2. Thus, the selection of a representative structure 463 should be improved and remains the next work for this project. 464 Ensemble screening with MSM models would open the possibility of uncovering novel 465 kinase inhibitors with diverse chemical scaffolds. It has advantages in addressing current 466 challenges in kinase inhibitor development for finding chemically diverse compounds. The 467 chemical diversity of kinase inhibitors could aid in overcoming the problem of drug resistance 468 generally caused by the mutation, a significant obstacle in kinase-targeted cancer therapies. 469 Additionally, it could increase chance to find hit compounds, not similar to the existing patents. 470 By exploring a diverse array of kinase inhibitors with accurately predicted structures, we would 471 be able to find inhibitors with different modes of action. This could lead to the development of 472 novel therapeutic strategies that are more robust in the face of drug resistance. Hence, our 473 approach could potentially be applied to more effective and precise kinase-targeted therapies. In 474 addition, the MSM method could be applied to important therapeutic targets such as GPCRs and 475 nuclear receptors. 476 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 24 477 Materials and Methods 478 Benchmark Dataset 479 To investigate whether our approach can generate a structural ensemble properly and 480 improve SBVS performance or not, we selected a kinase subset of DUD-E [16]. For each target, 481 a reference PDB structure for screening and a compound library composed of active and decoy 482 molecules are provided in the DUD-E set. The active compounds of DUD-E set are composed of 483 the molecules with affinity 1 μM or better were extracted from ChEMBL09 [5]. The reference 484 structures were selected by considering the resolution and enrichment of finding active 485 molecules by DOCK3.5. The decoys of the DUD-E set are constructed by gathering compounds 486 with similar characteristic to the active molecules such as logP and number of rotatable bonds 487 from ZINC database [27]. The kinase subset, which is used in this work, consists of 26 kinases 488 with include 205.6 actives and 12,830 decoys on average. Among the 26 targets in DUD-E 489 kinase subset, SRC kinase was removed from DUD-E benchmark set because the given reference 490 structure is not originated from human. 491 492 Kinase Structural State Annotation 493 In the active site of protein kinases, the activation loop, 20-30 residues long, is the most 494 important secondary structural element [28] for determining the structural state. The loop starts 495 from the conserved three-residue-long sequence, the DFG motif. In this work, the standalone 496 version of KinCoRe [11] (https://github.com/vivekmodi/Kincore-standalone, Accessed 497 4/14/2022) was employed to annotate the conformational state of all experimental and modeled 498 kinase structures. The program categorizes the conformational state into 12 classes by the .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 25 499 location of the activation loop and dihedral angles of the DFG motif. The spatial state of the 500 activation loop is defined by two distances: 1) a distance between Phe-ring of the DFG motif and 501 C atom of the fourth residue from the conserved Glu in the C-helix of N-lobe and 2) a distance 502 from the Phe-ring to the conserved C atom of conserved Lys in 3 strand of N-lobe (Fig 1). 503 Based on the distances, the activation loop location is classified into three classes: DFGin, which 504 is the Phe-ring located under C-helix, DFGout, the Phe-ring is moved into ATP binding pocket, 505 and DFGinter, an intermediate state between DFGin and DFGout. The program further classifies 506 the activation loop structural state by calculating dihedral angles: ,  backbone dihedral angles 507 of X-DFG (a residue before the DFG motif), Asp, and Phe of DFG motif, and 1 angle of DFG- 508 Phe. As a result, DFGin, the dominant class, has seven subclasses (BLAminus, BLAplus, 509 ABAminus, BLBminus, BLBplus, BLBtrans, and Unassigned), while DFGinter (BABtrans and 510 Unassigned) and DFGout (BBAminus and Unassigned) have only two subclasses. The three 511 letters after the activation loop states follows the region of Ramachandran map occupied by X, 512 D, and F residues: A, B, L for alpha, beta, and left-handed, respectively. The 1 angle of Phe is 513 indicated as plus (+60 degree), minus (-60 degree), and trans (180 degree). The last class is 514 Unassigned-Unassigned, the activation loop and DFG conformations cannot be determined. 515 516 Construction of Template Database for Each Structural State 517 To construct the structural template database for MSM, KLIFS [18], a database of 518 experimentally determined kinase structures was used. The database contains catalytic domain 519 structures of kinases, extracted from PDB and their inhibitors and provides the interaction 520 information between the protein and the compound. As of May 2023, the database is composed 521 of 6,344 structures (13,382 monomers). Among the kinase structures in KLIFS (Accessed .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 26 522 1/18/2023), we filtered out the non-human proteins and proteins that produced errors during 523 KinCoRe annotation, resulting in 11,106 monomer structures. To construct the template structure 524 database for each state, the crystal structures with the same annotation by KinCoRe were 525 gathered. 526 527 Standard AlphaFold2 Modeling 528 The kinase structure was modeled with the standard protocol of AF2 (v.2.3.1) to compare 529 with MSM models. Only the kinase domain was modeled from the full sequence of a protein. 530 The MSA for the kinase sequence was generated from BFD (v.3.2019), MGnify (v.5.2022), and 531 UniRef90 (v.2.2022) via HHblits (from HH-suite v3.3.0) and Jackhammer (from HMMER 532 v3.3.2). Four highest sequence identity proteins with 3D atomic coordinates were selected to 533 provide template structures. The number of recycles was set to three and the model relaxation 534 step was integrated into the procedure. As AF2 has five different trained models and runs all of 535 them independently in a single run, five structures were generated from a single run. The models 536 with the highest plDDT score out of the five predicted structures for the comparison, since 537 plDDT is a confidence measure of AF2 predicted models. 538 539 Multistate Modeling of Kinase using Structural Template 540 The workflow of MSM is given in Fig 3. From a given target kinase sequence to be 541 modeled, the templates were searched by MMseqs2 (release 11) easy-search (e-value cutoff: 1e- 542 3) [29] against all sequences in each structural state. For each structural state, the top five 543 templates, which were determined by e-value, were used for the modeling. To mimic the real 544 drug discovery process and check the influence of the template for virtual screening, we .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 27 545 generated two template sets for modeling: one set has the query protein which means 100% 546 sequence identity with the query sequence (TT set), and the other template sets without the query 547 protein (NT set). Modeling with each template was conducted independently. 548 Since AF2 produced five models per single run, 25 structures for each specific state were 549 generated in total. Among them, we selected one structure for benchmarking with two criteria: 550 conformational state and quality of the predicted model. First, we filtered out the models with 551 different KinCoRe annotations from the template structure classification. For example, to predict 552 a model the DFGin-ABAminus conformation of ABL1, five templates structure of TYR family 553 (LYN: 5XY1, EPHA2: 7KJB, 5NK3, 4TRL, and IGF1R: 3F5P) were selected. However, the 554 models were annotated as DFGin-BLAminus rather than DFGin-ABAminus conformation (S4 555 Fig). Thus, all predicted models were discarded. After filtering by the KinCoRe annotation, the 556 models with plDDT less than 70 were also removed. Among the remained structures, the model 557 with the highest plDDT score was finally selected. Other details for modeling are the same as 558 standard AF2 modeling. 559 560 Assessment of Model Quality 561 The TM-score [19] evaluates the structural similarity of protein structures. It is scaled 562 according to the size of the protein and exhibits a better sensitivity to the overall structural 563 alignment compared to the RMSD. The models with the same KinCoRe notation as the crystal 564 structures, extracted from the DUD-E set were chosen for comparison with the crystal structure 565 employing the TM-score. The TM-score's values range from 0 to 1, where 1 signifies a perfect 566 alignment. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 28 567 MolProbity [17] is a comprehensive validation tool for the structural integrity of proteins and 568 nucleic acids. MolProbity validates protein structures through hydrogen replacement, 569 comprehensive all-atom contact analysis, and evaluation of torsional angles. The MolProbity 570 score integrates various factors into a single metric to indicate the model's reliability, where a 571 lower score denotes a better model. The models from the standard AF2 and MSM were validated 572 using MolProbity implemented in Phenix software [30]. 573 574 Docking and Virtual Screening using AutoDock-GPU 575 AutoDock-GPU 1.5.3 [20], which is open-source and GPU-accelerated, was employed to 576 benchmark the virtual screening performance of kinase structures. We used AutoDockTools [31] 577 to convert the receptor PDB files to PDBQT and Meeko [32] based on RDKit [26] to convert the 578 ligand files into PDBQT format. 579 To define a docking pocket location, the model structure was superimposed with the reference 580 crystal structures provided in the DUD-E set using PyMOL alignment module [33]. Then, a 581 cubic box centered at the geometrical center position of the cognate ligand structure of the 582 reference protein structure was defined. Each dimension of the box has a size of 22.5 Å, a default 583 option of the program. The parameter nrun, the number of pose generations and searches in 584 AutoDock-GPU, was set to 50 to find the optimal AutoDock score between protein and ligand. 585 586 Scoring Schemes for Ensemble Docking 587 To get the representative score of a ligand that docked to the multiple receptor structures 588 in ensemble screening, we employed two scoring schemes: AutoDock best score (ADB) and 589 Boltzmann-weighted score (BW). After gathering all AutoDock scores of a compound docked to .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 29 590 the MSM structures, ADB scheme picks the lowest value as a representative score for the ligand. 591 For example, if a ligand is docked to five kinase structures with docking scores of -11 kcal/mol, 592 -12 kcal/mol, -8 kcal/mol, -7 kcal/mol, and -9 kcal/mol, then the compound has a score of -12 593 kcal/mol. 594 Instead of using the docking score from a single structure, the BW scheme calculates a weighted 595 average of the docking scores. We modified BW score of Shin et al. [1], which was originally 596 used to calculate a score of a protein with multiple ligand conformations, to apply a single ligand 597 to multiple protein conformations (Equation 1). 598 BW Score(𝑃,𝐿) = ∑ 𝑁𝑠𝑡𝑎𝑡𝑒 𝑃𝑠𝑡𝑎𝑡𝑒 𝐴𝑢𝑡𝑜𝐷𝑜𝑐𝑘(𝑃𝑠𝑡𝑎𝑡𝑒,𝐿) × exp[ ―β × 𝐴𝑢𝑡𝑜𝐷𝑜𝑐𝑘(𝑃𝑠𝑡𝑎𝑡𝑒,𝐿)] ∑ 𝑁𝑠𝑡𝑎𝑡𝑒 𝑃𝑠𝑡𝑎𝑡𝑒 exp[ ―β × 𝐴𝑢𝑡𝑜𝐷𝑜𝑐𝑘(𝑃𝑠𝑡𝑎𝑡𝑒,𝐿)] (1) where 𝛽 = 1, P and 599 L are protein and ligand, respectively. Pstate means target protein structure with a specific state. 600 AutoDock(Pstate, L) means AutoDock score of a ligand for the target protein with a specific state. 601 602 Evaluation Metrics for Docking and Screening 603 To evaluate the performance of cognate docking, RMSDs of the docked conformations 604 from the bound conformation of the crystal structure were calculated. Then we picked two 605 conformations: one with lowest AutoDock score and the other one is the lowest RMSD 606 conformation. We also measured the docking success rate of the 24 target proteins with the 607 RMSD cutoff of 2 Å, a widely used criterion for many docking studies [21-23]. 608 In order to compare the virtual screening performance of MSM model ensemble screening with 609 X-ray crystallography and standard AF2 structures, the EFX%, AUC, and BEDROC were 610 calculated. 611 The EFX% is a widely used metric to evaluate virtual screening methods. The enrichment factor .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 30 612 quantifies the extent to which active compounds are sampled in the top N% of compounds 613 relative to the total compound set (Equation 2). 614 𝐸𝐹𝑋% = Number of actives in the top X% / Number of compounds for X% Total number of actives in the library / Total number of compounds in the library (2) 615 We set X as 1, 5, and 10. A random selection of compounds makes the EF value 1. 616 One of the most popular measures for the discrimination problem is AUC. The true positive rate 617 in relation to the false positive rate was plotted to create a receiver operating characteristic curve. 618 In the case of virtual screening, the ratio of active chemicals represents the true positive rate, 619 while the ratio of decoy molecules represents the false positive rate. When a program detects all 620 active compounds before ranking any decoy compounds, AUC reaches 1.0, the maximum value 621 and AUC 0.5 means that the program performance is the same as the random selection. 622 Although AUC gives an overall performance discriminating power between actives and decoys 623 of SBVS, it has a problem called ‘early recognition’ [25]. In virtual screening, the highly ranked 624 compounds are passed to experiment, not all compounds. Thus, it is important to rank active 625 molecules within high rank. To solve this issue, Boltzmann-enhanced discrimination of receiver 626 operating characteristic (BEDROC) puts an exponential weight on the highly ranked active 627 compounds (Equation 3). 628 𝐵𝐸𝐷𝑅𝑂𝐶 = ∑𝑁 𝑖=1 𝑒 ―𝛼𝑟𝑖/𝑁 𝑅𝑎( 1 ― 𝑒𝛼 𝑒𝛼/𝑁 ― 1) × 𝑅𝑎 sinh (𝛼 2) cosh (𝛼 2) ― cosh (𝛼 2 ― 𝛼𝑅𝑎) + 1 1 ― 𝑒 𝛼(1―𝑅𝑎) (3) 629 N is the number of compounds, Ra is the ratio of active compounds in the library, i the index of 630 the active compounds, and ri is the rank of the active compound i. In this work, the weight, , is 631 set to 20. 632 633 Acknowledgements .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 31 634 All authors acknowledge the support from the Bio & Medical Technology Development 635 Program of the National Research Foundation (NRF) funded by the Korean government (No. 636 2022M3E5F3081268). WHS also acknowledges support from Korea University Grant (No. 637 K2327351). JL also acknowledges the support from NRF Grants funded by the Korean 638 government (MSIT) (Nos. 2022R1C1C1005080 and 2020M3A9G7103933) and Korea 639 Environment Industry & Technology Institute (KEITI) through “Advanced Technology 640 Development Project for Predicting and Preventing Chemical Accidents” Program, funded by 641 Korea Ministry of Environment (MOE) (RS-2023-00219144). 642 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 32 643 References 644 1. Shin WH, Christoffer CW, Wang J, Kihara D. PL-PatchSurfer2: Improved Local Surface 645 Matching-Based Virtual Screening Method That Is Tolerant to Target and Ligand Structure 646 Variation. J Chem Inf Model. 2016;56(9):1676-91. 647 2. Bordogna A, Pandini A, Bonati L. Predicting the accuracy of protein-ligand docking on 648 homology models. J Comput Chem. 2011;32(1):81-98. 649 3. Fan H, Irwin JJ, Webb BM, Klebe G, Shoichet BK, Sali A. Molecular docking screens using 650 comparative models of proteins. J Chem Inf Model. 2009;49(11):2512-27. 651 4. Santos R, Ursu O, Gaulton A, Bento AP, Donadi RS, Bologa CG, et al. A comprehensive map 652 of molecular drug targets. Nat Rev Drug Discov. 2017;16(1):19-34. 653 5. Gaulton A, Bellis LJ, Bento AP, Chambers J, Davies M, Hersey A, et al. ChEMBL: a large- 654 scale bioactivity database for drug discovery. Nucleic Acids Res. 2012;40(Database 655 issue):D1100-7. 656 6. McClendon CL, Kornev AP, Gilson MK, Taylor SS. Dynamic architecture of a protein kinase. 657 Proc Natl Acad Sci U S A. 2014;111(43):E4623-31. 658 7. Hari SB, Merritt EA, Maly DJ. Sequence determinants of a specific inactive protein kinase 659 conformation. Chem Biol. 2013;20(6):806-15. 660 8. Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate 661 protein structure prediction with AlphaFold. Nature. 2021;596(7873):583-9. 662 9. Baek M, DiMaio F, Anishchenko I, Dauparas J, Ovchinnikov S, Lee GR, et al. Accurate 663 prediction of protein structures and interactions using a three-track neural network. Science 664 2021;373(6557):871-6. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 33 665 10. Heo L, Feig M. Multi-state modeling of G-protein coupled receptors at experimental 666 accuracy. Proteins. 2022;90(11):1873-85. 667 11. Modi V, Dunbrack RL, Jr. Defining a new nomenclature for the structures of active and 668 inactive kinases. Proc Natl Acad Sci U S A. 2019;116(14):6818-27. 669 12. Meng Y, Lin YL, Roux B. Computational study of the "DFG-flip" conformational transition 670 in c-Abl and c-Src tyrosine kinases. J Phys Chem B. 2015;119(4):1443-56. 671 13. Meng Y, Pond MP, Roux B. Tyrosine Kinase Activation and Conformational Flexibility: 672 Lessons from Src-Family Tyrosine Kinases. Acc Chem Res. 2017;50(5):1193-201. 673 14. Haldane A, Flynn WF, He P, Vijayan RS, Levy RM. Structural propensities of kinase family 674 proteins from a Potts model of residue co-variation. Protein Sci. 2016;25(8):1378-84. 675 15. Carles F, Bourg S, Meyer C, Bonnet P. PKIDB: A Curated, Annotated and Updated Database 676 of Protein Kinase Inhibitors in Clinical Trials. Molecules. 2018;23(4). 677 16. Mysinger MM, Carchia M, Irwin JJ, Shoichet BK. Directory of useful decoys, enhanced 678 (DUD-E): better ligands and decoys for better benchmarking. J Med Chem. 2012;55(14):6582- 679 94. 680 17. Chen VB, Arendall WB, III, Headd JJ, Keedy DA, Immormino RM, Kapral GJ, et al. 681 MolProbity: all-atom structure validation for macromolecular crystallography. Acta 682 Crystallographica Section D. 2010;66(1):12-21. 683 18. Kanev GK, de Graaf C, Westerman BA, de Esch IJP, Kooistra AJ. KLIFS: an overhaul after 684 the first 5 years of supporting kinase research. Nucleic Acids Res. 2021;49(D1):D562-D9. 685 19. Zhang Y, Skolnick J. Scoring function for automated assessment of protein structure 686 template quality. Proteins 2004;57(4):702-10. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 34 687 20. Santos-Martins D, Solis-Vasquez L, Tillack AF, Sanner MF, Koch A, Forli S. Accelerating 688 AutoDock4 with GPUs and Gradient-Based Local Search. J Chem Theory Comput 2021;17(2): 689 1060-73. 690 21. Shin WH, Kim JK, Kim DS, Seok C. GalaxyDock2: protein-ligand docking using beta- 691 complex and global optimization. J Comput Chem. 2013;34(30):2647-56. 692 22. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new 693 scoring function, efficient optimization, and multithreading. J Comput Chem. 2010;31(2):455- 694 61. 695 23. Lee J, Seok C. A statistical rescoring scheme for protein-ligand docking: Consideration of 696 entropic effect. Proteins. 2008;70(3):1074-83. 697 24. Adasme MF, Linnemann KL, Bolz SN, Kaiser F, Salentin, S, Haupt VJ, et al. PLIP 2021: 698 expanding the scope of the protein-ligand interaction profiler to DNA and RNA. Nucleic Acids 699 Res 2021;49(W1):W530-4. 700 25. Truchon JF, Bayly CI. Evaluating virtual screening methods: good and bad metrics for the 701 "early recognition" problem. J Chem Inf Model. 2007;47(2):488-508. 702 26. Landrum G. RDKit: Open-source cheminformatics 2006. Accessed 2022. 703 27. Irwin JJ, Shoichet BK. ZINC--a free database of commercially available compounds for 704 virtual screening. J Chem Inf Model. 2005;45(1):177-82. 705 28. Steichen JM, Kuchinskas M, Keshwani MM, Yang J, Adams JA, Taylor SS. Structural basis 706 for the regulation of protein kinase A by activation loop phosphorylation. J Biol Chem. 707 2012;287(18):14672-80. 708 29. Mirdita M, Steinegger M, Soding J. MMseqs2 desktop and local web server app for fast, 709 interactive sequence searches. Bioinformatics. 2019;35(16):2856-8. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 35 710 30. Leibschner D, Afonine PV, Baker ML, Bukóczi G, Chen VB, Croll TI, et al. Macromolecular 711 structure determination using X-rays, neutrons and electrons: recent developments in Phenix. 712 Acta Cryst. D 2019;75(10):861-877. 713 31. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, et al. AutoDock4 714 and AutoDockTools4: Automated docking with selective receptor flexibility. J Comput Chem. 715 2009;30(16):2785-91. 716 32. Forli S. Meeko 2023 [Available from: https://github.com/forlilab/Meeko. 717 33. Schrodinger, LLC. The PyMOL Molecular Graphics System, Version 1.8. 2015. 718 .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 36 719 Supporting Information 720 S1 Table. Percentage of human kinase structural states in PDB database, standard AlphaFold2, 721 and average plDDT and MolProbity of predicted structures by multi-state modeling protocol. 722 S2 Table. The number of predicted models generated by MSM protocol. 723 S3 Table. Cognate docking result of 24 DUD-E kinase proteins. The values out of the 724 parentheses are RMSD (Å) of the best AutoDock score while in the parentheses are the 725 conformation of the lowest RMSD. 726 S4 Table. Percentage of retrieved interacting residues of the crystal structures from the failed 727 cognate MSM docking benchmark cases. 728 S5 Table. Virtual screening results of individual proteins. 729 S6 Table. Details in multi-state modeling of ABL1. 730 S7 Table. Details in multi-state modeling of KDR. 731 S1 Fig. Distribution of MolProbity score by modeling method. The average MolProbity score 732 was 1.04, 1.21, 1.24 for AF2, TT and NT respectively. 733 S2 Fig. Interacting residues identified by PLIP from the crystal structure of PRKCB (PDB ID: 734 2IOE, A) and docking results with the MSM model with TT (B). The interacting residues are 735 colored as blue while the docked ligands are shown in yellow. The RMSD between the 736 conformation is 2.79 Å. The hydrophobic interaction between the molecules is represented as 737 dashed lines while the hydrogen bonding is shown as blue solid lines. 738 S3 Fig. Predicted structures of ABL1 with DFGout-BBAminus state. TT model (template: 739 7HZ0) is shown in sky blue and NT model (template: 3PYY) is shown in gold. The activation 740 loop of each model is colored as green and orange for TT model and NT model, respectively. .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint 37 741 S4 Fig. Example of modeling result that the conformational state is not matched with template. 742 A. Overall structure comparison between the modeled structure (colored in beige) and the 743 template (PDB ID: 5XY1, colored in sky blue). B. Focused view for the binding site. C. 744 Ramachandran plot for modeled structure and the template. The modeled structure were 745 annotated as DFGin-BLAminus (red x), while the template was DFGin-ABAminus (green x). .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint .CC-BY 4.0 International licenseavailable under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint