Análisis de componentes principales: una herramienta para dilucidar los movimientos esenciales de las biomoléculas

El Análisis de Componentes Principales (PCA, por sus siglas en inglés) es un procedimiento que se utiliza para reducir el volumen de datos utilizado para describir el estado de un sistema, tratando de minimizar la información que se pierde al hacer esta reducción [1, 2]. El algoritmo fue propuesto a principios del siglo XX [3, 4], pero inicialmente tuvo escasas aplicaciones porque resulta muy tedioso realizar “a mano” los cálculos requeridos, especialmente si el sistema considerado es grande. Pero justamente, cuando el sistema es grande, es cuando puede ser útil o necesario aplicar PCA. Posteriormente, el advenimiento de las computadoras modernas permitió que PCA sea aplicado con relativa facilidad sobre sistemas verdaderamente grandes y complejos. En la actualidad, las implementaciones de PCA son increíblemente diversas y aparecen en casi todas las ramas de la ciencia y la tecnología. Por ejemplo, PCA se utiliza en algoritmos de reconocimiento facial [5, 6], en análisis genéticos [7] y en la compresión de imágenes [8]; pero también en las ciencias médicas [9] y ambientales [10, 11], e incluso en las ciencias sociales, para revelar correlaciones entre las diferentes variables que determinan los estados de un sistema.¹

Cuando se realiza una simulación de dinámica molecular de una biomolécula, se obtienen archivos enormes que contienen las coordenadas cartesianas (x, y, z) de los átomos del sistema ,² en cada instante en el que se tomaron muestras. El problema que se plantea entonces es como transformar ese volumen descomunal de numeritos en datos inteligibles. Y en el caso particular de simulaciones de biomoléculas, estamos interesados en determinar cuáles son sus movimientos más significativos, asumiendo que los mismos son requeridos para el cumplimiento de su función biológica. Justamente en este punto, en donde la utilización de PCA resulta fundamental.

La primera aplicación de PCA para caracterizar los movimientos de las proteínas fue realizada por García en 1992 [14], quien analizó una simulación de 240 ps de Crambina y observó que los movimientos de partes distantes la proteína estaba altamente correlacionados. Esta correlación demuestra que es posible reducir la cantidad de datos requeridos para describir las fluctuaciones del sistema.³ Poco después, Berendsen y sus colaboradores propusieron utilizar PCA para transformar las coordenadas cartesianas directamente obtenidas de las simulaciones, en nuevas coordenadas con una característica muy especial: solo un puñado de ellas describe la mayor parte de las deformaciones moleculares [15].

En esta sección presentaremos, en primer lugar, el procedimiento seguido para un análisis de Componentes Principales. A continuación, indicaremos lo que se consigue con cada uno de esos pasos (o sea, indicaremos cual es el resultado o la acción efectivamente realizada en cada uno). Por último, terminaremos esta sección discutiendo cuáles son las motivaciones que usualmente nos llevan a hacer PCA de simulaciones MD de biomoléculas.

Para presentar el procedimiento seguido en un PCA, consideraremos que el análisis se realiza sobre las coordenadas colectadas de una simulación de dinámica molecular (MD-PCA), dado que esta es la situación mas común. Por ejemplo, en simulaciones de proteínas, típicamente se emplean las coordenadas cartesianas de los carbonos alfa, o de los átomos del esqueleto proteico. Sin embargo, otras opciones son posibles y muchas veces son preferibles. Estas alternativas serán discutidas en la Sección 5.2. Asimismo, debemos señalar que un PCA también se puede realizar sobre coordenadas derivadas de otras fuentes. Por ejemplo, se puede hacer PCA de modelos moleculares generados experimentalmente, tales como la resonancia magnética nuclear (RMN) o la cristalografía de rayos X. En cualquier caso, los resultados de estos usos alternativos pueden comprenderse sin inconvenientes con base en la discusión proporcionada aquí para el MD-PCA.

Utilizaremos la letra N para indicar el número de estructuras de la macromolécula que se muestrearon de la simulación, mientras que denota a un vector que contiene las coordenadas cartesianas atómicas de la k-ésima; estructura. El n-ésimo elemento de será denotado como . Por lo tanto, el índice k va de 1 a N mientras que n va de 1 a 3 , donde es el número de átomos incluidos en el análisis. Usualmente no incluye a todos los átomos de la macromolécula sino a un conjunto reducido. Por ejemplo, en los estudios de proteínas se suelen incluir solo los Cα o los átomos del esqueleto proteico. Discutiremos este tema con mas detalle en la Sec. 5.2. Con base en las definiciones recién elaboradas, los pasos para realizar el PCA de una macromolécula son los siguientes.

