Guillermo Narváez ¹, Leonardo Lareo ², Javier Rincón ³

1 Facultad de Medicina, Maestría en Bioquímica, Universidad Nacional de Colombia, Bogotá, D.C., Colombia.

Los receptores ionotrópicos de glutamato son canales iónicos activados por ligandos expresados principalmente a nivel del sistema nervioso central, los cuales median la mayoría de los procesos de neurotransmisión excitatoria y plasticidad sináptica (18,19). De los tres subtipos de receptores ionotrópicos de glutamato conocidos, el que más interés biomédico ha tenido es el receptor de glutamato activado por N-metil-D-aspartato (NMDA), ya que la estimulación excesiva de estos receptores bajo condiciones patológicas conduce a la denominada "excitotoxicidad neuronal", y tal estimulación se ha relacionado con la pérdida de neuronas que acompaña la isquemia cerebral y ciertas enfermedades neuro-degenerativas (20). Debido a que la isquemia cerebral, ya sea por trombosis o apoplejía, es una de las principales causas de muerte a nivel mundial y la primera causa de incapacidad a largo término, este evento patológico ha sido uno de los principales blancos terapéuticos hacia los cuales se ha enfocado la investigación con el fin de poder frenar, empleando para ello antagonistas del receptor NMDA, la cascada de eventos excitotóxicos generada por la excesiva estimulación glutamatérgica.

El receptor de glutamato subclase NMDA emplea como coagonista a la glicina, ligando tan importante como el ácido glutámico para la apertura del canal (21). La investigación preclínica y clínica ha evidenciado que los agonistas parciales o los antagonistas del sitio de unión de glicina poseen mejores perfiles de efectos colaterales que los bloqueadores directos del canal, aun cuando desde el punto de vista químico ha sido un verdadero reto desarrollar antagonistas potentes que compitan con este pequeño aminoácido dentro de los límites normales de dosificación (22).

En el presente trabajo se generaron ecuaciones matemáticas que relacionan la estructura química de antagonistas de glicina con sus propiedades de afinidad por el correspondiente dominio de unión en el receptor ionotrópico de glutamato subclase NMDA. Adicionalmente, la información aportada por el modelo se utilizó para proponer nuevos antagonistas de glicina, cuya afinidad fue predicha teóricamente con las ecuaciones desarrolladas.

Para la construcción del modelo matemático se contó con una base de datos generada a partir de la revisión bibliográfica entre los años 1993 y 2004 (23), en la cual se tabularon aproximadamente 300 antagonistas del sitio de unión de glicina agrupados en 22 familias químicas, con sus correspondientes datos de afinidad expresados como CI50. El término CI50 (concentración inhibitoria 50) corresponde a la concentración molar del antagonista que inhibe en un 50% la unión de un ligando de referencia. Para el presente caso, cuanto menor sea el valor numérico de la CI50 para un antagonista, mayor será su afinidad por el dominio de unión a glicina. Como criterios de selección para los antagonistas que serían incluidos en el trabajo se tuvo en cuenta la homogeneidad experimental en la obtención de la CI50; el reporte de la precisión en la determinación experimental de la CI50, y variaciones estructurales sistemáticas sobre la región de interés del grupo farmacofórico.

Debido a que el modelo matemático tendría características semicuantitativas se definió el siguiente criterio de clasificación de los compuestos a ser empleados para construir el modelo:

Antagonistas positivos. Son aquellos con CI50 menor a 250 nM. Se considera que estos compuestos tienen una alta capacidad de unión al dominio de glicina del receptor NMDA y, por consiguiente, pueden competir y desplazar al ligando endógeno, antagonizando su efecto.

Antagonistas negativos. Son aquellos que tienen una CI50 mayor a 250 nM. Estos compuestos no presentan la suficiente capacidad de unión para ser utilizados farmacológicamente como anta-gonistas de glicina.

