Nuevos modelos de farmacóforos basados ​​en complejos covalentes y no covalentes de la proteasa principal (Mpro) del SARS-CoV-2 dilucidados mediante simulaciones MD de microsegundos

Simulaciones de dinámica molecular

Realizamos simulaciones MD de 500 ns para 15 complejos de proteína-ligando con posición de ligando experimental disponible (Tabla S1) y 38 complejos de proteína-ligando (Tabla S2, S3) con posición de ligando prevista a partir del acoplamiento molecular flexible. La desviación cuadrática media (RMSD) y la fluctuación cuadrática media por residuo (RMSF) sobre los átomos de carbono alfa para la etapa de producción MD de todos los sistemas se informan en Datos complementarios 1. Después de la simulación MD extensión de los complejos mejor clasificados de los cálculos de MM/PBSA, se obtuvieron los gráficos RMSD y RMSF (Figura S8, S9, S10, S11). Para algunos complejos, los gráficos RMSD mostraron posibles cambios conformacionales debido a fluctuaciones abruptas. Para comprender la flexibilidad de la proteína en detalle, la Fig. 3 presenta la MPro estructura cristalina del complejo 6LU7 en cintas y representación superficial coloreada con valores RMSF.

figura 3
figura 3

(A) Valores RMSF de simulación MD para el complejo 6LU7 aplicado en SARS-CoV-2 MPro Representación de cintas de dímero. (B) Valores RMSF de simulación MD para el complejo 6LU7 aplicado en SARS-CoV-2 MPro superficie del sitio activo.

Identificamos tres zonas importantes con valores más altos de RMSF. La primera zona se ubica alrededor de los residuos presentes en el rango 45-55, la segunda corresponde a los residuos ubicados en el rango 189-193 y la tercera zona corresponde al rango 200-300. Además, también se detectaron grandes movimientos en los residuos terminales. Las dos primeras zonas corresponden a los subsitios S2 (45-55) y S4 (189-193), por lo que la flexibilidad de estas áreas puede ser importante en la interacción proteína-ligando. Para destacar, el subsitio S2 presentó una pérdida gradual de estructura secundaria de hélice alfa a bucle para algunos complejos, como se muestra en la Figura S24 para el complejo 7JU7. Además, en algunos cristales de proteína como 6Y2F, observamos la falta de residuos en el subsitio S2, lo que podría estar relacionado con la flexibilidad de esta zona. Con respecto a los inhibidores covalentes y no covalentes, no se apreciaron diferencias notables entre estas categorías en los gráficos de RMSF o las inspecciones visuales de las trayectorias de MD en consideración de los cambios conformacionales más grandes que ocurren en los subsitios S2/S4.

Con el fin de realizar una revisión en profundidad de las trayectorias MD, llevamos a cabo un análisis de componentes principales (PCA) para carbonos alfa y análisis de enlaces H para algunos complejos. Los gráficos PCA para los complejos 6LU7, 6LZE, 6Y2F y 7K6D (Figura S12) mostraron que el SARS-CoV-2 MPro sufrió cambios conformacionales, especialmente entre los 500 y los 1000Â ns. Los autovectores correspondientes a los componentes principales se trazaron utilizando diagramas de puercoespín (Figura S13, S14, S15, S16) para ilustrar de manera visual cuáles son los movimientos de proteínas que ocurrieron en las simulaciones MD. Encontramos que estos movimientos pertenecen a los extremos terminales, dominio III y subsitios S2/S4 lo cual concuerda con las zonas flexibles encontradas en las gráficas RMSF, siendo los movimientos de extremos terminales y dominio III los más prominentes. Los movimientos observados se correlacionan con las rotaciones en el sentido de las agujas del reloj o en el sentido contrario a las agujas del reloj del dominio III y los subsitios S2/S4, lo que concuerda con los resultados de estudios MD similares.33. Además, la magnitud de estos movimientos varía entre los dos protómeros, lo que sugiere un comportamiento asimétrico en las unidades estructurales que constituyen el dímero. Paralelamente, se generaron superposiciones de los marcos de simulación MD y la estructura de rayos X para ciertos complejos proteína-ligando (Figura S22), lo que mostró que el dominio III en forma de dímero se mantiene a una distancia aproximadamente constante del dominio II, que es en línea con otros estudios34. Por otro lado, se estudiaron los enlaces H involucrados en las interacciones proteína-ligando para los complejos 6LU7, 6Y2F, 7K6D y 0026 (Figura S23). Los enlaces H con mayor frecuencia de aparición pertenecían a los residuos Glu-166, Phe-140, His-163, His-164, Cys-145 y Gly-143. Finalmente, aunque está fuera del alcance de este artículo, se informan resultados interesantes basados ​​en el análisis de enlaces H para residuos de dedos N entre protómeros en la sección Información complementaria denominada: Análisis de enlaces H de N-dedos.