1. Eliminar del conjunto de coordenadas colectadas , el efecto de la traslación y la rotación global de la macromolécula. Para ello, todas las estructuras del conjunto se alinean con una de ellas, típicamente la primera.⁴ Si nos “salteamos” este paso (o si lo realizamos en forma deficiente) la animación del primer autovector de PCA muestra a la molécula rotando y/o trasladándose. En ese caso, los vectores que describen las deformaciones internas de la molécula (aquellas que efectivamente queremos identificar con PCA), quedan relegados a posiciones de menor importancia, lo que puede arruinar todo el análisis. Los fundamentos de los métodos de alineamiento, junto con sus detalles técnicos se discuten en la Sec. 5.1.

2. Calcular la matriz de covarianza C del conjunto de coordenadas . Los elementos de C están dados por,

donde y son los valores medios de las coordenadas y , respectivamente. Por ejemplo,

3. Diagonalizar la matriz C. Esto significa que debemos encontrar la matriz R que transforma a C en la matriz diagonal Λ mediante,

donde es la transpuesta de la matriz R. Dado que C es una matriz simétrica, la matriz R siempre existe, es ortogonal y puede ser transformada en ortonormal. De hecho, los programas que calculan PCA de biomoléculas siempre producen matrices R que están normalizadas. Por lo tanto, , o lo que es equivalente,

indicando que las columnas de R forman una base de vectores ortonormales (o versores). Asimismo, la gran mayoría de los programas ordenan los autovectores siguiendo un orden decreciente de sus autovalores.

Es conveniente notar que la Ec. 3 es equivalente a resolver la ecuación de autovalores,

para n= 1 . . . N , donde el vector es la n-ésima columna de R, mientras que es el n-ésimo elemento diagonal de Λ. En principio, el número de autovectores con autovalores distintos de cero es , porque este es el número de grados de libertad internos del sistema. ⁵Sin embargo, pueden aparecer más autovalores nulos si el tamaño de la muestra es insuficiente. Esto ocurre porque el número de direcciones requeridas para dar cuenta de las fluctuaciones de un conjunto de N puntos es . Entonces, si realizamos PCA sobre un conjunto de N estructuras, solo obtendremos autovectores con autovalores no nulos. Aunque el sistema pueda tener otras deformaciones, estas no van a aparecer en nuestro análisis porque no le dimos al algoritmo información suficiente. En general, esta limitación no genera mayores inconvenientes porque solo un puñado de autovectores se utiliza en el análisis posterior, pero es un punto a tener en cuenta cuando se dispone de un número reducido de estructuras.

4. Definir el espacio esencial (ES) de la molécula. Para comprender lo que esto significa es necesario considerar el espectro de autovalores típico del PCA de una molécula biológica. La Fig. 1 muestra un ejemplo de estos espectros. La línea violeta de este grafico muestra los autovalores individuales en función de su índice. Se puede apreciar que la curva tiene una caída inicial muy abrupta, de manera tal que para índices mayores a 10, los autovalores son casi despreciables. Esta es una característica muy notable si consideramos que el sistema tiene un total de 354 autovalores no-nulos. El ES de la molécula se define como el conjunto formado por los primeros autovectores, siendo . En la Sec. 5.3 discutiremos los criterios utilizados para decidir el valor de .

5. Proyectar los vectores desplazamiento, , sobre los autovectores del ES. De esta manera se obtienen los Componentes Principales (PC) para cada tiempo k en el cual se tomó una muestra,

Los PCs son las nuevas coordenadas que describen la dinámica de la macromolécula en el ES. Usualmente, todo el análisis subsiguiente se realiza sobre estas coordenadas y el resto se descarta. En cierto sentido, este procedimiento implica perder parte de la información contenida en el conjunto de coordenadas . No obstante, esta acción se justifica por considerar que los movimientos que quedan afuera del ES solo aportan ruido térmico y por lo tanto complican en análisis de los movimientos relevantes.

