ES2440415A2

ES2440415A2 - Metodo y sistema de simulacion mediante dinamica molecular con control de estabilidad

Info

Publication number: ES2440415A2
Application number: ES201231189A
Authority: ES
Inventors: Roldán MARTÍNEZ ABÁN; Pablo ECHENIQUE ROBBA; Javier Sancho Sanz; Roberto MARTÍNEZ BENITO; Alfonso Campo Camacho; Álvaro LÓPEZ MEDRANO
Original assignee: Plebiotic S L; PLEBIOTIC SL
Current assignee: Plebiotic S L; PLEBIOTIC SL
Priority date: 2012-07-25
Filing date: 2012-07-25
Publication date: 2014-01-28
Anticipated expiration: 2032-07-25
Also published as: WO2014016447A1; ES2440415B1; ES2440415R1; US9727690B2; US20150278432A1; EP2879066A1

Abstract

Método y sistema de simulación mediante dinámica molecular con control de estabilidad. La presente invención tiene su aplicación dentro del campo de dinámica molecular, la cual consiste en métodos de computación para la predicción de la estructura y función de biomoléculas, y en particular de proteínas, mediante la simulación de los procesos de plegamiento de las mismas y de interacción con otras biomoléculas en un solvente. Más particularmente, la invención se refiere a un método y un sistema de control de estabilidad de las simulaciones y de elección del paso de tiempo usado en la integración numérica de las ecuaciones del movimiento. La invención logra reducir el tiempo de simulación de la dinámica molecular mediante la optimización de la elección paso de tiempo a través de un control adaptativo o permitiendo pasos de tiempo mayores corrigiendo las trayectorias en base a un control de la potencia.

Description

Método y sistema de simulación mediante dinámica molecular con control de estabilidad

Objeto de la invención

La presente invención tiene su aplicación dentro del campo de dinámica molecular, la cual consiste en métodos de computación para la predicción de la estructura y función de biomoléculas, y en particular de proteínas, mediante la simulación de los procesos de plegamiento de las mismas y de interacción con otras biomoléculas en un solvente. La Dinámica Molecular también permite simular sistemas más sencillos, como gases o fluidos, siendo la presente invención también aplicable. Más particularmente, la invención que aquí se describe, dentro de la Dinámica Molecular, se refiere a un método de control de estabilidad de las simulaciones y de elección del paso de tiempo usado en la integración numérica de las ecuaciones del movimiento.

El método de la invención se implementa en un ordenador o cualquier equipo electrónico programable adecuado, y aporta el efecto o ventaja de reducir el tiempo de simulación de la Dinámica Molecular mediante la optimización de la elección paso de tiempo a través de un control adaptativo o permitiendo pasos de tiempo mayores corrigiendo las trayectorias en base a un control de la potencia.

Antecedentes de la invención

La estructura tridimensional de una proteína determina su función biológica (véase Berg JM, Tymoczko JL, Stryer L., Biochemistry, 5th edition, W H Freeman, 2002). Por tanto, es de vital importancia la determinación de dicha estructura para, entre otros, el diseño de fármacos. Una forma de determinar la mencionada estructura es mediante técnicas de cristalización por rayos X o por resonancia magnética nuclear (véase P. Narayanan, Essentials of Biophysics, Second Edition, Anshan Ltd., 2010). No obstante, éstas técnicas no están exentas de imperfecciones y no permiten con facilidad el proceso opuesto, es decir, la determinación de la cadena de aminoácidos que permite obtener una cierta estructura (véase Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006). Estos hechos, unidos al aumento de la capacidad de procesamiento de los computadores, favorecido el desarrollo, en los últimos años, de la Dinámica Molecular que consiste en la simulación por ordenador del proceso de plegamiento de una proteína o de su interacción con otras moléculas.

Ha sido demostrado (véanse P. L. Freddolino, K. Schulten, Common Structural Transitions in Explicit-Solvent Simulations of Villin Headpiece Folding, Biophysical Journal Volume 97 October 2009 2338–2347, P. L. Freddolino,F. Liu, M. Gruebele, K.Schulten,Ten-Microsecond Molecular Dynamics Simulation of a Fast-Folding WW Domain, Biophysical Journal, Volume 94, Issue 10, L75-L77, 15 May 2008, S. K. Sadiq, G. De Fabritiis, Explicit solvent dynamics and energetics of HIV-1 protease flap-opening and closing, Proteins: Structure, Function, and Bioinformatics Volume 78, Issue 14, pages 2873–2885, 1 November 2010) que, la aproximación mediante la mecánica de Newton, es válida para modelar la dinámica (proceso de plegamiento o de interacción con otras moléculas) de grandes moléculas como las proteínas, la cual es conocida como Dinámica Molecular. Teniendo pues

en cuenta la segunda Ley de Newton y que, en ese caso cada átomo es considerado como un cuerpo de masa m, la

Fi

ecuación que determina la aceleración a a la que éste es sometido debido la existencia de fuerzas entre éste y

el resto de átomos de la propia proteína o del solvente, es la siguiente:

N

I

Fi = m · a

(1)

i=1

Siendo N el número de átomos. La aceleración a a la que está sometido un cuerpo es la derivada respecto del tiempo de la velocidad v

del mismo. .

Por tanto dicha velocidad puede ser obtenida de la integración respecto del tiempo de la mencionada aceleración a

De igual modo la posición x es obtenida como la integral respecto del tiempo de la velocidad v dado que ésta es la derivada de aquélla. De este modo la posición de cada átomo puede ser obtenida mediante integración a partir de la velocidad y ésta a su vez mediante integración de la aceleración que, conocidas las fuerzas que actúan sobre dicho átomo, nos proporciona la anterior ecuación 1.

De igual manera para similar las condiciones reales de las proteínas, y en general de las biomoléculas, se rodea a la molécula de un solvente, agua, simulando las condiciones de la célula. Este solvente es necesario para ciertos procesos biológicos, como el plegamiento de proteínas dónde el carácter hidrófilo o hidrofóbico es clave para determinar su estructura.

Esto requiere simular además de la proteína a estudio, miles de moléculas de agua que simulan el solvente, por lo que los sistemas a estudio cuentan con miles de átomos a simular.

Por motivos similares, para simular las condiciones experimentales de los sistemas macroscópicos, se introducen las llamadas condiciones periódicas de frontera (véase T. Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide. Interdisciplinary Applied Mathematics series, vol. 21. Springer: New York, pp272–6. 2002). De esta manera el sistema (biomoléculas más solvente) que queremos simular se rodea de copias periódicas de sí mismo en las 3 direcciones. Esta es una manera computacionalmente abordable de simular las condiciones de los sistemas, dónde se tienen un elevado número de moléculas (del orden del número de Avogadro).

