Modelos de predicción y refinamiento estructural para TYK2

El modelo estructural que utiliza el enhebrado para conectar los diferentes fragmentos cristalizados en los dominios proteicos recurre al programa predictor de estructura de pép-tidos y proteínas Phyre2 (^{Kelley et al., 2015}), el cual pertenece al grupo de bioinformática estructural del Imperial College de Londres, en donde se realiza un ensamblaje por homología a partir de una comparación con las estructuras cristalizadas ya reportadas y la secuencia primaria a la que se asemejan. Este programa unifica los fragmentos encontrados por homología, concibiendo el esqueleto atómico como un sistema flexible que obedece las leyes dinámicas de Newton en la conceptualización clásica (Kelley et al., 2015; Kelley & Sternberg, 2009).

La TYK2 mutada en la posición 684 del dominio pseudoquinasa 1 utiliza la misma metodología de modelamiento molecular para la enzyma wilde type, que genera el archivo en formato PDB. Como parte de su algoritmo de ensamblaje, el programa Phyre2 hace un refinamiento previo, sin embargo, se emplean programas externos para el refinamiento de los modelos hipotéticos, lo que evita los posibles solapamientos entre átomos y mejora la distribución espacial acorde con la dinámica de los campos de fuerza que afectan el modelo, para así constatar si sus parámetros tienden a cambiar mediante interacciones no covalentes del esqueleto carbonado en el modelo peptídico que se manifiesta en la estructura secundaria; para ello, se utilizan los programas de fragmentos guiados de dinámica molecular Modrefiner y FG-MD (^{Xu & Zhang, 2011}; Zhang et al., 2011).

Caracterización energética y estereoquímica de los modelos hipotéticos de TYK2

El modelo hipotético sin refinar obtenido cuenta con fragmentos que se crearon ab-initio, los cuales deben adoptar una distribución tridimensional que permita ilustrar la plausible funcionalidad de los dominios activos y no activos de la enzima TYK2. Para caracterizar y evaluar los modelos hipotéticos, se hizo una caracterización energética y estereoquímica que permitiera seleccionar el mejor modelo, sin refinar o refinado, con el fin de emplearlo en el estudio funcional de las variantes contrasentido. En la caracterización energética se hace uso del parámetro de caracterización QMEAN6 (^{Benkert et al., 2008}), el cual define un valor normalizado de la combinación lineal de términos relacionados con la estabilidad del modelo. La caracterización estereoquímica, por su parte, se hace con el programa RAMPAGE (Crystallography and Bioinformatics Group, 2017), herramienta del grupo de bioinformática y cristalografía que permite medir los ángulos diedrales Phi (φ) y Psi (ѱ) (que conectan el esqueleto carbonado de los enlaces peptídicos) en el análisis de la movilidad proteica.

Visualizadores estructurales - alineamiento tridimensional

Para estudiar las coordenadas espaciales y atómicas, se utilizaron varios visualizadores estructurales de los programas Deepview/Swiss-PdbViewer versión 3.7 (^{Guex & Peitsch, 1997}), los cuales permiten obtener un enfoque preliminar y topológico del estudio. Las estructuras se visualizaron en el programa U.C.S.F. Chimera versión 1.11 (^{Pettersen et al., 2004}), y cada estructura secundaria se modeló en un formato del tipo de cintas y estructura de densidad. El mejor modelo obtenido se alineó con los cristales reportados en el PDB a partir de un alineamiento tridimensional con la herramienta Matchmaker, la cual considera la estructura de referencia y el modelo hipotético, usando el algoritmo de alineamiento global Needleman-Wunsch y la matriz BLOSUM62 para la secuencia de aminoácidos, con lo que se obtuvo una representación tridimensional del modelo que permite observar si hay coincidencias o no a nivel espacial en esa comparación estructural, dado que en la literatura se reportan estructuras que conservan su actividad enzimática, por lo que se usa este alineamiento global para entender los efectos a nivel funcional (^{Carter et al., 2015}; ^{Zhao et al., 2015}). Estos alineamientos se cuantificaron a partir de la medición de la desviación estándar RMSD.

Análisis bioinformático de las variaciones en las proteínas - cálculo de fosforilación