La sección anterior describe como realizar PCA sobre la trayectoria de una macromolécula, pero no explica nada acerca de cuáles son los resultados de esos procedimientos. En otras palabras, ¿qué interpretación podemos dar a los autovectores, autovalores y Componentes Principales obtenidos con el procedimiento descrito? Esta sección está dedicada a discutir esta cuestión, a fin de que se comprendan las motivaciones para realizar PCA y se puedan estimar sus limitaciones.

Una forma conveniente de iniciar la discusión es señalar que la conformación de una biomolécula puede ser proporcionada de maneras alternativas. Una de ellas, quizás la más común, es usar el vector , cuyas componentes miden los desplazamientos cartesianos de los ´atomos desde su posición promedio

Alternativamente, podríamos emplear el vector , cuyos elementos son los Componentes Principales. Estas dos representaciones están relacionadas por cambio en la base utilizada para expresar al vector que señala la configuración del sistema dentro de s u e norme espacio de configuraciones. Veamos, con un ejemplo sencillo, qué queremos decir con esto.

En las secciones precedentes hemos presentado el procedimiento para hacer PCA. Además, explicamos lo que se logra con cada uno de sus pasos discutimos cuáles son las motivaciones para llevar a cabo dicho procedimiento. Muñidos de esta información, estamos ahora en condiciones de ver algunos “truquitos” que permiten sacar un mayor provecho de la aplicación de este algoritmo.

Hasta este punto, hemos considerado que el PCA se realiza sobre estructuras muestreadas de una única simulación MD. Sin embargo, una práctica muy extendida consiste en utilizar una trayectoria artificial obtenida al combinar simulaciones que emplean diferentes estados iniciales del mismo sistema. A este procedimiento también se lo llama “concatenación” de trayectorias porque los archivos de las mismas se ponen uno a continuación del otro. Las trayectorias combinadas podrían corresponder, por ejemplo, a las formas abierta y cerrada de un receptor, o a una enzima unida a diferentes sustratos, o una misma proteína en diferentes solventes. Este procedimiento fue propuesto en 1995 por Berendsen y colaboradores, quienes analizaron las diferencias entre las simulaciones de Termolisina en vacío y en agua. Desde entonces, el procedimiento se ha implementado en muchas situaciones. En la bibliografía, se conoce como Comb-ED por Combined Essential Dynamics, pero también aparece como PCA de trayectorias concatenadas. En su artículo original, Berendsen y colaboradores explicaron cómo interpretar los resultados del PCA de trayectorias concatenadas utilizando argumentos intuitivos. Mucho más recientemente, los autores de este artículo encontramos las expresiones analíticas que relacionan la matriz de covarianza de una trayectoria combinada con aquellas de las trayectorias individuales. Estas expresiones, proporcionan una base más sólida para comprender la información que proveen los componentes principales obtenidos por este procedimiento.

La matriz de covarianza de una trayectoria que combina n simulaciones
individuales de la misma molécula, , puede expresarse en funcion de las matrices de covarianza individuales, , como [37],

Aquí, es la matriz de covarianza formada por las estructuras promedio de las simulaciones individuales. Sus elementos están dados por,

donde es el promedio de la coordenada i en la simulación k, mientras que es el promedio de la coordenada i en la trayectoria combinada.

En palabras, la Ec. 9 nos dice que la matriz de covarianza de la trayectoria combinada se obtiene al sumar el promedio de las matrices de covarianza de las trayectorias individuales y la matriz de covarianza de sus estructuras promedio. Debemos observar que la matriz tiene solo n 1 autovectores con autovalores no nulos. Estos abarcan un espacio vectorial que contiene a las estructuras promedio de las simulaciones individuales.

Una forma sencilla de discutir qué se puede esperar de un PCA combinado es considerar un ejemplo particular. El más simple es el de concatenar dos simulaciones. De acuerdo con la Ec. 9, los elementos de la matriz de covarianza combinada para este caso particular están dados por:,

donde y son los elementos de las matrices de covarianza calculadas con las simulaciones 1 y 2, respectivamente, mientras que es la matriz de covarianza de sus estructuras promedio,

donde .

Como se discute en la Sección 3.2, los principales autovectores de las matrices y indican las direcciones de mayor fluctuación observada en cada una de las respectivas trayectorias mientras que sus autovalores miden la amplitud de esas fluctuaciones. La matriz , por otro lado, se calcula a partir de solo dos estructuras. Por lo tanto, hay una sola dirección de deformación: la que nos conduce de una estructura promedio a la otra pasando por el promedio global. Esta es la dirección del único autovector de con autovalor distinto de cero. Además, como se demostró en la Ref. [37], este autovalor está relacionado con la desviación cuadrática media entre las dos estructuras promedio de acuerdo con,