Una manera eficiente de calcular la fuerza que actúa sobre un átomo debido a su interacción con estas copias, es el llamado Particle Mesh Ewald (T. Darden, D. M.York, L. G. Pedersen , Particle mesh Ewald: An N log(N) method for Ewald sums in large Systems, J. Chem. Phys. (Communication) 98, 10089-10092 (1993)). Antes de entrar en detalle sobre los diferentes tipos de fuerza y la problemática que plantean, es conveniente hablar de la base de la Dinámica Molecular, la resolución de la segunda ley de Newton, ecuación 1. Se trata de una ecuación diferencial de segundo orden en la posición y por tanto su resolución pasa por una integración.

Esta ecuación es equivalente a un sistema de ecuaciones de primer orden, las denominadas ecuaciones de Hamilton, en lo que se llama formalismo hamiltoniano (véase B. Leimkuhler, S. Reich, Simulating Hamiltonian dynamics, Cambridge University Press 2004).Dado que estas soluciones no son resolubles analíticamente en casos como la Dinámica Molecular, es necesaria una resolución numérica de dichas ecuaciones diferenciales.

Dicha resolución, llamada integración numérica, se basa en el paso de un tiempo continuo a uno discreto. La separación entre cada uno de los instantes de tiempo discreto se le denomina paso de tiempo o timestep. Partiendo de unas condiciones iniciales de posición y velocidad, se trata de hallar sus valores en todos los instantes posteriores. A la función que, a partir del valor de posición y velocidad en un instante, permite obtener un valor para el instante siguiente, se le denomina integrador. Por tanto, el integrador y el paso de tiempo son los que determinan la calidad de la simulación.

Dado que al pasar de un tiempo continuo a uno discreto estamos realizando una aproximación, cuando mayor es el paso de tiempo, peor es la resolución de la ecuación diferencial.

Puede ocurrir que si elegimos un paso de tiempo demasiado grande la solución que obtenemos no tenga nada que ver con la real y diverja. En ese caso se dice que para ese paso de tiempo el integrador es inestable. El rango de valores de paso de tiempo para los que un integrador es estable se denomina zona de estabilidad (véase E. Hairer,

S. P. Nørsett, Gerhard Wanner, Solving ordinary differential equations I: Nonstiff problems, second edition, Springer

Verlag, Berlin, 1993) Las fuerzas Fi a las que cada átomo está sometido son de distinta naturaleza en cuanto a su

magnitud, su variabilidad en el tiempo y son función de la posición de los restantes átomos de la proteína o de los átomos del solvente.

El proceso de integración numérica anteriormente comentado está basado en ir avanzando la simulación a intervalos de tiempo hasta alcanzar el tiempo final de la misma. Debido a la altísima frecuencia (variabilidad en el tiempo) de

algunas de las fuerzas Fi estos intervalos de tiempo son pequeñísimos, del orden de los femtosegundos (10-15

segundos). Los tiempos totales de simulación requeridos para por ejemplo simular el plegamiento completo de una proteína van desde los microsegundos hasta varios minutos (David Sheehan, Physical Biochemistry, Second Edition, John Wiley & Sons, 2009) por tanto el número de intervalos de tiempo necesarios para completar la simulación van

desde los hasta los . En cada uno de estos intervalos de tiempo es necesario hacer complejos cálculos

computacionales, principalmente los correspondientes a las fuerzas Fi , y además éstos crecen exponencialmente

con la complejidad del sistema. Todo esto hace que, incluso los ordenadores más potentes actualmente, necesiten un tiempo de computación elevadísimo, del orden de días, meses e incluso años para completar la simulación de sistemas moleculares de un tamaño medio (véanse Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD,. Journal of Computational Chemistry, 26:1781-1802, 2005). Se deduce fácilmente de esto que las técnicas conocidas en el estado del arte de la Dinámica Molecular presentan el problema técnico de que no permiten realizar simulaciones con la rapidez que la industria farmacéutica o los organismos de investigación están demandando.

Entre las mencionadas Fi , las de enlace, modelan vibraciones de altísima frecuencia y por tanto necesitan

efectivamente unos pequeñísimos intervalos de tiempo para ser integradas: de 1 a 10 femtosegundos. Por el

contrario, entre estas mismas Fi se encuentran otras como las de Van Der Waals o las electroestáticas (Tamar

Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006), que son fuerzas de largo alcance comparadas con las anteriores de enlace, que tienen una frecuencia mucho más pequeña, y para las que, en principio, unos intervalos de tiempo mucho mayores (del orden de dos órdenes de magnitud más que los correspondientes a las de enlace) podrían ser suficientes para ser integradas. Hay que destacar que estas fuerzas de largo alcance, al implicar interacciones de un átomo con el resto de los de la proteína o molécula, así como con los del solvente, tienen un coste computacional altísimo comparado con las de enlace que sólo tienen en cuenta interacciones entre un átomo y su vecino más próximo. Así mientras la evaluación de éstas últimas tiene una complejidad de O(N), para aquéllas la complejidad es de O(N2) (Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006).

El integrador que por su simplicidad y excepcional estabilidad en simulaciones de larga duración como las requeridas en Dinámica Molecular es el método de Leapfrog/Verlet/Stormer (véase L.Verlet., Computer ‘experiment’ on classical fluids: I. Thermodynamical properties of Lennard-Jones molecules. Phys.Rev.,159 (1):98-103, July 1967.) habitualmente conocido de forma abreviada como de Verlet. Es un integrador explícito y de orden 2, con la particularidad de que sólo requiere evaluar las fuerzas una vez cada iteración, mientras que otros integradores del mismo orden (p.ej. Runge Kutta orden 2) requieren 2 evaluaciones, y por tanto prácticamente el doble de coste computacional.

Como integrador numérico, Verlet tiene importantes propiedades matemáticas como la simplecticidad y la simetría temporal que también tienen las ecuaciones diferenciales que trata de integrar. Por lo tanto representa más fielmente el sistema que se está simulando, manteniendo ciertas propiedades físicas como la conservación de la energía y la reversibilidad temporal, y matemáticas como la conservación del volumen físico, lo que redunda en mayor estabilidad. Existen variantes de este método como el velocity Verlet o el position Verlet, ampliamente utilizadas, que también poseen estas importantes propiedades.

