Cálculos de química Cuántica

Los cálculos se realizaron con el software ORCA en su versión Orca-5.0.3⁷. Esto se llevó a cabo utilizando el método DFT⁸. Inicialmente, la taspina fue preoptimizada con wB97M y el conjunto base def2-ZVP; luego, se optimizó su geometría con el funcional wB97M-D4 ^[9] y el conjunto base def2-TZVPP ¹⁰, todo en fase gaseosa. Los cálculos de los índices de Fukui también se realizaron con Orca-5.0.3, utilizando las coordenadas obtenidas en la optimización de la geometría. Para esto se empleó también el funcional wB97M-D4 y el conjunto base def2-TZVPP. Los inputs fueron las moléculas de taspina aniónica, catiónica y neutra. Para determinar el gradiente reducido de la densidad (RDG), se utilizaron dos softwares adicionales: VMD y Multiwfn, que es un programa de análisis de funciones de onda¹¹. Todos los cálculos mencionados anteriormente se realizaron usando una computadora personal AMD Ryzen 7 5800X3D, 3.40 GHz, 8 núcleos, 16 hilos de ejecución.

Acoplamiento molecular

Un paso previo antes de realizar el anclaje molecular fue la búsqueda de la enzima que interactúa con la taspina mediante el uso de SwissTargetPrediction ya que permitió identificar a las posibles proteínas dianas involucradas en su mecanismo de acción. La estructura elegida fue la acetilcolinesterasa (AChE-S) que se descargó del banco de proteínas (RCSB PDB) identificada con el código 4M0E, en la cual se simuló el acoplamiento molecular contra la taspina. La proteína se preparó utilizando el programa UCSF Chimera, donde se removieron las moléculas de agua. Luego, se identificaron los sitios activos de la AChE donde se simularon las interacciones del ligando y los aminoácidos terminales, empleando el servidor en línea CASTp que localizó la cavidad con los sitios activos. Seguido se utilizó la interfaz gráfica AutoDockTools para crear el archivo confing.txt, que contuvo el centro de la caja con las coordenadas en X= 1.154, Y= -41.283 y Z= 13.081, de la cavidad para el acoplamiento. luego, se creó una carpeta de trabajo donde se guardó el archivo config.txt, donde se modificó la exhaustividad a 32 y el número de conformaciones a 1, la taspina y la proteína AChE en formato pdbqt, y un script llamado docking.sh utilizado para realizar el anclaje molecular de forma repetitiva unas 50 veces. Una vez hecha la carpeta de trabajo se empleó el programa AutoDock Vina, desde la terminal de linux ejecutando los comandos chmod +x docking.sh y ./docking.sh. Por último, se realizó una simulación computacional más costosa y exhaustiva mediante el anclaje molecular flexible que nos permite seleccionar una cantidad limitada de los aminoácidos terminales de la cavidad y hacerlos flexibles para que cambien su conformación junto con el ligando durante la simulación.

Método computacional de Simulación de Dinámica Molecular (DM)

Las simulaciones de DM all-atom fueron realizadas considerando las coordenadas atómicas entre la proteína AChE-S y taspina usando Gromacs v5.1.4 con campo de fuerza AMBER99SB y modelo de agua TIP4P. La caja de simulación fue definida con una distancia de 10A° desde el centro del complejo proteína ligando. El sistema preparado fue solvatado con moléculas de agua y neutralizado con la adición de 0.15M de iones Na+ y Cl- ¹³. El proceso de minimización de energía involucró 50000 pasos por cada descenso seguido por gradientes conjugados. Las condiciones PBC fueron definidas por direcciones x, y, z¹⁴ y la simulación fue realizada a temperatura fisiológica de 300K. Se aplicó el algoritmo VERLET para restringir todos los enlaces involucrados con hidrógeno y fuerzas electrostáticas de largo alcance tratadas con PME (Particle mesh Ewald). El sistema fue equilibrado en dos pasos, NVT y NPT a 300K por periodos de 1000ps. Durante la simulación, se utilizaron el termostato V-rescale y Parrinello-Rahman¹⁵ para mantener la presión y la temperatura. Se utilizó el algoritmo LINCS para restringir los enlaces y ángulos. Las interacciones de Van der Waals fueron restringidas con el potencial LJ con un límite de 1.0 nm.

