Ecuaciones de fase para osciladores no lineales con ciclo límite estable

FACULTAD DE MATEMÁTICAS DEPARTAMENTO DE ECUACIONES DIFERENCIALES Y ANÁLISIS NUMÉRICO Trabajo Fin de Grado Ecuaciones de fase para osciladores no lineales con ciclo límite estable. Alejandra Sánchez Cañas ———————— Dirigido por: Soledad Fernández García 2 Índice general Resumen I Abstract III Introducción V 1. Modelos Matemáticos Básicos en Neurociencia 1 1.1. ConceptosPrevios ............................ 1 1.2. Modelo matemático para una Neurona: Integrate-and-Fire . . . . . . 2 1.2.1. Modelo de Neurona Integrate-and-Fire (IF) . . . . . . . . . . 3 1.2.2. Modelo de Neurona Leaky IF (LIF) . . . . . . . . . . . . . . . 4 2. Herramientas básicas 5 2.1. Funciones lipschitzianas . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2. Existencia y unicidad de solución . . . . . . . . . . . . . . . . . . . . 6 2.3. Órbitas................................... 8 2.4. Estabilidad ................................ 9 3. Bifurcaciones 15 3.1. Bifurcación Silla Nodo sobre una Órbita Periódica (SNIPER) . . . . . 16 3.2. Bifurcaciones de Hopf Supercrítica y Subcrítica . . . . . . . . . . . . 17 3.3. Bifurcación Homoclina . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4. Curvas de respuesta de fase (PRC) 23 4.1. Isócronas.................................. 23 4.1.1. Definición ............................. 23 4.1.2. Ejemplos.............................. 24 4.1.3. Aplicaciones............................ 26 4.2. Curvas de Respuesta de Fase (PRC) . . . . . . . . . . . . . . . . . . . 26 4.2.1. Definición y metodología . . . . . . . . . . . . . . . . . . . . . 26 4.3. Ecuacionesdefase ............................ 28 4.3.1. Reducciones de fase . . . . . . . . . . . . . . . . . . . . . . . . 28 4.3.2. CálculodelaPRC......................... 29 4.4. PRCs para modelos específicos . . . . . . . . . . . . . . . . . . . . . . 30 4.4.1. Curvas de respuesta de fase en un entorno de bifurcaciones de codimensión uno con impulsos periódicos . . . . . . . . . . . . 30 4.4.2. Modelos neuronales de dimensión uno . . . . . . . . . . . . . . 42 4.4.3. Precisión de las PRCs analíticas . . . . . . . . . . . . . . . . . 45 3 4ÍNDICE GENERAL Resumen La compresión del funcionamiento del sistema nervioso es uno de los grandes retos de la comunidad científica actual. En este trabajo, trataremos problemas clásicos en neurociencia desde el punto de vista de la matemática aplicada. En particular, utilizaremos las ecuaciones de fase para osciladores no lineales con ciclo límite estable para el análisis de modelos de tipo Integrate-and-Fire. Así, estudiaremos cómo los intervalos de reacción de las poblaciones de osciladores neuronales responden a los estímulos sensoriales, y conectaremos esto con los mecanismos de computación y modulación neuronal. Las simulaciones numéricas que aparecen en el trabajo, se han realizado con el programa MATLAB, y las gráficas cualitativas, con el programa de representación vectorial InkScape. i ii RESUMEN Abstract Understanding how the nervous system works, is one of the greatest challenges of nowadays scientific community. In this work, we will address classical problems in neuroscience from the point of view of applied mathematics. In particular, we will use the phase equations for non-linear oscillators with a stable limit cycle for the analysis of Integrate-and-Fire models. Thus, we will study how the reaction intervals of neural oscillator populations respond to sensory stimuli and we will connect this to the mechanisms of neural computation and modulation. The numerical simulations in this work have been carried out using MATLAB program, and the qualitative graphs, with the vector representation program InkScape. iii iv ABSTRACT Introducción La biomatemática (o matemática biológica) es una rama de la ciencia encargada de modelar los procesos biológicos mediante técnicas propias de las matemáticas. La neurociencia es el estudio del sistema nervioso, el cual está compuesto por el cerebro, la médula espinal y las redes de células nerviosas sensitivas o motoras, llamadas neuronas, en todo el cuerpo. El objetivo de la neurociencia es comprender cómo funciona el sistema nervioso para producir y regular emociones, pensamientos, conductas y funciones corporales básicas, incluidas la respiración y mantener el latido del corazón. La compresión del funcionamiento del sistema nervioso es uno de los grandes retos de la comunidad científica actual. En este trabajo, trataremos la neurociencia desde el punto de vista de la matemática aplicada. En particular, utilizaremos las ecuaciones de fase para osciladores no lineales con ciclo límite estable para el análisis de modelos de tipo Integrate-and-Fire. Así, estudiaremos cómo los intervalos de reacción de las poblaciones de osciladores neuronales responden a los estímulos sensoriales, y conectaremos esto con los mecanismos de computación y modulación neuronal. Para la mayoría de los modelos neuronales, observamos que la respuesta a un estímulo fijo es inversamente proporcional a la frecuencia de impulsos antes del estímulo. Las simulaciones numéricas que aparecen en el trabajo, se han realizado con el programa MATLAB, y las gráficas cualitativas, con el programa de representación vectorial InkScape. Este trabajo se organiza como explicamos a continuación. Para comenzar el trabajo, debemos introducir algunos conceptos básicos sobre neurociencia, como el concepto de neurona, su constitución y el proceso de la sinapsis neuronal. También presentaremos la conducta de un modelo general básico de una única neurona, el modelo de Integrate-and-Fire. Esto lo haremos en el Capítulo 1. A continuación, en el Capítulo 2, enunciaremos los resultados teóricos necesarios para llevar a cabo el estudio de las bifurcaciones y las curvas de respuesta de fase. Comenzaremos recordando el concepto y los resultados más destacados sobre ecuaciones diferenciales, como son las funciones lipschitzianas, que serán necesarias para el estudio de existencia y unicidad del Problema de Cauchy, el cual se desarrollará en la segunda sección de este capítulo. Además, veremos algunos resultados sobre órbitas y su estabilidad, como, por ejemplo la aplicación de Poincaré. Después, en el Capítulo 3, vamos a presentar algunas de las distintas bifurcaciones de codimensión 1 más relevantes que aparecen en modelos neuronales. Para ello, comenzaremos introduciendo el concepto de bifurcación y, posteriormente, estudiaremos la dinámica algunos sistemas dependientes de parámetros. Nos centraremos en las birfucaciones Silla-Nodo, de Hopf supercrítica y subcrítica, y homoclina. v 6CAPÍTULO 2. HERRAMIENTAS BÁSICAS espacio vertorial sobre R. En particular, todas las funciones constantes pertenecen aLip(x, ). Definición 2.2. Supongamos que es abierto. Se dice que f:≠æ Rnes una función localmente lipschitziana respecto de la variable xen , si para cada punto (t0,x 0)œexiste una bola abierta de centro dicho punto, contenida en , y tal que la restricción de fa dicha bola es globalmente lipschitziana respecto de la variable x. Es decir, fes localmente lipschitziana respecto de la variable xen , si para cada punto (t0,x 0)œexisten ‘(t0,x 0), y una constante L(t0,x 0)>0, tales que: B((t0,x 0),‘(t0,x 0)) µ, y para todo (t, x1),(t, x2)œB((t0,x 0),‘(t0,x 0)), |f(t, x1)≠f(t, x2)|ÆL(t0,x 0)|x1≠x2|. Si es abierto, se denota Liploc(x, )al conjunto de todas las funciones f: ≠æ Rnque son localmente lipschitzianas respecto de la variable xen . Se comprueba de manera inmediata que Liploc(x, ), con las operaciones usuales de suma de funciones y producto de un número real por una función, es un espacio vectorial sobre R,yLip(x, )µLiploc(x, )de forma estricta. En general, no es sencillo comprobar si una función es o no Lipschitziana. Por ello, es útil la siguiente observación: Observación 2.1. Si f:æRnes tal que existen todas las derivadas parciales de todas las componentes de frespecto de xj,j=1, ..., n y son continuas en , entonces fœLiploc(x, ). Todos los problemas que consideramos a lo largo de este trabajo, en los problemas de Cauchy asociados hay existencia y unicidad local de solución porque las funciones son CŒ, luego se cumple la Observación 2.1 y el Teorema de Picard, el cual enunciaremos en la siguiente sección. 2.2. Existencia y unicidad de solución Sea ™Rn+1. Consideremos el Problema de Cauchy, (2.1). El estudio de la existencia y unicidad de solución local, resulta del Teorema de Picard: Teorema 2.1 (de Picard).Sean ™Rn+1 un conjunto abierto no vacío, y f: ≠æ Rntal que: fœC(;Rn)flLiploc(x, ). Con estas condiciones, para cada (x0,t 0)œ, existe un ”>0, tal que si denotamos I”=[t0≠”, t0+”], existe una y solo una solución del (PC)en I”. 2.2. EXISTENCIA Y UNICIDAD DE SOLUCIÓN 7 Bajo las mismas condiciones, se tienen resultados de existencia y unidad de solución maximal. Definición 2.3. Sean ™Rn+1 un conjunto abierto no vacío y f:≠æ Rntal que: fœC(;Rn)flLiploc(x, ). Para cada (t0,x 0)œdenotemos por (I(t0,x 0),Ï(·;t0,x 0)) a la solución maximal del (PC). Se define el conjunto ={(t, t0,x 0)œRn+2;(t0,x 0)œ,tœI(t0,x 0)} y la función Ï:(t, t0,x 0)œ‘æ Ï(t, t0,x 0)œRn. A la función Ïasí definida se la denomina solución (maximal) del (PC)expresada en función de los datos iniciales. Para establecer el resultado de existencia y unicidad de solución global, se hace uso del siguiente resultado: Teorema 2.2 (Lema de Gronwall).Se tienen los siguientes resultados: Supongamos dados ≠Œ <x 0<x 1<+Œ, dos funciones u,kœC([x0,x 1]),y una constante hœR,tales que k(x)Ø0,y 1. u(x)Æh+⁄x x0 k(s)u(s)ds, ’xœ[x0,x 1]. En este caso, se satisface u(x)Æhesx x0k(s)ds,’xœ[x0,x 1]. 2. u(x)Æh+⁄x0 xk(s)u(s)ds, ’xœ[x1,x 0]. En este caso, se satisface u(x)Æhesx0 xk(s)ds,’xœ[x1,x 0]. Con una aplicación sencilla del Lema de Gronwall, obtenemos el siguiente resultado de unicidad global para el Problema de Cauchy. Teorema 2.3 (de unicidad Global).Sean ™Rn+1 un conjunto abierto no vacío yf:≠æ Rntal que: fœC(;Rn)flLiploc(x, ), y consideremos fijado un punto (t0,x 0)œ. En estas condiciones, (I1,Ï 1)e(I2,Ï 2) son dos soluciones del problema de Cauchy, entonces, Ï1(t)=Ï2(t),’tœI1flI2. 8CAPÍTULO 2. HERRAMIENTAS BÁSICAS Finalmente, presentamos el resultado principal de existencia y unicidad de solución global: Teorema 2.4 (de existencia y unicidad de solución global).Sean ™Rn+1 un conjunto abierto no vacío y f:≠æ Rntal que: fœC(;Rn)flLiploc(x, ). En estas condiciones, para cada (t0,x 0)œdado, existe una única solución maximal o global del problema de Cauchy (PC)que denotamos (I(t0,x 0),Ï(·;t0,x 0)). Además, el intervalo I(t0,x 0)de definición de solución global es abierto. Para introducir la teoría de las Curvas de Respueta de fase (PRCs, por sus siglas en inglés), objeto principal de este trabajo, nos basta con estudiar sistemas autónomos, es decir, la función fno depende explícitamente del tiempo, por lo que nos restringiremos a ese caso. Además, necesitamos conocer el concepto de órbita, que presentamos a continuación. 2.3. Órbitas Sean ™Rn+1 un conjunto abierto no vacío y f:≠æ Rntal que: fœC0(;Rn)flLiploc(x, ). Consideremos el sistema autónomo: ˙x=f(x)=Q c c c c a f1(x) f2(x) . . . fn(x) R d d d d b .(2.2) Analizaremos las propiedades cualitativas de las soluciones del sistema autónomo (2.2). Comenzaremos definiendo el concepto de órbita que pasa por un punto x0œ: Definición 2.4. Sea x0œy denotemos I(x0)=I(0,x 0). Se llama órbita del sistema autónomo (2.2) asociada (o, que pasa por) x0al conjunto “(x0)dado por: “(x0)={Ï(t;0,x 0):tœI(x0)}µ. De la definición se deduce que x0œ“(x0), y que, en general, “(x0)es una curva en Rn. Las órbitas más sencillas son los puntos de equilibrio, que definimos a continuación, y que se corresponden con las soluciones constantes. Definición 2.5. Sea x0œ. Se dice que x0es un punto crítico o punto de equilibrio para el sistema (2.2) si se tiene f(x0)=0. Los puntos de equilibrio de pueden calcular como la intersección de las curvas de tangencia horizontal y vertical, también llamadas nulclinas o isoclinas cero. 2.4. ESTABILIDAD 9 Definición 2.6. El conjunto ={xœ:fi(x)=0}se denomina nulclina de la variable xi. Observación 2.1. Un punto es crítico si y sólo si pertenece a la intersección de todas las nulclinas. Tras los puntos críticos, el siguiente tipo de órbita más interesante es la órbita periódica, que definimos a continuación. Definición 2.7. Sea x0œ, se dice que la órbita “(x0)del sistema (2.2) es cíclica, periódica o cerrada si I(x0)=Ry existe T>0tal que Ï(t+T;0,x 0)=Ï(t;0,x 0),’tœR. En este caso se dice que Tes un periodo de la órbita. Además, si una órbita periódica se reduce a un punto, se llamará degenerada. A las órbitas periódicas aisladas las denominamos ciclos límite. En un sistema dinámico puede haber dos o más puntos de equilibrio y, por tanto, también es posible que haya órbitas conectando dichos equilibrios. Definición 2.8. A las órbitas que conectan puntos de equilibrio se les llama separatrices, y pueden darse dos casos: 1. Si una órbita conecta dos puntos de equilibrio diferentes, es denominada órbita heteroclina. 2. Si una órbita conecta un punto de equilibrio consigo mismo, se la denomina órbita homoclina. Los equilibrios y las órbitas periódicas son elementos organizadores de la dinámica del sistema. Así, una vez identificados, es importante conocer su estabilidad. Intuitivamente, el concepto de estabilidad nos dice si el resto de órbitas del sistema permanecen cerca o se alejan de ellos, a medida que evoluciona el tiempo. En la siguiente sección, nos centraremos en introducir el concepto de estabilidad. 2.4. Estabilidad Definición 2.9. Un punto de equilibrio del sistema (2.2) se dice estable si para todo ‘>0, existe ”>0tal que, para cualquier xcon Îx≠x0Î<”, la solución Ï(t, x) del sistema (2.2) que pasa por xen t=0satisface la desigualdad ÎÏ(t, x)≠x0Î<‘, para todo tØ0. El punto de equilibrio xse dirá inestable si no es estable. Definición 2.10. Un punto de equilibrio x0se dice que es asintóticamente estable si es estable y, además, existe r>0tal que ÎÏ(t, x)≠x0Îæ0,cuando tæ+Œ, para todo xsatisfaciendo Îx≠x0Î<r. 10 CAPÍTULO 2. HERRAMIENTAS BÁSICAS A continuación, nos centramos en la estabilidad de los puntos de equilibrio en el caso más sencillo, los sistemas lineales. Consideremos el sistema lineal ˙x=Ax, (2.3) donde A=Q c c c c a a11 a12 ··· a1n a21 a22 ··· a2n . . .. . ..... . . an1an2··· ann R d d d d b , con aij œR. Teorema 2.5. Si todos los autovalores de la matriz de coeficientes Aen el sistema lineal (2.3) tienen parte real negativa, entonces el punto de equilibrio x0=0es asintóticamente estable. Más aún, existen constantes positivas Ky–tales que ÎeAtx0ÎÆKe≠–tÎx0Î,’tÆ0,x 0œRn. Si algún autovalor de la matriz Atiene parte real positiva, entonces el punto de equilibrio x0=0es inestable. Es evidente en las definiciones anteriores que el tipo de estabilidad de un punto de equilibrio es una propiedad local. Por consiguiente, es de esperar que bajo ciertas condiciones el tipo de estabilidad de x0se pueda aproximar por el campo vectorial fy sus derivadas, el cual es un campo vectorial lineal. Para ello, supongamos f= (f1,f 2,···,f n)TœC1()y sea la matriz Df(x)= Q c c c c c c c c c c c c c c a df 1 dx1(x)df 1 dx2(x)··· df 1 dxn(x) df 2 dx1(x)df 2 dx2(x)··· df 2 dxn(x) . . .. . ..... . . df n dx1(x)df n dx2(x)··· df n dxn(x) R d d d d d d d d d d d d d d b , la matriz Jacobiana de fen el punto x. Definición 2.11. Si x0es un punto de equilibrio del sistema (2.2) yfœC1, entonces la ecuación diferencial lineal ˙x=Df(x0)x, se denomina ecuación lineal variacional o linealización del campo vectorial fen el punto de equilibrio x0. Ahora enunciamos un resultado fundamental sobre la estabilidad asintótica de un punto de equilibrio a partir de la linealización del campo vectorial. 2.4. ESTABILIDAD 11 Teorema 2.6. Sea fœC1(). Si todos los autovalores de la matriz Jacobiana Df(x0)tienen parte real negativa, entonces el punto de equilibrio x0de la ecuación diferencial (2.2) es asintóticamente estable. Teorema 2.7. Sea fœC1(). Si algún autovalor de la matriz Jacobiana Df(x0) tienen parte real positiva, entonces el punto de equilibrio x0de la ecuación diferencial (2.2) es inestable. En la siguiente definición, señalamos ciertos equilibrios sobre los que podemos esperar capturar la dinámica local de los sistemas no lineales a partir de la linealización. Definición 2.12. Un punto de equilibrio x0del sistema (2.2) se dice hiperbólico si todos los autovalores de la matriz Jacobiana Df(x0)tienen parte real no nula. El siguiente resultado es un resultado central en el estudio de los sistemas dinámicos, que nos permite conocer el comportamiento del sistema en un entorno de los puntos de equilibrio a partir de su linealización, siempre que el equilibrio sea hiperbólico. Teorema 2.8 (Grobman-Hartman).Si x0es un punto de equilibrio hiperbólico de ˙x=f(x), entonces existe un entorno de x0en el cual fes topológicamente equivalente al campo vectorial ˙x=Df(x0)x. Para el estudio de las curvas de respuesta de fase, vamos a considerar sistemas de dos variables. Por ello, nos restringimos, a continuación, al estudio de sistemas planos lineales. Consideremos el sistema lineal (2.3), en el que f:R2≠æ R2. En concreto, consideremos los sistemas de la forma: Y ] [ ˙x1=a11x1+a12x2, ˙x2=a21x1+a22x2,(2.4) o equivalentemente, ˙x=Ax, donde A=Aa11 a12 a21 a22B con aij œR,1Æi, j Æ2. Antes de enunciar un resultado a partir del cual podemos clasificar topológicamente los puntos de equilibrio del sistema (2.4), recordemos algunos conceptos básicos de álgebra: Definición 2.13. Un número (real o complejo) ⁄se denomina autovalor de la matriz Asi existe un vector no nulo (real o complejo) vtal que Av =⁄v. El vector vse denomina autovector de Aasociado al autovalor ⁄. 12 CAPÍTULO 2. HERRAMIENTAS BÁSICAS En esta definición es necesario que el sistema lineal (A≠⁄I)=Osea homogéneo con solución no nula. Esto es equivalente, por el teorema de Rouché Frobenius, a que det(A≠⁄I)=0. Desarrollando esta expresión, nos queda el polinomio cuadrático en ⁄, ⁄2≠·⁄+=0.(2.5) Por tanto, los autovalores de Ason las raíces de este polinomio, que se denomina polinomio característico de A. Los autovalores son: ⁄1=·+Ô·2≠4 2,⁄ 2=·≠Ô·2≠4 2. La naturaleza y la estabilidad de los puntos de equilibrio de un sistema autónomo lineal se pueden describir a partir de la traza (·) y el determinante (), de la matriz A: Lema 2.1. Dados y·, determinante y traza, respectivamente, de la matriz Adel sistema (2.4), tenemos que: Si <0, los autovalores son reales y tienen signos opuestos; los puntos de equilibrio se denominan puntos de silla, y son inestables. Si >0y: •Si ·2≠4>0, los autovalores son reales y tienen el mismo signo; en este caso los puntos de equilibrio se denominan nodos y son: ¶Estables si ·<0. ¶Inestables si ·>0. •Si ·2≠4<0, los autovalores son complejos conjugados; en este caso los puntos de equilibrio se denominan focos o espirales y son: ¶Estables si ·<0. ¶Inestables si ·>0. ¶Centros si ·=0. •Si ·2≠4<0y los autovalores son: ¶Iguales y negativos, el punto de equilibrio es un nodo impropio asintóticamente estable. ¶Iguales y positivos, el punto de equilibrio es un nodo impropio inestable. Si =0, al menos un autovalor es cero. Entonces el punto de equilibrio no es un punto de equilibrio aislado. Hay una recta de puntos fijos, o un plano de puntos fijos, si A=0. Para terminar esta sección, realizaremos un estudio de las órbitas cercanas a una órbita periódica de un sistema autónomo plano y utilizaremos la aplicación de Poincaré para determinar su estabilidad. 2.4. ESTABILIDAD 13 Sea Ï(t, p)una solución periódica con periodo mínimo Tde la ecuación diferencial (2.2) y denotemos la órbita periódica correspondiente por . Escogemos ahora un vector vœR2tal que, vy el vector tangente f(p)de en psean linealmente independientes. De esta forma, podemos definir un segmento: L‘={xœR,x=p+av, 0Æ|a|Æ‘} que se denomina sección transversal de la órbita periódica en el punto p. A continuación, se define una aplicación en un subconjunto de L‘inducida por el flujo. Se elige un ‘tan pequeño que L‘intersecta a la curva en un solo punto p, y que todas las órbitas que crucen L‘lo hagan en la misma dirección. Como Ï(T,p)=py sus soluciones dependen de forma continua del valor inicial, hay un ”>0tal que, si x0œL”,entonces hay un primer instante T(x0)>0en el que Ï(T(x0),x 0)œL‘. Definición 2.14. La aplicación de Poincaré en un entorno de una órbita periódiva se define como: :L”≠æ L‘ x0‘≠æ Ï(T(x0),x 0) Los puntos de la sección transversal L‘tienen un orden natural: dos puntos x0=p+a0vyx1=p+a1vcumplen que x0Øx1en L”, si y solo si, a0Øa1. Por tanto, una aplicación de Poincaré se dice que es monótona si x0Øx1en L” implica que (x0)Ø(x1). Teorema 2.9. La aplicación de Poincaré tiene las siguientes propiedades: La aplicación de Poincaré cerca de la órbita periódica es una aplicación monótona C1. La órbita “(x0)de un punto x0œL”es una órbita periódica si y solo si es un punto fijo de la aplicación de Poincaré, es decir, (x0)=x0. La órbita periódica , con pœ, es asintóticamente estable si Õ(p)<1,e inestable si Õ(p)>1. Concluimos aquí este capítulo en el que hemos introducido herramientas básicas del estudio de las ecuaciones diferenciales. A continuación, nos centraremos en los cambios cualitativos que pueden darse en ecuaciones dependientes de parámetros, llamadas bifurcaciones. 14 CAPÍTULO 2. HERRAMIENTAS BÁSICAS Capítulo 3 Bifurcaciones Antes de estudiar las curvas de respuesta de fase, necesitamos introducir el concepto de bifurcación. Al considerar sistemas dependientes de parámetros, la dinámica de los sistemas puede cambiar a medida que modificamos los valores de los parámetros. Por ejemplo, el sistema puede pasar de tener un equilibrio estable, a uno inestable y un ciclo límite estable. Decimos que esos cambios en la dinámica del sistema son una bifurcación. En este capítulo nos centramos en su estudio. Para este estudio nos hemos apoyado principalmente en las referencias [6], [10] y [15]. Comencemos considerando el sistema dinámico, ˙x=f(x, –),(3.1) donde xœRny–œRmrepresentan las variables, que denominaremos de fase, y los parámetros, respectivamente. Cuando los parámetros varían, la fase también varía. Existen dos posibilidades: o bien el sistema permanece topológicamente equivalente al original, o bien su topología cambia. Definición 3.1. Una bifurcación es un cambio topológico del sistema al pasar sus parámetros por un cierto valor, llamado punto de bifurcación (crítico). Definición 3.2. Llamamos codimensión al número de condiciones independientes que determinan la bifurcación. Nos centraremos en algunas de las bifurcaciones más relevantes de codimensión 1. Consideremos el sistema continuo dependiente de un parámetro: ˙x=f(x, –),xœRn,–œR. Sea x=x0un equilibrio hiperbólico del sistema para –=–0. 15 22 CAPÍTULO 3. BIFURCACIONES Ahora que ya hemos introducido la bifurcaciones más relevantes, en el siguiente capítulo llegamos al objetivo principal del trabajo, y hablaremos sobre las curvas de respuesta de fase. Capítulo 4 Curvas de respuesta de fase (PRC) En este capítulo, comenzaremos introduciendo los conceptos de curva isócrona y de respuesta de fase, los cuales son necesarios para el estudio de las curvas de respuesta de fase para modelos específicos que realizaremos a continuación. En particular, se considerarán los modelos de tipo Integrate-and-Fire que se presentan en el Capítulo 1. Para la elaboración de este capítulo se han consultado las siguientes referencias [2], [6], [17], [18], [19] y [20]. 4.1. Isócronas A lo largo de esta sección daremos la definición de isócrona, así como un ejemplo y algunas aplicaciones. Este concepto resulta de interés, ya que será necesario para entender correctamente la definición de curva de respuesta de fase que estudiaremos a continuación. 4.1.1. Definición Consideremos un sistema autónomo de ecuaciones diferenciales ordinarias ˙x=f(x),xœRn,(nØ2), teniendo una órbita periódica hiperbólica Ï(t)con periodo T. Para cada punto x cerca de la atracción de la órbita periódica existe un único ◊(x)tal que, l´ım tæŒ|x(t)≠Ï(t+◊(x))|=0, donde x(t)es la solución del sistema dinámico anterior, con punto inicial x. Al valor ◊(x)se le denomina fase asintótica de x. Esta noción nos permitirá asignar una fase a cada punto en la base de atracción de una órbita periódica, por tanto, ◊(x)es una función que toma valores en [0,T). 23 24 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) Definición 4.1. La colección de todos los puntos en la cuenca de atracción de Ï con la misma fase asintótica se denomina isócrona. Se puede probar que existen isócronas para cualquier órbita periódica hiperbólica estable. En este caso las isócronas son variedades de codimensión uno, transversales a la órbita periódica Ï. 4.1.2. Ejemplos Este ejemplo, tomado de [17], ilustra algunas propiedades de las isócronas. Supongamos el siguiente sistema, en coordenadas polares: Y ] [ ˙ R=(1≠R)R2, ˙ „=R. (4.1) Este sistema tiene una órbita periódica atractiva Ï={R=1}. Llamemos ◊(„, R) a la función que define la fase asintótica. Las isócronas son conjuntos de curvas unidimensionales con la misma fase asintótica. A continuación vamos a deducir la expresión de la función ◊. Como el campo vectorial es simétrico bajo rotaciones alrededor del origen, la familia de isócronas debe ser invariante bajo esa simetría. Así, asumimos que ◊(„, R)es de la forma ◊(„, R)=„≠f(R),(4.2) para alguna función fa determinar. Entonces cada isócrona es una curva en coordenadas polares. En particular, la isócrona correspondiente a la fase asintótica ces el conjunto de puntos {„, R}que satisface c=„≠f(R). Sustituyendo el valor para el que se tiene la órbita periódica, R=1, en la segunda ecuación de (4.1), vemos que sobre el ciclo límite, ˙ „=1. Además, todos los puntos de una isócrona comparten la misma fase, y por tanto, ˙ ◊(x)=1. Teniendo esto en cuenta, de la expresión (4.2), deducimos d◊ dt =1=d„ dt ≠df dR dR dt . Y despejando de esta expresión obtenemos, df dR =˙ „≠1 ˙ R. Utilizando las expresiones del sistema (4.1) e integrando, obtenemos: f(R)=⁄R≠1 (1 ≠R)R2dR =≠⁄1 R2dR =1 R+C. 4.1. ISÓCRONAS 25 Así, f(R)= 1 R+C, con Cuna constante. Como la órbita periódica se tiene para R=1, tomamos C=≠1, y obtenemos que cada isócrona esta definida por „=c+1 R≠1. En la Figura 4.1, podemos ver ilustradas varias isócronas representativas de la órbita periódica, distinguida por el color rojo. En este ejemplo, las isócronas se definen en todo el plano, excepto en el origen, que es el único elemento del conjunto sin fases. Ahora, definamos una nueva coordenada lejos del origen, =„≠f(R), la cual coincide con la fase asintótica. Derivando con respecto a ty sustituyendo las expresiones de (4.1), tenemos que: d dt =d„ dt ≠df (R) dR ·dR dt =R≠˙ „≠1 ˙ R·(1 ≠R)·R2 =R≠R≠1 (1 ≠R)R2·(1 ≠R)·R2=1. Luego nuestro sistema tendrá la siguiente forma: Y ] [ dR dt =(1≠R)R2, d dt =1. Esto ilustra cómo se puede usar la fase asintótica para encontrar un cambio de coordenadas en el que la dinámica de la coordenada de fase se desacopla de las coordenadas restantes. Figura 4.1: Cuatro isócronas representativas de la órbita periódica (en rojo) considerada en el ejemplo. 26 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) A continuación, enumeraremos algunas aplicaciones del estudio de las isócronas. 4.1.3. Aplicaciones Las isócronas tienen diversas aplicaciones. Como se ha señalado en el ejemplo anterior, pueden utilizarse para reducir eficazmente la dimensión de la ecuación en la proximidad de una órbita periódica. También son útiles para: Extender la noción de fase de una órbita periódica en un entorno de la órbita periódica. Es habitual tomar la fase como un avance lineal en el tiempo, con valores en el intervalo [0,T),[0,2fi),oen[0,1). Definición de las curvas de respuesta de fase de los osciladores. Estudio del tiempo que pasa una trayectoria en diferentes regiones del espacio de fases. Con estas aplicaciones finalizamos el estudio de las isócronas. En la siguiente sección, introduciremos las curvas de respuesta de fase, con el objetivo de desarrollar, posteriormente, las curvas de respuesta de fase para modelos específicos. Entre otros, consideraremos los modelos de tipo Integrate-and-Fire presentados en el Capítulo 1. 4.2. Curvas de Respuesta de Fase (PRC) Para comenzar esta sección empezaremos introduciendo el concepto de las Curvas de Respuesta de Fase (PRC, por sus siglas en inglés) y posteriormente las clasificaremos. 4.2.1. Definición y metodología Las PRC describen cómo las neuronas responden a un estímulo. Es decir, una PRC es una curva que describe la relación entre un estímulo, como la exposición a la luz, y una respuesta, como un cambio en el ritmo cardíaco. Estimulando la neurona en diferentes fases, podemos medir la PRC, de la siguiente forma, PRC(◊)=◊n≠◊,(4.3) siendo ◊nla nueva fase después del estímulo y ◊la fase inicial. Lo podemos ver representado en la Figura 4.2. Suponemos que el período de funcionamiento libre y sin perturbaciones del oscilador es T0. El período de la oscilación se puede normalizar a2fi. Por lo tanto, cada punto de la oscilación se puede describir de forma única con una fase ◊œ(0,2fi]. Elegiendo un punto de referencia arbitrario en un momento t, podemos aplicar una perturbación arbitraria en una fase ◊=ts/T0, 4.2. CURVAS DE RESPUESTA DE FASE (PRC) 27 donde tses el tiempo transcurrido desde la fase de referencia. El ciclo que contiene la perturbación tiene una longitud T1(véase Figura 4.2). Figura 4.2: Voltaje de una neurona que se dispara periódicamente (en verde). Voltaje perturbado (en rojo). Definimos el reinicio de fase como, ◊=(T0≠T1)/T0. Si se acorta el período, se supone que la trayectoria se ha desplazado en la dirección del movimiento, provocando un incremento en la fase llamado avance de fase.Por otro lado, si el período se alarga, se supone que la trayectoria se ha desplazado en una dirección opuesta a la dirección del movimiento, lo que provoca una disminución en la fase, llamada retraso de fase. La teoría de los osciladores neuronales débilmente acoplados requiere el concep- to de PRC infinitesimal (iPRC), que es matemáticamente equivalente a la derivada parcial de la fase con respecto al voltaje (d◊/dV ), ya que generalmente solo se consideran perturbaciones en el voltaje. A continuación, clasificaremos dos tipos de PRCs, los cuales los utilizaremos para calcular las PRC para modelos específicos. 28 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) 4.3. Ecuaciones de fase para osciladores no lineales con ciclos límite atractivos 4.3.1. Reducciones de fase En primer lugar describimos un cambio de coordenadas a las variables de fase que simplificará el análisis que desarrollaremos. Nuestro punto de partida es un modelo general basado en la conducta de una sola neurona, descrito por las ecuaciones (1.1) en el Capítulo 1. Recordemos que, Ves la diferencia de voltaje a través de la membrana y nes un vector de dimensión N≠1. Además Ig(V,n)es la corriente de membrana asociada al vector n,yCla capacidad de la membrana. Por otro lado, Ibestablece la frecuencia del oscilador y corresponderá, a continuación, a un parámetro de bifurcación.La función I(V,t)representa las corrientes sinápticas procedentes de otras áreas cerebrales debido a la presentación de estímulos; a continuación, consideraremos I(V,t)=I(t). Como vimos, estas ecuaciones en la forma general venían dadas por la ecuación (1.2), donde F(x)es el campo vectorial de base yG(x, t)es el efecto del estímulo. Finalmente, recordemos que habíamos tomado, G(x, t)=I(I(t),0)T. Asumamos ahora que el oscilador neuronal de referencia (G©0) tiene un ciclo límite normalmente hiperbólico y atractivo, que denotaremos por “. Esta configuración se mantiene bajo pequeñas perturbaciones, y a partir de ahora supondremos que dicho ciclo límite siempre existe para cada neurona. El objetivo que nos planteamos es simplificar la ecuación (1.2) definiendo una variable escalar de fase ◊(x)œ[0,2fi)para todo xen alguna región Udel plano, cercana a“, tal que d◊(x) dt =Ê para todo xœUcuando G©0. Aquí, tomaremos Ê=2fi/T, (4.4) donde Tes el periodo de la ecuación (1.2) con G©0. Aplicando la regla de la cadena a la ecuación (1.2), se sigue que, d◊(x) dt =ˆ◊ ˆx(x)·F(x)+ˆ◊ ˆx(x)·G(x, t)=Ê+ˆ◊ ˆx(x)·G(x, t).(4.5) Esta ecuación define una ecuación en derivadas parciales de primer orden que debe satisfacer ◊(·). Evaluando el campo vectorial en el campo base x“(◊), que definimos como la intersección de “y las curvas de nivel ◊(x), a partir de (4.5) tenemos que, d◊(x) dt =Ê+ˆ◊ ˆx(x“(◊)) ·G(x“(◊,t)) + E, (4.6) donde Ees un error de O(|G|2), y donde el escalar |G|acota G(x, t)sobre todas las componentes, así como sobre xyt. 4.3. ECUACIONES DE FASE 29 Excluyendo este término de error, podemos reescribir la ecuación (4.6) como la ecuación de fase en dimensión uno, d◊ dt =Ê+ˆ◊ ˆx(◊)·G(◊,t),(4.7) la cual es válida en toda región Ucercana a “. 4.3.2. Cálculo de la PRC. En el caso de las ecuaciones del sistema (1.1), la única derivada parcial que debemos calcular para definir completamente la ecuación (4.7) es con respecto al voltaje, y definimos la curva de respuesta de fase como, ˆ◊ ˆV(◊)©z(◊).(4.8) Así, obtenemos, d◊ dt =Ê+z(◊)I(t)©v(◊,t).(4.9) A continuación, describimos el método directo para calcular la función z(◊),el cual usaremos en la siguiente sección para calcular las PRCs de las formas normales de la bifurcaciones que suelen aparecer en los modelos neuronales, y que hemos analizado en el Capítulo 3. Método directo: Este método es una forma sencilla y clásica de calcular z(◊)que es útil en estudios experimentales, numéricos y analíticos. Por definición, z(◊) = l´ım Væ0 ◊ V,(4.10) donde ◊=[◊(x“+(V,0)T)≠◊(x“)], es el cambio en ◊(x)resultante de la perturbación VæV+V del punto base x“en “(ver Figura 4.3). Ya que ˙ ◊=Êen toda vecindad de “, la diferencia ◊se preserva bajo el flujo de la fase base (G©0); por lo tanto, se puede medir en el límite cuando tæŒ, cuando la trayectoria perturbada ha regresado al ciclo límite “. Es decir, z(◊)se puede encontrar comparando las fases de las soluciones en el límite infinito de tiempo, partiendo de puntos que están a una distancia infinitesimal de los puntos base en el ciclo límite. 30 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) Figura 4.3: El método directo para calcular ˆ◊ ˆVen el punto indicado por ú, es tomar el límite de ◊/Vpara perturbaciones pequeñas en el voltaje, V. 4.4. PRCs para modelos específicos En esta última sección del capítulo y del trabajo proporcionaremos aproximaciones analíticas a las PRCs para sistemas multidimensionales con ciclos límite que surgen de las cuatro bifurcaciones de codimensión uno descritas en el Capítulo 2, las cuales son apropiadas para los modelos basados en la conductancia, como las ecuaciones (1.1). 4.4.1. Curvas de respuesta de fase en un entorno de bifurcaciones de codimensión uno con impulsos periódicos Tal y como vimos en el Capítulo 2, la teoría de bifurcaciones identifica cuatro bifurcaciones de codimensión uno, las cuales han sido estudiadas en dicho capítulo, en las cuales está involucrado un ciclo límite. Los cuatro tipos de bifurcación han sido identificados en modelos neuronales específicos con un parámetro. Aquí, la corriente de entrada de base es Ib. En esta sección, calculamos las PRC para los ciclos límite que surgen de las cuatro bifurcaciones. Esto se logra, cuando es posible, mediante el uso de ecuaciones en forma normal de una y dos dimensiones. Las formas normales se obtienen a través de cambios de variable, con el objetivo de eliminar tantos términos como sea posible, en un proceso que preserva la dinámica cualitativa del sistema. Para obtener la PRC en términos de las variables originales, es decir, ˆ◊ dV , en lugar de en términos de las variables de la forma normal con las PRCs asociadas ˆ◊ ˆxyˆ◊ dy , es necesario deshacer estas transformaciones de coordenadas. Sin embargo, como las transformaciones de coordenadas de la forma normal sólo afectan a los términos no lineales, obtenemos 4.4. PRCS PARA MODELOS ESPECÍFICOS 31 la sencilla relación, ˆ◊ ˆV=‹x ˆ◊ ˆx+‹y ˆ◊ ˆy+O(x, y),(4.11) donde ‹x=ˆx ˆV|x=y=0,‹ y=ˆy ˆV|x=y=0. El término restante en la ecuación (4.11) se supone que es pequeño cerca de la bifurcación de relevancia y se desprecia a continuación. Utilizamos el método directo para calcular los PRCs a partir de las ecuaciones en forma normal. Esto implica linealizar sobre la órbita periódica estable, lo cual es apropiado porque las perturbaciones Va considerar son pequeñas. La solución explícita de las ecuaciones en forma normal arroja ◊, y tomando límites, obtenemos la PRC. Sin pérdida de generalidad, la fase del voltaje (spike) se fija en ◊s=0, y se definen las coordenadas para que la fase aumente a un ritmo constante Êen ausencia de estímulos externos. Aquí calculamos explícitamente cómo varían los PRCs con la frecuencia del oscilador. Comenzamos con el caso de la bifurcación SNIPER. Bifurcación silla-nodo en una órbita periódica (SNIPER) Como comentamos en el Capítulo 2, la bifurcación SNIPER ocurre cuando tiene lugar una bifurcación silla-nodo de puntos de equilibrio sobre una órbita periódica. Al igual que en ese capítulo, consideremos el sistema: ˙x=÷+x2.(4.12) Para ÷>0la solución de la ecuación recorre un intervalo cualquiera en un tiempo finito; el período Tde la órbita puede aproximarse calculando el tiempo total necesario para que la solución de la ecuación (4.12) pase de x=≠Œ ax=+Œ y haciendo que la solución sea periódica restableciendo xa≠Œ cada vez que se dispara en x=Œ. Esto nos proporciona T=fi Ô÷, por lo tanto, sustituyendo esta expresión en la ecuación (4.4), obtenemos que Ê=2 Ô÷. Teniendo en cuenta la expresión (4.7) y recordando que habíamos supuesto G©0, se obtiene: ˆ◊ ˆx=ʈt ˆx=Ê ˆx ˆt ,(4.13) 38 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) Podemos ver representada la PRC, zH(◊), en la Figura 4.6. Figura 4.6: Aproximación analítica de la PRC, zH(◊). Las curvas azul, roja y amarilla, son las PRCs correspondientes resultantes de incrementar la frecuencia Ê, respectivamente: Ê=0,204,0,312,0,414. Para el caso de la bifurcación Bautin, obtenemos de manera similar que zB(◊)= d◊ dV =cB |Ê≠ÊSN|sen (Ê≠„B),(4.42) donde cB=S U≠2dÛ≠c 2f≠4gA≠c 2fB3/2T V f ckÒ‹2 x+‹2 y es una constante cuyo signo dependerá de fyg,y „B=tan ≠1A‹x ‹yB. Podemos ver representada la PRC, zB(◊), en la Figura 4.7 4.4. PRCS PARA MODELOS ESPECÍFICOS 39 Figura 4.7: Aproximación analítica de la PRC, zB(◊). Las curvas azul, roja y amarilla, son las PRCs resultantes de incrementar la frecuencia Ê, respectivamente: Ê=0,339,0,439,0,539. Así, se concluye el estudio para el caso de Hopf en los casos subcrítico y supercrítico. A continuación procedemos al caso de la bifurcación homoclínica. Bifurcación homoclínica Por último, supongamos que el modelo neuronal tiene un parámetro µtal que existe una órbita homoclínica a un punto de silla hiperbólico pcon autovalores reales en µ=0. Entonces habrá una órbita periódica “para µ>0, pero no para µ<0. Específicamente, asumimos un único autovalor inestable ⁄ude magnitud menor que el de todos los valores propios estables, ⁄u<|⁄s,j|, para que la órbita periódica bifurcada sea estable. Si los parámetros se eligen cerca de la bifurcación homoclínica, las soluciones cerca de la órbita periódica pasan la mayor parte del tiempo cerca de p, donde el campo vectorial está dominado por su linealización. Esto puede escribirse genéricamente en la forma diagonal: Y ] [ ˙x=⁄ux, ˙yj=⁄s,j yj,j=1,···,N ≠1,(4.43) donde los ejes xeyjson tangentes a la región inestable y a la estable de p, respectivamente, y ⁄s,j <0<⁄ uson los valores propios correspondientes. Para simplificar, suponemos aquí que los segmentos de los ejes mostrados en la Figura 4.8 están realmente contenidos en las respectivas regiones; esto siempre se puede conseguir localmente mediante un cambio de coordenadas regular. 40 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) Definamos la caja B=[0,]◊···◊[0,]que encierra a “durante la parte dominante de su período, pero dentro de la cual las expresiones de (4.43) siguen siendo una buena aproximación. No modelamos explícitamente “fuera de B,pero observamos que la trayectoria se reinyecta después de un tiempo despreciable (comparado con el transcurrido en B) a distancia ‘del conjunto, donde Xvaría con el parámetro de bifurcación µ(véase Figura 4.8). Así, las órbitas periódicas que se producen más cerca del punto de bifurcación corresponden a valores menores de ‘y tienen periodos mayores. Figura 4.8: Configuración para derivar el PRC para las oscilaciones cerca de una bifurcación homoclínica, mostrada con N=2. Reprentamos con el color azul la órbita periódica “. Aproximamos el periodo T(‘)como el tiempo que la coordenada xde “tarda en viajar desde ‘hasta según la primera expresión de (4.43): T(‘)= 1 ⁄u lnA ‘B.(4.44) Obsérvese que la coordenada xde “determina por sí sola T(‘). Fijamos ◊=0 en x=‘, y suponemos una reinyección instantánea, ◊=2fien x=. Entonces Ê=2fi T(‘), y como en la ecuacion (4.13), ˆ◊ ˆx=Ê dx dt =Ê ⁄ux(◊)=Ê ⁄u‘exp(≠⁄◊/Ê).(4.45) En la última igualdad, hemos usado la solución de la primera expresión del sistema (4.43), x(t)=‘exp(⁄ut),cont=◊ Ê. Dado que, como se ha comentado anteriormente, el movimiento en las direcciones yjno afecta a la fase de “, sólo los componentes de una perturbación Va lo largo del eje xcontribuyen a la curva de respuesta de fase; así, el PRC viene dado por zHC =ˆ◊ ˆV=‹x ˆ◊ ˆx,(4.46) donde ‹xes tal y como se ha definido en la ecuación (4.11). Usando la expresión (4.44), y teniendo en cuenta que T(‘)=2fi Ê, 4.4. PRCS PARA MODELOS ESPECÍFICOS 41 se obtiene, fácilmente, ‘=exp(≠2fi⁄u/Ê). Sustituyendo el valor de ‘en la ecuación (4.45), obtenemos, ˆ◊ ˆx=Ê ⁄uexpA2fi⁄u ÊBexpA≠⁄◊ ÊB Finalmente, sustituyendo este valor en (4.46), concluimos que, zHC(◊)= ˆ◊ ˆV=chcÊexpA2fi⁄u ÊBexpA≠⁄u ◊ ÊB,(4.47) donde chc =‹x ⁄ues una constante dependiente del modelo. Esta es una función exponencialmente decreciente de ◊con su máximo en: zm´ax =chcÊexpA2fi⁄u ÊB(4.48) y mínimo en zm´ın =zm´ax expA≠2fi⁄u ÊB=chcÊ.(4.49) Podemos ver representada la PRC para este modelo en la figura 4.9. Figura 4.9: Aproximación analítica de la PRC, para el modelo homoclínico. Las curvas azul, roja y amarilla, son las PRCs correspondientes a este modelo resultantes de incrementar la frecuencia Ê, respectivamente: Ê=0,0572,0,0972,0,0871. A partir de ahora tomaremos chc >0. Nótese que zHC es discontinua en el punto ◊s=2fi, lo que nos obliga a tomar un límite en la definición de las tasas de disparo promediadas por la población a continuación, pero esto no afecta por lo demás al siguiente análisis. 42 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) 4.4.2. Modelos neuronales de dimensión uno Los modelos integrate-and-fire generalizados son de la forma ˙ V=F(V)+G(V,t),(4.50) donde V(t)está limitado entre un voltaje de inicialización Vry un umbral Vth. Imponemos la siguiente dinámica de reajuste: si V(t)cruza Vth desde abajo, se produce un spike, y V(t)se reajusta a Vr. Aquí no se pierde nada en la transformación de la ecuación (4.7). En particular, el término de error de la ecuación (4.6) no se aplica. De hecho, como se señala en [7], ˆ◊ ˆVse puede obtener directamente de la ecuación (4.50) con G(V,t)©0: z(◊)= ˆ◊ ˆV=ʈt ˆV=Ê F(V(◊)),(4.51) donde ◊está definido de modo que ˙ ◊=Ê. A continuación, calculamos las PRCs para los dos modelos de Integrate-and-Fire que se presentaron en el Capítulo 1. Neurona Integrate-and-Fire (IF) En primer lugar, consideremos el caso más simple de modelo Integrate-and-Fire: C˙ V=Ib+I(t),(4.52) con Vr=0, y Vth =1, donde Ibes la corriente base, Ces la capacidad de la membrana, y G(V,t)=I(t). En adelante, fijaremos C=1para el modelo IF. La frecuencia angular de I(t)=0 es Ê=2fiIb, y de la ecuación (4.51) se tiene: zIF =Ê F(V(◊)) =Ê Ib©2fi.(4.53) Por lo tanto, la PRC para el modelo IF es constante en Êe independiente de la frecuencia. Podemos verla representada en la Figura 4.10. Figura 4.10: Aproximación analítica de la PRC para el modelo IF. 4.4. PRCS PARA MODELOS ESPECÍFICOS 43 Neurona Leaky IF (LIF) Ahora, consideramos el modelo LIF: C˙ V=Ib+gL(VL≠V)+I(t)(4.54) con Vr=0, y Vth =1<V L+VIb gL , donde Ibes la corriente base, gL>0yVLson la conducta de fuga y el potencial de reacción, Ces la capacidad, y G(V,t)=I(t). Al igual que antes, para este modelo fijamos C=1. Suponemos IbØgL(1 ≠VL)y tomamos I(t)=0,así ˙ V=dV dt =Ib+gL(VL≠V) por tanto, separando las variables e integrando ⁄1 gL(VL≠V)+Ib dV =⁄dt, es decir, ≠log (gL(VL≠V)+Ib) gL +c=t, (4.55) con cuna constante de integración. Teniendo en cuenta que V(0) = 0, se obtiene que, c=ln(gLVL+Ib) gL . Ahora, teniendo en cuenta que V(T)=1 ≠log (gL(VL≠1) + Ib)=AT≠ln(gLVL+Ib) gLBgL, obteniéndose así el periodo T: T= lnAIb+gLVL Ib+gLVL≠gLBg≠1 L. De este modo cuando I(t)=0, teniendo en cuenta que Ê=2fi T, la neurona se dispara periódicamente con frecuencia Ê=2figLCln AIb+gLVL Ib+gLVL≠gLBD≠1 .(4.56) Esta expresión muestra cómo Ibinterviene como parámetro de bifurcación, con Ib= gL(1 ≠VL)correspondiente al punto de bifurcación en el que Ê=0. Resolviendo la ecuación (4.54) para V(t)con la condición inicial V(0) = Vr=0, 44 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) y luego usando que ◊=Êty la expresión (4.51), se obtiene: zLIF (◊)= Ê gL31≠exp3≠2figL Ê44expAgL◊ ÊB.(4.57) Por lo tanto, la PRC para el modelo LIF es una función expenencialmente creciente de ◊, con un máximo que disminuye con Ê, zmax(Ê)= Ê gL3exp32figL Ê4≠14,(4.58) y mínimo zmin(Ê)=zmax exp3≠2figL Ê4=Ê gL31≠exp3≠2figL Ê44.(4.59) En la Figura 4.11 representamos zLIF (◊). Figura 4.11: Aproximación analítica de la PRC, para el modelo LIF. Las curvas azul, roja y amarilla, son las PRCs resultantes de incrementar la frecuencia Ê. LIF: Ê=0,419,0,539,0,639.gL=0,110. Recordemos que la PRC cerca de una bifurcación homoclínica es también una función exponencial pero con pendiente opuesta: esto se debe a que tanto la dinámica principal cerca de una bifurcación homoclina como la dinámica LIF son lineales, mientras que las trayectorias se aceleran en el caso homoclínico y se desaceleran en el LIF. Este es nuestro cálculo analítico final de la PRC. Resumimos los resultados obtenidos anteriormente en la tabla. 4.4. PRCS PARA MODELOS ESPECÍFICOS 45 Bifurcación z(◊)zmax zmin SNIPER csn Ê[1 ≠cos(◊)] 2csn Ê0 Hopf cH Ô|Ê≠ÊH|[sen(◊≠„H)] cH Ô|Ê≠ÊH|≠cH Ô|Ê≠ÊH| Bautin |cB| |Ê≠ÊSN |[sen(◊≠„B)] + O(1) |cB| |Ê≠ÊSN |+O(1) ≠|cB| |Ê≠ÊSN |+O(1) Homoclina chcÊexp 12fi⁄u Ê2exp(≠⁄u◊/Ê)chcÊexp12fi⁄u Ê2chcÊ IF 2fi2fi2fi LIF Ê gL11≠e≠2figL/Ê2egL◊/ÊÊ gL1e≠2figL/Ê≠12Ê gL11≠e≠2figL/Ê2 4.4.3. Precisión de las PRCs analíticas Las bifurcaciones SNIPER, Hopf, Bautin y Homoclina aparecen en modelos neuronales más complejos que los presentados en esta trabajo. En particular, la bifurcación SNIPER aparece en el modelo de Hindmarsch-Rose, la Hopf supercrítica en el modelo de FitzHugh-Nagumo, la Bautin en el de Hodgkin-Huxley, y la homoclina en el de Morris Lecar. Estos modelos neuronales se pueden encontrar en [9], [11], [13] y [14]. A modo de conclusión, debemos remarcar que el rango de parámetros en el que las PRCs de los modelos neuronales se aproximan bien mediante las expresiones analíticas, derivadas anteriormente, varía de un modelo a otro. Una limitación general es que los cálculos de la forma normal para la bifurcación de Bautin y Hopf supercrítica ignoran la naturaleza de relajación de la dinámica de los osciladores neuronales clásicos. Sin embargo, las PRC analíticas calculadas - ecuaciones (4.15), (4.41), (4.42) y (4.47) - son cualitativamente, correctas, y en muchos casos, cuantitativamente también (en las figuras 4.4, 4.6, 4.7, 4.9, 4.10 y 4.11, están representadas estas aproximaciones analíticas). Finalmente, como se puede observar en el artículo [3], los PRCs analíticos derivados aquí predecirán correctamente aspectos cualitativos clave de las respuestas de la población neuronal a los estímulos. 46 CAPÍTULO 4. CURVAS DE RESPUESTA DE FASE (PRC) Bibliografía [1] A. L. Hodgkin y A. F. Huxley, Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo. The Journal of physiology, 1952. [2] E. Brown, J. Moehlis, and P. Holmes, On the Phase Reduction and Response Dynamics of Neural Oscillator Populations, Neural Computation, 2004. [3] Departamento de Ecuaciones Diferenciales y Análisis Numérico, Apuntes de la Asignatura EDO, 2017/18, Universidad de Sevilla. [4] Departamento de Ecuaciones Diferenciales y Análisis Numérico, Apuntes de la Asignatura AED, 2019/2020, Universidad de Sevilla. [5] Departamento de Ecuaciones Diferenciales y Análisis Numérico, Apuntes de la Asignatura ANED, 2020/2021, Universidad de Sevilla. [6] Ermentrout, G. Type I membranes, phase resetting curves, and synchrony. Neural Computation, 1996. [7] Ermentrout, G., Phase locking of weakly coupled oscillators. J. Math. Biol., 1981. [8] J. Hale y H. Koçak, Dynamics and Bifurcations, New York: Springer-Verlag, 1996. [9] Hodgkin, A. L. y Huxley, A. F. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiology, 1952. [10] Kuznetsov, Yuri A., Elements of applied bifurcation theory, New York: Springer-Verlag, 2004. [11] Murray, J. Mathematical biology I: An introduction (3rd ed). New York: Springer-Verlag, 2002. [12] E. Pérez Castilla, Simulación numérica de modelos de ecuaciones diferenciales en Neurociencia., Trabajo de fin de Grado, Universidad de Sevilla, 2020. [13] Rinzel, J. and Ermentrout, G.B. Analysis of Neural Excitability and Oscillations MIT Press, Cambridge, 1989. [14] R. M. Rose and J. L. Hindmarsh The Assembly of Ionic Currents in a Thalamic Neuron I. The Three-Dimensional Models Royal Society, 1989. 47