Teniendo en mente las consideraciones previas, se pueden prever dos casos límite para la diagonalización de la matriz . Una posibilidad es que las dos simulaciones independientes correspondan a estructuras relativamente rígidas que muestrean regiones muy dispares del espacio configuracional. En ese caso, la distancia entre sus estructuras promedio es mucho mayor que las fluctuaciones o observadas en las trayectorias individuales. En consecuencia, la matriz hace la contribución dominante a en la Ec. 11. Si esta contribución supera ampliamente la de las matrices de covarianza individuales, el primer autovector de es casi paralelo al de . Eso implica que está prácticamente alineado con la recta que pasa por las dos estructuras promedio.

Un ejemplo de esta situación se presenta en la Fig. 3, que utiliza datos tomados de la Ref. 37 para mostrar la proyección de dos trayectorias de la albumina sérica humana sobre el primer vector obtenido al hacer un PCA combinado con ambas trayectorias. Una de las simulaciones fue realizada con la enzima unida al ácido láurico mientras que en la otra utilizó la forma apo. En la figura se puede apreciar que la distancia entre las estructuras promedio es mayor que las fluctuaciones de cada trayectoria individual. Para este ejemplo, el primer autovalor de la matriz es 4499.40 , mientras que el único autovalor de la matriz es 4361.59 . El producto escalar entre los dos autovectores se vale 0.999, demostrando que son prácticamente el mismo vector. Para hacer este PCA se utilizaron los 585 átomos de Cα de la enzima y si se reemplazan estos datos en la Ec. 13, se encuentra que el RMSD entre las estructuras promedio es de 5.46 , lo que coincide perfectamente con el obtenido por medios independientes. Obviamente, no hay ninguna necesidad de realizar todo el protocolo de PCA combinado para calcular el RMSD entre dos estructuras promedio, cuando todos los programas de visualización molecular permiten hacer el mismo calculo con solo dos “clicks” en el botón del “mouse”. Respecto al significado del resto de los autovectores de , no es posible hacer predicciones generales ya que los resultados dependen de los pesos relativos de las matrices y , que pueden cambiar de un caso a otro.

Uno de los principales objetivos de PCA es identificar los cambios conformacionales más prominentes de la biomolécula en estudio. Estas conformaciones se modifican cuando varían las distancias entre los átomos de la molécula, pero no cambian cuando la misma se traslada o rota como un cuerpo rígido. Sin embargo, las simulaciones MD que se realizan en coordenadas atómicas cartesianas mezclan estos tipos de movimiento. Por lo tanto, antes de realizar un PCA usando estas coordenadas, es necesario eliminar la influencia de la traslación y rotación global, tal como se indicó en la Sec. 3.1.

El procedimiento se basa en la siguiente idea. Si adosamos un sistema de ejes cartesianos a la molécula, la traslación y la rotación global de la misma pueden ser descritos por los movimientos de dicho sistema de ejes. De ahora en más, nos referiremos al mismo como el sistema de ejes FM por “Fijo en la Molécula”, para distinguirlo de un sistema de ejes que no se traslada ni rota, al que llamaremos FE por “Fijo en el Espacio”. De este modo, la traslación global de la molécula se describe por tres coordenadas cartesianas que señalan el origen del sistema FM desde el origen del sistema FE, mientras que las rotaciones se describen por tres ángulos de Euler que proveen suorientación [34]. La clave del procedimiento reside en definir qué prescripción o receta vamos a utilizar para decir cómo hay que adosar el sistema FM a cada estructura molecular. Por ahora supongamos que hemos resuelto este punto. Entonces, dado un conjunto de estructuras moleculares muestreadas de una trayectoria MD, debemos adosar el sistema FM a cada una de ellas utilizando siempre la misma prescripción. Finalmente, eliminamos la rotación y la traslación global de la molécula alineando el sistema de ejes FM de cada estructura del conjunto con el de una estructura que hayamos elegido como referencia. La trayectoria que se genera al hacer este alineamiento es la que habríamos obtenido si el sistema de ejes FM se hubiera quedado quieto. Así, con buena aproximación, las coordenadas cartesianas resultantes solo dan cuenta de las deformaciones internas de la molécula. Sin embargo, el algoritmo no es exacto y la calidad del resultado depende de los criterios seguidos para adosar el sistema FM a cada estructura molecular muestreada. Veamos entonces, cuáles son los criterios que se suelen aplicar para ello.