El uso de integradores simplécticos es típico en diversos campos de la Mecánica (véase B. Gladman, M. Duncan, J. Candy, Symplectic integrators for long-term integrations in celestial Mechanics, Celest. Mech. Dynam. Astron., 52, 221–240 1991), debido al descubrimiento de su mayor eficiencia respecto a integradores no simplecticos (véase M.

P. Calvo, J. M. Sanz-Serna, The developmentof variable-step symplectic integrators, with application to the two-body problem, SIAM J. Sci. Comput., 14, 936–952. 1993). Es también conocido que el uso de un paso de tiempo variable para un integrador simpléctico reduce su eficiencia (R. D. Skeel, C. W. Gear, Does variable step size ruin a symplectic integrator? Physica D, 60,311–313. 1992) y su estabilidad (véase R. D. Skeel, Variable step size destabilizes the Störmer/leapfrog/Verlet method, BIT, 33, 172–175. 1993). Debido a la complejidad y por tanto el comportamiento caótico de los sistemas que se simulan en Dinámica Molecular, con miles de grados de libertad, el seguimiento de un átomo en particular no es relevante. Una característica de los sistemas caóticos es su sensibilidad a las condiciones iniciales, de tal manera que una pequeña variación en una de las variables dinámicas lleva a una trayectoria completamente diferente. En un sistema molecular siempre tenemos una imprecisión inicial en las posiciones y velocidades iniciales del sistema, debido a problemas de medida de las posiciones de los átomos y a la temperatura, por lo que una trayectoria en particular no será representativa. Sin embargo estaremos interesados en variables estadísticas así como las configuraciones correspondientes a mínimos locales de la energía, que será las conformaciones hacia las que tiendan las trayectorias. En el caso de una proteína, corresponderá a su estado nativo. Por tanto en Dinámica Molecular, una alta precisión en las trayectorias no es necesario, lo que permite el uso de integradores no simplécticos, y el uso un paso de tiempo variable para mejorar la eficiencia. (véase S. J. Stuart, J. M. Hicks, M. T. Mury, An Iterative Variable-timestep Algorithm for Molecular Dynamics Simulations, Molecular Simulation, Volume 29, Issue 3 2003 , pages 177 – 186).

Tanto las técnicas, anteriormente comentadas y que veremos a continuación, de imposición de restricciones a las componentes de alta frecuencia de las fuerzas de enlace, como las basadas en el aumento del intervalo de integración de las fuerzas lentas, son modificaciones o variaciones del mencionado método Leapfrog/Verlet/Stormer.

Entre las técnicas de imposición de restricciones (ligaduras) para amortiguar las altas frecuencias se encuentran los conocidos como SHAKE (véase J.P. Ryckaert, G. Ciccotti and H.J.C. Berendsen., Numerical integration of ther Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys.,

23:327-341, 1977). RATTLE (véase H.C. Andersen. Rattle: a ‘velocity’ versión of the SHAKE algorithm for molecular dynamics calculations., J.Comput. Phys., 52:24-34, 1983.) y WIGGLE (Lee, Sang-Ho, K. Palmo, S. Krimm (2005). WIGGLE: A new constrained molecular dynamics algorithm in Cartesian coordinates, Journal of Computational Physics 210: 171–182 2005). Ambos se combinan con el método de Verlet o alguna de comentadas variantes. No obstante su éxito en cuanto a aumentar el intervalo de integración ha sido bastante limitado y tan sólo han aportado factores de aceleración de entre 2 y 5 (véase Michael A. Sherman et al., Method for large timesteps in molecular

modeling, US Patent 10053253), una mejora totalmente insuficiente dados los larguísimos tiempos de computación que, como se ha comentado, ofrece el estado del arte de la Dinámica Molecular.

La técnica más conocida basada en el aumento del intervalo de integración de las fuerzas lentas es conocida como MTS (Multiple- Timestep). En ellas los métodos tradicionales de integración en Dinámica Molecular (Verlet y sus variantes) son modificados evaluando las fuerzas lentas menos frecuentemente que las rápidas. El método MTS más popular y con diferencia más extendido en simuladores comerciales (GROMACS: H. J. C. Berendsen, D. van der Spoel and R. van Drunen, GROMACS: A message-passing parallel molecular dynamics implementation Computer Physics Communications Volume 91, Issues 1-3, 2 September 1995, Pages 43-56, AMBER: D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. DeBolt, D. Ferguson, G. Seibel, P. Kollman,

AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules, Computer Physics Communications Volume 91, Issues 1-3, 2 September 1995, Pages 1-41, CHARMM: B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, M. Karplus, CHARMM: A program for macromolecular energy, minimization, and dynamics calculations, Journal of Computational Chemistry Volume 4, Issue 2, pages 187– 217, Summer 1983, NAMD: J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD,. Journal of Computational Chemistry, 26:17811802, 2005,) es el conocido como r-Respa (M.E. Tuckerman, B.J.Berne, and G.J. Martyna., Reversible multiple time scale molecular dynamics. J.Chem.Phys., 97:1990-2001, 1992.), también llamado Verlet-I (véase H.Grubmüller, H.Heller, A.Windemuth, and K. Schulten. Generalized Verlet algorithm for efficient molecular dynamics simulations with long-range interactions., Mol. Sim.,6:121-142,1991.) y que fue propuesto de forma independiente en las dos anteriores referencias. La filosofía del método r-Respa es la habitual de los métodos MTS, es decir, la evaluación menos frecuente de las fuerzas rápidas, de modo que en cada una de estas evaluaciones se acumula la magnitud de la fuerza que no se ha ido teniendo en cuenta entre dichas evaluaciones, así estas fuerzas son actualizadas a impulsos. De nuevo, la mejora en el paso de integración de r-Respa es marginal, con unos tiempos típicos desde 0.5 fs, 2 fs, 4 fs (véase H. Grubmüller, P. Tavan, Multiple Time Step Algorithms for Molecular Dynamics Simulations o f Proteins: How Good Are They?, Journal of Computational Chemistry, Vol. 19, No. 13, 1998), hasta 1 fs, 2 fs y 5 fs (véase T. C. Bishop, R. D. Skeel, K. Schulten, Difficulties with Multiple Time Stepping and Fast Multipole Algorithm in Molecular Dynamics, Journal of Computational Chemistry, Vol. 18, No. 14, 1997) y por tanto muy lejos de lo que realmente se necesita en Dinámica Molecular. La razón principal de esta mejora marginal de r-Respa es que la actualización de las fuerzas lentas mediante impulsos tiene lugar de forma periódica, en particular a una frecuencia cercana a múltiplos de algunas de las frecuencias naturales del sistema. Esta circunstancia hace que el sistema entre en resonancia (Schlick, M. Mandziuk, R.D. Skeel, and K. Srinivas. Nonlinear resonance artifacts in molecular dynamics simulations. J.Comput. Phys., 139:1-29, 1998.) y por tanto que la simulación no sea viable. Para evitar estas resonancias se introdujo la familia de integradores MOLLY (véase B. Garcia-Archilla, J. M. Sanz- Serna & R. D. Skeel, Long-time-step methods for oscillatory differential equations, SIAM J. Sci. Comput. 20(1998), 930-963 y también J. M. Sanz-Serna, Mollified impulse methods for highly-oscillatory differential equations, SIAM J. Numer. Anal. 46(2008), 1040-1059). A pesar de ser una alternativa al uso de restricciones, el paso de tiempo sigue estando restringido a unos pocos femtosegundos, permitiendo solo un 50% de aumento en el timestep (J.A. Izaguirre, Q. Ma ,T. Matthey, J. Willcock, T. Slabach, B. Moore, G. Viamontes Overcoming instabilities in Verlet-I/r- RESPA with the mollified impulse method,).