El umbral de 250 nM se eligió porque experimental-mente se ha encontrado que la afinidad de la glicina en receptores NMDA heteroméricos se encuentra entre 200 y 300 nM (24,25). Cualquier compuesto que deba competir con este coagonista debe tener la suficiente afinidad para contrarrestar las con-centraciones de glicina presentes en la hendidura sináptica, de tal forma que el objetivo es buscar antagonistas que posean una afinidad expresada en términos de CI50, inferior a 250 nM (26).

Para generar la conformación molecular tridimensional energéticamente más favorable de cada antagonista seleccionado se utilizó el programa Molecular Modeling Pro CHEMSW® 2003, el cual emplea herramientas de mecánica molecular para llevar a cabo esta tarea.

Se utilizaron los algoritmos ampliamente descritos en la literatura (7-12). Los índices de conectividad utilizados para el presente estudio fueron los siguientes.

Indices de conectividad simples: 0X, 1X, 2X, 3X, 4X, 5X, 6X, 3Xc, 4Xpc.

Indices de conectividad de valencia:0Xv, 1Xv, 2Xv, 3Xv, 4Xv, 5Xv, 6Xv, 3Xvc, 4Xvpc.

Las CI50 y los diferentes índices de conectividad que codifican la estructura química de cada antagonista se utilizaron para generar el modelo matemático-estadístico mediante análisis de regresión simple y múltiple, expresando la afinidad como CI50, Ln CI50, 1/CI50, y evaluando la relación lineal, logarítmica, recíproca y multiplicativa a un nivel de probabilidad del 95%. Se utilizó el programa Excel de Microsoft Office 2003.

Con el fin de tener una metodología alterna que permitiera relacionar la afinidad con otro parámetro diferente a los índices de conectividad, se investigó si existía alguna correlación entre la CI50 y la energía de interacción de los antagonistas con su correspondiente dominio de unión. Para ello se utilizaron los segmentos S1S2 del dominio de unión de glicina del receptor de glutamato NMDA, cuya estructura terciaria, determinada mediante cristalografía de rayos X, fue recientemente publicada (27,28). La estructura tridimensional empleada tenía como ligando acoplado al ácido dicloroquinurénico un antagonista de glicina prototipo. En la figura 2 se observan los residuos aminoacídicos involucrados en la interacción con el ácido dicloroquinurénico. Para determinar la energía de interacción de los antagonistas seleccionados con el dominio de unión de glicina se utilizó el programa especializado ArgusLab 4.0.1 (29). Las dimensiones de la celda de unión para estos análisis fueron de 21,37Å x 18,67Å x 14,44Å; sin embargo, se hicieron ensayos de acople con celdas de hasta 15 Å x 14 Å x 14 Å. Se emplearon celdillas con una resolución de 0,4 Å, y se activó la opción de alta resolución de acople con un número máximo de interacciones de 450 poses.

Los análisis de energía de interacción o docking se utilizaron para validar los resultados obtenidos con el modelo matemático-estadístico desarrollado con base en índices de conectividad cuando éste se aplicara a moléculas nuevas cuyas propiedades de unión al dominio de glicina se deseara conocer. De esta forma el análisis de docking se convirtió en un control del modelo propuesto mediante los índices de conectividad.

Con base en el comportamiento de los índices de conectividad que más influencia mostraron sobre la afinidad, se diseñaron nuevos antagonistas del dominio de unión de glicina, prediciendo sus propiedades de unión a través de las ecuaciones matemáticas generadas. Los datos obtenidos con el modelo se confirmaron mediante el análisis de interacción bioquímica o docking.

Todos los programas computacionales mencionados fueron ejecutados en un equipo DELL™ XPS M1710 con procesador Intel® Core™ 2 Duo T7200 (2 GHz/ 667 MHz FSB) y memoria de 1GB SDRAM DDR2 a 667MHZ, 2 DIMM (1G2D).

Con base en los criterios mencionados se seleccionaron para el estudio 45 compuestos (30-32) que presentan modificaciones estructurales sistemáticas sobre el anillo bencénico de la dihidroquinoxalina-2,3-diona ( figura 3 , cuadro 1 ), región que de acuerdo con el grupo farmacofórico propuesto hasta la fecha (24,33), influye significativamente sobre la afinidad de este tipo de antagonistas.