El origen del sistema de ejes fijo a la molécula se elige siempre como el centro de masa de la molécula. Con esta elección, la energía cinética de la traslación global se separa de la energía de las rotaciones y de las vibraciones [34]. Sin embargo, no hay una manera única de seleccionar la orientación de los ejes, sino que existen diversas opciones y ninguna de ellas proporciona una separación exacta entre vibración y rotación. Para moléculas pequeñas y/o semirrígidas, el criterio más empleado es la segunda condición de Eckart. Esta alternativa asegura que las pequeñas deformaciones alrededor de la posición de equilibrio de la molécula no contribuyen al momento angular total. Además, esta elección minimiza el acoplamiento entre la energía rotacional y la vibracional [40].

Desafortunadamente, las macromoléculas en solución a temperatura ambiente no son cuerpos semirrígidos, sino que sus superficies de energía libre tienen muchos mínimos locales. Por lo tanto, la segunda condición de Eckart, que asume que hay una única estructura de equilibrio, no puede aplicarse. Una alternativa sería indicar que el sistema FM coincide con los ejes de inercia de la molécula. Sin embargo, en el análisis de simulaciones de MD, el enfoque más comúnmente utilizado es el ajuste de cuadrados mínimos rotacionales presentado por McLachlan [41]. De acuerdo con el mismo, se busca superponer cada una de las estructuras del conjunto con la de referencia, mediante una transformación combinada. Primero, se aplica una traslación que hace coincidir sus centros de masa. Luego, se hace una rotación que minimiza la desviación cuadrática media entre los átomos de ambas estructuras. Se ha demostrado que el ajuste de cuadrados mínimos proporciona la misma alineación que las condiciones de Eckart, si en ambos casos se emplea la misma estructura de referencia [42]. Sin embargo, como se mencionó anteriormente, las macromoléculas biológicas a temperatura ambiente no tienen un único mínimo que pueda seleccionarse inequívocamente como referencia. Para empeorar las cosas, el procedimiento de ajuste de cuadrados mínimos produce resultados ligeramente diferentes si se emplean distintas referencias. Así, la separación aproximada de la rotación funciona bien para moléculas relativamente estructuradas, pero empeora cuanto más flexible sea la molécula en estudio.

5.2 ¿Qué coordenadas me conviene utilizar?

La elección de las coordenadas que conviene utilizar para hacer PCA está estrechamente ligada al propósito del análisis. De lo discutido en la sección anterior se desprende que un PCA realizado en coordenadas cartesianas no es apropiado para ningún trabajo cuantitativo tal como la identificación de las conformaciones estables y metaestables de una biomolécula o la evaluación de sus pesos relativos. Sin embargo, para obtener una representación pictórica de los movimientos relevantes de la molécula se utilizan coordenadas cartesianas. Dicha representación podría ser una imagen que muestre los desplazamientos producidos por un determinado autovector, como el que se muestra en la Fig. 5, o podría ser una película con su animación, como la Película III presentada como Información Suplementaria. Por otra parte, debido al estrecho paralelismo entre PCA y NMA, la mayoría de los programas utilizados para representar o animar modos normales también pueden ser aplicados para ilustrar Componentes Principales. La extensión Normal Mode Wizard del programa VMD es probablemente la más ampliamente empleada para ese propósito [43].

Las dificultades causadas por el uso de coordenadas cartesianas en el PCA de moléculas flexibles han sido estudiadas en detalle en la Ref. 44. En ese estudio determinaron, de manera unívoca, que dichas dificultades surgen de la separación inapropiada entre las deformaciones internas de la molécula y su rotación global, lo que introduce ruido numérico en el análisis. En consecuencia, las superficies de energía libre calculadas a partir de PCA realizadas en coordenadas cartesianas tienen características artificiales que impiden la identificación de los estados estables y metaestables del sistema. Estos problemas se han encontrado incluso en el PCA de una proteína tan rígida como BPTI, cuando se utilizó una trayectoria muy larga [44]. Para superar estos inconvenientes se han ideado esquemas de superposición sofisticados [45] que mejoran el análisis en algunos casos, pero no en todos [44]. Alternativamente, todas las dificultades se pueden evitar si se lleva a cabo el PCA en coordenadas internas, tales como distancias y/o ángulos diedros. Como discutiremos a continuación, cada una de estas opciones se adapta a un tipo particular de problema y tiene sus propias limitaciones y deficiencias.