La caracterización empleada en el seguimiento de las variantes a nivel proteínico toma en cuenta la estimación de la fosforilación mediante el cálculo con la secuencia de la biblioteca nativa y la mutación Ile684Ser y considera sus puntos de fosforilación en el grupo hidroxilo de aminoácidos polares neutros, como la serina y la treonina, o los hidroxilos de grupos aromáticos, como el fenol en la tirosina. Las puntuaciones en 15 quinasas, en las que se resaltan PKA, PKC, PKE, RSK, EGFR o MAPK38, se evalúan para entender las posibilidades de fosforilación a lo largo de las posiciones en la proteína (^{Artimo et al., 2012}; ^{Blom et al., 2004}).

Modelación de las variantes con el método híbrido de mecánica cuántica y mecánica molecular (QM-MM)

En el análisis funcional de la variante en la enzima TYK2, se utilizó una simulación del modelo tridimensional obtenido para las enzimas silvestre y mutada. En los modelos se aplica el nivel de teoría acorde con la capacidad del hardware disponible, el tiempo computacional, y la cantidad de núcleos y de electrones presentes en el sistema poliatómico. En la construcción de la simulación se hace una parametrización del modelo en el archivo de entrada, saturando con hidrógenos las regiones vecinas al sitio de mutación, lo que en este caso se hizo en el análisis cuántico enfocado en la posición 684 considerando las hojas beta desde Val640 hasta Asp645 y desde Gly680 hasta Glu688. Los cambios se monitorean en las estructuras tridimensionales mediante la matriz topológica. Para estos sistemas macromoleculares se utiliza el método semiempírico con la base Austin Model 1 (AM1) del programa Spartan18 de Wavefunction de medición de cambios energéticos, distancias de enlace y cambios en los ángulos (^{Foresman & Frisch, 1996}; ^{Dewar et al., 1993}), ya que la base AM1 funciona muy bien en los métodos híbridos de estudio de sistemas enzimáticos (^{Cano et al., 2021}; ^{Cardona-Pemberthy et al., 2018}; ^{Grillo et al., 2020}; ^{López-Rivera et al., 2022}; ^{Peng et al., 2019}; ^{Soto-Ospina et al., 2021a}; Soto-Ospina et al., 2021b; ^{Tejada-Moreno et al., 2022}), que provienen de análisis cristalográficos y experimentales que favorecen la parametrización.

En la simulación se hace una optimización geométrica de los diversos confórmeros energéticos y se obtiene el de menor energía, el cual sirve de base para el enfoque analítico en una modelación superficial de la estructura electrónica utilizando la herramienta "superficie" para generar diversas superficies a manera de mapas de potenciales electrostáticos en rangos de -200kJ/mol a 200kJ/mol, y obtener la distribución electrónica y la posibilidad de fosforilar las vecindades de la región que sufre la mutación en escalas de ocho colores que representan el potencial (Wavefunction, 1991). El ion fosfato considerado como el responsable de la fosforilación en el dominio pseudoquinasa 1 debido al aminoácido Ser (S), se optimizó de manera geométrica con el método del funcional de densidad DFT en un ambiente acuoso con la segunda esfera de solvatación y el método wB97x-D con la base 6-31G*. Se midieron las superficies equipotenciales, los mapas de potencial electrostático y los orbitales HOMO-LUMO, con el fin de caracterizar el comportamiento del ion considerando las escalas de distribución del potencial.

Acoplamiento molecular

El estudio de la interacción entre el ion fosfato y la pseudoquinasa 1 de la enzima TYK2 wild type y la mutada Ile684Ser, se hizo con el programa Autodock Vina que permite el estudio conformacional y el acoplamiento receptor-ligando, y genera la calificación de la interacción y la desviación estándar entre los confórmeros. Se usaron los modelos construidos con los predictores de estructura y el ion fosfato se obtuvo de Pubchem en formato .SDF para realizar la simulación.

Modelo tridimensional - caracterización energética y estereoquímica

El modelo tridimensional se realizó a partir de la secuencia primaria de aminoácidos contemplando el total de 1.187 unidades. Esta secuencia primaria se caracterizó con el alineamiento T-Coffee 2.0 y se obtuvo un valor de 99,92 % de identidad y una desviación estándar RMSD de 2.970, producto de la variante de isoleucina (I) por serina (S) en la posición 684, como se muestra en la figura suplementaria 1,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341.