Clasificación de afinidad MM/PBSA de inhibidores de coronavirus

Realizamos un protocolo secuencial de cálculos de energía de unión de proteína-ligando de MM/PBSA basados ​​en trayectorias de MD, a partir de las cuales se creó una clasificación de afinidad (Tabla S15) que contiene 23 activos. Estos compuestos se dividen en 4 inhibidores no covalentes y 19 covalentes. La extensión de las trayectorias MD de 500 ns no mostró diferencias de clasificación apreciables, por lo tanto, los valores promedio de energía de unión proteína-ligando se mantuvieron aproximadamente constantes durante los 500 ns adicionales. Un conjunto de 10 inhibidores que carecen de pose cristalográfica pero que exhiben excelentes valores de energía de unión proteína-ligando (por debajo de -25 kcal/mol) pertenecen a la clasificación. Un punto a destacar es la diversidad estructural de los compuestos activos identificados, considerando que entre estos compuestos identificamos inhibidores con aceptor de Michael, cetonas, ojivas de aldehídos e inhibidores no covalentes con algunos andamios basados ​​en productos naturales.

El procedimiento de agrupamiento de PCA permitió la extracción de 181 poses de ligando de las 23 trayectorias MD (Tabla S16), y se generaron diagramas de interacción ligando-proteína con porcentajes de frecuencia de interacción para los 23 complejos proteína-ligando (Figura S36, S36, S37 , S38, S39, S40). Una comparación de los diagramas nos permitió encontrar los residuos con el mayor número de interacciones proteína-ligando y la naturaleza de estas interacciones (Fig. 4A). Así, clasificamos el SARS-CoV-2 MPro residuos del sitio activo en 6 zonas: agujero de oxianión S1′, hidrofóbico S1′, S1, S2, S3/S4 HBA/HBD y S4 hidrofóbico. La zona del agujero de oxianión S1′ está compuesta por Cys-145, Ser-144 y Gly-143, donde el ligando interactúa con los grupos del esqueleto N-H a través de enlaces de hidrógeno (Figura S44). Para el subsitio S1′, denominamos «S1′ hidrofóbico» a una zona compuesta por residuos His-41, Thr-25, Thr-26 y Leu-27 que interactúan con ligandos empleando contactos hidrofóbicos (Figura S44). Cerca del subsitio S1′ se ubica S1, compuesto por His-172, Leu-141, Asn-142, Phe-140, His-163 e His-164. La mayoría de los residuos de S1 interactúan con el ligando a través de enlaces de hidrógeno, excepto Leu-141 e His-172 que interactúan a través de contactos hidrofóbicos (Figura S42). El residuo de Asn-142 interactúa como HBA o HBD con el ligando a través de la cadena lateral, mientras que Phe-140 interactúa mediante el esqueleto. Los residuos His-163 e His-164 interactúan utilizando la cadena lateral y la columna vertebral, respectivamente. El subsitio S2, constituido por Met-49 y Met-165, interactúa con el ligando utilizando contactos hidrofóbicos (Figura S41). Además, encontramos que el residuo His-41 a veces se ubica cerca de Met-49, lo que ayuda en los contactos hidrofóbicos para S2. Para los subsitios S3 y S4, generamos dos clasificaciones: la primera se llama S3/S4 HBA/HBD y está compuesta por Glu-166, Gln-189, Gln-192 y Thr-190. El segundo se denomina hidrofóbico S4 y está constituido por Leu-167, Pro-168, Ala-191 y Ala-193. Los residuos Glu-166, Gln-189, Gln-192 y Thr-190 mostraron la capacidad de actuar como HBA y HBD. Glu-166 interactúa principalmente usando su columna vertebral, aunque ocasionalmente interactúa usando su cadena lateral (Figura S41). Gln-189 interactúa por su cadena lateral, mientras que Gln-192 y Thr-190 lo hacen tanto por su cadena lateral como por su columna vertebral (Figura S43). Finalmente, los residuos Leu-167, Pro-168, Ala-191 y Ala-193 interactúan con el ligando a través de contactos hidrofóbicos, aunque ocasionalmente están presentes enlaces de hidrógeno con la columna vertebral (Figura S43).