Una alternativa es aumentar la masa de los átomos de hidrogeno, ya que debido a su baja masa sus enlaces con otros átomos, son los responsables de las más altas frecuencias de oscilación. Estos métodos permiten aumentar el paso de tiempo para las fuerzas rápidas hasta 4 fs sin alterar la dinámica, en comparación con los 2 fs de MOLLY, y alcanzan hasta los 7 fs a costa de eliminar las oscilaciones más rápidas del sistema, pero sin modificar las propiedades de equilibrio termodinámico del sistema (véase K.A. Feenstra, B. Hess, H. J. C. Berendsen, Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich Systems, Journal of Comput. Chem. Volume 20, Issue 8, pages 786–798, June 1999)

Veamos en más detalle el modelo que tenemos en la Dinámica Molecular. La unidad básica de este modelo es el átomo. Para un sistema de N átomos se introduce un potencial que modela las interacciones entre los átomos. La razón de introducir un potencial en vez de directamente las fuerzas es tanto física, ya que la energía se puede escribir como la energía cinética más dicho potencial, como matemática, ya que el potencial es una función escalar.

A partir del potencial se puede obtener fácilmente la fuerza sin más que aplicar el gradiente, así si r1 , r2 , ... , ri , …

, rN son las coordenadas de los N átomos,

la ecuación de Newton se escribe a partir del potencial UTOTAL como

d 2r B

mi 2 i =-UTOTAL (r1,r2,.......,rN ) (2)dt Bri

Dadas los diferentes tipos interacciones que ocurren, este potencial es la suma de diferentes contribuciones:

U =U +U +U +U +U (3)

TOTAL BOND ANGLE DIHEDRAL VDW COULOMB

Siendo sus expresiones analíticas más comunes las siguientes

bond 2

U =Ik (r -r0) (4)

BOND i ii

enlaces angle 2

U =k (B -B 0) (5)

ANGLE i ii

Iangulos dihed

U =I Ik [1+cos( n ¢ -r )] (6)

DIHEDRAL i i i i n

diedros

i

ij 12 ij 6