Generación de la conformación tridimensional energéticamente más favorable

La generación de la conformación estructural tridimensional energéticamente más favorable se llevó a cabo sin mayores complicaciones. Como se puede observar en la figura 3 y en el cuadro 1 , los antagonistas seleccionados para el estudio tienen mucha restricción de rotación en torno a los enlaces de los átomos que conforman los anillos, lo cual disminuye significativamente la posibilidad de que tengan una amplia gama de conformaciones energéticamente favorables. Los resultados de esta etapa se emplearon posterior-mente durante los análisis de interacción energética o docking.

Los índices de conectividad simple y de valencia calculados para cada uno de los 45 antagonistas de estudio se presentan en los cuadro 2 y cuadro 3 . Con estos resultados se procedió a efectuar el modelaje matemático y estadístico.

Analizando el diagrama de dispersión para el índice de conectividad 4Xvpc se observó que había un comportamiento diferencial entre los antagonistas evaluados, generándose subgrupos de antagonistas que correlacionaban muy bien entre sí. Es así como en la figura 4 se observaron cuatro subgrupos de antagonistas (líneas punteadas) cuya correlación entre afinidad, expresada como Ln CI50, y su correspondiente índice 4Xvpc era significativamente alta. Advirtién-dose este comportamiento, se decidió dividir el total de los antagonistas analizados en dos grupos: aquellos que se ubicaban por debajo de la línea de regresión inicial (línea continua en la figura 4 ) (grupo A); y aquellos que se ubicaban por encima de esta línea (grupo B). Se encontró que en el grupo A se ubicaron preferencialmente aquellos compuestos que tenían valores de afinidad menores de 250 nM (15 de 23 compuestos, correspondiente al 65,2%), es decir, aquellos que sí presentaban una unión significativa al receptor; mientras que en el grupo B, el 81,8% (18 de 22) de los compuestos tenía una afinidad por encima de 250 nM, es decir, los considerados sin afinidad por el sitio de unión de glicina. Posteriormente, mediante regresión lineal simple se analizaron los dos grupos de antagonistas generados (A y B) para la variable 4Xvpc, obteniéndose los resul-tados ilustrados en la figura 5 . Para el grupo A de antagonistas ( figura 5 , panel izquierdo), el coeficiente de determinación (r2) y de correlación (r) para la regresión efectuada fue de 0,8589 y de 0,9267, respectivamente, cumpliéndose todos los supuestos de normalidad estadística. La ecuación de la curva de regresión para este modelo fue la siguiente (ecuación 1):

Para el grupo B de compuestos ( figura 5 , panel derecho), el coeficiente de determinación (r2) y de correlación (r) para la regresión fue de 0,8761 y de 0,9359, respectivamente, cumpliéndose los supuestos de normalidad estadística. Para este análisis estadístico se excluyó el compuesto 1ant28, el cual correspondía a un punto atípico de residuales. Este compuesto también podría considerarse como un "atípico químico", pues es el único que dentro de su estructura involucra un tioeter como sustituyente del anillo bencénico. La ecuación de la curva de regresión para este modelo fue la siguiente (ecuación 2):

Sin embargo, aun cuando la predicción de la afinidad mediante la ecuación 3 es correcta para todos los antagonistas, se tenía el problema de que la variable dummy "PosNeg" dependía del valor de la afinidad, valor éste que en la práctica no se tiene para una molécula nueva. Por consiguiente, el interés se centró en construir una función que se aproximara de la mejor manera posible a la diferenciación hecha con las anteriores variables dummys, pero que ya no dependiera de la afinidad. Esta nueva variable se construyó teniendo como referencia la línea de regresión entre las variables Ln(CI50) y el índice de conectividad 4Xvpc, y el umbral de afinidad de 250 nM. Como se observa en la figura 6 , a partir de un valor de 1,18 de la variable 4Xvpc, todos los compuestos son positivos (afinidad menor a 250 nM), y por debajo de 1,18, todos los compuestos son negativos (afinidad mayor a 250 nM), quedando solamente cuatro compuestos mal caracterizados. Con estos nuevos datos se decidió construir un nuevo modelo de regresión múltiple empleando para ello una nueva variable dummy denominada "Factor de corrección de grupo, FCG", la cual toma un valor de 1 si la variable 4Xvpc es mayor a 1,18 o de 0 si esta variable es menor o igual a 1,18. La regresión para las 43 observa-ciones arrojó un coeficiente de correlación múltiple (r) de 0,8109, un valor F de 38,43, un error típico de 1,66, siendo las dos variables regresoras (4Xvpc, FCG) significativas (p<0,05). Así mismo, la regresión cumplió con los supuestos de normalidad estadística. La ecuación de la curva de regresión para este modelo fue la siguiente:

