{"paper_id":"27a03e7a-1f7b-49bb-9411-34a1f1f458f7","body_text":"1\n1 Applying multi-state modeling using AlphaFold2 for \n2 kinases and its application for ensemble screening\n3\n4 Jinung Song1¶, Junsu Ha2¶, Juyong Lee1,2,3, Junsu Ko2, Woong-Hee Shin2,4*\n5 1 College of Pharmacy, Seoul National University, Seoul, Republic of Korea\n6 2 Arontier Co., Seoul, Republic of Korea\n7 3 Department of Molecular Medicine and Biopharmaceutical Sciences, Graduate School of \n8 Convergence Science and Technology, Seoul, Republic of Korea\n9 4 Department of Biomedical Informatics, Korea University College of Medicine, Seoul, Republic \n10 of Korea\n11\n12 * Corresponding author\n13 E-mail: whshin@korea.ac.kr (WHS)\n14\n15 ¶: These authors contributed equally to this work.\n16\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n2\n17 Abstract\n18 Structure-based virtual screening (SBVS) is a pivotal computational approach in drug \n19 discovery, enabling the identification of potential drug candidates within vast chemical libraries \n20 by predicting their interactions with target proteins. The SBVS relies on the receptor protein \n21 structures, making it sensitive to structural variations. Kinase, one of the major drug targets, is \n22 known as one of the typical examples of an active site conformation change caused by the type \n23 of binding inhibitors. Examination of human kinase structures shows that the majority of \n24 conformations have the DFGin state. Thus, SBVS using the structures might cause a favor of \n25 type of ligand type I inhibitors, bind to the DFGin state, rather than finding the diverse scaffolds. \n26 Recent advances in protein structure prediction, such as AlphaFold2 (AF2), offer promising \n27 solutions but may still be possibly influenced by the structural bias in existing templates.\n28 To address these challenges, we introduce a multi-state modeling (MSM) protocol for kinase \n29 structures. We apply MSM to AF2 by providing state-specific templates, allowing us to \n30 overcome structural biases and thus apply them to kinase SBVS. We benchmarked our MSM \n31 models in three categories: quality of predicted models, reproducibility of ligand binding poses, \n32 and identification of hit compounds by ensemble SBVS. The results demonstrate that MSM-\n33 generated models exhibit comparable or improved structural accuracy compared to standard AF2 \n34 models. We also show that MSM models enhance the accuracy of cognate docking, effectively \n35 capturing the interactions between kinases and their ligands.\n36 In virtual screening experiments using DUD-E compound libraries, our MSM approach \n37 consistently outperforms standard AF2 modeling. Notably, MSM-based ensemble screening \n38 excels in identifying diverse hit compounds for kinases with structurally diverse active sites, \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n3\n39 surpassing standard AF2 models. We highlight the potential of MSM in broadening the scope of \n40 kinase inhibitor discovery by facilitating the identification of chemically diverse inhibitors.\n41\n42 Author Summary\n43 One of the main problems with structure-based virtual screening is structural flexibility. \n44 Ensemble screening is one of the conventional approaches to solving the issue. Gathering \n45 experimental structures or molecular simulations could be used to compile the receptor \n46 structures. Recent developments in algorithms for predicting protein structures, like AlphaFold2, \n47 suggest that different receptor conformations could be produced. However, the prediction \n48 approaches produce biased structures because of the bias in the structure database. In order to \n49 solve the problem, we developed a protocol called multi-state modeling for kinases. Rather than \n50 supplying multiple sequence alignments as an input, we gave the AlphaFold2 a specific template \n51 structure and the sequence alignment between the template and query.\n52 Our findings imply that our technique can yield a particular structural state of interest with an \n53 enhanced or comparable structural quality to AlphaFold2 and predict highly accurate protein-\n54 ligand complex structures. Lastly, compared to the typical AlphaFold2 models, ensemble \n55 screening using the multi-state modeling approach improves the structure-based virtual screening \n56 performance, particularly for diverse active molecular scaffolds. \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n4\n57 Introduction\n58 Structure-based virtual screening (SBVS) is one of the most widely used computational \n59 drug discovery approaches to identify novel active compounds against the target from a virtual \n60 molecule library. By predicting the interaction between the ligand and the target protein, the \n61 technique ranks the ligands based on their scores that mimic the binding affinity between the \n62 molecules. The SBVS is a cost- and time-efficient way to explore the vast chemical space, by \n63 significantly narrowing down the number of drug candidates that need to be synthesized and \n64 tested in experiments. Calculating the interaction between the receptor and ligand can be done in \n65 various ways: molecular docking, molecular dynamics, fingerprint, pharmacophore matching, \n66 and so forth. As its name implies, the method requires the target protein structure, and the \n67 performance depends on the protein conformation. For targets whose experimental structures are \n68 not available, researchers should predict the three-dimensional structures using methods such as \n69 homology modeling.\n70 One of the major obstacles in the SBVS method, especially molecular docking, is caused by its \n71 static picture of a receptor structure. Proteins are flexible molecules, so they can change their \n72 shape depending on the binding partners. The structural change for receptor protein might lead to \n73 a failure of molecular docking [1-3]. One of the techniques for treating receptor flexibility is \n74 ensemble screening. The method uses pre-generated diverse receptor structures gathered from \n75 experimental structure databases or simulation trajectories such as molecular dynamics or normal \n76 mode analysis. Ligands in the screening library are docked to individual structures and ranked by \n77 their representative scores. There are various methods to get the representative score: arithmetic \n78 mean, harmonic mean, etc. Since the ensemble method uses multiple target structures, it is \n79 important to reflect the structural diversity when selecting the receptor ensemble. However, the \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n5\n80 crystal structures of the target would be biased to thermodynamically stable states or the major \n81 type of inhibitor-bound form. This might make the screening for obtaining diverse scaffolds or \n82 even lead to a failure in the SBVS even structural ensembles are used. \n83 Kinases are one of the typical examples having structural diversity. They play a key role \n84 in biological processes in phosphorylation transferring a phosphate group from ATP to proteins. \n85 Since the process modulates key cellular operations like cell cycle regulation, metabolism, and \n86 apoptosis, kinases have become attractive targets in drug discovery. According to Santos et al. \n87 [4], kinases belong to one of the four main target families, which are G protein-coupled receptors \n88 (GPCRs), ion channels, nuclear receptors, and kinases. In the ChEMBL database [5], the kinases \n89 are targeted by 14.2% of compounds.\n90 The kinase domain, a structural domain with catalytic function, has a highly conserved structure. \n91 It is composed of an N-lobe and a C-lobe, linked by a hinge region [6]. The N-lobe is structured \n92 with five β-strands and a single alpha helix, known as the C-helix, while the C-lobe is comprised \n93 of multiple alpha helices. The catalytic activity of kinase domains is regulated by two loops, \n94 namely the activation loop and the catalytic loop. The His-Arg-Asp (HRD) motif of the catalytic \n95 loop directly interacts with the hydroxyl group of the target protein residue (serine, threonine, or \n96 tyrosine) set for phosphorylation. The Asp-Phe-Gly (DFG) motif found in the N-terminal of the \n97 activation loop plays a crucial role in anchoring ATP to the active site. The conformational states \n98 of the active site of kinase domains are classified as DFGin, DFGinter, or DFGout, based on the \n99 orientation of the aspartic acid of the DFG motif relative to the site. Fig 1 illustrates two distinct \n100 structural states, DFGin (Fig 1A) and DFGout (Fig 1B), of BRAF. The DFGin state locates Phe \n101 of the DFG motif into the ATP binding pocket, allowing kinase to hold ATP, thus it is called the \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n6\n102 active state. On the other hand, DFGout conformation directs the Phe out of the ATP binding \n103 pocket, so it is classified as an inactive state.\n104\n105 Fig 1. The structure of kinase domain and key structural elements. The active state (DFGin-\n106 BLAminus) of BRAF (PDB ID: 6UAN, A) and the inactive state (DFGout-BBAminus) of BRAF \n107 (PDB ID: 8C7X, B). The active site is focused in the red box on the right of each panel. The \n108 activation loop and DFG motif changes conformation depending on the state.\n109\n110 The kinase inhibitors can be categorized depending on the binding site conformation of kinases \n111 and the place the compounds bind. Type I inhibitors are the ones that bind to the ATP-binding \n112 site and thus compete with ATP. Most type I inhibitors bind to the DFGin state. We observed \n113 that the majority of experimentally determined human kinase structures form DFGin state (87%, \n114 as of May 2023). On the other hand, type II inhibitors often associate with the kinase to the \n115 DFGout state. The type II compounds tend to occupy ATP-binding site partially and a \n116 hydrophobic pocket close to the ATP-binding site, which is opened when the activation loop \n117 forms the DFGout conformation. In general, type II inhibitors have more selectivity than type I \n118 inhibitors. According to Hari et al. [7], the type II inhibitors' selectivity is influenced by the \n119 inherent variations in kinase's capacity to adopt the DFGout conformation. Lastly, type III \n120 inhibitors, sometimes referred to as allosteric inhibitors, bind to a kinase site that is not directly \n121 connected to the ATP-binding site. Therefore, it is crucial to take into account as many different \n122 structural states of kinases as possible to find diverse hit molecules for SBVS. \n123 Recent advances in protein structure prediction using deep learning techniques, such as \n124 AlphaFold2 (AF2) [8] and RoseTTAFold [9], allow for accurate modeling of protein structures. \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n7\n125 However, these methods rely on pre-trained models from the PDB database, so the generated \n126 models might be affected by the conformational state distribution in the PDB if the protein can \n127 form diverse structures. Therefore, the methods could produce kinase structures similar to the \n128 DFGin state, and thus SBVS targeting kinases using the predicted models could potentially yield \n129 outcomes favoring type I inhibitors even using the modeled structures.\n130 To solve the issue, Heo and Feig [10] suggest a method called multi-state modeling (MSM) to \n131 predict GPCRs with high accuracy of the desired structural state. The authors found that \n132 predicted models by AF2 tend to have either active or inactive states depending on GPCR classes \n133 due to the small number of experimentally determined structures of GPCRs. Instead of providing \n134 multiple sequence alignment (MSA) as an input of AF2, the authors align a query sequence to a \n135 template sequence with the structural state of interest. The MSM technique showed an improved \n136 performance of modeling for GPCRs. Cognate docking using MSM models also showed \n137 enhanced accuracy with smaller root-mean-square-distance (RMSD) from the crystal structure.\n138 Inspired by the work, we established an MSM protocol for modeling the kinase structures \n139 to overcome the structural bias of kinases for SBVS by giving state-specific templates to AF2. \n140 All human kinase experimental structures were classified by the active site conformation using \n141 KinCoRe [11] to construct a state-specific template database. Our protocol was able to predict \n142 kinase conformations with the desired structural state with a high accuracy. Then we \n143 benchmarked cognate docking to observe the MSM models produced accurate binding modes \n144 with lower RMSD than standard AF2 models. Finally, we performed ensemble SBVS with \n145 generated models by MSM.  The ensemble method showed higher performance compared to \n146 screening with modeled structures using AF2, especially when the active molecules are diverse.\n147\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n8\n148 Results and Discussion\n149 Conformation Distributions of Experimental Kinase Structures and \n150 AlphaFold2 Predicted Models\n151 Proteins are flexible molecules, so they may change their structures when bind to their \n152 partners. This conformational change could act as an obstacle to SBVS. The active site of \n153 kinases also forms diverse conformations depending on the type of the binding compound. For \n154 example, the DFGin state binds to the type I inhibitor, while the DFGout state binds to the type II \n155 inhibitor. Dunbrack and his colleagues introduced a rule, called KinCoRe [11], to classify the \n156 kinase structure. The rule classifies structural states of kinase structures into 12 categories based \n157 on the spatial state of the activation loop and the dihedral angle of DFG motif. Details of the \n158 criteria can be found in the Materials and Methods section and Modi et al [11].\n159 Fig 2 shows a distribution of kinase conformational states annotated by KinCoRe scheme [11] in \n160 the PDB database (blue bar). Based on the KinCoRe notation, more than half (53.6%) of \n161 experimentally determined kinase structures have DFGin-BLAminus conformation. Details of \n162 the structure distribution are shown in Supplementary Material (S1 Table). Other than the major \n163 state, the other conformational states occupy less than 10% of PDB structures. \n164\n165 Fig 2. Distribution of crystal structures and standard AF2 structures for each kinase \n166 conformation. The X-axis is conformational states annotated by KinCoRe and the y-axis is the \n167 percentage of each state. All human crystal kinase structures deposited in RCSB-PDB, DUD-E \n168 target protein crystal structures, and predicted AF2 models with default parameters (standard \n169 AF2) are colored as blue, orange, and green, respectively.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n9\n170\n171 The preference for the DFGin state in experimentally determined structures might be influenced \n172 by the thermodynamic stabilities of the conformational states. Meng et al [12], studied the \n173 transition between the structural states of c-Abl and c-Src kinases using umbrella sampling and \n174 the potential of mean force. With a difference of 1.4 kcal/mol, the DFGin state of c-Abl has a \n175 lower Gibbs free energy than the DFGout state. The calculation shows the active conformation \n176 (DFGin) is the dominant population while the DFGout state occupied only 9% in their \n177 simulation. Similarly, for c-Src, the DFGin state is more favored than the DFGout state. The free \n178 energy difference between the states is calculated as 5.4 kcal/mol. The authors also studied a \n179 thermodynamic barrier of the transition between the states [13]. It is estimated as the order of 2-4 \n180 kcal/mol making the transition from highly populated the DFGin state to the DFGout state hard. \n181 Levy and his colleagues [14] collected 2,896 kinase structures and multiple sequence alignment \n182 and applied Potts model to predict structural propensities from sequences. With the statistical \n183 potential, most of the kinases are predicted to have preferences for the DFGin state. The penalty \n184 for forming the DFGout state reaches 2-3 kcal/mol for the extreme case.\n185 Not only for the thermodynamic preference, but the type of kinase inhibitor distribution might \n186 also affect the skewness of kinase structures. We counted the number of compounds in kinase \n187 inhibitor types from PKIDB (assessed August 2023) [15]. PKIDB is a curated database of kinase \n188 inhibitors in clinical trials. Out of 369 compounds in the database, only 84 molecules have their \n189 inhibition type annotation, because annotating the ligand type needs the complex structure. Type \n190 I inhibitor (DFGin bound) is the dominant form, occupying 66% (55 compounds) of the \n191 annotated inhibitors. In contrast, 17 molecules are labeled as type II (DFGout bound). Thus, \n192 DFGin conformation (active state) might have more chance to be crystallized than the DFGout \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n10\n193 state. The biased conformational states of kinases would make the discovery of chemically \n194 diverse kinase inhibitors harder. \n195 We also examined the conformational state distribution of 25 DUD-E targets [16], that are \n196 benchmarked throughout this study, in the PDB database (orange bar). DUD-E is a widely used \n197 benchmark dataset for evaluating virtual screening methods, composed of pharmaceutically \n198 important targets such as GPCRs, kinases, and nuclear receptors. Among the set, we extracted 25 \n199 kinase targets for benchmarking. The experimental structures are less biased than all PDB \n200 structures, about 30% of DUD-E proteins have DFGin-BLAminus conformation. The other \n201 states tend to have a higher proportion than the PDB structure distribution of the states. The \n202 difference in the state distribution between all kinases and DUD-E targets potentially means that \n203 the highly biased nature of kinase structures might not be suitable for finding diverse types of \n204 kinase inhibitor discovery.\n205 From an MSA of a given sequence, AF2 [8] extracts coevolution information by the MSA \n206 Transformer and predicts the three-dimensional structure based on the information and deep-\n207 learning models trained on existing protein structures. We modeled 25 kinase catalytic domains \n208 provided in the DUD-E kinase subset with the default parameters of AF2 resulting in 125 \n209 structures (five models per target), then assigned the conformational states of the models using \n210 KinCoRe. Throughout this paper, AF2 with default parameters is called standard AF2. The \n211 average plDDT and MolProbity [17] score of the predicted structures is 89.38 and 1.04, \n212 respectively. This implies that they are properly modeled. Out of 125 predicted models, 91 \n213 structures (72%) are annotated to have DFGin-BLAminus conformation (green bar). The \n214 distribution of predicted models is more skewed than all human kinase structures and the DUD-E \n215 set proteins. However, the other states have a lower proportion than the experimental structures. \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n11\n216 The standard AF2 did not produce DFGin-ABAminus, DFGin-BLBtrans, DFGinter-BABtrans, \n217 and DFGout-Unassigned conformations, which occupy a small portion of the human kinase PDB \n218 structures (7.3%, 3.0%, 0.3%, and 3.2%, respectively, S1 Table). Since the experimentally \n219 determined kinase structures are biased to the DFGin-BLAminus, predicted kinase structures by \n220 AF2 might have a chance to be biased to the conformation, which can be observed in previous \n221 research [10]. In addition to the biased trained models, the template selection process in AF2 \n222 modeling does not take the structural state of the kinase into account. \n223\n224 Predicting State-specific Kinase Structures using Multi-state \n225 Modeling Protocol\n226 AF2 with MSM protocol provides conformational state-specific structures as templates \n227 for AF2 models. All human kinase structures from the KLIFS database [18] and further \n228 classified by their structural state following KinCoRe rules [11]. For each state, the five highest \n229 sequence similar structures were selected as templates. Then AF2 was executed with the \n230 template information to predict five models for each template. The models with undesired \n231 conformational state were removed, and then the model with the highest plDDT was selected for \n232 our benchmark. The overall workflow is illustrated in Fig 3. Details of the kinase MSM protocol \n233 are elaborated in the Materials and Methods section.\n234\n235 Fig 3. Workflow of multi-state modeling of kinases\n236\n237 We prepared two template sets, a trivial template (TT) set containing 100% sequence identical \n238 structural template and a nontrivial template (NT) set not having identical protein. With our \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n12\n239 MSM protocol, AF2 was able to model structures in a specific state, producing 8.5 and 8.2 \n240 models with TT set and NT set on average, respectively. The number of predicted models for \n241 each target ranges from three (WEE1) to 11 (ABL1 and LCK, S2 Table). The average plDDT \n242 scores of structures are 89.63 and 88.08 for models from the TT set and NT set, respectively, \n243 showing comparable values to the standard AF2 predictions (89.38). As expected, TT set models \n244 have higher accuracy on average than NT set, but their difference is marginal. The distribution of \n245 the MolProbity score is given in S1 Fig. We also examined the quality of models for each \n246 structural state. The average plDDT values range from 86.62 (DFGin-BLAplus) to 94.48 \n247 (DFGin-BLAminus), meaning that the quality of the models does not depend much on the \n248 structural state of kinases. The highest plDDT score of DFGin-BLAminus might be caused by \n249 the highest frequency of the state in the PDB database (S1 Table). The Pearson’s correlation \n250 coefficient between the average plDDT and percentage in the PDB structure of each structural \n251 state is 0.739. Even though the MSM protocol provides a structural template for AF2, the model \n252 quality might be influenced by the pre-trained models of AF2, so the average plDDT of each \n253 model follows the distribution of the crystal structure. \n254 The TT set model and the NT set model have average MolProbity scores of 1.21 and 1.24, \n255 respectively. These figures represent a modest decline from the baseline AF2 models (1.04). We \n256 also examined each structural state's average MolProbity score. 1.06 (DFGin-ABAminus) is the \n257 lowest value, and 1.33 (DFGinter-Unassigned) is the highest. For DFGin-BLAminus, the highest \n258 populated states for both PDB and standard AF2, MSM models show the average MolProbity \n259 score as 1.08. The average MolProbity score and the distribution in the PDB structure have a -\n260 0.59 Pearson's correlation coefficient, indicating a weak correlation, and the high score of \n261 MolProbity might be caused by the states with the small PDB populated states.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n13\n262 To investigate the accuracy of the models, we measured the TM-Score [19] of predicted models \n263 to crystal structures given in the DUD-E set, called ‘reference structure’ throughout this paper. \n264 TM-Score assesses a structural similarity between two given proteins, ranging from zero (not \n265 similar) to one (identical). For MSM models, we used the predicted structures in the same \n266 structural state as the reference. The average TM-Scores are 0.92 and 0.90 for the models \n267 predicted using TT and NT sets, respectively. Compared with standard AF2 models (average \n268 TM-Score: 0.87), the MSM technique provided more similar models than the standard AF2 \n269 protocol, which is expected since the MSM provides structural templates. Models with TT sets \n270 generally have more accurate structures than those with NT sets as also expected.\n271 Consequently, by utilizing structures that represent a variety of states for the target kinase, the \n272 MSM is able to generate diverse structures as desired with high accuracy. Thus, it could provide \n273 a proper structure set for kinase ensemble SBVS.\n274\n275 Cognate Docking Accuracy of a Compound to the Multi-state \n276 Modelled Structures\n277 To examine whether modeled structures are suitable for molecular docking and thus \n278 structure-based virtual screening or not, we first conducted a cognate docking experiment on the \n279 predicted structures, both standard AF2 and MSM models. Ligands from complex crystal \n280 structures of the DUD-E kinase subset were used for this benchmark. For the standard AF2 \n281 model, the highest plDDT model for each protein, which is also used for performing virtual \n282 screening in the next section, was selected for evaluation. For evaluating the MSM protocol, we \n283 used the modeled structure with the same KinCoRe annotation as the reference PDB structure \n284 provided by DUD-E. Out of 25 kinases in DUD-E, IGF1R was removed, since the crystal \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n14\n285 structure was not assigned any of the structural states by KinCoRe, Therefore, the number of \n286 target proteins becomes 24. AutoDock-GPU [20] was employed to predict the 50 binding poses. \n287 The RMSDs of all predicted poses to the crystal binding pose were calculated. To analyze, we \n288 took three values: the RMSD of the best AutoDock score pose, that of the closest pose to the \n289 crystal binding mode, and the average RMSD of all 50 poses. \n290 Table 1 summarizes the docking accuracy evaluation results. Taking the best AutoDock score \n291 models, the standard AF2, MSM structures modeled with TT, and NT have average RMSD of \n292 2.74 Å, 2.15 Å, and 3.49 Å, respectively. In addition, the success cases with an RMSD cutoff of \n293 2.0 Å, a standard criterion for judging docking success [21-23], are 11 (standard AF2), 16 (with \n294 TT), and 9 (with NT) out of 24 receptors. Individual RMSD values are given in S3 Table. \n295\n296 Table 1. Docking accuracy benchmark result. The values are the average values of 24 proteins \n297 and the numbers in the parentheses are the number of success cases with RMSD < 2 Å.\nStructure Standard AF2 MSM with TT MSM with NT\nBest AutoDock \nScore\n2.74 (11) 2.15 (16) 3.49 (9)\nLowest RMSD 1.50 (20) 1.40 (19) 2.14 (14)\nAverage RMSD 2.93 2.49 3.34\n298\n299 Comparing the MSM models and AF2 models, the multi-state models with TT sets have the most \n300 accurate docking poses. As observed in the previous section, template information influenced the \n301 quality of the docking poses, i.e., the predicted docking poses of the models with TT sets have \n302 smaller RMSD and more successful cases than those with NT sets. One successful example is \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n15\n303 AKT2 (Fig 4). Although the TM-scores of MSM with TT set (0.98) and standard AF2 (0.96) are \n304 similar, the best AutoDock score models docked to the MSM using TT sets has more accurate \n305 binding pose than standard AF2, RMSD of 0.86 Å (magenta, Fig 4B) and 3.24 Å (cyan, Fig 4C), \n306 respectively.\n307\n308 Fig 4. Predicted cognate docking structures of AKT2. A. Superimposed structures of crystal \n309 structure (PDB ID: 3D0E, grey), MSM model (DFGin-BLAminus, magenta; TM-Score: 0.98), \n310 and standard AF2 (cyan; TM-Score: 0.96). B.  Predicted docking poses of the cognate ligand on \n311 the AF2 with MSM. The RMSD of the predicted pose is 0.86 Å. C. Predicted docking poses of \n312 the cognate ligand on the standard AF2. The RMSD of the predicted pose is 3.24 Å.\n313\n314 Although the docking is not successful (RMSD > 2 Å), MSM structure was able to retrieve the \n315 interaction between protein and ligand found in the crystal structure for some cases. Out of eight \n316 cases, the interacting residues in the reference structures were successfully captured more than \n317 50% in five proteins analyzed by PLIP [24] (S4 Table). For example, the predicted binding pose \n318 of the best AutoDock score conformation of PRKCB cognate docking was 2.79 Å when the \n319 MSM with TT structure was used as the receptor structure. However, out of eleven interacting \n320 residues identified in the reference structure, ten residues were retrieved in the MSM-ligand \n321 complex model. By analyzing the interaction pattern, the predicted docking pose to the multi-\n322 state modeled structure has hydrophobic interaction with L348, F353, V356, A369, K371, A483, \n323 and D484 and hydrogen bonding with T404, E421, and V423. These interactions are also \n324 observed in the crystal structure (S2 Fig). The cognate docking benchmark result would imply \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n16\n325 that MSM models could be used to predict binding poses of kinase-ligand complexes, thus they \n326 are suitable for virtual ensemble screening.\n327 Similar trends are observed for the lowest RMSD conformations and average RMSD of 50 \n328 conformations. Among the AF2 predicted structures, MSM with TT shows the most accurate \n329 models (average of the lowest RMSD pose: 1.40 Å, success case: 19) and the standard AF2 \n330 performs slightly lower (average of the lowest RMSD pose: 1.50 Å, success case: 20). The \n331 docking accuracy of MSM with NT sets is the worst among the receptor structure sets (average \n332 of the lowest RMSD pose: 2.14 Å, success case: 14). In most cases, the best AutoDock score \n333 conformations do not match to the lowest RMSD conformations, which means that the \n334 AutoDock score could not be able to find the optimal docking poses. The average RMSD of 50 \n335 conformations follows the same order as the other two metrics: MSM with TT is the smallest, \n336 and MSM with NT is the highest.\n337\n338 Virtual Screening Performance with Multi-state Models\n339 To investigate the advantage of using the MSM technique for SBVS, the compound \n340 library for each kinase protein from DUD-E docked to structures with diverse states generated by \n341 using our method and compared to the structures modeled with standard AF2. AutoDock-GPU \n342 was used for the benchmark. For ensemble docking using generated models by MSM, since a \n343 molecule was docked to a couple of receptor structures and each compound-structure pair had a \n344 docking score, we needed to decide representative scores of the compounds to rank the \n345 molecules. We employed two types of representative scores: the AutoDock best (ADB) and the \n346 Boltzmann-weighted (BW) scores. The ADB score picks the lowest AutoDock score among the \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n17\n347 docked results, while the BW score is a weighted average of the scores across all structures. The \n348 details of the BW score are given in the Materials and Methods section.\n349 Table 2 shows the performance of SBVS using the various receptor structure sets. To evaluate \n350 the performance, five metrics are used: enrichment factors at 1%, 5%, and 10% (EFX%), an area \n351 under ROC curve (AUC), and Boltzmann-enhanced discrimination of receiver operating \n352 characteristic (BEDROC). EF at X% indicates a capability of finding active molecules within the \n353 top X%, and AUC represents the discrimination power of a screening method between active and \n354 decoy molecules. BEDROC puts exponential weights to the early rank of molecules, thus it is \n355 able to solve the ‘early recognition problem’ caused in AUC [25]. Details of the metrics are \n356 illustrated in the Materials and Methods section. As observed in the cognate docking benchmark, \n357 MSM structures are better or equal to standard AF2 results, regardless of the scoring method or \n358 the template set for modeling. Details of individual results are provided in S5 Table. Comparing \n359 EF1% values target-by-target, the MSM performed better than or equal to standard AF2 in 18 \n360 proteins out of 25 targets (72%) using the ADB score screened to the structures modeled with the \n361 NT set. Even with the lowest EF1% combination, TT models with BW scoring, the MSM \n362 performed better than or equal to standard AF2 in more than half of proteins (13 proteins). \n363 Interestingly, although the TT set models with the same KinCoRe classification as the reference \n364 have more accurate structure and docking poses than the NT set models, the virtual screening \n365 performance is slightly worse in both scoring schemes. \n366\n367 Table 2. Performance of structure-based virtual screening on various receptor models.\nMSM with TT MSM with NT\nADB BW ADB BW\nStandard \nAF2\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n18\nDUD-E Kinase Subset (25 proteins)\nEF1% 7.2 7.2 8.2 8.0 6.6\nEF5% 3.5 3.7 3,7 3.7 3.5\nEF10% 2.7 2.7 2.6 2.7 2.6\nAUC 0.660 0.639 0.661 0.664 0.659\nBEDROC 0.193 0.196 0.183 0.189 0.183\nAverage Dissimilarity ≥ 0.7 (14 proteins)\nEF1% 5.3 5.6 6.4 6.9 4.6\nEF5% 3.0 3.1 3.2 3.2 2.8\nEF10% 2.4 2.4 2.2 2.3 2.1\nAUC 0.646 0.652 0.648 0.650 0.641\nBEDROC 0.166 0.171 0.154 0.161 0.146\nAverage Dissimilarity < 0.7 (11 proteins)\nEF1% 9.6 9.3 10.6 9.3 9.1\nEF5% 4.2 4.4 4.4 4.3 4.3\nEF10% 3.0 3.1 3.1 3.1 3.3\nAUC 0.678 0.681 0.678 0.681 0.681\nBEDROC 0.228 0.228 0.220 0.226 0.229\n368\n369 To elucidate the variance in performance between MSMs using TT and NT models, we \n370 identified kinases that exhibited notable differences in EF1% between the two sets. Of the kinases \n371 studied, both ABL1 and KDR demonstrated superior performance using the MSM with the NT \n372 set compared to the TT set, across both ensemble scoring methods (S5 Table). We measured \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n19\n373 TM-Scores of the individual models with the same structural state, those of templates used to \n374 predict the structures, and EF1% of the models (S6 and S7 Tables). Among the kinase states \n375 analyzed, the DFGout-BBAminus for both proteins produced a significant difference in both \n376 templates and models. When the kinase forms an active state, the activation loop forms an \n377 extended conformation to facilitate the catalytic function of the protein, thus the DFGin \n378 conformation has a conserved structure. On the other hand, for the DFGout conformation, the \n379 activation loop collapsed on the protein surface with high flexibility [11]. Therefore, although \n380 the protein structures have the same KinCoRe notation with DFGout, the activation loop can \n381 have different structures. The TM-Scores of the templates to model the DFGout-BBAminus state \n382 are 0.88 and 0.96 for ABL1 and KDR, respectively. The structural difference of templates \n383 influenced the predicted models, resulting in TM-Scores 0.89 for both proteins. S3 Fig shows a \n384 structural difference of ABL1 predicted models with DFGout-BBAminus conformation. The \n385 difference of activation loop location also might affect to the virtual screening performance, \n386 leading TT set (3.28 and 8.78 for ABL1 and KDR, respectively) has lower EF1% value than NT \n387 set (15.83 and 15.85 for ABL1 and KDR, respectively). Significantly, this difference of \n388 performance in the DFGout-BBAminus state impacted the overall ensemble docking results.\n389 One benefit of using the MSM models is that the predicted models are diverse, so it is potentially \n390 useful for discovering various scaffolds of hit chemicals. To examine whether the hypothesis is \n391 true or not, we divided the DUD-E kinase subset into two based on the diversity of active \n392 compounds in the screening library. The pairwise Tanimoto coefficients (Tc) between the active \n393 compounds were calculated by using RDKitFP fingerprint in RDKit [26]. Then the pairwise \n394 distances between the compounds are calculated as (1 – Tc). The diversity of active compounds \n395 is defined as the average distances of the active compounds (S5 Table). With a threshold of 0.7, \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n20\n396 the kinases were classified into two groups, 14 proteins with higher than or equal to the cutoff \n397 and the remaining 11 targets. \n398 For the kinases with more diverse active compounds, the ensemble screening with MSM models \n399 showed higher performance than standard AF2 models in all metrics. For EF1%, the MSM \n400 performed much better than standard AF2 (Table 2). Out of 14 proteins, the MSM with NT set \n401 and BW scoring has higher or equal EF1% values than standard AF2 models in 11 targets. On the \n402 other hand, when the active compounds become less diverse, MSM still performed better than \n403 the standard AF2, but the gap between them declined for EF1%. This implies that our approach \n404 would be powerful for discovering diverse molecular scaffolds. \n405 One of the diverse active hit examples is ABL1. The average dissimilarity of the active \n406 compounds is 0.73. Regardless of the template set and the scoring scheme, EF1% of MSM \n407 ensemble screening showed higher performance (TT models with ADB: 9.3, TT models with \n408 BW: 10.4, NT models with ADB: 15.3, and NT models with BW: 17.5) than standard AF2 \n409 model (6.0). We also examined the diversity of active compounds ranked within the top 1% \n410 ranked molecules. Our ensemble protocol tends to find diverse hits: 0.51, 0.53, 0.63, and 0.62 for \n411 TT models with ADB, TT models with BW, NT models with ADB, and NT models with BW, \n412 respectively. In contrast, the diversity of active compounds using the standard AF2 is 0.27.  \n413 CSF1R is another example showing MSM ensemble screening is able to find diverse hit \n414 compounds. The average dissimilarity of hit molecules is 0.73. The gap of EF1% between MSM \n415 models (5.4) and the standard AF2 model (3.6) is smaller than the case of ABL1. The diversities \n416 of discovered hit compounds within the top 1% are 0.70, 0.65, 0.73, and 0.73 in the order of TT \n417 models with ADB, TT models with BW, NT models with ADB, and NT models with BW, \n418 similar to the average dissimilarity of all active compounds. However, the hit compounds within \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n21\n419 the top 1 % identified using the standard AF2 model is 0.23. Fig 5 shows the docking pose of \n420 one of the active compounds, ChEMBL245377. The MSM ranked the active compound within \n421 the top 1% (10th using TT models and BW scoring, out of 12316 compounds) which is not \n422 ranked by the standard AF2 structure (167th). The compound is not successfully docked using the \n423 the standard AF2 model.\n424\n425 Fig 5. Predicted binding poses of ChEMBL245377 and CSF1R structures. A. The \n426 superimposed complex structures of ChEMBL245377 and CSF1R. The crystal structure of \n427 CSF1R (PDB ID:3KRJ), standard AF2 predicted structure, and MSM model (DFGout-\n428 BBAminus, with TT sets) are represented as ribbon diagrams colored as gold, pink, and orange, \n429 respectively, while the predicted docking poses of the compound are represented as sticks. B,C. \n430 Focused binding sites and docking poses of ChEMBL245377. The compound is ranked 10th by \n431 ensemble docking (B) and 167th in docking for the standard AF2 predicted structure (C), \n432\n433 Although the ensemble screening with MSM structures performed generally better than standard \n434 AF2, there is an issue with selecting the representative score of a compound. For instance, for \n435 FGFR1, the MSM model with TT set and BW score showed EF1% as 2.1. When we observed \n436 EF1% values of individual models, the DFGin-BLBplus state outperformed than any other \n437 structures including standard AF2 (Table 3). Thus, a proper method for selecting or calculating \n438 representative scores for a compound should be designed to achieve high performance for MSM \n439 ensemble screening.\n440\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n22\n441 Table 3. EF1% for FGFR1 structures of specific states. The MSM structures are built from TT \n442 set.\nStructural State EF1%\nDFGin-BLAminus 2.9\nDFGin-BLAplus 6.4\nDFGin-BLBminus 0.7\nDFGin-BLBplus 8.6\nDFGin-BLBtrans 2.1\nDFGout-BBAminus 3.6\nDFGout-Unassigned 1.4\nMSM with BW score 2.1\nStandard AF2 1.4\n443  \n444 Conclusion\n445 The receptor conformation affects SBVS performance. Like other proteins, kinase adjusts \n446 the conformation of its binding site in response to the binding ligand. Therefore, it is crucial to \n447 have adequate kinase structure to obtain inhibitors with the required mode of action or diversity. \n448 For human kinase structures that were identified through experiments, however, there is a clear \n449 bias toward the active state. The prediction of the AF2 structure could be influenced by the bias \n450 in the PDB database. We noticed that there is a bias toward the active state in the predicted \n451 structures with standard AF2 protocol. The results of the SBVS using the predicted structure \n452 would be compromised by this bias in receptor structure. To overcome the bias, we applied the \n453 MSM technique by providing a structural template to AF2 to generate structures with diverse \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n23\n454 states and using the models for ensemble docking. Compared with standard AF2 models, MSM \n455 protocol produced more accurate or comparable models although it did not give MSA as an \n456 input. Also, in cognate docking study, MSM models provided ligand docking poses close to the \n457 crystal structure. With the diverse predicted kinase structure, we performed ensemble screening. \n458 The ensemble method showed enhanced or comparable SBVS performance to the standard AF2 \n459 modeled structure result. We also observed that our method would be more suitable when the \n460 ligands are diverse, leading to the identification of a diverse range of kinase inhibitors. Even for \n461 the targets that ensemble docking method could not find active compounds, we found that some \n462 of the models outperformed the standard AF2. Thus, the selection of a representative structure \n463 should be improved and remains the next work for this project. \n464 Ensemble screening with MSM models would open the possibility of uncovering novel \n465 kinase inhibitors with diverse chemical scaffolds. It has advantages in addressing current \n466 challenges in kinase inhibitor development for finding chemically diverse compounds. The \n467 chemical diversity of kinase inhibitors could aid in overcoming the problem of drug resistance \n468 generally caused by the mutation, a significant obstacle in kinase-targeted cancer therapies. \n469 Additionally, it could increase chance to find hit compounds, not similar to the existing patents. \n470 By exploring a diverse array of kinase inhibitors with accurately predicted structures, we would \n471 be able to find inhibitors with different modes of action. This could lead to the development of \n472 novel therapeutic strategies that are more robust in the face of drug resistance. Hence, our \n473 approach could potentially be applied to more effective and precise kinase-targeted therapies. In \n474 addition, the MSM method could be applied to important therapeutic targets such as GPCRs and \n475 nuclear receptors. \n476\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n24\n477 Materials and Methods\n478 Benchmark Dataset\n479 To investigate whether our approach can generate a structural ensemble properly and \n480 improve SBVS performance or not, we selected a kinase subset of DUD-E [16]. For each target, \n481 a reference PDB structure for screening and a compound library composed of active and decoy \n482 molecules are provided in the DUD-E set. The active compounds of DUD-E set are composed of \n483 the molecules with affinity 1 μM or better were extracted from ChEMBL09 [5]. The reference \n484 structures were selected by considering the resolution and enrichment of finding active \n485 molecules by DOCK3.5. The decoys of the DUD-E set are constructed by gathering compounds \n486 with similar characteristic to the active molecules such as logP and number of rotatable bonds \n487 from ZINC database [27]. The kinase subset, which is used in this work, consists of 26 kinases \n488 with include 205.6 actives and 12,830 decoys on average. Among the 26 targets in DUD-E \n489 kinase subset, SRC kinase was removed from DUD-E benchmark set because the given reference \n490 structure is not originated from human.\n491\n492 Kinase Structural State Annotation\n493 In the active site of protein kinases, the activation loop, 20-30 residues long, is the most \n494 important secondary structural element [28] for determining the structural state. The loop starts \n495 from the conserved three-residue-long sequence, the DFG motif. In this work, the standalone \n496 version of KinCoRe [11] (https://github.com/vivekmodi/Kincore-standalone, Accessed \n497 4/14/2022) was employed to annotate the conformational state of all experimental and modeled \n498 kinase structures. The program categorizes the conformational state into 12 classes by the \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n25\n499 location of the activation loop and dihedral angles of the DFG motif. The spatial state of the \n500 activation loop is defined by two distances: 1) a distance between Phe-ring of the DFG motif and \n501 C atom of the fourth residue from the conserved Glu in the C-helix of N-lobe and 2) a distance \n502 from the Phe-ring to the conserved C atom of conserved Lys in 3 strand of N-lobe (Fig 1). \n503 Based on the distances, the activation loop location is classified into three classes: DFGin, which \n504 is the Phe-ring located under C-helix, DFGout, the Phe-ring is moved into ATP binding pocket, \n505 and DFGinter, an intermediate state between DFGin and DFGout. The program further classifies \n506 the activation loop structural state by calculating dihedral angles: ,  backbone dihedral angles \n507 of X-DFG (a residue before the DFG motif), Asp, and Phe of DFG motif, and 1 angle of DFG-\n508 Phe. As a result, DFGin, the dominant class, has seven subclasses (BLAminus, BLAplus, \n509 ABAminus, BLBminus, BLBplus, BLBtrans, and Unassigned), while DFGinter (BABtrans and \n510 Unassigned) and DFGout (BBAminus and Unassigned) have only two subclasses. The three \n511 letters after the activation loop states follows the region of Ramachandran map occupied by X, \n512 D, and F residues: A, B, L for alpha, beta, and left-handed, respectively. The 1 angle of Phe is \n513 indicated as plus (+60 degree), minus (-60 degree), and trans (180 degree). The last class is \n514 Unassigned-Unassigned, the activation loop and DFG conformations cannot be determined. \n515\n516 Construction of Template Database for Each Structural State\n517 To construct the structural template database for MSM, KLIFS [18], a database of \n518 experimentally determined kinase structures was used. The database contains catalytic domain \n519 structures of kinases, extracted from PDB and their inhibitors and provides the interaction \n520 information between the protein and the compound. As of May 2023, the database is composed \n521 of 6,344 structures (13,382 monomers). Among the kinase structures in KLIFS (Accessed \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n26\n522 1/18/2023), we filtered out the non-human proteins and proteins that produced errors during \n523 KinCoRe annotation, resulting in 11,106 monomer structures. To construct the template structure \n524 database for each state, the crystal structures with the same annotation by KinCoRe were \n525 gathered.\n526\n527 Standard AlphaFold2 Modeling\n528 The kinase structure was modeled with the standard protocol of AF2 (v.2.3.1) to compare \n529 with MSM models. Only the kinase domain was modeled from the full sequence of a protein. \n530 The MSA for the kinase sequence was generated from BFD (v.3.2019), MGnify (v.5.2022), and \n531 UniRef90 (v.2.2022) via HHblits (from HH-suite v3.3.0) and Jackhammer (from HMMER \n532 v3.3.2). Four highest sequence identity proteins with 3D atomic coordinates were selected to \n533 provide template structures. The number of recycles was set to three and the model relaxation \n534 step was integrated into the procedure. As AF2 has five different trained models and runs all of \n535 them independently in a single run, five structures were generated from a single run. The models \n536 with the highest plDDT score out of the five predicted structures for the comparison, since \n537 plDDT is a confidence measure of AF2 predicted models.\n538\n539 Multistate Modeling of Kinase using Structural Template\n540 The workflow of MSM is given in Fig 3. From a given target kinase sequence to be \n541 modeled, the templates were searched by MMseqs2 (release 11) easy-search (e-value cutoff: 1e-\n542 3) [29] against all sequences in each structural state. For each structural state, the top five \n543 templates, which were determined by e-value, were used for the modeling. To mimic the real \n544 drug discovery process and check the influence of the template for virtual screening, we \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n27\n545 generated two template sets for modeling: one set has the query protein which means 100% \n546 sequence identity with the query sequence (TT set), and the other template sets without the query \n547 protein (NT set). Modeling with each template was conducted independently.\n548 Since AF2 produced five models per single run, 25 structures for each specific state were \n549 generated in total. Among them, we selected one structure for benchmarking with two criteria: \n550 conformational state and quality of the predicted model. First, we filtered out the models with \n551 different KinCoRe annotations from the template structure classification. For example, to predict \n552 a model the DFGin-ABAminus conformation of ABL1, five templates structure of TYR family \n553 (LYN: 5XY1, EPHA2: 7KJB, 5NK3, 4TRL, and IGF1R: 3F5P) were selected. However, the \n554 models were annotated as DFGin-BLAminus rather than DFGin-ABAminus conformation (S4 \n555 Fig). Thus, all predicted models were discarded. After filtering by the KinCoRe annotation, the \n556 models with plDDT less than 70 were also removed. Among the remained structures, the model \n557 with the highest plDDT score was finally selected.  Other details for modeling are the same as \n558 standard AF2 modeling. \n559\n560 Assessment of Model Quality\n561 The TM-score [19] evaluates the structural similarity of protein structures. It is scaled \n562 according to the size of the protein and exhibits a better sensitivity to the overall structural \n563 alignment compared to the RMSD. The models with the same KinCoRe notation as the crystal \n564 structures, extracted from the DUD-E set were chosen for comparison with the crystal structure \n565 employing the TM-score. The TM-score's values range from 0 to 1, where 1 signifies a perfect \n566 alignment.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n28\n567 MolProbity [17] is a comprehensive validation tool for the structural integrity of proteins and \n568 nucleic acids. MolProbity validates protein structures through hydrogen replacement, \n569 comprehensive all-atom contact analysis, and evaluation of torsional angles. The MolProbity \n570 score integrates various factors into a single metric to indicate the model's reliability, where a \n571 lower score denotes a better model. The models from the standard AF2 and MSM were validated \n572 using MolProbity implemented in Phenix software [30].\n573\n574 Docking and Virtual Screening using AutoDock-GPU\n575 AutoDock-GPU 1.5.3 [20], which is open-source and GPU-accelerated, was employed to \n576 benchmark the virtual screening performance of kinase structures. We used AutoDockTools [31] \n577 to convert the receptor PDB files to PDBQT and Meeko [32] based on RDKit [26] to convert the \n578 ligand files into PDBQT format.\n579 To define a docking pocket location, the model structure was superimposed with the reference \n580 crystal structures provided in the DUD-E set using PyMOL alignment module [33]. Then, a \n581 cubic box centered at the geometrical center position of the cognate ligand structure of the \n582 reference protein structure was defined. Each dimension of the box has a size of 22.5 Å, a default \n583 option of the program. The parameter nrun, the number of pose generations and searches in \n584 AutoDock-GPU, was set to 50 to find the optimal AutoDock score between protein and ligand.\n585\n586 Scoring Schemes for Ensemble Docking\n587 To get the representative score of a ligand that docked to the multiple receptor structures \n588 in ensemble screening, we employed two scoring schemes: AutoDock best score (ADB) and \n589 Boltzmann-weighted score (BW). After gathering all AutoDock scores of a compound docked to \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n29\n590 the MSM structures, ADB scheme picks the lowest value as a representative score for the ligand. \n591 For example, if a ligand is docked to five kinase structures with docking scores of  -11 kcal/mol, \n592 -12 kcal/mol, -8 kcal/mol, -7 kcal/mol, and -9 kcal/mol, then the compound has a score of -12 \n593 kcal/mol.\n594 Instead of using the docking score from a single structure, the BW scheme calculates a weighted \n595 average of the docking scores. We modified BW score of Shin et al. [1], which was originally \n596 used to calculate a score of a protein with multiple ligand conformations, to apply a single ligand \n597 to multiple protein conformations (Equation 1).\n598 BW Score(𝑃,𝐿) =\n∑\n𝑁𝑠𝑡𝑎𝑡𝑒\n𝑃𝑠𝑡𝑎𝑡𝑒\n𝐴𝑢𝑡𝑜𝐷𝑜𝑐𝑘(𝑃𝑠𝑡𝑎𝑡𝑒,𝐿) × exp[ ―β × 𝐴𝑢𝑡𝑜𝐷𝑜𝑐𝑘(𝑃𝑠𝑡𝑎𝑡𝑒,𝐿)]\n∑\n𝑁𝑠𝑡𝑎𝑡𝑒\n𝑃𝑠𝑡𝑎𝑡𝑒\nexp[ ―β × 𝐴𝑢𝑡𝑜𝐷𝑜𝑐𝑘(𝑃𝑠𝑡𝑎𝑡𝑒,𝐿)]   (1) where 𝛽 = 1, P and \n599 L are protein and ligand, respectively. Pstate means target protein structure with a specific state. \n600 AutoDock(Pstate, L) means AutoDock score of a ligand for the target protein with a specific state. \n601\n602 Evaluation Metrics for Docking and Screening\n603 To evaluate the performance of cognate docking, RMSDs of the docked conformations \n604 from the bound conformation of the crystal structure were calculated. Then we picked two \n605 conformations: one with lowest AutoDock score and the other one is the lowest RMSD \n606 conformation. We also measured the docking success rate of the 24 target proteins with the \n607 RMSD cutoff of 2 Å, a widely used criterion for many docking studies [21-23].\n608 In order to compare the virtual screening performance of MSM model ensemble screening with \n609 X-ray crystallography and standard AF2 structures, the EFX%, AUC, and BEDROC were \n610 calculated. \n611 The EFX% is a widely used metric to evaluate virtual screening methods. The enrichment factor \n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n30\n612 quantifies the extent to which active compounds are sampled in the top N% of compounds \n613 relative to the total compound set (Equation 2).\n614 𝐸𝐹𝑋% =\nNumber of actives in the top X% / Number of compounds for X%\nTotal number of actives in the library / Total number of compounds in the library  (2)\n615 We set X as 1, 5, and 10. A random selection of compounds makes the EF value 1.\n616 One of the most popular measures for the discrimination problem is AUC. The true positive rate \n617 in relation to the false positive rate was plotted to create a receiver operating characteristic curve. \n618 In the case of virtual screening, the ratio of active chemicals represents the true positive rate, \n619 while the ratio of decoy molecules represents the false positive rate. When a program detects all \n620 active compounds before ranking any decoy compounds, AUC reaches 1.0, the maximum value \n621 and AUC 0.5 means that the program performance is the same as the random selection.\n622 Although AUC gives an overall performance discriminating power between actives and decoys \n623 of SBVS, it has a problem called ‘early recognition’ [25]. In virtual screening, the highly ranked \n624 compounds are passed to experiment, not all compounds. Thus, it is important to rank active \n625 molecules within high rank. To solve this issue, Boltzmann-enhanced discrimination of receiver \n626 operating characteristic (BEDROC) puts an exponential weight on the highly ranked active \n627 compounds (Equation 3). \n628 𝐵𝐸𝐷𝑅𝑂𝐶 =\n∑𝑁\n𝑖=1 𝑒\n―𝛼𝑟𝑖/𝑁\n𝑅𝑎(\n1 ― 𝑒𝛼\n𝑒𝛼/𝑁 ― 1) ×\n𝑅𝑎 sinh (𝛼\n2)\ncosh (𝛼\n2) ― cosh (𝛼\n2 ― 𝛼𝑅𝑎) +\n1\n1 ― 𝑒\n𝛼(1―𝑅𝑎) (3)\n629 N is the number of compounds, Ra is the ratio of active compounds in the library, i the index of \n630 the active compounds, and ri is the rank of the active compound i. In this work, the weight, , is \n631 set to 20.\n632\n633 Acknowledgements\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n31\n634 All authors acknowledge the support from the Bio & Medical Technology Development \n635 Program of the National Research Foundation (NRF) funded by the Korean government (No. \n636 2022M3E5F3081268). WHS also acknowledges support from Korea University Grant (No. \n637 K2327351). JL also acknowledges the support from NRF Grants funded by the Korean \n638 government (MSIT) (Nos. 2022R1C1C1005080 and 2020M3A9G7103933) and Korea \n639 Environment Industry & Technology Institute (KEITI) through “Advanced Technology \n640 Development Project for Predicting and Preventing Chemical Accidents” Program, funded by \n641 Korea Ministry of Environment (MOE) (RS-2023-00219144).\n642\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n32\n643 References\n644 1. Shin WH, Christoffer CW, Wang J, Kihara D. PL-PatchSurfer2: Improved Local Surface \n645 Matching-Based Virtual Screening Method That Is Tolerant to Target and Ligand Structure \n646 Variation. J Chem Inf Model. 2016;56(9):1676-91.\n647 2. Bordogna A, Pandini A, Bonati L. Predicting the accuracy of protein-ligand docking on \n648 homology models. J Comput Chem. 2011;32(1):81-98.\n649 3. Fan H, Irwin JJ, Webb BM, Klebe G, Shoichet BK, Sali A. Molecular docking screens using \n650 comparative models of proteins. J Chem Inf Model. 2009;49(11):2512-27.\n651 4. Santos R, Ursu O, Gaulton A, Bento AP, Donadi RS, Bologa CG, et al. A comprehensive map \n652 of molecular drug targets. Nat Rev Drug Discov. 2017;16(1):19-34.\n653 5. Gaulton A, Bellis LJ, Bento AP, Chambers J, Davies M, Hersey A, et al. ChEMBL: a large-\n654 scale bioactivity database for drug discovery. Nucleic Acids Res. 2012;40(Database \n655 issue):D1100-7.\n656 6. McClendon CL, Kornev AP, Gilson MK, Taylor SS. Dynamic architecture of a protein kinase. \n657 Proc Natl Acad Sci U S A. 2014;111(43):E4623-31.\n658 7. Hari SB, Merritt EA, Maly DJ. Sequence determinants of a specific inactive protein kinase \n659 conformation. Chem Biol. 2013;20(6):806-15.\n660 8. Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate \n661 protein structure prediction with AlphaFold. Nature. 2021;596(7873):583-9.\n662 9. Baek M, DiMaio F, Anishchenko I, Dauparas J, Ovchinnikov S, Lee GR, et al. Accurate \n663 prediction of protein structures and interactions using a three-track neural network. Science \n664 2021;373(6557):871-6.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n33\n665 10. Heo L, Feig M. Multi-state modeling of G-protein coupled receptors at experimental \n666 accuracy. Proteins. 2022;90(11):1873-85.\n667 11. Modi V, Dunbrack RL, Jr. Defining a new nomenclature for the structures of active and \n668 inactive kinases. Proc Natl Acad Sci U S A. 2019;116(14):6818-27.\n669 12. Meng Y, Lin YL, Roux B. Computational study of the \"DFG-flip\" conformational transition \n670 in c-Abl and c-Src tyrosine kinases. J Phys Chem B. 2015;119(4):1443-56.\n671 13. Meng Y, Pond MP, Roux B. Tyrosine Kinase Activation and Conformational Flexibility: \n672 Lessons from Src-Family Tyrosine Kinases. Acc Chem Res. 2017;50(5):1193-201.\n673 14. Haldane A, Flynn WF, He P, Vijayan RS, Levy RM. Structural propensities of kinase family \n674 proteins from a Potts model of residue co-variation. Protein Sci. 2016;25(8):1378-84.\n675 15. Carles F, Bourg S, Meyer C, Bonnet P. PKIDB: A Curated, Annotated and Updated Database \n676 of Protein Kinase Inhibitors in Clinical Trials. Molecules. 2018;23(4).\n677 16. Mysinger MM, Carchia M, Irwin JJ, Shoichet BK. Directory of useful decoys, enhanced \n678 (DUD-E): better ligands and decoys for better benchmarking. J Med Chem. 2012;55(14):6582-\n679 94.\n680 17. Chen VB, Arendall WB, III, Headd JJ, Keedy DA, Immormino RM, Kapral GJ, et al. \n681 MolProbity: all-atom structure validation for macromolecular crystallography. Acta \n682 Crystallographica Section D. 2010;66(1):12-21.\n683 18. Kanev GK, de Graaf C, Westerman BA, de Esch IJP, Kooistra AJ. KLIFS: an overhaul after \n684 the first 5 years of supporting kinase research. Nucleic Acids Res. 2021;49(D1):D562-D9.\n685 19. Zhang Y, Skolnick J. Scoring function for automated assessment of protein structure \n686 template quality. Proteins 2004;57(4):702-10.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n34\n687 20. Santos-Martins D, Solis-Vasquez L, Tillack AF, Sanner MF, Koch A, Forli S. Accelerating \n688 AutoDock4 with GPUs and Gradient-Based Local Search. J Chem Theory Comput 2021;17(2): \n689 1060-73.\n690 21. Shin WH, Kim JK, Kim DS, Seok C. GalaxyDock2: protein-ligand docking using beta-\n691 complex and global optimization. J Comput Chem. 2013;34(30):2647-56.\n692 22. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new \n693 scoring function, efficient optimization, and multithreading. J Comput Chem. 2010;31(2):455-\n694 61.\n695 23. Lee J, Seok C. A statistical rescoring scheme for protein-ligand docking: Consideration of \n696 entropic effect. Proteins. 2008;70(3):1074-83.\n697 24. Adasme MF, Linnemann KL, Bolz SN, Kaiser F, Salentin, S, Haupt VJ, et al. PLIP 2021: \n698 expanding the scope of the protein-ligand interaction profiler to DNA and RNA. Nucleic Acids \n699 Res 2021;49(W1):W530-4.\n700 25. Truchon JF, Bayly CI. Evaluating virtual screening methods: good and bad metrics for the \n701 \"early recognition\" problem. J Chem Inf Model. 2007;47(2):488-508.\n702 26. Landrum G. RDKit: Open-source cheminformatics 2006. Accessed 2022.\n703 27. Irwin JJ, Shoichet BK. ZINC--a free database of commercially available compounds for \n704 virtual screening. J Chem Inf Model. 2005;45(1):177-82.\n705 28. Steichen JM, Kuchinskas M, Keshwani MM, Yang J, Adams JA, Taylor SS. Structural basis \n706 for the regulation of protein kinase A by activation loop phosphorylation. J Biol Chem. \n707 2012;287(18):14672-80.\n708 29. Mirdita M, Steinegger M, Soding J. MMseqs2 desktop and local web server app for fast, \n709 interactive sequence searches. Bioinformatics. 2019;35(16):2856-8.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n35\n710 30. Leibschner D, Afonine PV, Baker ML, Bukóczi G, Chen VB, Croll TI, et al. Macromolecular \n711 structure determination using X-rays, neutrons and electrons: recent developments in Phenix. \n712 Acta Cryst. D 2019;75(10):861-877.\n713 31. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, et al. AutoDock4 \n714 and AutoDockTools4: Automated docking with selective receptor flexibility. J Comput Chem. \n715 2009;30(16):2785-91.\n716 32. Forli S. Meeko 2023 [Available from: https://github.com/forlilab/Meeko.\n717 33. Schrodinger, LLC. The PyMOL Molecular Graphics System, Version 1.8. 2015.\n718\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n36\n719 Supporting Information\n720 S1 Table. Percentage of human kinase structural states in PDB database, standard AlphaFold2, \n721 and average plDDT and MolProbity of predicted structures by multi-state modeling protocol.\n722 S2 Table. The number of predicted models generated by MSM protocol.\n723 S3 Table. Cognate docking result of 24 DUD-E kinase proteins. The values out of the \n724 parentheses are RMSD (Å) of the best AutoDock score while in the parentheses are the \n725 conformation of the lowest RMSD.\n726 S4 Table. Percentage of retrieved interacting residues of the crystal structures from the failed \n727 cognate MSM docking benchmark cases. \n728 S5 Table. Virtual screening results of individual proteins.\n729 S6 Table. Details in multi-state modeling of ABL1. \n730 S7 Table. Details in multi-state modeling of KDR.\n731 S1 Fig. Distribution of MolProbity score by modeling method. The average MolProbity score \n732 was 1.04, 1.21, 1.24 for AF2, TT and NT respectively.\n733 S2 Fig. Interacting residues identified by PLIP from the crystal structure of PRKCB (PDB ID: \n734 2IOE, A) and docking results with the MSM model with TT (B). The interacting residues are \n735 colored as blue while the docked ligands are shown in yellow. The RMSD between the \n736 conformation is 2.79 Å. The hydrophobic interaction between the molecules is represented as \n737 dashed lines while the hydrogen bonding is shown as blue solid lines.\n738 S3 Fig. Predicted structures of ABL1 with DFGout-BBAminus state. TT model (template: \n739 7HZ0) is shown in sky blue and NT model (template: 3PYY) is shown in gold. The activation \n740 loop of each model is colored as green and orange for TT model and NT model, respectively.\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n37\n741 S4 Fig. Example of modeling result that the conformational state is not matched with template. \n742 A. Overall structure comparison between the modeled structure (colored in beige) and the \n743 template (PDB ID: 5XY1, colored in sky blue). B. Focused view for the binding site. C. \n744 Ramachandran plot for modeled structure and the template. The modeled structure were \n745 annotated as DFGin-BLAminus (red x), while the template was DFGin-ABAminus (green x).\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint \n\n.CC-BY 4.0 International licenseavailable under a \nwas not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprint (whichthis version posted April 6, 2024. ; https://doi.org/10.1101/2024.04.04.588044doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}