Métodos de Simulación en Química y Bioquímica

1. Métodos de Simulación en Química y Bioquímica Tópicos de Fisicoquímica de Sistemas Biológicos 2007

2. Simulación moléculas pequeñas aisladas (fase gaseosa) ha llegado niveles de precisión muy altos. (QM) Sin embargo, procesos reactivos en solución y/o entornos complejos como biomoléculas son mucho más difíciles de investigar. ESE ES NUESTRO OBJETIVO!!!! Desarrollo metodologías Aplicaciones a reacciones en solución. Aplicaciones: reacciones en biomoléculas

3. Desafío: Introducir en las predicciones realizadas, los efectos del entorno y de la temperatura. Cuando se logre ese objetivo, la Química Computacional será capaz de realizar predicciones realmente relevantes.

4. METODOS DE SIMULACIONCOMPUTACIONAL (en Química y Bioquímica) • Herramientas informáticas (homologías, bases de datos, etc) • Herramientas fisicoquímicas implementadas en programas 1) Modelos: recetas para evaluar superficie energía potencial (SEP). (potenciales clásicos, QM, etc) 2) Esquemas para extraer información de SEP. Monte Carlo, Dinámica Molecular.

5. MODELOS CUANTICOS • Modelos basados en la mecánica cuántica: • Resuelven (o tratan!) ecuación Schrödinger. • Aproximación básica: Born-Oppenheimer • Distintos niveles aproximación: • semiempíricos, HF, DFT, etc • Indispensable para: procesos reactivos, estados excitados. • En principio: válido para cualquier sistema. • En la práctica: limitaciones computacionales en • cuanto tamaños.

6. Aproximación fundamental: Orbital • Función de onda: determinante de Slater. Equivale a asignar un electrón a cada orbital. • Matemáticamente: ecuaciones de Hartree-Fock. (SCF)

7. Orbitales: se expanden con funciones de base. • Número de funciones: determina calidad del resultado. • Base infinita: límite de HF (no es correcto) • Siglas usuales: 3-21G, 6-31G, 6-31G**, etc.

8. Costo computacional: • Crece con N4 (eso es mucho!!!) Qué se hace? Métodos aproximados: semiempíricos: Usa HF, pero desprecia algunas integrales, y mete resultados experimentales.

9. Métodos semiempíricos • Muy eficientes. • Problemas con metales, estados de transición. • AM1, PM3, ZINDO, etc…… • Siglas: distintos niveles de aproximación y parametrizaciones diferentes.

10. Métodos post-HF • Corrigen error en aproximación orbital.. • Son muy costosos computacionalmente. • MP2 (perturbativo) • CI (variacional). Sirve para estados excitados.

11. Método alternativo: DFT • Teoría del funcional de la densidad • Variable fundamental: densidad • Eficiente computacionalmente. Da muy bien. Preferido en biología….

12. QM: DFT y Exchange-correlation con:

13. Que cosas se calculan? • Distribución de cargas. • Factibilidad de reacciones (ΔE y barreras) • Propiedades estructurales (optimización de geometrías)

14. MODELOS CLASICOS • Modelos basados en la ideas clásicas: • (ley Coulomb, resortes, etc) • Distintos niveles aproximación según complejidad • función propuesta, ejemplo, potenciales polarizables. • Numéricamente muy eficientes. • En principio: no son válidos para procesos reactivos. • Parametrización es clave y requiere artesanía!

16. Química Cuántica Tradicional: optimizaciones de geometría Corresponde a analizar mínimos de la SEP (0 K) Cómo incluimos T? Existe una gran experiencia en métodos de simulación clásicos en los últimos 20 años, y más recientemente también con modelos cuánticos.

17. Ejemplos de fenómenos de interés químico y/o biológico que se pueden investigar mediante técnicas de simulación: • Reacciones químicas en solución. • Reacciones enzimáticas. QM • Interacciones de macromoléculas con ligandos. • Predicción de estructura de macromoléculas. • Interacción entre macromoléculas. MM

