Repositorio Institucional de Documentos

Modelos Matemáticos de Neuronas Nuria Begué Pedrosa Trabajo de fin del grado de Matemáticas Universidad de Zaragoza Agradecimientos Me llena de satisfacción poder llegar a escribir estas palabras. En primer lugar, agradezco a mi director, Roberto Barrio, su generosidad al haber compartido sus conocimientos conmigo. Y, mostrarse siempre dispuesto a ayudarme con cada una de mis dudas y aconsejarme en todo momento. Por otro lado, voy a aprovechar este escrito para dar las gracias a todas aquellas personas que durante este periodo de mi vida me han dado su apoyo incondicional confiando en mi plenamente. Finalmente, con la siguiente cita quiero alabar ese ansia de la ciencia y en particular, de las Matemáticas, de abrir nuevos horizontes sobre lo sombrío y dudoso, llenando de luz todas nuestras incógnitas. Y, elogiar, desde mi humildad, a todas aquellas personas que lo han realizado a lo largo de sus vidas. Lo que sabemos es una gota de agua; lo que ignoramos es un océano. Isaac Newton. III Índice general Agradecimientos III Summary VII 0.1. Hodgkin and Huxley Mathematical Neuronal Model ................. VII 0.2. Hindmarsh and Rose Mathematical Neuronal Model ................. XI 0.2.1. The 1982 Model ................................ XI 0.2.2. The 1984 Model ................................ XIII 0.3. Return to the reality ................................... XVI 1. Introducción 1 2. Modelo matemático neuronal de Hodgkin y Huxley 3 2.1. Objetivo a modelar: Acción potencial ......................... 3 2.2. Modelo matemático neuronal de Hodgkin y Huxley .................. 4 2.3. Circuito eléctrico equivalente .............................. 5 2.3.1. La conductancia del Potasio .......................... 7 2.3.2. La conductancia del Sodio ........................... 9 2.3.3. Obtención de las constantes de proporción ................... 10 2.3.4. Modelo completo ................................ 11 3. Estudio cualitativo de modelos matemáticos de neuronas 13 3.1. Modelo Matemático Neuronal de FitzHugh y Nagumo ................ 15 3.2. Modelo Matemático Neuronal de Hindmarsh y Rose ................. 15 3.3. Conceptos de la Teoría de Bifurcaciones ........................ 16 3.3.1. Bifurcación fold o silla-nodo .......................... 17 3.3.2. Bifurcación Andronov-Hopf .......................... 19 3.4. Modelo Matemático Neuronal de Hindmarsh y Rose (1982) .............. 22 3.5. Modelo Matemático Neuronal de Hindmarsh y Rose (1984) .............. 26 3.6. Vuelta a la realidad ................................... 33 4. APÉNDICE: Técnicas y herramientas utilizadas. 35 Bibliografía 37 V Summary 0.1. Hodgkin and Huxley Mathematical Neuronal Model All living cells have an electric voltage, or potential difference, between their inside and outside which is called membrane potential. At first, the cells have a resting potential but the membrane potential will change when the flux of ions acrosses the cell membrane, for example during an action potential. This is feasible due to some gated channels are closed or open. In the case of neuron, applying some stimuli will make the neuron excitable. If this stimuli is large enough such that a threshold value, the neuron will generate an action potential and will start to fire. On the other hand, if the stimuli is not enough, the neuron will return to its resting state. The action potential generate signals that are the core of transmission of nervous impulse. One of the most outstanding model in mathematical neuroscience is the Hodgkin and Huxley model for the propagation of an action potential along the squid’s giant axon. This choice was really useful thanks to the large diameter of the axon because tools, like electric microscope, and numerical tools did not exist at this time. Hodgkin and Huxley developed two experimental methods to separate the ion currents and compute how ion currents, mainly K+and Na+, depend on voltage. The real predictive power of the model became evident when these differencial equations could accurately reproduce the key of the biophysical properties of the action portential. For this important achievement, Hodgkin and Huxley were awared in 1963 with the Nobel Prize in Physiology and Medicine. Nowadays, we are interested in understanding this superb modelling exercise. Firstly, we can notice the electric properties of a neuron which allow us to understand the neuron as an electric circuit. This is a very useful way to describe the behavior of the membrane potential; which is commonly called the equivalent circuit model, see figure 1. The circuit consists on three components: conductors or resistors, representing the ion channels; batteries are the concentration gradients of the ions and a capacitor, showing the ability of the membrane to store charge. Thus, Kirchoff’s laws were applied to modelate the behavior of the potencial membrane, Vm. This change is described by the follow differential equation: CmdVm dt +Iion =0.(1) If an external current is applied, then sum (1) will equal the external or applied current, giving us, CmdVm dt =Iext −Iion (2) The current is related to the membrane voltage through Ohm’s law relationship of the form V=IR. Nevertheless, we consider this relation in terms of conductance rather than resistance, where the conductance gis the inverse of resistance, g=1/R. In relation to ion channels, the equilibrium potential or Nernst’s potential, Vi, must be introduced for each ion. This is the potential at which the current flowing to bring closer zero. The current is the product of the conductance and the difference between the membrane potencial, Vm, and the Nernst’s equilibrium, Vi. The total ionic current, Iion, is the algebraic sum of the individual contributions from all participating channels types, that is: VII VIII Capítulo 0. Summary Figura 1: Equivalent circuit model. Iion =INa +IK+IL=gNa(Vm−VNa) +gK(Vm−VK) +gL(Vm−VL).(3) Besides, the other aim is how the gates concern the model. Thus, Hodgkin and Huxley defined pias the fraction of gates are open, then 1−piis the fraction of the closed gates. This probability depends on the membrane potential. The rate at which gates transition from the non-permissive state to the permissive state is denoted by a variable αi(V), and βi(V)that represents the opposite transition. From the law of mass action, dpi dt =αi(V)(1−pi)−βi(V)pi.(4) If Vmis fixed as V, then the equation (4) will reach a stable value, that is, dpi dt =0 when t−→∞, giving us, pi∞(V) = αi(V) αi(V)+βi(V).(5) The course of the time for reaching this state of equilibrium is described by, τi(V) = 1 αi(V)+βi(V).(6) At this point, we wonder about which are the differences between IKconductances and INa ones. The response was given thanks to voltage-clamp. Hogkin and Huxley were able to isolate K+or Na+, respectively, to prove the fast activation and slow inactivation gates. On the one hand, they worked out IKconductance depends on a activacion variable, n. Thereby, the differential equation governing the state variable nis given by, dn dt =αn(V)(1−n)−βn(V)n.(7) Hodgkin-Huxley realized a better fit could be obtained if they considered the conductance to be proportional to a higher power of n. They concluded that a reasonable fit to the K+conductance data could be using an exponent r=4. Thus they arrived at a description for the K+conductance under voltage-clamp conditions given by, Modelos Matemáticos de Neuronas 0.1. Hodgkin and Huxley Mathematical Neuronal Model IX Figura 2: Activation and inactivation variables as solution of the Hodgkin-Huxley model, where green is n(t), activation of K, and m(t), red, represents activation gate of Na and blue is h(t)or inactivation of K. gK=gKn4.(8) where gKis a constant of standardisation which determines the maximum conductance when all gates are open. On the other hand, the strategy for modelling the sodium conductance is similar to the described above for the potassium conductance, except that the sodium conductance shows a more complex performance. The difference resides in the fact that sodium channels present two dynamics such as inactivation and activacion. For modelling this process, Hodgkin and Huxley postulated that the sodium had two types of gates, an activation gate, which they labeled m, and an inactivation gate, which they labeled h. dm dt =αm(V)(1−m)−βm(V)m,(9) dh dt =αh(V)(1−n)−βh(V)h.(10) Hodgkin and Huxley were able to fit the remaining parameters from the voltage clamp data. The sodium conductance gNa was modeled by an expression of the form, gNa =gNam3h.(11) In summary, the Hodgkin-Huxley model is a system of four differetianl equations; there is one equation for the membrane potential and three equations for the channel gating variables. In the case of a space-clamp squid axon. We must point out the existence of a leak current, IL, compound of ion Cl−which does not depend on the time. Thus, the model is giving by, Autor: Nuria Begué Pedrosa XVI Capítulo 0. Summary coincide with the stable lymit cicle. In other words, it happens when a lymit cicle hits a saddle-point, then the limit cycle vanishes. The other fact noticing from the homoclinic bifurcation theory is that the period of the lymit cycle grows logarithmically fast. This explains the increase of the interspike intervals by the end of bursts. We can realize that this bifurcation occurs before the inverse AH, this is the reason that the final of the burst presents this new configuration. The study of this bifurcation is very complicated due to this one is global, while the previous are local. At z=83/54, the curve has a fold bifurcation . Therefore, the dynamic system goes to the stable limit cycle. Thus, the slow motion finalizes. Therefore, the fold bifurcation and the homoclinic one commence and end the fast motion of the burst, respectively. 0.3. Return to the reality The qualitative study of these models have given us interesting conclusions about the comprehension of the neuronal dynamic. The bifurcations explain each sort of the neuronal activity during an action potential. Nevertheless, we would like to check if the models are real. If we analize the figure 10, we can notice that in the line athe figures FRB and IB are examples of the neurons who represent square wave bursting during an action potential. The more common neurons of the cortex describe a dynamic based on regular spiking. Figura 10: Real neuronal dynamic as an action potential. Modelos Matemáticos de Neuronas Capítulo 1 Introducción La introducción de la matemática en modelos biológicos nos permite una abstracción y comprensión más extensa de dichos procesos. Actualmente, la introducción de la matemática en el campo de la neurociencia tiene mucha fuerza. El primer capítulo tiene por objeto mostrar uno de los trabajos de modelización matemática más éxitosos tanto por su reconocimiento como por su aplicación, desarrollado por Hodgkin y Huxley a mediados del siglo pasado. Con el cual se pretende aplicar el trabajo de modelización matemática en el estudio del impulso nervioso de una neurona, el cual descansa sobre la acción potencial generada por movimientos de cargas a través de la membrana de la neurona. Luego, el modelo de HodgkinHuxley constituye todo un ejercicio multidisciplinar gracias al cual se pudo entender la generación y propagación del impulso nervioso. Se observa en el lenguaje la estrecha relación con la teoría de circuitos. Por tanto, no nos ha de extrañar que partiendo de conceptos propios de los circuitos llegamos a la expresión matemática de la acción potencial, en forma de un sistema de ecuaciones diferenciales no lineales. Sin embargo, la expresión matemática de la misma hace que el estudio cualitativo del sistema resulte complejo llevando a la necesidad de modelos más sencillos que recojan toda la información de la dinámica neuronal. En el siguiente capítulo se expondremos modelos matemáticos neuronales más sencillos en cuanto a estructura. La diferencia relevante con el capítulo anterior radica en que vamos a pensar la neurona como un sistema dinámico no lineal. Este nuevo enfoque nos permite introducir una serie de conceptos matemáticos relacionados con la teoría de bifurcaciones. De modo que a través de un estudio cualitativo del modelo podamos entender las distintas dinámicas neuronales que presentaremos en este trabajo. 1 Capítulo 2 Modelo matemático neuronal de Hodgkin y Huxley En este capítulo se desarrolla el modelo matemático neuronal de Hodgkin y Huxley. Considerado el esqueleto de la biomatemática moderna basada en modelos neuronales. Además vamos a explicar los mecanismos básicos que generan el impulso nervioso y su analogía con un circuito eléctrico. De modo que se introducen términos propios de la teoría de circuitos en el ejercicio de modelización, el cual nos conducirá a la obtención de un sistema de ecuaciones diferenciales no lineales capaz de recoger toda la dinámica neuronal. 2.1. Objetivo a modelar: Acción potencial Todas las células de los organismos presentan una diferencia de potencial, potencial de membrana, entre la parte interna y externa de la membrana, debida a la presencia de iones (partículas atómicas con carga), a ambos lados de ésta. En términos matemáticos, el potencial de membrana viene definido como: Vm=VIn −Vext Las membranas de las células presentan una serie de poros que permiten el flujo de iones a través de la membrana. Hay dos tipos de poros: con compuertas o sin éstas. Los primeros pueden estar abiertos o no, lo cual permite un flujo selectivo de iones. Mientras que los segundos siempre permanecen abiertos. Los principales iones que constituyen dicho flujo son: el sodio, Na+, el potasio, K+, y el cloro, Cl−. Donde el signo indica la carga que aportan, positiva o negativa. Otra de las características a destacar es el estado que presenta la célula según la dirección y los iones que constituyan la corriente eléctrica. Al inicio, la célula se halla en reposo. Durante este estado se dice que la célula presenta un potencial de reposo, que en el caso de las neuronas es aproximadamente de -70mV. Una corriente positiva hacia el interior, constituida por iones Na+, consigue que el potencial de membrana se vaya aproximando a cero. En este caso, se dice que la célula está despolarizada. Si, por el contrario, se sucede una corriente positiva de iones K+hacia el exterior y en sentido contrario una negativa de Cl−, se dice que la célula está hiperpolarizada. El potencial de membrana cambiará durante una acción potencial, por ejemplo. El potencial de acción es una onda de descarga eléctrica que viaja a lo largo de la membrana celular. En el caso de las neuronas, el potencial de acción es causado por un intercambio de iones a través de la membrana de la neurona como respuesta a un estímulo. Por tanto, el potencial de acción es la base de la transmisión del impulso nervioso. Si entendemos la dinámica de este proceso, entenderemos cómo se transmite la información en el organismo. 3 4Capítulo 2. Modelo matemático neuronal de Hodgkin y Huxley Hay diferentes formas de potencial de acción, pero todas tiene en común el efecto denominado todo o nada en la despolarización de la membrana. Es decir, si el voltaje no excede un valor crítico llamado umbral, no se iniciará la acción potencial y el potencial regresará a su nivel de reposo. Si, por el contrario, el voltaje supera dicho umbral, entonces la neurona se excitará dando lugar a un potencial de acción. 2.2. Modelo matemático neuronal de Hodgkin y Huxley Uno de los modelos más representativos y exitosos en neurociencia es el modelo neuronal diseñado por Alan Lloyd Hodkin y Andrew Huxley en 1952 [17], por el cual recibieron en 1963 el premio Nobel en Fisiología-Medicina. El objetivo de este modelo es describir la generación y propagación del potencial de acción. Hodgkin y Huxley formularon un modelo matemático para explicar el comportamiento de las células nerviosas del axón de un calamar gigante. Dicho axón tiene un diámetro de medio milímetro frente a un axón típico que es mil veces menor. En aquel momento, esta elección fue todo un éxito debido a la no existencia de microscopios eléctricos y de las simulaciones por ordenador. Además, se requirieron determinadas técnicas experimentales como: “Space Clamp”: es el mantenimiento de una distribución del voltaje uniforme en la membrana, sobre la zona donde se intenta medir la corriente, intensidad de la membrana. Para ello introdujeron un alambre delgado de plata a lo largo del axón, para simular lo ocurrido con un parche de membrana. “Voltage Clamp”: consiste en dejar que el potencial de la membrana se mantenga en un nivel de voltaje deseado. Gracias a las cuales consiguieron identificar de manera independiente la contribución individual de los iones. Demostrando que el sodio, Na+, y el potasio, K+, suponen una importante aportación a la corrientes iónicas que se generan durante una acción potencial. Además, Hodgkin y Huxley predijeron y, finalmente probaron, que la amplitud del potencial de acción depende de la concentración externa de sodio. Figura 2.1: Potencial de acción. Axón del tipo de calamar usado en los experimentos [19]. Modelos Matemáticos de Neuronas 2.3. Circuito eléctrico equivalente 5 2.3. Circuito eléctrico equivalente Las partículas que se mueven a través de la membrana determinan las propiedades eléctricas de las células y, en particular, de las neuronas. Por tanto, una manera bastante útil para describir el comportamiento del potencial de membrana es representar la célula en términos del circuito eléctrico equivalente, representado en la figura 2.2. El circuito consta de tres componentes: la membrana celular que puede ser vista como un condensador el cual tiene la capacidad de almacenar carga, las resistencias son representadas por los distintos tipos de canales iónicos incrustados en la membrana, mientras que las baterías quedan caracterizadas debido a los potenciales establecidos por las diferencias de concentración iónica entre la parte interna y externa. Figura 2.2: Circuito eléctrico equivalente propuesto por Hodgkin y Huxley para un segmento del axón del calamar [20]. En el circuito equivalente, mostrado en la figura 2.2, la corriente que atraviesa la membrana tiene dos componentes principales, una asociada con la capacidad eléctrica de la membrana y la otra asociada con el flujo de iones que circula por los canales. La intensidad de corriente, Ic, está definida como la tasa de cambio de la carga en la superficie de la membrana, es decir, Ic=dq/dt. La carga, q(t), está relacionada con el voltaje instantáneo de la membrana, Vm, y la capacidad de la membrana, Cm, por la relación q=CmVm. Por consiguiente, la intensidad corriente se reescribe como Ic=CmdVm/dt. La corriente iónica está subdividida en tres componentes, una corriente de sodio, INa, una corriente de potasio, IK, y una pequeña corriente, IL, debida a fugas constituida por iones de cloro, principalmente. Por consiguiente, usando las leyes de Kirchhoff1[18], el comportamiento del circuito eléctrico equivalente queda descrito por la siguiente ecuación diferencial: 1Las leyes de Kirchhoff son los axiomas sobre los que se asienta la teoría de circuitos. La primera ley de Kirchhoff se puede enunciar como: La suma algebraica de las intensidades de las corrientes que circulan por el conjunto de todos los elementos concurrentes en un punto, consideradas como entrantes en ese punto, es, en todo momento, cero. La segunda ley queda definida de la siguiente manera: La suma algebraica de las tensiones a lo largo de una línea cerrada, contabilizadas de acuerdo con un determinado sentido de circulación, es, en todo instante, cero. Autor: Nuria Begué Pedrosa 6Capítulo 2. Modelo matemático neuronal de Hodgkin y Huxley CmdVm dt +Iion =0,(2.1) Cuando una corriente externa es aplicada, entonces la suma anterior será igual a la intensidad aplicada, es decir: CmdVm dt =Iext −Iion,(2.2) donde Iext es la intensidad de corriente aplicada, que puede ser introducida a través de un electrodo intracelular. La ecuación presentada es fundamental para expresar el cambio en el potencial de membrana, Vm, debido al flujo de las corrientes que cruzan a través de la membrana. La intensidad de corriente iónica está relacionada con el voltaje de la membrana a través de la ley de Ohm2,V=IR. Sin embargo, en muchas ocasiones es más común presentarla en términos de la conductancia la cual se define: g=1/R. De modo que, obtenemos la relación I=gV. Otro concepto a tener en cuenta es el potencial de equilibrio de cada ion, Vi. El potencial de equilibrio se define como el potencial para el cual la corriente iónica neta que fluye a través de la membrana es cero. Este potencial recibe el nombre de potencial de Nernst 3[1]. La intensidad total, Iion, es la suma de cada una de las intensidades del circuito. Además, cada intensidad es el producto de la conductancia por la diferencia entre el potencial de membrana y el potencial de equilibrio, o potencial de Nernst. Se tiene que: Iion =INa +IK+IL=gNa(Vm−VNa) +gK(Vm−VK) +gL(Vm−VL).(2.3) Durante este minucioso trabajo de modelización, otro de los aspectos que fueron considerados son los canales que permiten el flujo iónico. Este movimiento iónico produce corrientes eléctricas que pueden llegar a producir el potencial de acción. El potencial de acción es debido a un aumento en la conductancia del sodio, gNa, provocando la entrada de iones Na+hacia el interior de la neurona haciendo, así, más positivo el interior (depolarización). Dicha conductancia varía en función del tiempo y empieza a disminuir aproximadamente hacia el máximo del potencial de acción, por lo que también depende del voltaje, es decir, gNa =F(t,V). Análogamente, el flujo de los iones de potasio conduce a la mismas conclusiones, por lo que la conductancia del potasio es una función del tiempo y del voltaje, es decir, gK=F(t,V). Luego, el problema a resolver era encontrar la función que modela las conductancias. En el modelo original propusieron que cada canal individual posee una compuerta la cual presenta dos estados permisivo y no permisivo. Por ello, definieron picomo la fracción de que las compuertas estuvieran abiertas. Aplicando la ley de acción de masas, se obtiene la transición entre los estados permisivo y no permisivo. La cual obedece la siguiente ecuación diferencial. dpi dt =αi(V)(1−pi)−βi(V)pi,(2.4) donde αi(V)yβi(V)son las constantes de proporción dependientes del voltaje, las cuales describen la transición entre los estados. 2La Ley de Ohm establece que "la intensidad de la corriente eléctrica que circula por un conductor eléctrico es directamente proporcional a la diferencia de potencial aplicada e inversamente proporcional a la resistencia del mismo". 3El potencial de Nernst viene dado por la siguiente fórmula: Vion =RT zF ln [ion]ext [ion]int . donde [ion]ext e[ion]int denotan las concentraciones externas e internas del ion respectivamente, Res la constante universal del gas, Tes la temperatura absoluta, Fes la constante de Faraday, y zes la carga del ion. Modelos Matemáticos de Neuronas 2.3. Circuito eléctrico equivalente 7 Como hemos comentado, la ecuación (2.4) corresponde a la ley de conservación de la masa recogida en el siguiente esquema: piαi ←−−1−pi,piβi −−→1−pi, donde αiyβiespecífica el número de transiciones a cada estado respectivamente. Si el voltaje de la membrana Vmse mantiene fijo para algún valor de V, entonces la fracción de compuertas abiertas, alcanzará un valor estacionario, es decir, dpi dt =0 cuando t−→∞, dado por: pi∞(V) = αi(V) αi(V)+βi(V),(2.5) El transcurso del tiempo para alcanzar este valor de equilibrio, está descrito por la constante del tiempo dada por: τi(V) = 1 αi(V)+βi(V).(2.6) Llevando las ecuaciones (2.5)y(2.6) a la ecuación (2.4) se cumple que: dpi dt =pi∞−pi τi.(2.7) Finalmente, cuando un canal se halla en estado permisivo contribuye en cierta medida a la conductancia total. La conductancia es proporcional al número de canales abiertos, los cuales son proporcionales, a su vez, a la probabilidad de que las compuertas asociadas se encuentren abiertas. Así, la conductancia total gjdebida a los canales de un tipo j, con compuertas de tipo i, es proporcional al producto de la posible aportación de cada compuerta individual, o pi: gj=gj r ∏ i=1 pi.(2.8) donde gjes una constante de normalización la cual determina la máxima conductancia posible cuando todos los canales están abiertos. 2.3.1. La conductancia del Potasio Para su caracterización realizaron experimentos basados en el “Voltage Clamp”. Y, así poder estudiar de manera independiente la conductancia del potasio sin la interacción del sodio. Observaron que la tendencia en los datos se debía a que un aumento en el voltaje aumentaba la corriente de iones. Una segunda tendencia es que la elevación de la conductancia es más rápida cuando se incrementa la despolarización. Hodgkin y Huxley supusieron la existencia de una partícula de activación del potasio n, que obedece la ecuación: dn dt =αn(V)(1−n)−βn(V)n,(2.9) En los experimentos se considera que el potencial de membrana comienza en el estado de reposo, Vm=0, el potencial es elevado a un nuevo potencial fijo Vc. Cuando la neurona se halla en reposo se cumple: n∞(0) = αn(0) αn(0)+βn(0),(2.10) Autor: Nuria Begué Pedrosa 8Capítulo 2. Modelo matemático neuronal de Hodgkin y Huxley Como Vmse mantiene fijo, sabemos por (2.5) que la variable nalcanzará un valor estacionario cuando t→∞dado por la siguiente ecuación: n∞(Vc) = αn(Vc) αn(Vc) +βn(Vc),(2.11) La solución de la ecuación diferencial dada por (2.9) describe la evolución en el tiempo de la conductancia del potasio: n(t) = n∞(Vc)−(n∞(Vc)−n∞(0))exp(−t τn(Vc)).(2.12) donde τn(Vc) = 1 αn(Vc) +βn(Vc).(2.13) Figura 2.3: Evolución de la variable de activación de los canales del potasio, n. (Hemos obtenido dicha gráfica al integrar en Octave el modelo completo dado por las ecuaciones (2.26),(2.27), (2.28), (2.29)). La ecuación (2.12) describe el transcurso de ncomo respuesta al cambio en el voltaje. Por otro lado, no tenemos que olvidar que nes una probabilidad, por tanto, toma valores entre 0 y 1, y debe ser multiplicada por una constante de normalización gK. Hodgkin y Huxley, gracias a los datos experimentales al fijar el voltaje en la ecuación (2.12), determinaron que el transcurso en el tiempo de la conductancia debía proporcionar una gráfica sigmoidal. Luego, determinaron que el valor de la curva de ajuste se obtenía para rn=4, en la ecuación (2.8). Entonces, la descripción de la conductancia del potasio viene dada por: gK=gKn4.(2.14) Modelos Matemáticos de Neuronas 2.3. Circuito eléctrico equivalente 9 2.3.2. La conductancia del Sodio Las observaciones experimentales llegaron a eludir las diferencias entre la conductancia del sodio respecto de la del potasio. El cambio se debe a que la conductancia del sodio es transitoria y solamente dura unos milisegundos. Por lo que la dinámica es más compleja. Sin embargo, la esencia del modelaje guarda una estrecha relación con el caso anterior. Para ajustar el comportamiento cinético, tuvieron que postular la existencia de dos partículas una de activación my otra de inactivación h. Dichas partículas son números adimensionales entre 0 y 1, es decir, reflejan la probabilidad de que las compuertas se hallen en dicho estado. dm dt =αm(V)(1−m)−βm(V)m,(2.15) dh dt =αh(V)(1−n)−βh(V)h.(2.16) Con una estrategia similar al caso anterior llegaron a una descripción para la conductancia del sodio dada por: gNa =gNam3h,(2.17) donde rm=3 y rh=1, son los exponentes que Hodgkin y Huxley determinaron para obtener la curva de ajuste que se aproximase a los datos experimentales sobre la conductancia. Figura 2.4: Dinámica de las variables de activación e inactivación de los canales del sodio solución del sistema dado por las ecuaciones (2.26),(2.27), (2.28), (2.29) que constituyen el modelo completo. Autor: Nuria Begué Pedrosa 16 Capítulo 3. Estudio cualitativo de modelos matemáticos de neuronas zdescribe la activación y/o inactivación de algunas corrientes. Dichas ecuaciones describen cualitativamente la dinámica presente en ciertos modelos neuronales basados en el formalismo del modelo de Hodgkin y Huxley. Las ecuaciones (3.3) y (3.4), con z=0, corresponden al modelo inicial de 1982 [22]. El cual fue mejorado en 1984 al incorporar la tercera ecuación, (3.5)[23]. No olvidemos que nuestra finalidad es el estudio cualitativo del modelo inicial, así como el modelo completo. Por lo tanto, en la siguiente sección vamos a desarrollar las teorías matemáticas necesarias para conseguirlo. 3.3. Conceptos de la Teoría de Bifurcaciones Como ya hemos comentado anteriormente, vamos a introducir los conceptos propios del estudio de bifurcaciones junto con los dos tipos de bifurcaciones que se van a analizar en este trabajo. Para más detalles ver [2]y[21]. Consideremos un sistema dinámico continuo en el tiempo dependiente de un vector de parámetros, α, ˙x=f(x,α),x∈Rn,α∈Rm,(3.6) donde fes suave con respecto a xyα. Al variar los valores del vector de parámetros puede suceder que el sistema siga siendo equivalente topológicamente al original, o no. Estas dos opciones quedan sintetizadas matemáticamente a través de las siguientes definiciones. Definición 3.3.1. Un sistema dinámico {T, Rn,ϕt}se dice topológicamente equivalente alrededor del equilibrio x0al sistema dinámico {T,Rn,ψt}alrededor del equilibrio y0si existe un homeomorfismo h:Rn−→Rnque cumple: (i) está definida en un pequeño entorno U ⊂Rnde x0; (ii) satisface y0=h(x0); (iii) las órbitas en el entorno U se corresponden a través del homeomorfismo con las de V =f(U)⊂ Rn, preservando la dirección del tiempo. Si Ues entorno abierto de x0, entonces Ves un entorno abierto de y0. Definición 3.3.2. La aparición de un retrato de fase topológicamente no equivalente bajo la variación de un parámetro recibe el nombre de bifurcación. Así, una bifurcación supone el cambio en el tipo de topología del sistema cuando el valor de sus parámetros atraviesan un valor crítico, o valor de bifurcación. Definición 3.3.3. Un diagrama de bifurcación de un sistema dinámico es una estratificación de su espacio de parámetros inducida por la equivalencia topológica, junto con los retratos de fase representativos de cada estrato. Las soluciones estables suelen representarse mediante líneas continuas, mientras que las soluciones inestables se representan con líneas punteadas. Una importante clasificación de las bifurcaciones viene dada por la siguiente definición, Modelos Matemáticos de Neuronas 3.3. Conceptos de la Teoría de Bifurcaciones 17 Definición 3.3.4. Una bifurcación local es una bifurcación que puede ser analizada puramente en términos de un cambio en la linealización alrededor del conjunto invariante o atractor. Una bifurcación global se produce cuando los cambios en las trayectorias en el espacio de fase no puede limitarse a un pequeño entorno. La aparición de bifurcaciones locales está asociado con la presencia de puntos críticos no hiperbólicos. Dado el sistema dinámico (3.6), existen genéricamente solamente dos maneras de que el punto de equilibrio deje de ser hiperbólico. Se puede dar el caso de que el valor propio se aproxime a cero y tengamos λ1=0 (figura 3.2), o un par de valores propios complejos que cruzan el eje imaginario y se tenga λ1,2=±iω0,ω0>0 (figuras 3.4 y3.3) para algún valor del parámetro. Dichas bifurcaciones serán presentadas matemáticamente de manera formal . Otro de los conceptos que aparece en la teoría de bifurcaciones es la codimensión, la cual se define intuitivamente: Definición 3.3.5. La codimensión de un punto de bifurcación es la diferencia entre la dimensión del espacio de parámetros y la correspondiente conjunto de bifurcaciones. Alternativamente, la codimensión es el número de parámetros necesarios para añadir de manera genérica la bifurcación al sistema. En nuestro estudio nos encontramos con los dos siguientes tipos de bifurcaciones locales de codimensión-1: la bifurcación fold y la bifurcación Andronov-Hopf. 3.3.1. Bifurcación fold o silla-nodo La bifurcación fold es el mecanismo básico de creación-destrucción de puntos fijos. Lo que sucede es que conforme un parámetro del sistema crece o decrece, los puntos de equilibrio del mismo se aproximan unos a otros, colisionan y desaparecen. Formalmente, se recoge la siguiente definición: Definición 3.3.6. La bifurcación asociada con la aparición de λ1=0se llama, genéricamente (en el caso de bifurcaciones de puntos críticos), bifurcación fold, o silla-nodo. Para entender mejor esta bifurcación vamos a considerar el siguiente sistema dinámico uni-dimensional dependiente de un parámetro de control, α: ˙x=α+x2≡f(x,α). Figura 3.2: Bifurcación de “fold” o silla-nodo [2]. Autor: Nuria Begué Pedrosa 18 Capítulo 3. Estudio cualitativo de modelos matemáticos de neuronas En la figura 3.2, se observa como al variar αel retrato de fases cambia drásticamente ya que pasamos de tener dos puntos de equilibrio a presentar uno, y finalmente ninguno. Por tanto, el cambio en el número de puntos de equilibrio nos lleva a concluir que estamos ante una bifurcación. En este caso el valor de bifurcación o valor crítico se localiza para α=0. Teorema 3.3.7. (Forma normal topológica para la bifurcación de fold). Sea un sistema unidimesional, dx dt =f(x,α),x∈R1,α∈R1, con f “suave”, tal que para α=0posee un equilibrio en x =0, y λ=fx(0,0) = 0. Si el sistema satisface las siguientes condiciones de no-degeneración: (H.1) a(0) = fxx(0,0)6=0, (H.2) fα(0,0)6=0, entonces es topológicamente equivalente a ˙ η=β±η2. Para el caso multi-dimensional se recoge el siguiente teorema [26], Teorema 3.3.8. Sea el siguiente sistema dinámico, dx dt =f(x,α),x∈Rn,α∈R1, donde n ≥2, tal que para α=0posee un equilibrio en el origen tal que la matriz Jacobiana, A, evaluada en el origen tiene un único valor propio, λ1=0, nsvalores propios con Re(λj)<0y ns valores propios con Re(λj)>0, con ns+nu+1=n. Entonces, de acuerdo con el teorema de la variedad central, existe una familia invariante unidimensional de variedades Wc αcerca del origen. El sistema restringido sobre Wc αen uni-dimensional. Por lo tanto, es equivalente a la forma normal del caso uni-dimensional. Además, bajo las condiciones (H.1) y (H.2) (del teorema 3.3.7), el sistema n-dimensional es topológicamente equivalente cerca del origen a ˙y=β±y2, ˙ys=−ys, ˙yu= +yu. donde y ∈R1, ys∈Rns, yu∈Rnu. Definición 3.3.9. El coeficiente cuadrático, a(0), se define para n ≥1. Se considera el desarrollo de Taylor de f(x,0)en x =0como f(x,0) = A0x+1 2B(x,x)+O(kxk3), donde B(x,y) es la función multilineal con componentes: Bj(x,y) = n ∑ k,l=1 ∂2f(ξ,0) ∂ξk∂ξlξ=0 xkyl, Modelos Matemáticos de Neuronas 3.3. Conceptos de la Teoría de Bifurcaciones 19 donde j = 1,2,..., n. Sea q ∈Rnun vector tal que A0q=0. Sea p ∈Rn,el vector adjunto tal que AT 0p=0, donde p,q=1. Entonces, a(0) = 1 2p,B(q,q). 3.3.2. Bifurcación Andronov-Hopf La bifurcación de Andronov-Hopf está asociada con la aparición de órbitas o ciclos límite estables o inestables, gracias a la cual podemos garantizar la existencia de ciclos límites en sistemas de dimensión mayor que el caso planar, el cual queda recogido en el Teorema de Poincaré-Bendixon. Por tanto, la importancia del estudio de esta bifurcación radica, también, por constituir una herramienta para la justificación de la aparición de ciclos límite en sistemas de dimensión mayor. Definición 3.3.10. La bifurcación correspondiente con la presencia de λ1,2=±iω,ω>0,recibe el nombre de bifurcación de Hopf, o Andronov-Hopf. En este tipo de bifurcación se produce una pérdida de estabilidad debido a que los autovalores complejos pasan de tener parte real negativa, a tener parte real positiva, atravesando por tanto el eje imaginario. Se distinguen al igual que antes dos tipos: bifurcación de Hopf supercrítica y subcrítica. Vamos a definir el siguiente concepto conocido como coeficiente de Lyapunov el cual juega un papel trascendente para este tipo de bifurcación. Definición 3.3.11. (El primer coeficiente de Lyapunov). Se considera el desarrollo de Taylor de f(x,0) en x = 0 como f(x,0) = A0x+1 2B(x,x)+ 1 6C(x,x,x)+O(kxk4), donde B(x,y) y C(x,y,z) son funciones multilineales con componentes: Bj(x,y) = n ∑ k,l=1 ∂2f(ξ,0) ∂ξk∂ξlξ=0 xkyl, Cj(x,y,z) = n ∑ k,l,m=1 ∂3f(ξ,0) ∂ξk∂ξl∂ξmξ=0 xkylzm, donde j = 1,2,..., n. Sea q ∈Cnun autovector complejo de A0correspondiente al autovalor iωital que A0q=iω0q. Sea p ∈Cn,el autovector adjunto tal que AT 0p=−iω0p, p,q= 1. Entonces, el primer coeficiente de Lyapunov se define l1(0) = 1 2ω0Re[p,C(q,q,q)−2p,B(q,A−1 0B(q,q))+p,B(q,(2iω0In−A0)−1B(q,q))]. Notar que el valor, pero no el signo, de l1(0)depende de la norma del vector q. La normalización q,q=1 es una de las opciones para eliminar la ambigüedad. Autor: Nuria Begué Pedrosa 20 Capítulo 3. Estudio cualitativo de modelos matemáticos de neuronas Figura 3.3: Bifurcación de Hopf supercrítica (ciclos límite estables), σ<0. Podemos observar como inicialmente el sistema presenta un punto crítico estable que tras la bifurcación pierde la estabilidad surgiendo, además, un ciclo límite estable [2]. Retomando el estudio de la bifurcación de Andronov-Hopf se presenta el siguiente teorema para el caso bidimesional el cual caracteriza dicha bifurcación, así como su tipo. Teorema 3.3.12. (Forma normal topológica para la bifurcación de Hopf). Sea un sistema bidimensional, dx dt =f(x,α),x∈R2,α∈R, con f “suave”, tal que para |α| pequeño posee un equilibrio en x=0 con valores propios de la matriz Jacobiana de f en x =0, λ1,2(α) = µ(α)±iω(α), donde µ(0) = 0,ω(0) = ω0>0. Si el sistema satisface las condiciones de no degeneración: (H.1) µ0(0)6=0, (H.2) l1(0)6=0,donde l1es el primer coeficiente de Lyapunov. entonces es topológicamente equivalente a ˙y1=αy1−y2+σy1(y2 1+y2 2), ˙y2=y1+αy2+σy2(y2 1+y2 2), siendo σ=sign(l1(0)). Nota.-Se observa que el signo del primer coeficiente de Lyapunov nos va a determinar qué tipo de bifurcación de Hopf sufre el sistema. Las figuras 3.3 y3.4 nos muestran, respectivamente, los retratos de fases que surgen bajo la variación del parámetro de control, α, de los sistemas de ecuaciones diferenciales que constituyen la forma normal de la bifurcación de Hopf (resultado del teorema, 3.3.12). En ambos casos, el valor crítico se tiene para α=0. Modelos Matemáticos de Neuronas 3.3. Conceptos de la Teoría de Bifurcaciones 21 Figura 3.4: Bifurcación de Hopf subcrítica (ciclos límite inestables), σ>0. Podemos notar que inicialmente el sistema presenta un ciclo límite inestable asociado a un punto crítico estable.Y,después de la bifurcación el punto crítico pierde su estabilidad convirtiéndose en un foco inestable [2]. Definición 3.3.13. Para el caso bidimensional se obtiene una expresión del primer coeficiente de Lyapunov dada por: l1(0) = 1 2ω2 0 Re(ig20g11 +ω0g21), donde g20 =p,B(q,q),g11 =p,B(q,q),g02 =p,B(q,q). Una vez que nos hemos introducido en el mundo de las bifurcaciones, nuestro propósito reside en aplicar los conocimientos adquiridos para entender los modelos objeto de nuestro estudio. Nos vamos a centrar en localizar el valor crítico para el cual se da la bifurcación, así como su tipo y, finalmente, cómo afecta a la dinámica neuronal correspondiente al impulso nervioso. Autor: Nuria Begué Pedrosa 22 Capítulo 3. Estudio cualitativo de modelos matemáticos de neuronas 3.4. Modelo Matemático Neuronal de Hindmarsh y Rose (1982) En 1982, Hindmarsh y Rose proponen un modelo constituido por las ecuaciones (3.3) y (3.4), donde z=0. Este modelo presenta dos puntos de equilibrio adicionales al modelo de Fitzhugh-Nagumo. El primer punto para el estado de reposo o potencial de reposo y el segundo asociado al ciclo límite estable. Vamos a realizar un estudio de la estabilidad local de los puntos críticos del modelo en función del conjunto de parámetros. En primer lugar, definimos las ecuaciones que satisfacen las isoclinas, ˙x=0, ˙y=0, cuyas soluciones constituyen los puntos de fijos del sistema. Por tanto, los puntos críticos vienen dados por las raíces de la siguiente ecuación; x3+tx2=s(3.7) donde t=1 a(d β−b),s=1 a(I+c β), Para nuestro estudio, suponemos que I=0, es decir, que la intensidad aplicada es nula. Otra de las hipótesis que establecemos es la no negatividad de los parámetros que constituyen este modelo. Sabemos que el modelo presenta tres puntos de equilibrio. Por tanto, la condición para que existan se traduce en que el discriminante sea positivo. En este caso se tiene la siguiente desigualdad: 27ca2 β<4d β−b3 ,(3.8) Debido a la condición de no negatividad de los parámetros, se sigue inmediatamente que: d β>b.(3.9) Sea (xeq,yeq), punto de equilibrio. Consideramos la aproximación lineal del sistema que nos permite realizar un estudio local de la estabilidad, siendo la matriz Jacobiana, A(xeq,yeq), la que gobierna el sistema dinámico linealizado. El tipo del punto de equilibrio xeq queda determinado por los signos del determinante y la traza de A(xeq,yeq).El determinante y la traza vienen dados por: det(A(xeq,yeq)) = 3aβx2 eq +2(d−bβ)xeq. traza(A(xeq,yeq)) = −3ax2 eq +2bxeq −β. En la siguiente tabla se recoge la estabilidad local del punto de equilibrio según a la región que pertenezca la coordenada xdel punto. Modelos Matemáticos de Neuronas 3.4. Modelo Matemático Neuronal de Hindmarsh y Rose (1982) 23 xeq signo(det(A)) signo(tr(A)) tipo de equilibrio x<−2(d β−b) 3a+ - nodo estable −2(d β−b) 3a<x< 0 - - punto de silla 0 < x<b−pb2−3aβ 3a+ - foco estable b−pb2−3aβ 3a<x<b+pb2−3aβ 3a+ + foco inestable x>b+pb2−3aβ 3a+ - foco estable Concluimos que el potencial de reposo queda representado por un nodo estable que pertenece a la primera fila. En la segunda fila, el punto de silla corresponde al umbral del potencial para el cual se da la acción potencial. Finalmente, el periodo de estallidos queda reflejado por el ciclo límite estable asociado al foco inestable (fila quinta). Para un valor genérico de I, consideramos el sistema bidimensional dado por las ecuaciones (3.3) y (3.4), con z=0. El cálculo analítico de los puntos de equilibrio del modelo con SAGE nos devuelve expresiones complejas, por lo que obtener conclusiones resulta una tarea ardua. Por ello, vamos a escoger un conjunto de valores estándar para los parámetros (a=1, b=3, c=1, d=5 y β=1)[16]. Entonces el modelo queda definido por el siguiente conjunto de ecuaciones diferenciales, ˙x=−x3+3x2+y+I,(3.10) ˙y=1−5x2−y.(3.11) donde Iconstituye nuestro parámetro de control. En primer lugar, si consideramos I=0, el conjunto de parámetros cumple la desigualdad (3.8). Luego, obtenemos tres puntos fijos: (−1 2(√5+1),−1 2(5√5+13)), tal que −1 2(√5+1)∈(−∞,−4 3). Entonces, según la tabla es un nodo estable. (−1,−4), tal que −1∈(−4 3,0). Entonces, el punto crítico corresponde a un punto silla. (1 2(√5−1),1 2(5√5−13)), tal que 1 2(√5−1)∈(1 3(3−√6),1 3(3+√6)). Entonces, el punto es un foco inestable asociado a un ciclo límite estable. Para comprobarlo, hemos estudiado la estabilidad local de cada punto a través de la matriz Jacobiana, asociada al sistema (3.10), (3.11), y las conclusiones son idénticas. Autor: Nuria Begué Pedrosa 24 Capítulo 3. Estudio cualitativo de modelos matemáticos de neuronas Figura 3.5: Soluciones al modelo descrito por las ecuaciones (3.10), (3.11); con I=0. Las gráficas de la izquierda muestran las isoclinas junto con el ciclo límite estable, mientras que las de la derecha dibujan el potencial de acción como función del tiempo. Otro aspecto a señalar es que, para este valor del parámetro de control, estamos ante un caso de biestabilidad. Es decir, dependiendo de las condiciones iniciales el sistema se mueve hacia uno de los dos atractores (nodo estable y ciclo límite estable). Los sistemas que presentan multi-estabilidad son muy frecuentes en modelos biológicos. La figura 3.5 demuestra que el modelo para I=0 se comporta de una manera u otra dependiendo de las condiciones iniciales. En el primer caso se sobrepasa el umbral, lo que da lugar a la transmisión del impulso nervioso. En términos cualitativos se traduce en que el sistema entra en el ciclo límite estable produciéndose un conjunto de espigas todas ellas regulares (“spiking” regular). En el segundo, el sistema evoluciona hacia el nodo estable, el cual representa el estado de reposo de la neurona. En segundo lugar, este modelo nos va a facilitar nuestro objetivo, la comprensión cualitativa del potencial de acción. El parámetro de control escogido es la intensidad aplicada a la neurona. Por tanto, vamos a estudiar para distintos valores de Icuál es el comportamiento dinámico de la neurona cuando tiene lugar el impulso nervioso. La figura 3.6 nos ayuda a entender el sistema cualitativamente. En esta gráfica (y todas las posteriores), el diagrama de bifurcaciones ha sido calculado con MAPLE, así como la localización de los Modelos Matemáticos de Neuronas 3.4. Modelo Matemático Neuronal de Hindmarsh y Rose (1982) 25 Figura 3.6: Diagrama de bifurcaciones para el modelo neuronal de Hindmarsh-Rose donde Int (I) es el parámetro de control. La curva está constituida por los puntos críticos del modelo dado por las ecuaciones (3.10)y(3.11) al variar el parámetro de control, I. Además, la figura señala las bifurcaciones que tienen lugar en el modelo. Las etiquetas AH y SN hacen referencia a Andronov-Hopf y saddle-node (bifurcación de “fold”), respectivamente. puntos de bifurcación. La curva está constituida por los puntos de equilibrio del sistema al variar el parámetro de control, además de señalar los valores de bifurcación junto con el tipo de bifurcación sucedida. En primer lugar, los dos puntos de tangencia representan el momento en el que tiene lugar la bifurcación de tipo nodo-silla o “fold” (SN), entre las dos bifurcaciones se localizan un conjunto de puntos silla. Por otro lado, la rama inferior está formada por nodos estables. Finalmente, la rama superior presenta un comportamiento más complejo sucediéndose la bifurcación de Andronov-Hopf (AH). Dicha bifurcación se localiza para I≃−0,9265. Para demostrarlo, consideramos el sistema de ecuaciones diferenciales dadas por (3.10), (3.11). Se cumple que para el valor de I≃ −0,9265, la matriz Jacobiana evaluada en x≃0,183503, y≃0,92647378 tiene como valores propios, λ1,2≃±i0,9143302193 (teniendo valores de la parte real positivos antes y negativos después de la bifurcación), es decir, la parte real es casi nula. Además, el sistema satisface las condiciones de no degeneración, l1≃−5,340144708, de donde σ=sign(l1)<0. Entonces, por el teorema 3.3.12, el sistema sufre una bifurcación de Andronov-Hopf supercrítica. Por lo tanto, tenemos demostrada la siguiente proposición: Proposición 3.4.1. Sea el sistema dinámico dado por las ecuaciones (3.10), (3.11), con z =0. Entonces, para el valor de I ≃−0,9265 el sistema sufre una bifurcación de Andronov-Hopf supercrítica cuyo retrato de fase se corresponde al de la figura 3.3. Por tanto, tras la bifurcación, el foco estable pierde su estabilidad, apereciendo junto a él un ciclo límite estable. Si estudiamos el modelo anterior para valores de I≥0, podemos afirmar que al aplicar una corriente externa con una amplitud suficiente (despolarización), se alcanza un punto de bifurcación para el valor del parámetro I=5/27, que corresponde con una bifurcación “fold” o silla-nodo. Para localizar el punto de bifurcación, consideramos las isoclinas del sistema dado por las ecuaciones (3.10), (3.11). Así, se obtiene: 0=f(x,I) = −x3−2x2+1+I Autor: Nuria Begué Pedrosa 32 Capítulo 3. Estudio cualitativo de modelos matemáticos de neuronas Figura 3.14: Retrato de fases para el valor de z=1,0856. El sistema presenta tres puntos de equilibrios, foco inestable (verde), un punto silla (negro), y nodo estable (azul). Las variedades estable e inestable del punto de silla dividen el retrato de fase en dos regiones atractoras, producidas debido al nodo estable y al ciclo límite estable asociado al foco inestable. A medida que el valor de zaumenta, el ciclo límite se aproxima al punto de silla y, cuando toquen se sucederá la bifurcación. En general, la bifurcación homoclínica es un mecanismo de destrucción o construcción de ciclos límite. Así, cabe señalar que el signo de σdetermina la estabilidad del ciclo límite: nace estable si σ<0 e inestable en el otro caso. Otro de los rasgos que caracteriza a este tipo de bifurcaciones se centra en la ley de escalado asociada a bifurcaciones de ciclos límite de modelos bidimensionales, [17]. En el caso de la bifurcación homoclínica, el periodo del ciclo crece en escala logarítmica. Este hecho explica que aumente el intervalo en que se suceden los estallidos. Con objeto de que entendamos esta bifurcación incluimos la figura 3.15, que muestra ejemplos de órbitas homoclínicas en el espacio. Definición 3.5.1. Una órbita homoclínica es una trayectoria, ϕ, que conecta la variedad estable e inestable de un punto silla. En nuestro caso, la trayectoria que conecta las variedades viene dada por el ciclo límite. Tras la bifurcación, el ciclo límite estable desaparece de modo que el sistema evoluciona hacia el nodo estable, dando lugar a la fase de recuperación del “bursting”. Al moverse por la rama inferior el sistema sufre una bifurcación tipo “fold” para el valor z=83/54. Para ese valor del parámetro, el sistema presenta el punto de bifurcación y el ciclo límite estable asociado al foco inestable. De este modo, el sistema evoluciona hacia la variedad rápida, teniendo lugar el conjunto de estallidos. Modelos Matemáticos de Neuronas 3.6. Vuelta a la realidad 33 Figura 3.15: Ejemplos de órbitas homoclínicas [2]. La figura de la izquierda se corresponde con una órbita homoclínica para un punto silla. La figura de la derecha muestra una órbita homoclínica para un punto silla-foco. En ambos casos, la órbita homoclínica Γ0conecta la variedad estable, Ws, con la variedad inestable, Wu, del punto silla. Entonces, los mecanismos que inician y finalizan la variedad rápida del “bursting” quedan descritos por la aparición de dos bifurcaciones, la bifurcación tipo “fold” y la bifurcación homoclínica, respectivamente. Observemos como este tipo de “bursting” está señalado en la tabla de la figura 3.7. Si contrastamos las figuras 3.8 y3.9, la diferente forma del estallido se debe precisamente a la existencia de la bifurcación homoclínica, que tiene como consecuencia la desaparición del ciclo límite de forma repentina (debido a la conexión con el punto de equilibrio), causando que el sistema evolucione hacia el estado de reposo. Si trazásemos una recta vertical atravesando la bifurcación, percibiríamos el corte vertical del estallido, perdiendo la forma de paraboloide del primero. 3.6. Vuelta a la realidad En este capítulo hemos elaborado el estudio cualitativo de distintos modelos que tienen por objeto recoger las distintas dinámicas en las que se presenta la acción potencial. Al tratar cada modelo como un sistema dinámico, podemos introducir toda una teoría matemáticamente fundamentada que nos permite entender cada actividad en función de las bifurcaciones sucedidas. Sin embargo, nos suscita interés saber si toda esta abstracción nos permite obtener conclusiones reales. La respuesta afirmativa nos aporta una grata satisfacción, ya que a través de teorías fundamentadas conseguimos reproducir con gran exactitud la dinámica neuronal que tiene lugar durante una acción potencial. Notemos que si recordamos la figura 3.1, podemos observar que está compuesta de una serie de gráficas donde se aprecia que los potenciales de acción se asemejan al “Square-Wave bursting”. Por ejemplo, en la fila alas gráficas etiquetadas con el nombre de FRB, IB. En la fila blas 2, 3, 4 son ejemplos de este tipo de “bursting”. Las neuronas más comunes de la corteza presentan una dinámica basada en disparos regulares o “spiking” regular. Nos gustaría comentar que nuestro organismo está constituido por distintos tipos de neuronas dependiendo del lugar que ocupan y la función que realizan. Esta diferenciación está también presente en el potencial de acción que reproducen a la hora de transmitir el impulso nervioso, por ello las gráficas que recogen la acción potencial gozan de variedad y riqueza. Autor: Nuria Begué Pedrosa Capítulo 4 APÉNDICE: Técnicas y herramientas utilizadas. 1) Métodos de resolución de EDO’s. En este trabajo la formulación del problema viene dado como un sistema de ecuaciones diferenciales no lineales ordinarias. Por lo tanto, ha sido necesario el uso de integradores numéricos para resolver dichas ecuaciones. Hemos usado los integradores ode23 yode45 de tipo RungeKutta que ofrece Octave y el programa Dynamics Solver con el método Dormand −Prince8(5,3), destinado para la resolución de ecuaciones diferenciales ordinarias. Este método pertenece también a la familia de métodos Runge-Kutta. 2) Métodos de análisis cualitativo. Todas las gráficas que recogen los diagramas de bifurcaciones de los distintos modelos, así como el estudio del retrato de fases para valores críticos del parámetro de control escogido las hemos obtenido gracias al programa Maple, el cual está orientado al cálculo de problemas matemáticos. Por otro lado, la localización de los puntos de equilibrio del sistema y su estabilidad la hemos estudiado a través de Maple ySAGE. Este último constituye un software libre. 35 Bibliografía [1] G.B. Ermentrout, D.H. Terman, Mathematical Foundations of Neuroscience, Interdisciplinary Applied Mathematics, Springer, 35, (2010). [2] Y. Kuznetsov, Elements of applied bifurcation theory, Springer, (2004). [3] A. Shilnikov, M.Kolomiets, Methods of qualitive theory for Hindmarsh-Rose model: A case study. A tutorial, International Journal of Bifurcation and Chaos, Vol. 18, No. 8 2141-2168 (2008). [4] R. Barrio, A. Shilnikov, Parameter-sweeping techniques for temporal dynamics of neuronal systems: case study of Hindmarsh-Rose model, Journal of Mathematical Neuroscience, 10.1186/2190-8567-1-6, (2011). [5] G. Innocenti, A. Morelli, R. Genesio, A. Torcini, Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos, Chaos 17, 043128, (2007). [6] M. Izhikevich, Bursting, Scholarpedia, (2006). http://www.scholarpedia.org/article/Bursting [7] R. Zillmer, Simple Neuron Models: FitzHugh-Nagumo and Hindmarsh-Rose, INFN, Sezione di Firenze. [8] J. Macuk, Introducción a las Ecuaciones Diferenciales en Maple, Facultad de Ingeniería, Universidad Diego Portales, (2003). [9] J. M. Oviedo, E. López, D. Busignani, Manual Introductorio a Maple 8.0: Álgebra Lineal Avanzada y Optimización, Departamento de Estadística y Matemática, Facultad de Ciencias Económicas, Universidad Nacional de Córdoba, (2005). [10] N. Carrillo, Modelización de la actividad neuroeléctrica, tesis Docotral, Universidad nacional autónoma de México, (2010). [11] R. Huerta, Un método para predecir el número de potenciales de acción producidos por el modelo clásico de Hodgkin Huxley cuando la corriente aplicada es constante, tesis Doctoral, (2004). [12] F. Verhulst, Periodic solutions and slow manifolds,International Journal of Bifurcation and Chaos, Vol. 17, No. 8, 2533-2540, (2007). [13] M. Hana, H. Zhua, The Loop Quantities and Bifurcations of Homoclinic Loops, Journal of differential equations, Vol. 234, 339-359, (2007). [14] S. Kahan, Bifurcaciones homoclínicas en el circuito de Chua, tesis Doctoral, Universidad de la República-Uruguay, (1997). [15] J.M. Valiente, Manual de iniciación a GNU Octave, Trabajo realizado dentro de un Proyecto Fin de Carrera en la E.U. Politécnica de Teruel, (2006). 37 38 BIBLIOGRAFÍA [16] E.Steur, Parameter Estimation in Hindmarsh-Rose Neurons, Traineeship report, Universidad técnica de Eindhoven, (2006). [17] A.L. Hodgkin, A.F. Huxley A Cuantitative Description of Membrane Current and its Application to Conduction and Excitation in Nerve. J Physil. 117, 500-544 (1952). [18] A. Pastor, J. Ortega, V.M. Parra, Á. Pérez, Circuitos eléctricos. Volumen I, Universidad Nacional de Educación a Distancia, 978-84-362-4981-1, (2007). [19] P. Lamberti, V. Rodríguez, Desarrollo del modelo matemático de Hodgkin y Huxley en neurociencias,(2007). http://electroneubio.secyt.gov.ar/Lamberti-Rodriguez_Hodgkin-Huxley.htm. [20] K.L. Anderson, J. Chism, Q. Hale, P. Klockenkemper, C. Pinkett, C. Smith, D. Badamdorj, Mathematical Modeling Action Potential in Cell Processes, (2013). http://www.tnstate.edu/mathematics/mathreu/filesreu/Action_Potential_Report.pdf [21] Steven H. Strogatz, Nonlinear Dynamics and Chaos, Westview Press, (1994). [22] J.L. Hindmarsh, R.M. Rose, A model of the nerve impulse using two first-order differential equations, Nature, 269:162-164, (1982). [23] J.L. Hindmarsh, R.M. Rose, A model of the nerve impulse using three coupled first-order differential equations, Proc. R. Soc. Lond. B, Biol. (1984). [24] M. Steriade, Neocortical cell classes are flexible entities, Nature Reviews Neuroscience 5, 121134 (2004). [25] P. A. Tipler, G. Mosca, Física para la ciencia y la tecnología. Volumen 2. Electricidad y magnetismo. Luz., Reverte, (2010). [26] Y. Kuznetsov, Saddle-node bifurcation, Scholarpedia, (2006). Modelos Matemáticos de Neuronas