Modelos efectivos para sistemas con una dinámica estocástica subyacente

FACULTAD DE FÍSICA DEPARTAMENTO DE FÍSICA ATÓMICA, MOLECULAR Y NUCLEAR Modelos efectivos para sistemas con una dinámica estocástica subyacente UNIVERSIDAD DE SEVILLA Esther Vega Estévez Curso 2022-2023 Doble Grado en Física y Matemáticas Tutorizado por Antonio Prados Montaño y Carlos Alberto Plata Ramos Resumen En gran variedad de sistemas físicos, como las biomoléculas o los sistemas de baja dimensionalidad como el grafeno, se ha observado experimentalmente la aparición de un fenómeno de metaestabilidad. Esto es, coexisten estados con propiedades elásticas distintas (plegado/desplegado en las biomoléculas, rizado/combado en el grafeno). La existencia de metaestabilidad puede explicarse como consecuencia del acoplamiento entre los modos elásticos del sistema con ciertos grados de libertad internos. En este trabajo revisamos un modelo sencillo propuesto recientemente, con un solo modo elástico y grados de libertad internos modelados como pseudoespines con interacciones a primeros vecinos. Encontramos que aparece metaestabilidad en un cierto rango de temperatura y fuerza aplicadas, lo que está de acuerdo—desde un punto de vista cualitativo—con los resultados experimentales en biomoléculas sencillas como las horquillas de ADN/ARN. Además, analizando la dinámica del sistema encontramos saltos estocásticos entre el estado estable y el metaestable a fuerza constante, así como un fenómeno de histéresis al variar la fuerza a un ritmo constante, de nuevo en acuerdo cualitativo con los experimentos. Los resultados analíticos han sido comparados con simulaciones numéricas construyendo códigos numéricos propios en MATLAB, encontrando siempre un buen acuerdo con la teoría desarrollada. 1 Abstract In a wide variety of physical systems, such as biomolecules or low dimensional systems like graphene, the occurrence of a metastability phenomenon has been observed experimentally. That is, states with different elastic properties coexist (folded/unfolded in biomolecules, rippled/buckled in graphene). The existence of metastability can be explained as a consequence of the coupling between the elastic modes of the system with certain internal degrees of freedom. In this work we review a recently proposed simple model with a single elastic mode and internal degrees of freedom modeled as pseudospins with first-neighbor interactions. We find that metastability appears over a certain range of applied temperature and force, which is in agreement—from a qualitative point of view—with experimental results on simple biomolecules such as DNA/RNA hairpins. Furthermore, analyzing the dynamics of the system we find stochastic jumps between the steady state and the metastable state at constant force, as well as a hysteresis phenomenon when varying the force at a constant rate, again in qualitative agreement with the experiments. Our analytical results have been compared with numerical simulations carried out using in-house built MATLAB codes, and a good agreement has been found in all cases. 2 Índice general 1. Introducción 4 1.1. Biofísica y materiales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2. Transicionesdefase ................................. 8 1.3. TeoríadeKramers.................................. 13 2. Modelo. Propiedades del equilibrio 16 3. Dinámica del modelo 26 4. Resultados numéricos 32 4.1. Equilibrio ....................................... 32 4.2. Dinámica ....................................... 35 4.2.1. Ciclodehistéresis.............................. 35 4.2.2. Tiempos de permanencia medios . . . . . . . . . . . . . . . . . . . . . 37 4.2.3. Exploración del diagrama de fases . . . . . . . . . . . . . . . . . . . . 39 5. Conclusiones 45 Anexos: Códigos de simulación 47 Bibliografía 51 3 Capítulo 1 Introducción 1.1. Biofísica y materiales La física es una ciencia que trata de describir y explicar los fenómenos naturales y, a partir de estas descripciones, hacer predicciones. Para ello, utilizamos el lenguaje matemático. Resulta interesante constatar que fenómenos naturales muy diferentes a primera vista son equivalentes matemáticamente. El uso de modelos sencillos hace posible el estudio de la naturaleza, a la vez que ayuda a construir un marco teórico que puede utilizarse en múltiples contextos, más allá de la motivación original del modelo. A lo largo de este trabajo desarrollaremos un modelo motivado por la biestabilidad observada en la respuesta elástica de biomoléculas y de ciertos materiales de baja dimensionalidad. El ADN y el ARN son biomoléculas complejas esenciales para la vida, ya que contienen la información genética de cada individuo. Están compuestas de bases nitrogenadas: guanina (G), adenina (A), citosina (C) y timina (T) (o uracilo (U) en el ARN). Estas bases se concatenan dando lugar a cadenas que, a su vez, se unen entre ellas siguiendo la regla de apareamiento C-G o A-T (A-U en el ARN); formándose así la doble hélice de ADN o de ARN. Estas cadenas están unidas entre ellas por enlaces débiles como son los puentes de hidrógeno. Otras biomoléculas complejas que van a ser objeto de interés en este trabajo son las proteínas. Estas están formadas por aminoácidos, que se unen en una cadena y se pliegan en una forma determinada. La forma que adopta cada proteína es muy importante, ya que es esencial para el desarrollo de su función característica. En las últimas décadas, gracias al desarrollo técnico de dispositivos que permiten medir fuerzas del orden del pN y longitudes del orden del nm, se han realizado múltiples experimentos con estas biomoléculas en los que se ha podido investigar con gran precisión su respuesta elástica al aplicar una fuerza sobre la misma [1–3]. En el caso de los ácidos nucleicos, en estos experimentos se aplica una fuerza en ambos extremos de lo que se denomina una horquilla de ADN o ARN, que es una doble hélice en forma de A, como se 4 CAPÍTULO 1. INTRODUCCIÓN 5 state system in an external field at finite temperature (solid line). From this analysis, P5ab’s unfolding free energy (⌬G)is193⫾6kJ mol ⫺1 (12,15). A second, independent, measure of P5ab’s ⌬Gis the average area under the reversible folding/unfolding plateau, which equals the potential of mean force of folding and yields a ⌬Gof 157 ⫾20 kJ mol ⫺1 . After correction for the free energy reduction of the unfolded state due to tethering (calculated to be 44 ⫾10 kJ mol ⫺1 )(12), these values compare well with the predicted ⌬Gof unfolding untethered P5ab calculated with the mfold free energy–minimization method (⌬G sigmoid ⫽ 149 ⫾16; ⌬G 具area典 ⫽113 ⫾30; ⌬G mfold ⫽147 kJ mol ⫺1 )(12,16). Several force-extension traces showed the molecule’s extension jumping between two values when the force was within ⬃1pNofthe unfolding plateau (Fig. 2A, left inset). We investigated this bi-stability by imposing a constant force on the molecule with feedback-sta- bilized optical tweezers capable of maintaining a preset force within ⫾0.05 pN by moving the beads closer or further apart. Then, the end-to- end distance of the P5ab hairpin hopped back and forth by ⬃18 nm, signaling the repeated folding and unfolding of a single RNA molecule. As in the pulling experiments, transitions between the two states were unresolvably fast (⬍10 ms) and without intermediates. By increasing the pre-set force, it was possible to tilt the foldedNunfolded equilibrium toward the unfolded state and thus directly to control the thermodynamics and kinetics of RNA folding in real time (Fig. 2C). As the force was increased, the molecule spent more time in the extended open form and less time in the short folded form. Whether hopping can be observed with a particular type of RNA depends on the time resolution of the instrument, its drift rate, and the kinetic barrier to folding/unfolding as determined by the potential energy surface of the molecule (12). In the instrument we used, hopping could be observed for rates between approximately 0.05 Hz and 20 Hz. A ratio of the average lifetimes of the molecule in the two states yields the equilibrium constant K(F) for folding/unfolding at that force (Fig. 2D). Linear extrapolation of K(F) to zero force, and correction for free energy reduction due to tethering (as above), yields a ⌬Gof 156 ⫾8kJmol ⫺1 , which coincides with the ⌬Gvalues obtained from stretching and the predicted value. Therefore, three different methods of measuring P5ab’s unfolding ⌬G give similar results: (i) the fit to the distribution of opening forces, (ii) the average area under the folding/unfolding plateau, and (iii) the ratio of folded and unfolded lifetimes. The sensitivity of RNA hopping to external force is determined by the force-depen- dent length difference between the unfolded and folded forms, ⌬x(F). In particular, an expression analogous to the van’t Hoff formula holds: dln K(F)/dF ⫽⌬x(F)/k B T(17). Indeed, the slope of the ln Kversus Fplot (Fig. 2D) multiplied by k B Tis 23 ⫾4 nm, and the ⌬x(F 1/2 ) value thus obtained is within experimental error of the value from the Fig. 1. (A) Sequence and secondary structure of the P5ab, P5abc⌬A, and P5abc RNAs. The five green dots represent magnesium ions that form bonds (green lines) with groups in the P5c helix and the A-rich bulge (3). (B) RNA molecules were attached between two 2-␮m beads with ⬃500–base pair RNA:DNA hybrid handles. Fig. 2. (A) Force-extension curves of the RNA-DNA handles without an insert (red) and with the P5ab RNA (black) in 10 mM Mg 2⫹ . Stretching and relaxing curves are superimposed. Inset, detail of force-extension trace showing hopping. Right inset, force-extension curves for the RNA hairpin without Mg 2⫹ .(B) Probability of opening versus force in Mg 2⫹ was obtained by summing a normalized histogram of hairpins opened versus force. Data are from 36 consecutive pulls of one molecule. Solid line, probability p(E) of a two-state system: p(E)⫽1(1⫹e E/k B T ). Best-fit (least squares) values, ⌬G(F 1/2 )⫽193 ⫾6kJmol ⫺1 ,⌬x⫽22 ⫾1nm(12). (C) Length versus time traces of the RNA hairpin at various constant forces in 10 mM Mg 2⫹ .(D) The logarithm of the equilibrium constant in Mg 2⫹ plotted as a function of force (error bar ⫽1 SD). (E) Detail of the stretching (blue) and relaxing (green) force-extension curves of the P5abc⌬A molecule taken at low and high loading rates in 10 mM Mg 2⫹ . REPORTS 27 APRIL 2001 VOL 292 SCIENCE www.sciencemag.org734 on April 11, 2011www.sciencemag.orgDownloaded from Figura 1.1: A. Estructura de distintos tipos de horquillas de ARN. B. Esquema del experimento realizado, sujetando las horquillas de ARN por sus extremos. Imagen tomada de la referencia [4]. ilustra en la figura 1.1. Para cada valor de fuerza se realizan medidas de la separación de las cadenas durante un tiempo determinado. Experimentos muy similares se han realizado con proteínas, en los que se aplica una fuerza que produce que la proteína se elongue e incluso se despliegue para pasar a ser una cadena monodimensional, obteniéndose resultados análogos. Al medir la elongación de las biomoléculas se ha encontrado que, en un cierto rango de valores de la fuerza (rango que puede depender de la temperatura), la molécula fluctúa entre dos estados de elongaciones distintas. Es decir, se observa metaestabilidad. Por otro lado, realizando ciclos de extensión y compresión variando la fuerza se ha observado que se produce una histéresis que depende del ritmo de cambio de la fuerza: a mayor ritmo, mayor histéresis [4]. Ilustramos estos comportamientos con la figura 1.2, tomada de la referencia [4]. En el panel superior (A), se muestra la existencia de metaestabilidad, donde se comprueba que en un cierto rango de fuerzas (aquí entre aproximadamente 14 y 14.6 pN) existen dos estados estables y el sistema fluctúa entre ambos. En el panel inferior (B), se reportan dos ciclos de histéresis. A la izquierda el ciclo es menor porque el ritmo al que varía la fuerza es lento, mientras que a la derecha es mayor, ya que se ha aumentado el ritmo de cambio. Por otro lado, algunos experimentos con materiales cuasi-bidimensionales como el grafeno dan lugar a resultados similares, en el sentido de la emergencia de biestabilidad de dos estados, denominados rizado y combado. Estos experimentos se realizan, en general, con el método de microscopía de efecto túnel (STM) [5], en el que la punta del microscopio se sitúa sobre la muestra y se aplica un voltaje entre estas dos. Este voltaje produce, por un lado, que la muestra se caliente y, por otro, una interacción electroestática entre CAPÍTULO 1. INTRODUCCIÓN 6 Figura 1.2: A. Elongación de una horquilla de ARN frente al tiempo para distintas fuerzas aplicadas. B. Ciclos de extensión (azul) y compresión (verde) de fuerzas de la molécula P5abc∆A (triple hélice de ARN) para la fuerza cambiando a ritmo bajo (izquierda) y alto (derecha). Imagen tomada de la referencia [4]. la punta y la muestra. Los resultados muestran, como vemos en la figura 1.3, que al aumentar el voltaje se da una transición del estado rizado al combado, que es reversible para intensidades bajas de la corriente de efecto túnel. Sin embargo, para corrientes altas la transición no es reversible y, al disminuir el voltaje, el grafeno se mantiene en el estado combado. En la figura 1.4 se representa el diagrama de fases que se obtiene en las referencias [6,7] al considerar un modelo con grados de libertad internos similar al que consideraremos en este trabajo. En las zonas I y II solo existe un estado estable (rizado y combado respectivamente), mientras que en las zonas IIIa y IIIb se da metaestabilidad entre los dos estados, siendo el combado el estable en IIIa y el metaestable en IIIb. A lo largo de este trabajo vamos a intentar dar una explicación cualitativa a los resultados experimentales antes descritos. Propondremos un modelo general que sea aplicable a las biomoléculas anteriormente mencionadas, así como a materiales que tengan propiedades similares. En los experimentos con biomoléculas el parámetro que se mantiene fijo es la elongación o, mediante dispositivos más elaborados, la fuerza. Por simplicidad del tratamiento que involucra, consideraremos que el parámetro de control es la fuerza mientras que estudiamos su efecto en la elongación. En el capítulo 2 introduciremos el modelo y analizaremos su estado de equilibrio, com- CAPÍTULO 1. INTRODUCCIÓN 7 Figura 1.3: Altura de la punta del microscopio en la lámina de grafeno frente al voltaje aplicado para distintos valores de la corriente de efecto túnel. Las curvas están ligeramente desplazadas verticalmente por claridad. Imagen tomada de la referencia [5]. Figura 1.4: Diagrama de fases del modelo para la elasticidad de sistemas de baja dimensión analizado en las referencias [6,7]. Se presentan las distintas fases del sistema en el plano (κ, θ), donde κes un parámetro que modela la interacción repulsiva entre electrones y θes la temperatura (ambos convenientemente adimensionalizados). En la zona I solo existe la fase rizada, en la II solo la combada y en la IIIa y IIIb coexisten la rizada y la combada como estados estables. En la zona IIIa el estado estable es el combado y el metaestable, el rizado, y viceversa en la IIIb. Imagen tomada de la referencia [6]. CAPÍTULO 1. INTRODUCCIÓN 14 Figura 1.6: Forma cualitativa de un potencial biestable U(x). Los puntos xRyxLcorresponden a mínimos locales del potencial biestable, siendo xRel punto de equilibrio estable—ya que es el mínimo absoluto de U(x)—y xLel punto de equilibrio metaestable. El punto xUcorresponde a un máximo local y es por tanto un punto de equilibrio inestable. CAPÍTULO 1. INTRODUCCIÓN 15 en cuenta, la podemos reemplazar por ∫∞ −∞ exp[U(xL) θ−∣U′′(xL)∣ 2θ(x′′ −xL)2]dx′′ =√2πθ U′′(xL)exp[−U(xL) θ].(1.24) Con estas aproximaciones, llegamos a que el tiempo medio de transición de xLaxRviene dado por τL=2π √U′′(xL)∣U′′(xU)∣exp[U(xU)−U(xL) θ].(1.25) Vemos que esta expresión sigue la forma de Arrhenius, que constituye el resultado principal de la teoría de Kramers [10]. Observamos que la exponencial en la ecuación (1.25) depende de la diferencia entre el máximo y el mínimo de potencial, es decir, de la barrera de potencial. Esto se debe a que el sistema fluctúa en torno a uno de los puntos de equilibrio hasta que, para una fluctuación suficientemente grande, adquiere una energía lo suficientemente alta como para traspasar la barrera de potencial. Entonces, el sistema cambia de estado y pasa a fluctuar en torno al otro punto de equilibrio estable. El factor ∣U′′(xU)∣−1/2indica que una barrera plana es más difícil de sobrepasar que una puntiaguda de la misma altura. Por otro lado, el factor U′′(xL)−1/2indica que en un pozo de fondo plano el sistema llega a los alrededores de la barrera menos veces, de forma que es más difícil de sobrepasar. Un análisis análogo se puede realizar para el tiempo de transición de xRaxL. De esta forma hemos obtenido una expresión para los tiempos de transición válida para sistemas sobreamortiguados. Un análisis similar puede realizarse en sistemas subamortiguados, pero es más complejo matemáticamente y no lo presentamos aquí. El resultado para el tiempo de transición sigue teniendo un factor tipo Arrhenius, pero el prefactor depende del coeficiente de fricción y adquiere una estructura más complicada [10]. Para los objetivos de este trabajo, lo esencial es que el tiempo de transición tiene siempre un factor tipo Arrhenius, y es por ello que hemos presentado la deducción de su expresión en el caso más sencillo—el límite sobreamortiguado. Capítulo 2 Modelo. Propiedades del equilibrio El modelo que estudiamos en este trabajo fue introducido originalmente en el contexto de la biofísica para describir las curvas fuerza-extensión de biomoléculas [11]. En particular, consideraremos una cadena monodimensional de longitud L0con Nsitios equiespaciados, cuyos estados se caracterizan por una variable que puede tomar dos valores. Por tanto, modelaremos estos espacios como pseudoespines de Ising con la variable σi=±1, con condiciones de contorno periódicas (es decir, σN+1=σ1), siendo σ={σ1, σ2,...,σN} la configuración de todos los espines.1Así, la distancia entre espines es d0=L0/(N−1). Consideraremos que los espines están en contacto con un baño térmico a temperatura T. Supongamos que la cadena se deforma, de modo que su longitud pasa a ser L=L0+∆. Asumimos que los espines siguen estando equiespaciados y, por tanto, tendremos que ahora la distancia entre ellos será d=d0+∆/(N−1). En estas condiciones, consideramos una energía potencial de la forma V(∆,σ)=1 2mω2∆2+J(∆)N ∑ i=1 σiσi+1,(2.1) donde J(∆)=J0−µ∆(2.2) es una función de ∆, con J0yµpositivas. Este potencial tiene un término macroscópico correspondiente a un potencial elástico, 1 2mω2∆2, y un término de energía procedente de los espines debido a la interacción a primeros vecinos. La constante de acoplamiento de los espines Jdepende linealmente de la elongación de la cadena y se iguala a J0para la longitud L0. Esta elección es razonable siempre que ∆≪L0. Hemos tomado la interacción entre espines como una a primeros vecinos para modelar de forma sencilla la interacción 1Hablamos de pseudoespines porque no son realmente momentos angulares intrínsecos, sino grados de libertad internos que supondremos por simplicidad descritos por una variable dicotómica. Para simplificar el lenguaje, de aquí en adelante emplearemos el término espín en lugar de pseudoespín. 16 CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 17 de corto alcance entre los grados de libertad internos que se da en muchos sistemas físicos, por ejemplo en ciertas biomoléculas, como los ácidos nucleicos o proteínas descritos en la introducción. Definimos ahora x=∆−J0 µ,(2.3) de forma que J(x=0)=0. De hecho, basándonos en esta variable xque acabamos de definir, somos capaces de clasificar la interacción entre espines fácilmente. Para x<0la interacción entre espines es ferromagnética, lo que favorece situaciones de menor elongación (estado plegado), y para x>0la interacción es antiferromagnética, por lo que se favorece situaciones de mayor elongación (estado desplegado). Por simplicidad, denominaremos, indistintamente, de ahora en adelante a xelongación reducida, simplemente elongación, o posición cuando no haya lugar a dudas. Escribiendo el potencial en función de x, V(x, σ)=1 2mω2x2+Fcx−µx ∑N i=1σiσi+1, Fc≡mω2J0 µ,(2.4) excepto por una constante aditiva irrelevante. El parámetro Fctiene dimensiones de fuerza y lo denominaremos fuerza crítica, por razones que quedarán claras más adelante. El último término de V(x, σ)se puede interpretar como una interacción a primeros vecinos con una constante de acoplamiento Jeff =µx. (2.5) Si ahora aplicamos una fuerza externa F, el potencial se modificará como V(x, σ)=1 2mω2x2−Hx −µx ∑N i=1σiσi+1, H ≡F−Fc.(2.6) La distribución de equilibrio para este sistema viene dada por la distribución canónica, Peq(x, p, σ)=1 Ze−βH(x,p,σ),(2.7) siendo β=(kBT)−1, con kBla constante de Boltzmann, Zla función de partición Z=∫+∞ −∞ dx ∫+∞ −∞ dp ∑ σ e−βH(x,p,σ)(2.8) yH(x, p, σ)el hamiltoniano del sistema. Podemos considerar que este hamiltoniano es H(x, p, σ)=H0(x, p)+Hint(x, σ),(2.9a) H0(x, p)=p2 2m+1 2mω2x2−Hx, (2.9b) CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 18 Hint(x, σ)=−Jeff(x)N ∑ i=1 σiσi+1.(2.9c) Aquí, H0(x, p)es el hamiltoniano del modo elástico sin acoplamiento de espines y Hint(x, σ) el de interacción entre el oscilador y los espines. Este último se corresponde a la interacción entre espines a primeros vecinos con la constante de acoplamiento Jeff definida en la ecuación (2.5). A partir de la ecuación (2.7) podemos obtener la probabilidad marginal Peq(x, p), que nos será útil para estudiar el comportamiento del modo elástico, lo que es nuestro objetivo principal. Esta probabilidad viene dada por Peq(x, p)=∑ σPeq(x, p, σ)=1 Ze−βH0(x,p)Z0(x),(2.10) donde Z0(x)=∑ σ e−βHint(x,σ)=[2cosh(Jeff kBT)]N=[2cosh(µx kBT)]N (2.11) es la función de partición del modelo de Ising en una dimensión con interacción a primeros vecinos, véanse por ejemplo las referencias [12, 13]. De esta forma, podemos escribir la probabilidad marginal Peq(x, p)como Peq(x, p)=1 Ze−β[T(p)+Veff(x)],(2.12) con T(p)la energía cinética del sistema y Veff(x)=1 2mω2x2−Hx −kBTln(Z0)=1 2mω2x2−Hx −NkBTlncosh(Jeff kBT),(2.13) salvo una constante. Veff(x)es el potencial efectivo total que actúa sobre el oscilador, que incluye un término debido al acoplamiento del grado de libertad elástico con los espines. A continuación vamos a redefinir las variables de forma que nuestras ecuaciones de movimiento sean adimensionales. Así podemos tratarlas de forma más fácil y cómoda, además de identificar los parámetros adimensionales que controlan el comportamiento del sistema. A partir de la ecuación (2.6), teniendo en cuenta que los σipueden tomar los valores ±1, vemos que el último término de la derecha de la ecuación (2.6) puede ser, a lo sumo, µN/m. Si tomamos este valor como un indicativo de la dimensión de este término, para que todos los términos de esta ecuación sean comparables, llegamos a que mω2[x]2=m[x] [t]2[x]=µN[x],(2.14) CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 19 de donde podemos deducir que las escalas características de xy de tson [x]=µN mω2,(2.15a) [t]=ω−1.(2.15b) Con estas escalas características podemos deducir que [p]=m[x] [t]=µN ω.(2.16) Podemos obtener también la escala característica de la energía, a partir de la energía cinética: [E]=[p]2 m=µ2N2 mω2.(2.17) Ahora, sabiendo que [E]=NkB[T]como podemos deducir de la ecuación (2.13), llegamos a que la unidad natural de temperatura del sistema es [T]=µ2N mω2kB≡Tc,(2.18) que de ahora en adelante llamaremos temperatura crítica. A partir de esta temperatura definimos la temperatura adimensional del sistema como θ=T/Tc. Este será un parámetro muy importante en futuros análisis porque, en el caso H=0, para θ>1(o T>Tc) tendremos que el único punto de equilibrio del sistema es x=0, mientras que para θ<1(o T<Tc) tendremos dos puntos de equilibrio igualmente probables, x=±x0, con x0>0. De esta forma, para θ=1se dará una transición de fase y, por ello, hemos llamado a Tctemperatura crítica. Por otro lado, como [F]=[E]/[x], llegamos a que la escala característica de la fuerza es [F]=µN. (2.19) Ahora, redefiniremos las variables para llegar a las ecuaciones adimensionales, de forma que para cualquier magnitud M, definimos M∗=M/[M]. Vemos que la fuerza crítica adimensional tiene una expresión sencilla en función de la temperatura crítica, F∗ c=Fc Nµ =mω2J0 Nµ2=J0 kBTc .(2.20) Con las nuevas variables definidas de esta forma, sustituyendo en (2.6), llegamos a que el CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 20 potencial adimensional sigue la expresión V∗(x∗,σ)=x∗2 2−H∗x∗−x∗ N N ∑ i=1 σiσi+1.(2.21) Sustituyendo estas mismas variables en (2.13) vemos que el potencial efectivo adimensional se puede escribir como V∗ eff(x∗)=1 2x∗2−H∗x∗−θlncosh(x∗ θ).(2.22) De la misma forma, sustituyendo en la ecuación (2.12) podemos obtener la probabilidad marginal adimensional. Teniendo en cuenta que la energía cinética adimensional es T∗(p∗)=p∗2 2,(2.23) llegamos a que la probabilidad marginal es P∗ eq(x∗, p∗)=1 Z∗exp{−N θ[1 2p∗2+V∗ eff(x∗)]}.(2.24) En posteriores análisis prescindiremos de los asteriscos para facilitar la notación, aunque, salvo que se explicite lo contrario, en todo caso estaremos haciendo alusión a las variables adimensionales. Con esto hemos sentado las bases del modelo con el que vamos a trabajar a partir de ahora. Analizaremos la situación de equilibrio y la dinámica de este sistema y comprobaremos que el comportamiento del mismo es consistente con el observado en las biomoléculas y los materiales de baja dimensionalidad que presentamos en la introducción. De esta manera, habremos conseguido describir, de modo cualitativo, con un modelo efectivo y sencillo el comportamiento de estos sistemas físicos. Vamos a comenzar por estudiar el estado de equilibrio del modelo que acabamos de describir. La situación de equilibrio es aquella que maximiza la probabilidad marginal Peq(x, p). Es decir, teniendo en cuenta la ecuación (2.24), el equilibrio sucede cuando Veff(x)alcance un mínimo, ya que este potencial efectivo corresponde a una especie de energía libre, como vimos anteriormente en la sección 1.3, y xaquí hace el papel de parámetro de orden. Así, derivando la ecuación (2.22), llegamos a que la xen el equilibrio debe verificar que xeq =H+tanh(xeq θ),(2.25) que es una ecuación para la posición de equilibrio del oscilador xeq en función de Hy de θ. La figura 2.1 presenta Hcomo función de xeq para una temperatura θ<1, en la que aparece metaestabilidad. En los próximos párrafos analizaremos la ecuación (2.25) de CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 21 Figura 2.1: Hfrente a xeq para θ=0.7. La línea continua azul e ha obtenido a partir de la ecuación (2.25). Para valores de Hcomprendidos entre −HbyHbhay tres valores para xeq, que están representados en la figura como xR,xUyxLen la línea discontinua. Tanto xRcomo xL representan un mínimo del potencial, mientras que xUrepresenta un máximo y es, por tanto, inestable. Como el valor de Hdestacado por una línea discontinua rosa verifica que H>0, el estado estable es xR>0y el metaestable, xL<0. Para H=0los estados de equilibrio ±x(0) eq son simétricos. Para valores de ∣H∣por encima de Hbhay un único punto de equilibrio. Este análisis cualitativo es el mismo para cualquier θ<1. modo exhaustivo, para entender el comportamiento que se observa en la figura. Si consideramos la situación H=0estaremos en las condiciones del modelo utilizado en la referencia [14]. Para este caso, tenemos que xeq −tanh(xeq θ)=0,(2.26) donde vemos que xeq =0es siempre una solución para cualquier valor de θ. Para cualquier θ>1esta es la única solución, que corresponde a un mínimo de Veff, por lo que es un equilibrio estable. Sin embargo, para θ≤1aparecen dos nuevas soluciones y xeq =0pasa a ser un máximo de Veff, es decir, un equilibrio inestable. Las dos nuevas soluciones que aparecen para θ≤1son simétricas y las denotamos por ±x(0) eq , tomando x(0) eq >0. Así CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 22 vemos que θ=1(o T=Tc) define un punto crítico, donde se da un cambio de fase en el que pasamos de tener un único estado de equilibrio, xeq =0a tener 3 posibles estados de equilibrio. Tomando el límite θ→1−podemos aproximar tanh(xeq θ)≃xeq θ−x3 eq 3θ3,(2.27) y obtenemos, utilizando la ecuación (2.26), que x(0) eq ≃√3(1−θ).(2.28) Vemos que la posición es el parámetro de orden de una transición de fase de tipo pitchfork, de forma similar a lo que vimos en la sección 1.2, concretamente en la figura 1.5. Por otro lado, tomando el límite θ→0+, llegamos a que tanh(xeq θ)≃1y obtenemos que x(0) eq ≃1, que es el valor máximo para xeq cuando H=0. En ese caso, para H=0, el potencial visto en (2.22) es una función simétrica de x, de forma que las dos posiciones de equilibrio estable ±x(0) eq son equiprobables, ya que corresponden a dos mínimos del potencial igualmente profundos. Si consideramos ahora el caso H≠0yθ<1, tenemos de nuevo que pueden existir (dependiendo del valor de H) tres extremos de Veff(x), uno de ellos máximo (xeq =0) y los otros dos mínimos. Sin embargo, ahora el término −Hx en la ecuación (2.22) rompe la simetría entre x>0yx<0, de forma que el potencial ya no es una función par de x. Esto hace que los mínimos de potencial no sean equiprobables y, por tanto, uno de ellos será estable mientras que el otro será metaestable. Viendo cómo afecta el término Hx a la ecuación (2.22) deducimos que se verán favorecidos los estados en los que Hx >0, ya que esto se traduce en una contribución negativa al potencial. De esta forma, para H>0, es decir, F>Fcesperamos encontrar el sistema en el estado desplegado x>0y si H<0 (o equivalentemente F<Fc), esperaremos encontrar el sistema en el estado plegado x<0. Es por esto que consideramos Fcla fuerza crítica, ya que marca la diferencia entre que se favorezca el estado plegado o desplegado. A partir de la ecuación (2.25) podemos obtener ∂Hxeq. Teniendo en cuenta que ∂Hxeq = (∂xeq H)−1, se llega a ∂Hxeq =[1−1 θsech2(xeq θ)]−1 .(2.29) Por otro lado, si derivamos dos veces Veff(x)con respecto a xutilizando la ecuación (2.22), obtenemos ∂2 xVeff∣x=xeq ∂Hxeq =1.(2.30) El potencial efectivo tiene un mínimo local, es decir, la posición de equilibrio correspon- CAPÍTULO 2. MODELO. PROPIEDADES DEL EQUILIBRIO 23 diente es estable, cuando ∂Hxeq >0, mientras que para ∂Hxeq <0tenemos un máximo del potencial efectivo. Para xeq =0, utilizando la ecuación (2.29), se tiene que ∂Hxeq =θ/(θ−1). Así podemos comprobar que para temperaturas por encima de la crítica (θ>1)xeq =0es un equilibrio estable, mientras que para temperaturas por debajo de la crítica (θ<1), es inestable. A partir de la ecuación (2.29) podemos estudiar la monotonía de xeq frente a H. Vemos que, para θ>1,xeq es creciente como función de H. Sin embargo, para θ<1, tenemos que Hno es monótona con respecto a xeq. De hecho, podemos ver dónde se dan máximos y mínimos de Hen función de xeq teniendo en cuenta la ecuación (2.29) y que ∂Hxeq =(∂xeq H)−1. De esta forma, obtenemos que ∂xeq H=1−1 θsech2(xeq θ),(2.31) de modo que podemos hallar los mínimos y los máximos de H(xeq)igualando la ecuación (2.31) a cero. Así, obtenemos que los extremos se dan, respectivamente, en los puntos ±xb que siguen la expresión cosh2(xb θ)=1 θ,(2.32) donde tomamos xb>0. El máximo y el mínimo de Hse dan, respectivamente, en ±Hb, con Hbcumpliendo Hb=xb−tanh(xb θ).(2.33) En la figura 2.1 vemos representada la función H(xeq). Como podemos observar, para H>HbyH<−Hbse tiene que Hes una función creciente de xeq. Sin embargo, para valores de Hentre −HbyHbencontramos que xeq es una función multivaluada de Hcon tres posibles valores, uno de ellos entre −xbyxb, que es inestable. Los otros dos valores que obtenemos son estables. El mínimo absoluto del potencial se da cuando Hxeq >0y para el valor de xeq hay un mínimo relativo del potencial, de forma que el sistema en esta posición se encuentra un estado metaestable. Así, vemos que encontramos una región de metaestabilidad o coexistencia cuando se cumple que θ<1y∣H∣<Hb, como se ha ilustrado en el diagrama de fases en la figura 2.2. Para estos valores, si H>0(zona CD) el estado estable es el desplegado y el metaestable, el plegado (mínimos global y local de Veff respectivamente). En el caso H<0(zona CP), la estabilidad relativa de dichos estados se invierte. Si θ>1(zona U) encontramos que hay un solo estado estable para el sistema, que verifica que Hxeq >0y varía de forma continua con (θ, H). En las zonas D y P (θ<1,H>HboH<−Hbrespectivamente) solo hay un estado estable, correspondiente al desplegado y al plegado respectivamente. La línea continua en la figura 2.2 define una transición de primer orden porque el estado de equilibrio global cambia discontinuamente CAPÍTULO 3. DINÁMICA 30 Recordemos que el potencial Veff viene definido por la ecuación (2.22). Observamos que la dependencia explícita con Hes lineal. Sin embargo, la posición de los puntos de equilibrio depende de Hcomo vemos en (2.25). Por tanto, el potencial efectivo tiene una dependencia no lineal con H. A continuación, realizaremos una serie de aproximaciones de forma que transformemos esta dependencia en una lineal y veremos en qué condiciones esta aproximación es razonable. Consideremos H=F−Fc≪1y hagamos el desarrollo de Taylor de xeq(H), utilizando la ecuación (2.25). Llegamos a que xR,L =±x(0) eq +χH +O(H2),(3.13) donde x(0) eq viene dada por la ecuación (2.26) (es decir, es la solución de equilibrio para H=0) y χes la “susceptibilidad” para H=0: χ=∂Hx(0) eq ∣H=0=⎛⎝1−1 θ+x(0)2 eq θ⎞⎠ −1 .(3.14) De esta forma, despreciando los términos cuadráticos en H, obtenemos que xR−xL=2x(0) eq . Por otro lado, la posición inestable varía como xU=−(1 θ−1)H+O(H2).(3.15) Aplicando las ecuaciones (2.22) y (3.13) podemos aproximar el potencial efectivo en xRy en xL, despreciando los términos de orden cuadrático en H, por Veff(xR,L)≃Veff (x(0) eq )∓Hx(0) eq .(3.16) Por otro lado, aplicando la ecuación (3.15) llegamos a que Veff(xU)=O(H2),(3.17) por lo que de ahora en adelante despreciaremos este término. Con estas aproximaciones, llegamos a que la barrera es BR,L ≃B(0)±NHx(0) eq .(3.18) Sustituyendo esto en la ecuación (3.10) llegamos a que τR,L ≃2π N√V′′ eff(xR,L)∣V′′ eff(xU)∣exp⎧⎪⎪⎨⎪⎪⎩±B(0)+NHx(0) eq θ⎫⎪⎪⎬⎪⎪⎭=τcexp⎧⎪⎪⎨⎪⎪⎩±Nx(0) eq θH⎫⎪⎪⎬⎪⎪⎭,(3.19) CAPÍTULO 3. DINÁMICA 31 donde τc=2π N√V′′ eff(xR,L)∣V′′ eff(xU)∣exp{±B(0) θ}.(3.20) Con esto, hemos llegado a una aproximación del tiempo de permanencia medio, que sigue una exponencial que depende linealmente de H. Recordemos que esta aproximación es válida únicamente cuando H≪1. Comprobaremos en el siguiente capítulo que esta aproximación es adecuada en un cierto rango de Hpequeño, viendo que los tiempos de permanencia medios en los pozos derecho e izquierdo en función de Hsiguen, aproximadamente, la ecuación (3.19). Capítulo 4 Resultados numéricos Hemos realizado una serie de simulaciones con el programa MATLAB para el modelo que hemos introducido. Los resultados de estas simulaciones nos permitirán comprobar que el sistema se comporta de modo consistente con los análisis que hemos realizado en los capítulos anteriores. Hemos tomado un número de espines N=500,α=4y una temperatura adimensional θ=0.9en todos los casos.1Para esta temperatura, según las ecuaciones (2.32) y (2.33), tenemos que xb=0.2947,Hb=0.0215 yx(0) eq =0.5254. Esperamos que para valores de ∣H∣por debajo de Hbhaya un estado estable y otro metaestable para el sistema, mientras que para valores por encima de Hbhaya un único estado estable. Comenzamos estudiando el sistema en el equilibrio, comprobando que se da el fenómeno de metaestabilidad para ∣H∣<Hby que se sigue la curva de equilibrio dada por la ecuación (2.25). Después pasamos a estudiar la dinámica del sistema, centrándonos primero en investigar el comportamiento del sistema cuando variamos Hde forma cuasiestática. A continuación, comprobamos cómo varían los tiempos de permanencia medios con Hy constatamos que conforman un proceso poissoniano, siguiendo una distribución exponencial a Hfijo. Por último, exploramos el diagrama de fases del sistema, realizando dos tipos de simulaciones: variando Hpara un θ<1fijo (en este caso, θ=0.9) y variando θpara un Hfijo. Esta última es la única de las simulaciones en la que tomaremos θ≠0.9. 4.1. Equilibrio En este apartado comenzamos estudiando la aparición de metaestabilidad en nuestro sistema. Para ver este comportamiento hemos tomado H=3×10−3<Hb, de forma que el estado estable será, según la ecuación (2.25), xR=0.5402 y el metaestable, xL=−0.5095. Como podemos comprobar en la figura 4.1, el sistema en equilibrio fluctúa entre el estado 1Hay un único caso en el que hemos variado θ, como detallamos más adelante en la sección 4.2.3. 32 CAPÍTULO 4. RESULTADOS NUMÉRICOS 33 Figura 4.1: Evolución de la posición del sistema en equilibrio frente al tiempo para H=3×10−3. Las líneas discontinuas señalan la posición de los mínimos del potencial. Se observa la situación de metaestabilidad, ya que xoscila entre la posición de equilibrio estable y la metaestable. Como H>0, en el estado de equilibrio estable se da que x>0. Se comprueba que el tiempo de permanencia es mayor en el estado estable que en el metaestable. CAPÍTULO 4. RESULTADOS NUMÉRICOS 34 Figura 4.2: Valores medios de xpara cada H. Los puntos rojos son los valores que hemos obtenido en la simulación y la línea continua es la función teórica H(xeq)que sigue la expresión dada por (2.25). Los puntos siguen en buena aproximación la curva teórica, aunque los que corresponden a equilibrios metaestables se alejan ligeramente debido a que el tiempo de permanencia medio en el estado metaestable es menor, de modo que la muestra a partir de la que hacemos la media tiene poca estadística. Aún así, la desviación con respecto a la curva teórica es razonable. estable y el metaestable, de acuerdo con nuestras predicciones teóricas. Vemos también que el tiempo de permanencia en el estado desplegado es mayor que en el estado plegado, lo que coincide con lo que esperábamos a partir de la ecuación (3.19). Hemos repetido esta simulación para distintos valores de Hvariando desde H=−0.2 hasta H=0.2, tomando más valores en el intervalo en que ∣H∣<Hb, ya que es el intervalo de mayor interés. De esta forma, hemos tomado los valores medios de la posición en equilibrio para cada valor de Hy, en caso de haber un estado metaestable, hemos medido también la posición media para el estado metaestable. En la figura 4.2 podemos ver que los valores que hemos obtenido siguen en buena aproximación la curva que esperábamos para Hen función de xeq. Observamos que existe una cierta desviación de la curva teórica en los valores de xque son metaestables. Esto puede deberse a que, al no ser el estado estable, el tiempo de permanencia medio del sistema en esos estados es muy corto, de forma que al realizar la media no tenemos una muestra lo suficientemente significativa de puntos. Una posible solución a esto es dejar al CAPÍTULO 4. RESULTADOS NUMÉRICOS 35 sistema evolucionar durante más tiempo, de forma que tengamos más puntos y podamos realizar una mejor medida. A pesar de la limitación de tiempo, vemos que los puntos medidos se ajustan razonablemente a la curva teórica. Por tanto, el sistema en equilibrio se comporta de acuerdo a lo que habíamos predicho en la ecuación (2.25). 4.2. Dinámica En esta sección pasamos a realizar simulaciones para estudiar la dinámica del sistema. En primer lugar observamos un fenómeno que podemos interpretar como un ciclo de histéresis, que ocurre debido a la existencia de un estado metaestable. A continuación, presentamos simulaciones que nos permiten obtener los tiempos de permanencia en el estado estable y el metaestable. De esta forma, podemos comparar estos resultados con nuestra predicción teórica y comprobar que las aproximaciones que hemos hecho en el anterior capítulo son adecuadas, para un cierto rango de H. Por último realizamos una exploración del diagrama de fases del sistema, comprobando cómo aparece y desaparece el fenómeno de metaestabilidad tanto en función de Hcomo de θ. 4.2.1. Ciclo de histéresis Vamos a comenzar por estudiar la curva de fuerza-extensión de este sistema, es decir, vamos a comenzar desde el equilibrio para un Hmin <0, que en nuestro caso es Hmin =−0.2<−Hb, de forma que comenzamos con un único estado estable, el plegado. En cada paso de tiempo ∆tvamos a ir aumentando Huna cantidad ∆Hhasta alcanzar Hmax =−Hmin =0.2. Tomaremos el ritmo de cambio ∆H/∆tde tal forma que Hvaríe tan lentamente como para que el sistema no se separe de la curva de equilibrio. Tras cada uno de estos pasos medimos la posición xy, a continuación, volvemos a variar H. Tras alcanzar Hmax realizaremos este mismo proceso a la inversa, disminuyendo Hen una cantidad ∆H cada intervalo de tiempo ∆t. Hemos realizado dos procesos como los descritos anteriormente, tomando dos valores distintos para el ritmo al que variamos H, mostrando los resultados en la figura 4.3. En ambos casos el intervalo de tiempo ha sido ∆t=1 Nα =5×10−4. En primer lugar, hemos tomado un ritmo de cambio rápido (dentro de las restricciones que nos hacen mantenernos en el equilibrio), en el que ∆H=1.33 ×10−7, de forma que el ritmo es ∆H/∆t=2.67 ×10−4. Con este ritmo nos mantenemos en la curva de equilibrio y a medida que vamos aumentando Hnos acercamos a la zona de metaestabilidad. Cuando alcanzamos H=−Hbentramos en la zona de metaestabilidad, pero el equilibrio estable sigue siendo el negativo siempre que H<0. En el momento en que Hcambia de signo, el estado estable pasa a ser aquel en el que x>0, pero partimos de un equilibrio metaestable. CAPÍTULO 4. RESULTADOS NUMÉRICOS 36 Figura 4.3: Variación de la posición de equilibrio del sistema al variar Hde forma cuasiestática. Las curvas rojas representan curvas en las que hemos aumentado Hun valor ∆Hcada intervalo ∆t=5×10−4, mientras que las azules punteadas representan curvas en las que hemos disminuido Hal mismo ritmo. Las curvas continuas se han realizado con un ritmo de cambio mayor, ∆H= 1.33 ×10−7, de forma que al sistema no le da tiempo a adquirir suficiente energía térmica para saltar la barrera cuando el estado estable pasa de ser el plegado al desplegado (o viceversa). En ambos casos, se recorre toda la curva de metaestabilidad antes de pasar al estado estable. Las curvas punteadas se han obtenido variando Ha un ritmo más lento ∆H=2.67 ×10−9, de forma que el sistema puede evolucionar hacia el estado estable cuando este pasa de ser el plegado al desplegado (o viceversa). El salto en este caso se produce muy cerca de H=0. Esto hace que el sistema tenga que superar una barrera de potencial para poder pasar al equilibrio estable, como vimos en el capítulo 3. Si el ritmo de cambio es suficientemente rápido, como es este primer caso, el sistema no tiene tiempo para acumular suficiente energía térmica como para saltar la barrera de potencial, de forma que permanece en el equilibrio metaestable hasta que la barrera es lo suficientemente pequeña o hasta que H>Hby desaparece la metaestabilidad. De la misma forma, al realizar el proceso inverso, vemos el mismo fenómeno, aunque esta vez el sistema se mantiene en el estado desplegado hasta que desaparece y el único estado estable es el plegado. Este fenómeno puede observarse con claridad en la figura 4.3. Las líneas rojas son en las que hemos ido aumentando H, mientras que las líneas azules son en las que hemos ido disminuyéndolo. Al aumentar Hobservamos como el sistema se CAPÍTULO 4. RESULTADOS NUMÉRICOS 37 mantiene en el equilibrio metaestable hasta que este desaparece (línea roja continua), al igual que al disminuirlo (línea azul continua). Hemos realizado una segunda simulación, esta vez para un ritmo más lento, en el que ∆H=2.67 ×10−9, de forma que el ritmo de cambio es ∆H/∆t=5.33 ×10−6. Para este ritmo, el sistema tiene suficiente tiempo para llegar al estado de equilibrio estable para cada valor de H. Este fenómeno podemos observarlo de nuevo en la figura 4.3, en las líneas punteadas roja y azul. En ambos casos el sistema pasa del estado plegado al desplegado (o viceversa) prácticamente inmediatamente después de alcanzar H=0. Si tomáramos un ritmo de cambio aún más lento podríamos observar varios saltos entre un estado y otro antes de que el sistema se mantuviera únicamente en el estado estable, como vimos en la sección anterior. 4.2.2. Tiempos de permanencia medios Procedemos ahora a comprobar que la linealización del potencial efectivo en función de Hque realizamos en el capítulo anterior es válida cuando H≪1. Para ello, hemos realizado simulaciones en las que hemos medido el tiempo de permanencia medio tanto en el estado plegado como en el desplegado, para distintos valores de H, con ∣H∣<3×10−3≪1. Si aplicamos logaritmo a los resultados que obtenemos, debemos observar que estos siguen aproximadamente una recta en función de H. Además, si tenemos en cuenta la ecuación (3.19), comprobaremos que las rectas de mejor ajuste deben tener como pendientes ±Nx(0) eq /θ≃±292. En la figura 4.4 hemos representado los logaritmos de los tiempos de permanencia medios en función de H(puntos para el τLy asteriscos para el τR) junto a sus rectas de mejor ajuste y al logaritmo del tiempo medio que obtenemos de la ecuación (3.10) antes de realizar la linealización. Hemos tomado todos estos datos sin tener en cuenta la ordenada en el origen, normalizando con el valor para H=0, ya que por otras aproximaciones realizadas (por ejemplo tomar el límite sobreamortiguado para llegar a la ecuación (1.22)) no coinciden perfectamente y nos centramos en estudiar las pendientes de las rectas y la validez de la aproximación lineal. Como podemos comprobar, la aproximación lineal es válida en estas escalas de H, ya que la curva teórica se asemeja mucho al ajuste lineal de los resultados obtenidos. Además, como podíamos esperar, las pendientes de las rectas son aproximadamente las mismas en valor absoluto. Como comprobación, en la tabla 4.1 comparamos estas pendientes con la pendiente predicha por la aproximación lineal, dada por la ecuación (3.19). Confirmamos que, tanto en el caso plegado como en el desplegado, las pendientes están muy cercanas a las esperadas teóricamente, con un error relativo menor del 10%. Pasamos ahora a estudiar la distribución de tiempos de permanencia del sistema, ya CAPÍTULO 4. RESULTADOS NUMÉRICOS 38 Figura 4.4: Logaritmo de los tiempos de permanencia medios frente a H. Las líneas continuas corresponden al logaritmo de los tiempos de permanencia medios que siguen la ecuación (3.10). Las líneas discontinuas son las rectas de mejor ajuste a los logaritmos de los tiempos de permanencia medios obtenidos en la simulación. El tiempo de permanencia en el estado plegado está representado por puntos y el del estado desplegado, por asteriscos. En todos los casos hemos omitido la ordenada en el origen. La linealización para estas escalas de Hes válida a la vista de esta comparación y hay un buen acuerdo entre la pendiente esperada y la que hemos obtenido con la recta de mejor ajuste. Las pendientes de ambas rectas son aproximadamente iguales y de signo contrario, como predice la ecuación (3.19). Tabla 4.1: Comparación entre las pendientes de las rectas de mejor ajuste y la obtenida en la aproximación lineal. Pendientes ±Nx(0) eq /θ ϵr Plegado −278 ±14 −292 5% Desplegado 275 ±8+292 6% CAPÍTULO 4. RESULTADOS NUMÉRICOS 39 que al ser un proceso estocástico el sistema no permanecerá siempre el mismo tiempo en cada estado. Vamos a comprobar a través de una simulación que la distribución de tiempos de permanencia sigue una distribución exponencial de la forma fR,L(t)=1 τR,L exp(−t τR,L ),(4.1) donde τR,L es el tiempo de permanencia medio en el estado desplegado o plegado, respectivamente. La forma exponencial es consecuencia de que, como estamos en el estado de equilibrio, la probabilidad de transición por unidad de tiempo no depende del instante de tiempo. Para hacer esta comprobación, hemos realizado una simulación en la que hemos tomado H=3×10−3y hemos dejado el sistema evolucionar un tiempo t=2×106, coleccionando los datos de instantes de tiempo para los que se daba una transición. Para estos valores de los parámetros del sistema, como el estado estable es el desplegado, sabemos que el tiempo de permanencia medio es mayor para el estado desplegado xRque para el plegado xL. De esta forma, conseguimos una muestra de tiempos de permanencia de longitud NT= 849 y, a partir de esta muestra, hemos realizado un histograma con los datos obtenidos. A partir de los datos de estos histogramas podemos obtener algunos valores fR,L(t)para ciertos valores de t. En la figura 4.5 hemos representado lnfR,L(t)frente a los valores de tiempo correspondientes y hemos calculado la recta de mejor ajuste. Comprobamos que, efectivamente, los puntos se adaptan en buena aproximación a una recta, de forma que confirmamos que la distribución de tiempos de permanencia es exponencial, conformando un proceso de Poisson. Además, observando las escalas de los ejes de abscisas en la figura 4.5, comprobamos que el tiempo de permanencia medio es menor en el estado plegado (gráfica inferior) que en el desplegado (gráfica superior) para nuestro caso, ya que hemos tomado H>0. Podemos ahora comparar las pendientes de estas rectas con las predichas por la ecuación (4.1), utilizando la medida del tiempo medio que hemos obtenido a partir de las simulaciones. La comparación de los valores obtenidos la tenemos en la tabla 4.2, donde comprobamos que las pendientes de las rectas de mejor ajuste y −1/τR,L son muy similares, con un error relativo que ronda el 10%, así que podemos concluir que los resultados son satisfactorios. 4.2.3. Exploración del diagrama de fases En esta sección vamos a presentar resultados numéricos para explorar el diagrama de fases del sistema, realizando simulaciones para distintos valores de Hy de θ. En la figura 4.6 hemos representado con una línea continua azul la variable ±Hbcomo función CAPÍTULO 5. CONCLUSIONES 46 las cercanías de la transición plegado-desplegado, H≪1. En el capítulo de resultados numéricos hemos comprobado que las predicciones teóricas deducidas en los capítulos anteriores se verifican con bastante buena aproximación en las simulaciones realizadas. En todos los casos y propiedades medidas hemos obtenido un resultado satisfactorio, que concuerda con la teoría anteriormente desarrollada. En particular, hemos verificado que se da metaestabilidad con el modelo propuesto, en el rango de valores de Hyθque habíamos predicho. También hemos hecho un estudio de los tiempos medios de permanencia en los estados de equilibrio en las simulaciones, encontrando un acuerdo razonable con la teoría, a pesar de las aproximaciones realizadas. Además, hemos reproducido el ciclo de histéresis experimental al subir y bajar la fuerza aplicada en nuestro sistema, en completa analogía con lo observado en la referencia [4]. Finalmente, podemos concluir que el modelo propuesto es satisfactorio, ya que explica de forma sencilla y cualitativa el comportamiento que queríamos describir. Posibles formas de mejorar el modelo serían considerar más grados de libertad en los espines o múltiples estados entre desplegado y plegado (como una cadena que se despliega poco a poco). Sería interesante analizar estos modelos más complejos en el futuro, para tratar de conseguir una descripción no solo cualitativa, sino también cuantitativa, de los resultados observados en los experimentos. Anexos: Códigos de simulación En este anexo incluimos el código de la simulación que da lugar a la figura 4.1. A partir de este código se pueden realizar modificaciones sencillas para realizar las otras simulaciones que hemos llevado a cabo en el capítulo 4. 1 2% Fijamos primero los parametros del sistema : N espines , H la fuerza ( constante ) , alfa la unidad de tiempo , theta la temperatura adimensional , tmax el tiempo de simulacion , delta el tiempo de un salto entre una medida de x y otra y x0 el valor de x para tiempo 0, alrededor del x de equilibrio para estos datos . 3 4N = 500; 5H = 3e -3; 6theta = 0.9; 7alfa = 4; 8delta = 1/( alfa*N); 9tmax = 25000; 10 M = tmax/delta; 11 x0 = 0.56; 12 gamma = @(y) tanh (2*y/theta ); 13 14 % Definimos un vector de valores para x, desde x0 hasta el valor de x en t = tmax , que iremos modificando a medida que obtenemos valores para x resolviendo la ecuacion diferencial , pero que comienza siendo nulo. 15 16 x = zeros (1,M +2); 17 x(1) = x0; 18 x(2) = x0; 47 CÓDIGOS DE SIMULACIÓN 48 19 20 % A continuacion , generamos un vector de pseudoespines aleatorios sigma : 21 22 sigma = randi ([0 1], 1, N); 23 24 for k = 1: length ( sigma ) 25 if sigma (k) == 1 26 sigma (k) = 1; 27 else 28 sigma (k) = -1; 29 end 30 end 31 32 % Creamos un bucle for en el que iremos modificando sigma de forma aleatoria , siguiendo la probabilidad de transicion por unidad de tiempo , a la vez que x para ver como varia. 33 % Para la resolucion de la ecuacion diferencial hemos utilizado el metodo de Verlet. 34 % Los pasos que daremos en este bucle for son M y damos pasos de tiempo de intervalo delta. 35 36 % Definimos la parte del hamiltoniano que depende de los pseudoespines como E0: 37 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 38 39 % Como hemos tomado una distribucion aleatoria para los pseudoespines y hemos tomado un valor inicial de x cercano al equilibrio pero no exactamente en el equilibrio , dejamos el sistema evolucionar durante un tiempo suficiente como para que alcance el equilibrio . En este caso , basta con dejar evolucionar al sistema 10 unidades de tiempo , tomando pasos de tiempo de delta : 40 41 for k = (1:(10/ delta )) 42 i = randi (N ,1) ; 43 r = rand (1 ,1); CÓDIGOS DE SIMULACIÓN 49 44 if i == 1 45 p1 = 0.5*(1 - 0.5* gamma (x(2) )* sigma (1) *( sigma (2) + sigma (N))); 46 if r <= p1 47 sigma (1) = - sigma (1); 48 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 49 end 50 elseif i == N 51 pN = 0.5*(1 - 0.5* gamma (x(2))*sigma (N)*( sigma (N -1) + sigma (1))); 52 if r <= pN 53 sigma (N) = - sigma (N); 54 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 55 end 56 else 57 pi = 0.5*(1 - 0.5* gamma (x(2))*sigma (i)*( sigma (i -1) + sigma (i+1))); 58 if r <= pi 59 sigma (i) = - sigma (i); 60 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 61 end 62 end 63 x(3) = 2*x (2) - x (1) + delta ^2*( -x(2) + H + E0); 64 x(1) = x(2); 65 x(2) = x(3); 66 end 67 68 % Ahora tenemos tanto los pseudoespines como los dos valores iniciales de x en el equilibrio , asi que podemos comenzar a guardar los valores de x que obtenemos para ver la evolucion del sistema : 69 70 for k = 1:M 71 i = randi (N ,1) ; 72 r = rand (1 ,1); 73 if i == 1 74 p1 = 0.5*(1 - 0.5* gamma (x(k+1))* sigma (1) *( sigma (2) + CÓDIGOS DE SIMULACIÓN 50 sigma (N))); 75 if r <= p1 76 sigma (1) = - sigma (1); 77 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 78 end 79 elseif i == N 80 pN = 0.5*(1 - 0.5* gamma (x(k +1))* sigma (N)*( sigma (N -1) + sigma (1) )); 81 if r <= pN 82 sigma (N) = - sigma (N); 83 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 84 end 85 else 86 pi = 0.5*(1 - 0.5* gamma (x(k +1))* sigma (i)*( sigma (i -1) + sigma (i+1))); 87 if r <= pi 88 sigma (i) = - sigma (i); 89 E0 = 1/N*[ sigma (N) sigma ]*[ sigma 0] '; 90 end 91 end 92 x(k +2) = 2*x(k+1) - x(k) + delta ^2*( -x(k+1) + H + E0); 93 end 94 95 % Definimos la variable t para poder representar una grafica de x frente a t y dibujamos esta grafica : 96 97 t = (0: delta :tmax); 98 99 figure (1) 100 plot (t, x(2:M +2)) 101 ax=gca; 102 ax. FontSize =15; 103 yticks((-1:0.25:1)) 104 yticklabels({'-1','','-0.5 ','','0','','0.5','','1'}) 105 xlabel('$t$','interpreter','latex') 106 ylabel('$x$','interpreter','latex') Bibliografía [1] C. Danilowicz, V. W. Coljee, C. Bouzigues, D. K. Lubensky, D. R. Nelson, and M. Prentiss, Proceedings of the National Academy of Sciences 100, 1694 (2003). [2] G. Hummer and A. Szabo, Proceedings of the National Academy of Sciences 107, 21441 (2010). [3] F. Ritort, Journal of Physics: Condensed Matter 18, R531 (2006). [4] J. Liphardt, B. Onoa, S. B. Smith, I. Tinoco, and C. Bustamante, Science 292, 733 (2001). [5] J. K. Schoelz, P. Xu, V. Meunier, P. Kumar, M. Neek-Amal, P. M. Thibado, and F. M. Peeters, Physical Review B 91, 045413 (2015). [6] M. Ruiz-García, L. L. Bonilla, and A. Prados, Physical Review E 96, 062147 (2017). [7] M. Ruiz-García, L. L. Bonilla, and A. Prados, Physical Review B 94, 205404 (2016). [8] L. D. Landau and E. M. Lifshitz, Statistical Physics Part I, (Elsevier, Amsterdam, 1980), 3rd edition. [9] N. G. Van Kampen, Stochastic Processes in Physics and Chemistry, (Elsevier, Amsterdan, 2007), 3rd edition. [10] P. Hänggi, P. Talkner, and M. Borkovec, Review of Modern Physics 62, 251 (1990). [11] A. Prados, A. Carpio, and L. L. Bonilla, Physical Review E 86, 021919 (2012). [12] H. A. Kramers and G. H. Wannier, Physical Review 60, 252 (1941). [13] R. P. Feynman, Statistical Mechanics: A Set Of Lectures, (CRC Press, New York, 1998). [14] A. Prados, L. L. Bonilla, and A. Carpio, Journal of Statistical Mechanics: Theory and Experiment P06016 (2010). [15] R. J. Glauber, Journal of Mathematical Physics 4, 294 (2004). 51