18. ¿Cómo extraigo información de la SEP? V= V(R1, R2, R3, …..,Rn) Termodinámica Estadística Promedio temporal propiedad A: Promedio configuracional propiedad A Promedio sobre ensamble de sistemas: Densidad probabilidad Probabilidad encontrar configuración R

19. Ensambles: • Microcanónico: NVE constantes • Canónico: NVT constantes • Otros algo menos usuales: • Isotérmico-isobárico: NPT • Gran canónico: potencial quimico, V, T Corresponden a situaciones experimentales diferentes y se pueden implementar en esquemas de simulación.

20. A N,T y V constantes: Función de partición Energía cinética Energía Potencial (SEP)

21. En casos normales V depende sólo de R Término K Siempre igual! Término V, depende del sistema

22. En principio si conozco Q en forma analítica, puedo obtener todas las propiedades termodinámicas de mi sistema! Por ejemplo Esto sólo se puede hacer para gases ideales, además suponiendo modelos simplificados para grados de libertad internos (oscilador armónico, rotor rígido, etc),

23. Volviendo a las ecuaciones básicas: HIPOTESIS ERGODICA : En principio equivalentes.

24. MC Ecuación muy elegante, pero como no conozco Q, no conozco  Manera aproximada de sacar cosas de esa ecuación: Esquema Monte Carlo (MC), Metropolis et al 1952 MD Manera aproximada de usar esta otra: Dinámica Molecular (MD)

25. DINAMICA MOLECULAR Leyes Newton: F=m a Conociendo superficie de energia potencial: saco F=Grad E Se usan métodos de diferencias finitas, expandiendo en Serie de Taylor: r(t+h)=r(t) + h v(t) + ½ h2 a(t) + …… r(t-h)=r(t) – h v(t) + ½ h2 a(t) + …… Sumando da: ALGORITMO VERLET r(t +h) = 2 r(t) –r(t-h) + h2 a(t)

26. Receta para calcular posiciones a tiempo t+h, conociendo las aceleraciones a tiempos t-h , t y la aceleracion= Fuerza/masa para cada átomo del sistema. Se parte de r(t), y r(t-h), y velocidades salen de: v(t)=[r(t+h)-r(t-h)]/(2h) Alternativa más usada: Velocity Verlet r(t+h)=r(t)+hv(t) +1/2 h2 a(t) v(t+h)=v(t)+1/2 h [a(t)+a(t+h)] Mejor, condiciones iniciales posiciones y velocidades.

27. MD: muestrea configuraciones a E constante (ensamble microcanónico). Se puede demostrar fácilmente. Como se si hago las cosas bien? Se debe conservar E total. Típicamente 1 en 104 K(t)= (pi(t))2/(2mi) = kT(3N-Nc)/2 Definicion microscópica de T V(t)= dada por SEP, E(t)=V(t)+K(t) 3N-Nc: grados de libertad internos del sistema Como influimos?. Elección de time step h

28. Ver escalas. Fluctuaciones E total Deben ser menores a 1/104

29. Movimiento lento h 0.1-05 fs h 1-2 fs Movimiento rápido

30. Time step: 0.1-0.2 fs vibraciones enlaces covalentes 1-2 fs doblamientos, torsiones, rotaciones, etc. Tiene que ver con escalas de tiempo de movimientos. Uniones fuertes: movimientos rápidos , menor h preciso para integrar bien. Truco: usar uniones rígidas: constraint dynamics Metodo usual: SHAKE, Ciccotti et al 1977. Se agrega una fuerza extra, debido al constraint (ejemplo unión HO =1 angstrom) Matemáticamente: se usan multiplicadores de Lagrange

31. MD: tambien se puede adaptar para muestrear • en el ensamble Canónico (mas parecido a • situación experimental) T constante • Esquemas más simplificados: • Escalar velocidades directamente! • K(t)= mi vi(t)2/2 = kT(3N-Nc)/2 • Multiplicar velocidades por x=sqrt(Tdeseada/T(t)) • Eso se hace en equilibración, pero no en corridas • de producción.