El análisis de las energías de interacción entre antagonistas de glicina y su correspondiente sitio de unión o docking se llevó a cabo para los 45 compuestos seleccionados utilizando las estructuras moleculares con la conformación tridimensional energéticamente más favorable generadas según lo descrito. Para evaluar el desempeño del programa de docking en el equipo donde fue instalado se hicieron controles de acople con proteínas y ligandos no relacionados (transportador de ácido retinoico y receptor de estradiol). Los RMSD promedio obtenidos fueron de 0,42 Å y 0,55 Å, respectivamente, con una formación de puentes de hidrógeno idéntica a la del ligando de referencia contenido en la correspondiente cristalografía de rayos X (datos no mostrados). Una vez obtenidas las energías de interacción para los antagonistas de glicina en estudio, se construyó el diagrama de dispersión ilustrado en la figura 7 . En esta figura se observa la generación de dos grupos de puntos separados por aproximadamente 400 calorías/mol, lo cual sugiere claramente que para obtener una afinidad inferior a 250 nM se requiere de una energía de interacción menor a –8,0 Kcal/mol. Desde el punto de vista estadístico, estos resultados se analizaron mediante la prueba de Fisher, con la cual se obtuvo un valor de p<0,05, indicando así que la energía de acople o de docking está estadísticamente asociada con la afinidad. Por otro lado, el coeficiente de correlación de Spearman obtenido fue de 0,999, con un valor de p<0,05 para los 32 compuestos analizados. Esto significa que el 99% de los compuestos cuya energía de interacción es menor a –8,0 Kcal/mol tendrán un valor de afinidad (Ln CI50) menor a 5,52, es decir menor a 250 nM.

Es importante mencionar que para poder catalogar el análisis de docking como correcto se tuvo en cuenta el siguiente criterio: los anillos base de cada uno de los antagonistas analizados deben haber quedado superimpuestos a los anillos base del ácido dicloroquinurénico (ADCQ) de referencia, y en la posición correcta, es decir, el anillo bencénico de la quinolina del ADCQ debe superponerse al anillo bencénico de la dihidro-quinoxalina-2,3-diona, anillo base de todos los antagonistas analizados. De igual forma, los dos anillos heterocíclicos deben haberse super-impuesto. Sólo si se observaba esta condición, se tomaba en cuenta la energía de interacción que reportaba el programa. Con base en este criterio se obtuvo una posición correcta para 32 de los 45 antagonistas de estudio, de tal forma que la figura 7 fue construida con estos 32 puntos.

De acuerdo con las anteriores condiciones, se sugirieron para la evaluación de afinidad a través del modelo matemático propuesto, las nuevas moléculas enumeradas en el cuadro 4 . En este punto es donde entran a jugar un papel decisivo los resultados obtenidos con los análisis de energías de interacción ( docking), metodología esta que se utilizó para validar los resultados obtenidos mediante el modelo de índices de conectividad. Es así como los datos generados por el modelo en cuanto a las características de unión de estas nuevas moléculas, se confirmaron a través de la determinación de la energía de interacción entre ellas y el dominio de unión de glicina para ratificar o no de esta forma los resultados obtenidos con el modelo matemático generado (cuadro 4 ).

