Especial Premio Nobel de Química 2013.
Las metodologías desarrolladas a partir de los trabajos de los tres premiados nos han dotado del equivalente a un «microscopio in silico» que nos permite descifrar las claves íntimas del funcionamiento de las macromoléculas biológicas. Hay pocas dudas que nuestra visión sobre los sistemas biomacromoleculares sería mucho más limitada que lo es en la actualidad sin los trabajos pioneros de Karplus, Levitt & Warshel.
La historia de la biología estructural computacional comenzó a finales de los años sesenta en el laboratorio de Shenior Lifson (un visionario científico fallecido en 2001) en el Weizmann Institute en Israel. En un momento en el que todos los químicos teóricos miraban hacia la mecánica cuántica, Lifson pensó que era posible una alternativa: emplear la mecánica clásica con potenciales ajustados para reproducir la realidad experimental o cálculos cuánticos precisos. Un ejemplo que suelo explicar en mi clase de «Estructura de Macroléculas» de la Universidad de Barcelona ilustra bien la idea básica de Lifson: consideremos una molécula de hidrógeno, si queremos representarla mediante mecánica cuántica simplemente definiremos dos núcleos de hidrógeno y dos electrones, describiremos las zonas por las que se pueden mover los segundos mediante funciones más o menos complejas y dejaremos que la resolución de la ecuación de Schroedinger nos encuentre la mejor distribución relativa de núcleos y electrones. Del cálculo convergido podremos derivar todas las propiedades de equilibrio de la molécula de hidrógeno, incluso determinar cómo reaccionaría al intentar alargar o acortar su longitud de enlace.
Está claro que la química cuántica nos permite conocer con todo detalle la molécula de hidrógeno. Ahora bien, ¿está justificado hacer un cálculo extremadamente complejo cada vez que queremos representar una molécula de hidrógeno en un entorno diferente? ¿Cómo abordaremos el estudio cuando la molécula de hidrógeno está encapsulada en un sistema con miles de átomos, imposible de tratar cuánticamente por su elevado coste computacional? ¿Por qué no aprovechar nuestro conocimiento experimental sobre la molécula de hidrógeno para reducir el coste de representarla teóricamente?
Lifson se formuló y contestó mentalmente a estas preguntas. Se dio cuenta que, especialmente en sistemas biológicos, las distorsiones de los términos de enlace son pequeñas y el entorno es poco agresivo, por lo que para muchos (ciertamente no para todos) los sistemas de interés podemos asumir que la molécula de hidrógeno se puede asimilar a dos partículas de van der Waals unidas por un muelle con una distancia de equilibrio igual a la experimentalmente conocida del enlace H-H, y con una constante de fuerza igual a la que justifica la frecuencia de vibración experimental de esta molécula. Simplificamos así la complejidad de la curva de disociación/fusión del hidrógeno, cuyo estudio requeriría un tratamiento cuántico complejo, asumiendo que nos movemos cerca de la posición de equilibrio, que los efectos cuánticos son moderados y que el potencial de Morse se puede ajustar a una simple parábola. Nótese, que en esta aproximación los electrones no se explicitan pero su efecto se introduce indirectamente en los «parámetros» (por ejemplo en el caso de la molécula de hidrógeno la distancia de equilibrio y la fuerza del enlace) que definen las ecuaciones de energía (el «campo de fuerzas, o force-field«). Nos movemos pues de un Hamiltoniano cuántico basado en primeros principios y dependiente de posiciones de núcleos i electrones a uno clásico, dependiente únicamente de las posiciones de los núcleos, mucho más eficiente computacionalmente, pero que precisa de input experimental en su definición.
A finales de los sesenta Lifson había delineado el primer «campo de fuerzas», el consistent force-field (CFF) que él pensaba podía servir para estudiar incluso las estructuras de proteínas y ácidos nucleicos que en esos momentos se empezaban a describir por cristalografía de rayos X. Fue en este periodo (en torno a 1967) cuando dos estudiantes de doctorado (Michael Levitt y Arieh Warshel) empezaron a codificar el force-field (Figura 1) en un programa FOTRAN que calculaba, no solo la energía de un sistema, sino sus derivadas primeras y segundas. En 1968 se publicó la primera minimización energética de una proteína mediante el uso de un campo de fuerza, algo que hubiera resultado impensable en ese momento (y posiblemente aún ahora) mediante metodología cuántica. Levitt y Warshel abandonaron Israel, el primero para continuar su tesis en Cambridge y el segundo para realizar un postdoc en el laboratorio del tercer laureado Martin Karplus. En aquel momento Karplus estaba estableciendo su grupo en Harvard, después de sus estudios de doctorado con Linus Pauling, su postdoc con Charles Coulson (Oxford), y estancias en Illinois y Columbia. Martin Karplus era una estrella emergente en el campo de la química teórica, uno de los primeros en ver el potencial que tenía estudiar los sistemas biológicos desde la perspectiva de la química física. La estancia de Warshel sirvió para que Karplus entrara en contacto con CFF y a partir de él, con la ayuda de Bruce Gelin, un programador experto proveniente del ejército, y Andy McCammon (un brillante estudiante de doctorado), desarrollaran un programa mucho más eficiente, que con el tiempo evolucionaría en CHARMM (Chemistry at HARvard Molecular Mechanics), un programa aún hoy ampliamente usado en todo el mundo, y que fue el origen de muchos de los códigos actuales (NAMD, GROMOS, GROMACS, AMBER,…) de dinámica y mecánica molecular.
A finales de los setenta se publicarían los dos artículos más influyentes en el campo, el primero (Warshel,A.; Nature, 1976, 260,679; citado en 389 publicaciones posteriores (hasta Marzo 2014)) en el que Warshel demostraba que mediante simulación molecular era posible entender la primera etapa de un fenómeno tan complejo como la visión, y el segundo (McCammon,J.A.; Gelin,B.R.; Karplus,M. Nature, 1977, 267, 679; citado en 2099 publicaciones posteriores) donde Karplus y colaboradores presentaban la primera simulación de dinámica molecular de una pequeña proteína. El impacto de estas dos publicaciones en el campo fue rápido y contundente, conduciendo a que muchos otros grupos en Estados Unidos y Europa entraran a trabajar en el desarrollo de esta tecnología emergente.
La idea general de un cálculo de dinámica no puede ser más sencilla (véase Figura 2), aunque su implementación práctica diste de serlo y los algoritmos computacionales que lo hacen sean unas obras complejas de ingeniería. El cálculo empieza con una conformación inicial del sistema, derivada de datos experimentales o de una simulación previa. Para esa estructura se calculan un conjunto de velocidades iniciales para cada átomo, dichas velocidades se seleccionan al azar pero en conjunto deben hacer cumplir la ecuación térmica de estado. La incorporación de la definición topológica del sistema (es decir definir el esquema de enlaces covalentes del sistema) y la aplicación del force-field a las coordenadas iniciales nos permiten obtener un primer estimado de la energía del sistema del cual por derivación para las tres coordenadas de cada átomo obtenemos la fuerza que actúa sobre ellos. Derivar aceleraciones a partir de fuerzas es trivial, pero integrar esas aceleraciones es mucho más complicado y requiere del uso de métodos numéricos que implican realizar integraciones en escalas temporales muy pequeñas (alrededor del femtosegundo), de las que se obtienen las nuevas velocidades y posiciones. En este punto se vuelven a calcular energías y fuerzas a fin de realizar otro paso de la trayectoria (Figura 2), repitiéndose el proceso hasta que se determinan el número deseado de etapas de integración. Remarcar una vez más el esquema dibujado aquí debe de considerarse únicamente a título ilustrativo, en tanto que los algoritmos reales permiten introducir diferentes etapas de integración, incorporar restricciones de presiones y/o temperatura, sesgar la trayectoria por fuerzas externas, o bien obligarla a adaptarse a descriptores experimentales y son por ello mucho más complejos.
Como se comentó en el párrafo anterior, desde su formulación inicial la dinámica molecular ha vivido una evolución constante, que no parece haber llegado a su fin. Así, durante los últimos años de la década de los setenta y todos los ochenta el campo vivió en ebullición formal, con aportaciones constantes que permitieron ampliar el rango de aplicación de la técnica de dinámica molecular e hicieron posible considerar sistemas cada vez más realistas. Los trabajos de Berendsen, Jorgensen, Van Gunsteren y muchos otros permitieron considerar las proteínas y ácidos nucleicos como sistemas hidratados en condiciones isotérmicas e isobáricas que se aproximaban a las fisiológicas. Grupos como los de Karplus, Kollman, Case, Brooks, van Gunsteren y muchos otros abrieron campos insospechados de aplicación de la dinámica molecular en el mundo de la biología durante la década de los ochenta, atrayendo por primera vez la atención de los grupos experimentales.
Los noventa alumbraron nuevos campos de fuerza, que manteniendo el formalismo básico del CFF de Lifson fueron calibrados con mucho más cuidado. Gracias a ello se consiguió poder predictivo y popularizar, por ejemplo, el uso de cálculos de dinámica molecular para refinar modelos estructurales de rayos X o RMN. Son estos campos de fuerza los que con refinados posteriores están todavía en uso, recogiendo en muchos casos millares de citas (como ejemplo solo el modelo TIP3P de agua desarrollado por el grupo de W. Jorgensen reunía casi 15000 citas a finales de marzo de 2014). Es una década en la que es imprescindible reconocer el trabajo ingente de grupos como los de van Gunsteren, Berendsen, Karplus, Scheraga, Kollman o Jorgensen (entre otros). Ellos fueron los que trasladaron las ideas iniciales en una metodología madura y accesible para toda la comunidad.
A partir de finales de los noventa quedó claro que aun disponiendo de campos de fuerza extremadamente precisos la dinámica molecular encontraría dos límites: i) relacionados con las aproximaciones implícitas a la mecánica clásica, ii) relacionados con la dificultad para representar correctamente el espacio de fases de sistemas tan complejos como las proteínas. A fin de mejorar la representación del espacio de fases surgieron tres aproximaciones diferentes, aunque complementarias. Por un lado, autores como McCammon, Jorgensen, Karplus, Kollman, Jarzynski, Schulten, van Gunsteren, Parrinello y muchos otros implementaron métodos para muestrear adecuadamente la parte más interesante del espacio de fases. Por otro, autores como Pande, Hummer, Best, Okamoto y otros crearon métodos que permiten definir un espacio de fase complejo, a partir de procesar datos de trayectorias individuales cortas, mucho más fáciles de obtener que las trayectorias largas. Finalmente, otros grupos optaron por focalizar su esfuerzo en optimizar los códigos de dinámica molecular, a fin de hacerlos más eficientes en el aprovechamiento de las nuevas tecnologías computacionales. Muy destacable en este campo el esfuerzo sostenido del grupo de K.Schulten (los desarrolladores del programa NAMD), y del consorcio GROMACS, liderado por E.Lindhal y D.Van der Spoel, que han proporcionado a la comunidad magníficas piezas de ingeniería de software detrás de las que se esconden centenares de hombres/años de trabajo y que permiten optimizar el uso de supercomputadores paralelos. Muy recientemente el esfuerzo de estos y otros grupos (por ejemplo el del consorcio AMBER liderado en este área por R.Walker) se ha focalizado en adaptar los programas a nuevos procesadores de uso específico como son las tarjetas gráficas (Graphical Processing Units, GPUs). Especialmente relevante ha sido en estos últimos años el trabajo de la empresa creada por un millonario americano (y antiguo catedrático de arquitectura de computadores), D.E. Shaw, que ha desarrollado una familia de ordenadores (llamados Anton en honor del inventor del microscopio) específicamente creados para realizar cálculos de dinámica molecular. Gracias a Anton Shaw Co. se ha conseguido realizar simulaciones que pasan la barrera mítica del milisegundo en longitud de trayectoria (pensemos que la integración de las ecuaciones de Newton se realiza cada femtosegundo, es decir que obtener 1 milisegundo de cálculo implica realizar 1012 cálculos de energías y fuerzas para sistemas con centenares de miles de grados de libertad), pudiendo así estudiar una gama enorme de procesos dinámicos, incluidos movimientos alostéricos, cambios conformacionales en receptores de membrana, unión flexible de ligandos, o plegamiento de proteínas.
Las limitaciones intrínsecas derivadas del uso de mecánica clásica fueron visibles ya en los inicios del desarrollo de la dinámica molecular. Karplus, Levitt, Warshel y muchos otros se percataron que la aproximación puramente clásica limitaba, a menudo en exceso, el rango de aplicación de estos métodos a sistemas bioquímicos. Excluía, por ejemplo, algo que a los tres les fascinaba: el estudio de la reactividad enzimática. Esto llevó a modificar el formalismo básico a fin de considerar el sistema dividido en dos partes, una pequeña (por ejemplo, el centro activo de un enzima), donde ocurren procesos cuya representación requiere un tratamiento cuántico del sistema (por ejemplo, una reacción química), y otra mucho más extensa, donde una descripción cuántica no es necesaria, ya que tratamos con interacciones débiles que se pueden representar razonablemente a partir de la mecánica clásica (Figura 3). Esta idea básica es la que ha dado lugar a los métodos híbridos QM/MM (QM por quantum mechanics y MM por molecular mechanics) que son los de elección para estudiar sistemas bioquímicos en los que se produzcan fenómenos de reactividad química.
A nivel práctico el uso de métodos QM/MM presenta algunos problemas que han sido resueltos de manera diferente por distintos grupos. Uno, sutil, pero tremendamente complejo emerge de la propia complejidad de representar el movimiento del electrón en la partición cuántica del sistema. En general, esto se ha solucionado trabajando con trayectorias en el régimen de Born-Openheimer, es decir, integrando el movimiento de los núcleos y asumiendo que para cada uno de los movimientos nucleares la nube electrónica ha alcanzado su estructura de equilibrio. Un segundo problema complejo es la definición de un acoplamiento eficiente entre las particiones cuántica y clásica. Aquí han pervivido dos líneas de pensamiento, la de Warshel que se inspira en el modelo “valence-bond” y la que sugirió originalmente Karplus y que fue seguida después por muchos otros autores entre ellos M. Field, J. Gao, o W. Thiel que utilizaba el modelo de orbitales moleculares acoplado a cálculos de campo autocoherente (MO-SCF). Según este modelo el acoplamiento clásico-cuántico se introduce mediante un operador de perturbación adicionado al Hamiltoniano molecular que describe la partición cuántica. Un último escollo, aun no resuelto totalmente, viene dado por el coste computacional intrínseco del cálculo cuántico de la energía y fuerzas asociadas. El coste dificulta, a menudo impide, que se puedan integrar trayectorias suficientemente largas en el régimen QM/MM empleando un nivel elevado de teoría en la partición cuántica.
Ya a mediados de los ochenta, como una consecuencia natural de la evolución del campo QM/MM empezaron a presentarse las primeras aproximaciones que intentaban obtener una representación cuántica de todo el sistema dentro de un cálculo de dinámica molecular. Es lo que se ha venido a denominar como cálculos de dinámica molecular ab initio. El formalismo más extendido es el de Carr-Parrinello, en el que se integran tanto los grados de libertad nucleares como electrónicos (estos de manera más aproximada) dentro del contexto de un cálculo de ondas planas y un formalismo del funcional de la densidad muy simple. Otras aproximaciones como el algoritmo del consorcio español SIESTA emplea orbitales moleculares en lugar de ondas planas y aceleradores para la diagonalización de la matriz de Fock. Todos estos métodos tienen aún un rango de aplicación moderado en bioquímica, por el coste computacional que implican. Recientemente, se han implementado por Gao y colaboradores aproximaciones de “dividir-y-conquistar” para partir un cálculo cuántico (en el régimen Born-Openheimer) en un sistema muy grande en varios cálculos más pequeños que se realizan de forma casi-independiente, mejorando con ello la eficiencia del cálculo. Con esta aproximación y el uso de niveles bajos de teoría (Hamiltonianos semiempíricos) se ha conseguido simular trayectorias en la escala del sub-nanosegundo de proteínas hidratadas. Es aún pronto para determinar el campo de aplicación para la dinámica ab initio, pero posiblemente sea relativamente limitado en bioquímica, donde el tamaño y falta de simetría del sistema hacen estos cálculos prohibitivos.
Vista la historia, ¿cuál podemos esperar que sea su futuro?
Podemos esperar mejoras continuas en los force-fields, que se adaptarán para incorporar términos energéticos (como la polarización o la transferencia de carga) hoy ignorados. Es difícil, no obstante, que en términos generales se llegue a una descripción puramente cuántica de los sistemas bioquímicos. Se esperan grandes avances en la formulación de force-fields de baja granularidad, que a expensas de una reducción en el nivel de resolución, permitirán estudiar sistemas de tamaño muy grande aproximándonos al paradigma de la simulación molecular a nivel celular, o al menos a nivel de una ruta concreta.
Es previsible que los códigos actuales continúen su mejora de manera sostenida en el tiempo, adaptándose a los nuevos ordenadores. Es, no obstante, poco realista pensar que se consiga, en general, mejorar sustancialmente su paralelismo. Sí que podemos pensar que la iniciativa de D.E. Shaw se popularice y aparezcan otros ordenadores específicos para realizar cálculos de dinámica molecular, lo que permitirá que el acceso a la escala del milisegundo se extienda y que el uso de la técnica se convierta en rutinario en proyectos de diseño racional de fármacos. No obstante creo que el avance más significativo en este campo vendrá por la popularización del uso de técnicas de reconstrucción de espacios de muestreo a partir de simulaciones cortas, muy fáciles de obtener en el contexto del “cloud computing”.
Aventuro, no obstante, que el salto más importante en la simulación biomolecular, y en concreto en la dinámica molecular, vendrá de la integración efectiva de datos experimentales en la simulación. Creo que en una década los grupos teóricos de simulación molecular irán incorporando expertise experimental, hasta acabar derivando en grupos híbridos, capaces de combinar sin discontinuidades resultados teóricos y experimentales. Paralelamente, los grupos hoy puramente experimentales se familiarizarán con el uso de información dinámica experimental y las bases de datos estructurales, como el protein data bank (PDB), se expandirán para introducir parámetros de flexibilidad derivados de cálculos teóricos.
Podemos visualizar un escenario futuro donde la simulación molecular, iniciada por Karplus-Levitt-Warshel se combine con simulaciones más globales, propias de la biología de sistemas, e incluso con técnicas mecánicas de simulación de órganos, para dotarnos con una capacidad global de simulación de los sistemas vivos complejos, aproximándonos al sueño de crear avatares in silico de los seres vivos.
Referencias:
- Warshel,A.; Nature, 1976, 260,679
- McCammon,J.A.; Gelin,B.R.; Karplus,M. Nature, 1977, 267, 679