32. 2) Algoritmo de Berendsen Se usa mucho en proteínas. Es algo más riguroso. El cambio en temperatura entre dos pasos es: T=1 + h/ (Tdeseada/T(t) –1), con  parámetro acoplamiento Entre el termostato y el sistema.

33. Esquema correcto para hacer MD en el conjunto canónico: se denomina algoritmo de Nose (1984) Se agrega grado de libertad para el termostato o Reservorio y es algo más complicado, pero más riguroso. Se agrega K y V para el reservorio a E total. Esquemas escalamiento tipo Berendsen: problemas: no equilibran (distribuyen E entre distintos grados de libertad). Pueden haber modos calientes y modos fríos.

34. Propiedades estructurales distancias, angulos, fluctuaciones, etc. Propiedades dependientes del tiempo ¿Qué se puede aprender de esto? Propiedades termodinámicas. valores medios E, cv. Cambios G.

35. Predicciones que se pueden realizar dependen del modelo empleado para la descripción de la superficie de energía potencial. Reacciones químicas: necesito Mecánica Cuántica. Muy costoso!! Procesos que no involucran ruptura o formación uniones covalentes: Pueden usar potenciales clásicos. Hay muchos trabajos en la Literatura en Bioquímica y Ciencias de Materiales.

36. Ejemplo: Bovine Pancreatic Trypsin Inhibitor (BPTI). Proteína pequeña: 536 átomos Muy estudiada

37. Medidas estructurales globales:

39. Otras medidas estructurales: funciones de correlación radial g(r). Por ejemplo para ver interacciones específicas (uniones H). Probabilidad de encontrar un átomo a una distancia r de otro átomo, comparada con una distribucion ideal. Dist. Ideal: Volumen de capa de ancho h es: 4/3(r + h)3 - 4/3(r)3  4r2h ; Nro partículas está dado por: 4r2h , donde  densidad (N/V)

40. Se calcula haciendo histogramas de toda la simulación. Ej: H2O, g(HO), se toman todas las distancias HO en cada paso, y se acumulan en un histograma. Luego se normaliza con el factor de volumen ideal.

41. Gráfico tiene este aspecto: Primer pico: interacciones directas. Sólido: muchos picos. Gas, estructura se pierde rápido. H...O(nitrato) O ...N(nitrato) O...O(nitrato)

43. Ejemplo de proceso dinámico: estudian la dinámica del plegado (Folding) y comparan la estructura plegada con la experimental obtenida por RMN

44. Método de Monte Carlo • A partir de una configuración de partida, se generan • configuraciones nuevas, moviendo las coordenadas de • cada átomo según: xi = x i-1 +  w, donde  es un • parámetro fijo, y w es un número generado al azar • entre 0 y 1. (Idem para y, z). • Para cada configuración i-sima se evalúa • p= exp(-E)/kT), donde E=Ei - E i-1 • Si p > 1, se acepta la movida. • Si p< 1, se compara p con otro número al azar • entre 0 y 1, ww. Si p< ww, se acepta la movida.

45. ¿Qué quiere decir aceptar movida? • Si acepto movida, uso esa configuración para la suma. • Si rechazo movida, vuelvo a poner en la suma la • configuración anterior. • Notar que si hice M movidas, la suma tendrá M • términos, independientemente del número de • aceptaciones.

46. Detalles Técnicos: • Elijo cuán lejos me muevo con el parámetro  • Si me muevo mucho, voy a rechazar casi todo, y • eso no es bueno. • Habitualmente, se elige , como para que se acepte • alrededor del 50% de las movidas. • Muestreo en el ensamble canónico (NVT)

47. Más detalles de MC: • no necesito programar fuerzas (derivadas) • Idea es consistente con que los sistemas visitan con mayor probabilidad configuraciones de menor energía, pero que eventualmente pueden visitar otras de mayor energía. • Configuraciones de muy baja energía darán origen a muchos rechazos, y por eso se cuentan muchas veces en suma • Problema: difícil que ocurran movimientos concertados

48. Esquema MC convencional: útil en sistemas sencillos. Ejemplo: cálculo solvatación en solutos orgánicos.

50. Proteína corismato mutasa