Con el consenso de las secuencias primarias se construyó el modelo tridimensional utilizando el programa de predicción de estructura Phyre2, el cual se refinó utilizando algunas herramientas de medición de las interacciones no covalentes de los esqueletos carbonados a partir de campos de fuerza intrínsecos en el algoritmo de cada programa y alterando el plegamiento tridimensional para evaluar el efecto con respecto a los no refinados. Cada modelo encontrado en el estudio se caracterizó con herramientas energéticas y estereoquímicas, cuyos valores se encuentran reportados en la tabla 1.

Tabla 1 Estimaciones del modelo hipotético propuesto con métodos de caracterización estructural 

Los valores reportados que se obtuvieron de la caracterización de los modelos refinados y sin refinar, se lograron utilizando el QMEAN6 en gráficos de Ramachandran de los ángulos diedrales, de los carbonos alfa y de los enlaces peptídicos; estos evidenciaron que el mejor modelo para la enzima nativa se determinó con el Phyre2 sin refinar, con un valor definido dentro de un rango normalizado, la estimación obtenida en QMEAN6 fue de 0,666, en tanto que los de los ángulos diedrales Phi (φ) y Psi (ѱ), cuya constitución ubica los residuos en los cuadrantes según su estructura secundaria, los valores generados fueron: (95,4 %) en la región favorable, (3,2 %) en la región permitida y (1,4 %) en la región prohibida, estos valores se ubican en las regiones de mayor importancia para los residuos que constituyen el total de la quinasa.

También se desarrolló el análisis de la variante, en la cual el mejor modelo obtenido se determinó con el Phyre2 sin refinar, con un valor de QMEAN6 de 0,662 y un gráfico de Ramachandran con características de (95,8 %) en la región favorable, (2,8 %) en la región permitida y (1,4 %) en la región prohibida, lo que generó un modelo de características muy parecidas en cuanto a la predicción estructural y la caracterización energética y estereoquímica. El predictor de estructura Phyre2 hace un refinamiento previo y las herramientas de refinación no mejoran los parámetros energéticos y estereoquímicos, por esa razón, se utilizan los modelos construidos por Phyre2 mejor calificados y sin refinamiento mediante los programas Modrefiner o FGMD.

El mejor modelo silvestre (Figura1A) y el mejor de la variante (Figura1B) se alinearon a nivel tridimensional para visualizar solo el modelo silvestre y compararlo con las plantillas cristalizadas (Figura1C).

Figura 1 Representación estructural de los modelos obtenidos por predicción estructural. A) Modelo hipotético de TYK2 silvestre. B) Modelo hipotético de la variante Ile684Ser en la enzima TYK2. C) Alineamiento tridimensional con plantillas cristalizadas (Inferior izquierda); Plantillas cristalizadas (Inferior derecha) 

Luego de su caracterización, los modelos obtenidos a nivel tridimensional se alinearon con las plantillas de las regiones cristalizadas, lo que mostró poca desviación estándar entre ellos y mantuvo la topología espacial, validando así el modelo artificial al unificarlo a partir del enhebrado para contemplar la estructura completa de la enzima TYK2. El mismo tratamiento se aplicó para la variante encontrada en el dominio pseudoquinasa 1, que también tuvo un muy buen alineamiento tridimensional de la estructura, como se muestra en el material suplementario (Figura suplementaria 2,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341).

Cuando se determina una variante en los exones que codifican para la enzima TYK2, se hace una comparación del modelo a nivel tridimensional, esto para detectar cambios en la topología con respecto a la variante y así entender cómo se afecta la funcionalidad de la enzima a través del algoritmo Needleman Wunschman para alineamientos globales entre ambos modelos (Figura 2A). Asimismo, se hizo el cálculo de los enlaces de hidrógeno en el programa teniendo en cuenta los parámetros de detección de enlaces de hidrógeno intramolecular e intermolecular a lo largo de la estructura, con una restricción de 2.7 Å en la distancia y 20° en la torsión angular; estos enlaces de hidrógeno se representan en la estructura con líneas rojas dentro del alineamiento de los modelos y del fragmento pseudoquinasa 1 donde se encontró la mutación relacionada Ile684Ser, como se muestra en la figura 2B.

Figura 2 Alineamiento tridimensional de los modelos predichos. A) Alineamiento estructural global en representación en cintas. B) Cálculo de enlaces de hidrógeno (líneas en rojo) para los modelos en representación en cintas: la región pseudoquinasa 1 donde ocurre la mutación para ambos modelos se amplificó en los recuadros. 