Figura 4
Figura 4

(A) Residuos clave para la interacción proteína-ligando obtenidos a través de simulaciones MD y cálculos MM/PBSA. HBA (rojo), HBD (verde), HBA/HBD (azul) e hidrofóbico (amarillo). (B) Perfil de contribución del subsitio de la energía de unión proteína-ligando para los complejos (A) 6LU7, (B) 7K6D, (C) 6Y2F y (D) 0026 normalizado por el número de residuos. S4 hidrofóbico (Leu-167, Pro-168, Ala-191 y Ala-193), S3/S4 HBA/HBD (Glu-166, Gln-189, Thr-190 y Gln-192), S2 (Met-49 y Met-165), S1 (Phe-140, Leu-141, Asn-142, His-163, His-164 y His-172), S1′ oxyanion hole (Gly-143, Ser-144, Cys-145) , S1′ hidrofóbico (Thr-25, Thr-26, Leu-27 y His-41).

Para reconocer qué residuos son los más energéticamente relevantes en la interacción proteína-ligando, se realizó una descomposición de la energía de unión proteína-ligando para algunos complejos con un enfoque especial en las zonas ya mencionadas (Figura S25). Se sumaron los aportes por residuos pertenecientes a las 6 zonas (Tabla S18 y Figura S46) para obtener los aportes energéticos de las zonas (Figura S45). La comparación energética entre zonas puede ser desigual ya que el número de residuos que integran las zonas no es el mismo. Para compensarlo, el los valores de energía se normalizaron por el número de residuos, lo que resultó en la Fig. 4B. La zona que más contribuye a la energía de unión proteína-ligando es S2, donde Met-165 es el residuo con la energía de unión proteína-ligando más alta para todos los complejos. Por otro lado, la segunda zona más contribuyente es S4 hidrofóbica, siendo Pro-168 el residuo con los valores de energía de enlace más bajos (mayor contribuyente). Las zonas restantes parecen tener contribuciones similares excepto S1′ hidrofóbica que es la zona con la contribución más baja a la energía de enlace. Apuntando a encontrar la relación entre el movimiento de los carbonos alfa y el promedio de contribución de unión de proteína-ligando, se generaron gráficos de contribución de unión de energía libre de proteína-ligando para los sistemas 6LU7, 7K6D, 6Y2F y 0026 junto con los valores RMSF para los carbonos alfa representados en el formato de un mapa de calor (Figura S26). Estas parcelas reafirmaron que los subsitios S2 y S4 son los que presentaron mayor movilidad y aporte a la energía libre (normalizado por el número de residuos).

Una vez identificados los residuos y zonas más relevantes del sitio activo, se evaluó la estabilidad de los complejos mediante el uso de gráficos de tiempo y energía de enlace. Inicialmente, se obtuvieron gráficos de energía de unión frente a tiempo para los 23 complejos (Figura S27, S28, S29, S30). Estos gráficos muestran que la energía de enlace permaneció oscilando alrededor de valores constantes para algunos sistemas, mientras que para otros se observaron cambios. A pesar de estos cambios energéticos, las energías de estos sistemas permanecieron en valores negativos, lo que indica que los complejos continuaron exhibiendo una buena estabilidad. Estos cambios se debieron a variaciones de las poses de los ligandos, que se confirmaron después de una revisión visual de las trayectorias de MD. Debido a que los perfiles de contribución indican un promedio de la contribución energética a lo largo del tiempo, se generaron gráficas de contribución vs tiempo para los residuos de las 6 zonas (Figura S31, S32, S33, S34). Esto se hizo con el fin de inspeccionar la estabilidad de cada interacción individual. Como resultado, se puede apreciar que la mayoría de las interacciones se mantuvieron estables, a diferencia de algunos residuos pertenecientes al subsitio S4/S1′, que mostraron una pérdida de interacciones en ciertos momentos durante la simulación MD.

Modelado de farmacóforos