La predicción semicuantitativa mediante la ecuación 4 del modelo matemático-estadístico desarrollado estuvo acorde con las correspon-dientes energías de interacción para 11 de las 12 nuevas moléculas evaluadas, lo cual sustenta sólidamente la validez del modelo propuesto. De acuerdo con este modelo, los compuestos P5, P6 y P9 son moléculas promisorias, candidatas para síntesis y evaluación biológica. Si bien las moléculas P3 y P4 tienen una energía de docking favorable y el valor de su índice de conectividad 4Xvpc cumple el requisito de ser superior a 1,18, éste se sale de la línea de regresión del grupo de compuestos analizados, lo cual no hace tan confiable la predicción de su afinidad como sí lo es para las moléculas P5, P6 y P9, cuyos valores del índice mencionado caen sobre la línea de regresión mostrada en la figura 6 .

Existe además el inconveniente de que muchos reportes experimentales se publican sin la suficiente rigurosidad estadística que permita evaluar la calidad de los datos presentados y, por tanto, sobre las conclusiones basadas en ellos. Es así como gran cantidad de moléculas no fueron incluidas en el estudio, pues carecían de una desviación estándar que permitiera evaluar la precisión en la obtención de la CI50. Este aspecto es importante, ya que no se le puede exigir al modelo matemático generado que prediga datos de CI50 con mayor precisión y exactitud que los que experimentalmente pueden obtenerse. El promedio de las desviaciones estándar reportadas para la determinación de la CI50 in vitro de los antagonistas de estudio ( cuadro 1 ) fue de alrededor de 16%, siendo los valores mínimos y máximos de 5% y 26%, respectivamente. Esta deficiente precisión en la evaluación de la afinidad, propia de la metodología experimental actual, puede explicar el hecho de que la estructura de un compuesto, codificada en los índices de conectividad, no explique por sí sola un mayor porcentaje de la variación en la afinidad. Pese a lo anterior, se considera que el modelo mate-mático generado es de una buena calidad estadística. Adicionalmente, el hecho de que se haya logrado seleccionar 45 moléculas que cumplieron los criterios de selección dio una base confiable para el desarrollo del presente trabajo.

La comparación del índice de conectividad 4Xvpc con su homólogo 4Xpc, el cual, a diferencia del primero, no codifica información acerca de la presencia de heteroátomos e insaturaciones, arrojó resultados estadísticos bastante pobres. Lo anterior sugiere que las insaturaciones, y en especial la presencia de heteroátomos, son factores decisivos para la adecuada interacción entre esta región y los aminoácidos presentes en el dominio de unión del antagonista en el receptor NMDA. El índice 4Xvpc es especialmente útil para codificar información acerca del número de sustituyentes sobre anillos bencénicos, los patrones de sustitución y su longitud (hasta tres enlaces), así como el tipo de sustituyente heteroatómico. Adicionalmente, este índice codifica información acerca de la orientación de los sustituyentes en el anillo: en cuanto mayor sea la adyacencia de los sustituyentes, mayor es el valor numérico del índice (9). Lo anterior puede observarse por ejemplo para los compuestos 1ant38 y 1ant35, ambos con dos halógenos como sustituyentes, donde el valor del índice 4Xvpc es mayor para el compuesto que tiene los halógenos vecinales. El incremento en el número de sustituyentes sobre el anillo aromático también incrementa el valor de este índice, lo cual se evidencia en los compuestos 1Ant36, 1Ant37, 1Ant38, 1Ant39 ( cuadros 1 y cuadro3 ).

Adicionalmente, el índice 4Xvpc tiene la capacidad de discriminar entre compuestos cuya única diferencia es el tipo de heteroátomo sustituyente sobre el anillo aromático. El índice de conectividad simple homólogo 4Xpc carece de esta capacidad (ver los compuestos 1Ant15, 1Ant32 y 1Ant27 en los cuadros 1 y cuadro3 ). Es importante anotar que cuanto mayor sea el delta de valencia para los heteroátomos en la molécula, menor es el valor numérico del índice 4Xvpc. De la misma forma, el índice 4Xvpc puede diferenciar entre sustituyentes que contienen los mismos átomos, pero en diferente posición (ver compuestos 1Ant17, 1Ant18). De allí la importancia de este descriptor para codificar la estructura química en esta clase de compuestos.