En el análisis cuantitativo del alineamiento estructural global, se encontró que la RMSD del alineamiento por pares fue de 1.575 Å entre los dos modelos de predicción, el silvestre y el mutado. Al observar la mutación y encontrar las posibles diferencias funcionales, se analizó el resultado del cálculo del enlazamiento de hidrógeno para determinar cómo al centrarse en la región pseudoquinasa 1, donde se localiza la mutación para la posición Ile684Ser, se esperaba que el cambio de un grupo funcional hidrocarbonado alifático por un grupo hidroxilo modificara el perfil de enlazamiento de hidrógeno en las vecindades de la posición 684.

Los resultados permiten observar que la comparación cuantitativa del perfil de enlazamiento dentro de las láminas beta evidenció la aparición de dos enlaces de hidrógeno que no estaban en la enzima silvestre. Dichos cambios se deben a las interacciones entre los aminoácidos N683 y S684, los cuales conciben distancias de enlace entre el donador y los átomos donadores, como el nitrógeno y el oxígeno, de 4.886 Å y 4.580 Å, respectivamente (Tabla suplementaria 1,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341).

Fue posible corroborar la aparición de dos líneas de enlace de hidrógeno que produjeron un cambio topológico en la estructura secundaria del dominio pseudoquinasa 1, lo que podría generar cambios en los diagramas de superficie de distribución de cargas e, incluso, en la estructura electrónica, y permitiría a la subunidad pseudoquinasa 1 quedar más expuesta a una posible fosforilación o interacción con la subunidad FERM y el dominio atípico de SH₂.

Con el fin de determinar los posibles cambios en las estructuras topológica y electrónica (Figura suplementaria 3,https://www.raccefyn.co/index.php/raccefyn/article/ view/1694/3341), se hizo luego una optimización geométrica a partir del método híbrido de QM-MM, utilizando una metodología semiempírica con la base Austin Model 1 (AM1) para el cálculo de la estructura optimizada, y un congelamiento a partir de una descripción de mecánica molecular para el resto de la región de la enzima, enfocando el cálculo de superficies y mapas de potencial electrostático de la región seleccionada silvestre y de la mutada Ile684Ser, con el fin de no comprometer recursos computacionales ni niveles de teoría muy altos. Además, con las estructuras optimizadas se hizo el cálculo de superficie para entender cómo ocurre el cambio del mapa de potencial electrostático y las regiones de alta o baja densidad electrónica. Los valores de potencial se muestran en la figura 3A.

Figura 3 Diagramas de superficie para la estructura electrónica de la enzima silvestre y mutada Ile684Ser. A) Mapa de potencial electrostático para la vecindad de la región 684 en la pseudoquinasa 1 (izquierda) y para la mutación Ile684Ser (derecha). B) Acoplamiento de TYK2 wild type con el ion PO₄ ^3- (izquierda) y TYK2 Ile684Ser con el ion PO₄ ^3- (derecha) 

Después de considerar los cambios en la estructura electrónica, se hizo la representación de los modelos silvestre y mutado de TYK2 en la superficie de densidad y la distribución de cargas para comparar los cambios relacionados con la variante Ile584Ser en lo que respecta a la topología de la enzima y el acceso a los diversos dominios funcionales. En la figura suplementaria 4A,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341, puede observarse que tanto la enzima wild type como la enzima mutada presentaron cambios en la tridimensionalidad de sus dominios pseudoquinasa 1.

Ahora, con el fin de entender el posible acoplamiento entre el ión fosfato responsable de la fosforilación y la región en las cercanías de la posición 684, se hizo la simulación de la distribución de su densidad electrónica a partir de un método que contempla la correlación electrónica como teoría de densidad funcional; este método genera información más precisa y sensible que el método semiempírico, pero tiene una limitante para sistemas poliatómicos y la infraestructura computacional, pues los tiempos de simulación son elevados. Se estudió también el ion en el mapa de potencial en forma de superficie sólida y transparente, como se muestra en la figura suplementaria 4B,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341.