Método de gradiente reducido de la densidad (RDG)

La gradiente reducida de la densidad (RDG) es una herramienta clave en la visualización de interacciones no covalentes ^[12]. El análisis del RDG revela diversas interacciones atractivas, de Van der Waals, y repulsivas que influyen en la estabilidad y comportamiento de la molécula. El gradiente reducido de densidad (s), un campo escalar de la densidad electrónica, se describe como:

La taspina, al ser una molécula con cuatro anillos aromáticos fusionados, así como grupos éster, éter y una amina primaria, exhibe una red compleja de interacciones intra e intermoleculares. Estas interacciones pueden visualizarse mediante el índice RDG al analizar los valores del parámetro λ₂ρ, donde la densidad electrónica y su gradiente indican la naturaleza de las interacciones. Donde H-B (atracciones fuertes: tipo enlace puente de hidrógeno), VdW (Interacciones de Van der Waals) y E.E. (Repulsiones fuertes efectos estéricos).

Figura 1(a) Isosuperficie de RDG de la taspina - Figura 1(b) Gráfico 2D RGD de la taspina 

En la figura 1(a), no se observa una isosuperficie de color azul, lo que sugiere la ausencia de enlaces por puente de hidrógeno. Esto concuerda con la estructura química de la taspina, la cual no posee átomos dadores de puentes de hidrógeno, pero sus átomos de nitrógeno y oxígeno pueden ser aceptores de puentes de hidrógeno. Sin embargo, se observa otras interacciones, como las de tipo Van der Waals (VdW) y repulsiones estéricas (E.E.), estas se representan mediante isosuperficies de color verde y rojo, respectivamente.

Por otro lado, en la figura 1(b), se observa que los valores asociados a las interacciones fuertes (H-B) se encuentran entre -0.35 < sign( 𝜆 2 )𝜌 < -0.25 a.u. Mientras que las interacciones VdW están en el rango de -0.20 < 𝑠𝑖𝑔𝑛( 𝜆 2 )𝜌 < 0 a.u., lo que indica la presencia de interacciones Van der Waals en la región ecuatorial de la molécula principalmente. Finalmente, las isosuperficies de color rojo que representan las repulsiones estéricas se localizan principalmente en el centro de los cuatro anillos aromáticos de la taspina, los valores de estas interacciones oscilan entre 0.10 < 𝑠𝑖𝑔𝑛( 𝜆 2 )𝜌 < 0.20 a.u.

Descriptores Globales y Locales

En la tabla 1 se observan los descriptores globales, estos indican la reactividad general del sistema. La diferencia de energía entre el HOMO y el LUMO da un ΔE = 0.0304 a.u., lo cual sugiere una moderada estabilidad frente a la transferencia electrónica. El valor negativo de la dureza global (η) indica que la molécula es flexible frente a cambios electrónicos, lo que aumenta su reactividad; el valor de la blandura global (s) es el inverso de la dureza, su valor negativo confirma lo mencionado. El valor negativo de la electronegatividad (χ) sugiere que la molécula tiene una baja tendencia a atraer electrones. Por último, el valor de electrofilicidad global (Ψ) sugiere que la taspina tiene baja propensión a aceptar electrones, indicando que se comporta más como nucleófilo que como electrófilo.

En el análisis de la función de Fukui, al considerar los estados de carga neutro, catiónico (+1) y aniónico (-1), se obtienen los valores de la función de Fukui condensada para ataques nucleofílicos (𝑓⁻), electrofílicos (𝑓⁺) y radicalarios (𝑓⁰), los cuales proporcionan información valiosa sobre los sitios más reactivos de una molécula. Los valores positivos de 𝑓⁺ y 𝑓⁻ revelan tendencias significativas en los centros carbonados, oxígenos y nitrógenos, indicando su propensión a sufrir ataques nucleofílicos y electrófilos, respectivamente. En la tabla 2, se presentan los descriptores locales calculados para identificar sitios reactivos dentro de la molécula. Los carbonos C6 y C4 presentan valores elevados de 𝑓⁻ de 0.097 y 0.068, respectivamente, lo que indica que estas posiciones son favorables para recibir ataques electrófilos. Por otro lado, los átomos de oxígeno O15 y O7 presentan altos valores de 𝑓⁺ de 0.065 y 0.054, lo que sugiere que estas posiciones son susceptibles a sufrir ataques nucleofílicos.