UVDW =II4Eij [(( ) -(( )] (7)

i

rr

j i iij ij

qq

U =IIji (8)

COULOMB ij i 4 E0rij

UBOND es el denominado potencial de distancia de enlace, siendo ri la distancia a la que están los átomos que bond

forman el i-ésimo enlace, ri0 es la distancia de equilibrio para dicho enlace, y ki es una constante que cuantifica

la fuerza de dicho enlace, en particular del potencial armónico utilizado. Este potencial, y la fuerza que se deriva de él, modelizan el hecho que átomos enlazados están a una cierta distancia. Por ejemplo, la distancia de dos átomos de carbono, unidos por un enlace simple, se sitúa en torno a 1,54 amstrongs, que sería la distancia de equilibrio para este modelo molecular UBOND es el denominado potencial ángulo de enlace, siendo Bi formado por 3 átomos enlazados, Bi0 es el ángulo

angle

de equilibrio para dicho ángulo, y ki es una constante que cuantifica magnitud de dicha fuerza, y que esta

relacionada con la fuerza y características de las fuerzas de enlace entre los 3 átomos que forman el ángulo. Modeliza el ángulo formado por átomos enlazados. Por ejemplo para la molécula de agua, el ángulo de equilibrio para el ángulo formado por el átomo de oxigeno y los de hidrogeno se sitúa en torno a 105º.

UDIHEDRAL es el denominado potencial de torsión o de ángulos diedros, siendo ¢i el ángulo que forman los planos

definidos por cuatro átomos enlazados, ri es ángulo de equilibrio para dicho ángulo de torsión, ni es la llamada dihed

multiplicidad, ya que dado el carácter periódico de esta fuerza pueden existir diversos ángulos de equilibrio y ki es una constante que cuantifica dicha fuerza. Este potencial modelizan la planaridad de ciertas moléculas. UVDW es el potencial que modela las interacciones de tipo Van der Waals entre átomos no enlazados. Su origen esta en las interacciones electromagnéticas entre átomos debida a la estructura atómica, en este caso la dependencia funcional es el llamado potencial de Lennard-Jones, dónde rij es la distancia entre los átomos i-ésimo

y j-ésimo, y Eij e (ij son dos constantes que cuantifican la intensidad de la interacción entre ese par de átomos, que tiene carácter atractivo a larga distancia, y repulsivo a muy corta distancia. UCOULOMB es el potencial de Coulomb, que modela las interacciones electrostáticas entre átomos no enlazados,

dónde rij es la distancia entre los átomos i-ésimo y j-ésimo, qi y qj son sus respectivas cargas efectivas, y E0 es

la constante que cuantifican la intensidad de la fuerza electrostática, llamada en este caso permitividad eléctrica del medio.

Las diversas constantes son obtenidas mediante medios semi-empíricos, las llamadas constantes del campo de fuerzas (véase J.W Ponder, D.A. Case, Force fields for protein simulations. Adv. Prot. Chem. 66 27-85 (2003)). Reflejan las propiedades observadas experimentalmente para las moléculas. Dado que el enlace atómico y la estructura atómica están regidos por la mecánica cuántica, cálculos basados en dinámica cuántica y el los datos experimentales son utilizados para obtener este modelo Newtoniano de la dinámica molecular. De esta manera, cada campo de fuerzas lleva tabuladas el conjunto de valores necesario para modelizar la molécula.

Existen otras posibilidades para la dependencia funcional de los distintos potenciales. (véase capítulo 8, Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006).

UBOND , U ANGLE , UDIHEDRAL modelan la estructura de la molécula debida al enlace atómico. Su efecto es el de dar

las distancias y los ángulos entre átomos enlazados su valor experimental. Dada la fuerza del enlace covalente se trata de fuerzas rápidas, y con oscilaciones de alta frecuencia. El hecho de que se calculen como un sumatorio sobre el número de enlaces, ángulos y ángulos diedros (torsiones), cuyo número el aproximadamente proporcional N, es lo que determina la dependencia O(N) del coste computacional.

Por su lado, UVDW y UCOULOMB modelan las interacciones electromagnéticas entre átomos no enlazados. El doble

sumatorio es el causante de la dependencia O (N2) del coste computacional. A estos potenciales hay que sumarles el efecto de las condiciones periódicas de frontera, que calculadas mediante

el método Particle Mesh Ewald anteriormente comentado, que aportan un potencial extra UEWALD .

Para la integración numérica correcta de la ecuación de Newton, que relaciona la variación de las coordenadas con la fuerza y la velocidad, el ritmo de variación de las fuerzas es crucial. Ello es debido a que cuanto mayor es la variación de la fuerza entre dos instantes de tiempo, mayor será el error numérico.

Un integrador numérico nos permite hallar el valor de las coordenadas de los átomos en un instante a partir de un instante anterior, siendo la separación de tiempos entre ambos instantes el paso de tiempo h.

Si llamamos qi al vector formado por las posiciones y velocidades de los átomos del sistema, tenemos

q = I (q ,h) (9)

ii+1

Siendo I (q , h) el integrador. La dependencia en las coordenadas es debida

i

esencialmente al cálculo de fuerzas, es decir

q = I (F(q),h) (10)

ii+1

Cada integrador viene caracterizado por un coste computacional debido al cálculo de las fuerzas CF y al cálculo de las nuevas coordenadas CC. En el caso de integradores explícitos, este último es despreciable. No así en el para integradores implícitos, dónde este coste puede ser superior. Para una buena integración es necesaria que la variación de la fuerzas en el intervalo h sea pequeña. Para el caso de un paso de tiempo variable, predecir esta variación es clave en la elección del paso de tiempo.

Para simular las condiciones de experimentales se introducen los llamados termostatos (temperatura constante) y barostatos (presión constante), condiciones que se dan en los sistemas reales. Típicamente son usados los métodos de Berendsen (Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. (1984). Molecular-Dynamics with Coupling to an External Bath, Journal of Chemical Physics 81 (8): 3684–3690 (1984)), Nose-Hoover (véase Nose, S, A unified formulation of the constant temperature molecular-dynamics methods, Journal of chemical physics 81 (1): 511–519 (1984) y también GJ Martyna, DJ Tobias, ML Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys 101(5), 1994) y especialmente los métodos basados en dinámica de Langevin (véase T. Schlick,. Molecular Modeling and Simulation. Springer. pp. 435–438 (2002) y también SE Feller, Y Zhang, RW Pastor and BR Brooks, Constant pressure molecular dynamics simulation: The Langevin piston method, J. Chem. Phys. 103(11), 1995). Un integrador que integra la dinámica de Langevin es el denominado BBK (véase A. Brunger, C. L. Brooks, M. Karplus, Stochastic boundary conditions for molecular dynamics simulations of ST2 water, Chem. Phys. Lett., 105 (1984), pp. 495–500.) La mencionada dinámica de Langevin añade un término estocástico a la ecuación de Newton, por lo que la energía deja de conservarse y la dinámica no es reversible pero simulando las propiedades termodinámicas de un sistema a temperatura constante. Además se introduce un término de frenado proporcional a la velocidad. Este término estabiliza la dinámica, y permite aumentar el paso de tiempo. Sin embargo para que la dinámica sea realista el factor de frenado no debe superar un valor, y el paso de tiempo que se alcanza no debe superar los 14 fs para las fuerzas lentas (véase J. A. Izaguirre, D. P. Catarello, J. M. Wozniak, R. D. Skeel, Langevin stabilization of molecular dynamics, Journal of Chem. Phys. Volume 114 number 5 1 (2001)) o incluso 16 fs (véase Q. Ma and J. A. Izaguirre, Targeted mollified impulse -- a multiscale stochastic integrator for long molecular dynamics simulations). Sin embargo estos métodos emplean un elevado término de frenado (dumping), lo que ralentiza la simulación, y por tanto implica más iteraciones para alcanzar el estado deseado, además de no alcanzar el paso de tiempo que podría permitir la dinámica de las fuerzas lentas.

A pesar de los recientes avances computacionales, que permiten simular milisegundos de una proteína en un tiempo razonable (véase D. E. Shaw, R. O. Dror, J. K. Salmon, J.P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M.

M. D., B. Batson, K. J. B., 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, S. Piana, Y. Shan, B. Towles, Millisecond-Scale Molecular Dynamics Simulations on Anton," Proceedings of the Conference on High Performance Computing, Networking, Storage and Analysis (SC09), New York, NY: ACM Press, 2009), estos están basados en mejoras únicamente computacionales. Sin embargo no logran superar la escala del femtosegundo para el paso de integración por lo que siguen necesitando ingentes recursos computacionales al alcance de muy pocos.

Descripción de la invención.

La presente invención resuelve la problemática anteriormente comentada mediante la materia definida en las adjuntas reivindicaciones independientes.

Esta invención presenta un método para simular de manera eficiente un sistema compuesto por átomos mediante Dinámica Molecular. Esta simulación esta basada en la integración de las leyes de Newton, para lo cual la elección del paso de tiempo requerido en el integrador numérico es crucial. Dada la complejidad de estos sistemas, la elección del paso de tiempo óptimo no es trivial.

Típicamente estás simulaciones son realizadas a un valor de paso de tiempo fijo, y suficientemente pequeño para garantizar una buena integración. Sin embargo, este valor no siempre es el máximo que podría utilizarse en la integración, con lo cual hay una pérdida de eficiencia. La presente invención se basa en usar un paso de tiempo variable que se adapte continuamente a la dinámica del sistema de átomos a simular, buscando el máximo valor posible sin comprometer la estabilidad de la simulación.

Para ello, en el método de la invención se calcula la potencia de cada átomo. La potencia de un átomo nos dice si un átomo va a ser acelerado o frenado, y en que magnitud. Esto nos da información sobre el estado dinámico del sistema, además de permitir predecir su comportamiento futuro.

Un primer efecto de conocer en todo momento el estado dinámico del sistema, es que podemos determinar un paso de tiempo que se adapta a la dinámica del mismo. Además el carácter predictivo de este criterio, clave para aumentar la ganancia, permite evitar los problemas debidos a átomos que van a ser excesivamente acelerados. Para ello basta reducir su velocidad de manera que su potencia entre en valores normales. Asimismo este criterio es aplicable a dinámicas que utilicen algoritmos como ligaduras, dinámica de Langevin o con diferentes pasos de tiempo para diferentes fuerzas (Multistep) dónde otros criterios no son fácilmente aplicables.

Por lo tanto un primer aspecto de la invención se refiere a un método implementado por ordenador para simular sistemas de átomos mediante dinámica molecular, que comprende el empleo de al menos un integrador numérico para integrar las ecuaciones de movimiento de los átomos de dicho sistema. El método además comprende controlar la potencia de al menos parte de los átomos del sistema a simular, mediante la modificación adaptativa del valor del paso de tiempo de dicho integrador numérico, manteniendo la potencia dentro de un rango de estabilidad de potencia dentro del cual se considera que la simulación del sistema es estable.

Mediante dicho ordenador se integran las ecuaciones de movimiento de los átomos en los tiempos marcados por el paso de tiempo elegido para el integrador, que gracias al control de potencia que mantiene la simulación estable, será máximo sin perder estabilidad, lo que produce el efecto de que el número de integraciones que tiene que calcular el ordenador se reduce y por lo tanto la simulación se completa en menor tiempo.

La invención logra el efecto o ventaja técnica de realizar simulaciones complejas de dinámica molecular con gran rapidez, pero utilizando ordenadores con una capacidad computacional convencional.

Otro aspecto de la invención se refiere a un equipo o sistema para simular sistemas de átomos mediante dinámica molecular, que comprende al menos un ordenador programado para integrar las ecuaciones de movimiento de los átomos del sistema mediante al menos un integrador numérico. Dicho al menos un ordenador está programado para controlar la potencia de al menos parte de los átomos del sistema a simular, modificando adaptativamente el valor del paso de tiempo de dicho integrador numérico, de forma que se mantiene la potencia dentro de un rango de estabilidad de potencia dentro del cual se considera que la simulación del sistema es estable.

El ordenador comprende medios de entrada de datos, medios de salida de datos, un procesador y un programa de ordenador almacenado en la memoria del ordenador. Preferentemente, dicho programa está adaptado para la selección de fármacos, y el sistema de átomos a simular se corresponde con una biomolécula, de modo que dicho

programa al ser ejecutado por dicho procesador, genera un dato de salida sobre si el fármaco satisface un criterio preestablecido.

El programa se lleva a cabo un bucle de control de potencia en el que se calcula la potencia de al menos parte de los átomos del sistema a simular, y se modifica adaptativamente el paso de tiempo de dicho integrador numérico en función de la potencia de los átomos calculada, de modo que si la potencia es menor que un umbral inferior del rango de estabilidad se aumenta el paso de tiempo del integrador, y de modo que si la potencia es mayor que un umbral superior del rango de estabilidad se reduce el paso de tiempo del integrador. Mediante dicho ordenador se integran las ecuaciones de movimiento de los átomos en los tiempos marcados por el paso de tiempo elegido en cada instante para el integrador. Por lo tanto, la invención logra optimizar el uso de la capacidad de procesamiento del ordenador, ya que reduce el número de integraciones que dicho ordenador tiene que calcular para completar la simulación, y reduciendo por lo tanto el tiempo de simulación.

La invención aplicada al campo del diseño de fármacos, posibilita ensayar un mayor volumen de compuestos químicos, y por lo tanto aumenta la posibilidad de encontrar compuestos útiles para el tratamiento de enfermedades.

La invención también se refiere a un programa de ordenador almacenado en un soporte informático, estando el programa adaptado para la selección de fármacos, de modo que cuando se ejecuta en un ordenador es capaz de implementar el método de la invención.

Descripción de las figuras

La figura 1 muestra un diagrama general de las fases de un método de simulación de dinámica molecular. El control de potencia se realiza en la actualización de coordenadas (posiciones y velocidades) del integrador.

La figura 2 muestra una gráfica sobre la dependencia de la potencia con el paso de tiempo, donde se muestran la potencia promedio y la máxima de los átomos para diversos valores de paso de tiempo.

La figura 3 muestra una gráfica donde se observa pi para las primeras 260000 iteraciones de la simulación de la

proteína ubiquitina, siendo el paso de integración 1 femtosegundo (1 fs).

La figura 4 muestra una gráfica sobre la evolución del paso de tiempo (en picosegundos) para un tipo de control con

2 -3

umbral 6000000 au · A · ps .

La figura 5 muestra una gráfica del comportamiento del paso de tiempo para 3 simulaciones con 3 criterios diferentes para el control del paso de tiempo (en pico segundos).En azul un control basado en la conservación del energía, criterio mencionado anteriormente. En verde el criterio es la alineación entre variación de posición y velocidad. En rojo el criterio aquí presentado. Se ha introducido un algoritmo de control de ligaduras para permitir un aumento del paso de tiempo.

La figura 6 muestra una gráfica sobre la combinación con algoritmos de ligaduras (WIGGLE) y repartición de la masa de los átomos de hidrógeno. En la gráfica se muestra el paso de tiempo en cada iteración en pico segundos.

Realización preferente de la invención

En la figura 1 muestra el esquema general de un simulador de dinámica molecular. El primer paso es leer las variables físicas del sistema que queremos simular. Entre ellas están las posiciones y velocidades iniciales de los átomos a simular, los parámetros del campo de fuerzas para dichos átomos, las condiciones a simular, por ejemplo temperatura y presión, y diversos parámetros de configuración para el sistema, como pueden ser el número de pasos a simular y diversos parámetros de configuración internos del programa.

Dichas variables deben ser llevadas a variables internas del programa. La siguiente parte es el núcleo del simulador, el integrador, que como hemos comentado, para cada paso de la simulación es el encargado de calcular las posiciones, velocidades y fuerzas de los átomos del sistema Es justamente en la parte del integrador donde está localizada la presente invención, en particular en la parte de actualización de coordenadas (posiciones y velocidades), modificando las velocidades o el paso de tiempo.

En una simulación estable, los valores de potencia tanto total como individual están dentro de un rango de estabilidad de valores que podemos considerar normal. De igual manera, si los valores de potencia están dentro de este rango de estabilidad, la simulación es estable.

Esto nos permite controlar la estabilidad de la simulación, que es algo a priori difícil, a partir del control de un parámetro, la potencia de los átomos.

La invención logra usar el máximo paso de tiempo posible, manteniendo siempre la simulación estable, que se puede conseguir como hemos mencionado, manteniendo la potencia dentro un rango de valores que podemos denominar rango de estabilidad.

Existen diferentes posibilidades para ello. En cada iteración i de la simulación, que corresponde a un instante de tiempo discreto, el estado de cada átomo

aa

a viene descrito por su posición ri y velocidad vi . A partir de las posiciones de cada átomo, podemos calcular

aa

la fuerza Fi y la aceleración ai que sobre él actúa.

a

Se define la potencia pi del átomo a en el instante i como el producto escalar de la velocidad de un átomo por la

fuerza a que esta sometido:

a aa

p =F ·v (11)

i ii

a

Definimos pi al mayor de todos los pi , es decir el valor de potencia del átomo de mayor potencia.

En una dinámica Newtoniana, la potencia es una medida de la variación de energía cinética para una partícula, y por tanto en el caso que nos ocupa, por un átomo. En el caso de potencia positiva un átomo gana energía cinética, mientras que si la potencia es negativa la pierde, el átomo es frenado. Dado que estamos realizando una integración numérica de la ecuación de Newton, es de esperar que una buena integración mantenga esta propiedad física.

aa

Podemos definir también la potencia total Pi como la suma de la potencia pi para todos los N átomos del

sistema:

a

Pi =Ip (12)

Ni

a

De igual manera definimos la variación de la potencia pi como la diferencia de potencia para un átomo entre dos

iteraciones consecutivas:

a aa

p=p -p (13)

i ii-1

y Pi como la variación de potencia total entre dos iteraciones

P=P -P (14)

i ii-1

a

Definimos pi como el mayor de todos los pi .

Para estas variables dinámicas del sistema podemos definir unos umbrales u , U , 0 y de tal manera que en

condiciones de buena simulación, las variables pi , Pi , pi , Pi tengan siempre valores menores que su

correspondientes umbrales.

Este umbral dependerá de dos características de la simulación:

-: el sistema a simular: tanto el campo de fuerzas utilizado como el tipo de átomos del sistema puede afectar a estos umbrales. Para las cantidades totales también afecta el tamaño del sistema -condiciones de la simulación: la temperatura, que está relacionada con la energía cinética influirá en los valores normales de potencia. El uso de ligaduras también afectará, al eliminar ciertas fuerzas de enlace. Asimismo podemos ser más exigentes en la calidad de la simulación, deseando una buena conservación de la energía, lo que impone que el paso de tiempo sea menor y los valores de potencia y por tanto los umbrales también lo sean.

La dependencia de la potencia con el paso de tiempo puede verse en la figura 2. Dónde se muestran la potencia promedio y la máxima de los átomos para diversos valores de paso de tiempo. El sistema utilizado es una molécula de la proteína ubiquitina solvatada en agua El integrador usado es Velocity Verlet y con el termostato tipo Berendsen.

2 -3

En el eje vertical tenemos representada la potencia en au · A · ps , Siendo au la unidad de masa atómica, A la unidad de distancia en amstrongs, y ps picosegundos. En el eje horizontal representamos el paso de tiempo, en femtosegundos (1 fs = 0.001 ps)

La manera de determinar estos parámetros para unas condiciones de simulación es obtener la estadística de las diferentes variables a partir de una corta simulación en dichas condiciones. Por inspección puede obtenerse un umbral a partir del máximo de dichas variables dinámicas.

El umbral se escogerá como el máximo más un margen del 10 % del valor de potencia máximo para ese sistema. Una vez determinado un umbral para un tipo de sistema (gases, fluidos, nano-dispositivos, biomoléculas) y en unas condiciones de temperatura y presión, este nos servirá para sistemas y condiciones similares para un mismo simulador.

Las primeras iteraciones de la dinámica pueden usarse como etapa de inicialización de los umbrales, automatizando el proceso.

En la figura 3 puede observarse pi para las primeras 260000 iteraciones de la simulación de la proteína ubiquitina,

siendo el paso de integración 1 femtosegundo (1 fs).

Se da el caso de que los átomos que presentan mayor potencia son los átomos de hidrógeno, y son generalmente ellos los que comprometen la estabilidad del sistema. Un recurso que ya hemos citado anteriormente es aumentar su masa, a expensas de los átomos pesados a los que están ligados.

Una vez determinados estos valores umbral, existen varias estrategias compatibles entre sí para que la potencia se mantenga en su rango normal, manteniendo la simulación estable.

1-Control del paso de tiempo h:

Si la potencia es menor que el umbral el sistema es estable, y es posible aumentar el paso de tiempo. Si es mayor entramos en la zona de inestabilidad y hay que reducir el paso de tiempo para tratar de recuperarla. Un algoritmo sencillo basado en este principio es el siguiente.

En cada iteración:

Si variable < umbral entonces h = a•h Si no, entonces h = b•h

Con a > 1 y b < 1, números reales positivos. Es conveniente elegir b < 1/a (condición conservadora), de manera que el sistema puede estabilizarse. Este tipo de control es complementado, cuando conviene, con vueltas al pasado:

En cada iteración:

Si variable < umbral entonces h = a•h Si no entonces

h = b•h

si variable >> umbral entonces nuevo estado i = estado i-k fin si fin si

Siendo k un número entero. En la práctica con el criterio aquí presentado, la vuelta al pasado sólo resulta necesaria, cuando se elije un umbral excesivamente elevado. Este tipo de algoritmos ha sido probado ya para la Dinámica Molecular, basado en un criterio diferente, la conservación de la energía (véase S. J. Stuart, J. M. Hicks, M. T. Mury, An Iterative Variable-Timestep Algorithm for Molecular Dynamics Simulations) Este control puede realizarse a partir de cualquiera de los pares (variable, umbral) basados en la potencia:

(: pi , u ) (15)

(: Pi ,U ) (16)

(: pi , 0 ) (17)

( Pi , ) (18)

De los 4 pares el aconsejado es ( pi , u ), al permitir un control más fino de la potencia, que redunda en mayor

ganancia.

A continuación vemos en la figura 4 la evolución del paso de tiempo (en picosegundos) para este tipo de control, con

2 -3

umbral 6000000 au · A · ps

Apreciar de nuevo la dependencia de la potencia y por tanto del umbral recomendable, con el paso de tiempo de la simulación.

Reseñar que una simulación a 2.5 femto segundos (0.0025 ps) sin control de potencia sería inestable, por lo que este método permite extender el paso de tiempo.

La ventaja de este control adaptativo basado en el control de la potencia es que permite más ganancia. En la figura 5 podemos ver como se comporta el paso de tiempo para 3 simulaciones con 3 criterios diferentes para el control del paso de tiempo (en pico segundos).

En azul un control basado en la conservación del energía, criterio mencionado anteriormente. En verde el criterio es la alineación entre variación de posición y velocidad. En rojo el criterio aquí presentado. Se ha introducido un algoritmo de control de ligaduras para permitir un aumento del paso de tiempo.

La poca variación del paso de tiempo, y la ganancia respecto al resto de criterios, es indicio que este es óptimo para el control del paso de tiempo basado en bucles de control.

2-Control de la potencia individual

Una manera de mantener la potencia de cada átomo dentro de los valores normales es re-escalar la velocidad del átomo si su potencia ha superado el valor umbral:

u

a aa

Si > u entonces = (19)

pi vi a vi

pi

aa

Siendo pi y vi la potencia y velocidad respectivamente del átomo a en la iteracíon i.

Esta corrección de la potencia de átomos problemáticos permite aumentar de manera artificial el paso de tiempo. Junto con el control adaptativo del paso de tiempo esto permite alcanzar valores aún mayores para el paso de tiempo.

Combinados con algoritmos de ligaduras (WIGGLE) y repartición de la masa de los átomos de hidrógeno, esto permite valores superiores a 18 femtosegundos, como puede observarse en la figura 6, que muestra el paso de tiempo en cada iteración (en pico segundos):

3-Particularización por fuerzas

Es posible particularizar estas variables para cada tipo de fuerza o grupo de fuerzas. Así la potencia debida al tipo de fuerza o al grupo de fuerzas k , en la iteración i , para el átomo a será:

a ,k a ,k a

p = F ·v (20)

ii

Siendo Fi a ,k donde el la suma de la fuerza total que actúa sobre el átomo a tipo de fuerza o al grupo de fuerzas

k .

De manera análoga es posible definir el resto de variables y de umbrales. Esto puede particularizarse para cada tipo de fuerza, por ejemplo Van der Waals, ángulos; o por grupos de fuerzas, por ejemplo ligadas, no ligadas. Esta separación por fuerzas nos permite controlar de manera independiente el paso de tiempo para una integración Multiple-Timestep (MTS), como es el caso/RESPA, dónde los diferentes pasos de tiempo son integrados a ritmo diferente, siguiendo un criterio idéntico al expuesto en el apartado 1.

Todos estos algoritmos permiten aumentar el paso de tiempo manteniendo la estabilidad de la simulación. Como hemos comentado, la elección del umbral va a depender de las condiciones de simulación. En unas condiciones muy conservadoras, sin uso de ligaduras ni re-particionado de masas de hidrógenos y usando como integrador BBK a 300 K de temperatura el paso de tiempo estable es 1.1 fs, para un sistema formado por una proteína similar a la ubiquitina, la DHFR (Dihydrofolate reductasa) solvatada en agua. Dado que no se usan ligaduras ni re-particionado de la masa de hidrógenos, son estos átomos los que causan los problemas en la integración. Para cada paso de la simulación, el control de potencia limita la velocidad de los átomos según la ecuación (19).

Con este control de potencia es posible alcanzar un paso de tiempo de 2.1 fs sin modificar la dinámica global, siendo la ganancia un 90% respecto a un BBK convencional en las misma condiciones.

2 -3

El umbral elegido en este caso fue 3100000 au · A · ps , lo que esta en concordancia con los valores que teníamos para la ubiquitina, sistema muy similar, a 1 fs Dado que el comportamiento dinámico de ciertos átomos puede ser muy diferente, es posible usar un umbral diferente para cada átomo

ua aa aa

Si pi > u entonces vi = vi (21)

pi a

Típicamente puede ocurrir que hay gran variación de masa de un átomo a otro, como ocurre con los hidrógenos, por lo que dicho umbral puede ser simplemente una función de la masa de cada átomo

a aa a um a

Si pi > um entonces vi = vi (22)

pi a

Los pasos de tiempo estables son menores que para el integrador anterior, velocity Verlet y termostato Berendsen, especialmente debido al termostato, ya que Berendsen hace un re-escalado de las velocidades, que por tanto las limita, pero frena la dinámica global. BBK hace uso por el contrario de un termostato tipo Langevin, que además reproduce más fielmente las condiciones termodinámicas que se quieren simular.

En el caso de querer usar un umbral por tipo de fuerza o tipo de átomo, la manera de determinarlo es análoga a la de determinar un único umbral. Se estudia el correspondiente comportamiento de la correspondiente potencia para una simulación a paso de tiempo muy conservador, y se toma el máximo más un margen.

Claims

REIVINDICACIONES

1. Método implementado por ordenador para simular sistemas de átomos mediante dinámica molecular, que comprende el empleo de al menos un integrador numérico para integrar las ecuaciones de movimiento de los átomos de dicho sistema, caracterizado porque además comprende controlar la potencia de al menos parte de los átomos del sistema a simular, modificando adaptativamente el valor del paso de tiempo de dicho integrador numérico, manteniendo la potencia dentro de un rango de estabilidad de potencia dentro del cual se considera que la simulación del sistema es estable, y mediante dicho ordenador integrar las ecuaciones de movimiento de los átomos al ritmo marcado por el paso de tiempo elegido para el integrador.
2.- Método según la reivindicación 1 caracterizado porque se lleva a cabo un bucle de control de potencia en el que se calcula la potencia de al menos parte de los átomos del sistema a simular, y se modifica adaptativamente el paso de tiempo de dicho integrador numérico en función de la potencia de los átomos calculada, de modo que si la potencia es menor que un umbral inferior del rango de estabilidad se aumenta el paso de tiempo del integrador, y de modo que si la potencia es mayor que un umbral superior del rango de estabilidad se reduce el paso de tiempo del integrador, y mediante dicho ordenador integrar las ecuaciones de movimiento de los átomos en los tiempos marcados por el paso de tiempo elegido en cada instante para el integrador.
3.- Método según la reivindicación 1 o 2 caracterizado porque además comprende reducir la velocidad de los átomos cuya potencia supere un umbral superior del rango de estabilidad, con objeto de mantener la estabilidad de la simulación.
4.- Método según la reivindicación 3 caracterizado porque la reducción de la velocidad de los átomos se realiza mediante el escalado de la velocidad del átomo, de modo que la nueva potencia del átomo se iguala a la del umbral superior del rango.
5.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque el umbral superior se define como cien veces el valor máximo de las variables dinámicas del sistema, para una simulación a paso de tiempo constante de 0.5 femtosegundos, y un umbral inferior se define como cien veces el valor mínimo de estas variables dinámicas para la misma simulación.
6.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque la potencia de los átomos se calcula como el producto escalar de la velocidad del átomo por la fuerza que actúa sobre él.
7.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque comprende el empleo de un integrador numérico que incluye algoritmos de ligaduras, seleccionados entre: Shake, Rattle, or Wiggle, que se aplican a variables que representan a las fuerzas y velocidades de los átomos del sistema, y donde dichos algoritmos de ligaduras se emplean para la imposición de restricciones a las componentes de alta frecuencia de las fuerzas de enlace del sistema de átomos, y la repartición de la masa de los átomos de hidrógeno.
8.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque se definen dos o más rangos de estabilidad, cada uno asociado a un átomo distinto del sistema.
9.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque la integración de las ecuaciones de movimiento de los átomos consiste en una integración con múltiples pasos de tiempo (MTS), y se controla de manera independiente los diferentes pasos de tiempos en base a la potencia calculada con las fuerzas y velocidades para cada paso de la integración.
10.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque el sistema de átomos a simular se corresponde con: un gas, un líquido, un nano-dispositivo, biomoléculas, o proteínas.
11.- Método según cualquiera de las reivindicaciones anteriores caracterizado porque el método está adaptado para la selección de fármacos, y porque el sistema de átomos a simular se corresponde con una biomolécula, y porque el método comprende introducir en dicho ordenador datos correspondientes al fármaco a simular, y una vez completada la simulación, generar datos referentes a la dinámica de las moléculas que comprenden: posiciones, velocidades, y energías.
12.- Sistema para simular sistemas de átomos mediante dinámica molecular, que comprende al menos un ordenador programado para integrar las ecuaciones de movimiento de los átomos del sistema mediante al menos un integrador numérico, caracterizado porque además dicho al menos un ordenador está programado para controlar la potencia de al menos parte de los átomos del sistema a simular, modificando adaptativamente el valor del paso de tiempo de dicho integrador numérico.
13.-Sistema según la reivindicación 12 caracterizado porque el ordenador comprende medios de entrada de datos, medios de salida de datos, un procesador y un programa de ordenador almacenado en la memoria del ordenador, donde dicho programa está adaptado para la selección de fármacos, y porque el sistema de átomos a simular se corresponde con una biomolécula, y porque dicho programa al ser ejecutado por dicho procesador, está adaptado a

5 generar datos referentes a la dinámica de las moléculas que comprenden: posiciones, velocidades, y energías.
14.- Sistema según la reivindicación 12 o 13 caracterizado porque el programa se lleva a cabo un bucle de control de potencia en el que se calcula la potencia de al menos parte de los átomos del sistema a simular, y se modifica adaptativamente el paso de tiempo de dicho integrador numérico en función de la potencia de los átomos calculada,

10 de modo que si la potencia es menor que un umbral inferior del rango de estabilidad se aumenta el paso de tiempo del integrador, y de modo que si la potencia es mayor que un umbral superior del rango de estabilidad se reduce el paso de tiempo del integrador, y mediante dicho ordenador integrar las ecuaciones de movimiento de los átomos en los tiempos marcados por el paso de tiempo elegido en cada instante para el integrador.

15 15.- Programa de ordenador adaptado para el estudio de fármacos, configurado de forma que cuando se ejecuta en un ordenador es capaz de implementar el método definido en cualquiera de las reivindicaciones 1 a 11.