Se hizo el acoplamiento del ion fosfato y el fragmento de estructura de la pseudoquinasa 1 y se encontró una mejor interacción entre el ion fosfato PO₃ ^4- y el aminoácido serina en la posición 684, lo que podría suponer una interacción entre ion y dipolo con el dipolo del grupo hidroxilo, mientras en la región pseudoquinasa 1 wild type se dio una posible interacción con isoleucina por dispersión de London, desfavoreciendo el acoplamiento y la plausible fosforilación, como se observa en la figura 3B.

Los resultados obtenidos permitieron constatar la existencia de un efecto en la estructura electrónica (Figura 3A), lo que lleva a inferir que las láminas beta anti-paralelas externas quedan expuestas a su superficie. La mutación favorece que la distribución electrónica en el grupo hidroxilo genere mayor exposición para una posible interacción con alguna quinasa, permitiendo su fosforilación o una eventual interacción como la que ocurre con el dominio FERM y el atípico SH₂.

Con la simulación se obtuvieron los valores energéticos de la optimización y se observó que la estructura nativa tuvo un valor de -2764k J/mol y la mutación Ile684Ser en la región adyacente a la posición 684 uno de -2982,48 kJ/mol, lo que evidenció que la mutación generó mayor estabilidad, por lo que su valor fue de menor energía que el de la silvestre. En el análisis superficial con el mapa de potencial electrostático, el grupo funcional isoleucina en la posición 684 de la enzima TYK2 silvestre no presentó potencial, lo que se explica por su naturaleza hidrocarbonada, pero cuando la mutación se produjo por la vía del aminoácido serina, esta presentó un momento dipolar neto establecido por la diferencia de electronegatividad entre el oxígeno y el carbono al que está enlazado, lo que produjo un aumento de densidad electrónica localizado en el átomo de oxígeno del grupo hidroxilo. Al cuantificar las propiedades de superficie se logró obtener un valor del volumen de distribución de carga de 465.65 ³ y un área superficial electrónica de 94.27 ². Al analizar los resultados de la superficie, también se pudo evidenciar un efecto de bloqueo por impedimento estérico dada la cadena hidrocarbonada presente, que en la mutación presenta menos bloqueo y una mayor porción de acceso a la región electrónica disponible para reaccionar como una fosforilación plausible.

Con las herramientas bioinformáticas del programa Netphos 3.1, se calculó la fosforilación para la enzima silvestre y para la mutada. Se encontró que la posición 684 de la silvestre no presentaba ninguna posibilidad de fosforilación, pero al sufrir la mutación en Ser684, la región pseudoquinasa 1 quedó sujeta a una fosforilación según el orden de probabilidad por acción de las quinasas, con un valor normalizado de (0,455) GSK3, (0,447) CaM-II, (0,439) CK-II y (0,431) DNAPK. Teniendo en cuenta los valores que no superaron el umbral de 0,5, fueron cuatro los puntos que se acercaron a superar el umbral. Al evaluar la enzima silvestre en esa posición no hubo una respuesta importante en cuanto a los valores de las quinasas, como se muestra en el recuadro de la figura suplementaria 5,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341.

Esto evidencia una posible hiperactivación de la enzima que puede desencadenar la respuesta patológica debido a la posible fosforilación favorecida por la formación del aducto con el grupo hidroxilo del aminoácido serina (Ser684). En la figura suplementaria 4,https://www.raccefyn.co/index.php/raccefyn/article/view/1694/3341 se resumen los cambios estructurales encontrados: la mutación produjo una alteración en las distribuciones de carga y un efecto sobre el dominio FERM y el atípico de SH₂, evidenciados en los cambios en el dominio pseudoquinasa 1 y en aquellos ocurridos en la superficie de densidad, que podrían afectar la función enzimática de la TYK2.

En la literatura no se encontraron reportes experimentales que confirmen la fosforilación de esta posición para la enzima TYK2, pero sí se sabe que la inhibición de esta enzima favorece el tratamiento contra la inflamación en la psoriasis, dado que es mediadora de la inflamación en la producción de IFN-y, IL-17 e IL-23. Por ello se hace una propuesta de simulación computacional como mecanismo de activación desregulada de la señalización de la enzima TYK2 a partir de una autofosforilación, lo que hasta el momento no se ha hecho a nivel experimental, por lo que este modelo explica su hiperactivación (^{Ghoreschi et al., 2021}; Zhi^{Li et al, 2013}; ^{Min et al., 2015}).