Para clasificar las 181 poses de ligando en varias conformaciones de proteína-ligando, las estructuras de proteína-ligando correspondientes a cada una de las 181 poses de ligando se extrajeron de las trayectorias de MD y se realizó un análisis de componentes principales en las coordenadas atómicas C-α del sitio activo. . A continuación, los tres componentes con el valor propio más alto se utilizaron para un proceso de agrupamiento con el método k-means (Figura S7), y finalmente se obtuvieron 18 agrupaciones de conformación de sitio activo. Las estructuras representativas de cada grupo se muestran en la Figura S47.

Después de dilucidar, agrupar y fusionar 181 modelos de farmacóforos iniciales (Figura S48), obtuvimos 18 representantes para modos de interacción no covalentes (Figura S49). Las funciones con una frecuencia de aparición baja o una contribución baja a la energía de enlace se configuraron como funciones opcionales. Además, una característica farmacofórica (HBA) cerca del agujero de oxianión se cambió a la característica del punto de unión de residuos (Figura S50). Esto está respaldado por el hallazgo de diferencias notables entre estos inhibidores covalentes y no covalentes en la flexibilidad del sitio activo de las trayectorias de MD. Los modelos que tienen la característica HBA ubicada en el orificio de oxianión se denominaron modelos no covalentes y los modelos que tienen la característica del punto de unión residual ubicada en el orificio de oxianión se denominaron modelos covalentes. La función de punto de unión de residuos detecta cabezas electrofílicas (ojivas) como aldehídos, cetonas, nitrilos o aceptores de Michael. Los 18 modelos covalentes y no covalentes se sometieron a un proceso de validación para medir la capacidad de los modelos para discriminar entre compuestos activos y señuelo utilizando métricas como ROC-AUC (área bajo la curva ROC) y BACC (precisión equilibrada). Los 23 activos de los cálculos de MM/PBSA se clasificaron en inhibidores covalentes y no covalentes según la presencia o ausencia de fragmentos reactivos como aceptores de Michael, aldehídos o cetonas en la posición P1′ del inhibidor. Dado que estos fragmentos reactivos son capaces de aceptar HB, los 23 activos se clasificaron en 23 ligandos no covalentes y 19 covalentes. El proceso de validación estuvo compuesto por dos pasos: (1) reducción de la rigidez de los modelos aumentando las tolerancias de los esferoides para filtrar más fácilmente los modelos con mal desempeño, y (2) disminuyendo la tolerancia de los esferoides y optimizando los valores de peso para generar modelos más ajustados. Los resultados de la primera validación se informan en la Tabla S19 y la Tabla S20. Los cinco modelos que presentaron el mejor desempeño fueron sometidos a la segunda validación como se ilustra en la Tabla S21 y la Tabla S22.

Con respecto a la ubicación de la característica para los modelos de farmacóforos, observamos que las características tienen una gran similitud con las zonas de interacción que se encuentran en la Fig. 4. Las características ubicadas en S2 son principalmente hidrofóbicas considerando que los inhibidores presentan anillos aromáticos o grupos alifáticos que pueden posicionarse en el subsitio S2. Para las características del subsitio S1, encontramos que son en su mayoría donantes/aceptores de HB, como consecuencia de la existencia de fragmentos de lactama con funciones carbonilo y amina que pueden actuar como aceptor y donante de HB respectivamente. En el caso de los subsitios S3/S4, notamos la presencia de características hidrofóbicas y de donantes/aceptores de HB. En última instancia, para el subsitio S1′, encontramos HBA y características hidrofóbicas ubicadas cerca del agujero de oxianión e His-41, respectivamente. Comparando con las contribuciones de unión de proteína-ligando de MM/PBSA (Tabla S18, Fig. 3B, S45, S46) las características más relevantes para la energía de unión con porcentajes de alta frecuencia (Figura S35, S36, S37, S38, S39, S40) son: el grupo hidrofóbico en S2 (Met-49, Met-165), el grupo hidrofóbico en S4 (Pro-168), HBA o punto de unión residual en el agujero de oxianión (Cys-145), HBD asociado a Phe-140 y HBD /HBA relacionado con Glu-166. Comparación de nuestros modelos de farmacóforos con otros estudios de modelado de farmacóforos35 encontramos la presencia de algunas características comunes relacionadas con los residuos de agujeros de oxianión, residuos de S2, His-164 y Glu-166.

