Análisis numérico del modelo de Hodgkin-Huxley

Traballo Fin de Grao Análisis numérico del modelo de Hodgkin-Huxley Alba Rois González 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Análisis numérico del modelo de Hodgkin-Huxley Autora: Alba Rois González Tutores: Óscar López Pouso y Carmen Rodríguez Iglesias Julio, 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Conocimiento: Matemática Aplicada Título: Análisis numérico del modelo de Hodgkin-Huxley Breve descripción del contenido Este modelo, propuesto por Hodgkin y Huxley en 1952, consiste en un sistema de cuatro ecuaciones diferenciales no lineales que describen la iniciación y transmisión del potencial de acción en células excitables, como las neuronas o los miocitos cardíacos. El objetivo de este trabajo es resolver las ecuaciones que componen el modelo mediante distintos métodos numéricos e implementarlos en Matlab R . Se prestará también atención a las propiedades numéricas de los distintos esquemas. Recomendaciones Otras observaciones iii Índice general Resumen viii Introducción xi 1. Descripción del modelo de Hodgkin-Huxley 1 1.1. Fisiología de la neurona. Potencial de acción . . . . . . . . . . . . . . . . . . 1 1.2. Modelomatemático................................ 8 1.2.1. Circuito eléctrico equivalente . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2. Modelocompleto............................. 13 1.2.3. Consideraciones acerca del modelo . . . . . . . . . . . . . . . . . . . 14 2. Resolución del modelo bajo hipótesis de homogeneidad en espacio 17 2.1. Existencia y unicidad de solución del PVI . . . . . . . . . . . . . . . . . . . 18 2.2. Métodos numéricos de resolución de EDO . . . . . . . . . . . . . . . . . . . 22 2.2.1. Eulerexplícito .............................. 22 2.2.2. Runge-Kutta ............................... 23 2.3. Análisis comparativo de los métodos numéricos . . . . . . . . . . . . . . . . 25 2.4. Análisis de los resultados de la simulación numérica . . . . . . . . . . . . . . 30 3. Resolución del modelo sin hipótesis de homogeneidad en espacio 37 3.1. Métodos numéricos de resolución de EDP . . . . . . . . . . . . . . . . . . . 40 3.1.1. Implícito estándar . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1.2. Crank-Nicolson.............................. 46 3.2. Consideraciones acerca de los métodos numéricos . . . . . . . . . . . . . . . 48 3.3. Análisis de los resultados de la simulación numérica . . . . . . . . . . . . . . 49 4. Conclusiones y nuevas líneas de trabajo 53 Índice de parámetros del modelo 57 v vi ÍNDICE GENERAL Apéndice I: Códigos de Matlab R 59 I.I. Modelo bajo hipótesis de homogeneidad en espacio . . . . . . . . . . . . . . . 59 I.I.I.Entradadedatos .............................. 59 I.I.II.Eulerexplícito ............................... 60 I.I.III. Función intrínseca ode23 . . . . . . . . . . . . . . . . . . . . . . . . . 62 I.I.IV. Runge-Kutta clásico de 4 etapas . . . . . . . . . . . . . . . . . . . . . 64 I.I.V.FunciónF.................................. 66 I.I.VI. Función corriente externa . . . . . . . . . . . . . . . . . . . . . . . . 67 I.I.VII. Funciones Alpha y Beta . . . . . . . . . . . . . . . . . . . . . . . . . 67 I.II. Modelos sin hipótesis de homogeneidad en espacio . . . . . . . . . . . . . . . 69 I.II.I.Entradadedatos.............................. 69 I.II.II. Implícito estándar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 I.II.III.Crank-Nicolson.............................. 72 Apéndice II: Sobre los paramétros de activación e inactivación de Na y K 77 Apéndice III: Efecto de la temperatura 81 Bibliografía 83 2CAPÍTULO 1. DESCRIPCIÓN DEL MODELO DE HODGKIN-HUXLEY Figura 1.1: Estructura de una neurona. Recuperado de https://medium.com/soldai/ inspiracion-biologica-de-las-redes-neuronales-artificiales-9af7d7b906a Antes de presentar el modelo propuesto por Hodgkin y Huxley, es necesario, no obstante, entender qué iones son los que se mueven a través de la membrana celular y cómo es este movimiento. En una célula, las especies iónicas más frecuentes son el ion sodio (Na`), el ion potasio (K`) y, en menor medida, aniones orgánicos. El movimiento de estos iones a través de la membrana celular es la base del inicio y propagación del potencial de acción, pero este movimiento no es arbitrario, sino que obedece a las características estructurales de la membrana celular. Generalidades acerca de la membrana celular La membrana celular es una bicapa lipídica que aísla el medio intracelular del medio extracelular, y una de sus funciones es regular el paso de sustancias que entran y salen de la célula. Se compone fundamentalmente de fosfolípidos, que son altamente hidrofóbicos, junto con otros tipos de lípidos, proteínas y carbohidratos. Su composición determina que actúe como una barrera selectiva que permite el paso de algunas moléculas pequeñas y no cargadas eléctricamente, pero impide el paso de moléculas hidrosolubles y cargadas eléctricamente (iones). Para permitir el paso de estas últimas, existen unas estructuras proteicas transmembrana, esto es, que atraviesan la membrana 1.1. FISIOLOGÍA DE LA NEURONA. POTENCIAL DE ACCIÓN 3 celular, y que son los canales ybombas iónicas. Ambos son selectivos para uno o varios tipos de iones, de manera que, a través de mecanismos físico-químicos, son capaces de discriminar las distintas especies iónicas y dejar pasar solo algunas de ellas [19]. Se distinguen dos tipos de canales iónicos: los canales de fuga, que están siempre abiertos, y otros que se abren o se cierran en respuesta a estímulos químicos (canales regulados por ligandos) o a cambios en el voltaje transmembrana (canales regulados por voltaje). Tal y como se comentó anteriormente, los canales son selectivos para ciertos tipos de iones, de tal forma que se hablará de canales de una cierta especie iónica, tales como canales de sodio o potasio, para referirse a canales altamente selectivos (permeables) para esa especie. Los canales, a diferencia de las bombas, no utilizan energía para el movimiento de los iones, sino que, cuando están abiertos, los iones fluyen a través de ellos en el sentido que determina su gradiente electroquímico. El gradiente electroquímico de un cierto ion es el resultado de la suma de dos tipos de fuerzas: la fuerza electrostática, que tiende a impulsar al ion hacia el medio donde predominan las partículas de carga contraria, y el gradiente de difusión o concentración, que determina el movimiento de un ion desde un lugar donde hay más concentración a otro donde hay menos. Por otra parte, las bombas iónicas utilizan energía para transportar iones a través de la membrana celular en contra de su gradiente electroquímico. Especialmente importante es la llamada bomba de sodio-potasio, que transporta 3 iones de sodio al exterior por cada 2 iones de potasio que bombea al interior de la célula, lo que determina que el medio extracelular sea más rico en Na`. Las diferentes concentraciones de los iones dentro y fuera de la célula determinan que exista una diferencia de potencial entre el interior y el exterior celular que se denomina potencial de membrana. El estado de equilibrio para un ion es aquel en el que el flujo neto de ese ion a través de la membrana es nulo. Esto sucede cuando existe un flujo debido al gradiente de difusión y otro flujo debido a la fuerza electrostática de sentido contrario y ambos son de la misma magnitud. Para cada especie iónica existe un valor del potencial de membrana en el que dicha especie se encuentra en estado de equilibrio. Este valor se conoce como potencial de equilibrio del ion. 4CAPÍTULO 1. DESCRIPCIÓN DEL MODELO DE HODGKIN-HUXLEY Figura 1.2: Gráfico que muestra la membrana celular en estado de reposo. Los canales de fuga (canales pasivos) permanecen abiertos, permitiendo el flujo de iones K`y, en menor medida, Na`a través de la membrana celular, y la bomba de sodio-potasio también regula el paso de iones. Recuperado de http://www.isipedia.com/enfermeria/fisiologia-humana/ propiedades-de-la-membrana-en-reposo Potencial de membrana en reposo En ausencia de un estímulo nervioso, es decir, en una neurona en reposo, el medio extracelular es rico en sodio (Na`), mientras que en el medio intracelular se encuentran fundamentalmente potasio (K`) y aniones orgánicos. Esta distribución asimétrica de los iones determina que exista una diferencia de carga entre el interior y exterior celular, que, en el caso de las neuronas, es del orden de ´70 mV . A esta diferencia de potencial se le llama potencial de membrana en reposo o, de forma abreviada, potencial de reposo, y se dice que la membrana está polarizada. Existe una serie de factores que determinan que, en reposo, el potencial de membrana sea negativo: Los canales de fuga, que son los que están abiertos en reposo, son mucho más permeables al ion K`que al ion Na`, lo que favorece la salida de K`del interior celular, donde es más abundante, provocando que el interior celular sea más negativo. Los aniones orgánicos presentes en el interior de la neurona, que tienen carga negativa, suelen estar unidos a proteínas, constituyendo moléculas voluminosas que no pueden atravesar la membrana celular. 1.1. FISIOLOGÍA DE LA NEURONA. POTENCIAL DE ACCIÓN 5 La bomba de sodio-potasio también contribuye a hacer el interior celular más negativo, ya que extrae 3 cargas positivas por cada 2 que bombea al interior. Todo ello determina que el potencial de reposo de una neurona sea negativo y, además, debido a la mayor permeabilidad al K`, este potencial de reposo, ´70 mV , es similar al potencial de equilibrio del K`,´82 mV . De hecho, si la membrana resultase permeable a un único ion, el potencial de membrana en reposo coincidiría con el potencial de equilibrio de ese ion. El potencial de reposo se mantiene constante, ya que el flujo neto de iones a través de la membrana es nulo cuando la neurona está en reposo, esto es, cuando no emite señal eléctrica. Esto se debe, por una parte, a la fuerza electrostática y al gradiente de difusión y, por otra, a la acción de la bomba de sodio-potasio. Así, se tiene lo siguiente: Ion potasio: Dado que en el interior hay una mayor concentración de K`, el gradiente de difusión determina que los iones de K`se muevan hacia el exterior de la célula. Por otra parte, como el interior de la célula está cargado negativamente respecto al exterior, la acción de la fuerza electrostática determina que los iones se muevan hacia el interior. De esta forma, las fuerzas se compensan y el flujo neto de K`a través de la membrana celular es nulo. Ion sodio: La concentración de Na`es mayor en el exterior que en el interior, por lo que los iones de Na`se mueven hacia el interior de la célula por acción del gradiente de difusión. Además, la fuerza electrostática determina, también, que los iones de Na`se muevan hacia el interior. De esta manera, existe un flujo neto no nulo de Na`hacia el interior celular y este flujo provocaría que el potencial de membrana se desplazase del potencial de reposo. Sin embargo, la actividad de la bomba de sodiopotasio evita que esto suceda, ya que esta proteína emplea energía para, en contra de su gradiente de concentración, transportar 3 iones de sodio al exterior por cada 2 iones de potasio que bombea al interior de la célula, compensando, así, el flujo de Na`que existe hacia el interior celular. De esta manera, cuando la neurona está en reposo, el flujo neto de iones a través de la membrana celular es nulo. Recordemos que la diferencia de cargas entre el interior y el exterior celular se llama potencial de membrana y que, en caso de reposo, se dice que la membrana está polarizada. Diremos que la membrana se despolariza cuando el potencial de membrana se vuelve más 6CAPÍTULO 1. DESCRIPCIÓN DEL MODELO DE HODGKIN-HUXLEY positivo y que la membrana se repolariza cuando el potencial de membrana se vuelve más negativo. La despolarización y repolarización de la membrana celular vienen determinadas por el flujo de iones a través de dicha membrana. En este sentido, uno de los descubrimientos más importantes de Hodgkin y Huxley fue el hecho de que existen canales independientes para los distintos iones, por los cuales estos atraviesan la membrana celular. De esta forma, en un momento dado, la membrana puede ser más permeable a un determinado ion y menos permeable a otro en función del número de canales de cada ion que se encuentren abiertos. Así, por ejemplo, si la membrana aumenta su permeabilidad al Na`y mantiene la permeabilidad al K`, entrarán más iones sodio al interior, de manera que el potencial de membrana será más positivo; la membrana se despolariza. Transmisión del potencial de acción El potencial de acción constituye la forma de transmisión de un estímulo nervioso y se desencadena como respuesta a una despolarización de la membrana celular que provoca que el potencial de membrana alcance un umbral de voltaje [12] determinado, que, en el caso de las neuronas, es de aproximadamente ´55 mV . En este sentido, el potencial de acción responde, de acuerdo a Hodgkin y Huxley, a la ley del todo o nada: si el estímulo eléctrico no es suficiente para que el potencial de membrana rebase el umbral, la respuesta no se genera y, si llega al umbral, el potencial de acción se desarrolla completamente. De esta forma, siempre que se supere el umbral, la amplitud y duración del potencial de acción es la misma, independientemente de la magnitud del estímulo aplicado. Por otra parte, si la magnitud del estímulo es lo suficientemente grande, en lugar de un único potencial de acción puede producirse una secuencia de varios potenciales de acción (disparos neuronales) separados por períodos de inactividad. Así, aunque la intensidad del estímulo no afecta a la forma y velocidad del potencial de acción, produce un aumento proporcional de la frecuencia de disparos neuronales [19]. De esta manera, el sistema nervioso utiliza la frecuencia de los impulsos para codificar la intensidad del estímulo. De acuerdo a los estudios de Hodgkin y Huxley, la transmisión del potencial de acción se puede describir cualitativamente como sigue: 1. Un estímulo eléctrico de suficiente magnitud, procedente de otra neurona u otro tipo de célula, despolariza la membrana, de manera que el potencial de membrana alcanza el umbral de voltaje. 1.1. FISIOLOGÍA DE LA NEURONA. POTENCIAL DE ACCIÓN 7 Figura 1.3: Gráfica potencial de membrana frente a tiempo, obtenida al resolver numéricamente el modelo de Hodgkin-Huxley con Matlab R . 2. Una vez alcanzado el umbral, la membrana aumenta su permeabilidad al Na`, esto es, se abren canales de sodio regulados por voltaje. Dado que la concentración de Na`es mayor en el exterior que en el interior de la célula y debido, también, a las fuerzas electrostáticas, el Na`entra en la neurona y esto provoca que la membrana se despolarice aún más. De esta manera, el interior de la célula pasa a tener más carga que el exterior y el potencial de membrana se vuelve positivo. 3. A medida que avanza la despolarización de la membrana se van abriendo canales de potasio regulados por voltaje. Debido a su gradiente de difusión y a las fuerzas electrostáticas, el K`se mueve del interior al exterior celular. 4. La despolarización continúa hasta que el potencial de membrana se aproxima al potencial de equilibrio del Na`,`40 mV , el potencial de membrana que hace que el flujo neto de Na`sea nulo. En ese momento, se cierran los canales de sodio y comienza la repolarización de la membrana. 5. Los canales de potasio están abiertos y, por tanto, el K`continúa saliendo de la célula. Esto provoca que disminuya la carga en el interior, de manera que el potencial de 8CAPÍTULO 1. DESCRIPCIÓN DEL MODELO DE HODGKIN-HUXLEY membrana va disminuyendo hasta volverse negativo. 6. La repolarización continúa hasta que el potencial de membrana alcanza el potencial de equilibrio del K`, que es ligeramente más negativo que el potencial de reposo. En esta situación, se dice que la membrana está hiperpolarizada. 7. Finalmente, la bomba de sodio-potasio y los canales de fuga (que son más permeables al K`) equilibran la concentración de iones, llevando nuevamente el potencial de membrana al potencial de reposo. Durante la transmisión de un potencial de acción, existe un intervalo de tiempo conocido como período refractario, en el cual los canales regulados por voltaje aún no han recobrado su estado original y la neurona no responde de la misma forma a un estímulo. Comprende dos fases: El período refractario absoluto, que coincide con la fase de despolarización y parte de la fase de repolarización y en el cual los canales de Na`regulados por voltajes están inactivados. En este período no puede generarse un nuevo potencial de acción, independientemente de la intensidad del estímulo aplicado. Esto se debe a que, a diferencia de los canales en estado cerrado, los canales inactivados no pueden abrirse aunque reciban estímulos. Para permitir su apertura, los canales de sodio tienen que reactivarse, esto es, pasar de estado inactivado a estado cerrado y este hecho sucede durante la repolarización. El período refractario relativo, que coincide con la fase en que la membrana está hiperpolarizada. En este período los canales de Na`regulados por voltaje ya están cerrados y puede generarse un nuevo potencial de acción. Sin embargo, para generar este segundo potencial de acción es necesario que se aplique un estímulo de gran magnitud, ya que la despolarización necesaria para alcanzar el umbral es mayor que en el caso del primer potencial de acción. 1.2. Modelo matemático Como ya se comentó en la sección anterior, a través de sus experimentos con el axón de calamar gigante, Hodgkin y Huxley lograron mostrar que las corrientes iónicas a través de la membrana neuronal son las responsables de la transmisión del impulso nervioso. Además, observaron que el comportamiento de la membrana celular puede ser descrito en términos de un circuito eléctrico. De esta manera, a través de las ecuaciones que modelan dicho circuito, lograron formalizar la transmisión del potencial de acción en un sistema de 1.2. MODELO MATEMÁTICO 9 ecuaciones diferenciales. Antes de presentar el modelo, hay algunas consideraciones que cabe tener en cuenta: Para poder medir las corrientes iónicas, Hodgkin y Huxley emplearon una técnica conocida como voltage clamp (fijación de voltaje), introducida por Kenneth Cole y George Marmont en 1947. Mediante el uso de dos electrodos y un circuito de retroalimentación, esta técnica permite mantener fijo el potencial de membrana. De esta forma, se puede estudiar el comportamiento de los canales iónicos para distintos potenciales de membrana, previamente fijados. Las ecuaciones del modelo de Hodgkin y Huxley fueron derivadas del axón de calamar gigante, un axón no mielinizado. Un axón mielinizado está recubierto por una capa proteica conocida como vaina de mielina, que actúa como un aislante eléctrico. Debido a la presencia de canales en la membrana celular, en la transmisión del impulso nervioso hay una fuga de corriente hacia el medio circundante y esto provoca que la amplitud de la señal se disipe. Para evitar esto, algunas neuronas presentan axones recubiertos con mielina. Este recubrimiento no es continuo, sino que hay segmentos mielinizados intercalados con otros no mielinizados. Estos últimos se llaman nodos de Ranvier (ver figura 1.1) y es donde están concentrados los canales iónicos, de tal forma que, debido a la gran resistencia que ofrece la membrana en los segmentos mielinizados, la transmisión del potencial de acción tiene lugar en estos nodos. Se dice, en este caso, que la conducción del impulso nervioso es saltatoria y la velocidad de transmisión es mayor que en el caso no mielinizado [5]. El modelo de Hodgkin-Huxley se puede adaptar para el caso de axones mielinizados, pero en este trabajo nos centraremos en el caso de axones no mielinizados, siguiendo lo estudiado por Hodgkin y Huxley. 1.2.1. Circuito eléctrico equivalente Las propiedades eléctricas de una neurona pueden representarse mediante un circuito eléctrico simple formado por un condensador, resistores y baterías [12]. La membrana celular, que separa el interior y el exterior celular, actúa como un condensador. Un condensador está formado por dos placas conductoras y un medio aislante que las separa, de tal manera que, cuando se someten a una diferencia de potencial, una de las placas adquiere carga positiva y la otra negativa. Así, se dice que un condensador 10 CAPÍTULO 1. DESCRIPCIÓN DEL MODELO DE HODGKIN-HUXLEY ` ´ENa RNa Cm RK ` ´ EK RL ` ´ EL Medio extracelular Medio intracelular Figura 1.4: Adaptación del circuito original propuesto por Hodgkin y Huxley [12] para representar el comportamiento eléctrico de la membrana del axón de calamar gigante. almacena cargas eléctricas y está caracterizado por su capacitancia, que se mide en faradios (F). La carga en un condensador se calcula como q“CV donde Ces la capacitancia y V la diferencia de potencial entre las placas. Se considera la membrana celular, pues, como un condensador de capacitancia Cma. Su carga será q“CmVm, en donde Vmes el potencial de membrana. Los distintos tipos de canales iónicos que se encuentran en la membrana actúan como resistores caracterizados por su resistencia, que se mide en ohmios (Ω). Hay un resistor para cada tipo de canal iónico: uno para el sodio, de resistencia RNa; otro para el potasio, de resistencia RKy otro de escape para el resto de iones (fundamentalmente Cl´), con resistencia RL. Los canales de escape se suponen constantes; sin embargo, el número de canales abiertos para el sodio y el potasio depende del potencial de membrana. Es- to se traduce en que RNa yRKson funciones dependientes del tiempo y del potencial de membrana, Vm, lo que se indica en el circuito con una flecha transversal sobre los resistores. La inversa de la resistencia es la conductancia, gion “1 Rion , que es una medida de la permeabilidad de la membrana a un determinado ion. Denotaremos por gNa “gNapt, Vmq, aCmdenota la capacitancia específica, es decir, el valor por unidad de área, y se mide en µF {cm2. 1.2. MODELO MATEMÁTICO 11 gK“gKpt, VmqygLlas conductancias para el sodio, el potasio y de escape, respectivamente, donde la última se supone constante. Por otra parte, la diferencia de potencial entre el potencial de membrana y el potencial de equilibrio para cada ion induce una fuerza electromotriz, que se representa por una batería en el circuito. Así, se tienen baterías distintas para el sodio, el potasio y el resto de iones, con voltajes respectivos ENa,EKyEL. Cuando una corriente Ipx, tqse inyecta en la célula, esta puede aumentar la carga en la membrana (condensador) o bien fluir a través de los canales iónicos (resistencias). De esta forma, la corriente total a través de la membrana, Ipx, tqb, tiene dos componentes: la corriente capacitiva, ICpx, tq, que se almacena en la membrana, y la corriente iónica, Iionpx, tq, asociada al flujo de iones a través de los resistores: Ipx, tq “ ICpx, tq ` Iionpx, tq “ ICpx, tq ` INapx, tq ` IKpx, tq ` ILpx, tq.(1.1) La corriente capacitiva, ICpx, tq, se define como la tasa de cambio de la carga almacenada en la membrana , IC“Bq Bt. Además, teníamos que la carga en la membrana se calculaba como q“CmVm. Por tanto, reescribimos la corriente capacitiva como ICpx, tq “ Cm BVm Bt. Por otra parte, por la ley de Ohm,I“V R, lo que, en términos de conductancia, sería I“V g. Ahora, como Vmes el voltaje total a través de la membrana celular y Ekes el voltaje de la batería correspondiente al ion k, el voltaje a través del resistor correspondiente será Vm´Ek, donde Ekes el potencial de equilibrio para el ion k. Así, aplicando la ley de Ohm, tenemos: INa “gNapVm´ENaq, IK“gKpVm´EKq, IL“gLpVm´ELq. (1.2) Sustituyendo en (1.1) las expresiones anteriores y la expresión para la corriente capabIpx, tqdenota la densidad de corriente total a través de la membrana y se mide en µA{cm2. 18 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... tenemos el siguiente problema de valor inicial (PVI): Hallar ~ Xptq,0ďtďT, para algún T, donde d~ X dt “~ Fpt, ~ Xqcon ~ Xp0q “ ~ X0. En relación a la condición inicial, en tiempo t“0suponemos que la neurona no está sometida a ningún estímulo y, por tanto, el potencial de membrana coincide con el potencial de reposo; Vmp0q “ Vr“ ´70. Además, los parámetros mp0q,hp0qynp0qpueden aproximarse por sus valores en estado estacionario m8pVmp0qq,h8pVmp0qq yn8pVmp0qq, tal y como puede verse en el Apéndice II. El sistema de ecuaciones diferenciales ordinarias (EDO) anterior es un sistema complejo y no lineal, por lo que no puede resolverse de forma analítica. Por tanto, es necesario considerar métodos numéricos de resolución de EDO para obtener soluciones aproximadas. 2.1. Existencia y unicidad de solución del PVI Antes de abordar la resolución numérica de un PVI es conveniente plantearse si dicho problema tiene solución y, en caso de que la tenga, si esta es única. Con este objetivo presentaremos un teorema que, bajo ciertas hipótesis, proporciona existencia y unicidad de solución de un PVI. Teorema 2.1 (de Picard-Lindelöf).Sean I“ rt0, T sy~ f:IˆRnÑRnuna función continua. Consideremos el PVI: $ & % d~xptq dt “~ fpt, ~xptqq tPI ~xpt0q “ ~x0 Si ~ fes lipschitziana respecto de la segunda variable, esto es, DLą0tal que k~ fpt, ~x1q ´ ~ fpt, ~x2qkďLk~x1´~x2k@~x1, ~x2PRny@tPI, entonces existe solución única,~xptq, del PVI en el intervalo I. Además, ~xptq P C1pIq. Observación 2. Si ~ f:IˆRnÑRn, con I“ rt0, Ts, es continua y todas las derivadas parciales Bfi Bxj con 1ďi, j ďnexisten y están acotadas en Ω“˚ IˆRn, entonces, por el Teorema del Valor Medio,~ fes lipschitziana respecto de la segunda variable con constante 2.1. EXISTENCIA Y UNICIDAD DE SOLUCIÓN DEL PVI 19 de Lipschitz: L“sup pt,~xqPΩ    B~ fpt, ~xq B~x   ă `8 Nótese que la derivada parcial Bfi Bxj es la entrada pi, jqde la matriz jacobiana B~ f B~x y que la norma considerada, }¨}, es la norma matricial subordinada a la norma }¨}de Rn. Corolario 2.2. Sea I“ rt0, Tsy sea ~ f:IˆRnÑRntal que ~ fPC1pIˆRnq. Entonces el PVI: $ & % d~xptq dt “~ fpt, ~xptqq tPI ~xpt0q “ ~x0 tiene solución única ~x PC1pIq. A continuación, se discutirá si el PVI que define el modelo de Hodgkin-Huxley en hipótesis de homogeneidad en espacio admite solución única. Para ello, se comprobará que se satisfacen las hipótesis del corolario anterior. En primer lugar, veamos que ~ Fpt, ~ Xqes continua en todo punto pt, ~ Xq P r0, Ts ˆ R4, ya que cada una de sus componentes es continua: F1pt, ~ Xq “ Iptq{Cm´1{Cm`ĚgNaX23X3pX1´ENaq ` ĎgKX44pX1´EKq ` ĎgLpX1´ELq˘ es claramente continua en todo pt, ~ Xq P r0, T s ˆ R4. F2pt, ~ Xq “ αmpX1qp1´X2q ´ βmpX1qX2,F3pt, ~ Xq “ αhpX1qp1´X3q ´ βhpX1qX3y F4pt, ~ Xq “ αnpX1qp1´X4q ´ βnpX1qX4son claramente continuas en todo pt, ~ Xq P r0, TsˆR4si αm, βm, αh, βh, αnyβnson continuas en la variable X1“Vm. En efecto, βm, αh, βhyβnson continuas, pero αmyαnpresentan una discontinuidad evitable en Vm“ ´45 yVm“ ´60, respectivamente. Podemos, no obstante, para garantizar la continuidad, considerar, en lugar de las expresiones propuestas por Hodgkin y Huxley, sus extensiones continuas: αmpVmq “ $ ’ ’ ’ & ’ ’ ’ % 0,1pVm`45q 1´exp ˆ´Vm`45 10 ˙si Vm‰ ´45 1 si Vm“ ´45 (2.5) 20 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... αnpVmq “ $ ’ ’ ’ & ’ ’ ’ % 0,01 Vm`60 1´exp ˆ´Vm`60 10 ˙si Vm‰ ´60 0,1 si Vm“ ´60 (2.6) A continuación, se debe comprobar que se satisface la segunda hipótesis del corolario (2.2), esto es, que las derivadas parciales existen y están acotadas. Para ello, se verá que las derivadas parciales son continuas en cualquier compacto conteniendo a la condición inicial y, por el Teorema de Weierstrass, se concluirá que están acotadas. De acuerdo a (2.1) el sistema de ecuaciones diferenciales acopladas que define el modelo puede escribirse como una ecuación diferencial vectorial: d~ X dt “~ Fpt, ~ Xq, donde ~ X“¨ ˚ ˚ ˚ ˚ ˝ X1 X2 X3 X4 ˛ ‹ ‹ ‹ ‹ ‚ y~ Fpt, ~ Xq “ ¨ ˚ ˚ ˚ ˚ ˝ F1pt, ~ Xq F2pt, ~ Xq F3pt, ~ Xq F4pt, ~ Xq ˛ ‹ ‹ ‹ ‹ ‚ . La matriz jacobiana asociada al campo vectorial ~ Fes: ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˝ BF1pt, ~ Xq BX1 BF1pt, ~ Xq BX2 BF1pt, ~ Xq BX3 BF1pt, ~ Xq BX4 BF2pt, ~ Xq BX1 BF2pt, ~ Xq BX2 0 0 BF3pt, ~ Xq BX1 0BF3pt, ~ Xq BX3 0 BF4pt, ~ Xq BX1 0 0 BF4pt, ~ Xq BX4 ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‚ , con BF1pt, ~ Xq BX1 “ ´ 1 Cm´ĚgNaX23X3`ĎgKX44`ĎgL¯, BF1pt, ~ Xq BX2 “ ´ 3 CmĚgNaX22X3pX1´ENaq, BF1pt, ~ Xq BX3 “ ´ 1 CmĚgNaX23pX1´ENaq, 2.1. EXISTENCIA Y UNICIDAD DE SOLUCIÓN DEL PVI 21 BF1pt, ~ Xq BX4 “ ´ 4 CmĎ gKX43pX1´EKq, BF2pt, ~ Xq BX1 “α1 mpX1qp1´X2q ´ β1 mpX1qX2“ “ $ ’ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ ’ % 2 9exp ˆ´X1`70 18 ˙`0.1p1´X2q 1´exp ˆ´X1`45 10 ˙´ 0.01 pX1`45qexp ˆ´X1`45 10 ˙p1´X2q ˆ1´exp ˆ´X1`45 10 ˙˙2si X1‰ ´45 0.05p1´X2q ` 4 18 exp ˆ´25 18˙si X1“ ´45 BF2pt, ~ Xq BX2 “ ´αmpX1q ´ βmpX1q “ $ ’ ’ ’ ’ & ’ ’ ’ ’ % ´0.1pX1`45q 1´exp ˆ´X1`45 10 ˙´4 exp ˆ´X1`70 18 ˙si X1‰ ´45 ´1´4 exp ˆ´25 18˙si X1“ ´45 BF3pt, ~ Xq BX1 “α1 hpX1qp1´X3q ´ β1 hpX1qX3“ ´0.0035 exp ˆ´X1`70 20 ˙p1´X3q ´ 0.1 exp ˆ´X1`40 10 ˙X3 ˆ1`exp ˆ´X1`40 10 ˙˙2, BF3pt, ~ Xq BX3 “ ´αhpX1q ´ βhpX1q “ ´0.07 exp ˆ´X1`70 20 ˙´1 1`exp ˆ´X1`40 10 ˙, BF4pt, ~ Xq BX1 “α1 npX1qp1´X4q ´ β1 npX1qX4“ “ $ ’ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ ’ % 0.016 exp ˆ´X1`70 80 ˙X4`0.01p1´X4q 1´exp ˆ´X1`60 10 ˙´ 0.001pX1`60qexp ˆ´X1`60 10 ˙p1´X4q ˆ1´exp ˆ´X1`60 10 ˙˙2si X1‰ ´60 0.005p1´X4q ` 0.016 exp ˆ´10 80˙si X1“ ´60 BF4pt, ~ Xq BX4 “ ´αnpX1q ´ βnpX1q “ “$ ’ ’ ’ ’ & ’ ’ ’ ’ % ´0.01 X1`60 1´exp ˆ´X1`60 10 ˙´0.125 exp ˆ´X1`70 80 ˙si Vm‰ ´60 ´0.1´0.125 exp ˆ´10 80˙si Vm“ ´60 Resulta sencillo comprobar que las derivadas parciales anteriores son continuas en IˆR4, por lo que, en particular, son continuas en cualquier entorno compacto KĂIˆR4 22 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... de pt0, ~x0q. Por tanto, se tiene que ~ FPC1pKq, con Kun compacto conteniendo a la condición inicial, y, por el corolario (2.2), se concluye que el PVI considerado admite solución única. 2.2. Métodos numéricos de resolución de EDO A continuación se presentarán distintos métodos numéricos para resolver EDO aplicándolos a las ecuaciones del modelo. 2.2.1. Euler explícito Es el método más sencillo de implementar. Geométricamente, se basa en aproximar la curva solución por una poligonal que se construye a través de rectas tangentes. Sea el PVI: $ & % d~ X dt “~ Fpt, ~ XqtP r0, Ts ˆ Rn, ~ Xp0q “ ~x0. Comenzamos considerando una discretización del intervalo r0, Tscon N`1puntos equiespaciados: 0“t0ăt1ă ¨ ¨ ¨ ă tN´1ătN“T. Denotamos por h“T´0 Nel paso de discretización, esto es, la distancia entre dos nodos consecutivos de la discretización. Supongamos que ~ XPC2pr0, Tsq y consideremos el polinomio de Taylor centrado en tncon la expresión integral del resto: ~ Xptn`1q “ ~ Xptnq ` ~ X1ptnqptn`1´tnq ` żtn`1 tn ptn`1´sq~ X2psqds Observación 3. Se tiene la siguiente acotación:   żtn`1 tn ptn`1´sq~ X2psqds  ď   ~ X2  C0żtn`1 tn ptn`1´sqds “1 2    ~ X2  C0ptn`1´tnq2 Ahora, dado que h“tn`1´tny que tiene que satisfacerse la EDO ~ X1“~ Fpt, ~ Xqen todos los puntos de la discretización, tenemos: ~ Xptn`1q “ ~ Xptnq ` h~ Fptn,~ Xptnqq ` Oph2q lo que sugiere el esquema numérico: #~ X0“~x0, ~ Xn`1“~ Xn`h~ Fptn,~ Xnqn“0,1,¨ ¨ ¨ , N ´1, 2.2. MÉTODOS NUMÉRICOS DE RESOLUCIÓN DE EDO 23 donde ~ Xnes la aproximación de la solución exacta en el nodo n-ésimo, ~ Xn«~ Xptnq. Observación 4. Este método numérico es explícito y converge con orden 1. Además, es un esquema de 1 paso, ya que, para hallar la aproximación en tiempo tn`1, basta conocer la aproximación en tiempo tn. 2.2.2. Runge-Kutta Sea el PVI: $ & % d~ X dt “~ Fpt, ~ XqtP r0, Ts, ~ Xp0q “ ~x0. Este problema se puede reescribir en términos de una ecuación integral como: ~ Xptq ´ ~x0“żt 0 ~ X1psqds ô~ Xptq ´ ~x0“żt 0 ~ Fps, ~ Xpsqq ds Consideramos una discretización del intervalo r0, Tscon N`1puntos equiespaciados 0“t0ăt1ă ¨ ¨ ¨ ă tN´1ătN“Tyh“T´0 Nel paso de discretización. Los métodos de Runge-Kutta son una familia de métodos de un paso que pueden ser implícitos o explícitos. En su formulación más general, un método de Runge-Kutta de s etapas se basa en, conocida ~ Xn, aproximación de la solución exacta en el nodo tn, hallar ~ Xn`1con tn`1“tn`hempleando una fórmula de cuadratura de snodos para aproximar la solución de la ecuación integral ~ Xptn`1q ´ ~ Xptnq “ ştn`1 tn ~ Fpt, ~ Xptqqdt. Así, dado ~ Xn, para calcular ~ Xn`1se emplean spuntos auxiliares tn,1,¨ ¨ ¨ , tn,s con tn,i “tn`hci,i“1,¨ ¨ ¨ , s y se tiene que ~ Xn`1“~ Xn`hřs i“1bi~ kn,i donde hbison los pesos de la fórmula de cuadratura y ~ kn,i “~ Fptn,i,~ Xn`hřs j“1aij~ kn,jq,i“1,¨ ¨ ¨ , s. Los coeficientes aij ,biycicon i, j “1,¨ ¨ ¨ , s que definen el esquema suelen presentarse en una tabla conocida como tabla de Butcher: c1a11 a12 ¨ ¨ ¨ a1s c2a21 a22 ¨ ¨ ¨ a2s . . .. . .. . .. . .. . . csas1as2¨ ¨ ¨ ass b1b2¨ ¨ ¨ bs 24 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... Finalmente, para un Runge-Kutta de setapas se tiene el siguiente esquema: $ ’ ’ & ’ ’ % ~ X0“~x0, ~ Xn`1“~ Xn`hřs i“1bi~ kn,i n“0,1,¨ ¨ ¨ , N ´1, con ~ kn,i “~ Fptn,i,~ Xn`hřs j“1aij~ kn,jq, tn,i “tn`hci, i “1,¨ ¨ ¨ , s. Observación 5. Resulta sencillo comprobar que si la matriz paijqi,j“1,¨¨¨ ,s es triangular inferior con todos los elementos de la diagonal principal iguales a cero, aij “0@jěi, esto es, estrictamente triangular inferior, entonces el esquema de Runge-Kutta correspondiente es explícito. Observación 6. Geométricamente, un método de Runge-Kutta de setapas emplea como dirección de avance, para aproximar la solución en un nodo, la dada por una recta cuya pendiente es una ponderación de las pendientes calculadas mediante spuntos auxiliares. Runge Kutta clásico de 4 etapas Los coeficientes que definen este método vienen dados por la tabla de Butcher siguiente: 0 0 0 0 0 1/2 1/2 0 0 0 1/2 0 1/2 0 0 1 0 0 1 0 1/6 1/3 1/3 1/6 Así, ~ X0“~x0y, para 0ănăN´1, se tiene que ~ Xn`1“~ Xn`h 4 ÿ i“1 bi~ kn,i “ ~ Xn`h 6p~ kn,1`2~ kn,2`2~ kn,3`~ kn,4q, donde $ ’ ’ ’ ’ & ’ ’ ’ ’ % ~ kn,1“~ Fptn,1,~ Xnq “ ~ Fptn,~ Xnq, ~ kn,2“~ Fptn,2,~ Xn`h1 2~ kn,1q “ ~ Fptn`h1 2,~ Xn`h1 2~ kn,1q, ~ kn,3“~ Fptn,3,~ Xn`h1 2~ kn,2q “ ~ Fptn`h1 2,~ Xn`h1 2~ kn,2q, ~ kn,4“~ Fptn,4,~ Xn`h~ kn,3q “ ~ Fptn`h, ~ Xn`h~ kn,3q. Observación 7. El método de Runge-Kutta clásico de 4 etapas es explícito y converge con orden 4. 2.3. ANÁLISIS COMPARATIVO DE LOS MÉTODOS NUMÉRICOS 25 2.3. Análisis comparativo de los métodos numéricos Con el objetivo de evaluar la bondad de la aproximación que proporcionan los distintos métodos numéricos, resulta interesante aplicarlos a un caso particular más sencillo para el cual es posible obtener la solución exacta. Si se consideran las conductancias para Na y K iguales a 0 y la intensidad de corriente externa, I, constante a lo largo del tiempo, la ecuación: Iptq “ Cm dVm dt `Ě gNam3hpVm´ENaq ` Ď gKn4pVm´EKq ` Ď gLpVm´ELq(2.7) se reduce a I“Cm dVm dt `Ď gLpVm´ELq(2.8) Por tanto, se tendría el siguiente PVI: $ & % dVm dt “1 Cm rI´ĎgLpVm´ELqs , Vmp0q “ ´70, donde la ecuación diferencial puede expresarse como una ecuación en variables separadas: 1 I´Ď gLpVm´ELqdVm´1 Cm dt “0ô1 I´Ď gLpVm´ELqdVm“1 Cm dt En primer lugar, se integran ambos lados de la ecuación: ż1 I´Ď gLpVm´ELqdVm“ż1 Cm dt ô ´ 1 Ď gL ln |I´Ď gLVm`Ď gLEL| “ 1 Cm t`k A continuación, se sustituye la condición inicial, Vmp0q“´70, obteniéndose el valor de la constante k: ´1 Ď gL ln |I`ĎgL70 `ĎgLEL| “ k. Finalmente, sustituyendo el valor de k, se tiene que: ´1 ĎgL ln |I´ĎgLVm`ĎgLEL| “ 1 Cm t´1 ĎgL ln |I`ĎgL70 `ĎgLEL| ô exp ˆ´1 Cm t`1 Ď gL ln |I`Ď gL70 `Ď gLEL|˙“exp ˆ1 Ď gL ln |I´Ď gLVm`Ď gLEL|˙ô expp´ 1 Cm tqpI`ĎgL70 `ĎgLELq1{ĎgL“ pI´Ď gLVm`Ď gLELq1{ĎgLô expp´ ĎgL Cm tqpI`ĎgL70 `ĎgLELq “ I´ĎgLVm`ĎgLELô 26 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... Vm“1 Ď gL r´ expp´ Ď gL Cm tqpI`Ď gL70 `Ď gLELq ` I`Ď gLELs Conocida la solución analítica para el problema particular planteado, se puede calcular el error absoluto para las aproximaciones que proporcionan los distintos métodos numéricos. Los errores máximos obtenidos se muestran en la siguiente tabla a modo de comparativa: Método Error máximo Orden de convergencia Euler explícito 1.34 ˆ10´11 Runge-Kutta clásico de 4 etapas 5.93 ˆ10´12 4 Tabla 2.1: Comparación métodos numéricos aplicados al modelo de Hodgkin-Huxley en el intervalo [0,150] con conductancias de Na y K nulas y corriente externa constante I=4. Se consideran 15000 pasos de tiempo. A la vista de la tabla anterior, resulta claro que el método Runge-Kutta clásico de 4 etapas es más preciso que el Euler explícito. Además, en este caso particular, con conductancias de Na y K nulas y corriente externa constante, en el que no se desarrollan potenciales de acción (ver figura (2.1)), los métodos explícitos considerados resultan numéricamente estables para pasos de tiempo de tamaño considerable, del orden de 1ms. En efecto, en la siguiente tabla se presentan, a modo de ejemplo, los errores máximos obtenidos aplicando el método Runge-Kutta clásico de 4 etapas para distintos pasos de tiempo: Runge-Kutta clásico de 4 etapas Paso de tiempo Error máximo 1ms 7.60 ˆ10´4 0,1ms 6.09 ˆ10´8 0,01 ms 5.94 ˆ10´12 Se puede observar que el método es numéricamente estable para los pasos de tiempo considerados y que la reducción del error al disminuir el paso de tiempo es la esperada para un método con orden de convergencia 4. 2.3. ANÁLISIS COMPARATIVO DE LOS MÉTODOS NUMÉRICOS 27 Figura 2.1: Respuesta neuronal para el caso particular con conductancias de Na y K nulas e intensidad constante I“4. No se generan potenciales de acción. Sin embargo, si se considera el PVI general (2.1), pueden ser necesarios pasos de tiempo pequeños para garantizar la estabilidad numérica de los métodos considerados, tal y como veremos a continuación. Análisis de la rigidez para el PVI general La gran velocidad con la que se desarrolla un potencial de acción determina que se produzca una rápida variación en el potencial de membrana, dando lugar a zonas de alto gradiente en la solución que coinciden con los intervalos de tiempo donde se producen disparos neuronales. Este hecho fuerza a tener que considerar pasos de tiempo muy pequeños para garantizar la estabilidad numérica de los métodos numéricos considerados, aumentando con ello el coste computacional. En este sentido, aun cuando no sea necesaria demasiada exactitud en la solución, métodos explícitos como los considerados precisan de pasos de tiempo pequeños, del orden de 0,01 ms, para que quede garantizada su estabilidad numérica. Se dice, en este caso, que el PVI presenta rigidez. De esta forma, un problema rígido (en inglés, stiff ) se puede definir como aquel para el cual determinados métodos numéricos resultan inestables a menos que se consideren pasos de tiempo muy pequeños. A la vista de lo anterior, puede resultar de interés considerar métodos de paso adaptativo que permitan tomar pasos de tiempo más pequeños en las zonas donde se desarrollan potenciales de acción y más grandes en otras. Con este objetivo, se emplean las funcio- 34 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... En la figura 2.5 se observa que, para un estímulo discontinuo de intensidad 3 µA{cm2 con período 10 ms y duración de los impulsos 10 msb, se genera una ráfaga de potenciales de acción. Sin embargo, si se disminuye el período de los impulsos hasta un valor inferior a 5 ms, se tiene que, manteniendo la magnitud y duración de los impulsos, se desarrolla un único potencial de acción al inicio (ver figura 2.6). Esto se debe a que el tiempo entre impulsos, 3msc, es inferior a la duración del período refractario relativo (figura 2.4) que sucede a un potencial de acción. Durante este período, el umbral de voltaje necesario para que se desencadene un potencial de acción es mayor que en situación de reposo, por lo que para que se generen nuevos potenciales de acción es necesario aumentar la magnitud del estímulo. Este hecho se muestra en la figura 2.7 en la que se aplica un estímulo de intensidad I“4µA{cm2, desencadenando, de esta forma, una ráfaga de potenciales de acción. Figura 2.5: Respuesta neuronal a una corriente externa discontinua de intensidad 3µA{cm2, duración de los impulsos 10 ms y período 10 ms. b Iptq “ #3µA{cm2tP r10i, 10 `10iq, i “0,2,4,6,¨ ¨ ¨ 0µA{cm2tP r10i, 10 `10is, i “1,3,5,7,¨ ¨ ¨ c Iptq “ #3µA{cm2tP r13i, 10 `13iq, i “0,1,2,3,¨ ¨ ¨ 0µA{cm2tP r10 `13i, 13 `13is, i “0,1,2,3,¨ ¨ ¨ 2.4. ANÁLISIS DE LOS RESULTADOS DE LA SIMULACIÓN NUMÉRICA 35 Figura 2.6: Respuesta neuronal a una corriente externa discontinua de intensidad 3µA{cm2, duración de los impulsos 10 ms y período 3 ms. Figura 2.7: Respuesta neuronal a una corriente externa discontinua de intensidad 4µA{cm2, duración de los impulsos 10 ms y período 3 ms. 36 CAPÍTULO 2. RESOLUCIÓN DEL MODELO BAJO HIPÓTESIS... Capítulo 3 Resolución del modelo sin hipótesis de homogeneidad en espacio Con el objetivo de estudiar cómo es la propagación del potencial de acción a lo largo del axón de una neurona, en este capítulo se considera el modelo de Hodgkin y Huxley sin el supuesto de fijación espacial. Si se supone que el axón es cilíndrico, puede hacerse uso de la llamada ecuación del cable, tal y como se vio en el capítulo 1, resultando el siguiente sistema de ecuaciones en derivadas parcialesa: Cm BV Bt“a 2R B2V Bx2´ĚgNam3hpV´ENaq ´ ĎgKn4pV´EKq ´ ĎgLpV´ELq, Bm Bt“αmpVqp1´mq ´ βmpVqm, Bh Bt“αhpVqp1´hq ´ βhpVqh, Bn Bt“αnpVqp1´nq ´ βnpVqn, (3.1) donde aes el radio del axón (en cm) y Rla resistividad axial específica (en Ω¨cm), que se consideran constantes a lo largo del axón. En concreto, se tomarán los valores que emplearon Hodgkin y Huxley [12]: a“0.0238 cm yR“35.4 Ω ¨cm. Por medio de este sistema es posible estudiar no solo la variación del potencial de membrana con respecto al tiempo, sino también cómo varía en función de la distancia longitudinal a lo largo del axón. Para ello, se considera una región de estudio de longitud aPara aligerar la notación, en este capítulo se denota el potencial de membrana por V, en lugar de Vm. 37 38 CAPÍTULO 3. RESOLUCIÓN DEL MODELO SIN HIPÓTESIS... Ly un intervalo temporal r0, Tsb. Para que el sistema de EDP esté bien definido se debe proporcionar el valor de V,m,h ynen cada punto del axón en tiempo t“0(condición inicial) y especificar, también, cómo es el potencial de membrana en los extremos de la región de estudio del axón (condiciones de contorno). Condición inicial En el instante inicial, la neurona se encuentra en reposo, por lo que V“Vr“ ´70, y los parámetros de activación y desactivación se aproximan por los correspondientes al estado estacionario (ver Apéndice II). Así, se tiene que: Vpx, 0q“´70,(3.2) mpx, 0q “ m8p´70q “ αmp´70q αmp´70q ` βmp´70q,(3.3) hpx, 0q “ h8p´70q “ αhp´70q αhp´70q ` βhp´70q,(3.4) npx, 0q “ n8p´70q “ αnp´70q αnp´70q ` βnp´70q,(3.5) con xP r0, Ls. Condiciones de contorno Dado que la primera ecuación del sistema es una EDP parabólica, que involucra a la segunda derivada espacial, es necesario proporcionar dos condiciones de contorno, una por cada extremo de la región de estudio del axón. Una primera opción consiste en especificar el valor del potencial de membrana en los extremos del intervalo espacial, esto es, considerar condiciones de tipo Dirichlet: Vp0, tq “ V0ptq,(3.6) VpL, tq “ VLptq.(3.7) bEn caso de que se quiera estudiar solamente la variación espacial, se puede considerar la neurona en un estado estacionario. Así, la única variable independiente del sistema es x, por lo que se trata de un sistema de EDO, y las funciones m,hynse reemplazan por las correspondientes al estado estacionario m8,h8yn8. 39 Este tipo de condiciones de contorno son las más sencillas de implementar. Sin embargo, en este trabajo se considerarán condiciones de contorno de tipo Neumann, en las cuales se especifica el valor de la derivada espacial, ya que tienen una interpretación más natural desde el punto de vista de la fisiología de la neurona. En efecto, dado que BV Bxes proporcional a la corriente longitudinal que circula a través del axón, imponer condiciones Neumann equivale a proporcionar el valor de la corriente longitudinal en los extremos de la región de estudio. Si se supone que en x“0se inyectan I0µA y en x“Lno hay salida de corriente (extremo aislado), se tienen las siguientes condiciones de contorno: BV Bxp0, tq“´RI0 πa2,(3.8) BV BxpL, tq “ 0.(3.9) Veamos cómo se obtiene la primera de estas condiciones de contorno [7]. aIlong ∆x Im Figura 3.1: Representación del axón de una neurona donde aes el radio, Ilong es la corriente que circula a través del axón e Imes la corriente a través de la membrana. Dada una sección del axón, la resistencia total que ejerce el axón en esa región es directamente proporcional a la longitud de la sección e inversamente proporcional al área transversal, siendo la resistividad axial específica, R, la constante de proporcionalidad. Así, para una sección de longitud ∆x(ver figura 3.1), se tiene que: RT“R∆x πa2.(3.10) Ahora, por la ley de Ohm, la variación de potencial en la sección es igual al producto de la intensidad de corriente que atraviesa esa sección por la resistencia, esto es: Vpx`∆x, tq ´ Vpx, tq“´IlongRT“ ´IlongR∆x πa2,(3.11) 40 CAPÍTULO 3. RESOLUCIÓN DEL MODELO SIN HIPÓTESIS... donde el signo menos se debe al convenio de que el sentido positivo de la corriente es el dado por el flujo de cargas positivas. Finalmente, cuando ∆xÑ0, se obtiene la siguiente expresión: BV Bxpx, tq“´IlongR πa2.(3.12) A la vista de la expresión (3.12), si en el extremo izquierdo de la región de estudio se inyecta una corriente de intensidad I0, como esta corriente es longitudinal se tiene: BV Bxp0, tq“´RI0 πa2,(3.13) que es la condición de contorno que se considera en x“0. Se tiene, entonces, el siguiente problema de contorno: $ ’ ’ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ ’ ’ % BV Bt´a 2RCm B2V Bx2“ ´ 1 Cm ¯gNam3hpV´ENaq ´ 1 Cm ¯gKn4pV´EKq ´ 1 Cm ¯gLpV´ELq, Vpx, 0q“´70, x P r0, Ls; BV Bxp0, tq“´RI0 πa2,BV BxpL, tq “ 0, t P r0, T s. $ ’ ’ ’ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ ’ ’ ’ % Bm Bt“αmpVqp1´mq ´ βmpVqm, mpx, 0q “ m8p´70q, x P r0, Ls; Bh Bt“αhpVqp1´hq ´ βhpVqh, hpx, 0q “ h8p´70q, x P r0, Ls; Bn Bt“αnpVqp1´nq ´ βnpVqn, npx, 0q “ n8p´70q, x P r0, Ls. El sistema de EDP anterior es no lineal, por lo que, en general, no puede ser resuelto de forma analítica. Por ello, se buscará aproximar su solución empleando métodos numéricos de resolución de EDP. 3.1. Métodos numéricos de resolución de EDP En la resolución de EDP se emplean fundamentalmente dos tipos de métodos: los métodos de diferencias finitas y los de elementos finitos. En este trabajo se emplearán métodos 3.1. MÉTODOS NUMÉRICOS DE RESOLUCIÓN DE EDP 41 de diferencias finitas debido a su mayor facilidad de implementación. Los métodos de diferencias finitas se basan en la idea de transformar el problema de resolver una EDP en el de resolver un sistema de ecuaciones algebraicas. Para ello, se discretiza el intervalo temporal y el dominio espacial y se emplean fórmulas de derivación numérica de tipo interpolatorio polinómico para aproximar las derivadas. Observación 10. Para hallar la solución discreta del sistema se emplearán las siguientes fórmulas de derivación numérica: Fórmula centrada para la derivada primera: f1pαq “ fpα`hq ´ fpα´hq 2h`Oph2q(3.14) Fórmula regresiva para la derivada primera: f1pαq “ fpαq ´ fpα´hq h`Ophq(3.15) Fórmula progresiva para la derivada primera: f1pαq “ fpα`hq ´ fpαq h`Ophq(3.16) Fórmula centrada para la derivada segunda: f2pαq “ fpα´hq ´ 2fpαq ` fpα`hq h2`Oph2q(3.17) Consideremos una discretización del intervalo espacial r0, Lscon Inodos equiespaciados, x1“0ăx2ă ¨ ¨ ¨ ă xI´1ăL“xI, y paso de discretización ∆x“L I´1. Consideremos, asimismo, t1“0ăt2ă ¨ ¨ ¨ ă tN´1ătN“T, discretización del intervalo espacial r0, Tscon Nnodos equiespaciados y paso de discretización ∆t“T N´1. A continuación, se presentan dos posibles esquemas en diferencias finitas para resolver el sistema. 3.1.1. Implícito estándar En este método se emplea una fórmula regresiva de primer orden para aproximar la derivada temporal (3.15) y una fórmula centrada de segundo orden para aproximar la se- 42 CAPÍTULO 3. RESOLUCIÓN DEL MODELO SIN HIPÓTESIS... gunda derivada espacial (3.17). Considérese la ecuación que rige el comportamiento del potencial de membrana. Se trata de una EDP parabólica, que puede reescribirse como: BV Bt´DB2V Bx2“Fpx, tq,(3.18) donde D“a 2RCm yF“ ´ 1 CmĚ gNam3hpV´ENaq ´ 1 CmĎ gKn4pV´EKq ´ 1 CmĎ gLpV´ELq. Así, dicha ecuación se discretiza como sigue: Vpxi, tn`1q ´ Vpxi, tnq ∆t«BV Btpxi, tn`1q “ DB2V Bx2pxi, tn`1q ` Fpxi, tn`1q « «DVpxi`1, tn`1q ´ 2Vpxi, tn `1q ` Vpxi´1, tn`1qq p∆xq2`Fpxi, tn`1q, lo que sugiere el esquema numérico: Vn`1 i´Vn i ∆t´DVn`1 i`1´2Vn`1 i`Vn`1 i´1 p∆xq2´Fn`1 i“0,1ďnďN´1, 2ďiďI´1. (3.19) Observación 11. Nótese que Vn ies la aproximación de la solución exacta en el punto pxi, tnqde la malla, esto es, Vn i«Vpxi, tnq. Asimismo, Fn idenota la aproximación del lado derecho de la ecuación, Fn i«Fpxi, tnq. Observación 12. Este esquema converge incondicionalmente con orden 1 en tiempo, Op∆tq, y orden 2 en espacio, Opp∆xq2q. i´1, n `1 i, n i`1, n `1i, n `1 Figura 3.2: Molécula computacional para el esquema implícito estándar. 3.1. MÉTODOS NUMÉRICOS DE RESOLUCIÓN DE EDP 43 Implementación del método implícito estándar A continuación veremos cómo se implementa el esquema anterior. En primer lugar, empleando la condición inicial, se tiene que V1 i“V0pxiq“´70, i“1,¨ ¨ ¨ , I. Ahora, por tratarse de un esquema implícito, para cada nP t1,¨ ¨ ¨ , N ´1ufijo, conocidos Vn 1,¨ ¨ ¨ , V n I, se obtienen las aproximaciones Vn`1 1,¨ ¨ ¨ , V n`1 Iresolviendo un sistema de ecuaciones. Para hallar las ecuaciones que conforman dicho sistema, se reordenan términos en la ecuación (3.19), obteniéndose: ˆ1 ∆t`a RCmp∆xq2`1 CmĚ gNapmn`1 iq3hn`1 i`1 CmĎ gKpnn`1 iq4`1 CmĎ gL˙Vn`1 i` ˆ´a 2RCmp∆xq2˙Vn`1 i´1`ˆ´a 2RCmp∆xq2˙Vn`1 i`1“ “1 ∆tVn i`1 CmĚ gNapmn`1 iq3hn`1 iENa `1 CmĎ gKpnn`1 iq4EK`1 CmĎ gLEL. De forma más reducida, la expresión anterior puede escribirse como: ´αV n`1 i´1`ˆ1 ∆t`2α`Φn`1 i˙Vn`1 i´αV n`1 i`1“1 ∆tVn i`Ψn`1 i,(3.20) donde α“a 2RCmp∆xq2, Φn`1 i“1 Cm“ĚgNapmn`1 iq3hn`1 i`ĎgKpnn`1 iq4`Ď gL‰, Ψn`1 i“1 Cm“ĚgNapmn`1 iq3hn`1 iENa `ĎgKpnn`1 iq4EK`ĎgLEL‰. Esta expresión solamente es válida para i“2,¨ ¨ ¨ , I ´1, pero se verá a continuación que, empleando las condiciones de contorno, se obtienen ecuaciones similares para i“1e i“I. Comenzamos consideramos la condición de contorno BV Bxp0, tq “ ´RI0 πa2y la discretizamos espacialmente empleando una fórmula centrada para la primera derivada ((3.14)). 50 CAPÍTULO 3. RESOLUCIÓN DEL MODELO SIN HIPÓTESIS... Para apreciar cómo sucede esta propagación, en la figura 3.5, obtenida tomando como referencia la figura 25 de [2], se representa el potencial de membrana en distintos puntos a lo largo del axón. Se puede observar que la forma y duración del potencial de acción permanecen invariables a lo largo del axón y, además, este viaja a una velocidad constante. En efecto, a la vista de la figura, para los datos considerados, el potencial de acción avanza a razón de 7cm cada 5ms, lo que, mediante unos sencillos cálculos, se traduce en una velocidad de propagación aproximada de 14 m{s. Figura 3.5: Propagación del potencial de acción a lo largo del axón para una corriente externa de intensidad 6µA{cm2aplicada durante el primer ms. Por otra parte, una vez estudiada la propagación, resulta de interés comprobar que, para un punto fijo del axón, el comportamiento del potencial de membrana dado por este modelo es el mismo que el proporcionado por el modelo con fijación espacial. Para ello, en la figura 3.6, se considera el potencial de membrana en el cono axónico y se observa que la amplitud y duración del potencial de acción coincide con la dada por el modelo con fijación espacial. 3.3. ANÁLISIS DE LOS RESULTADOS DE LA SIMULACIÓN NUMÉRICA 51 Figura 3.6: Potencial de acción en el cono axónico. 52 CAPÍTULO 3. RESOLUCIÓN DEL MODELO SIN HIPÓTESIS... Capítulo 4 Conclusiones y nuevas líneas de trabajo Conclusiones En este trabajo se intentó dar respuesta a la cuestión de cómo sucede la propagación del potencial de acción en una neurona, a través de la resolución del modelo propuesto por Hodgkin y Huxley. Este modelo resulta fundamental, ya que plantea el nexo existente entre los canales iónicos y la transmisión del potencial de acción, sentando las bases de todos los estudios posteriores en relación a la dinámica neuronal. Además, introduce un enfoque novedoso, otorgando igual importancia a la descripción cualitativa y cuantitativa del fenómeno, que se intentó mantener a lo largo de este trabajo. Siguiendo esta idea, se comenzó describiendo las propiedades del potencial de acción, que después fueron observadas en la simulación numérica. La simulación numérica realizada permite extraer algunas conclusiones generales en relación al modelo: A pesar de su relativa simplicidad, el modelo refleja de una manera bastante precisa la dinámica de los canales iónicos y del potencial de acción. El modelo muestra que el potencial de acción es un proceso de todo o nada, que se desarrolla solamente si el estímulo aplicado resulta suficiente para que el potencial de membrana supere el umbral de voltaje. Además, la intensidad y frecuencia de los estímulos son determinantes en su desarrollo. La simplificación realizada por Hodgkin y Huxley, al considerar que la propagación del potencial de acción ocurría de acuerdo a la ecuación de ondas, resulta ser acertada. 53 54 CAPÍTULO 4. CONCLUSIONES Y NUEVAS LÍNEAS DE TRABAJO Aunque esta simplificación parece ser fruto de las limitaciones computacionales de la época, a través de la simulación se comprobó que, en efecto, el potencial de acción se propaga en forma de una onda que viaja a velocidad constante. Por otra parte, el estudio realizado nos permitió alcanzar una mayor comprensión de los métodos numéricos en general, más allá de su aplicación particular en el modelo. Estas aportaciones pueden resumirse en los siguientes puntos: No existe un método universalmente mejor que otro en todos las situaciones y solo el conocimiento del problema a resolver puede llevarnos al método más adecuado. Además, en la elección del método, cabe tener en cuenta la precisión que se busca en la solución. Siempre resulta interesante emplear diferentes métodos para resolver un mismo problema, no solo con el objetivo de buscar el método más adecuado, sino también como una forma de validar los métodos, comparando las soluciones que proporcionan. Esto resulta especialmente relevante cuando no es posible construir soluciones manufacturadas con solución exacta no trivial que permitan validar el método. La rigidez es un fenómeno complejo y los métodos de paso adaptativo, aunque pueda esperarse lo contrario, no siempre permiten resolver problemas rígidos. En efecto, aunque estos métodos controlan el error local en cada paso, si las zonas de alto gradiente son estrechas, pasos grandes pueden dar lugar a un error pequeño, obviándose, por tanto, estas zonas de alto gradiente. El hecho de emplear métodos incondicionalmente numéricamente estables no siempre es garantía de obtener buenas aproximaciones, ya que se debe prestar también atención al condicionamiento del problema. Más en particular, se pueden extraer algunas conclusiones acerca de los métodos empleados para resolver el modelo: En el modelo con fijación espacial, las funciones intrínsecas de Matlab R ode23 y ode45 no permiten resolver el problema, ya que no detectan las zonas de alto gradiente que surgen cuando se desarrollan potenciales de acción. No obstante, limitando el tamaño máximo de paso posibilitamos que las funciones resolviesen el problema. Además, la solución proporcionada por estas no dista demasiado, en términos de precisión, de la dada por la función ode15s, que está diseñada específicamente para problemas rígidos. 55 En la propagación del potencial de acción, el método de Crank-Nicolson presenta oscilaciones espurias para pasos de discretización temporal a partir de un cierto tamaño. Este hecho determina que, aun cuando Crank-Nicolson ofrece orden 2 en espacio frente al orden 1 del esquema implícito estándar, se prefiere este último cuando se precisa utilizar pasos de discretización temporal de tamaño considerable. En caso de necesitar mucha precisión en la solución, es preferible, sin embargo, emplear Crank- Nicolson, tomando un paso de discretización temporal pequeño para evitar posibles oscilaciones. Nuevas líneas de trabajo El trabajo realizado admite un análisis más profundo. En este sentido, puede resultar de interés modificar algunos de los parámetros del modelo, como las conductancias iónicas. Así, sería posible simular, por ejemplo, cómo afectan las canalopatías (trastornos que causan disfunción de cierto tipo de canales iónicos) o ciertas enfermedades que alteran la permeabilidad iónica, a la transmisión del impulso nervioso. Por otra parte, puede llevarse a cabo un análisis similar al realizado en este trabajo para una posible generalización del modelo para axones mielinizados. En este caso, dado que la mielina actúa como un aislante eléctrico, la resistencia de la membrana aumenta con respecto al modelo original, mientras que la capacitancia disminuye; resultando en una mayor velocidad de transmisión del potencial de acción. Además, en lugar de estudiar la propagación del potencial de acción en una única neurona, el modelo puede generalizarse al estudio de una red neuronal, incorporando la sinapsis neuronal. A través de estas generalizaciones, surge un nuevo horizonte de posibles aplicaciones, especialmente en el ámbito de la medicina, como, por ejemplo: El estudio del efecto de las enfermedades desmielinizantes, tales como la esclerosis múltiple, en la transmisión el impulso nervioso. El estudio del efecto de la electroestimulación en determinadas deficiencias o trastornos, tales como la electroestimulación del nervio auditivo en personas sordas [16] o la estimulación de la corteza cerebral en personas que sufren de neuralgia del nervio trigémino [13]. El estudio de la transmisión del impulso nervioso bajo el efecto de determinadas sustancias, como alcohol o sedantes [1]. 56 CAPÍTULO 4. CONCLUSIONES Y NUEVAS LÍNEAS DE TRABAJO Por último, cabe señalar que el modelo de Hodgkin y Huxley también es aplicable, con ligeras modificaciones, al estudio de la transmisión del potencial de acción en otras células excitables, como los miocitos cardíacos. No obstante, la transmisión del potencial de acción en estas células es ligeramente más compleja, pues están implicadas más especies iónicas. Las numerosas aplicaciones e investigaciones que surgen alrededor del modelo nos dan una idea de su relevancia y vigencia. A través de este trabajo se quiso, precisamente, poner de manifiesto la importancia de este modelo que, aunque con ciertas limitaciones, especialmente en la descripción de la dinámica de los canales iónicos, consiguió dar una explicación satisfactoria en términos cuantitativos al fenómeno de la transmisión del potencial de acción. Índice de parámetros del modelo Parámetro Descripción Valor Unidades ENa potencial de equilibrio Na`45 mV EKpotencial de equilibrio K`-82 mV ELpotencial de equilibrio resto de iones -59.4 mV gNa conductancia máxima Na`120 mS{cm2 gKconductancia máxima K`36 mS{cm2 gLconductancia máxima resto de iones 0.3 mS{cm2 Cmcapacitancia de la membrana 1 µA{cm2 Vrpotencial de membrana en reposo -70 mV aradio del axón 0.0238 cm Rresistividad axial específica 35.4 Ω¨cm 57 58 ÍNDICE DE PARÁMETROS Apéndice I: Códigos de MatlabR  I.I. Modelo bajo hipótesis de homogeneidad en espacio I.I.I. Entrada de datos fprintf( ofi , ’\n∗∗∗∗∗ DATOS ∗∗∗∗∗ ’) ; %INTERVALO DE TIEMPO fprintf( ofi , ’\n∗Extremo izquierdo : iT = %´.15E. ’ , 0) T = 150; %EXTREMO DERECHO INTERVALO ( en ms) fprintf( ofi , ’\n∗Extremo derecho : T = %´.15E. ’ , T) ; %NÚMERO DE NODOS DISCRETIZACION N=15000; fprintf( ofi , ’\n\n∗Número de pasos de tiempo : N = %´i.’, N) ; %CAPACITANCIA MEMBRANA C=1; %muF/cm^2 fprintf( ofi , ’\n∗Capacitancia de la membrana : C = %´.15E. ’ , C) ; %POTENCIALES DE EQUILIBRIO ENa=45; %mV fprintf( ofi , ’\n∗Potencial e q u i l i b r i o Na: ENa = %´.15E. ’ , ENa) ; EK=´82; %mV fprintf( ofi , ’\n∗Potencial e q u i l i b r i o K: EK = %´.15E. ’ , EK) ; EL=´59.4; %mV fprintf( ofi , ’\n∗Potencial e q u i l i b r i o de fuga : EL = %´.15E. ’ , EL) ; %CONDUCTANCIAS IÓNICAS MÁXIMAS gbarNa=120; %mS/cm^2 fprintf( ofi , ’\n∗Conductancia máxima Na: gNa = %´.15E. ’ , gbarNa 59 66 APÉNDICE I legend(’m’ ,’n ’ ,’h ’ ) ; %CÁLCULO ERROR i f i s o l==1 for i =1: length (t) err ( i )=abs(V_m( i )´sol_ex ( t ( i ) , I ( t ( i ) ) ) ) ; end [ max_err , pos]=max( err ) ; fprintf( ofi , [ ’\n∗Máximo error : max_err =’ ,num2str (max_err) , ’ en la p osi c i ón ’ ,num2str ( pos ) ] ) ; figure(3) plot ( t ,V_m, ’r’) hold on plot ( t , sol_ex ( t , I ( t ) ) , ’g ’ ) t i t l e (’ Soluci ón exacta vs . aproximada ’ ) legend(’ Sol . aprox ’ ,’ Sol . exacta ’ ) hold off end rk_4.m I.I.V. Función F function z = f ( t ,X) % dX/dt=F( t ,X) % Función que calcula el valor del campo ve c t o r i al F para un instante t y %X=(Vm,m, h , n) dados . global C ENa EK EL gbarNa gbarK gbarL V_m = X(1 ,1) ; m = X(2 ,1) ; h = X(3 ,1) ; n = X(4 ,1) ; gNa=gbarNa∗m^3∗h ; gK=gbarK∗n^4; gL=gbarL ; 67 INa=gNa∗(V_m´ENa) ; IK=gK∗(V_m´EK) ; IL=gL∗(V_m´EL) ; % Lado derecho ecuaci ón d i f e r e n c i a l z = [(1/C) ∗( I ( t )´(INa+IK+IL ) ) ; a_m(V_m)∗(1´m)´b_m(V_m)∗m; a_h( V_m)∗(1´h)´b_h(V_m)∗h ; a_n(V_m)∗(1´n)´b_n(V_m)∗n ] ; end Función f.m I.I.VI. Función corriente externa function y= I ( t ) %CORRIENTE EXTERNA % Constante y=4; % Pulso constante de intensidad 4 muA/cm^2 aplicado desde el milisegundo 50 % hasta e l f i n a l . % i f t>=50 % y=4; % els e % y=0; % end end Función I.m I.I.VII. Funciones Alpha y Beta function a=a_m(V) % funci ón alpha_m a=0.1∗(V+45)./(1´exp(´(V+45)/10) ) ; 68 APÉNDICE I end Función a_m.m function b=b_m(V) % funci ón beta_m b=4.0∗exp(´1/18∗(V+70)) ; end Función b_m.m function a=a_h(V) % funci ón alpha_h a=0.07∗exp(´1/20∗(V+70)) ; end Función a_h.m function b =b_h(V) % funci ón beta_h b=1/(1+exp(´1/10∗(V+40)) ) ; end Función b_h.m function a=a_n(V) % funci ón alpha_n a=0.01∗(V+60)./(1 ´exp(´(V+60)/10) ) ; end Función a_n.m function b=b_n(V) % funci ón beta_n b=0.125∗exp(´(V+70)/80) ; end Función b_n.m 69 I.II. Modelos sin hipótesis de homogeneidad en espacio I.II.I. Entrada de datos fprintf( ofi , ’\n∗∗∗∗∗ DATOS ∗∗∗∗∗ ’) ; %INTERVALO DE TIEMPO fprintf( ofi , ’\n∗Tiempo i n i c i a l : iT = %´.15E. ’ , 0) T = 50; % ms fprintf( ofi , ’\n∗Tiempo f i n a l : T = %´.15E. ’ , T) ; %PASO DE DISCRETIZACIÓN TEMPORAL dt =0.005; fprintf( ofi , ’\n\n∗Paso de d i s c r e t i z a c i ón temporal : dt = %´i.’, dt ) ; %DOMINIO ESPACIAL fprintf( ofi , ’\n∗Extremo izquierdo dominio espacia l : a= %´.15E. ’,0) ; L=50; % cm fprintf( ofi , ’\n∗Extremo derecho dominio e spa cial : b=%´.15E. ’ ,L ) ; %PASO DE DISCRETIZACIÓN ESPACIAL dx=0.05; %CAPACITANCIA MEMBRANA C=1; %muF/cm^2 fprintf( ofi , ’\n∗Capacitancia de la membrana : C = %´.15E. ’ , C) ; %POTENCIALES DE EQUILIBRIO ENa=45; %mV fprintf( ofi , ’\n∗Potencial e q u i l i b r i o Na: ENa = %´.15E. ’ , ENa) ; EK=´82; %mV fprintf( ofi , ’\n∗Potencial e q u i l i b r i o K: EK = %´.15E. ’ , EK) ; EL=´59.4; %mV fprintf( ofi , ’\n∗Potencial e q u i l i b r i o de fuga : EL = %´.15E. ’ , EL) ; %CONDUCTANCIAS IÓNICAS MÁXIMAS gbarNa=120; %mS/cm^2 fprintf( ofi , ’\n∗Conductancia máxima Na: gNa = %´.15E. ’ , gbarNa 70 APÉNDICE I ) ; gbarK=36; %mS/cm^2 fprintf( ofi , ’\n∗Conductancia máxima K: gK = %´.15E. ’ , gbarK) ; gbarL =0.3; %mS/cm^2 fprintf( ofi , ’\n∗Conductancia máxima resto de iones : gL = %´.15 E. ’ , gbarL ) ; %RADIO AXÓN a=0.0238; % cm fprintf( ofi , ’\n∗Radio del axón : a=%´.15E. ’ ,a) ; % RESISTIVIDAD AXIAL ESPECÍFICA R=35.4e´3; % ohmiocm fprintf( ofi , ’\n∗Resistividad axi al espec í f i c a : R= %´.15E. ’ ,R) ; datos.m I.II.II. Implícito estándar % Esquema impl í cito est ándar para el modelo de Hodgkin y Huxley sin % f i j a c i ón espac ial clear a l l clo se a l l o f i=fopen(’ output4 . txt ’ ,’wt ’ ) ; datos ; % DISCRETIZACIÓN N=T/dt+1; t=0+dt ∗[ 0 :N´1]; I=L/dx+1; x=0+dx ∗[ 0 : I ´1]; %CONDICIÓN INICIAL V_0=´70.∗ones (1 , length (x) ) ; m_0=(a_m(V_0) ) /(a_m(V_0)+b_m(V_0) ) ; 71 h_0=(a_h(V_0) ) /(a_h(V_0)+b_h(V_0) ) ; n_0=(a_n(V_0) ) /(a_n(V_0)+b_n(V_0) ) ; alpha=a/(2∗R∗C∗dx^2) ; V=zeros( I ,N) ; % matriz que en la columna j´ésima contiene la s % aproximaciones de V en lo s nodos de la d i s c r e t i z a c i ón % espacial V( : , 1 )=V_0; m=zeros(length (x) , length (t)); h=zeros(length (x) , length (t)); n=zeros(length (x) , length (t)); m( : , 1 )=m_0; h ( : , 1 )=h_0; n ( : , 1 )=n_0; %CÁLCULO APROXIMACIONES for k=2: length(t) m( : , k)=dt∗((1´m( : , k´1)) .∗a_m(V( : , k´1))´m( : , k´1) .∗b_m(V( : , k´1)) ) +m( : , k´1) ; % m_i^k h ( : , k)=dt∗((1´h ( : , k´1)) .∗a_h(V( : , k´1))´h ( : , k´1) .∗b_h(V( : , k´1)) ) +h ( : , k´1) ; % h_i^k n ( : , k)=dt∗((1´n ( : , k´1)) .∗a_n(V( : , k´1))´n ( : , k´1) .∗b_n(V( : , k´1)) ) +n ( : , k´1) ; % n_i^k phi=(1/C) ∗( gbarNa∗m( : , k) .^3.∗h ( : , k)+gbarK∗n ( : , k ).^4+gbarL ) ; psi =(1/C) ∗(gbarNa∗m( : , k) .^3.∗h ( : , k) .∗ENa+gbarK∗n ( : , k) .^4.∗EK+ gbarL∗EL) ; d=sparse(1: I , 1 : I ,1/ dt+2∗alpha+phi , I , I ) ; % diagonal pr inc ipal u=sparse(1: I ´1 ,2: I ,[ ´2∗alpha ´alpha∗ones (1 , I ´2) ] , I , I ) ; % diagonal superior l=sparse (2 : I , 1 : I ´1,[´alpha∗ones (1 , I ´2) ´2∗alpha ] , I , I ) ; % diagonal i n f e r i o r A=f u l l (d+u+l ) ; b=(1/dt )∗V( : , k´1)+psi ; b(1)=b(1)+I_0( t (k) ) /(C∗dx∗pi ∗a) ; V( : , k)=A\b ; end %REPRESENTACIÓN GRÁFICA % Representaci ón 3D 72 APÉNDICE I figure(1) [X,T]=meshgrid(x , t ) ; surf (X,T,V’ , ’ edgecolor ’ ,’ none ’ ) ; % las opciones ’ edgecolor ’ ’none ’ son necesarias para que pueda verse e l mapa de colores , ya que sino , al ser la malla muy fina , l as l í neas de malla ( que son negras ) no permiten ver los col ores . t i t l e (’ Propagación del potencial de acci ón ’ ) ; xlabel (’x (cm) ’ ) ; ylabel (’ t (ms) ’ ) ; zlabel(’ Potencial de membrana (mV) ’ ) ; figure(2) x_0=find(x==8); x_1=find(x==15) ; x_2=find(x==22) ; x_3=find(x==30) ; x_4=find(x==37) ; x_5=find(x==45) ; plot ( t ,V( 1 , : ) , ’g ’ ) ; hold on plot ( t ,V(x_0 , : ) ) ; plot ( t ,V(x_1 , : ) , ’r’) ; plot ( t ,V(x_2 , : ) , ’b ’ ) ; plot ( t ,V(x_3 , : ) ) ; plot ( t ,V(x_4 , : ) ) ; plot ( t ,V(x_5 , : ) , ’y ’ ) ; legend(’x=0 ’ ,’x=8 ’ ,’x=15 ’ ,’x=22 ’ ,’x=30 ’ ,’x=37 ’ ,’x=45 ’ ) ; xlabel (’ t (ms) ’ ) ; ylabel (’V (mV) ’ ) ; hold off implicito.m I.II.III. Crank-Nicolson % Esquema Crank´Nicolson para el modelo de Hodgkin y Huxley sin 73 % f i j a c i ón espac ial clear a l l clo se a l l o f i=fopen(’ output5 . txt ’ ,’wt ’ ) ; datos ; % DISCRETIZACIÓN N=T/dt+1; t=0+dt ∗[ 0 :N´1]; I=L/dx+1; x=0+dx ∗[ 0 : I ´1]; %CONDICIÓN INICIAL V_0=´70.∗ones (1 , length (x) ) ; m_0=(a_m(V_0) ) /(a_m(V_0)+b_m(V_0) ) ; h_0=(a_h(V_0) ) /(a_h(V_0)+b_h(V_0) ) ; n_0=(a_n(V_0) ) /(a_n(V_0)+b_n(V_0) ) ; alpha=a/(2∗R∗C∗dx^2) ; V=zeros( I ,N) ; % matriz que en la columna j´ésima contiene la s % aproximaciones de V en lo s nodos de la d i s c r e t i z a c i ón % espacial V( : , 1 )=V_0; m=zeros(length (x) , length (t)); h=zeros(length (x) , length (t)); n=zeros(length (x) , length (t)); m( : , 1 )=m_0; h ( : , 1 )=h_0; n ( : , 1 )=n_0; %CÁLCULO APROXIMACIONES for k=2: length(t) m( : , k)=dt∗((1´m( : , k´1)) .∗a_m(V( : , k´1))´m( : , k´1) .∗b_m(V( : , k´1)) ) +m( : , k´1) ; % m_i^k h ( : , k)=dt∗((1´h ( : , k´1)) .∗a_h(V( : , k´1))´h ( : , k´1) .∗b_h(V( : , k´1)) ) +h ( : , k´1) ; % h_i^k 74 APÉNDICE I n ( : , k)=dt∗((1´n ( : , k´1)) .∗a_n(V( : , k´1))´n ( : , k´1) .∗b_n(V( : , k´1)) ) +n ( : , k´1) ; % n_i^k phi=(1/C) ∗( gbarNa∗m( : , k) .^3.∗h ( : , k)+gbarK∗n ( : , k ).^4+gbarL ) ; psi =(1/C) ∗(gbarNa∗m( : , k) .^3.∗h ( : , k) .∗ENa+gbarK∗n ( : , k) .^4.∗EK+ gbarL∗EL) ; f=´(1/C) .∗( gbarNa∗m( : , k´1) .^3.∗h ( : , k´1) . ∗(V( : , k´1)´ENa)+gbarK .∗ n ( : , k´1) .^4.∗(V( : , k´1)´EK)+gbarL .∗(V( : , k´1)´EL) ) ; d=sparse(1: I , 1 : I ,1/ dt+2∗alpha+1/2∗phi , I , I ) ; % diagonal principal u=sparse(1: I ´1 ,2: I ,[ ´2∗alpha ´alpha∗ones (1 , I ´2) ] , I , I ) ; % diagonal superior l=sparse (2 : I , 1 : I ´1,[´alpha∗ones (1 , I ´2) ´2∗alpha ] , I , I ) ; % diagonal i n f e r i o r A=f u l l (d+u+l ) ; b=zeros(length (x) ,1) ; b(1) =(1/dt´alpha ) .∗V(1 ,k´1)+alpha .∗V(2 ,k´1)+psi (1)+I_0( t (k) ) /(C ∗dx∗pi ∗a)+f (1) ; b (2 : I ´1)=alpha /2.∗V( 1: I ´2,k´1)+(1/dt´alpha ) .∗V(2: I ´1,k´1)+alpha /2.∗V(3: I , k´1)+1/2∗psi ( 2: I ´1)+1/2∗f ( 2: I ´1) ; b( I )=alpha .∗V( I ´1,k´1)+(1/dt´alpha ) .∗V( I , k´1)+1/2∗psi ( I )+1/2∗f ( I ) ; V( : , k)=A\b ; end %REPRESENTACIÓN GRÁFICA figure(1) [X,T]=meshgrid(x , t ) ; surf (X,T,V’ , ’ edgecolor ’ ,’ none ’ ) ; % las opciones ’ edgecolor ’ ’none ’ son necesarias para que pueda verse el mapa de colores , ya que sino , al ser una malla muy fina , l as l í neas de malla ( que son negras ) no permiten ver los col ores . t i t l e (’ Propagación del potencial de acci ón ’ ) ; xlabel (’x (cm) ’ ) ; ylabel (’ t (ms) ’ ) ; zlabel(’ Potencial de membrana (mV) ’ ) ; figure(2) 75 x_0=find(x==8); x_1=find(x==15) ; x_2=find(x==22) ; x_3=find(x==30) ; x_4=find(x==37) ; x_5=find(x==45) ; plot ( t ,V( 1 , : ) , ’g ’ ) ; hold on plot ( t ,V(x_0 , : ) ) ; plot ( t ,V(x_1 , : ) , ’r’) ; plot ( t ,V(x_2 , : ) , ’b ’ ) ; plot ( t ,V(x_3 , : ) ) ; plot ( t ,V(x_4 , : ) ) ; plot ( t ,V(x_5 , : ) , ’y ’ ) ; legend(’x=0 ’ ,’x=8 ’ ,’x=15 ’ ,’x=22 ’ ,’x=30 ’ ,’x=37 ’ ,’x=45 ’ ) ; xlabel (’ t (ms) ’ ) ; ylabel (’V (mV) ’ ) ; hold off crank_nicolson.m 82 APÉNDICE III biológico como consecuencia de un incremento de 10˝Cen la temperatura. El resto de parámetros del modelo fueron escogidos independientes a la temperatura. De esta forma, el efecto en el potencial de acción es el siguiente: dado que la velocidad de apertura y cierre de los canales aumenta con la temperatura, la despolarización es más rápida y la repolarización se inicia antes, por lo que el máximo del potencial de membrana es menor. Así, conforme aumenta la temperatura, los potenciales de acción son más rápidos y, también, más débiles. Bibliografía [1] Abiva, J.; Joseph, E.; Mikaelian, A. y Rogers, C., Alcohol’s effect on neuron firing, Repositorio de la Universidad Estatal de Arizona, 2005. [2] Amin, M., Action potential propagation in neurons [Master Thesis, Utrech University], Repositorio Institucional, 2019. [3] Amsallem, D. y Nordström, J., High-order accurate difference schemes for the Hodgkin- Huxley equations, Journal of Computational Physics, 252 (2013), 573–590. [4] Börgers, C., An Introduction to Modeling Neuronal Dynamics, 23–29, Springer International Publishing AG, 2017. [5] Boron, W. F. y Boulpaep, E. L., Fisiología médica, 3aed., 173–182, Elsevier, 2017. [6] Calderón, J. S., Análisis del modelo neuronal de Hodgkin-Huxley sobre la transmisión de potenciales de acción en las neuronas [Monografía de Pregrado, Universidad de los Andes], Repositorio Institucional, 2016. [7] Ermentrout, G. B. y Terman, D. H., Mathematical foundations of neuroscience, 11–16, Springer, New York, 2010. [8] Gerstner, W.; Kistler, W. M.; Naud, R. y Paninski, L., Neuronal Dynamics: From single neurons to networks and models of cognition, 28–57, Cambridge University Press, Cambridge, 2014. [9] Gopinathan, A. y Mathew, J., Solving the Hodgkin-Huxley equations using the compact difference scheme somadendrite, SCTIMST, Kerala, 2013. [10] Hodgkin, A. L.; Huxley, A. F. y Katz, B., Measurement of current-voltage relations in the membrane of the giant axon of Loligo, Journal of Physiology 116 (1952), 424–448. 83 84 BIBLIOGRAFÍA [11] Hodgkin, A. L. y Huxley, A. F., Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo, Journal of Physiology 116 (1952), 449–472. [12] Hodgkin, A. L. y Huxley, A. F., A quantitative description of membrane current and its application to conduction and excitation in nerve, Journal of Phisiology 117 (1952), 500–544. [13] Khodashenas, M.; Baghdadi, G. y Towhidkhah, F., A modified Hodgkin-Huxley model to show the effect of motor cortex stimulation on the trigeminal neuralgia network, Journal of Mathematical Neuroscience 9(2019). [14] Kristensen, T. y MacNearney, D., Simulating spiking neurons by Hodgkin-Huxley model [Sesión de conferencia], International Conference on Bioinformatics and Computational Biology (BIOCOMP), Las Vegas, EEUU, 22-25 julio 2013. [15] Mascagni, M. y Sherman, A., Numerical methods for neuronal modeling, Koch, C. y Segev, T. (Eds.), Methods of neuronal modeling: from ions to networks, 2aed., 569–606, MIT Press, Massachusetts, 1996. [16] Morton, J., Numerical solution of Hodgkin-Huxley’s partial differential system for ner- ve conduction [Master of Science Thesis, Oregon State University], Repositorio Institucional, 1967. [17] Motz, H. y Rattay, F., A study of the application of Hodgkin-Huxley and the Frankenhaeuser-Huxley model for electrostimulation of the acoustic nerve, Journal of Neuroscience 18 (1986), 699–712. [18] Nelson, M. E., Electrophysiological models, Koslow, S. y Subramaniam, S. (Eds.), Databasing the Brain: From Data to Knowledge, 285–301, Wiley, New York, 2004. [19] Torres, M., Modelos neuronales memresistivos [Tesis de Maestría, Instituto Potosino de Investigación Científica y Tecnológica (IPICYT)], Repositorio Institucional- IPICYT, 2016.