Autoorganización en sistemas compuestos: sincronización, Turing, fluctuaciones

ESCUELA DE DOCTORADO INTERNACIONAL DE LA USC David García Selfa Tesis doctoral Autoorganización en sistemas compuestos: sincronización, Turing, fluctuaciones Santiago de Compostela, 2022 Programa de doctorado en Ciencia de Materiales TESIS DE DOCTORADO AUTOORGANIZACIÓN EN SISTEMAS COMPUESTOS: SINCRONIZACIÓN, TURING, FLUCTUACIONES David García Selfa ESCUELA DE DOCTORADO INTERNACIONAL DE LA UNIVERSIDAD DE SANTIAGO DE COMPOSTELA PROGRAMA DE DOCTORADO EN CIENCIAS DE MATERIALES SANTIAGO DE COMPOSTELA AÑO 2022 AUTORIZACIÓN DEL DIRECTOR/TUTOR DE LA TESIS D./Dña. Alberto Pérez Muñuzuri En condición de: Tutor/a y director/a Título de la tesis: Autoorganización de sistemas compuestos: sincronización, Turing, fluctuaciones INFORMA: Que la presente tesis, se corresponde con el trabajo realizado por D/Dña David García Selfa, bajo mi dirección/tutorización, y a utorizo su presentación, considerando que reúne l os r equisitos exigidos en el R eglamento de Estudios de Doctorado de la USC, y que como director/tutor de esta no incurre en las causas de abstención establecidas e n la Ley 40/2015. En Santiago de Compostela, 10 de octubre de 2022 Firma electrónica DECLARACIÓN DEL AUTOR/A DE LA TESIS D./Dña. David García Selfa Título da tese: Autoorganización de sistemas compuestos: sincronización, Turing, fluctuaciones Presento mi tesis, siguiendo el procedimiento adecuado al Reglamento y declaro que: 1) La tesis abarca los resultados de la elaboración de mi trabajo. 2) De ser el caso, en la tesis se hace referencia a las colaboraciones que tuvo este trabajo. 3) Confirmo que la tesis no incurre en ningún tipo de plagio de otros autores ni de trabajos presentados por mí para la obtención de otros títulos. 4) La tesis es la versión definitiva presentada para su defensa y coincide la versión impresa con la presentada en formato electrónico. Y me comprometo a presentar el Compromiso Documental de Supervisión en el caso que el original no esté depositado en la Escuela. En Santiago de Compostela, 10 de octubre de 2022. Firma electrónica A Sandra y a mis padres David García Selfa IV Agradecimientos La realización de una tesis doctoral es un camino que, por mucho que ya se haya dicho, es duro. También es gratificante, y mucho, por el conocimiento que se adquiere, pero también por las personas que me han acompañado en esa empresa. Y uno no puede dejar de lado un merecido reconocimiento a las personas buenas y queridas que me he encontrado en este camino y a las que han estado conmigo desde hace mucho (o desde siempre). Empezando por las personas que me acompañan y guían en el mundo académico, debo agradecer que hayan hecho posible esta tesis a mi director/tutor Alberto Pérez Muñuzuri, que me dio la oportunidad de hacer la tesis y, no menos importante, su amistad y la de su familia. También a los otros investigadores principales del Grupo de Física No Lineal: Vicente Pérez Muñuzuri y Gonzalo Míguez. Y, como no, a los compañeros de fatigas (y alegrías) en esta aventura, desde los más veteranos a los más nuevos: Darío, Dani, Carlos, Jorge, Mariamo, Irma, Alberto, Damián, Martín, Sara, Álex, Ismael, Santi, Pablo... y seguro que me dejo algunos (que me perdonen). También desde el mundo académico, quienes han permitido mi modesta colaboración y con los que he trabajado y aprendido sin límites: Juan Pérez Mercader, admirado por mí desde hace mucho tiempo y con quien es un privilegio trabajar, Christian Bick, que me acogió en la Universidad de Exeter y me abrió su casa y su familia, David Simakov y Gourab Ghoshal (thank you very much!). Asimismo, a los geógrafos Ángel Miramontes y a José Balsa. No puedo dejar de agradecer tampoco, a mis compañeras y comV David García Selfa pañeros de trabajo en el CESGA: Andrés Gómez, Carlos Mouriño, Lorena Fernández, Aurelio Rodríguez y al resto de la plantilla. Da gusto trabajar con ellos. Yendo al plano familiar y más personal, debo estar agradecido a mi compañera, esposa y motivo principal de querer estar en este mundo y hacer esta tesis: Sandra. Sin ella nada es posible y con ella todo es posible. También a las mejores personas que conozco, mis padres, David y Loli, y al resto de mi familia: mi hermano Antonio, hermanas Irene y Mariló, sobrinas y sobrinos Adrián, Jose, Martín, Irene, Cecilia, Alejandra y Valentina, cuñados y cuñada. Por supuesto, al resto de mi familia. También es importante agradecer el haber contado, y contar, con mis amigas y amigos, aunque la mayoría estén (geográficamente) lejos: os echo de menos. Y, por fin, un recuerdo agradecido a los que nos han ido dejando y ya no están entre nosotros, salvo en la memoria (siempre). Gracias a todas y a todos... incluso a los que me he dejado (involuntariamente) en el tintero. VI Resumen Contexto del trabajo La naturaleza, de forma espontánea, se autoorganiza dando lugar a estructuras altamente complejas que son responsables de muchas funciones importantes para la vida. En muchos de estos sistemas la estructuración aparece como consecuencia de comportamientos cooperativos de sus partes. Existen varios mecanismos que son relevantes para estos procesos, englobando mecanismos de sincronización, mecanismo de Turing, etc. Por otro lado, el análisis de cómo interaccionan las partes implica considerar diferentes tipos de acoplamientos y de conectividades entre las diferentes partes, es decir, redes de conexiones que se pueden modelar mediante redes complejas. Esta tesis se enmarca dentro de este contexto y plantea el estudio sistemático de la formación de estructuras de relevancia en la naturaleza en sistemas constituidos por individualidades más pequeñas. Se considerarán diferentes tipos de interacción entre las diversas partes y la presencia de fluctuaciones tanto en los propios parámetros de cada parte como en la interacción con factores externos. Resumen del trabajo La memoria de esta tesis se distribuye de la siguiente manera: VII David García Selfa Capítulo 1. Introducción En este capítulo se exponen los antecedentes del trabajo: la complejidad, consecuencia de sistemas cuyas partes interactúan entre sí de forma que surgen comportamientos dinámicos colectivos que no tendrían lugar si se consideran las partes aisladas entre ellas, es decir, el comportamiento colectivo no viene dado únicamente por el comportamiento individual de sus partes, sino por la interacción entre ellas. Se introducen los conceptos básicos de la teoría de bifurcaciones y la teoría de redes complejas. También se presentan en este capítulo los sistemas de reacción difusión, como paradigma de sistema complejo, y de algunos de sus comportamientos complejos que llevan a la autoorganización. Por otra parte, se describe la reacción de Belousov-Zhabotinsky, que va a ser la reacción en la que se basan los sistemas complejos utilizados a lo largo de este trabajo, así como su modelo matemático básico: el oregonator. Por último, introducimos los comportamientos complejos fundamentalmente estudiados en la tesis: resonancia (a partir de fluctuaciones externas), sincronización y patrones de Turing. Capítulo 2. Métodos numéricos y experimentales Se detalla cómo se han construido los osciladores químicos que se estudian en la tesis, tanto experimentalmente (en el Capítulo 5), como teórica y numéricamente en la mayoría de capítulos (excepto en el Capítulo 3). La población de osciladores químicos se ha construido mediante cuentas de resina que permiten intercambiar especies químicas con un medio activo. Las superficies de las cuentas llevarán fijado el catalizador de una reacción de Belousov-Zhabotinsky, mientras que el medio activo será una solución de dicha reacción libre de catalizador. Por otra parte, se presentan las bases teóricas de las técnicas numéricas empleadas en el estudio de los diferentes sistemas tratados y sus comportamientos: análisis de estabilidad y factores de crecimiento para la inestabilidad de Turing descrito mediante redes complejas, los análisis de bifurcación y continuidad y, por último, el modelo de aproximación de fases para sistemas VIII Resumen sincronizados. Capítulo 3. Oregonator no isotermo Se presenta un nuevo modelo de oregonator en el cual se ha incorporado la temperatura como variable a partir de un balance energético en un reactor agitado de flujo continuo con una reacción de Belousov-Zhabotinsky con refrigeración y fuente externa de calor. En este sistema podremos ver nuevos comportamientos, hasta ahora no estudiados, gracias al efecto de poder modular y controlar el sistema mediante la temperatura. En concreto, vamos a estudiar en profundidad un efecto resonante que tiene lugar cuando hay fluctuaciones de la temperatura gracias a una fuente externa de calor, lo que permite establecer un mecanismo de control del sistema. Este modelo, sin embargo, abre la puerta a otros tipos de comportamientos complejos si introducimos otras fluctuaciones, tanto deterministas como estocásticas, en la fuente externa de calor y/o en el refrigerante del reactor. Capítulo 4. Inestabilidad de Turing en osciladores químicos acoplados a través de un medio activo Este capítulo inicia una serie de capítulos en los que usaremos una población de osciladores químicos en un medio activo tal como se ha descrito en el Capítulo 2. Aquí vamos a usar como herramienta matemática para describir el sistema un modelo de red compleja en el que tanto el medio activo como cada cuenta de resina viene representada por un nodo de la red. Como todas las cuentas interactúan con el medio activo, existirá un enlace (conexión) entre el nodo que representa el medio activo y todos los nodos que representan a las cuentas. Por otra parte, la difusión vendrá descrita a partir de la matriz laplaciana (versión discreta del operador laplaciano). Aunque en un sistema de este tipo, es decir, de osciladores químicos basados en la reacción de Belousov-Zhabotinsky, nunca se han podido detectar hasta ahora patrones de Turing, se ha sido capaz IX David García Selfa de demostrar su existencia teniendo en cuenta una descripción más detallada de la difusión a partir de términos de difusión cruzada. Capítulo 5. Sincronización de dos osciladores químicos diferentes en un medio activo Se hace un trabajo experimental en el que se estudia y analiza cómo son los estados de sincronización entre sólo dos osciladores químicos en un medio activo en función de la distancia entre ellos. Además, a diferencia de la mayoría de trabajos en los cuales los osciladores suelen ser idénticos o casi idénticos, en este caso se trabaja con dos osciladores químicos de distinto tamaño y con distinta carga de catalizador y, por tanto, mostrarán distinto periodo de oscilación. Se observará un interesante cambio de comportamiento en la sincronización conforme aumenta la distancia entre osciladores hasta la pérdida de sincronización para distancias grandes. Capítulo 6. Osciladores químicos acoplados a través de un medio oscilante Se da un paso más en la descripción de los distintos comportamientos observados en el sistema formado por una población de osciladores químicos en un medio activo. En este capítulo, gracias a considerar el papel fundamental que tiene el carácter oscilatorio del medio activo (si bien de distinta naturaleza que el resto de osciladores), se explican por primera vez, desde un punto de vista teórico y empleando teoría de bifurcaciones, estados de sincronización detectados anteriormente tanto mediante simulación numérica como mediante experimentación, especialmente la transición entre un estado de sincronización y otro bautizado como de súpersincronización. Así, usando el medio activo como un oscilador, se ha construido un modelo reducido del sistema que ha permitido reconstruir matemáticamente y explicar la dinámica anteriormente observada. Además, se ha usado el modelo de aproximación de fases para describir la discontinuidad en los periodos observada en X Resumen la transición mencionada entre estos estados así como para calcular dichos periodos. Capítulo 7. Dos reactores acoplados con retardo temporal Por último, se da un paso más en la complejidad aportando un nueva interacción: el acoplamiento, con retardo temporal, entre dos poblaciones de osciladores como los estudiados en capítulos anteriores. Ahora, mediante simulación numérica y el estudio de la dinámica del modelo, se ponen de manifiesto nuevos comportamientos y estados de sincronización en las poblaciones de osciladores químicos. Este trabajo da pie tanto a la explicación de fenómenos de sincronización entre distintas poblaciones observados en la naturaleza y abre la puerta al establecimiento de mecanismos de control de los estados de sincronización mediante la regulación de la densidad de las poblaciones y del flujo de medio activo entre las mismas. Capítulo 8. Conclusiones y perspectivas Se exponen las conclusiones generales obtenidas del trabajo realizado en esta tesis, así como las perspectivas y el trabajo futuro que se puede seguir. XI David García Selfa XII Resumo Contexto do traballo A natureza de forma espontánea autoorganízase dando lugar a estruturas altamente complexas que son responsables de moitas funcións importantes para a vida. En moitos destes sistemas a estruturación aparece como consecuencia de comportamentos cooperativos das partes. Existen varios mecanismos que son relevantes para estes procesos, englobando mecanismos de sincronización, mecanismo de Turing, etc. Doutra banda, a análise de como interaccionan as partes implica considerar diferentes tipos de axustes e de conectividades entre as diferentes partes, é dicir, redes de conexións que podense modelar mediante redes complexas. O presente proxecto de tese enmárcase dentro deste contexto e expón o estudo sistemático da formación de estruturas de relevancia na natureza en sistemas constituídos por individualidades máis pequenas. Consideraranse diferentes tipos de interacción entre as diversas partes e a presenza de fluctuacións tanto nos propios parámetros de cada parte como na interacción con factores externos. Resumo do traballo La memoria de esta tesis se distribuye de la siguiente manera: XIII David García Selfa a way that collective dynamic behaviors arise and this behaviors would not take place if the parts are considered isolated from each other, that is, the collective behavior is not given only by the individual behavior of its parts, but by the interaction between them. Diffusion-reaction systems are also presented in this chapter, as a complex system paradigm, and some of their complex behaviors that lead to self-organization. Basic concepts of bifurcation theory and complex network theory are introduced. On the other hand, the Belousov-Zhabotinsky reaction is described. This work is based on this reaction, as well as its basic model: the oregonator. Finally, we introduce the complex behaviors fundamentally studied in the thesis: resonance (from external fluctuations), synchronization and Turing patterns. Chapter 2. Experimental and numerical methods It is described how the chemical oscillators studied in the thesis have been built, both experimentally (in Chapter 5), and theoretically and numerically in most chapters (except in Chapter 3). The population of chemical oscillators has been built using resin beads that allow chemical species to be exchanged with an active medium. The surfaces of the beads will carry the catalyst of a Belousov-Zhabotinsky reaction fixed, while the active medium will be a solution of said reaction free of catalyst. On the other hand, the theoretical bases of the numerical techniques used in the study of the different treated systems and their behaviors are presented: stability analysis and growth factors for Turing instability described by means of complex networks, bifurcation and continuity analyses, and finally, the phase approximation model for synchronized systems. Capítulo 3. Non-isothermal Oregonator A new model of oregonator is presented in which temperature has been incorporated as a variable from an energy balance in a continuous flow stirred reactor with a Belousov-Zhabotinsky reaction XX Summary with cooling and external heat source. In this system we will be able to see new behaviors, hitherto unstudied, thanks to the effect of being able to modulate and control the system through temperature. Specifically, we are going to study in depth a resonant effect that occurs when there are fluctuations in temperature thanks to an external source of heat, which allows establishing a control mechanism for the system. This model, however, opens the door to other types of complex behavior if we introduce other fluctuations, both deterministic and stochastic, in the external heat source and/or in the coolant of the reactor. Chapter 4. Turing instability in chemical oscillators coupled via an active medium This chapter begins a series of chapters in which we will use a population of chemical oscillators in an active medium as described in Chapter 2. Here we will use as a mathematical tool to describe the system a complex network model in which both the active medium as each resin bead is represented by a network node. As all the accounts interact with the active medium, there will be a link (connection) between the node that represents the active medium and all the nodes that represent the accounts. On the other hand, diffusion will be described from the Laplacian matrix (a discrete version of the Laplacian operator). Although in such a system, that is, of chemical oscillators based on the Belousov-Zhabotinsky reaction, it has never been possible to detect Turing patterns up to now, it has been possible to demonstrate their existence by taking into account a more detailed description of diffusion from cross-diffusion terms. Chapter 5. Synchronization of two different chemical oscillators in an active medium An experimental work is done, in which the synchronization states between only two chemical oscillators in an active medium are studied and analyzed as a function of the distance between them. In addition, unlike most works in which the oscillators are usually XXI David García Selfa identical or almost identical, in this case we work with two chemical oscillators of different sizes and with different catalyst loads. An interesting change in synchronization behavior will be observed as the distance between oscillators increases until the loss of synchronization for large distances. Chapter 6. Chemical oscillators coupled via an oscillating medium A further step is taken in the description of the different behaviors observed in the system formed by a population of chemical oscillators in an active medium. In this chapter, thanks to considering the fundamental role of the oscillatory character of the active medium (although of a different nature than the rest of the oscillators), previously detected synchronization states are explained for the first time both by numerical simulation and by experimentation, especially the transition between a state of synchronization and another baptized as super-synchronization. Thus, using the active medium as an oscillator, a reduced model of the system has been built that has allowed mathematical reconstruction and explanation of the previously observed dynamics. In addition, the phase approximation model has been used to describe the discontinuity in the periods observed in the mentioned transition between these states as well as to calculate said periods. Chapter 7. Two delayed coupled oscillators Finally, a further step is taken in complexity by providing a new interaction: the coupling, with time delay, between two populations of oscillators like those studied in previous chapters. Now, through numerical simulation and the study of the model’s dynamics, new behaviors and synchronization states are revealed in the populations of chemical oscillators. This work gives rise both to the explanation of synchronization phenomena between different populations observed in nature and opens the door to the establishment of control mechanisms of the states of synchronization by regulating the XXII Summary density of the populations and the flow of active medium between them. Chapter 8. Conclusions and outlook The general conclusions obtained from the work carried out in this thesis are presented, as well as the perspectives and future work. XXIII David García Selfa XXIV Objetivos Presentemos brevemente los objetivos perseguidos en esta tesis. Como acabamos de describir en el resumen, la interacción entre las partes que forman un sistema complejo, así como las interacciones con el entorno (el ambiente), hacen que surjan nuevos e inesperados comportamientos en el sistema. Estos comportamientos se manifiestan cuando el sistema se lleva lejos del equilibrio termodinámico y da lugar a la autoorganización del mismo mostrando diferentes patrones espacio-temporales tales como patrones de Turing o sincronización. Así, en los sucesivos capítulos de esta tesis pretendemos estudiar un sistema formado por osciladores químicos basados en la reacción de Belousov-Zhabotinsky. Esta reacción es altamente no lineal y hace que la riqueza de comportamientos sea aún mayor. Además, no nos vamos a conformar con la reacción en sí y su interacción con el ambiente, sino que además vamos a establecer sistemas compuestos cuyas partes son estos osciladores, es decir, vamos a describir cómo se construyen y a estudiar estas colectividades de osciladores. E incluso, para terminar, usaremos un sistema compuesto por dos de estos colectivos acoplados entre sí teniendo en cuenta el retardo temporal entre ambos colectivos. En el estudio no sólo vamos a limitarnos a la descripción físicoquímica de los sistemas compuestos por osciladores químicos, sino que vamos a tratar de explicar y anticipar su comportamiento matemáticamente desde el punto de vista de la dinámica de sistemas. Para ello, aplicaremos técnicas y modelos que van desde el estudio de la dinámica, incluyendo estabilidad, modelos de aproximación de fases y análisis de series temporales, hasta topologías de redes XXV David García Selfa complejas. Un objetivo principal de este trabajo va a ser poner en relieve cómo cierta parte del sistema, el medio activo en el que están inmersos los osciladores, juega un papel fundamental en las estructuras espacio-temporales descritas y que no podrían ser explicadas sin tener en cuenta el papel de esta parte (el medio activo). El resultado de este estudio se relaciona con sistemas naturales, especialmente sistemas biológicos, que manifiestan estos tipos de autoorganización, y que pueden describirse mediante modelos similares de reacción-difusión. XXVI Lista de publicaciones derivadas de la tesis La elaboración de esta tesis ha derivado en la publicación de varios artículos científicos y en la preparación de otros. 1. García-Selfa, D., Muñuzuri, A. P., Pérez-Mercader, J., and Simakov, D. S. A. (2019). Resonant Behavior in a Periodically Forced Nonisothermal Oregonator. The Journal of Physical Chemistry A, 123(38):8083–8088. Factor de impacto (JCR): 2.600 (2019). Cuartil (JCR): Q2 (Physics, Atomic, Molecular and Chemical). 2. Mussa Juane, M., García-Selfa, D., and Muñuzuri, A. P. (2020). Turing instability in nonlinear chemical oscillators coupled via an active medium. Chaos, Solitons and Fractals, 133. Factor de impacto (JCR): 5.944 (2020). Cuartil (JCR): Q1 (Physics, Mathematical). 3. García-Selfa, D., Ghoshal, G., Bick, C., Pérez-Mercader, J., and Muñuzuri, A. P. (2021). Chemical oscillators synchronized via an active oscillating medium: Dynamics and phase approximation model. Chaos, Solitons and Fractals, 145. Factor de impacto (JCR): 9.922 (2021). Cuartil (JCR): Q1 (Physics, Mathematical). 4. Carballosa, A., Balsa-Barreiro, J., Garea, A., Garcia-Selfa, D., Miramontes, A., and Munuzuri, A. P. (2021). Risk evaluation at municipality level of a COVID-19 outbreak incorporating XXVII David García Selfa relevant geographic data: the study case of Galicia. Scientific reports, 11(1). Factor de impacto (JCR): 4.997 (2021). Cuartil (JCR): Q2 (Multidiscicplinary Sciences). 5. Paramés-Estévez, S., Carballosa, A., García-Selfa, D., and Muñuzuri, A. P. (2023). Artificial intelligence techniques used to extract relevant information from complex social networks. Entropy, 25, 507. Factor de impacto (JCR): 2.738 (2021) - No se dispone del índice del año actual. Cuartil (JCR): Q2 (Physics, Multidiscicplinary). 6. García-Selfa, D., Rey-Devesa, P., and Muñuzuri, A. P. Synchronization of two different chemical oscillators in an active medium. En preparación. 7. García-Selfa, D., Bick, C., Muñuzuri, A. P., and Pérez-Mercader, J. Delayed coupling of two populations of chemical oscillators. En preparación. En 1 y 3, soy el autor principal. En 2 soy coautor principal, junto a M. Mussa-Juane (dicha coautora no era doctora en la fecha de publicación, aunque sí lo es en la actualidad). En 4, mi contribución principal es el cálculo de los factores de crecimiento y el contenido del artículo no se emplea en esta tesis, si bien dichos cálculos derivan de contenido de esta tesis; A. Carballosa y A. Garea son coautores no doctores. En 5, he contribuido en el tratamiento de datos, análisis formal, investigación y validación; el contenido del artículo no se emplea en esta tesis, si bien algunas técnicas y conceptos teóricos empleados se desarrollan en esta tesis; S. Paramés-Estévez y A. Carballosa son coautores no doctores. Las autorizaciones editoriales correspondientes al uso del contenido de algunos de estos artículos (1, 2 y 3) se incluyen al final de la tesis y ningún contenido de estos artículos se ha empleado en otras tesis. Asimismo, la elaboración de la tesis ha dado lugar al registro de software destinado a la detección de esferas en imágenes usando XXVIII Lista de publicaciones derivadas de la tesis técnicas de inteligencia artificial, concretamente de aprendizaje automático. El mismo tiene número de asiento registral 03/2022/523, se titula “Sphere detection” y los autores y propietarios de la propiedad intelectual somos Miguel Cruces Fernández, Alberto Pérez Muñuzuri y David García Selfa. XXIX David García Selfa del equilibrio termodinámico: oscilaciones, excitabilidad, multiestabilidad, patrones espacio-temporales (patrones de Turing, ondas, sincronización, etc.) Un sistema de reacción-difusión es, en sentido estricto, un sistema que implica una reacción química, que transforma sus componentes de forma local, y una difusión, que desplaza espacialmente estos componentes [Kuramoto,1984;De Wit,2007; Epstein and Pojman,1998]. Y, aunque estos sistemas en origen y por definición describen un proceso físico-químico, sus modelos matemáticos se pueden extender a otros ámbitos tan dispares como la biología [Murray,1993], la sociolingüística [Vidal-Franco et al., 2019] o la epidemiología [Carballosa et al.,2021] entre otros. De hecho, en el artículo seminal de Alan Turing de 1952, propone un sistema químico de reacción-difusión química para explicar la morfogénesis (cómo surgen patrones biológicos durante el crecimiento) y hace el primer estudio de estabilidad del sistema [Turing,1952]. El modelo básico para describir un sistema de reacción-difusión sin considerar fuerzas externas ni fluctuaciones se puede formular mediante ∂x ∂t =f(x;α) + D∇2x;x∈Rn, α ∈Rm,(1.3.1) siendo xlas concentraciones de las especies químicas de las reacciones, αun conjunto de parámetros fenomenológicos (como las constantes de velocidad de las reacciones ki, i = 1, ..., n), f(x;α) las expresiones asociadas a la ley de acción de masas para las reacciones químicas y Dla matriz diagonal de difusión (sus elementos diagonales Dison los coeficientes de difusión, considerados aquí uniformes), con lo que el segundo término del lado derecho de la ecuación es una expresión asociada a la ley de conservación de la masa en el transporte de las especies químicas por difusión. En este modelo básico suponemos que no hay efectos cruzados de difusión y se cumple la ley de Fick. Los términos de reacción de las ecuaciones anteriores (f) son típicamente no lineales, pues conllevan productos de las concentraciones de las especies químicas y los parámetros, las constantes de velocidades de las reacciones, dependen de la temperatura (ley de 6 Capítulo 1. Introducción Arrhenius). Sin embargo, el término de difusión es típicamente lineal en una primera aproximación. Así, es en el término de reacción en el que recae en mayor parte la complejidad del sistema y estos comportamientos complejos se estudian, principalmente, a través del análisis de las bifurcaciones (ramificación de comportamientos en sistemas topológicamente equivalentes que tiene lugar al variar los parámetros del sistema partiendo un punto de equilibrio y alejándonos de éste), del análisis de estabilidad y otras técnicas de la dinámica no lineal. También es interesante señalar que la simulación numérica es fundamental para este tipo de estudios (y así lo haremos a lo largo de este trabajo) [Strogatz,2018;Nicolis,1995; Guckenheimer and Holmes,1983]. El comportamiento del sistema va a depender de si posee, o no, grados de libertad espaciales: Si no posee dependencia espacial, el sistema no es más que un sistema de ecuaciones diferenciales ordinarias acopladas (no lineal, eso sí) cuya solución ya puede presentar distintos comportamientos como oscilaciones, excitabilidad, multiestabilidad, etc. En cambio, si posee dependencia espacial, tenemos un sistema de ecuaciones diferenciales en derivadas parciales acoplado de tipo parabólico y es capaz de generar un amplio abanico de soluciones en las que encontramos numerosos tipos de patrones espacio-temporales cuyas escalas espaciales y temporales vienen dadas por los valores de ki(de dimensiones [T]−1) y los valores de Di(de dimensiones [L]2[T]−1): la escala temporal por k−1 iy la escala espacial por √Dk−1. Ejemplos de patrones que pueden generar estas soluciones son [De Wit,2007]: •Sincronización, caos espacio-temporal y agrupamiento (clustering): por la forma de acoplamiento espacial, debido a la importancia relativa de los términos de reacción y difusión, de osciladores localmente periódicos o caóticos. 7 David García Selfa •Patrones de Turing: patrones espaciales estacionarios surgidos de una inestabilidad por ruptura de la simetría de un estado homogéneo (Figuras 1.2 (b), (c) y (d)). •Frentes de ondas: por el acoplamiento espacial de dos estados estacionarios con diferente estabilidad (la propagación va del más al menos estable). Si se dan oscilaciones o excitabilidad locales surgen ondas más exóticas, como ondas espirales u ondas cilíndricas (Figuras 1.2 (a), (e) y (f)). •Patrones compuestos: surgen de la interacción entre distintos tipos de inestabilidad. 1.3.1. Estabilidad y oscilaciones Uno de los elementos con los que vamos a trabajar en esta tesis son los osciladores químicos y, para explicar el comportamiento básico de un oscilador químico, vamos a usar una reacción con un mecanismo relativamente sencillo e introducida por Lengyel et al. en 1990 [Lengyel et al.,1990]: la reacción CDIMA (Chlorine Dioxide-IodineMalonic-Acid). La introducción de este modelo permitió describir la primera reacción con la cual se obtuvieron estructuras de Turing (además de otros patrones espacio-temporales, como oscilaciones, si se toman valores adecuados de los parámetros). Partimos de un sistema de reacción-difusión (1.3.1) en el que, por estar bien agitado, eliminamos el término difusivo y nos quedamos sólo con la reacción dx dt =f(x;α) ; x∈R2.(1.3.2) En este caso, el sistema (convenientemente adimensionalizado con x1→x,x2→y) queda, para el modelo CDIMA,        dx dt =f1(x, y, a) = a−x−4xy 1 + x2 dy dt =f2(x, y, b) = bx 1−y 1 + x2.(1.3.3) 8 Capítulo 1. Introducción (a) (b) (c) (d) (e) (f) Figura 1.2: Algunos patrones espacio-temporales. (a) Estructuras de Hopf (simulación numérica de una reacción-difusión autocatalítica cúbica o CARD). (b) Patrones de Turing obtenidos experimentalmente en una reacción CDIMA. (c) y (d) Distintos patrones de Turing (simulación numérica de una reacción CARD). (e) y (f) Ondas obtenidas experimentalmente para una reacción de Belousov-Zhabotinsky. Imágenes cedidas por el Dr. Alberto P. Muñuzuri. 9 David García Selfa Para estudiar su comportamiento dinámico, representamos en el espacio de fases (o espacio de estados1) la isoclinas nulas del sistema, esto es,      dx dt = 0 →y=(a−x)(1 + x2) 4x dy dt = 0 →y= 1 + x2 .(1.3.4) El sistema tiene un punto de equilibrio en el corte de las dos isoclinas nulas, esto es, x∗=a 5,y∗= 1 + a 52. Haciendo un estudio de estabilidad se obtiene que, si el determinante y la traza del jacobiano del sistema en el punto de equilibrio son positivos, se puede aplicar el teorema de Poincaré-Bendixson [Strogatz,2018] y obtener un ciclo límite (trayectoria cerrada en el espacio de las fases), lo que implica que el sistema será oscilante. Esto se cumple para b < bc, con bc=3a 5−25 a(Figura 1.3 (a), (b)). Para b > bc, la trayectoria en el espacio de las fases es una espiral hacia el punto de equilibrio de forma que el sistema acaba en equilibrio (Figura 1.3 (c), (d)). Para este valor de bctenemos una bifurcación de Hopf (en este caso, supercrítica) [Strogatz,2018]. Con este modelo se pueden obtener más bifurcaciones, en este caso variando el parámetro a, de forma que podemos encontrar que el sistema es biestable y se obtienen ondas viajeras u ondas de transición. Sin embargo, este sencillo ejemplo ya cumple su propósito de ilustrar una bifurcación. 1.4. La reacción de Belousov-Zhabotinsky: el “oregonator” Introducimos, ahora, la reacción con la que vamos a trabajar fundamentalmente en esta tesis. Esta reacción (Belousov-Zhabotinsky) y su modelo matemático (el oregonator) nos van a permitir estudiar 1A lo largo de esta tesis, especialmente cuando hablemos de sincronización y de la fase de un oscilador, usaremos con frecuencia el término equivalente de espacio de estados para evitar confusiones. 10 Capítulo 1. Introducción 0246810 x 0 5 10 15 y (a) 0 20 40 60 80 100 t 0 1 2 3 4 5 x (b) 0246810 x 0 5 10 15 y (c) 0 20 40 60 80 100 t 0 1 2 3 4 5 x (d) Figura 1.3: Simulaciones numéricas para la reacción CDIMA. Para un valor del parámetro a= 10, la bifurcación de Hopf tiene lugar para bc= 3.5. (a) Trayectoria en el espacio de las fases y (b) evolución temporal de xpara b<bc, esto es, ciclo límite y oscilaciones. (c) Trayectoria en el espacio de las fases y (d) evolución temporal de xpara b > bc, esto es, caída al estado de equilibrio. Nótese que las isoclinas nulas no cambian, en este caso, al variar el parámetro b, pues dependen exclusivamente del parámetro ay éste no cambia en esta bifurcación. 11 David García Selfa comportamientos más ricos y complejos que la introducida anteriormente. Además, se basa en remedar artificialmente un proceso biológico clave. En todos los organismos aeróbicos existe un mecanismo llamado ciclo de Krebs (o ciclo del ácido cítrico o ciclo de los ácidos tricarboxílicos) que es clave en la llamada ruta metabólica, que forma parte de la respiración celular, en la que se genera energía, además de producir compuestos esenciales para la vida como algunos aminoácidos (Figura 1.4). Figura 1.4: Ciclo de Krebs [https://es.wikipedia.org/wiki/Ciclo_de_ Krebs]. En 1951, mientras Boris P. Belousov buscaba un equivalente inorgánico del ciclo de Krebs, llegó a observar oscilaciones químicas en una solución acuosa de ácido cítrico, bromato acidulado y sulfato de cerio: la solución oscilaba cambiando de incolora a color amarillo. Cuando intentó publicar su trabajo, éste fue rechazado por considerarse que violaba las leyes de la termodinámica: la creencia, entonces, era que las reacciones químicas siempre evolucionaban de forma monótona hasta el equilibrio [Epstein and Pojman,1998]. No 12 Capítulo 1. Introducción fue hasta 1959, ocho años después, cuando se publicaron en las actas de un congreso [Belousov,1959]. En 1964, Anatol Zhabotinsky fue capaz de reproducir el experimento sustituyendo el ácido cítrico por ácido malónico y explicando que el cambio de color se debía a las oscilaciones en la concentración de Ce4+ [Zhabotinsky,1964]. Durante mucho tiempo se asumió que las oscilaciones en sistemas químicos cerrados homogéneos eran imposibles, como se comprobaba con un balance energético exhaustivo cerca del equilibrio termodinámico. Sin embargo, con la reacción de Belousov-Zhabotinsky (BZ) y su complejo mecanismo, quedó probado que las oscilaciones químicas eran muy probables lejos del equilibrio termodinámico. Figura 1.5: Mecanismo FKN de la reacción de Belousov-Zhabotinsky [Field et al., 1972]. En 1972, Field, Körös y Noyes describieron el mecanismo de la reacción BZ con 10 procesos claves a partir de detallados análisis termodinámicos y cinéticos (Figura 1.5) [Field et al.,1972]: (R1) - (R3): El ion bromato (BrO− 3) se reduce a bromo (Br2) vía ácido bromoso (HBrO2) e hipobromoso (HBrO) debido a la acción reductora del ion bromuro (Br−), cuya concentración decae. 13 David García Selfa (R4) - (R6): Cuando la concentración el ion bromuro llega a estar por debajo de un valor crítico ([Br−]<[Br−]cr), el ácido bromoso (HBrO2) comienza a competir con el ion bromuro (Br−) para reducir el ion bromato (BrO− 3). La producción autocatalítica de ácido bromoso (HBrO2) viene acompañada de la oxidación del ion metálico catalítico (Ce3+ →Ce4+) produciendo el repentino cambio de color. Además, se frena el crecimiento exponencial de la concentración de ácido bromoso (HBrO2). (R7) - (R10): Se cierra el ciclo retroalimentado reduciendo el catalizador (Ce4+ →Ce3+) y produciendo iones bromuro (Br−) cambiando, ahora lentamente, el color. De nuevo, Field y Noyes simplificaron el modelo FKN en 1974 en la Universidad de Oregón, que fue bautizado como oregonator [Field and Noyes,1974]. Este modelo se basa en la cinética de las reacciones y en una aproximación al estado estacionario (algunos componentes de la reacción apenas cambian tras unos cuantos ciclos, mientras que otros lo hacen muy rápidamente. Así, el mecanismo central de la reacción queda reducido a 5 pasos (Figura 1.6)). La reacción mostrada en la Figura 1.6 presenta los siguientes 5 pasos fundamentales del oregonator: (P1) A+Y→P+X (P2) X+Y→2P (P3) A+X→2X+ 2Z (P4) 2X→A+P (P5) B+Z→fY (1.4.1) (P1): producción del activador. (P2): inhibición. (P3): producción autocatalítica del activador. (P4): decaimiento del activador. 14 Capítulo 1. Introducción BrO3 − HBrO2 Ce4+ 2BrO2, Ce3+ CH2(COOH)2Br− 2HOBr [Br−]<[Br−]cr [Br−] > [Br−]cr ciclo auto catalítico par redox oxidación de sustrato orgánico inhibición del ciclo auto catalítico Figura 1.6: El oregonator. Esquema del modelo simplificado del mecanismo FKN [Field and Noyes, 1974]. (P5): producción del inhibidor. Donde tenemos los químicos de las reacciones cuyas concentraciones serán la variables del modelo: Activador: X - HBrO2. Inhibidor: Y - Br−. Catalizador metálico: Z - Ce4+. Y los que proporcionarán los parámetros del modelo: Reactivos: •A - BrO− 3. •Sustrato orgánico: BCH2(COOH)2. •Factor estequimétrico:f. Producto: P - HOBr. 15 David García Selfa Si la red es dirigida, se definen el grado de entrada y el grado de salida de un nodo, cuya suma es el número de conexiones total del mismo, esto es, el grado. También es interesante trabajar con el grado medio o conectividad media de la red, ⟨k⟩, que es el promedio de los grados de sus nodos. Asimismo, se define la centralidad de un nodo a partir de su número de conexiones: el nodo de mayor grado es el más central. Otro concepto importante es la longitud de camino mínima media: si definimos la distancia entre nodos lij como el número de enlaces que conecta iyjpor el camino más corto, la longitud de camino mínima media se define como ⟨l⟩=1 N(N−1) X i=j lij . También es importante el agrupamiento (clustering) de un nodo, C(i), que se define como el número de conexiones mutuas entre vecinos próximos, esto es, el cociente entre el número de enlaces que realmente tienen los vecinos entre sí, niy el número máximo que podría existir: C(i) = ni ki(ki−1)/2. Aunque existen y se siguen proponiendo y formulando nuevos tipos de redes complejas, pasamos a enumerar algunos tipos de redes clasificadas según sus propiedades [Albert and Barabási,2002]: Redes regulares. Son redes cuyos nodos forman una malla regular. A pesar de ser poco realistas por su topología, permiten describir algunas propiedades básicas de la dinámica de los sistemas que representan. Un ejemplo clásico, que permite solución analítica, es el modelo de Ising. Redes aleatorias. Propuestas por Erdös y Rényi en 1959, esta red se define a partir de Nnodos de forma que los enlaces entre ellos tiene la misma probabilidad pde existir, por lo 22 Capítulo 1. Introducción que el grado de cada nodo es parecida al grado medio de la red y de forma que la distribución de sus grados es de Poisson [Erdös and Rényi,1960]. En la Figura 1.8 podemos ver un ejemplo de este tipo de red con N= 70 nodos y probabilidad de conexión p= 0.2. En esta red, si el número de nodos es suficientemente grande (N≫1), el número de enlaces es aproximadamente pN(N−1)/2. Estas redes, por su regularidad, tampoco parecen ser muy realistas, pero el modelo mostró propiedades interesantísimas desde el punto de vista de la mecánica estadística: en este tipo de redes se dan transiciones de fase, es decir, al variar poco a poco la probabilidad p, algunas propiedades de la red cambian drásticamente cuando pasan cierto valor umbral (un punto crítico). 0 5 10 15 20 Degree 0 2 4 6 8 10 12 Frequency Figura 1.8: Red aleatoria con N= 70 nodos y probabilidad de conexión p= 0.2. A la derecha podemos ver representada la distribución de los grados de la red. Red de mundo pequeño. Estas redes fueron ideadas por Watts y Strogatz en 1998 inspirándose en redes sociales [Watts and Strogatz,1998]. La red se construye a partir de una red unidimensional con condiciones de contorno periódicas (circular) de Nnodos en la cual cada nodo de conecta con los m vecinos más cercanos. A partir de ahí, cada nodo se reconecta con probabilidad pde forma aleatoria, de forma que el grado medio no cambia: ⟨k⟩=m. Sin embargo, incluso para valores 23 David García Selfa pequeños de p, la longitud de camino mínima promedio se acorta considerablemente, de forma que resulta fácil conectar nodos alejados en pocos pasos, de ahí su nombre pues este tipo de red se inspira en redes sociales en las que se exploraba el número de personas intermedias que tienen relaciones comunes entre dos personas cualesquiera, surgiendo a menudo la expresión “ ¡qué pequeño es el mundo! ”. En la Figura 1.9 se muestra la construcción de una red pequeño mundo con N= 40 ym= 4 en la cual hemos ido tomando p= 0,p= 0.3 yp= 1. Red libre de escala. Son un modelo de red más realista que debemos a Albert y a Barabási y que añaden la característica de ir creciendo en el tiempo como ocurre en muchas redes reales [Albert and Barabási,2002]. En esta red se van añadiendo nodos en cada paso de forma que la probabilidad de enlazarse a los nodos ya existentes es proporcional a los grados de dichos nodos en ese instante, es decir, los nuevos nodos tienen mayor probabilidad de conectarse a los nodos centrales (los más conectados). Estas redes presentan una distribución de grados en forma de ley de potencias, esto es, P(k)∼k−γ, con 2< γ < 3, así como propiedades de red de mundo pequeño. Al cumplir los grados una ley de potencias, la red adquiere naturaleza fractal siendo, por tanto, libre de escala. En la Figura 1.10 se puede observar un ejemplo de este tipo de red. En esta tesis trabajaremos con redes aleatorias que nos servirán para simular las interacciones entre osciladores químicos en un tanque continuamente agitado. 1.7. Inestabilidad de Turing Anteriormente, hablamos de la morfogénesis, esto es, del proceso biológico que establece la forma que va a tomar un organismo vivo. Este proceso se da en la biología a muchos niveles: desde la 24 Capítulo 1. Introducción p=0 01234 Degree 0 5 10 15 20 25 30 35 40 Frequency p=0 p=0.3 0 1 2 3 4 5 6 Degree 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Frequency p=0.3 p=1 02468 Degree 0 2 4 6 8 10 Frequency p=1 Figura 1.9: Red de mundo pequeño con N= 40 ym= 4 en la cual hemos ido tomando p= 0,p= 0.3yp= 1. A la derecha podemos ver representadas las distribuciones de los grados de la red que, aunque vas cambiando, su grado medio (⟨k⟩= 4) no cambia. 25 David García Selfa 100101102 Degree 100 101 102 Frequency Figura 1.10: Red libre de escala: los nodos nuevos tienen mayor probabilidad de conectarse a los nodos con mayor grado de los ya existentes. A la derecha podemos ver representada la distribución de los grados de la red (en escala logarítmica), en la cual se puede apreciar que la distribución sigue una ley de escala. estructura de una sola célula hasta la formación colectivos multicelulares ordenados y tejidos y, a partir de estos, órganos y organismos completos. Por supuesto, para comprender los diferentes mecanismos implicados en la vida, se necesita de la biología molecular y la bioquímica, especialmente de la genética. Sin embargo, aunque la genética controle la creación de patrones, los genes no pueden producirlos por sí mismos: sólo proveen una receta, pues no poseen información sobre en qué parte del organismo (o del futuro organismo) se encuentran. Se necesita, entonces, otro mecanismo adicional, y es aquí donde Turing propone que agentes químicos puedan reaccionar y difundirse de una forma determinada, y bajo ciertas condiciones, para establecer patrones espaciales estacionarios heterogéneos de concentraciones de estos agentes químicos (morfogenes) [Turing,1952;Murray,2003]. En la naturaleza existen vistosos y famosos ejemplos, como los patrones únicos en las pieles de ciertos mamíferos (cebras, leopardos, peces, etc.) Es notorio cómo animales con una información genética casi idéntica, como una madre de cebra y su cría, pueden mostrar patrones muy diferentes si durante el proceso de morfogénesis de la cría se da cualquier variación en los parámetros (bifurcación) del mecanismo de reacción-difusión que 26 Capítulo 1. Introducción haga surgir patrones distintos a los habituales (ver Figura 1.11). Figura 1.11: Anomalía en los patrones de la piel de una cebra. Foto cortesía del Dr. J. Carlos Mouriño Gallego, tomada en Kenia en 2019. Aclaremos qué son los patrones de Turing y qué propiedades tienen. Para ello, nos vamos a basar en el artículo del equipo de De Kepper, en el que dieron a conocer la primera prueba experimental de dichos patrones [Castets et al.,1990]. Así, citando literalmente, estas estructuras son patrones estacionarios de concentraciones debidos únicamente al acoplamiento de procesos de reacción y difusión, dejando de lado, por tanto, cualquier estructura de tipo hidrodinámica. Los patrones de Turing son resultado de fenómenos de ruptura de la simetría espontáneos asociados a las bifurcaciones de estados estacionarios (autoorganización) y corresponden a soluciones estacionarias estables de un sistema de ecuaciones de reacción-difusión (ec. 1.3.1). Por otra parte, presentan propiedades como que se caracterizan por una longitud de onda característica que no depende de los parámetros geométricos del sistema (salvo en casos en los que el sistema sea del orden de unas pocas longitudes de onda, pues tiene que ajustarse a las condiciones de contorno). Y, por último, mencionar condiciones necesarias para que se den estos patrones: (a) que la cinética de las reacciones incorpore realimentación positiva de las especies activadoras (como la autoca27 David García Selfa tálisis) y procesos inhibitorios, y (b) la difusión del inhibidor debe ser más rápida que la del activador (hay excepciones, como cuando los patrones se propician por perturbaciones finitas y se deben a bifurcaciones secundarias). En el Capítulo 2 (Métodos numéricos y experimentales) haremos un estudio de estabilidad lineal para determinar cuándo se dan los patrones de Turing. 28 2. Métodos En este capítulo se detalla cómo se han construido los osciladores químicos que se estudian en la tesis, tanto experimentalmente como teórica y numéricamente. La población de osciladores químicos se ha construido mediante cuentas de resina que permiten intercambiar especies químicas con un medio activo. Las superficies de las cuentas llevarán fijado el catalizador de una reacción de Belousov-Zhabotinsky, mientras que el medio activo será una solución de dicha reacción libre de catalizador. Por otra parte, se presentan las bases teóricas de las técnicas numéricas empleadas en el estudio de los diferentes sistemas tratados y sus comportamientos: análisis de estabilidad y factores de crecimiento para la inestabilidad de Turing descrito mediante redes complejas, los análisis de bifurcación y continuidad y, por último, el modelo de aproximación de fases para sistemas sincronizados. 2.1. Osciladores químicos en medio activo En esta tesis se ha trabajado en varios capítulos con colectivos de osciladores químicos formados por cuentas de resina cargadas con catalizador y sumergidas en una reacción de Belousov-Zhabotinsky (BZ) libre se catalizador. En la mayoría de capítulos, el estudio ha sido numérico, pero basado en este sistema que se ha llevado a la experimentación en numerosos estudios, como en nuestro Capítulo 5. En esta sección describimos cómo se construyen estos osciladores en el laboratorio. Para obtener 20 mL de solución BZ libre de catalizador se usó la receta que pasamos a detallar: se mezclaron las concentraciones específicas de ácido sulfúrico, bromuro sódico, bromato sódico, ácido 29 David García Selfa malónico y agua como se muestra en la Tabla 2.1. Solución stock Cantidad Concentración final [H2SO4] = 5 M 2.6 mL [H2SO4] = 0.65 M [NaBr] = 0.5M 3.5 mL [NaBr] = 0.0875 M [NaBrO3] = 2 M 4.1 mL [NaBrO3] = 0.41 M [CH2(COOH)2] = 1.5M 1.8 mL [CH2(COOH)2] = 0.135 M H2O8.0 mL Tabla 2.1: Receta para obtener 20 mL de solución BZ libre de catalizador. Por otra parte, cargamos con ferroína ([Fe(o−phen)3]2+) dos grupos de cuentas de resina esféricas de tamaños distintos: Dowex®50wx4, forma hidrógeno-hidrógeno, 50-100 mesh (diámetros aproximados de entre 300 y 150 µm): cargadas con 1 µmol de ferroína por 1 g de cuentas de resina. Dowex®50wx4, forma hidrógeno-hidrógeno, 100-200 mesh (diámetros aproximados de entre 150 y 75 µm): cargadas con 5µmol de ferroína por 1 g de cuentas de resina. Para cargar las cuentas de resina con ferroína se procedió de la siguiente manera: En un vaso de precipitados se añaden 10 g de cuentas de resina y 50 mL de una concentración: •[[Fe(o−phen)3]2+] = 2 ×10−5M (para obtener 1 µmol de ferroína por 1 g de cuentas de resina). •[[Fe(o−phen)3]2+] = 1 ×10−4M (para obtener 5 µmol de ferroína por 1 g de cuentas de resina). Se mezclan durante 4 horas con un agitador magnético. Se filtran al vacío y se limpian con agua ligeramente acidulada (con 306 mL de [H2SO4] = 0.1M). 30 David García Selfa Se dejan secar. Todos estos preparados se realizan a temperatura ambiente de 23 ±1oC. 2.2. Análisis de estabilidad y factores de crecimiento en patrones de Turing En el Capítulo 4 de esta tesis vamos a tratar con patrones de Turing usando una topología de red compleja. Así, vamos a ver cómo se realiza el análisis de estabilidad y poder determinar si se dan dichos patrones (inestabilidad de Turing). Partimos de un sistema de reacción-difusión de dos variables (el activador uy el inhibidor v):      ∂u ∂t =f(u, v) + Du∇2u ∂v ∂t =g(u, v) + Dv∇2v .(2.2.1) Suponiendo que existe un punto de equilibrio estable homogéneo (u∗, v∗), esto es, (f(u∗, v∗)=0 g(u∗, v∗)=0 (2.2.2) y, dado que es estable, la matriz jacobiana del sistema homogéneo J=fufv gugv cumple que det(J)>0y Tr(J)<0en el punto de equilibrio (u∗, v∗). Si introducimos, ahora, en pequeña perturbación en el punto de equilibrio δx=δu δv=u−u∗ v−v∗, podemos linealizar la ecuación de reacción difusión como dδx dt =Jδx+D∇2δx,(2.2.3) 31 David García Selfa 38 3. Resonancias en un “oregonator” no isotermo forzado Vamos a estudiar el comportamiento de un oscilador químico sometido al acoplamiento de una fluctuación externa. En concreto, acoplaremos una reacción de Belousov-Zhabotinsky (BZ) a una fuente de calor externa oscilante. Hasta ahora la reacción BZ ha venido siendo estudiada bajo condiciones isotermas, pero nosotros vamos a tener que formular un modelo no isotermo para poder ver cómo se comporta frente a fluctuaciones de temperatura debido a fuentes de calor. Consideraremos, entonces, un modelo de oregonator no isotermo partiendo del modelo de tres variables e incorporando la temperatura como cuarta variable (no como parámetro), añadiendo un balance energético al sistema de ecuaciones. El efecto de la temperatura en las velocidades de la reacción se incluye a través de la ley de Arrhenius (constantes de velocidad de reacción dependiente de la temperatura). Para poder modelizar una situación realista en un entorno de laboratorio, el sistema contará con refrigeración externa y estará sometido a las fluctuaciones de un forzamiento externo por radiación infrarroja. Mediante simulaciones numéricas y estudios paramétricos, encontramos que el sistema muestra resonancias debido a oscilaciones inducidas. Hemos descubierto que una fuente externa de calor (por ejemplo, un diodo emisor de luz) en condiciones resonantes puede ser usado para inducir una bifurcación de Hopf en un reactor BZ experimental. 3.1. Introducción Ya hablamos en la introducción de la reacción BZ como ejemplo primordial de dinámica química no lineal fuera del equilibrio que proporcionó los primeros osciladores químicos [Belousov,1959; 39 David García Selfa Zhabotinsky,1964]. Los osciladores naturales están sometidos normalmente a fluctuaciones, tanto estocásticas como periódicas, de las variables físicas externas, tales como la radiación luminosa o la temperatura ambiente, que pueden modificar, modular o incluso acompasar su comportamiento y han sido estudiadas tanto en modelos matemáticos [Simakov and Pérez-Mercader,2013b;Serna et al.,2017] como experimentalmente para una BZ fotosensible [Simakov and Pérez-Mercader,2013a]. La compensación de la temperatura en relojes circadianos es un buen ejemplo en la naturaleza [Ruoff et al.,1999;Bodenstein et al.,2012;François et al.,2012]. En este capítulo, partimos de un oregonator de tres variables y lo modificamos adecuadamente añadiendo una cuarta variable, la temperatura, para estudiar el efecto de la modulación en la misma. La incorporación de la temperatura se hace incluyendo un balance energético que nos proporciona la evolución del calor. En estudios anteriores, tanto teóricos como experimentales, se ha explorado el forzamiento externo de una reacción BZ mediante modulación de la temperatura, de forma que se llega a proponer que la reacción BZ posee un mecanismo de compensación de la temperatura [Ruoff, 1995;Masia et al.,2001;Bánsági et al.,2009;Novak et al.,2011]. Esto tiene gran importancia, pues al ser la reacción BZ un modelo simplificado del ciclo de Krebs, se ha tomado como modelo de oscilador químico para predecir y controlar sistemas tan complejos como el músculo cardiaco [Krinsky,1984]. Sin embargo, en los estudios precedentes, la temperatura siempre ha sido considerada como un parámetro forzado externamente, no como una variable dependiente, por lo que no han investigado de forma explícita la evolución del calor. Sabemos que los efectos del calor, como el producido en las reacciones, pueden afectar a la dinámica de un sistema químico de forma muy significativa, sobre todo en sistemas concentrados, tal como se ha descrito en varios estudios sobre sistemas heterogéneos a altas temperaturas, donde se autoinducen pulsos oscilantes de temperatura, por ejemplo, en la oxidación catalítica de monóxido de carbono [Flytzani-Stephanopoulos et al.,1980;Scheintuch,1981]. Pero el efecto de la evolución dinámica del calor en osciladores ho40 Capítulo 3. Oregonator no isotermo mogéneos a baja temperatura, como la reacción BZ, considerando la temperatura como una variable dependiente (no un parámetro) no ha sido estudiado hasta ahora. Así, para el propósito descrito, consideraremos un modelo de oregonator para una reacción BZ que trabaja en un reactor de flujo y que incluye un balance de energía para reflejar las condiciones no isotermas con términos de calor de reacción, de intercambio de calor a través de un refrigerante y de radiación infrarroja (IR), este último como forzamiento externo. Esta descripción completa del sistema químico nos permite incorporar fluctuaciones periódicas indirectamente en la temperatura de la mezcla de la reacción remedando condiciones realistas. Además, demostraremos que la modulación periódica externa a través de radiación IR puede hacer surgir comportamientos no triviales como oscilaciones pseudodeterministas y resonancias. IR LED coolant IN coolant coolant OUT thermocouple REDOX electrode insulation stirring bar outlet (waste) inlet (feed solutions) Figura 3.1: Esquema del reactor agitado de flujo continuo y refrigerado considerado en el capítulo. 41 David García Selfa 3.2. Modelo de “oregonator” no isotermo Consideremos un reactor agitado de flujo continuo (CSTR por sus siglas en inglés: continuously stirred tank reactor), bien mezclado, rodeado de un refrigerante (refrigeración externa), en condiciones estacionarias y a escala de laboratorio (como el usado experimentalmente en [Simakov and Pérez-Mercader,2013b]) y cuyo esquema se muestra en la Figura 3.1. Haciendo un balance de masas para la reacción BZ y tomando el esquema de oregonator presentado en la introducción (1.4.2) más un término de flujo (recordemos que estamos alimentando el reactor)[Zhong and Xin,2000], llegamos a              dX dt =k1AY −k2XY +k3AX −2k4X2−kfX dY dt =−k1AY −k2XY +fk5BZ −kfY dZ dt = 2k3AX −k5BZ −kfZ ,(3.2.1) donde tenemos las tres variables X= [HBrO2],Y= [Br−]yZ= [Ce4+], así como los parámetros A= [BrO− 3],B= [CH2(COOH)2], fyP= [HOBr], las constantes cinéticas correspondientes a las reacciones ki, con i= 1, ..., 5, y la velocidad de flujo alimentación kf. La dependencia de las velocidades de reacción kicon la temperatura viene dada por la ley de Arrhenius ki(T) = ki0exp −Ei R1 T−1 T0 ,(3.2.2) siendo ki0=ki(T0)las constantes de velocidad de reacción a la temperatura de referencia T0,Rla constante de los gases ideales y Eilas energías de activación de las reacciones. La evolución del calor en el CSTR sometido a radiación externa (la fuente de IR que, como propusimos, puede implementarse mediante LED) y con refrigeración externa se modela mediante el 42 Capítulo 3. Oregonator no isotermo balance de energía VRρwCpw dT dt =−VRX i ∆HRiri+QfX j CfjHfj−Hj+˙ Q+VRRIR , (3.2.3) donde VRes el volumen del reactor, ρwes la densidad del agua, Cpw el la capacidad calorífica a presión constante del agua, ∆HRies la variación de entalpía de la reacción i: ∆HRi= ∆H0 Ri+ZT Tr ∆CpidT ≈∆H0 Ri, pues ∆Cpi≈0, ries la velocidad de la reacción i,Qfes el caudal volumétrico del flujo de alimentación, Cfjes la concentración de la especie jde la alimentación, Hfjes la entalpía específica de la especie j,˙ Qes la potencia calorífica transferida, Ses el área de la superficie de transmisión de calor, Tces la temperatura del refrigerante, Ues el coeficiente de transmisión térmica de la superficie y RIR es la potencia térmica (infrarroja) irradiada por unidad de volumen. Por tanto, tenemos un término que representa la generación de calor mediante las reacciones químicas −VRX i ∆HRiri, un término que representa el intercambio de calor a través del refrigerante que rodea al reactor y que viene dado por la ley de enfriamiento de Newton ˙ Q=US (Tc−T) y un término que representa el calentamiento por radiación infrarroja VRRIR. 43 David García Selfa Las velocidades de reacción vienen dadas por r1(T) = k1AY H2 r2(T) = k2XY H r3(T) = k2AX r4(T) = k4X2 r5(T) = k5BZ, con ki(T)cumpliendo la ley de Arrhenius como se ha dicho más arriba. Asumimos, ahora, que el término QfPjCfjHfj−Hjes despreciable y sustituimos las velocidades de las reacciones y la ley de enfriamiento de Newton en la ecuación 3.2.3: VRρwCpw dT dt =−VR(k1AY H2∆HR1+k2XY H∆HR2+k3AXH∆HR3 +k4X2∆HR4+k5BZ∆HR5) +US(Tc−T) + VRRIR . (3.2.4) Usando las misma concentraciones y tiempo adimensionales que en 1.4.3 (escalado de Tyson) y la temperatura adimensional θ=T Tc , llegamos a la ecuación adimensional para la evolución de la temperatura γdθ dt =ψ(x, y, z, θ) + (1 −θ) + ω, (3.2.5) 44 Capítulo 3. Oregonator no isotermo siendo γ(θ) = VRρwCpwk5B US ω=VRRIR USTc ψ(x, y, z, θ) = −VR USTcα1y+α2xy +α3x+α4x2+α5z,con α1(θ) = k1k3A2H2 k2 ∆HR1 α2(θ) = k2 3A2H2 2k4 ∆HR2 α3(θ) = k2 3A2H2 2k4 ∆HR3 α4(θ) = k2 3A2H2 4k4 ∆HR4 α5(θ) = k2 3A2H2 k4 ∆HR5. Combinando las ecuaciones del oregonator 3.2.1 adimensionalizada y la ecuación 3.2.5, obtenemos, por fin, las ecuaciones de nuestro modelo de oregonator no isotermo                      ϵdx dτ =x(1 −x) + y(q−x)−ϵκx δdy dτ =fz −y(q+x)−δκy dz dτ =x−z−κz γdθ dt =ψ(x, y, z, θ) + (1 −θ) + ω ,(3.2.6) donde κ=kf k5B. Es importante resaltar que, como las constantes de velocidad de reacción (ki(T)) dependen de la temperatura, todos los parámetros adimensionales de la ecuación 3.2.6, esto es, κ,ϵ,δ,qyγtienen 45 David García Selfa una dependencia no lineal de la temperatura a través de la ley de Arrhenius. Podemos ver inmediatamente, entonces, que hay dos mecanismos distintos para producir fluctuaciones en la temperatura del sistema: (a) uno de carácter aditivo a partir de la intensidad de la radiación IR (ω) y (b) otro de carácter multiplicativo a través de la temperatura del refrigerante (1−θ). En este trabajo sólo vamos a considerar el primero, pero tenemos como trabajo futuro considerar el segundo. 3.3. Resultados y discusión Para llevar a cabo las simulaciones, integramos numéricamente las ecuaciones 3.2.6 usando un método de Runge-Kutta de cuarto orden, con los valores iniciales x(0) = 1,y(0) = 0.1,z(0) = 0 y θ(0) = 1, los valores de los parámetros de la Tabla (3.1), usado en el trabajo [Pullela et al.,2009], además de: Constante de los gases y temperatura de referencia: R= 8.314 J mol K,T0= 298 K. Geometría del reactor y propiedades del fluido: VR= 0.05 R, S= 1 dm2,ρw= 1kg L,Cpw= 4.186 kJ kgK. Parámetros de la reacción: A0= 0.24 M, B0= 0.18 M, f= 1, kf= 5.556 ×10−4s−1. Usaremos como parámetros de control la concentración de H2SO4 de la solución (H0), la temperatura del refrigerante (Tc) y el coeficiente de transmisión de calor (U). En primer lugar, analizamos el comportamiento del sistema sin forzamiento externo (ω= 0 en 3.2.6), con un coeficiente de transmisión de calor U= 10 J m−2K−1s−1. En el rango de valores empleados para H0(0.07 −0.16 M), tiene lugar una bifurcación de Hopf para Tc= 283 −313 K (10 −40 oC), como podemos ver en la Figura 3.2. Los valores de estos parámetros se han elegido de forma que se acerquen a los valores usados en experimentos reales. 46 Capítulo 3. Oregonator no isotermo Reacción 1 2 3 4 5 ∆HRi[kJ mol−1] 43 -71 98 -114 0 Ei[kJ mol−1] 60 25 60 75 70 ki02[M−3s−1]106[M−2s−1]10 [M−2s−1]2×103[M−1s−1]1[M−1s−1] Tabla 3.1: Entalpías, energías de activación y constantes de velocidad a la temperatura de referencia (obtenidas de [Pullela et al., 2009]). 47 David García Selfa de la frecuencia normalizada, obteniéndose las curvas de resonancia características. Las ramas de estrella rojas de los paneles (b) y(d) representan las oscilaciones deterministas a baja frecuencia obtenidas cuando se cruza la bifurcación de Hopf (tal como sucede en la Figura 3.7b). Las filas de la Figura 3.8 corresponden a distintos valores del coeficiente de transmisión superficial del calor (U). Disminuyendo el valor de U, el pico de resonancia se hace menos pronunciado. También se puede observar cómo la temperatura correspondiente al pico de resonancia se hace mayor conforme disminuye U, al ser peor la transmisión del calor entre el refrigerante y la reacción. En la Figura 3.9, en la que comparamos las curvas de resonancia para distintos valores de U, se puede observar este fenómeno con más claridad. En el recuadro de la Figura 3.9 representamos la frecuencia normalizada de resonancia frente al coeficiente de transmisión superficial del calor. Vamos a analizar los tres casos presentados en la Figura 3.7. En el primer caso (fila superior: (a) y (b)), tenemos una frecuencia por debajo de la frecuencia natural, el comportamiento es trivial y el sistema muestra oscilaciones cuando se cruza la bifurcación de Hopf. En el último caso (fila inferior: (e) y (f)), tenemos otro caso trivial y se alcanza un estado estacionario (no oscilante) tras un estado transitorio, pues la dinámica del sistema no es lo suficientemente rápida y sólo es capaz de “percibir” la intensidad promedio de la radiación infrarroja que hace que la temperatura del reactor posicione el sistema más allá de la bifurcación de Hopf. El caso no trivial (fila intermedia: (c) y (d)), es el más interesante, pues para frecuencias de forzamiento cercanas a la frecuencia natural del sistema (que corresponde a la que tiene en las proximidades de la bifurcación de Hopf) el sistema presenta oscilaciones cuando éste debiera tender a un estado estacionario por estar mucho más allá de la bifurcación de Hopf. Si analizamos, ahora, la Figura 3.8 vemos cómo la temperatura máxima decrece mientras que la mínima crece conforme se aumenta la frecuencia de forzamiento, hasta que ambas convergen para frecuencias altas. Esto pone de manifiesto la capacidad calorífica 54 Capítulo 3. Oregonator no isotermo 10-2 10-1 100 Normalized forcing frequency, f f/f0 290 300 310 320 T (K) 10-2 100 Normalized forcing frequency, f f/f0 10-4 10-2 100 xmax-xmin 10-2 10-1 100 Normalized forcing frequency, f f/f0 290 300 310 320 T (K) 10-2 100 Normalized forcing frequency, f f/f0 10-4 10-2 100 xmax-xmin 10-2 10-1 100 Normalized forcing frequency, f f/f0 290 300 310 320 T (K) 10-2 100 Normalized forcing frequency, f f/f0 10-4 10-2 100 xmax-xmin (a) U=10 J m-2K-1s-1 (b) U=10 J m-2K-1s-1 (c) U=7 J m-2K-1s-1 (d) U=7 J m-2K-1s-1 (e) U=5 J m-2K-1s-1 (f) U=5 J m-2K-1s-1 Figura 3.8: Curvas de resonancia obtenidas para distintos valores de l coeficiente de transmisión superficial del calor (U). En la columna izquierda se muestra la temperatura del reactor frente a la frecuencia normalizada (la línea punteada corresponde a la temperatura a la cual tiene lugar la bifurcación de Hopf: T= 296.3 K). En la columna derecha se muestra la amplitud de la concentración adimensional del activador frente a la frecuencia normalizada. Los valores empleados de los parámetros son: H0= 0.1M, Tc= 283 K, U= 10 J m−2K−1,RIR0= 0.028 kW/L yAf= 0.018 kW/L. La frecuencia natural es f0= 0.00459 Hz[García-Selfa et al., 2019]. 55 David García Selfa 10 0 Normalized forcing frequency, ff/f0 10 -4 10 -3 10 -2 10 -1 xmax -x min U=10.3 J m -2K-1 s-1 U=8 J m -2K-1 s-1 U=7 J m -2K-1 s-1 U=5 J m -2K-1 s-1 0.7 0.8 0.9 1 (ff/f0)peak 4 6 8 10 U (J m-2K-1s-1) Figura 3.9: Curvas de resonancia para distintos valores del coeficiente de transmisión superficial del calor (U). El recuadro muestra la frecuencia de resonancia normalizada frente al coeficiente de transmisión superficial del calor: la curva punteada muestra el ajuste que cumple la ley de potencias U= 5.817 (ff/f0)5.817 peak + 4.463. Los valores empleados de los parámetros son: H0= 0.1M, Tc= 283 K, U= 10 J m−2K−1,RIR0= 0.028 kW/L y Af= 0.018 kW/L. La frecuencia natural es f0= 0.00459 Hz [García-Selfa et al., 2019]. 56 Capítulo 3. Oregonator no isotermo finita de la reacción (ecuación 3.2.3) que afecta a la velocidad de la respuesta dinámica del sistema frente a fluctuaciones externas. En los trabajos previos a éste, en los que la temperatura era un parámetro y no una variable, este comportamiento no se podía observar y, aparentemente, la temperatura se podía aumentar incluso a altas frecuencia, lo que es físicamente imposible. Por tanto, el hecho de considerar la temperatura como una variable dependiente (esto es, al incorporar un balance energético) nos permite una representación más precisa de la dinámica del sistema, especialmente a altas frecuencias, cuando el efecto de la capacidad calorífica cobra importancia. En la columna derecha podemos ver la amplitud de la concentración del activador (como diferencia entre los valores máximo y mínimo) frente a la frecuencia de forzamiento normalizada. Para frecuencias relativamente bajas, la amplitud de la concentración oscila entre dos valores, tal como hace la temperatura. Este comportamiento es trivial: las oscilaciones deterministas (ramas superiores, representadas por estrellas rojas, como se ven en la Figura 3.8 (b) y (d)) tienen lugar cuando la temperatura está por debajo de la bifurcación de Hopf (como en la Figura 3.7 (a) y (b)). Por encima de ciertas frecuencias, las oscilaciones no tienen lugar, pues la frecuencia de forzamiento es mucho más rápida que la dinámica del sistema. La variación de amplitud debida al forzamiento externo (círculos azules) es no nula porque los cambios periódicos en la temperatura del reactor obliga a las variables del sistema a ajustarse a las condiciones (que cambian periódicamente). Pero el comportamiento no trivial se observa en las proximidades de la frecuencia natural (ff/f0≈1), esto es, la frecuencia de resonancia, donde se observan importantes oscilaciones en la amplitud de la concentración del activador (Figura 3.8 (b)) a pesar de estar más allá de la bifurcación de Hopf (Figura 3.8 (a)). Estas oscilaciones, a pesar de tener amplitudes relativamente bajas, son significativamente mayores que las oscilaciones pasivas que se observan a otras frecuencias, como se puede ver el pico que aparece para ff/f0≈1. Sin embargo, podemos ver que para la frecuencia a la que tiene lugar el pico (resonancia), la temperatura no oscila y el sistema se encuentra 57 David García Selfa siempre en el dominio no-oscilatorio, más allá de la bifurcación de Hopf. Por tanto, son oscilaciones de naturaleza no determinista y es claramente un fenómeno resonante que propicia la aparición de oscilaciones macroscópicas. Téngase en cuenta que el oregonator es un modelo de la reacción BZ, que es un oscilador químico no lineal. Por otra parte, señalar que nuestro oscilador está siempre amortiguado debido al la refrigeración (ecuación 3.2.3). Así, la resonancia observada no es idéntica a la observada en en los osciladores lineales clásicos, en los cuales la resonancia es catastrófica y aparece para ff/f0= 1 en ausencia de amortiguación. En la Figura 3.9 se muestra el comportamiento resonante para distintos valores de U. Cada máximo de las curvas de amplitud (xmáx. −xmín.) corresponde a una resonancia. Como ocurre con los osciladores lineales forzados amortiguados clásicos, las resonancias tienen lugar a menores frecuencias normalizadas conforme aumenta la amortiguación (recuadro interior de la Figura 3.9). Así, para valores de Uinferiores, la refrigeración es menos eficaz y la temperatura del reactor aumenta, desplazando al sistema hacia la región de equilibrio, lejos de la bifurcación de Hopf. Este comportamiento coincide con el observado en los osciladores lineales forzados amortiguados clásicos. 3.3.2. Leyes de potencias Por otra parte, es interesante subrayar la ley de potencias que sigue el ajuste de Ufrente a las frecuencias de resonancia normalizadas: U= 5.817 ff f05.817 peak + 4.463 .(3.3.2) Para los valores del coeficiente de transmisión del calor representados hemos obtenido, además de las frecuencias de resonancia normalizadas ((f/f0)peak), los valores medios de las temperaturas entre los cuales el reactor oscila simétricamente y, también, la diferencia entre esta temperatura media y la temperatura del refrigerante (∆T= (Tmáx. −Tmín.)/2−Tc), con Tc= 283 K. Los resultados los podemos ver en la Tabla 3.2. 58 Capítulo 3. Oregonator no isotermo U[Jm−2K−1s−1] (f/f0)peak T= (Tmáx. −Tmín.)/2[K] ∆T[K] 5.0 0.6667 311.15 28.15 7.0 0.8627 303.00 20.00 8.0 0.9216 300.50 17.50 10.3 1.0000 296.60 16.60 Tabla 3.2: Resultados obtenidos de las frecuencias de resonancia normalizadas, temperaturas medias del reactor y diferencias entre las temperaturas medias y la del refrigerante para distintos valores de U. 0.04 0.06 0.08 0.1 U S (J K-1s-1 ) 0.4 0.6 0.8 1 (f/f0)peak y=-0.007541x-1.414+1.188 r2=0.9999 Figura 3.10: Representación de los valores de las frecuencias de resonancia normalizadas frente a la tasa de transferencia de calor (circulitos rojos) y de la curva de ajuste (línea discontinua). Para los cálculos se ha tomado un área S= 1 dm2. Obsérvese el magnífico ajuste a una ley de potencias. En la Figura 3.10 podemos ver cómo las frecuencias de resonancia (producidas por una fuente de calor oscilante externa) sigue una ley de potencias con el intercambio de calor ambiental. En la Figura 3.11 podemos observar cómo las diferencias de temperatura entre el reactor y el ambiente (refrigerante) también sigue una ley de potencias con la tasa de calor transferido al ambiente. Para 59 David García Selfa 0.04 0.06 0.08 0.1 U S (J K-1s-1) 10 20 30 40 " T (K) y=-1.173x-1.051+0.8134 r2=1 Figura 3.11: Representación de las diferencias entre las temperaturas medias y la del refrigerante frente a la tasa de transferencia de calor (circulitos azules) y de la curva de ajuste (línea discontinua). Para los cálculos se ha tomado un área S= 1 dm2. Obsérvese el magnífico ajuste a una ley de potencias. los cálculos de ha tomado una superficie de área S= 1 dm2. Por otra parte, podemos ver que las diferencias de temperatura y las frecuencias de resonancia también siguen una ley de potencias: ∆T=−29.52 f f01.668 peak + 43.15 , r = 0.9998 (3.3.3) Este resultado se puede relacionar con algunos comportamientos biológicos como, por ejemplo, el ritmo metabólico de animales endotermos (de sangre caliente). Los endotermos son capaces de equilibrar su transferencia de calor con el ambiente con su ritmo metabólico, de forma que su temperatura se mantiene constante cuando varía la temperatura ambiente [Fristoe et al.,2015]. Este balance energético se modela mediante la ley de Scholander-Irving: Tb−Ta=B C,(3.3.4) 60 Capítulo 3. Oregonator no isotermo donde Tbes la temperatura corporal del endotermo, Taes la temperatura del ambiente, Bes el ritmo de calor metabólico producido yCes el coeficiente de transferencia de calor. En nuestro sistema, podemos considerar T−Tc=B+RIRVR U S ,(3.3.5) donde, como hemos venido usando hasta ahora, Tes la temperatura del reactor, Tces la temperatura del refrigerante, Ues el coeficiente de transferencia de calor, Sel área de la superficie de transferencia de calor, RIRVRrepresenta la tasa de calor transmitida por radiación infrarroja. Y Bpuede considerarse como el ritmo de calor metabólico del reactor (calor que se produce mediante las reacciones químicas). Nótese que el calor metabólico depende de la temperatura a través de la ley de Arrhenius [Gillooly et al.,2001]. Así, podemos tomar B= ∆T U S −RIRVR.(3.3.6) Con los valores que hemos venido usando hasta ahora (RIR0VR= 0.0014 Js−1oscilando con una amplitud de 0.0009 Js−1) y los valores mostrados en la Tabla 3.2, obtenemos B≈1.4J s−1. Así, como tenemos una temperatura ambiente constante (la temperatura del refrigerante) y una fuente de calor externa oscilante, nuestro sistema mantiene el ritmo de calor metabólico cambiando el coeficiente de transferencia de calor. Como nuestro sistema modela (de una forma muy simplificada) cierto tipo de comportamiento biológico, pues el origen de la reacción BZ está en el ciclo de Krebs, también sería interesante estudiar la dependencia del tamaño del reactor con el ritmo de calor metabólico con el objetivo de hallar algún tipo de ley alométrica [West et al.,2002;West and Brown,2005]. 3.4. Conclusiones En este capítulo hemos estudiado un modelo de oregonator de tres variables al que le hemos añadido la temperatura como cuarta va61 David García Selfa riable a partir de la incorporación de un balance energético (calor). El balance energético incluía los calores generados por la reacción, un término de refrigeración y un término de forzamiento externo a partir de radiación infrarroja. El acoplamiento entre el calor y el balance másico se ha hecho a través de las constantes de velocidad de las reacciones mediante la ley de Arrhenius. El modelo resultante de oregonator no isotermo de cuatro variables no sólo tiene en cuenta explícitamente las variaciones de temperatura de la reacción, sino que permite investigar la respuesta dinámica del sistema frente a fluctuaciones externas en flujos de calor. Tomando los parámetros cinéticos experimentales hallados en la literatura, así como un conjunto de parámetros adecuados y ajustados a la descripción del reactor a escala real de laboratorio, hemos podido predecir la bifurcación de Hopf a temperaturas realistas (10 −40 oC). Hemos sometido al sistema a fluctuaciones aditivas para controlar su dinámica (radiación infrarroja modulada periódicamente). Este parámetro de control se ha elegido por ser de fácil implementación en el laboratorio, ya sea mediante LED o mediante láser. Como trabajo futuro, se pretende estudiar el efecto de fluctuaciones multiplicativas (a través del refrigerante), así como usar fluctuaciones estocásticas, tanto aditivas como multiplicativas. Hemos demostrado, a través de simulaciones numéricas, que el sistema presenta resonancia y oscilaciones mantenidas incluso cuando, para los parámetros seleccionados, la configuración más probable es un estado de equilibrio no oscilatorio. Las resonancias se han observado para frecuencias de forzamiento cercanas a la frecuencia natural (correspondiente a oscilaciones cercanas a la bifurcación de Hopf). Debido al efecto de la amortiguación en la ecuación del calor, no se han observado resonancias de orden superior. Podríamos decir que la periodicidad en flujos de calor ambientales pueden inducir un comportamiento periódico en un oscilador químico. Este hallazgo puede abrir el camino para controlar osciladores químicos y abre la puerta a extender este mecanismo de control a osciladores bioquímicos y, por tanto, a relojes biológicos. Asimismo, gracias a algunas leyes de potencias que se mani62 Capítulo 3. Oregonator no isotermo fiestan en los resultados obtenidos y que relacionan el intercambio de calor con el ambiente y las frecuencias de resonancia de las oscilaciones, podríamos estar ante la base teórica de algunas leyes alométricas relacionadas con el ritmo metabólico celular. 63 David García Selfa (a)(b) (c) (d) 0 500 1000 node i 0 500 1000 ki 0 500 1000 node i 0 500 1000 ki Figura 4.1: Distribución de la conectividad en la red. Representación de la red con (a) p= 0.01 y(b) p= 0.5(sólo se representan el 10 % de los nodos para facilitar la observación). Se muestra la conectividad de cada nodo para (c) p= 0.01 y(d) p= 0.5[Mussa Juane et al., 2020]. 70 Capítulo 4. Inestabilidad de Turing en osciladores acoplados competitiva es mayor. En nuestro caso, los coeficientes de difusión cruzada serán negativos. Completamos el modelo con las ecuaciones para el catalizador: dzi dt =h(xi, yi, zi)para i≥1 zo= 0 para i= 0 .(4.2.6) 4.3. Análisis de estabilidad lineal Veamos cómo llevamos a cabo el análisis de estabilidad para el sistema de (n+ 1) ecuaciones 4.2.5 y4.2.6. Tomamos una perturbación del sistema considerado de forma que (xi, yi, zi) = (x∗ i, y∗ i, z∗ i) + (δxi, δyi, δzi), siendo (x∗ i, y∗ i, z∗ i)el estado estacionario. Expresamos esta perturbación como                δdxi dt =fxδxi+fyδyi+fzδzi+dxiKex n X j=0 Lij((1 + a12y∗ i)δxj+a12x∗ iδyj) δdyi dt =gxδxi+gyδyi+gzδzi+dyiKex n X j=0 Lijδyj δdzi dt =hxδxi+hyδhi+hzδzi . (4.3.1) Definiendo x, de dimensión (3n+ 2), como x= (x0, y0, x1, y1, z1, ..., xn, yn, zn)⊺, expresamos la perturbación como δdx dt = (δx0, δy0, δx1, δy1, δz1, ..., δxn, δyn, δzn)⊺. Así, el sistema 4.3.1 queda reducido a δdx dt = (J+L)δx,(4.3.2) 71 David García Selfa donde Jes una matriz (3n+ 2) ×(3n+ 2) que contiene la (n+ 1) matrices jacobianas de los nodos J=     J00 J1... 0Jn      (4.3.3) yLcontiene los elementos del término de la red L=      d0KexL00 (1 + a12 y∗ 0)d0L00a12 x∗ 0··· d0KexL0n(1 + a12 y∗ 0)d0L0na12 0 0d0KexL00 ··· 0d0Kex L0n0 d1KexL10 (1 + a12 y∗ 1)d1L10a12 x∗ 1··· d1KexL1n(1 + a12 y∗ 1)d1L1na12 0 0d1KexL10 ··· 0d1Kex L1n0 0 0 ··· 0 0 0 . . . . . . . . . . . . . . . . . . . . . dnKexLn0(1 + a12 y∗ n)dnLn0a12x∗ n··· dnKexLnn(1 + a12 y∗ n)dnLnna12 0 0dnKexLn0··· 0dnKex Lnn 0 0 0 ··· 0 0 0       (4.3.4) y la perturbación se puede desarrollar a partir de un conjunto de (3n+ 2) ondas planas δx=AαeλαtΦα,(4.3.5) siendo λαlos (3n+ 2) autovalores de la matriz (J+L)y siendo Φαlos (3n+ 2) autovectores asociados. El conjunto λαse conoce como factores de crecimiento y miden la rapidez con la que evoluciona un nodo tras la perturbación. Para que se pueda dar una inestabilidad de Turing se tiene que dar que, al menos, algún λα>0[Nakao and Mikhailov,2010]. 4.4. Resultados y discusión Una vez definido el modelo del sistema, llevamos a cabo simulaciones numéricas para demostrar que se alcanzan inestabilidades de Turing (recordemos que se han llegado a encontrar inestabilidades de Turing en redes con otros modelos, como el de Murray-Mimura [Nakao and Mikhailov,2010]), pero no con osciladores químicos con reacción BZ. Señalar, también, que dado que la red no posee métricas espaciales, los patrones de Turing que encontramos no vienen 72 Capítulo 4. Inestabilidad de Turing en osciladores acoplados caracterizados por una longitud de onda características: son patrones de Turing deslocalizados y, así, aparecen dos estados bien diferenciados de forma que los nodos evolucionan hacia un estado o hacia el otro [McCullen and Wagenknecht,2016;Wolfrum, 2012]. Por otra parte, con el fin de validar el modelo, veremos cómo se manifiestan patrones conocidos del sistema para valores de los parámetros adecuados del modelo tales como sincronización, quimeras, etc. En las simulaciones realizadas usamos un método de integración de Runge-Kutta de orden 4 y el número de osciladores químicos ha sido n= 1000 (además de la solución externa) y, además, se ha realizado el análisis de estabilidad descrito en el apartado anterior. En la Figura 4.2 se muestran los resultados de la simulación para valores de los parámetros que cumplen las condiciones de inestabilidad de Turing según el criterio descrito en el análisis de estabilidad de la sección anterior (algunos valores de los factores de crecimiento positivos). Inicialmente (paneles a, b), los estados son aproximadamente iguales, salvo una pequeña perturbación alcanzada mediante fluctuaciones estocásticas (ruido uniforme). Después de un tiempo (paneles c, d), las perturbaciones crecen hasta alcanzar un estado estacionario final con dos valores de las concentraciones claramente diferenciados (paneles e, f), de forma que unos nodos permanecen con valores cercanos a los iniciales y otros evolucionan a valores mucho mayores. En la Figura 4.3 se pueden ver los histogramas con el número de nodos en cada estado. Indicar que las variables yy zvarían de forma similar a la x. En la Figura 4.4 podemos ver la evolución de x(t)para cada nodo. En la Figura 4.5 y en la Figura 4.6 se muestra cómo surgen diferentes estados según los conjuntos de parámetros seleccionados para las ecuaciones 4.2.5 y4.2.6. Los parámetros p(probabilidad de colisión entre cuentas de resina que nos remite a la cantidad de osciladores interactuando) y Kex (constante de intercambio: intensidad de la interacción de los osciladores). Para los valores de parámetros seleccionados en la Figura 4.5 (h= 7,ϵ= 10,δ= 1000), aparecen muerte de oscilaciones (zona roja del diagrama de fases mostrado 73 David García Selfa (a) (b) (c) (d) (e) (f) Figura 4.2: Evolución de la concentración del activador (x(t)) bajo condiciones de inestabilidad de Turing. En la columna de la izquierda se muestran los grafos de la red para tres tiempos: (a) t= 0 u. t., (c) t= 80 u. t.y (e) t= 1000 u. t. (u. t. = unidades de tiempo adimensionales). El color de capa nodo se corresponde con la concentración según el mapa de la derecha. La columna de la derecha ((b),(d) y (f)) representa el valor de xpara cada nodo, correspondiendo a los mismos tiempos que los representados en la columna derecha (con el mismo código de color). A t= 0u. t. el sistema con n+1 nodos está en estado estacionario (a, b) con una pequeña perturbación de hasta el 10 % (mediante un ruido distribuido uniformemente). Este estado perturbado evoluciona en el tiempo a dos estados claramente diferenciados (dos últimas filas de la figura). Una vez se alcanza la situación (e, f), el sistema permanece estable. Valores de los parámetros: q= 0.000001,ϵ= 10,δ= 1000, h= 7,Kex = 1,⟨Vn⟩= 0.0001,Vs= 0.1,a12 =−30,a21 = 0,p= 0.01, n= 1000 [Mussa Juane et al., 2020]. 74 Capítulo 4. Inestabilidad de Turing en osciladores acoplados (a) t= 0 u. t. (b) t= 40 u. t. (c) t= 80 u. t. (d) t= 1000 u. t. Figura 4.3: Histograma con los valores de xde cada nodo bajo condiciones de inestabilidad de Turing para: (a) t= 0 u. t., (b) t= 40 u. t., (c) t= 80 u. t. y (d) t= 1000 u. t. (estado estacionario final). Valores de los parámetros: q= 0.000001, ϵ= 10,δ= 1000,h= 7,Kex = 1,⟨Vn⟩= 0.0001,Vs= 0.1,a12 =−30,a21 = 0, n= 1000 yp= 0.01 [Mussa Juane et al., 2020]. 75 David García Selfa Figura 4.4: Evolución temporal de x(t)para las 1000 cuentas de resina (nodos) con el mismo mapa de color empleado en Figura 4.2. (Con los mismos valores de los parámetros que en la Figura 4.2) [Mussa Juane et al., 2020]. en el panel (a) y evolucion en el panel (b)) e inestabilidad de Turing (zona amarilla, con valores de Kex bajos y de paltos en el panel (a) y evolución en el panel (c)). Para los valores de la Figura 4.6 (h= 0.7,ϵ= 0.001,δ= 0.1) observamos cuatro regiones en el diagrama de fases (panel (a)): (I), en rojo, corresponde a un estado estacionario homogéneo, asociado a la evolución del panel (b);(II), en verde, corresponde a un estado sincronizado con todas las cuentas de resina oscilando con la misma amplitud y el medio, en fase, pero con distinta amplitud, como vemos en el panel (c); III, en azul, tenemos un estado quimera (subconjuntos de osciladores, cada uno con un estado diferente: sincronización, caos, muerte de oscilaciones, etc. [Tinsley et al.,2012]), como se puede ver en el panel (d);IV, en magenta, tenemos estados similares a los de (I), pero con otros valores de los estados estacionarios, como vemos en el panel (e). Además de los diferentes valores del parámetro estequiométrico h, las diferencias fundamentales de las dinámicas descritas en las Figuras 4.5 y4.6 la establecen los parámetros ϵyδ, que fijan las escalas temporales: para valores menores de estos pa76 Capítulo 4. Inestabilidad de Turing en osciladores acoplados rámetros, la dinámica es más rápida y predominan las oscilaciones, impidiendo que se formen inestabilidades de Turing, cuyos estados finales corresponden a estado de equilibrio. La Figura 4.7 representamos los factores de crecimiento, calculados según lo expuesto en la sección anterior. Estos factores miden la rapidez con la que un nodo evoluciona a un nuevo estado estacionario después de haber sido perturbado en su estado de equilibrio inicial. Los valores positivos de los factores de crecimiento están asociados a inestabilidades de Turing: el estado estacionario inicial se hace inestable y el sistema evoluciona hasta otro estado estable distinto. 4.5. Conclusiones Se ha modelizado matemáticamente una población de osciladores químicos en un medio activo (solución externa) y la colisión entre cuentas mediante una red compleja en la que cada cuenta es un nodo (cuya conexión con otro nodo se corresponde a colisiones entre ellos) y el medio activo es otro nodo conectado con todos los demás. Hemos introducido, también, una difusión generalizada que incluye difusión cruzada, siendo ésta crucial para que se puedan producir inestabilidades de Turing. A partir de estas herramientas hemos demostrado la existencia de patrones de Turing no localizados que pueden explicar ciertos fenómenos observados en la naturaleza [Miller,2010;Drescher et al., 2009;Mackie et al.,1988;Haddock et al.,2005;Alvarez-Maubecin et al.,2000]. Para poder poner de manifiesto tanto las inestabilidades de Turing como otros comportamientos descritos con anterioridad (sincronización, quimeras, etc.) se han realizado simulaciones numéricas intensivas y se ha analizado el efecto de los parámetros del modelo. Hemos hallado que el número de colisiones entre las cuentas de resina, determinado por p, ha de ser suficientemente alto. También hemos visto que la dinámica del sistema, determinada por ϵyδha de ser suficientemente lenta para suprimir oscilaciones. Se han calculado los factores de crecimiento con el fin de determinar 77 David García Selfa (a) 10-4 10-3 10-2 10-1 100 p 10 -2 10-1 100 101 102 K ex (I) (II) (b) (I) (c) (II) t (t.u.) 01000 t (t.u.) 0 1000 xi(t) xi(t) Figura 4.5: (a) Distintos estados en el espacio de parámetros Kex vs. p(log-log). Cada estado observado viene descrito por la evolución temporal de x(t)mostrados en (b) y (c). (b) Muerte de oscilaciones: todos los nodos dejan de oscilar y acaban en un mismo estado estacionario final y la solución externa (nodo 0) muestra el mismo comportamiento. (c) Patrones de Turing no localizados: en el sistema aparecen dos subconjuntos de cuentas con estados finales diferentes. Las líneas continuas representan la evolución de x(t)de las cuentas y las líneas discontinuas la de la solución externa. Valores de los parámetros: q= 0.000001,ϵ= 10,δ= 1000, h= 7,⟨Vn⟩= 0.0001,Vs= 0.1,a12 =−30,a21 = 0,n= 1000 [Mussa Juane et al., 2020]. 78 Capítulo 4. Inestabilidad de Turing en osciladores acoplados (b) (c) (d) (e) (a) 10 10 10 10 10 -4 -3 -2 -1 0 p -2 -1 0 1 10 10 10 10 10 2 K ex (I) (II) (III) (IV) (I) (II) (III) (IV) 0 3000 0 3000 0 3000 0 3000 t (t.u.) t (t.u.) t (t.u.) t (t.u.) xi(t)xi(t) xi(t)xi(t) Figura 4.6: (a) Distintos estados en el espacio de parámetros Kex vs. p(log-log). Cada comportamiento viene descrito por la evolución de x(t)en (b), (c), (d) y (e). (b) Muerte de oscilaciones. (c) Osciladores sincronizados (todas las cuentas sincronizadas, pero con distinta amplitud para la solución externa). (d) Estado quimera (algunos nodos sincronizados, otros caóticos, otros estacionarios, etc. (e) Muerte de oscilaciones con diferente estado final para la solución externa. Valores de los parámetros: q= 0.000001,ϵ= 0.001,δ= 0.1,h= 0.7,⟨Vn⟩= 0.0001, Vs= 0.1,a12 =−30,a21 = 0,n= 1000 [Mussa Juane et al., 2020]. 79 David García Selfa 5.2.2. Simulaciones numéricas Con el fin de realizar simulaciones numéricas que permitan respaldar los experimentos realizados, se usa el modelo de reaccióndifusión dado por el siguiente sistema de ecuaciones diferenciales                ∂x(r, t) ∂t =1 ϵ(r)(x(1 −x) + y(q(r)−x)) + Dx(r)∇2x ∂y(r, t) ∂t =1 δ(r)(2h(r)z−y(q(r) + x)) + Dy(r)∇2y ∂z(r, t) ∂t =C(r)(x−z) + Dz(r)∇2z , (5.2.4) donde las variables x,yyzrepresentan las concentraciones adimensionalizadas de activador, inhibidor y catalizador, respectivamente, y los parámetros dependen de la posición r; de hecho, hay tres dominios en los que los parámetros tienen valores distintos: dentro de la cuenta mayor (tipo A), dentro de la cuenta menor (tipo B) y en el resto (solución externa). C(r)vale 1en las cuentas (donde hay fijado catalizador) y 0en el medio externo (libre de catalizador). Se resuelven las ecuaciones usando un método de Runge-Kutta de cuarto orden con los valores de los parámetros: A: radio de 6 a 10 píxeles, ϵ= 0.01,δ= 0.001,q= 0.0018, h= 0.27. B: radio de 1 a 5 píxeles, ϵ= 0.01,δ= 0.001,q= 0.004, h= 0.35. Medio externo: ϵ= 0.01,δ= 0.015,q= 0.003,h= 0.64. 5.3. Resultados y discusión Analizamos lo obtenido en las dos series de experimentos: caracterización de los osciladores químicos y sincronización de los osciladores químicos. 86 Capítulo 5. Sincronización de dos osciladores químicos diferentes 5.3.1. Caracterización de los osciladores químicos En esta serie de experimentos el único, e imprescindible, objetivo fue el de caracterizar los osciladores, esto es, medir sus tamaños y periodos. Tras una larga serie de experimentos, se consigue tener una cantidad suficiente de medidas para poder caracterizar de forma sifnificativa los osciladores. En la Figura 5.2 se muestra la distribución de periodos de los osciladores de tipo Ay de tipo B(en el panel (a) y(b), respectivamente). Los histogramas en las partes superior y derecha de cada panel representan, respectivamente, las distribuciones de diámetros y de periodos de los osciladores. En la parte central de las gráficas de cada panel se observa la clara correlación entre diámetros y periodos. Como se sabe por la literatura, el periodo es menor para diámetros mayores de las cuentas de resina [Yoshikawa et al.,1998]. Además, la concentración de ferroína ligada a la cuenta de resina no parece reflejar dependencia con el periodo de las oscilaciones, ya que el periodo de las oscilaciones depende de la velocidad a la que se difunden las especies químicas activas a través de la cuenta la resina hacia el medio activo y esto depende de la razón entre el área de la superficie de la cuenta y su volumen, pero la ferroína está ligada a la superficie de la cuenta de resina, por lo que no presenta difusión. En cualquier caso, demostramos con esta figura que el procedimiento empleado para cargar las cuentas de resina ha permitido obtener dos grupos claramente diferenciados por los tamaños y con un rango de periodos también claramente diferenciado. 5.3.2. Sincronización de los osciladores químicos Una vez caracterizados, estudiamos el acoplamiento y la sincronización de los osciladores químicos en la segunda serie de experimentos midiendo las distancias, tamaños, amplitudes, fases y diferencias de fase de cada par de cuentas de resina sumergidas en el medio activo (siempre una cuenta de cada tipo). En primer lugar, como se muestra en la Figura 5.3, estudiamos experimentalmente los periodos de los osciladores de tipo A y de ti87 David García Selfa (a) (b) 80 90 100 110 120 130 140 150 diameter ( m) 40 50 60 70 80 90 period (s) 140 160 180 200 220 240 260 280 300 diameter ( m) 36 38 40 42 44 46 48 50 52 period (s) Figura 5.2: Distribución de periodos de los osciladores químicos. (a) Cuentas de entre 150 y 300 µm de diámetro, aproximadamente, cargadas con 1 µmol de ferroína por gramo de cuentas de resina. (b) Cuentas de entre 75 y 150 µm de diámetro, aproximadamente, cargadas con 5 µmol de ferroína por gramo de cuentas de resina. Podemos ver que el periodo de las oscilaciones depende del tamaño de las cuentas de forma aproximadamente lineal: el periodo es menor para cuentas más grandes. 88 Capítulo 5. Sincronización de dos osciladores químicos diferentes (a) (b) (1) (3)(2) (1) (3)(2) Figura 5.3: Datos experimentales. (a) Periodos y (b) razón entre los periodos (TA TB) de las oscilaciones en función de la distancia de separación de los centros de las cuentas de resina. En (a), los círculos azules se corresponden con los osciladores tipo A(las más grandes) y las estrellas naranja con los de tipo B(las más pequeñas). (a) (1) (3)(2) (b) (1) (3)(2) Figura 5.4: Datos se las simulaciones numéricas. (a) Periodos y (b) razón entre los periodos (TA TB) de las oscilaciones en función de la distancia de separación de los centros de las cuentas de resina. En (a), los círculos azules se corresponden con los osciladores tipo A(las más grandes) y las estrellas naranja con los de tipo B(las más pequeñas). 89 David García Selfa Sync 1:2 No Sync Sync 1:1 t=97.6 s t=91.6 s t=103.6 s t=115.6 st=109.6 s t=100.6 s t=91.6 s t=109.6 s t=127.6 st=118.6 s t=10 s t=8 s t=12 s t=16 st=14 s (a) (b) (c) Figura 5.5: Distintos casos de sincronización obtenidos experimentalmente: (a) Sincronización 1:1, con los osciladores en contacto. (b) Sincronización 1:2, con los osciladores a distancias intermedias. (c) No sincronizados, con los osciladores alejados. 90 Capítulo 5. Sincronización de dos osciladores químicos diferentes (a) (b) (c) Sync 3:2 No Sync Sync 1:1 Figura 5.6: Distintos casos de sincronización obtenidos numéricamente: (a) Sincronización 1:1, con los osciladores en contacto. (b) Sincronización 3:2, con los osciladores a distancias intermedias. (c) No sincronizados, con los osciladores alejados. 91 David García Selfa po B (Figura 5.3 (a)) y la razón entre los mismos, esto es, TA TB(Figura 5.3 (b)) en función de la distancia entre los centros de las cuentas. Se distinguen tres zonas, según el acoplamiento entre osciladores: en la zona (1) el acoplamiento es más fuerte y los osciladores están completamente sincronizados, teniendo ambos el mismo periodo; en la zona (2) el acoplamiento es más débil y aparecen formas más complejas de sincronización, desde osciladores que sincronizan con un periodo que es el doble que el de su pareja (sincronización 1:2), u osciladores con el mismo o con distinto periodo, hasta parejas que tratan de sincronizarse pero no acaban de sincronizarse por completo; en la zona (3) los osciladores están tan separados que no hay acoplamiento entre ellos. En segundo lugar, como se muestra en la Figura 5.4, estudiamos los periodos de los osciladores de tipo A y de tipo B (Figura 5.4 (a)) y la razón entre los mismos, esto es, TA TB(Figura 5.4 (b)) en función de la distancia entre los centros de las cuentas, pero usando las simulaciones numéricas basadas en el modelo dado por la ec. 5.2.4. Se distinguen igualmente tres zonas, según el acoplamiento entre osciladores: en la zona (1) el acoplamiento es más fuerte y los osciladores están completamente sincronizados, teniendo ambos el mismo periodo (sincronización 1:1); en la zona (2) el acoplamiento es más débil y aparecen formas más complejas de sincronización (sincronización 1:1, sincronización 2:3, sincronización 3:4, etc.); en la zona (3) los osciladores están tan separados que no hay acoplamiento entre ellos. En las Figuras 5.5 y5.6 se muestran unos fotogramas, experimentales y simulados numéricamenbte, respectivamente, en los cuales se pueden observar distintos tipos de sincronización asociados a los casos correspondientes a las zonas (1),(2) y(3). En las Figuras 5.7,5.8 y5.9 se muestran análisis de una selección, de entre todos los experimentos, de oscilaciones correspondientes a los casos (1),(2) y(3), respectivamente. En estas figuras podemos ver la evolución temporal de la concentración del catalizador de los osciladores químicos a partir del cambio de color periódico debido a la oxidación-reducción de la ferroína (paneles (a) y(c)). Las imá92 Capítulo 5. Sincronización de dos osciladores químicos diferentes time (s) 0400 368.9 mm 112.5 mm (a) (b) (c) (d) Figura 5.7: Caso (1). Los osciladores están muy próximos entre ellos o en contacto. El acoplamiento es fuerte y los osciladores oscilan sincronizadamente. (a) Series temporales de las oscilaciones en la concentración del catalizador: las líneas y puntos azules corresponden a osciladores de tipo Ay los naranja a los de tipo B.(b) Evolución temporal del parámetro de orden de Kuramoto. (c) Evolución temporal de una sección de los fotogramas que atraviesa diametralmente las cuentas de resina. (d) Evolución temporal de la diferencia de fase entre los dos osciladores. 93 David García Selfa time (s) (a) (b) (c) (d) time (s) 0500 273.3 mm 136.7 mm 66.5 mm Figura 5.8: Caso (2). Los osciladores están separados distancias intermedias entre ellos. El acoplamiento es débil y, en este ejemplo, las oscilaciones tienen aproximadamente un periodo el doble la una de la otra. (a) Series temporales de las oscilaciones en la concentración del catalizador: las líneas y puntos azules corresponden a osciladores de tipo Ay los naranja a los de tipo B.(b) Evolución temporal del parámetro de orden de Kuramoto. (c) Evolución temporal de una sección de los fotogramas que atraviesa diametralmente las cuentas de resina. (d) Evolución temporal de la diferencia de fase entre los dos osciladores. 94 Capítulo 5. Sincronización de dos osciladores químicos diferentes (a) (b) (c) (d) time (s) 0300 272.1 mm 122.1 mm 244.3 mm Figura 5.9: Caso (3). Los osciladores están muy alejados entre ellos. No hay acoplamiento y los osciladores no presentan sincronización. (a) Series temporales de las oscilaciones en la concentración del catalizador: las líneas y puntos azules corresponden a osciladores de tipo Ay los naranja a los de tipo B.(b) Evolución temporal del parámetro de orden de Kuramoto. (c) Evolución temporal de una sección de los fotogramas que atraviesa diametralmente las cuentas de resina. (d) Evolución temporal de la diferencia de fase entre los dos osciladores.. 95 David García Selfa y el medio activo. Y los parámetros ϵ,δ,qyhestán asociados a las concentraciones iniciales, las velocidades de las reacciones y la estequiometría. Para el medio activo, la dinámica viene dada por las concentraciones de activador e inhibidor                ϵ∂xs(t) ∂t =xs(t) (1 −xs(t)) + ys(t) (q−xs(t)) +⟨V⟩n Vs Kex nbeads X i=1 (xi(t)−xs(t)) δ∂ys(t) ∂t =−ys(t) (q+xs(t)) + ⟨V⟩n Vs Kex nbeads X i=1 (yi(t)−ys(t)) ,(6.2.2) donde ⟨V⟩nes el volumen medio de las cuentas y Vses el volumen total de la solución externa. Como la solución es libre de catalizador, no se incluye su concentración en estas ecuaciones directamente, pero el medio sí contiene a todas las cuentas cargadas con catalizador por lo que puede presentar oscilaciones. A partir de ahora, nos vamos a centrar en la transición entre el estado sincronizado (todas cuentas y el medio oscilando con la misma frecuencia, pero con distinta amplitud) y el supersincronizado (todas las cuentas y el medio oscilando en fase con la misma frecuencia y con la misma amplitud). Dado que la transición de interés implica que todas las cuentas están sincronizadas, consideramos todos los osciladores, excepto el medio activo, son idénticos y, por tanto, podemos considerar la variedad de sincronización donde todos estos osciladores son idénticos y xi=xb,yi=yb,zi=zbpara todo i∈1, ..., n, denotando el índice ba las cuentas. Así, reducimos el sistema descrito por las (n+ 1) ecuaciones 6.2.1 y6.2.2 a un 102 Capítulo 6. Osciladores químicos acoplados a través de un medio modelo de 5 dimensiones descrito por el sistema de 5 ecuaciones                                  ϵ∂xb(t) ∂t =xb(t) (1 −xb(t)) + yb(t) (q−xb(t)) −Kex (xb(t)−xs(t)) δ∂yb(t) ∂t = 2hzb(t)−yb(t) (q+xb(t)) −Kex (yb(t)−ys(t)) ∂zb(t) ∂t =xb(t)−zb(t) ϵ∂xs(t) ∂t =xs(t) (1 −xs(t)) + ys(t) (q−xs(t)) +ρKex (xb(t)−xs(t)) δ∂ys(t) ∂t =−ys(t) (q+xs(t)) + ρKex (yb(t)−ys(t)) ,(6.2.3) siendo ρ=n⟨V⟩n Vsla densidad del sistema. Ahora, usando el software de bifurcación y continuidad Matcont [Dhooge et al.,2008] y AUTO [Ermentrout,2002], vamos a obtener el diagrama del espacio de estados para nuestro sistema de 5 dimensiones. En la Figura 6.1 representamos el diagrama del espacio GH CPC LPC HH+ LPC HH+ A B (a) (b) Figura 6.1: Diagrama del espacio de estados. (a) Distintos comportamientos observados en el sistema cuando varían Kex yρ. (b) el mismo diagrama del espacio de estados pero tridimensional, al representar la concentración xben el eje vertical [García-Selfa et al., 2021]. de estados correspondiente al sistema de ecuaciones 6.2.3. En el diagrama del espacio de estados mostramos las bifurcaciones mediante curvas continuas. En el panel (a) se puede ver el espacio de 103 David García Selfa V s K ex 0 0.05 0.10 0.15 0.20 0.25 0.005 0.010 0.020 0.030 Figura 6.2: Diagrama del espacio de estados obtenidos experimental y numéricamente en [Ghoshal et al., 2016]. Adaptado con permiso de los autores. estados al variar la constante de intercambio Kex y la densidad ρy podemos comprobar, en primer lugar, que se reproduce el comportamiento obtenido en [Ghoshal et al.,2016] tanto numérica como experimentalmente (Figura 6.2). Asimismo, en el panel (b) vemos cómo se vuelven a reproducir los resultados de [Ghoshal et al.,2016] en la representación tridimensional del mismo diagrama obtenida al incorporar la variable xben el eje vertical. En estos diagramas podemos ver las distintas bifurcaciones implicadas en las transiciones analizadas. En ambas representaciones podemos ver la bifurcación de Hopf generalizada (GH) que separa las dos ramas de bifurcación de Hopf supercrítica (H−, donde el primer coeficiente de Lyapunov es negativo), la bifurcación de Hopf subcrítica (H+, donde el primer coeficiente de Lyapunov es positivo) y la bifurcación de punto de silla de órbitas periódicas (LCP, donde el sistema tiene un ciclo límite no hiperbólico único con multiplicador de Floquet +1). Por otra parte, coincidiendo con el punto GH, tenemos un punto de cúspide de ciclos (CPC) en el cual se bifurcan el comportamiento supercrítico del subcrítico. El único comportamiento no descrito 104 Capítulo 6. Osciladores químicos acoplados a través de un medio 1 0.5 xb H 0 0.5 0.2 zb 0.4 0.75 LPC LPC LPC LPC LPC LPC 10 A B (a) (b) Figura 6.3: (a) Ciclos límites correspondientes a la transición entre el punto de equilibrio A (muerte de oscilaciones) al punto B de estado de supersincronización pasando por el estado de sincronización entre la bifurcación de Hopf supercrítica, H, y la bifurcación de punto se silla de órbitas periódicas ,LPC. La línea discontinua A-B corresponde a lo mostrado en la Figura 6.1a.(b) Frecuencias correspondientes a estas transiciones [García-Selfa et al., 2021]. en esta figura, respecto a lo descrito en la Figura 2-d de [Ghoshal et al.,2016], es la ausencia de sincronización cuando Kex es muy pequeña, pues en nuestro modelo reducido todos los osciladores son idénticos. En la Figura 6.3 (a) podemos observar las órbitas periódicas en la transición desde el equilibrio (equivalente a la muerte de oscilaciones del sistema de 3n+ 2 dimensiones) del punto Aal estado de sincronización del punto B, pasando a través de una bifurcación de Hopf supercrítica Hy la bifurcación de puntos de silla de órbitas periódicas LPC. Obsérvese que el cambio rápido de ciclo límite es una explosión Canard, que se da incluso el modelo de oregonator simple [Bo Peng et al.,1991;Brøns and Bar-Eli,1991;Krupa and Szmolyan,2001]. En el panel (b) podemos ver el descenso drástico de frecuencias en la explosión Canard, por lo que esta explosión puede observarse también en el colectivo de oregonators sincronizados a través del acoplamiento con un medio activo, además de en un simple oregonator. 105 David García Selfa 6.3. Modelo de aproximación de fases Para obtener ls periodos de las oscilaciones en los estados de sincronización y supersincronización, vamos a usar el modelo de aproximación de fases descrito en la sección 2.4. Para un oscilador dado, su estado puede describirse mediante su posición a lo largo de su ciclo límite, esto es, su fase. Si tenemos dos osciladores desacoplados, sus fases se pueden describir sobre la superficie de un toro. Incluso si dos osciladores con ciclos límites estables están débilmente acoplados, el toro persiste y se puede describir su estado mediante sus fases [Nakao,2016;Pietras and Daffertshofer,2019]. Si el acoplamiento es fuerte, no podemos asegurar esta descripción. Resolvemos numéricamente (con XPPAUTO Ermentrout [2002]) nuestro problema a partir del sistema formado por un oscilador que representa a todas las cuentas de resina (pues todas son idénticas y están sincronizadas, con densidad ρ) desacoplado del resto del sistema (el medio activo). Ahora, para estudiar el acoplamiento entre cuentas y medio activo, tomamos dos copias idénticas del sistema desacoplado                          ∂xbi(t) ∂t =1 ϵ(xbi(t) (1 −xbi(t)) + ybi(t) (q−xbi(t))) ∂ybi(t) ∂t =1 δ(2hzbi(t)−ybi(t) (q+xbi(t))) ∂zbi(t) ∂t =xbi(t)−zbi(t) ∂xsi(t) ∂t =1 ϵ(xsi(t) (1 −xsi(t)) + ysi(t) (q−xsi(t)) + ρKex (xbi(t)−xsi(t))) ∂ysi(t) ∂t =1 δ(−ysi(t) (q+xsi(t)) + ρKex (ybi(t)−ysi(t))) , (6.3.1) donde i= 1,2. Para el conjunto de parámetros correspondiente a un régimen oscilatorio descrito más arriba, las cuentas y el medio activo de cada copia del sistema oscila con la misma frecuencia natural ω0. Ahora, acoplamos las cuentas de una copia del sistema con el medio activo de la otra copia del sistema. Concretamente, acopla106 Capítulo 6. Osciladores químicos acoplados a través de un medio mos la cuenta del sistema i= 1 con el medio activo del sistema i= 2 a partir de G1=      −1 ϵ(xb1−xs2) −1 δ(yb1−ys2) 0 0 0       ,(6.3.2) con una intensidad de acoplamiento dada por la constante K= Kex 2(1+ρ). Esto no permite obtener la descripción mediante las fases si esta constante es pequeña. Por simplicidad, introducimos la constante intensidad de acoplamiento dentro de la función Hide forma que la fase del oscilador 1evoluciona como dθ1(t) dt =ω0+H1(θ2−θ1).(6.3.3) Ahora, calculamos los periodos de las oscilaciones tanto en el estado sincronizado como en el estado supersincronizado y comparamos el modelo completo no lineal y la aproximación de fases. Para la dinámica de fases sabemos que las cuentas y el medio activo oscilan en fase (diferencia de fase ϕ0= 0 = 2π) con frecuencia Ωy periodo T=2π Ω, luego dθ1(t) dt =Ω=2π T=ω0+H1(0) .(6.3.4) Los parámetros del modelo correspondiente al sistema completo son: q= 0.002 , ϵ = 0.01 , δ = 0.015 , h = 0.70 .(6.3.5) En el análisis que realizamos, se han tomado como parámetros de control Kex (la constante de tasa de intercambio de especies químicas) y ρ(la densidad de cuentas en el medio). En la Figura 6.4 se muestran las curvas de respuesta de fase (PRC) Z(θ)para una cuenta con densidad ρ= 1.2. 107 David García Selfa 0/2 3 /2 2 0 3 6 9 Zxb ( ) 0/2 3 /2 2 -4 -2 0 2 Zyb () (a) (b) Figura 6.4: Curvas de respuesta de fases (PRC) para una cuenta con ρ= 1.2. Conjunto de parámetros químicos usados: q= 0.002 ,ϵ= 0.01 ,δ= 0.015, h= 0.70 [García-Selfa et al., 2021]. 108 Capítulo 6. Osciladores químicos acoplados a través de un medio 0/2 3 /2 2 -4 -2 0 2 4 H1( ) Kex=0.077 0/2 3 /2 2 -3.9 -3.8 -3.7 -3.6 -3.5 H1( ) Kex=0.033 (a) (b) Figura 6.5: Dos casos de H1(ϕ)con densidad ρ= 1.2correspondientes a (a) estado sincronizado con Kex = 0.033 y (b) estado supersincronizado con Kex = 0.077. Obsérvese la notable diferencia en las escalas del eje yde ambas figuras [García-Selfa et al., 2021]. 109 David García Selfa Las funciones de interacción H1(ϕ)correspondientes a los estados de sincronización y supersincronización con ρ= 1.2se muestran en la Figura 6.5. Al ir variando los valores de los parámetros, vamos obteniendo los periodos de las oscilaciones colectivas en buena aproximación. En la Figura 6.6 podemos ver los periodos calculados para cada valor de Kex, con ρ= 1.2. Estos periodos han sido calculados usando al integración numérica del sistema 6.2.3 y usando el modelo de aproximación de fases. Podemos ver la transición entre sincronización y supersincronización como una discontinuidad debido a la bifurcación descrita entre ciclos límite de distinta naturaleza. 0 0.03 0.06 0.09 0.12 Kex 0 5 10 15 T (t.u.) =1.2 Phase approximation Numerical integration Figura 6.6: Periodos de los osciladores en función de Kex con densidad ρ= 1.2. La discontinuidad en el periodo muestra la transición entre los estados de sincronización y supersincronización [García-Selfa et al., 2021]. Los cambios cualitativos entre los estados de sincronización y supersincronización también se pueden evidenciar mediante los cambios en la función de interacción de fases [Hesse et al.,2017]. De hecho, los cambios en esta función son un indicador de la bifurcación subyacente. Para mostrar este efecto en el sistema de osciladores químicos, desarrollamos la función de interacción en series de Fourier y observamos cómo cambian los modos de Fourier al variar los valores de los parámetros. 110 Capítulo 6. Osciladores químicos acoplados a través de un medio Así, desarrollamos la función de interacción en series de Fourier de senos y cosenos: H1(ϕ) = a0 2+ ∞ X k=1 [akcos (kϕ) + bksin (kϕ)] .(6.3.6) En la Figura 6.7 se muestran los coeficientes del desarrollo de Fou0.01 0.04 0.07 0.10 Kex -80 -40 0 40 ak (1/2)a0 a1 a2 a3 a4 a5 0.01 0.04 0.07 0.10 Kex -30 -20 -10 0 10 bk b0=0 b1 b2 b3 b4 b5 Figura 6.7: Coeficientes de la serie en senos y cosenos de Fourier. Debido a la diferente naturaleza de los ciclos límite en los estados de sincronización y supersincronización, tenemos diferentes modo en la función de interacción para esos estados [García-Selfa et al., 2021]. rier obtenidos para distintos valores de Kex.H1(ϕ)se ha obtenido con el método mencionado anteriormente para cada valor de Kex. En el estado sincronizado podemos ver que H1(ϕ)∼a0 2+a1cos (ϕ) + b1sin (ϕ),(6.3.7) 111