Los mejores 5 modelos covalentes y no covalentes se muestran en la Fig. 5, S51 y S52. Identificamos que el modelo de farmacóforo no covalente (NCM-1) con el mejor rendimiento tiene un valor ROC-AUC de 0,75 y un valor BACC de 0,74. Por otro lado, el mejor modelo de farmacóforo covalente (CM-1) presentó un valor ROC-AUC de 0,93 y un BACC de 0,86. El CM-1 y el NCM-1 comparten 10 características farmacofóricas idénticas: 4 HBA, 4 HBD y 3 grupos hidrofóbicos (Fig. 6). Entre estas 10 características farmacofóricas, 5 características se configuran como opcionales, a saber, 3 HBA, 1 HBD y 1 grupo hidrofóbico. Además, para probar la influencia de las características en los valores ROC-AUC, se realizaron dos pruebas de los modelos en la Fig. 5: en la primera se eliminaron todas las características configuradas como opcionales (Tabla S23, S24), y en la segunda se eliminaron las características pertenecientes al agujero de oxianión (Tabla S25, S26). Como resultado de la primera prueba, se determinó que la eliminación de las características opcionales provoca una pérdida total de la capacidad de los modelos para detectar moléculas, con valores ROC-AUC entre 0,5 y 0,0. En cambio, para la segunda prueba se obtuvo una disminución de los valores de ROC-AUC, aunque no tan drástica como en la primera prueba. Esto significa que los modelos conservan cierta capacidad predictiva sin la presencia de la característica del agujero de oxianión.

Figura 5
Figura 5

Modelos de farmacóforos basados ​​en proteínas y ligandos y curvas ROC. Los modelos están ordenados de mejor a peor según los valores ROC-AUC y BACC. Modelos no covalentes (NCM) a la izquierda y modelos covalentes (CM) a la derecha. HBD (verde), HBA (rojo), hidrófobo (amarillo), punto de unión de residuos (naranja) y función opcional

. Cada modelo de farmacóforo corresponde a un grupo de conformación de sitio activo, las convenciones se pueden encontrar en la Figura S49 y la Figura S50.
Figura 6

figura 6

Los mejores modelos de CM-1 covalente y NCM-1 no covalente según los valores ROC-AUC y BACC. HBD (verde), HBA (rojo), hidrófobo (amarillo), punto de unión de residuos (naranja) y función opcional . Los residuos informados en los modelos fueron identificados por LigandScout. Con respecto a la detección de farmacóforos contra Nirmatrelvir, los resultados se pueden encontrar en las Tablas S27 y S28. Los modelos covalente y no covalente mostraron un buen rendimiento con funciones de puntuación de farmacóforo superiores a 70 en todos los casos, detectando Nirmatrelvir como un verdadero positivo. Además, el modelo CM-2 pudo replicar la pose experimental del Nirmatrelvir (Fig. 7). El grupo nitrilo estaba ubicado en S1′, la lactama en S1, el fragmento bicíclico hidrofóbico en S2, el esqueleto peptídico en S3/S2 y el trifluorometilo en S4. Las interacciones relacionadas con las características del farmacóforo corresponden a las observadas en el conjunto de ligandos utilizados, que son: S1 (Phe-140 e His-164), S1′ (Gly-143), S2 (Met-49 y Met-165) , y S3 (Glu-166). Monitoreo de los fragmentos electrofílicos correspondientes a los inhibidores covalentes de nuestro conjunto de ligandos, en conjunto con el Nirmatrelvir molécula, encontramos diferenciación entre las categorías de aceptor de aldehído, nitrilo, cetona y Michael. Los inhibidores basados ​​en aldehídos y nitrilos no satisfacían la característica hidrofóbica de S1′ debido a la ausencia de grupos hidrofóbicos en la posición P1′ del inhibidor, lo que no es el caso de las cetonas y los aceptores de Michael.

Figura 7
figura 7

Pose de ligando de Nirmatrelvir obtenida a través de CM-2 a la izquierda. Pose de ligando experimental de Nirmatrelvir a la derecha. (PDB 7SI9) El grupo nitrilo en la pose experimental no aparece correctamente ya que está unido covalentemente a la proteasa. HBD (verde), HBA (rojo), hidrófobo (amarillo), punto de unión residual (naranja) y volumen de exclusión (gris).

Deja un comentario