Para un conjunto de N puntos, hay distancias entre pares, mientras que el número de grados de libertad internos, o sea los grados de libertad que describen las deformaciones de la nube de puntos, es . Por lo tanto, para , hay más distancias que coordenadas internas independientes. Además, el número de distancias crece muy rápidamente con N y pronto se convierte en un número enorme. ⁹ Esto hace que la diagonalización de la matriz de covarianza, requerida para hacer PCA, se vuelva computacionalmente muy costosa o incluso sea imposible de realizar. Por lo expuesto, antes de hacer un PCA basado en distancias, uno debe elegir cuidadosamente cuáles se van a incluir en el análisis. Las opciones típicas son utilizar las distancias entre carbonos alfa o las distancias entre los átomos más próximos de cada par de residuos. Sin embargo, incluso este conjunto suele ser demasiado grande para una implementación numérica eficiente [46], por lo que se suelen agregar otros criterios para reducirlo aún más. Para tal fin, el criterio más utilizado es emplear solo las distancias entre residuos que están en contacto en la estructura nativa [47]. Además, las distancias entre residuos que están separados por menos de cuatro residuos en la secuencia proteica también se descartan.

Por su parte, para el PCA basado en diedros, se pueden utilizar los ángulos del esqueleto de la molécula o también los de las cadenas laterales. La primera opción es, por mucho, la más utilizada, porque permite revelar los cambios conformacionales más significativos de la molécula. En cambio, si se hace un PCA de diedros de las cadenas laterales, lo que se observa son variaciones de los contactos entre residuos, pero el PCA basado en distancias de contacto es claramente una herramienta más adecuada para este fin [47].

Antes de realizar un PCA basado en ángulos diedros, uno debe tener en cuenta los problemas que pueden originarse de la periodicidad de los mismos. Debido a esta periodicidad, por ejemplo, puntos muy cercanos como 179° y -179° están separados artificialmente cuando el rango angular utilizado es (-180°:180°]. Esta separación artificial persiste, aunque se cambie el rango utilizado, porque dicho cambio solo modifica la región en la que se encontraran las dificultades. Por ejemplo, al usar [0°:360°), no hay problemas alrededor de 180° porque los mismos se han desplazado a alrededor de 0°/360°. A veces ocurre que los ángulos utilizados no muestrean todo el intervalo de 360°.¹⁰ En tales circunstancias, uno puede elegir el intervalo angular apropiado para cada caso, de manera tal que las discontinuidades aparezcan en las regiones que no son muestreadas por el sistema. En esta idea se basa el método llamado de dPCA+, desarrollado por Gerard Stock y colaboradores [48]. Una ilustración de cómo se realiza el cambio en el intervalo angular se proporciona en la Fig. 3 de la Ref. 48. Sin embargo, la estrategia de dPCA+ es engorrosa cuando uno debe lidiar con muchos ángulos diferentes, ya que requiere determinar el rango apropiado para cada uno de ellos. Además, para algunas moléculas y algunos ángulos en particular, puede suceder que todo el intervalo angular sea efectivamente visitado. Esto hace que dPCA+ sea inútil.

El procedimiento más ampliamente empleado para realizar PCA de ángulos diedros se llama dPCA y fue desarrollado por Stock varios años antes que dPCA+ [49]. El algoritmo propone duplicar el número de variables mediante la transformación,

