Dinâmica Molecular

Vista previa del material en texto

Caṕıtulo 17
Dinámica molecular
Juan A. Bueren-Calabuig
17.1. Introdución
Para poder comprender con detalle la estructura, la función y los mecanismos de reacción de sistemas
biológicos y qúımicos es necesario emplear métodos teóricos de modelado molecular y de simulación
que complementen los resultados obtenidos mediante técnicas experimentales. A pesar de que la cris-
talograf́ıa de rayos X permite estudiar las propiedades f́ısicas y estructurales de las moléculas, la mera
visualización de una única estructura molecular no es capaz de resolver problemas más complejos como,
por ejemplo, el movimiento de una protéına o su unión a un ligando. En estos casos es necesario acudir
a la simulación de tales procesos. El principal objetivo del modelado molecular aplicado en distintas
áreas como la bioqúımica o la farmacoloǵıa, consiste en describir las interacciones biomoleculares en
base a las leyes generales de la qúımica y de la f́ısica [15, 19].
Hoy en d́ıa, el modelado molecular está ı́ntimamente ligado al uso de computadoras y gráficos inter-
activos. El espectacular desarrollo de los ordenadores y de las técnicas computacionales ha permitido
el uso de modelos teóricos para el estudio de todos aquellos problemas relacionados con la geometŕıa
y la enerǵıa de las moléculas. De hecho, se podŕıa llegar a modelar con exactitud cualquier sistema
biológico con suficiente poder de cálculo y un alto nivel de teoŕıa. Sin embargo, en la práctica, esto no
es siempre posible debido a limitaciones tanto de tiempo como de recursos materiales. Por este motivo,
es necesario introducir aproximaciones para mantener la viabilidad del experimento aplicando distintos
niveles de teoŕıa en función del sistema que va a ser analizado [29]. Como norma general, a medida que
aumenta el nivel teórico aumenta la calidad de los resultados obtenidos, pero también se incrementa
notablemente el coste computacional.
La descripción más rigurosa de un sistema viene definida por la mecánica cuántica (Quantum Me-
chanics, QM) que tiene en cuenta expĺıcitamente los electrones en sus cálculos, haciendo posible el
estudio de la estructura, las propiedades que dependen de la distribución electrónica y la reactividad
qúımica (formación y ruptura de enlaces). Sin embargo, su aplicabilidad está restringida a sistemas
con centenares de átomos como máximo, siendo dif́ıcilmente viable para aquellos constituidos por mi-
les de átomos en ausencia de recursos de supercomputación. Por otro lado, la Mecánica Molecular o
clásica (MM) ignora los movimientos electrónicos y calcula la enerǵıa de una molécula o conjunto de
moléculas únicamente en función de la disposición de los núcleos atómicos por lo que resulta el nivel
de teoŕıa generalmente escogido para estudiar macromoléculas biológicas como el DNA o las protéınas.
425
La Dinámica Molecular (DM) se basa en los principios de la mecánica clásica y constituye uno de los
métodos más utilizados para la simulación de macromoléculas biológicas como el DNA o las protéınas.
Ambas técnicas, no obstante, se pueden combinar en aproximaciones de Mecánica Cuántica / Mecánica
Molecular (QM/MM). Con este método h́ıbrido una pequeña parte del sistema implicada, por ejemplo,
en una reacción qúımica, se estudia mediante QM mientras que el resto de los átomos se considera
mediante MM.
17.2. Mecánica molecular
Los cálculos de MM se basan en la aproximación de Born-Oppenheimer que permite separar los mo-
vimientos del núcleo y de los electrones. Se considera que, debido a que la masa del núcleo es muy
superior a la de los electrones, éstos pueden adaptarse rápidamente a cualquier cambio en las posiciones
de los núcleos. Por lo tanto, la enerǵıa de una molécula en su estado basal, puede considerarse como
una función de las coordenadas de los núcleos atómicos. Esta función se la denomina Campo de Fuerzas
o Force Field. Los cambios que se producen en la enerǵıa potencial de un sistema pueden representarse
como una superficie, denominada superficie de enerǵıa potencial [11]. Uno de los objetivos en modelado
molecular es encontrar los puntos mı́nimos en la superficie energética que corresponden a estructuras
moleculares optimizadas. También es importante encontrar los “puntos de silla” (puntos de pendiente
igual a cero en cualquier dirección). Estos puntos se consideran barreras de mı́nima enerǵıa en los
caminos que conectan los distintos mı́nimos y que corresponden a los estados de transición.
17.2.1. El campo de fuerzas
Los primeros campos de fuerzas para simulaciones biomoleculares fueron desarrollados por primera vez
hacia 1970 [16, 41, 42, 44]. Desde entonces un gran número de campos de fuerza emṕıricos para MM han
sido desarrollados para la simulación de protéınas, ácidos nucleicos, ĺıpidos y otras moléculas biológicas
[40]. Es importante diferenciar entre los programas empleados en simulación y los parámetros de MM
desarrollados para ellos dado que en muchas ocasiones los nombres coinciden. Entre los programas
más usados para las simulaciones de DM por ordenador de moléculas biológicas destacan AMBER [9],
CHARMM [5], GROMOS [38], NAMD [32] y TINKER [34]. Un campo de fuerzas está formado por la
función de enerǵıa potencial y por los parámetros emṕıricos usados por cada uno de los términos. Una
caracteŕıstica importante de un campo de fuerzas es su capacidad de poder ser transferido, es decir, una
misma serie de parámetros puede ser empleada para el estudio de distintas moléculas relacionadas entre
śı [21]. La mayoŕıa de los campos de fuerzas empleados para sistemas moleculares se pueden definir
mediante una ecuación con dos componentes principales que describen las interacciones enlazantes y
no enlazantes del sistema:
Etotal = Eenlace + Eángulo + Etorsional︸ ︷︷ ︸
Enlazantes
+Evdw + Eelec︸ ︷︷ ︸
No enlazantes
(17.1)
Los términos enlazados incluyen contribuciones debidas a los enlaces covalentes, ángulos de valencia
y ángulos torsionales propios e impropios. Los términos no enlazados se definen por un término de
atracción-repulsión de tipo Lennard-Jones para las fuerzas de van der Waals y un término Coulómbico
para las interacciones electrostáticas. En la Figura 17.1 se esquematizan ambos tipos de interacciones.
A continuación se describen los componentes individuales más importantes que forman parte del campo
de fuerzas en MM.
Figura 17.1: Esquema de las interacciones enlazantes y no enlazantes que se tienen en cuenta en un
campo de fuerzas de MM.
17.2.2. Términos enlazados
Término de enlace
El término de enlace (bond stretching) se encarga de mantener las longitudes de enlace cercanas a los
valores de equilibrio medidos experimentalmente. En términos generales la enerǵıa potencial para un
enlace covalente se define mediante una función de Morse representada en rojo en la Figura 17.2.
Sin embargo, el potencial de Morse no se emplea generalmente en los campos de fuerzas dado que, en
los cálculos de MM, raras veces los enlaces se desv́ıan significativamente de sus valores de equilibrio.
Por ello se emplean expresiones más simples como la ley de Hooke por la cual la enerǵıa vaŕıa en función
del desplazamiento desde la longitud de referencia del enlace req:
Eenlace =
Kr
2
(r − req)2 (17.2)
Término de ángulo
La variación de los ángulos (angle bending) con respecto a sus valores de referencia también se repre-
senta mediante un potencial armónico de Hooke:
Eángulo =
Kθ
2
(θ − θeq)2 (17.3)
Término de torsión
El término de torsión (torsional term) describe la variación de la enerǵıa asociada a la rotación alre-
dedor de un enlace B-C dentro de una serie de cuatro átomos A-B-C-D donde A-B, B-C, y C-D están
Figura 17.2: Forma del potencial de enlace según la representación de Morse y según una representación
armónica.
unidos. Se caracteriza por presentar una periodicidad en el ángulo φ: si el enlace rota 360º la enerǵıadebe volver al mismo valor. Su perfil energético se expresa como una serie de Fourier :
Etorsional =
Vn
2
[1 + cos(nφ− γ)] (17.4)
La constante Vn determina la altura de la barrera de torsión alrededor del enlace B-C, n describe la
multiplicidad (el número de mı́nimos en la función cuando se rota 360º, φ el ángulo de torsión y γ el
ángulo de fase (indica en qué punto pasa la torsión por mı́nimo energético).
17.2.3. Términos no enlazados
Las moléculas y los átomos independientes interaccionan entre ellos a través de interacciones que no
dependen de una relación espećıfica de enlace entre átomos. Estos términos se agrupan en dos grupos:
interacciones de van der Waals e interacciones electrostáticas.
Interacciones de van der Waals
La interacción de van der Waals entre dos átomos se origina a partir de un balance entre fuerzas
atractivas y repulsivas. Esta enerǵıa de interacción vaŕıa en función de la distancia entre ambos átomos
como muestra la Figura 17.3. La enerǵıa de interacción es cero a una distancia interatómica infinita (e
incluso despreciable a distancias relativamente cortas). Al reducirse la distancia, la enerǵıa disminuye
hasta llegar a un mı́nimo. Después, la enerǵıa crece rápidamente al continuar disminuyendo la distancia.
Figura 17.3: Potencial de Lennard-Jones.
Las interacciones de atracción y repulsión entre átomos y moléculas pueden ser calculadas usando
mecánica cuántica. Sin embargo, ya que en muchos de los sistemas a estudiar existe un gran número
de interacciones de van der Waals, se hace necesario emplear una función que las calcule de manera
más eficiente. La función más conocida es la función de Lennard-Jones:
Evdw(rAB) =
aAB
r12
AB
− bAB
r6
AB
(17.5)
Donde a y b son dos constantes espećıficas del par de átomos A y B.
Interacciones electrostáticas
La distribución de la carga en una molécula se puede representar como una ordenación de las cargas
puntuales. Estas cargas reproducen las propiedades electrostáticas de la molécula. En el caso de que
las cargas estén centradas en los núcleos, se las denomina cargas atómicas parciales. La interacción
electrostática de una molécula se calcula, por tanto, como la suma de las interacciones entre pares de
cargas puntuales según la ley de Coulomb:
Eelec(rAB) =
qAqB
4πε0rAB
(17.6)
qA y qB son las cargas puntuales de cada átomo, rAB la distancia entre ellos y ε0 la constante dieléctrica
del medio que las separa. Si las cargas puntuales de dos átomos son contrarias éstos se atraerán entre
śı, pero si las cargas son del mismo signo se repelerán (Figura 17.4).
Figura 17.4: Potencial Coulómbico.
17.2.4. Parametrización del campo de fuerzas
Los campos de fuerzas contienen un gran número de parámetros, incluso si están diseñados para modelar
sistemas pequeños. El objetivo a la hora de desarrollar estos parámetros es encontrar un modelo capaz
de aproximarse lo mejor posible a los valores experimentales. Por ello, una fase importante en la
parametrización consiste justamente obtener estos parámetros de datos experimentales. En MM, esta
información proviene de datos estructurales, energéticos o electrónicos. Además, para complementar
estos datos experimentales, se recurre a los cálculos cuánticos ab initio que son capaces de reproducir
resultados experimentales en muchos sistemas. El gran inconveniente del uso de estos métodos para la
obtención de parámetros es que se requieren muchos recursos computacionales para asegurarse de que
los datos ab initio son suficientemente precisos.
17.2.5. Minimización de enerǵıa
Como se dijo anteriormente, la superficie de enerǵıa potencial de un sistema viene definida por el modo
en el que la enerǵıa de las moléculas vaŕıa en función de sus coordenadas. En modelado molecular,
es necesario estudiar los puntos mı́nimos en la superficie de enerǵıa potencial que corresponden a
los estados estables del sistema. Cualquier cambio en esta configuración produciŕıa un incremento de
la enerǵıa potencial [21].En muchos casos puede haber un gran número de mı́nimos en la superficie
energética. Aquel punto con el mı́nimo energético más bajo se le conoce como mı́nimo global mientras
que los demás puntos se los denomina mı́nimos locales (Figura 17.5). En un punto mı́nimo, la primera
derivada de la función de potencial f con respecto a las coordenadas cartesianas o internas (xi) es igual
a cero mientras que las segundas derivadas son positivas (Ecuación 17.7).
El punto energético más elevado en el camino entre dos mı́nimos se le conoce como punto de silla y
corresponde al estado de transición. El proceso que permite identificar las geometŕıas del sistema que
corresponden a los puntos de mı́nima enerǵıa potencial se denomina algoritmo de minimización.
Los métodos más empleados para la minimización de enerǵıa en modelado molecular son los basados
en las derivadas de la función de potencial, especialmente el de “descenso más pronunciado” (steepest
descent) y el de “gradiente conjugado” (conjugate gradient). Ambos modifican las coordenadas de los
átomos mientras desplazan el sistema hacia el del punto de mı́nima enerǵıa. El primero realiza una
búsqueda siguiendo la máxima pendiente y se suele emplear cuando la estructura a estudiar está muy
alejada del mı́nimo. Al aproximarse a este punto mı́nimo, el método de steepest descent se vuelve
demasiado oscilante y por ello puede sustituirse con el método de gradiente conjugado que evita estas
oscilaciones.
∂f
∂xi
= 0 :
∂2f
∂x2
i
> 0 (17.7)
Figura 17.5: Esquema de la superficie energética.
Existe un gran número de algoritmos de minimización además de los expuestos anteriormente. Dado
que este caṕıtulo, al estar destinado a ofrecer una visión general de las técnicas de simulación, no ofrece
los detalles de los distintos algoritmos, se recomienda la consulta de bibliograf́ıa más detallada [17, 21].
17.3. Simulaciones de dinámica molecular
Los métodos de minimización generan configuraciones individuales de mı́nima enerǵıa que, en mu-
chos casos, pueden ser suficientes para predecir ciertas propiedades de un sistema. Sin embargo, estos
métodos dif́ıcilmente pueden ser aplicados cuando queremos estudiar las propiedades estructurales y
termodinámicas de sistemas macromoleculares que contienen numerosos puntos mı́nimos de enerǵıa.
Los métodos de simulación computacional son capaces de obtener una muestra representativa de las
configuraciones de estos sistemas macroscópicos a una determinada temperatura. Existen dos grandes
técnicas de simulación computacional: el método de Monte Carlo y las simulaciones de Dinámica Mole-
cular (DM). Mientras que el primero se basa en analizar la enerǵıa de distintas coordenadas generadas
de manera aleatoria, la DM permite generar una trayectoria de puntos que evolucionan con el tiempo
siguiendo la segunda ecuación de Newton. Se trata por tanto de un método determinista, es decir,
el estado de un punto de la trayectoria permite predecir el estado del siguiente. En este caṕıtulo, se
describirán los principios y aplicaciones de las simulaciones de dinámica molecular.
17.3.1. Cálculo de las fuerzas
Como se dijo al inicio del caṕıtulo, La DM está basada en los principios de la mecánica molecular :
los núcleos atómicos son suficientemente pesados como para ser considerados como part́ıculas clásicas
cuya dinámica puede ser estudiada mediante la segunda ecuación de Newton, F = ma, cuya forma
diferencial se escribe:
d2y
dt2
=
Fxi
mi
(17.8)
siendo Fxi la fuerza aplicada a la part́ıcula i en la dimensión x y t el tiempo. A partir de unas
coordenadas (x0) y velocidades iniciales (v0) se determina la enerǵıa potencial y la enerǵıa cinética del
sistema, siendo la enerǵıa total la suma de ambas:
Etotal = Epotential(x0) + Ecinética(v0) (17.9)
La evolución temporal del sistema se puede seguir aplicando métodos de integración numérica. De este
modo se obtienen pequeñasetapas sucesivas separadas en el tiempo por un intervalo fijo ∂ti (tiempo
de integración). La fuerza que actúa sobre cada part́ıcula en un instante de tiempo t se determina
mediante la derivada de la enerǵıa potencial con respecto a las coordenadas:
F = −
∂Epotencial
∂x
= 0 (17.10)
El valor de la fuerza se asume constante durante cada paso de integración. Sumando sobre cada núcleo
las fuerzas de interacción con las demás part́ıculas del sistema se obtiene la fuerza total a partir de la
cual podemos determinar las aceleraciones de las part́ıculas (∂
2y
∂t2
=
Fxi
mi
) que se combinan a continuación
con las posiciones y velocidades a tiempo t para calcular las posiciones y velocidades a tiempo t+ ∂t.
17.3.2. Integración de las ecuaciones de movimiento
El método más simple para integrar las ecuaciones de movimiento es mediante algoritmos que aplican
las series de Taylor. Empleando esta aproximación, la posición de una part́ıcula a tiempo t + ∂t se
expresa en función de su posición, velocidad y aceleración:
xi(t = ∂t) = xi(t) +
∂xi
∂t
+
1
2
∂2xi
∂t2
(17.11)
Uno de los métodos más utilizados para integrar las ecuaciones de movimiento en simulaciones de DM
es el algoritmo de Verlet [37]. Éste usa las posiciones y aceleraciones a tiempo t y las posiciones en la
etapa previa, xi(t− ∂t), para calcular las nuevas posiciones a t+ ∂t, xi(t+ ∂). El tiempo de integración
se selecciona teniendo en cuenta que un tiempo de integración demasiado grande dará lugar a ines-
tabilidades por solapamiento de altas enerǵıas pero un tiempo demasiado corto estudiará únicamente
una región muy limitada del espacio conformacional. En el caso de biomoléculas como protéınas y
DNA, los movimientos más rápidos provienen de las vibraciones de los enlaces entre átomos pesados
e hidrógeno (X-H) que tienen lugar en la escala del femtosegundo. Sin embargo, estos movimientos
rápidos contribuyen poco al comportamiento global de la molécula y supondŕıan emplear más recursos
computacionales. Para evitarlo, se utiliza el algoritmo SHAKE [35] que mantiene constantes las longi-
tudes de los enlaces X-H en sus valores de equilibrio y permite emplear un tiempo de integración mayor
(2 fs) reduciendo el coste computacional.
17.3.3. Condiciones de ĺımite periódico
El agua juega un papel esencial en todos los organismos vivos. Por ejemplo, muchas enzimas necesitan
moléculas de agua para poder desarrollar sus funciones biológicas [27]. Las simulaciones de DM, por
tanto, deben intentar reproducir las interacciones de un solvente acuoso con el sistema biológico a
estudiar, bien de manera expĺıcita o impĺıcita. Durante una simulación, el sistema biológico incluido
en una caja de agua incluye átomos que se localizan en el borde del área de simulación. Para evitar
anomaĺıas y asegurar una inmersión completa del soluto en el solvente, se emplean condiciones de ĺımite
periódico de manera que el sistema se considera rodeado por réplicas iguales en todas las direcciones
y el cálculo de la simulación se desarrolla como si el sistema fuera infinito en el espacio (Figura 17.6).
La forma y tamaño de las cajas depende de la geometŕıa del sistema, usándose cajas octaédricas y
rectangulares para el estudio de ácidos nucleicos y protéınas (Figura 17.7).
Figura 17.6: Condiciones de ĺımite periódico en dos dimensiones.
17.3.4. Cálculo de las interacciones no enlazantes
El cálculo de los términos no enlazantes en una DM es el proceso más costoso computacionalmente.
Una aproximación muy utilizada para acelerar este proceso consiste en anular las interacciones entre
Figura 17.7: Caja rectangular y octaedro truncado.
pares de átomos que están separados por una distancia mayor a una distancia determinada o distancia
ĺımite. Es importante tener en cuenta que la distancia ĺımite debe ser inferior a la distancia entre
cualquier punto de la caja hasta cualquiera de sus copias vecinas.
El uso de esta distancia ĺımite, sin embargo, afecta especialmente al término electrostático debido a
que esta interacción es de largo alcance. Usando el sumatorio de Ewald [14] se consigue computar la
interacción electrostática total sumando las interacciones con infinitos sistemas que son réplicas del
original.
17.3.5. Preparación y ejecución de una DM
En esta sección discutiremos algunas etapas del desarrollo de las simulaciones de DM. Un ejemplo de
un protocolo para ejecutar simulaciones de DM se observa en la Figura 17.8.
Figura 17.8: Protocolo general para la ejecución de DM.
En primer lugar, es necesario establecer una configuración inicial del sistema que puede provenir de
datos experimentales (cristalograf́ıa de rayos X, NMR), de modelos teóricos o de una combinación
de ambos. Una vez optimizada la estructura, es necesario asignar velocidades iniciales a cada uno
de los átomos. Esto se realiza de manera aleatoria a partir de una distribución Maxwell-Boltzmann
a una determinada temperatura. Finalmente, es preciso mantener fijas algunas de las condiciones de
simulación: número de part́ıculas (N), volumen (V), temperatura (T), presión (P) o enerǵıa total del
sistema (E). Al combinarlas se obtienen simulaciones microcanónicas (NVE), isotérmico-isobáricas
(NPT) y canónicas (NVT) siendo NPT y NVT las más empleadas en DM de protéınas y ácidos
nucleicos.
Una vez que el sistema ha sido preparado y optimizado y las velocidades han sido asignadas, puede
comenzar la simulación propiamente dicha. Las simulaciones de DM se componen de dos etapas: una fase
de equilibrado y una fase de producción. El objetivo del equilibrado es llevar al sistema a un estado de
equilibrio a partir de la configuración inicial. Durante esta fase se monitorizan varios parámetros como la
enerǵıa potencial, la temperatura y la densidad hasta que se estabilizan. Para conseguir un equilibrado
óptimo, en ocasiones se aplican restricciones al sistema, liberándolas a continuación lentamente para
permitir su adaptación a las condiciones deseadas. Este proceso suele durar entre 200 y 500 picosegundos
aunque en ciertos casos conviene equilibrar el sistema durante varios nanosegundos.
Después de un correcto equilibrado comienza la producción en la cual permitimos la evolución del
sistema obteniendo una trayectoria que será analizada a continuación. Cuanto mayor sea el tiempo de
simulación, se explorará más espacio conformacional y por lo tanto los resultados serán más precisos.
En general se estima que el tiempo de simulación debeŕıa ser al menos diez veces más largo que la escala
temporal del proceso a estudiar. En la práctica esto significa que con los ordenadores actuales sólo se
pueden simular procesos que tienen lugar en la escala de unos pocos nanosegundos aunque gracias a
la mejora en los algoritmos de paralelización, la capacidad de almacenamiento en discos y el uso de
supercomputadores se ha podido pasar de simulaciones de 10 ns a las de 1 microsegundo [29].
17.3.6. Simulaciones de dinámica molecular de macromoléculas biológicas
Los primeros métodos de simulación de DM fueron desarrollados por Alder y Wainwright en 1957
empleando un modelo de esferas sólidas en el cual los átomos interactuaban únicamente mediante
interacciones perfectas ([2]. En 1964, Rahman realizó la primera simulación con potenciales continuos
para reproducir las interacciones atómicas con mayor realismo (Rahman 1964). Ya durante los años
70, las simulaciones de DM se empezaron a aplicar a sistemas más complejos y en 1976 se realizó la
primera simulación biomolecular investigando las isomerización del retinal [43] y la primera simulación
de DM de una protéına, el inhibidor de la tripsina pancreática bovina (BPTI), por McCammon en 1977
[24]. Esta primera DM teńıa una extensión de 9.2 ps e inclúıa 500 átomos. En la actualidad, gracias al
desarrollo de las teoŕıas de campos de fuerza y al incremento en el poder de cálculo, somos capaces de
estudiar sistemas que contienen hastaun millón de átomos en la escala del micro y del milisegundo. De
hecho, el uso de superordenadores especialmente diseñados para simulaciones de DM como el llamado
“Anton” permiten realizar investigaciones que alcanzan una escala de tiempo en la cual tienen lugar
muchos procesos biológicos [39].
El número de publicaciones relacionadas con la teoŕıa de DM y su aplicación en el estudio sistemas
biológicos está creciendo extraordinariamente [1]. Un ejemplo es la evolución desde 1970 hasta la actua-
lidad del número de art́ıculos que investigan las protéınas y la DM incluidos en PubMed (Figura 17.9).
Por supuesto, los experimentos tienen un papel fundamental a la hora de validar los métodos de si-
mulación [19]. A la hora de desarrollar las investigaciones es imprescindible una acción interdisciplinar
en la cual experimentos teóricos complementen a los experimentales y viceversa. Gracias a esta mutua
colaboración se enriquecen las investigaciones y se generan resultados de mayor relevancia [7].
Las simulaciones de DM se emplean en la actualidad para estudiar prácticamente cualquier tipo de
macromolécula con un interés biológico o médico (protéınas, ácidos nucleicos, carbohidratos) [4]. Entre
las aplicaciones más importantes de las simulaciones de dinámica biomolecular destacan el estudio de
los cambios conformacionales (Elber 2005) o simulaciones de plegamiento y desplegamiento de protéınas
[12]. La DM es también un método muy apropiado para estudiar los canales iónicos o la membrana de
las protéınas cuyo estudio mediante métodos experimentales resulta complejo [3, 36].
Figura 17.9: Número de art́ıculos incluidos en PubMed que resultan de una búsqueda de los términos
“Molecular Dynamics & Proteins”.
Otra ventaja de las simulaciones de DM es que permiten refinar estructuras obtenidas mediante cristalo-
graf́ıa de rayos X o por NMR [6, 10] (ver Sección 7.7). Además, la simulación de sistemas biomoleculares
nos permiten calcular la enerǵıa libre de unión de pequeños ligandos a sus moléculas diana [20, 25] o
la enerǵıa de activación necesaria para las reacciones catalizadas por enzimas [40].
La DM se ha convertido en una herramienta esencial en el estudio del DNA aportando gran información
acerca de sus caracteŕısticas estructurales y dinámicas. El DNA se caracteriza por ser una molécula
muy flexible pudiendo presentar distintas conformaciones [30]. Los métodos experimentales son muchas
veces insuficientes para analizar estos cambios conformacionales y por ello se recurre a la simulación
para estudiar de la flexibilidad y estructura del DNA. Un ejemplo es el estudio de la desnaturalización
de una doble hélice de DNA simulado por DM demostrando que se trata de un proceso muy complejo
siguiendo distintas rutas [31]. Investigaciones similares fueron también desarrolladas para estudiar el
mecanismo de acción de fármacos antitumorales que interaccionan con el DNA. Se demostró mediante
DM que los fármacos covalentemente unidos al DNA evitan la fusión completa de la doble hélice
inducida por alta temperatura justificando su acción estabilizadora del DNA [7] (Figura 17.10).
Uno de los problemas de la DM como algoritmo de muestreo de la superficie de enerǵıa potencial de
los sistemas biomoleculares es que la trayectoria puede quedarse estancada explorando únicamente un
mı́nimo local [40]. Distintos métodos han sido desarrollados para superar esta dificultad y ampliar el
espacio conformacional a estudiar, por ejemplo, la aplicación de un potencial tipo “umbrella” para
forzar el muestreo siguiendo una determinada coordenada de reacción [33] o el acoplamiento de dis-
tintas trayectorias permitiendo que intercambien distintos parámetros entre ellas, como por ejemplo la
temperatura [28].
Gracias al uso de técnicas que amplifican el poder atomı́stico de las simulaciones de DM, se puede ex-
tender el margen de problemas biológicos a estudiar mediante modelado molecular. En las simulaciones
Figura 17.10: Modelo Molecular obtenido por simulación de DM del entrecruzamiento intercatenario
formado entre el fármaco antitumoral Mitomicina C y las dos hebras de una molécula de DNA. Adaptado
de Bueren-Calabuig et al. 2012 [8].
de modelos aproximados, llamados de “grano grueso” o coarse grained, ideados por Levitt y Warshel
[22, 23], un grupo de átomos se modela como un único sitio de interacción de manera que, al bajar el
detalle atomı́stico, se consiguen alcanzar escalas de tiempo biológicamente más relevantes para estudiar
procesos más lentos y sistemas más grandes [4, 18].
En 2013, Martin Karplus (Universidad de Estrasburgo y Universidad de Harvard), Michael Levitt (Uni-
versidad de Standford) y Ariel Warshel (Universidad de Southern California) fueron galardonados con el
Premio Nobel de Qúımica “por el desarrollo de modelos multiescala para sistemas qúımicos complejos”.
En los años 70, contribuyeron a desarrollar técnicas computacionales para modelar reacciones qúımicas
o el plegamiento de protéınas. Al aplicar las técnicas de QM/MM, consiguieron estudiar reacciones
enzimáticas con un gran nivel de detalle (ver Figura 17.11).
Figura 17.11: En 2013, El Premio Nobel de Qúımica reconoció las técnicas de modelado molecular por
ordenador que muestran, por ejemplo, cómo la lisozima hidroliza los enlaces glicośıdicos (en amarillo).
Los autores combinan métodos QM muy precisos estudiar la reacción (derecha), con métodos de MM
que describen el resto de la protéına (izquierda) [43].
17.4. Métodos h́ıbridos QM/MM
Esta técnica permite combinar potenciales cuánticos y mecánico-moleculares en un potencial h́ıbrido
QM/MM [13, 21, 26]. En este tipo de aproximación, una pequeña fracción del sistema, por ejemplo,
aquella implicada en una reacción qúımica, se analizan mediante una función QM mientras que el
potencial de los demás átomos del sistema se examina por cálculos clásicos de MM. Este método
combina la simplicidad y velocidad del tratamiento MM con el potencial de la qúımica cuántica que
permite el modelado de la formación y rotura de enlaces, aśı como la inclusión de la polarización
electrónica. Al dividir el sistema, somos capaces de realizar cálculos significativamente más grandes
que los que hubieran sido posibles únicamente con una aproximación QM pura y al mismo tiempo
podemos estudiar reacciones qúımicas que no son susceptibles de un modelado riguroso mediante MM.
En el estudio de un mecanismo de acción enzimático, la región QM correspondeŕıa normalmente al
sitio activo incluyendo los grupos reactivos de la enzima mientras que la región MM incluiŕıa la parte
mayoritaria de la protéına, aquella que no incluye los residuos implicados en la reacción (Figura 17.12).
La enerǵıa total ETOT para este tipo de sistemas se puede escribir de la siguiente forma:
ETOT = EQM + EMM + EQM/MM (17.12)
donde EQM y EMM corresponden a la enerǵıa de aquellas partes del sistema tratadas exclusivamente
con QM y MM, respectivamente. EQM/MM es la enerǵıa de interacción entre las partes mecánico-
cuánticas y mecánico-moleculares.
Se han implementado distintos métodos con el fin de estudiar un sistema mediante QM/MM, los cuales
difieren entre śı por el nivel de teoŕıa utilizado para la QM (semiemṕırico, ab initio, enlace de valencia o
funcional de densidad), por el modelo de MM y por el modo de representar el disolvente (expĺıcitamente
o usando un modelo simplificado). Otra diferencia importante es el modo de tratar la unión entre las
regiones QM/MM. En general, es preferible cortar enlaces no polares (como enlaces sencillos C–C) que
cortar enlaces insaturados o polares. El método más utilizado incluye simplemente un “link atom” o
átomo enlazante (normalmente hidrógeno) que asegura la conservación de la valencia.
Figura 17.12: Ejemplo de un esquema QM/MM utilizado en el estudio del mecanismo de acción de la
trans-Sialidasa de Trypanosoma cruzi [33]. El centro activo (en azul) seráestudiado por un método QM
mientras que los demás aminoácidos de la enzima y el solvente (en rojo) serán tratados mediante MM
clásica.
17.5. Programas y tutoriales de dinámica molecular
En la Tabla 17.1 se listan algunos de los programas más utilizados para realizar simulaciones de dinámica
molecular junto con tutoriales y enlaces a sus páginas web.
Programa/Tutorial URL
AMBER http://ambermd.org
Simulación de un fragmento de DNA http://ambermd.org/tutorials/basic/tutorial1
Uso de VMD con AMBER http://ambermd.org/tutorials/basic/tutorial2
Simulación de un fármaco http://ambermd.org/tutorials/basic/tutorial4b
Análisis de trayectorias con Ptraj http://ambermd.org/tutorials/basic/tutorial5
Dinámica Molecular dirigida http://ambermd.org/tutorials/advanced/tutorial10
QM/MM http://ambermd.org/tutorials/advanced/tutorial2
CHARMM http://www.charmm.org
GROMOS http://www.gromos.net
NAMD http://www.ks.uiuc.edu/Research/namd
TINKER http://dasher.wustl.edu/ffe
Tabla 17.1: Listado de enlaces a programas y tutoriales para realizar simulaciones de dinámica molecular.
http://ambermd.org
http://ambermd.org/tutorials/basic/tutorial1
http://ambermd.org/tutorials/basic/tutorial2
http://ambermd.org/tutorials/basic/tutorial4b
http://ambermd.org/tutorials/basic/tutorial5
http://ambermd.org/tutorials/advanced/tutorial10
http://ambermd.org/tutorials/advanced/tutorial2
http://www.charmm.org
http://www.gromos.net
http://www.ks.uiuc.edu/Research/namd
http://dasher.wustl.edu/ffe
17.6. Bibliograf́ıa
[1] S. A. Adcock and J. A. McCammon. Molecular dynamics: survey of methods for simulating the activity of proteins.
Chemical reviews, 106(5):1589–1615, 2006.
[2] B. J. Alder and T. E. Wainwright. Phase transition for a hard sphere system. The Journal of chemical physics,
27(5):1208, 1957.
[3] O. Beckstein and M. S. P. Sansom. A hydrophobic gate in an ion channel: the closed state of the nicotinic acetylcholine
receptor. Phys Biol, 3(2):147–159, 2006.
[4] D. W. Borhani and D. E. Shaw. The future of molecular dynamics simulations in drug discovery. J. Comput. Aided
Mol. Des., 26(1):15–26, 2012.
[5] B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. Charmm: A program
for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem., 4(2):187–217, 1983.
[6] A. T. Brunger and P. D. Adams. Molecular dynamics applied to x-ray structure refinement. Acc. Chem. Res.,
35(6):404–412, 2002.
[7] J. A. Bueren-Calabuig. Modelado Molecular de la Interacción de Fármacos Antitumorales y Nucleasas con el ADN.
Universidad Alcalá, 2011.
[8] J. A. Bueren-Calabuig, A. Negri, A. Morreale, and F. Gago. Rationale for the opposite stereochemistry of the major
monoadducts and interstrand crosslinks formed by mitomycin c and its decarbamoylated analogue at cpg steps in
dna and the effect of cytosine modification on reactivity. Org. Biomol. Chem., 10(8):1543, 2012.
[9] D. A. Case, T. E. Cheatham III, and T. Darden. The amber biomolecular simulation programs. Journal of . . . , 2005.
[10] J. Chen, H.-S. Won, W. Im, H. J. Dyson, and C. L. Brooks. Generation of native-like protein structures from limited
nmr data, modern force fields and advanced conformational sampling. J. Biomol. NMR, 31(1):59–64, 2005.
[11] C. J. Cramer. Essentials of Computational Chemistry: Theories and Models. Wiley, 2 edition, 2004.
[12] V. Daggett and A. Fersht. The present view of the mechanism of protein folding. Nat. Rev. Mol. Cell Biol.,
4(6):497–502, 2003.
[13] G. de M Seabra, R. C. Walker, M. Elstner, D. A. Case, and A. E. Roitberg. Implementation of the scc-dftb method
for hybrid qm/mm simulations within the amber molecular dynamics package. J. Phys. Chem. A, 111(26):5655–5664,
2007.
[14] P. P. Ewald. Die berechnung optischer und elektrostatischer gitterpotentiale. Ann. Phys., 369(3):253–287, 1921.
[15] F. Gago. Métodos computacionales de modelado molecular y diseño de fármacos. Monograf́ıas de la Real Academia
Nacional de Farmacia, 2009.
[16] A. T. Hagler, E. Huler, and S. Lifson. Energy functions for peptides and proteins. i. derivation of a consistent force
field including the hydrogen bond from amide crystals. Journal of the American Chemical, 1974.
[17] F. Jensen. Introduction to Computational Chemistry. Wiley, 2 edition, 2006.
[18] S. C. L. Kamerlin, S. Vicatos, A. Dryga, and A. Warshel. Coarse-grained (multiscale) simulations in studies of
biophysical and chemical systems. Annu Rev Phys Chem, 62:41–64, 2011.
[19] M. Karplus and J. A. McCammon. Molecular dynamics simulations of biomolecules. Nature structural biology,
9(9):646–652, 2002.
[20] P. A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, M. Lee, T. Lee, Y. Duan, W. Wang, O. Donini,
P. Cieplak, J. Srinivasan, D. A. Case, and T. E. Cheatham. Calculating structures and free energies of complex
molecules: combining molecular mechanics and continuum models. Acc. Chem. Res., 33(12):889–897, 2000.
[21] A. Leach. Molecular Modelling: Principles and Applications (2nd Edition). Prentice Hall, 2 edition, 2001.
[22] M. Levitt. A simplified representation of protein conformations for rapid simulation of protein folding. Journal of
molecular biology, 104(1):59–107, 1976.
[23] M. Levitt and A. Warshel. Computer simulation of protein folding. Nature, 253(5494):694–698, 1975.
[24] J. A. McCammon, B. R. Gelin, and M. Karplus. Dynamics of folded proteins. Nature, 267(5612):585–590, 1977.
[25] J. Michel, N. Foloppe, and J. W. Essex. Rigorous free energy calculations in structure-based drug design. Mol. Inf.,
29(8-9):570–578, 2010.
[26] A. J. Mulholland. Chemical accuracy in qm/mm calculations on enzyme-catalysed reactions. Chem Cent J, 2007.
[27] D. R. Nutt and J. C. Smith. Molecular dynamics simulations of proteins: can the explicit water model be varied? J.
Chem. Theory Comput., 3(4):1550–1560, 2007.
[28] A. Patriksson and D. van der Spoel. A temperature predictor for parallel tempering simulations. Phys Chem Chem
Phys, 10(15):2073–2077, 2008.
[29] A. Pérez. Estudio de mecanismos de interacción macromolecular. Universidad de Barcelona, 2008.
[30] A. Pérez, F. J. Luque, and M. Orozco. Frontiers in molecular dynamics simulations of dna. Acc. Chem. Res.,
45(2):196–205, 2012.
[31] A. Pérez and M. Orozco. Real-time atomistic description of dna unfolding. Angewandte Chemie (International ed,
49(28):4805–4808, 2010.
[32] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, and
K. Schulten. Scalable molecular dynamics with namd. J. Comput. Chem., 26(16):1781–1802, 2005.
[33] G. Pierdominici-Sottile, N. A. Horenstein, and A. E. Roitberg. Free energy study of the catalytic mechanism of try-
panosoma cruzitrans-sialidase. from the michaelis complex to the covalent intermediate. Biochemistry, 50(46):10150–
10158, 2011.
[34] J. W. Ponder and F. M. Richards. An efficient newton-like method for molecular mechanics energy minimization of
large molecules. J. Comput. Chem., 8(7):1016–1024, 1987.
[35] J. P. Ryckaert, G. Ciccotti, and H. Berendsen. Sciencedirect.com - journal of computational physics - numerical
integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.
Journal of Computational Physics, 1977.
[36] M. S. P. Sansom, P. J. Bond, S. S. Deol, A. Grottesi, S. Haider, and Z. A. Sands. Molecular simulations and lipid-
protein interactions: potassium channels and other membrane proteins. Biochem. Soc. Trans., 33(Pt 5):916–920,
2005.
[37] D. Schiff and L. Verlet. Ground state of liquid helium-4 and helium-3. Phys. Rev., 160(1):208–218, 1967.
[38] W. R. P. Scott, P. H. Hünenberger, I. G. Tironi, A. E. Mark, S. R. Billeter, J. Fennen, A. E. Torda, T. Huber,
P. Krüger, and W. F. van Gunsteren. The gromos biomolecular simulation program package. J. Phys. Chem. A,
103(19):3596–3607, 1999.
[39] D. E. Shaw, K. J. Bowers, E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson,
K.Lindorff-Larsen, P. Maragakis, M. A. Moraes, R. O. Dror, S. Piana, Y. Shan, B. Towles, J. K. Salmon, J. P.
Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. Deneroff, and B. Batson. Proceedings of the conference
on high performance computing networking, storage and analysis, 2009.
[40] M. W. van der Kamp and A. J. Mulholland. Computational enzymology: insight into biological catalysts from
modelling. Nat. Prod. Rep., 25(6):1001, 2008.
[41] P. K. Warme, F. A. Momany, R. S V, R. W. Tuttle, and S. H. A. Computation of structures of homologous proteins.
alpha-lactalbumin from lysozyme. Biochemistry, 13(4):768–782, 1974.
[42] A. Warshel and M. Karplus. Calculation of ground and excited state potential surfaces of conjugated molecules. 1.
formulation and parametrization. Journal of the American Chemical Society, 16:5612–5625, 1972.
[43] A. Warshel and M. Levitt. Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization
of the carbonium ion in the reaction of lysozyme. J Mol Biol, 103(2):227–49, 1976.
[44] A. Warshel, M. Levitt, and S. Lifson. Consistent force field for calculation of vibrational spectra and conformations
of some amides and lactam rings. Journal of Molecular Spectroscopy, 33(1):84–99, 1970.