Con el fin de evaluar el efecto de la protonación de la histidina y de los estados rotaméricos en el rendimiento predictivo de los receptores, realizamos un cribado virtual (VS) para la enzima RmlC de Mtb basándonos en los resultados de un estudio previo de cribado de alto rendimiento (HTS). A continuación, examinaremos en primer lugar las interacciones típicas del ligando co-cristalino TRH para sondear la dependencia de la pose del ligando en la protonación de la histidina. Además, contextualizamos el análisis del rendimiento del enriquecimiento y el poder predictivo de varios modelos de receptor, analizando las interacciones con el receptor para mostrar el efecto de los diferentes estados de protonación de la histidina en la VS. Por último, comparamos los valores de pKa predichos calculados por varios paquetes de cálculo de pKa comunes con los estados de protonación del receptor con el mejor poder de predicción.
Docking de TRH
Se llevó a cabo el docking del ligando de co-cristal TRH en 36 modelos de receptor para mostrar la dependencia de la pose, o la orientación del ligando en relación con el receptor, de la protonación de la histidina y los estados rotaméricos. Los patrones de enlace de hidrógeno químicamente intuitivos para las coordenadas cristalinas de His62 e His119, mostrados en la Fig. 2b, implican la importancia potencial de los enlaces de hidrógeno en el acoplamiento de TRH. El docking de este ligando permitió examinar preliminarmente la dependencia de la pose de las posibles redes de enlaces de hidrógeno con el receptor.
La variación de los estados de protonación de la histidina tiene un claro efecto en la predicción de la pose para el ligando co-cristal determinado. La RMSD de la pose de acoplamiento de la TRH autoacoplada en las coordenadas del cristal para diferentes estados de protonación y rotaméricos de His62 e His119 varió de 2,91 a 5,44 Å. El estado de protonación de ambas histidinas con la mejor RMSD media es HIE, que coincide con los estados de protonación más probables de las coordenadas del cristal de la TRH. Además, en todos los casos el algoritmo de acoplamiento predice correctamente la posición del pirofosfato del ligando, pero la gran desviación respecto a las coordenadas cristalinas se debe principalmente a la inversión de las moléculas de timidina y ramnosa alrededor del pirofosfato, lo que da lugar a diferentes patrones de enlace de hidrógeno entre la TRH y las dos histidinas. Esto indica la importancia de las redes de enlace de hidrógeno con His62 e His119 en la predicción de la pose del ligando co-cristal TRH. Por lo tanto, después de examinar la dependencia de la pose de los enlaces de hidrógeno proporcionados por dos histidinas, ampliamos nuestro estudio para examinar sistemáticamente la clasificación de los compuestos en VS y cómo se ve afectada por la protonación y los estados rotoméricos de las histidinas.
Se realizó un cribado virtual
El acoplamiento molecular para examinar la dependencia de la clasificación de los compuestos en la protonación y los estados rotoméricos de las histidinas. El conjunto de ligandos incluía diez activos y 2.000 inactivos seleccionados al azar de un HTS. Observamos que las puntuaciones de Tanimoto indican que la mayoría de nuestros señuelos tienen una baja similitud con los activos. Este conjunto de señuelos presenta un reto menor para el algoritmo de acoplamiento y el rendimiento predictivo del propio VS puede verse afectado cuando se utilizan señuelos con mayor similitud a los activos. Sin embargo, este estudio pretendía examinar no el rendimiento predictivo del algoritmo de acoplamiento per se, sino cómo los estados de protonación de la histidina afectan al rendimiento relativo en el VS.
Ligandos activos acoplados y el producto análogo fueron examinados inicialmente para caracterizar las interacciones importantes en el sitio de unión de RmlC. En todos los modelos de receptores, las interacciones hidrofóbicas de apilamiento pi-pi contribuyen significativamente a la puntuación de acoplamiento de los compuestos activos dentro del sitio activo de RmlC. El compuesto inicial de HTS, SID7975595, ocupa una posición alta en la mayoría de los modelos de receptores, entre el 8º y el 51º puesto en 26 de los 36 receptores. Aunque la similitud estructural entre SID7975595 y el ligando co-cristalino TRH es limitada, el anillo tricíclico de SID7975595 sustituye fácilmente a la fracción de timidina de TRH, mientras que el anillo de benzimidazolona sustituye a la fracción de ramnosa, proporcionando la base estructural de la inhibición. Como se muestra en la Fig. 3, la interacción hidrofóbica entre los activos y el receptor suele implicar a las Tyr132 y Tyr138 de la cadena A y a la Phe26 de la cadena B (nótese que una parte de la cadena B se inmiscuye en el sitio activo de la cadena A). Al interactuar con los residuos esenciales del sitio de unión y evitar que las moléculas de agua accedan a Phe26 y Tyr132, los activos proporcionan abundantes contactos hidrofóbicos para lograr la alta afinidad de unión. Como se discute en Sivendran et al, la sustitución del grupo etilo unido al nitrógeno del anillo tricíclico de SID7975595 por un grupo alilo (por ejemplo, el compuesto activo 77074) mejora aún más la afinidad de unión al formar un sello hidrofóbico aún más estrecho. En comparación, la sustitución de este grupo por un grupo metilo más pequeño o un átomo de hidrógeno da lugar a una menor afinidad de unión . Además de los contactos hidrofóbicos descritos anteriormente, algunos de los activos también forman enlaces de hidrógeno con Ser51, Arg59 y Arg170. Una figura que describe las interacciones de los activos acoplados puede encontrarse en el Recurso Online 3.
Interacción prevista del compuesto de éxito inicial SID7975595 con HID62 y HIP119 volteados en el modelo de receptor 23. En general, los activos no tienen interacciones fuertes con His62 o His119, pero la variación de los estados de protonación de la histidina tiene un efecto profundo en los resultados clasificados. Se observan interacciones favorables con otros residuos del sitio de unión, como Tyr132 y Tyr138, como se muestra aquí
Interesantemente, los activos generalmente no logran interacciones polares con His62 e His119. Como se muestra en la Fig. 3, el oxígeno del carbonilo y los dos nitrógenos de la benzimidazolona de SID7975595 se orientan en dirección contraria a His62 e His119. La dirección de los hidrógenos aromáticos de los activos suele ser incapaz de participar en redes de enlace de hidrógeno con las dos histidinas. No obstante, los diferentes estados de protonación y rotaméricos de estas histidinas sí afectan a los resultados de la VS a través de sus interacciones con los señuelos.
Evaluación de las diferencias en la clasificación
No es infrecuente que sólo el 1 % superior de los compuestos cribados pueda probarse experimentalmente en un estudio de la VS, debido a la limitación de recursos. Por lo tanto, la métrica del factor de enriquecimiento (EF)1%, que refleja el rendimiento del enriquecimiento de la base de datos en el 1 % superior (20 compuestos acoplados) de una biblioteca, adquiere especial relevancia para evaluar el poder predictivo de la VS. El EF1% oscila entre 0 y 80 para 36 modelos de receptores (Tabla 1), lo que indica que los resultados de la VS son sensibles a los estados de protonación y rotaméricos de His62 e His119 de RmlC. Sin embargo, 28 de los 36 receptores clasifican más de ocho activos dentro del 10 % superior en la VS, como refleja el EF10% (Tabla 1), lo que sugiere que la mayoría de los receptores son capaces de distinguir los activos y los señuelos cuando se considera una porción mayor (10 %) de la base de datos. Los resultados de EF también sugieren que los modelos de receptores con HIP62 o HIP119 tienden a tener un rendimiento de enriquecimiento pobre, probablemente debido a las extensas redes de enlace de hidrógeno con los señuelos, como se discute más adelante.
Se evaluó el área bajo la curva característica de funcionamiento del receptor (AUC) para cada modelo de receptor para informar del rendimiento de enriquecimiento de los modelos sobre diferentes estados de protonación y rotaméricos de His62 e His119. Como se muestra en la Fig. 4a y en la Tabla 1, los valores de AUC de todos los modelos de receptores oscilan entre 0,868 y 0,996, lo que indica un buen rendimiento predictivo general (un AUC de 0,5 corresponde a la no diferenciación entre los activos y los señuelos). En general, el resultado del AUC es complementario a la evaluación del EF para el rendimiento predictivo del receptor. Resumiendo la Tabla 1, la Fig. 4c muestra cómo el rango del rendimiento del receptor depende de los dos estados de protonación de la histidina y de los rotaméricos. Considerando el rango del 25-75 % de las AUCs (Fig. 4c, indicado por las líneas más gruesas), los modelos His62 muestran una mayor variación a través de los estados His119. Los modelos de His119, por otro lado, tienen un rendimiento más consistente independientemente de los estados de protonación de His62, con la excepción del estado HIP. Esto indica que los diferentes estados de protonación de His62 tienen una menor influencia que los de His119 en el rendimiento del receptor en nuestro cribado.
a Valores AUC de 36 modelos de receptor. Los estados de protonación y rotaméricos están marcados para cada histidina. Los estados invertidos se marcan con la letra F. Un color más oscuro indica un mayor AUC y un mejor rendimiento predictivo del modelo de receptor correspondiente. b Porcentaje medio de enlaces de hidrógeno del 1 % de los mejores compuestos en 36 ejecuciones de VS. Los estados de protonación y rotamérico están marcados para cada histidina. El color más claro indica un mayor porcentaje de enlaces de hidrógeno, con la unidad de % para la barra de color. El R2 para la correlación entre los AUC y el porcentaje medio de enlaces de hidrógeno para cada ejecución de VS es de 0,42 (véase el Recurso en línea 4 para el gráfico de dispersión). c Dependencia del rendimiento del receptor en His62 (arriba) e His119 (abajo). La mediana de los valores AUC de cada estado de protonación se muestra con una línea horizontal grande. Las pequeñas marcas en cada modelo de histidina marcan seis estados de protonación diferentes de la otra histidina. Las líneas verticales más gruesas representan el rango 25-75 % de las AUC. Los mejores modelos de receptor se muestran explícitamente con los estados de protonación de los modelos
Se observa una mayor dependencia del enriquecimiento de los estados de protonación de His119 en los modelos HIE62 e HIP62. Con el estado HIE62, los modelos con HIP119 volteado (modelo 6) y HIE119 volteado (modelo 2) producen el mayor rendimiento del receptor. Los modelos 3 y 5 con HID119 y HIP119, respectivamente, conducen al peor enriquecimiento. Al examinar por qué el estado HIE62 tiene la mayor variación en las AUC, se encuentra que His62 tiene apilamiento pi-pi o no tiene interacciones con los ligandos, y hace sólo unos pocos enlaces de hidrógeno con los señuelos de alto rango. Por tanto, el rendimiento del receptor depende de la interacción de His119 con los señuelos. Esto también se observa al examinar el amplio rango de rendimiento de las AUC de los modelos HIP62. Las redes de enlace de hidrógeno con los señuelos se discutirán más adelante en la siguiente sección.
Para evaluar la significación estadística de la diferencia de los valores AUC entre un par de modelos de receptor, realizamos una prueba p de dos caras al nivel del 95 % sobre la hipótesis nula de que el par tiene valores AUC estadísticamente comparables, contra la hipótesis alternativa de que su diferencia en los valores AUC y el poder predictivo es estadísticamente significativa. El resultado se muestra en la tabla de recursos en línea 1, donde se destacan los valores p inferiores a 0,05. En promedio, los receptores tienen más de 16 valores de p inferiores a 0,05, lo que demuestra la sensibilidad de la VS sobre la protonación de la histidina y los estados rotaméricos. Como era de esperar, los receptores con las diferencias más significativas corresponden a los modelos con los valores AUC más altos (modelo 6) o más bajos (modelos 3, 29 y 5). El modelo 6 es estadísticamente mejor en la clasificación de los activos sobre los señuelos que otros 26 receptores del conjunto. Los modelos 3, 29 y 5 son notablemente peores en la clasificación de los activos que otros 29, 25 y 31 receptores, respectivamente.
Se llevó a cabo un análisis cuantitativo de las interacciones de enlace de hidrógeno para el 1 % superior (20 compuestos acoplados) de cada resultado de VS para dar cuenta de las abundantes interacciones de enlace de hidrógeno con los residuos del sitio de unión que se observan a menudo con los señuelos. Los resultados indican una correlación inversa entre la contribución del enlace de hidrógeno y el rendimiento del receptor. La figura 4b muestra el porcentaje medio de enlaces de hidrógeno de cada modelo de receptor para el 1 % superior de compuestos acoplados. El porcentaje de enlace de hidrógeno se define como la parte del término de enlace de hidrógeno de Glide XP en la puntuación total de acoplamiento. La comparación de la Fig. 4a, b revela la relación inversa entre el porcentaje de enlace de hidrógeno y el AUC con un R2 de 0,42 (y = -56,18x + 67,95, la correlación se traza en el Recurso Online 4). La relación inversa se observa comúnmente con los modelos con HIP119, HIP119 volteado, o HID62, donde el alto porcentaje de enlace de hidrógeno resultó en un pobre enriquecimiento. Por ejemplo, el modelo de receptor 29 con HIP62 y HIP119, en el que ambas histidinas presentan donantes de enlaces de hidrógeno orientados hacia el sitio activo, tiene uno de los peores AUC debido al alto porcentaje de enlaces de hidrógeno en los primeros aciertos.
En particular, el potencial de enlace de hidrógeno de His119 a menudo determina el rendimiento del receptor. Por ejemplo, el modelo con HID62 y HIP119 fue un valor atípico entre los modelos HID62 en la Fig. 4c, con un enriquecimiento notablemente bajo en comparación con el buen rendimiento general de los otros cinco modelos HID62. Los modelos HID62 tienen una mediana de AUC alta de 0,989, a pesar de los frecuentes enlaces de hidrógeno con los señuelos de HID62. Esto se debe a que los estados His119 logran pocas interacciones de enlace de hidrógeno con los señuelos. Sólo con el estado HIP119 el modelo HID62 realiza enlaces de hidrógeno con varios señuelos, lo que da lugar a un AUC relativamente bajo. Esta observación concuerda con la mayor dependencia del rendimiento del receptor en los estados de protonación de His119, como se ha comentado anteriormente. El Recurso en línea 4 describe la distribución del AUC y el porcentaje de enlaces de hidrógeno junto con la dirección del donante o aceptor de enlaces de hidrógeno de dos histidinas frente al receptor.
Los análisis anteriores ponen de manifiesto el efecto obstaculizador de los enlaces de hidrógeno a los señuelos en el poder predictivo de VS, debido a las diversas coordenadas de dos histidinas con diferentes estados de protonación y rotaméricos. La dispersión de la correlación observada con el R2 de 0,42 se atribuye probablemente a varias causas, incluyendo la naturaleza química del conjunto de datos de los señuelos, así como las ligeras diferencias en la geometría de cada receptor al minimizarlo en la preparación inicial de la proteína. Al mostrar claramente la sensibilidad de los resultados del cribado virtual sobre los diferentes estados de protonación y rotaméricos de las histidinas en el sitio activo, enfatizamos que se debe tener cuidado al preparar las coordenadas atómicas de un receptor para la VS, particularmente considerando las propiedades generales de los ligandos que se están cribando. Esto incluye tener en cuenta el enlace de hidrógeno con el ligando del cocristal y su efecto en la preparación de la proteína, así como un análisis exhaustivo de las redes de enlace de hidrógeno proximales. Esto suele lograrse examinando los resultados de paquetes de software de predicción de pKa ampliamente utilizados, y hasta este punto, hemos comparado los resultados de diferentes paquetes en relación con nuestros resultados de VS y los discutimos más adelante.
Docking de los señuelos
Varios factores conducen a diferencias en la clasificación entre los receptores, particularmente con respecto a los señuelos. Por lo general, los señuelos que se clasificaron más alto que los activos eran de alto peso molecular y tenían más potencial para tener enlaces de hidrógeno con el receptor. En esta sección, analizamos más a fondo los patrones de interacción frecuentes observados entre los señuelos y el receptor, centrándonos en los modelos de receptores con escaso enriquecimiento.
Los señuelos tienden a tener un mayor peso molecular y más estructuras de anillo que los activos (Tabla 2). Esto hace que los señuelos tengan una mayor clasificación, debido a las interacciones hidrofóbicas en ausencia de enlaces de hidrógeno con el receptor. La figura 5a muestra las interacciones hidrofóbicas logradas a través del gran compuesto inactivo 16952387 en el modelo de receptor 19. Este compuesto suele estar entre los cinco primeros en muchas ejecuciones de cribado virtual por sus importantes interacciones de apilamiento pi-pi con Phe26, Tyr132 y Tyr138. Esta tendencia se observa con frecuencia en el cribado virtual, donde las moléculas más grandes se clasifican mejor como resultado de las extensas interacciones con el receptor .
a Interacción del compuesto inactivo 16952387 con HID62 y HIE119 volteados en el modelo 19 del receptor. El compuesto no tiene interacción con ninguna de las dos histidinas. Las interacciones de apilamiento Pi-pi con Phe26 de la cadena B, Tyr132 y Tyr138 contribuyen a su alto rango, junto con los enlaces de hidrógeno con Arg23, Arg59, Arg170 y Ser51 (no mostrado). b Interacción del compuesto inactivo 17388064 con HIE62 y HID119 en el modelo de receptor 3. Ambas histidinas proporcionan enlaces de hidrógeno al compuesto
El rendimiento de enriquecimiento es particularmente bajo para los receptores que proporcionan abundantes redes de enlaces de hidrógeno a los señuelos. Las interacciones a través de His62 e His119 no se observaron ampliamente para los activos, y por lo tanto los compuestos con mayores contribuciones entálpicas se clasifican erróneamente más favorablemente. Un ejemplo mostrado en la Fig. 5b representa la interacción del compuesto inactivo 17388064 en el modelo de receptor 3 (AUC 0,868), clasificado como primero. En este receptor, que es el peor en la clasificación de compuestos basada en el AUC, el compuesto 17388064 forma dos enlaces de hidrógeno con dos histidinas, uno entre su hidrógeno de hidroxilo y el δ-nitrógeno de HIE62 y el otro entre su oxígeno de hidroxilo y el hidrógeno del δ-nitrógeno de HID119. Este compuesto tiene cinco donantes de enlaces de hidrógeno y nueve aceptores, un gran número comparado con las medias respectivas de los de los señuelos y los activos (Tabla 2). Por lo tanto, con una alta contribución de enlace de hidrógeno a la puntuación total de 34,7 ± 6,62 %, este compuesto señuelo se observa con frecuencia para formar al menos un enlace de hidrógeno con cualquiera de las dos histidinas, logrando así altos rangos en múltiples ejecuciones de VS.
Otros dos modelos de receptores, el modelo 29 con HIP62 y HIP119 y el modelo 5 con HIE62 y HIP119, muestran patrones de interacción con los señuelos similares al modelo 3. Estos tres modelos tienen los valores AUC más bajos, con una media de 0,870 entre ellos. Como se ha comentado anteriormente, sus valores de AUC difieren significativamente de los de otros receptores, lo que refleja la sutil relación entre los enlaces de hidrógeno logrados a través de His62 e His119 y el escaso enriquecimiento. Una figura adicional que describe las redes de enlaces de hidrógeno entre los señuelos y los receptores se proporciona en el Recurso Online 5.
Predicción de pKa para His62 e His119
Nuestros resultados demuestran claramente la sensibilidad del cribado virtual sobre la protonación de la histidina y los estados rotaméricos. En muchos estudios biofísicos computacionales, los estados de protonación de los residuos titulables se determinan utilizando varios programas de predicción de pKa. Para evaluar el rendimiento de estos programas para identificar el modelo de receptor con el mejor poder de predicción en docking, comparamos los resultados de predicción de pKa de His62 e His119 de PROPKA, Maestro, H++ y MCCE, como se muestra en la Tabla 3, con los valores de pKa calculados.
En primer lugar, PROPKA 3.1 predice que tanto His62 como His119 son neutros independientemente de la presencia de TRH durante la preparación. El programa, sin embargo, no puede asignar estados rotaméricos de las histidinas. Por lo tanto, se debe determinar manualmente un estado de HID, HID volteado, HIE o HIE volteado. Al igual que PROPKA, el programa H++, que utiliza una electrostática continua de estructura única, también encuentra que ambas histidinas son neutras, aunque los valores de pKa predichos son diferentes a los de PROPKA. El programa MCCE, que se basa en la electrostática continua multiconformación, predice que la His62 es neutra mientras que la His119 está protonada.
A continuación, utilizamos el Asistente de Preparación de Proteínas en Maestro para calcular el pKa de la His62 y la His119 con y sin TRH. Nótese que Maestro es capaz de variar los estados rotaméricos, mientras que PROPKA no puede. Una actualización reciente permite a Maestro emplear PROPKA en su predicción de pKa en lugar de Epik. Con Epik, Maestro predice tanto His62 como His119 en estados doblemente protonados, independientemente de la presencia de TRH. Curiosamente, el modelo de receptor que corresponde a este estado de multihistidina tiene el peor poder predictivo con un AUC de 0,869. Cuando se utiliza PROPKA, se predicen HIP62 y HIE119 para el complejo proteína-TRH y HIE62 y HIE119 para la proteína apo. Estas dos predicciones de PROPKA en Maestro corresponden a los modelos de rendimiento de enriquecimiento moderado, con AUCs de 0,971 para el modelo 25 (HIP62 y HIE119) y 0,942 para el modelo 1 (HIE62 y HIE119), respectivamente.
Dado que las predicciones anteriores realizadas por diferentes programas informáticos varían significativamente entre sí, se debe tener precaución al utilizar estos resultados como guía para preparar una proteína para el cribado virtual. Sin un conocimiento íntimo del verdadero estado de protonación del receptor, así como de los ligandos que se están cribando, es difícil abordar este problema. Por lo tanto, sugerimos que un análisis a pequeña escala, como el realizado en este estudio, y la comparación con los datos experimentales, si están disponibles, podrían proporcionar una descripción más precisa de los estados de protonación y rotaméricos de los residuos titulables en los receptores de proteínas para futuros cribados a mayor escala. Alternativamente, puede valer la pena un modelo que incluya la incorporación explícita de estados alternativos de protonación y rotamería de las cadenas laterales durante el acoplamiento, potencialmente con información almacenada en la rejilla como existe para los hidroxilos y tioles rotatorios en Glide. El examen de los resultados con respecto a los estados de protonación y las clasificaciones basadas en las interacciones con las histidinas debería examinarse cuidadosamente antes de proceder a las pruebas experimentales.
Además, la flexibilidad del receptor probablemente afectará a los estados de protonación de los residuos ionizables. Aunque esto no se estudió explícitamente aquí, aparte de la minimización de cada receptor después de asignar los estados de protonación, la flexibilidad de la proteína es claramente importante para el diseño y desarrollo de fármacos . Considerar el espacio conformacional y de protonación en conjunto se vuelve rápidamente intratable con métodos físicos como los descritos aquí, pero los métodos de muestreo mejorados se muestran prometedores para abordar tales dificultades . Esto incluye simulaciones de dinámica molecular a pH constante, para las que el pH es una variable termodinámica externa, utilizada para la predicción ciega de los valores de pKa de los residuos titulables . La aplicación efectiva de los resultados de estas simulaciones al diseño molecular es un área de interés actual. Los conjuntos de equilibrio de estas simulaciones pueden utilizarse junto con el docking como una aplicación del esquema complejo relajado, en el que el cribado virtual se realiza con un conjunto de estructuras protonadas de forma diferente, para mejorar los resultados de enriquecimiento . Si se tiene en cuenta la flexibilidad del receptor en la preparación de la diana, se obtendrá un muestreo más amplio del espacio conformacional y de protonación, mejorando así el rendimiento del VS.