Adicionalmente, los descriptores duales confirman regiones de ambivalencia química en sitios como C17 y C2, donde los valores de 𝑓⁺ y 𝑓⁻ son muy similares de 0.00001 a.u. y 0.00002 a.u., sugiriendo su capacidad para actuar como nucleófilos o electrófilos dependiendo del contexto de la reacción. Finalmente, estos descriptores locales complementan la caracterización global del sistema molecular al proporcionar información detallada de los puntos más reactivos y sus posibles interacciones, permitiendo una predicción precisa de las rutas de reacción más probables.

Estudios de anclaje molecular

La búsqueda de las posibles Dianas mediante el SwissTargetPrediction¹⁷, dio como resultado un listado de enzimas candidatas, dentro de las cuales se encuentran los receptores acoplados a proteínas G, kinasas, proteasas e hidrolasas. El programa elige en función de la similitud molecular y la alta probabilidad estadística dada por la herramienta bioinformática que seleccionó a la acetilcolinesterasa (AChE-S) clasificada como hidrolasa (figura 2).

Figura 2 Interfaz del servidor SwissTargetPrediction donde se muestra a la taspina con los diferentes tipos de Dianas seleccionadas. 

El análisis de la estructura de la AChE empleando el servidor en línea CASTP¹⁸, evidenció una cavidad cuyo sitio activo presentó un volumen de 933 ų y una superficie de 599.6 Ų (Figura 3). En esta región con el programa UCSF Chimera¹⁹, se identificó a los aminoácidos terminales activos que aparecen de color verde en la figura. Estos aminoácidos fueron ASN 233, PRO 235, GLU 313, ILE 316, ASN 317, VAL 367, VAL 370, HIS 405, CYS 409, PRO 410, ALA 412, GLN 413, GLY 416, ARG 417, TYR 503, ALA 505, GLN 508, LEU 524, ALA 526, CYS 529, ALA 530, TRP 532, ASN 533, ARG 534, LEU 536, PRO 537, LYS 538 y LEU 540.

Figura 3 (a) Cavidad seleccionada por el CASTP de la AChE de color rojo. Figura 3 (b) Secuencia completa de aminoácidos de la AChE-S, donde los de color verde corresponden a los aminoacidos terminales activos.  

Respecto a las simulaciones por anclaje molecular que se realizaron 50 veces entre la taspina y la enzima acetilcolinesterasa presentaron en promedio una energía mínima de unión -7.1 Kcal/mol, y al realizar el anclaje molecular flexible la energía llegó hasta los -8.71 Kcal/mol, lo cual significa que la formación del complejo taspina/acetilcolinesterasa es estable, favorable y espontánea. El análisis post docking se realizó con el Discovery Studio Visualizer V24 ^[20], con el cual se identificaron interacciones más importantes en el sitio activos correspondientes a los aminoácidos PRO 290, ARG 296, VAL 239, GLY 240, PRO 235 y HIS 405 (figura 4). Así como la estructura de la superficie hidrofóbica (figura 5) donde se aprecia la cavidad en la cual se fijó la taspina durante el acoplamiento molecular.

Figura 4 En la imagen de la izquierda se muestran las interacciones de los aminoácidos terminales más representativos con la taspina. En la imagen de la derecha aparecen los diferentes tipos de interacciones no covalentes más importantes en el momento del acoplamiento. 

En el programa también se identificaron interacciones no covalentes que influyeron en el valor de la energía mínima de unión en el complejo (Figura 4). Las más importantes fueron el convencional puente de hidrógeno, enlace de carbono-hidrógeno, alquil y pi-alquil. En la tabla 2, se registraron junto con las distancias a la cual actuaron frente al ligando.

Figura 5 Superficie hidrofóbica de la AChE donde se muestra la cavidad donde se acopló la taspina. 

Estudios de simulación de Dinámica Molecular (DM)