donde n toma valores desde 1 hasta el número total de ángulos considerados. Esta alternativa tiene la clara desventaja de aumentar artificialmente el tamaño de las matrices con las que uno tiene que lidiar, pero también tiene diversas ventajas que explican su “popularidad”. Principalmente, con dPCA se evitan las discontinuidades que surgen de la periodicidad de los ángulos y se elimina el problema de los promedios mal definidos. Este último inconveniente afecta gravemente el cálculo de la matriz de covarianza. Para entender su origen, consideraremos dos ángulos, uno que oscila alrededor de 0°, entre -10° y 10°; y el otro que oscila alrededor de 180°, entre -170° y 170°. Dado que el promedio se sitúa en 0° en ambos casos, las desviaciones respecto del promedio, que se utilizan para calcular la matriz de covarianza, resultan ser mucho mayores en el segundo que en el primer caso. Sin embargo, esta diferencia no refleja la verdadera situación física, ya que en los dos casos los desplazamientos son de 10° respecto al promedio. Afortunadamente, el método dPCA elimina naturalmente todos estos inconvenientes. Además, provee un mapeo biunívoco de la distribución angular original. Esto significa que no crea ni suprime mínimos en las superficies de energía libre de la molécula. Finalmente, el procedimiento puede ser fácilmente sistematizado y constituye una de las opciones disponibles para realizar PCA en el marco de la gran mayoría de los programas que analizan simulaciones MD.

La experiencia con el análisis de diferentes sistemas indica que dPCA es el método de elección cuando se pretenden identificar los estados estables y metaestables de moléculas que experimentan grandes cambios conformacionales. Esto se aplica, por ejemplo, a estudios de plegamiento de proteínas, así como a cambios conformacionales de ARNs [50]. Por otro lado, si uno tiene como objetivo caracterizar los estados de proteínas estructuradas que llevan a cabo sus movimientos funcionales, el PCA basado en distancias resulta ser el más conveniente [47]. En cualquier caso, siempre se deben probar las diferentes alternativas, comparar sus resultados e intentar evaluar la calidad de sus predicciones mediante procedimientos independientes. Para concluir, señalamos que la Ref. 51 ofrece una discusión muy clara y detallada de las ventajas y desventajas de las diferentes coordenadas utilizadas para llevar a cabo PCA de macromoléculas. Recomendamos encarecidamente la lectura de ese artículo a aquellos que buscan profundizar en los temas presentados en esta sección.

5.3 ¿Cómo se determina el tamaño del espacio esencial?

Hasta este punto, hemos estado hablando acerca del ES, una pequeña subespecie vectorial que contiene las direcciones de las fluctuaciones más importantes de la molécula bajo estudio. Sin embargo, no hemos indicado aun ningún criterio que nos permita decidir cuántos autovectores deben incluirse en dicho espacio. De hecho, se pueden aplicar diversos criterios y la elección del mismo depende del propósito que se persigue al hacer el PCA.

Como se mencionó anteriormente, los autovalores obtenidos por PCA son las fluctuaciones cuadráticas medias de los Componentes Principales. Asimismo, al realizar el PCA en coordenadas cartesianas, la suma de todas estas fluctuaciones equivale a la suma de las fluctuaciones cuadráticas de los átomos incluidos en el análisis. Así, uno de los criterios más simples para decidir el tamaño del ES consiste en incluir los primeros n autovectores de PCA, de modo que estos den cuenta de una fracción dada de las fluctuaciones cuadráticas totales. Por lo tanto, si f es la fracción de las fluctuaciones que se desea tener en cuenta, y el PCA se realizó utilizando las coordenadas cartesianas de átomos, el número n de autovectores en el ES se elige de manera que,

Así, por ejemplo, para los datos mostrados en la Fig. 1, elegir f = 0,9 conduce a un ES con 12 autovectores mientras que establecer f = 0,75 incluye solo 5 autovectores.

En este artículo hemos discutido cuál es la motivación para realizar un PCA con los datos recopilados de una simulación de MD. Además, hemos indicado cómo interpretar los resultados. Se han presentado las principales deficiencias del MD-PCA y se han proporcionado varios consejos para mitigar el efecto de las mismas. Para finalizar, presentamos a continuación una lista de consejos que recomendamos seguir para lograr un uso provechoso de esta técnica.

* El primer paso para hacer un MD-PCA es seleccionar qué partes de la molécula se van a incluir en el análisis. Todas aquellas regiones que realizan grandes movimientos de carácter arbitrario deben descartarse. Los N- y C- terminales de las proteínas son ejemplos típicos de esto.

* El PCA de coordenadas cartesianas debe usarse principalmente con fines ilustrativos, tales como observar visualmente las deformaciones producidas por los principales autovectores.

* Al estudiar moléculas que experimentan cambios conformacionales significativos, tales como ARN, proteínas intrínsecamente desordenadas, o para simular procesos de plegamiento de proteínas, cualquiera de las variantes de PCA de ángulos diedros del esqueleto molecular representa la mejor opción.