Afin d’évaluer l’effet de la protonation de l’histidine et des états rotamériques sur la performance prédictive des récepteurs, nous avons effectué un criblage virtuel (VS) pour l’enzyme Mtb RmlC sur la base des résultats d’une étude précédente de criblage à haut débit (HTS). Ci-dessous, nous examinerons d’abord les interactions typiques du ligand co-cristallin TRH pour sonder la dépendance de la pose du ligand à la protonation de l’histidine. Nous contextualisons ensuite l’analyse de la performance d’enrichissement et du pouvoir prédictif de divers modèles de récepteurs, en discutant les interactions avec le récepteur pour montrer l’effet des différents états de protonation de l’histidine sur le VS. Enfin, nous comparons les valeurs de pKa prédites calculées par plusieurs paquets de calcul de pKa communs aux états de protonation du récepteur avec le meilleur pouvoir prédictif.
Docking de TRH
Docking du ligand co-cristallin TRH en retour dans 36 modèles de récepteur a été effectué pour montrer la pose, ou l’orientation du ligand par rapport au récepteur, la dépendance à la protonation de l’histidine et aux états rotamériques. Les modèles de liaison hydrogène chimiquement intuitifs pour les coordonnées cristallines de His62 et His119, présentés dans la figure 2b, impliquent l’importance potentielle des liaisons hydrogène dans le docking de la TRH. Le docking de ce ligand a permis un examen préliminaire de la dépendance de la pose sur les réseaux de liaison hydrogène possibles avec le récepteur.
La variation des états de protonation de l’histidine a un effet clair sur la prédiction de la pose pour le ligand cocristallin déterminé. La RMSD de la pose de docking de la TRH auto-dockée dans les coordonnées cristallines pour différents états de protonation et de rotation de His62 et His119 a varié de 2,91 à 5,44 Å. L’état de protonation des deux histidines avec la meilleure RMSD moyenne est HIE, ce qui correspond aux états de protonation les plus probables des coordonnées cristallines de la TRH. De plus, dans tous les cas, l’algorithme d’arrimage prédit correctement la position du pyrophosphate du ligand, mais la grande déviation par rapport aux coordonnées cristallines provient principalement de l’inversion des fragments de thymidine et de rhamnose autour du pyrophosphate, ce qui entraîne des modèles de liaison hydrogène différents entre la TRH et les deux histidines. Cela indique l’importance des réseaux de liaison hydrogène avec His62 et His119 dans la prédiction de la pose du ligand co-cristallin TRH. Par conséquent, après avoir examiné la dépendance de la pose aux liaisons hydrogène fournies par deux histidines, nous avons élargi notre étude pour examiner systématiquement le classement des composés dans le VS et comment il est affecté par les états de protonation et de rotomérie des histidines.
Criblage virtuel
Le docking moléculaire a été effectué pour examiner la dépendance du classement des composés aux états de protonation et de rotomérie des histidines. L’ensemble de ligands comprenait dix actifs et 2 000 inactifs choisis au hasard dans un HTS. Nous notons que les scores de Tanimoto indiquent que la plupart de nos leurres ont une faible similarité avec les actifs. Un tel ensemble de leurres présente un plus petit défi pour l’algorithme d’amarrage et la performance prédictive de VS lui-même peut être affectée lorsque des leurres ayant une plus grande similarité avec les actifs sont utilisés. Cependant, cette étude visait à examiner non pas la performance prédictive de l’algorithme de docking en soi, mais comment les états de protonation de l’histidine affectent la performance relative dans VS.
Ligands actifs dockés et l’analogue du produit ont été examinés initialement pour caractériser les interactions importantes dans le site de liaison de RmlC. Dans tous les modèles de récepteurs, les interactions hydrophobes d’empilement pi-pi contribuent de manière significative au score de docking des composés actifs dans le site actif de RmlC. Le composé initial de l’étude HTS, SID7975595, est bien classé dans la plupart des modèles de récepteurs, entre le 8ème et le 51ème rang dans 26 des 36 récepteurs. Bien qu’il n’y ait qu’une similarité structurelle limitée entre le SID7975595 et le ligand co-cristallin TRH, le cycle tricyclique du SID7975595 remplace facilement la partie thymidine de la TRH, tandis que le cycle benzimidazolone remplace la partie rhamnose, fournissant ainsi la base structurelle de l’inhibition. Comme le montre la figure 3, l’interaction hydrophobe entre les actifs et le récepteur implique souvent Tyr132 et Tyr138 de la chaîne A et Phe26 de la chaîne B (notez qu’une partie de la chaîne B fait intrusion dans le site actif de la chaîne A). En interagissant avec les résidus essentiels du site de liaison et en empêchant les molécules d’eau d’accéder à Phe26 et Tyr132, les actifs fournissent des contacts hydrophobes abondants pour obtenir une affinité de liaison élevée. Comme discuté dans Sivendran et al, la substitution du groupe éthyle attaché à l’azote du cycle tricyclique de SID7975595 par un groupe allyle (par exemple, le composé actif 77074) augmente encore l’affinité de liaison en formant un joint hydrophobe encore plus étanche. En comparaison, la substitution de ce groupe par un groupe méthyle plus petit ou un atome d’hydrogène entraîne une affinité de liaison plus faible. En plus des contacts hydrophobes décrits ci-dessus, certains des actifs forment également des liaisons hydrogène avec Ser51, Arg59 et Arg170. Une figure décrivant les interactions des actifs dockés peut être trouvée dans la ressource en ligne 3.
Intéressant, les actifs ne réalisent généralement pas d’interactions polaires avec His62 et His119. Comme le montre la figure 3, l’oxygène du carbonyle et les deux azotes de la benzimidazolone de SID7975595 sont orientés à l’opposé de His62 et His119. La direction des hydrogènes aromatiques des actifs est souvent incapable de participer à des réseaux de liaison hydrogène avec les deux histidines. Néanmoins, les différents états de protonation et de rotation de ces histidines affectent les résultats de VS par leurs interactions avec les leurres.
Évaluation des différences de classement
Il n’est pas rare que seul le top 1 % des composés criblés puisse être testé expérimentalement dans une étude VS, en raison des ressources limitées. Par conséquent, le facteur d’enrichissement (EF)1%, qui reflète la performance d’enrichissement de la base de données dans le top 1 % (20 composés dockés) d’une bibliothèque, devient particulièrement pertinent pour évaluer le pouvoir prédictif du VS. Le EF1% varie de 0 à 80 pour 36 modèles de récepteurs (Tableau 1), ce qui indique que les résultats de la VS sont sensibles aux états de protonation et de rotaméricité de His62 et His119 de RmlC. Néanmoins, 28 des 36 récepteurs classent plus de huit actifs dans les 10 % supérieurs de la VS, comme le reflète le EF10% (Tableau 1), ce qui suggère que la plupart des récepteurs sont capables de distinguer les actifs et les leurres lorsqu’une plus grande partie (10 %) de la base de données est considérée. Les résultats EF suggèrent également que les modèles de récepteur avec HIP62 ou HIP119 ont tendance à avoir une mauvaise performance d’enrichissement, probablement en raison des réseaux de liaison hydrogène étendus avec les leurres, comme discuté plus tard.
L’aire sous la courbe caractéristique d’exploitation du récepteur (AUC) pour chaque modèle de récepteur a été évaluée pour rendre compte de la performance d’enrichissement des modèles sur différents états de protonation et de rotamère de His62 et His119. Comme le montrent la figure 4a et le tableau 1, les valeurs de l’AUC de tous les modèles de récepteurs vont de 0,868 à 0,996, ce qui indique une bonne performance prédictive globale (une AUC de 0,5 correspond à une absence de différenciation entre les actifs et les leurres). En général, le résultat de l’AUC est complémentaire de l’évaluation de l’EF pour la performance prédictive des récepteurs. En résumant le tableau 1, la figure 4c montre comment la gamme des performances du récepteur dépend des deux états de protonation et de rotation de l’histidine. Si l’on considère la plage de 25 à 75 % des AUC (Fig. 4c, indiquée par les lignes plus épaisses), les modèles His62 montrent une plus grande variation entre les états His119. Les modèles His119, en revanche, ont une performance plus cohérente quels que soient les états de protonation de His62, à l’exception de l’état HIP. Cela indique que les différents états de protonation de His62 ont une influence plus faible que ceux de His119 sur la performance des récepteurs dans notre criblage.
Une plus forte dépendance de l’enrichissement aux états de protonation de His119 est observée dans les modèles HIE62 et HIP62. Avec l’état HIE62, les modèles avec HIP119 retourné (modèle 6) et HIE119 retourné (modèle 2) donnent les meilleures performances du récepteur. Les modèles 3 et 5 avec HID119 et HIP119, respectivement, conduisent au pire enrichissement. En examinant pourquoi l’état HIE62 a la plus grande variation dans les AUC, on constate que His62 a soit un empilement pi-pi, soit aucune interaction avec les ligands, et ne fait que quelques liaisons hydrogène avec les leurres de haut rang. Par conséquent, la performance du récepteur dépend de l’interaction de His119 avec les decoys. C’est également ce que l’on constate en examinant la large gamme de performances des AUC des modèles HIP62. Les réseaux de liaison hydrogène avec les decoys seront discutés plus loin dans la section suivante.
Afin d’évaluer la signification statistique de la différence des valeurs AUC entre une paire de modèles de récepteurs, nous avons effectué un test p bilatéral au niveau de 95 % sur l’hypothèse nulle que la paire a des valeurs AUC statistiquement comparables, contre l’hypothèse alternative que leur différence dans les valeurs AUC et le pouvoir prédictif est statistiquement significative. Le résultat est présenté dans le tableau 1 des ressources en ligne, les valeurs p inférieures à 0,05 étant soulignées. En moyenne, les récepteurs ont plus de 16 valeurs p inférieures à 0,05, ce qui démontre la sensibilité de la VS sur la protonation de l’histidine et les états rotamériques. Comme on pouvait s’y attendre, les récepteurs présentant les différences les plus significatives correspondent aux modèles ayant les valeurs d’AUC les plus élevées (modèle 6) ou les plus faibles (modèles 3, 29 et 5). Le modèle 6 est statistiquement meilleur pour classer les actifs par rapport aux leurres que les 26 autres récepteurs de l’ensemble. Les modèles 3, 29 et 5 sont distinctement moins bons pour classer les actifs que 29, 25 et 31 autres récepteurs, respectivement.
Une analyse quantitative des interactions de liaison hydrogène a été réalisée pour les 1 % supérieurs (20 composés dockés) de chaque résultat VS pour tenir compte des interactions de liaison hydrogène abondantes avec les résidus du site de liaison souvent observées avec les leurres. Les résultats indiquent une corrélation inverse entre la contribution des liaisons hydrogène et la performance des récepteurs. La figure 4b montre le pourcentage moyen de liaison hydrogène de chaque modèle de récepteur pour les 1 % de composés les mieux ancrés. Le pourcentage de liaison hydrogène est défini comme la part du terme de liaison hydrogène de Glide XP dans le score total d’amarrage. La comparaison des figures 4a et b révèle une relation inverse entre le pourcentage de liaison hydrogène et la SSC avec un R2 de 0,42 (y = -56,18x + 67,95, la corrélation est tracée dans la ressource en ligne 4). La relation inverse est généralement observée avec les modèles comportant HIP119, HIP119 inversé ou HID62, où le pourcentage élevé de liaisons hydrogène entraîne un mauvais enrichissement. Par exemple, le modèle de récepteur 29 avec HIP62 et HIP119, où les deux histidines présentant des donneurs de liaison hydrogène faisant face au site actif, a l’une des pires AUC en raison du pourcentage élevé de liaisons hydrogène dans les meilleurs résultats.
Notamment, le potentiel de liaison hydrogène de His119 détermine souvent la performance du récepteur. Par exemple, le modèle avec HID62 et HIP119 était une aberration parmi les modèles HID62 dans la figure 4c, avec un enrichissement remarquablement faible par rapport à la bonne performance globale des cinq autres modèles HID62. Les modèles HID62 ont une AUC médiane élevée de 0,989, malgré les fréquentes liaisons hydrogène avec les leurres de HID62. Cela est dû au fait que les états His119 réalisent peu d’interactions de liaison hydrogène avec les leurres. Ce n’est qu’avec l’état HIP119 que le modèle HID62 établit des liaisons hydrogène avec un certain nombre de leurres, ce qui explique l’AUC relativement faible. Cette observation concorde avec la plus forte dépendance de la performance du récepteur aux états de protonation de His119, comme discuté ci-dessus. La ressource en ligne 4 décrit la distribution de l’AUC et le pourcentage de liaison hydrogène ainsi que la direction du donneur ou de l’accepteur de liaison hydrogène à partir de deux histidines faisant face au récepteur.
Les analyses ci-dessus mettent en évidence l’effet d’entrave de la liaison hydrogène aux leurres sur le pouvoir prédictif de la VS, en raison des diverses coordonnées de deux histidines avec des états de protonation et de rotamérisation différents. La dispersion de la corrélation observée avec le R2 de 0,42 est probablement attribuée à plusieurs causes, notamment la nature chimique de l’ensemble de données des leurres, ainsi que les légères différences de géométrie de chaque récepteur lors de la minimisation dans la préparation initiale de la protéine. En montrant clairement la sensibilité des résultats du criblage virtuel sur les différents états de protonation et de rotation des histidines dans le site actif, nous soulignons qu’il faut être prudent lors de la préparation des coordonnées atomiques d’un récepteur pour le VS, en particulier si l’on considère les propriétés générales des ligands à cribler. Cela inclut la prise en compte de la liaison hydrogène avec le ligand co-cristallin et son effet sur la préparation de la protéine, ainsi qu’une analyse complète des réseaux de liaison hydrogène proximaux. Ceci est généralement réalisé en examinant les résultats de progiciels de prédiction de pKa largement utilisés, et jusqu’à présent, nous avons comparé les résultats de différents progiciels par rapport à nos résultats VS et nous les discutons plus loin.
Docking des leurres
Divers facteurs conduisent à des différences de classement entre les récepteurs, en particulier en ce qui concerne les leurres. En général, les leurres qui se sont classés plus haut que les actifs avaient un poids moléculaire élevé et avaient plus de potentiel pour avoir des liaisons hydrogène avec le récepteur. Dans cette section, nous analysons plus en détail les schémas d’interaction fréquents observés entre les leurres et le récepteur, en nous concentrant sur les modèles de récepteur avec un faible enrichissement.
Les leurres ont tendance à avoir un poids moléculaire plus important et plus de structures cycliques que les actifs (tableau 2). Il en résulte un classement plus élevé des leurres, en raison des interactions hydrophobes en l’absence de liaisons hydrogène avec le récepteur. La figure 5a montre les interactions hydrophobes obtenues par le grand composé inactif 16952387 dans le modèle de récepteur 19. Ce composé est souvent classé parmi les cinq premiers dans de nombreux criblages virtuels pour ses importantes interactions d’empilement pi-pi avec Phe26, Tyr132 et Tyr138. Cette tendance est fréquemment observée dans le criblage virtuel où les plus grandes molécules se classent mieux en raison des interactions étendues avec le récepteur .
La performance d’enrichissement est particulièrement faible pour les récepteurs fournissant des réseaux de liaisons hydrogène abondants aux leurres. Les interactions par His62 et His119 n’ont pas été largement observées pour les actifs, et donc les composés avec des contributions enthalpiques plus importantes se classent à tort plus favorablement. Un exemple illustré à la figure 5b montre l’interaction du composé inactif 17388064 dans le modèle de récepteur 3 (AUC 0,868), classé en premier. Dans ce récepteur, qui est le plus mauvais pour classer les composés en fonction de l’AUC, le composé 17388064 forme deux liaisons hydrogène avec deux histidines, l’une entre son hydrogène hydroxyle et l’azote δ de HIE62 et l’autre entre son oxygène hydroxyle et l’hydrogène sur l’azote δ de HID119. Ce composé possède cinq donneurs de liaisons hydrogène et neuf accepteurs, un nombre important comparé aux moyennes respectives de celles des leurres et des actifs (Tableau 2). Par conséquent, avec une contribution élevée des liaisons hydrogène au score total de 34,7 ± 6,62 %, on observe fréquemment que ce composé leurre forme au moins une liaison hydrogène avec l’une ou l’autre des deux histidines, ce qui lui permet d’obtenir des rangs élevés dans de multiples passages VS.
Deux autres modèles de récepteurs, le modèle 29 avec HIP62 et HIP119 et le modèle 5 avec HIE62 et HIP119, présentent des schémas d’interaction avec les leurres similaires à ceux du modèle 3. Ces trois modèles présentent les valeurs d’AUC les plus faibles, avec une moyenne de 0,870 entre eux. Comme nous l’avons vu plus haut, leurs valeurs AUC diffèrent considérablement de celles des autres récepteurs, ce qui reflète la relation subtile entre les liaisons hydrogène réalisées par His62 et His119 et le faible enrichissement. Une figure supplémentaire décrivant les réseaux de liaisons hydrogène entre les leurres et les récepteurs est fournie dans la ressource en ligne 5.
prédiction du pKa pour His62 et His119
Nos résultats démontrent clairement la sensibilité du criblage virtuel sur la protonation de l’histidine et les états rotamériques. Dans de nombreuses études biophysiques computationnelles, les états de protonation des résidus titrables sont déterminés à l’aide de divers programmes de prédiction du pKa. Afin d’évaluer la performance de ces programmes pour identifier le modèle de récepteur avec le meilleur pouvoir prédictif dans le docking, nous avons comparé les résultats de prédiction de pKa de His62 et His119 de PROPKA, Maestro, H++, et MCCE, comme indiqué dans le tableau 3 avec les valeurs de pKa calculées.
Premièrement, PROPKA 3.1 prédit que His62 et His119 sont neutres indépendamment de la présence de TRH pendant la préparation. Le programme, cependant, ne peut pas attribuer d’états rotamériques aux histidines. Par conséquent, un état de HID, HID inversé, HIE ou HIE inversé doit être déterminé manuellement. Comme PROPKA, le programme H++, qui utilise une électrostatique continue à structure unique, trouve également que les deux histidines sont neutres, bien que les valeurs de pKa prédites soient différentes de celles de PROPKA. Le programme MCCE, qui est basé sur l’électrostatique continue multi-conformation, prédit que l’His62 est neutre tandis que l’His119 est protonée.
Puis, nous avons utilisé l’assistant de préparation des protéines dans Maestro pour calculer les pKa de l’His62 et de l’His119 avec et sans TRH. Notez que Maestro est capable de faire varier les états rotamériques, alors que PROPKA ne le peut pas. Une mise à jour récente permet à Maestro d’utiliser PROPKA dans sa prédiction de pKa au lieu d’Epik. Avec Epik, Maestro prédit à la fois His62 et His119 dans des états doublement protonés, indépendamment de la présence de TRH. Il est intéressant de noter que le modèle de récepteur qui correspond à cet état multi-histidine a le pire pouvoir de prédiction avec une AUC de 0,869. Lorsque PROPKA est utilisé, HIP62 et HIE119 sont prédits pour le complexe protéine-TRH et HIE62 et HIE119 pour la protéine apo. Ces deux prédictions par PROPKA dans Maestro correspondent aux modèles de performance d’enrichissement modérée, avec des AUC de 0,971 pour le modèle 25 (HIP62 et HIE119) et de 0,942 pour le modèle 1 (HIE62 et HIE119), respectivement.
Compte tenu du fait que les prédictions ci-dessus effectuées par différents logiciels varient considérablement les unes des autres, il convient de faire preuve de prudence lors de l’utilisation de ces résultats comme ligne directrice pour préparer une protéine pour le criblage virtuel. Sans une connaissance intime du véritable état de protonation du récepteur ainsi que des ligands criblés, il est difficile d’aborder ce problème. Par conséquent, nous suggérons qu’une analyse à petite échelle, comme celle réalisée dans cette étude, et une comparaison avec des données expérimentales, si elles sont disponibles, pourraient fournir une description plus précise de la protonation et des états rotamériques des résidus titrables dans les récepteurs protéiques pour les futurs criblages à plus grande échelle. Alternativement, un modèle qui inclut l’incorporation explicite d’états alternatifs de protonation et de rotation de la chaîne latérale pendant le docking, potentiellement avec des informations stockées dans la grille comme cela existe pour les hydroxyles et les thiols rotatifs dans Glide, peut valoir la peine d’être poursuivi. L’examen des résultats en ce qui concerne les états de protonation et les classements basés sur les interactions avec les histidines devrait être examiné attentivement avant de passer aux tests expérimentaux.
En outre, la flexibilité du récepteur affectera probablement les états de protonation des résidus ionisables. Bien que cela n’ait pas été explicitement étudié ici, à part la minimisation de chaque récepteur après l’attribution des états de protonation, la flexibilité des protéines est clairement importante pour la conception et le développement de médicaments . Considérer conjointement l’espace conformationnel et l’espace de protonation devient rapidement intraitable avec des méthodes physiques telles que celles décrites ici, mais des méthodes d’échantillonnage améliorées sont prometteuses pour résoudre ces difficultés. Cela inclut les simulations de dynamique moléculaire à pH constant, pour lesquelles le pH est une variable thermodynamique externe, utilisées pour la prédiction aveugle des valeurs de pKa des résidus titrables . L’application efficace des résultats de ces simulations à la conception moléculaire est un domaine d’intérêt permanent. Les ensembles d’équilibre issus de ces simulations peuvent être utilisés conjointement avec le docking comme une application du schéma complexe relaxé, où le criblage virtuel est effectué avec un ensemble de structures différemment protonées, pour améliorer les résultats de l’enrichissement . La prise en compte de la flexibilité des récepteurs dans la préparation de la cible conduira à un échantillonnage plus large de l’espace conformationnel et de protonation, améliorant ainsi les performances du VS.