Se utilizó la mejor posición del análisis de docking molecular y se realizó la simulación de DM all-atoms en un periodo de 15ns a temperatura fisiológica. Las interacciones de la AChE y taspina se analizaron mediante la investigación de los gráficos de RMSD, Rg, SASA y RMSF, como se muestran en las figuras 6-9. La Figura 6, muestra que el RMSD del complejo AChE-Taspina presenta un equilibrio estable a aproximadamente 0.7 - 0.9 nm, estos valores pueden deberse a factores como flexibilidad intrínseca de la proteína, o cambios conformacionales en esta. El backbone de la proteína con RMSD en 0.2 - 0.3 nm indica que la estructura es bastante estable a lo largo de la simulación, lo que nos indica que la proteína no ha sufrido grandes cambios conformacionales, al menos en su cadena principal. El ligando con un RMSD 0.1 - 0.2 nm indica que ha tenido una desviación mínima durante la simulación, lo que sugiere que el ligando ha mantenido una posición fija con respecto a la posición inicial, mostrando una fuerte interacción entre el ligando y la proteína.

Figura 6 Valores de RMSD del complejo, backbone y ligando. 

En la figura 7 , se determinó el radio de giro del sistema a lo largo de 15 ns, donde nos muestra que ambas curvas estabilizan desde el inicio y no exhiben grandes fluctuaciones, lo que nos indica que la proteína y el complejo han alcanzado un equilibrio sin presentar grandes cambios conformacionales, demostrando flexibilidad dinámica de la proteína y su interacción con el ligando. En la proteína, el radio de giro se encuentra entre 2.25 - 2.3 nm, lo que indica una estructura bastante compacta. Con respecto al complejo, se observa que presenta un valor ligeramente más alto, con un valor de 2.3 - 2.35 nm, lo que sugiere que el ligando aumenta el tamaño del sistema, pero sin una contribución significativa por un buen acoplamiento a la proteína.

Figura 7 Valores de Radio de Giro (Rg) del complejo y backbone de la proteína. 

La figura 8, nos muestra el análisis del Área Superficial Accesible al Solvente (SASA) que mide el área de la superficie de la molécula que está expuesta al solvente¹⁶. Las fluctuaciones que se observan son típicas en una simulación de dinámica molecular que corroboran que no hubo cambios drásticos en la exposición al solvente a lo largo del tiempo de simulación. Este gráfico denota que después de los 5ns tanto el complejo como el backbone alcanzan un equilibrio, con fluctuaciones menores, lo que indica que la proteína y el ligando han alcanzado un estado estable en términos de interacción con el solvente. Adicionalmente, observa que el SASA del complejo es un poco menor al backbone, esto sugiere que el ligando está bien acoplado a la proteína y posiblemente protegido por la estructura proteica, esto es un indicio favorable entre el ligando y proteína, denotando que el ligando no está demasiado expuesto al solvente.

Figura 8 Análisis SASA del complejo y backbone de la proteína. 

En la figura 9, se representa el análisis RMSF que nos permite evaluar la flexibilidad de cada residuo de la proteína a lo largo de la simulación. Ambas curvas muestran valores de RMSF bajos, lo que indica que la proteína ha permanecido bastante estable durante la simulación, denotando ser un buen indicador de estabilidad en la dinámica. Se observan dos picos (el primero y último) con elevados valores de RMSF que alcanzan los 0.5 - 0.6 nm. Siendo esto común en las proteínas, ya que los extremos suelen ser más flexibles y no están restringidos estructuralmente. Aunque las dos curvas están superpuestas, los valores de curva roja tienden a ser ligeramente más altos, lo que indica que la presencia del ligando puede estar contribuyendo a fluctuaciones adicionales en la estructura de la proteína. Esto una vez más nos muestra que la simulación de dinámica molecular ha sido bien equilibrada y que la proteína ha mantenido su estructura durante la simulación.

Figura 9 Análisis RMSF del complejo y backbone de la proteína. 

Estos análisis de la simulación de Dinámica Molecular revelan en conjunto una estabilidad estructural entre el AchE y taspina, lo que sugiere que podría ser un buen inhibidor potencial de la enzima, ya que podría estar bien acoplada al sitio de unión.