Download Modelos matemáticos para el estudio de la activación de la corteza

Modelos matemáticos para el estudio de la
activación de la corteza cerebral
Juan Alberto Escamilla Reyna
1999
Introducción
La modelación matemática, en el área de las neurociencias, ha tenido
éxito en los distintos niveles de organización del cerebro, baste recordar uno
de los éxitos más famosos, el logrado por los neurofisiólogos A. Hodgkin y A.
Huxley (H-H), que en el año 1952 [57], construyeron un modelo matemático
para explicar la dinámica del potencial de membrana de una neurona, en
función de las corrientes iónicas subyacentes. El modelo se toma todavı́a
como ejemplo para la elaboración de modelos “realistas” de neuronas [21],
[99], [100].
El modelo de H-H consta de cuatro ecuaciones diferenciales no lineales
–otros modelos constan de más ecuaciones–. El estudio matemático de estos
modelos es muy complejo y se ha abordado usando diversas herramientas
matemáticas : métodos numéricos, teorı́a cualitativa de ecuaciones diferenciales, teorı́a de bifurcaciones, métodos de pequeño parámetro, etc. [16], [21],
[28], [97], [100].
Otros éxitos famosos, son los logrados en el área de las redes neuronales,
desde el trabajo pionero de Mc Culloch-Pitts [80], en el cual se construye
una red neuronal, con neuronas artificiales muy simplificadas, pero basada
en algunas propiedades de las neuronas reales, hasta los trabajos de Hopfield
[59], [60]. Los trabajos de Hopfield han sido muy importantes para el estudio
de la memoria asociativa.
Actualmente, el área de las redes neuronales es muy fructı́fera, se desarrolla fundamentalmente en dos lı́neas de investigación: la primera, tiene
como objetivo construir redes neuronales que resuelvan problemas en diversas áreas, computación, matemáticas, etc., independientemente de que si
las propiedades de las neuronas y los parámetros que intervienen en la red
tienen alguna relación significativa con las neuronas reales. La segunda, intenta construir redes neuronales “realistas”, que ayuden a explicar fenómenos
fisiológicos reales. Para la primera lı́nea de investigación se pueden ver los
trabajos [10], [12], [43], para la segunda [24], [40], [41], [75].
En esta área de investigación también se presentan problemas
matemáticos complejos que se han abordado con métodos numéricos, ecuaciones diferenciales, teorı́a de bifurcaciones y técnicas de pequeño parámetro.
En el área de la electroencefalografı́a, se han elaborado modelos considerando al cerebro como un medio conductor de la electricidad, y con
una aproximación cuasi-estática, para tratar el problema de localización de
I
fuentes, para una distribución de potencial en el cuero cabelludo. Este problema no tiene solución única y usualmente se agregan hipótesis, sugeridas
por la anatomı́a y fisiologı́a del cerebro, a las posibles fuentes que originan
la distribución de potencial dada sobre el cuero cabelludo. Estos modelos
plantean problemas matemáticos que se ubican en el área matemática de los
problemas inversos. También se auxilian con los métodos matemáticos de la
teorı́a del potencial.
La principal limitante de los modelos de medio conductor cuasi-estáticos
es que no pueden informar sobre el origen fisiológico de los patrones temporales rı́tmicos observados en el electroencefalograma.
En cuanto al estudio dinámico de la actividad cortical y su relación con la
dinámica de los electroencefalogramas (EEG), se han elaborado diversos modelos. Estos se pueden clasificar, grosso modo, en modelos locales, modelos
globales y modelos locales-globales.
En los modelos locales se trata de explicar la dinámica de la actividad
cortical en función de lo que ocurre a nivel sináptico [33], [34], [75], [76], [77],
[95]. En los modelos globales, se considera que la dinámica de la actividad
cortical se puede explicar usando propiedades macroscópicas de la propia
corteza [87], [88], [89], [90], [92] y en los modelos locales-globales, se mezclan
ambos niveles [64], [66], [68], [104], [105].
Entre los modelos globales, se encuentra el elaborado por Nuñez [89],
[92]. Es un modelo global lineal de interacción neuronal que describe la
propagación de ondas de frecuencias bajas en la corteza cerebral, sus resultados coinciden –cualitativamente– con datos experimentales relacionados con
el ritmo alfa y con algunos otros fenómenos electroencelográficos [89], [92].
Sin embargo, su modelo no incluye otros fenómenos fisiológicos, como los
correspondientes a los diversos estados del sueño.
El modelo relaciona la actividad sináptica –excitatoria e inhibitoria– en
cada elemento de corteza x en el tiempo t, con el número de potenciales
de acción que se “generan” en los demás elementos de corteza en tiempos
anteriores, teniendo en cuenta las conexiones sinápticas y las velocidades de
propagación de los potenciales de acción.
El objetivo de esta tesis es elaborar y estudiar modelos que describen la
dinámica de la actividad cortical, utilizando variables y parámetros que se
puedan medir y ası́ poder correlacionarlos con las mediciones reales de dicha
actividad. Los modelos que elaboramos desarrollan las ideas de Nuñez en al
II
menos dos direcciones: se introducen nuevos parámetros y se toma en cuenta
la no homogeneidad de la corteza cerebral.
Además, el enfoque matemático usado –ecuaciones diferenciales parciales
ordinarias, técnicas de pequeño parámetro, etc.– nos permitió estudiar modelos no lineales de una manera menos complicada. Ver capı́tulos 2 y 3 de
esta tesis.
La tesis fue dividida en tres capı́tulos. En el capı́tulo 1, se hace un recuento breve sobre la estructura y funcionamiento del cerebro en sus distintos
niveles, los diversos métodos y técnicas para obtener información y comentarios sobre los diversos modelos matemáticos que existen a nivel neurona,
redes y estructuras completas. En el capı́tulo 2, describimos los fundamentos anatómicos y fisiológicos de los modelos matemáticos que construimos.
Analizamos dichos modelos y los comparamos con el modelo de Nuñez. En
el tercer capı́tulo describimos los resultados fundamentales que obtuvimos.
III
Contenido
Capı́tulo 1
El cerebro : su estructura, funciones y diferentes
formas de medir y modelar su actividad
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
El cerebro : un sistema complejo.......................................................
Algunas estructuras importantes del cerebro.....................................
La neurona y su fisiologı́a..................................................................
Bases fı́sicas del potencial de reposo de la membrana y del potencial
de acción............................................................................................
Técnicas para obtener información del cerebro..................................
Algunos modelos de estructuras del cerebro......................................
Redes neuronales................................................................................
Modelos de medio conductor..............................................................
Conclusiones.......................................................................................
1
2
4
8
12
19
22
30
34
Capı́tulo 2
Modelación global de la actividad sináptica
de la corteza cerebral
2.1 Introducción.........................................................................................
2.2 Bases fisiológicas para la modelación de la actividad sináptica en la
corteza cerebral....................................................................................
2.3 Descripción y simplificaciones en el modelo global de actividad
sináptica en una banda de corteza.......................................................
2.4 Algunos modelos alternativos al modelo de Nuñez..............................
2.5 Desarrollo y estudio de un modelo de la actividad sináptica de la
corteza cerebral....................................................................................
2.6 Transformación del sistema de ecuaciones integrales en un sistema de
ecuaciones diferenciales parciales.........................................................
2.7 La ecuación de activación para la corteza cerebral..............................
2.8 Construcción de un modelo compartimentado de activación de la
corteza cerebral....................................................................................
IV
35
36
44
51
56
60
65
68
Capı́tulo 3
Análisis matemático del modelo de activación
de la corteza cerebral
3.1 Un modelo matemático para una rebanada de corteza cerebral..........
3.2 Controlabilidad de los estados fisiológicos de la corteza cerebral........
3.3 Existencia de ondas viajeras en la corteza cerebral.............................
3.4 Ondas estacionarias en la corteza cerebral..........................................
82
89
106
126
Conclusiones ..............................................................................................
Bibliografı́a ................................................................................................
129
131
V
Capı́tulo 1
EL CEREBRO: SU ESTRUCTURA, FUNCIONES Y DIFERENTES FORMAS DE
MEDIR Y MODELAR SU ACTIVIDAD
1.1 EL CEREBRO: UN SISTEMA COMPLEJO
El cerebro es un sistema extremadamente complejo. Está compuesto
de aproximadamente 1011 neuronas densamente interconectadas. Existen
alrededor de 1015 conexiones entre las neuronas. Cada neurona se conecta
con alrededor de 104 neuronas y recibe conexiones de aproximadamente 104
neuronas. Además las neuronas están rodeadas por un número igual o mayor
de células, llamadas células gliales.
Estudiar el diagrama de conexiones del cerebro es una tarea enorme, a
pesar de que en muchas áreas del cerebro se ha logrado establecer que las
interconexiones entre las neuronas no se hacen al azar, sino que “siguen
leyes” bastante definidas. Pero aunque lográramos establecer perfectamente
el diagrama de conexiones del cerebro, éste serı́a extremadamente difı́cil de
manejar.
Tarea igual de compleja o más es determinar cómo funcionan las distintas
estructuras del cerebro y cómo se acoplan entre ellas. El cerebro procesa la
información en diferentes escalas espacio-temporales. Se dan interacciones a
nivel sináptico, entre grupos de neuronas de diversos tamaños, entre áreas
cerebrales, etc. de las cuales surgen fenómenos que no se pueden explicar en
1
una sola escala espacio-temporal y que requieren un tratamiento cuidadoso
en varias de estas escalas.
1.2 ALGUNAS ESTRUCTURAS IMPORTANTES DEL
CEREBRO
La información que se presenta en esta sección, fue obtenida de los libros
[17], [71], [106]. Para una información más detallada sobre las distintas
estructuras del cerebro y sus funciones recurrir a ellos.
Se ha logrado establecer que el cerebro se encuentra dividido en muchas
estructuras que interactúan entre sı́. De algunas de ellas se conoce aproximadamente su estructura anatómica y las funciones cerebrales en las que están
involucradas.
La corteza cerebral, la parte exterior del cerebro y, desde el punto de vista
evolutivo, la parte más reciente del cerebro, es un órgano muy complejo de
aproximadamente dos centı́metros de espesor y una superficie aproximada
de 1600 cm2 . En la cavidad craneana, la corteza se encuentra densamente
plegada, formando giros y surcos. Se divide en cuatro regiones, los lóbulos
frontal, parietal, occipital y temporal. Cada uno de ellos cumple con funciones distintas. Se considera que la corteza se encuentra involucrada en
las siguientes actividades : percepción sensorial, actividades motoras voluntarias, almacenamiento de distintos tipos de memoria, actividad conciente,
etc.
A pesar de su complejidad – gran cantidad de circunvoluciones y aproximadamente 1010 neuronas densamente interconectadas– la corteza tiene un
conjunto de caracterı́sticas anatómicas y funcionales que facilitan la elaboración de modelos que nos permiten estudiarla teóricamente.
Algunas de estas caracterı́sticas son:
1. Es un órgano estratificado que consta de seis capas horizontales de
distintas densidades, además la corteza es bastante uniforme en el plano
determinado por cada capa.
2
2. Tiene una organización columnar, es decir las neuronas se agrupan
por columnas de diversos tamaños - minicolumnas, columnas y macrocolumnas corticales - que se traslapan y se orientan perpendicularmente
a la superficie cortical.
3. La corteza cerebral está mucho más densamente interconectada consigo
misma que con otras estructuras del cerebro. La estructura que tiene
más conexiones con la corteza es el tálamo y éste contribuye con el 1%
de las conexiones que tiene la corteza.
En el siguiente capı́tulo presentaremos mayor información cualitativa y
cuantitativa que justifique los modelos teóricos que estudiaremos. Por ahora
terminaremos diciendo que en la corteza cerebral humana se procesa la información en múltiples escalas espacio temporales. Se dan interacciones entre
neuronas, entre minicolumnas, entre columnas corticales, entre macrocolumnas, entre áreas corticales e interacciones que involucran todas estas escalas.
El tálamo es una estación de relevo y centro integrador de las entradas
sensoriales a la corteza. Se cree que también regula el nivel de conciencia y
los aspectos emocionales de las sensaciones.
El hipotálamo se considera como la interfaz cerebral con los sistemas
hormonal y autónomo que controlan la homeóstasis interna del cuerpo.
El tálamo y el hipotálamo constituyen conjuntamente el diencéfalo o cerebro intermedio.
El sistema lı́mbico ocupa una porción cortical y una porción subcortical.
Se encarga de las funciones relacionadas con la emoción, la motivación y
ciertos tipos de memoria. Dos subestructuras del sistema lı́mbico son el
hipocampo y la amı́gdala. El hipocampo juega un papel importante en la
memoria. La amı́gdala coordina las acciones de los sistemas autónomo y
endocrino.
El cerebelo tiene una estructura muy sistematizada y homógenea en todo
el órgano. La función más importante del cerebelo es la coordinación de la
actividad motora y de la “postura”.
Los ganglios basales son cinco núcleos altamente interconectados: caudado, putamen, globo pálido, núcleo subtálamico y sustancia nigra. Los
3
ganglios basales se considera que están implicados en aspectos cognitivos
superiores del control motor y en otras muchas funciones.
Por último, los hemisferios cerebrales están formados por la corteza y
tres estructuras profundas: los ganglios basales, el hipocampo y el núcleo
amigdalino. Aunque cada hemisferio tiene funciones especializadas, ambos
trabajan en asocación en lo que se refiere a funciones perceptivas, cognitivas
y motoras superiores.
Existen muchas otras estructuras que no mencionamos en esta sección.
Tampoco mencionamos las distintas interrelaciones que existen entre las diversas estructuras. Para mayor información se puede recurrir a los libros
mencionados al inicio de esta sección.
1.3 LA NEURONA Y SU FISIOLOGÍA
Una exposición detallada de los temas expuestos en esta sección y la
siguiente se puede encontrar en [17], [20], [71], [106], [116], [117].
Las neuronas del cerebro presentan una gran variedad de formas y
tamaños, (Fig. #1) , pero en todas ellas se pueden identificar cuatro partes
(Fig. #2):
1. El cuerpo celular o soma.
2. Las dendritas.
3. El axón.
4. Los terminales axónicas o sinápticas.
En el funcionamiento de la neurona, cada una de estas partes tiene un
papel especı́fico:
1. El cuerpo celular es el centro metabólico de la neurona y el lugar donde
se procesa la información. Su diámetro varı́a de 5 a 100µ.
4
2. Las dendritas son el centro receptor de la información de la neurona.
Por lo general, la neurona recibe la información (en su mayor parte de
tipo eléctrico) a través de las dendritas, aunque algunas veces la recibe
directamente a través del cuerpo celular.
3. El axón es una fibra nerviosa que nace del cuerpo celular y es el encargado de transportar las señales que genera la neurona, en ocasiones
transportando la señal a distancias considerables. El diámetro de los
axones en humano, varı́a de 0.1µ a 20µ. Frecuentemente los axones
están envueltos en una vaina grasa, denominada mielina. La vaina está
dividida por espacios regulares sin mielina, llamados nodos de Ranvier. La misión de la vaina de mielina es aumentar la velocidad de
propagación de las señales generadas por la neurona.
4. Los terminales axónicos o sinápticos constituyen los elementos de transmisión de la neurona. A través de ellos una neurona contacta y transmite la información a la zona receptora de otra neurona. La zona
de contacto se llama sinápsis y por lo general la zona postsináptica se
ubica en las dendritas, aunque también ocurren contactos sinápticos en
el cuerpo celular de la neurona y algunas veces en el axón de la neurona.
Existen aproximadamente 104 contactos sinápticos por neurona.
En condiciones de reposo, toda la neurona (cuerpo celular, soma, dendritas) está polarizada de manera tal que el interior tiene un potencial de
-70 mV , considerando que la superficie externa de la membrana es el nivel
de referencia, es decir, tiene potencial cero. Este potencial de la neurona es
llamado potencial de reposo de la membrana. El potencial de reposo puede
variar, dependiendo del tipo de neurona, pero el interior de la membrana de
la neurona queda siempre cargado negativamente, con respecto al exterior
de la misma. Además este potencial puede cambiar en respuesta a diversos
estı́mulos como temperatura, PH, concentraciones extracelulares de iones,
corrientes eléctricas, etc.
Como ya mencionamos anteriormente la neurona recibe la información
(de otras neuronas o información producida por medios artificiales) a través
de las dendritas. Esta información se presenta como un cambio en el potencial eléctrico en las dendritas; que llamaremos potenciales postsinápticos.
Estos potenciales, tienen una magnitud entre 0.1 − 10 mV , son de naturaleza
5
local, pasivos y disminuyen en intensidad progresivamente y no se detectan
más allá de 1 ó 2 mm del sitio de origen. El cuerpo celular los suma espacial
y temporalmente y cuando se alcanza un cierto umbral, la neurona “genera”
una señal, llamada potencial de acción. Esta señal es un impulso eléctrico
que se propaga por el axón hasta los terminales axónicos, donde la despolarización produce la liberación de neurotransmisores que se propagan a través
del espacio sináptico y afectan a la membrana postsináptica. Esta afectación
modifica el estado eléctrico de la membrana y puede ser de tipo excitatorio
despolarizándola y ası́ ayudar a la neurona receptora a que “genere” más
fácilmente un potencial de acción, o de tipo inhibitorio que tiene el efecto
contrario. Estos potenciales postsinápticos serán llamados potenciales excitatorios postsinápticos (PEPS) si excitan a la neurona y potenciales inhibitorios
postsinápticos (PIPS) si la inhiben.
La señal que genera una neurona, es decir, el potencial de acción tiene las
siguientes propiedades :
1. Se propaga activamente a lo largo del axón.
2. No disminuye su intensidad en función de la distancia.
3. Es de naturaleza todo o nada.
4. Es semejante en todas las neuronas. La amplitud del potencial de
acción es de unos 100 mV y dura de 0.5 a 2 milisegundos. Su velocidad
de propagación depende del diámetro del axón y de su mielinización.
5. La mayorı́a de las neuronas generan un potencial de acción cuando el
estı́mulo aplicado es despolarizante. La magnitud del estı́mulo debe
ser mayor que un cierto valor crı́tico, llamado el umbral de disparo de
la neurona. Este umbral depende de la neurona y aún en la misma
neurona puede cambiar.
También existen neuronas que generan un potencial de acción cuando
el estı́mulo aplicado es hiperpolarizante y de cierta duración.
6. Si el mismo estı́mulo aplicado a la neurona, se mantiene durante un
cierto tiempo, el umbral de disparo de la neurona disminuye.
6
7. Una vez que se genera un potencial de acción, la neurona es incapaz,
durante cierto tiempo, de generar un nuevo potencial de acción, independientemente de la magnitud del estı́mulo que se le aplique, en este
caso se dice que la neurona está en estado refractario. Al tiempo que
dura la neurona en estado refractario se le llama perı́odo refractario
absoluto. Existe un intervalo de tiempo, llamado perı́odo refractario
relativo, durante el cual el umbral de disparo de la neurona aumenta y
con esto dificultar la generación de nuevos potenciales de acción.
El perı́odo refractario absoluto de una neurona determina la frecuencia
de generación de potenciales de acción.
8. Un potencial de acción se puede generar mediante estı́mulos de corriente
que varı́en linealmente con el tiempo, con tal que la pendiente sea mayor
que un cierto valor crı́tico. Si la pendiente es menor que este valor
crı́tico, no se genera ningún potencial de acción, sin importar qué tan
grande sea la despolarización de la membrana. A esta propiedad de la
membrana se le llama acomodación.
El mecanismo de propagación del potencial de acción se puede describir
esquemáticamente de la siguiente manera: el estı́mulo inicial en un punto
del axón (una despolarización de la neurona de cierta magnitud) provoca
corrientes locales que fluyen de manera pasiva a través de la membrana,
causando una variación en el potencial de reposo en zonas cercanas al punto
estimulado.
Esta despolarización en las zonas cercanas “genera” un mecanismo en la
membrana que genera un voltaje muchas veces mayor que el inicial, el cual
a su vez produce corrientes locales que causan un cambio en el potencial en
las zonas adyacentes y ası́ sucesivamente.
Está claramente demostrado que el potencial de acción no se puede propagar por el axón como si éste fuera un cable (en forma pasiva) ya que en este
caso el potencial decrecerı́a en forma exponencial y la rapidez de decrecimiento estarı́a en función de la resistencia de la membrana axónica y de la
resistencia a las corrientes que fluyen longitudinalmente a través del axón
7
1.4 BASES FÍSICAS DEL POTENCIAL DE REPOSO DE LA
MEMBRANA Y DEL POTENCIAL DE ACCIÓN
Toda neurona contiene en su interior varios tipos de iones, entre los que
se encuentra el sodio (N a+ ), el potasio (K + ) y el cloro (Cl− ) en distintas
concentraciones, además de algunas otras moléculas como aminoácidos, proteinas, etc. La membrana de una neurona tiene distintas permeabilidades con
respecto a cada tipo de ion. La difusión de los iones a través de la membrana
celular se realiza a través de zonas especiales restringidas de la membrana,
llamados “poros” o “canales” y existe selectividad de estos canales para los
distintos tipos de iones. La membrana tiene también mecanismos activos
para “bombear” algunos de los tipos de iones (principalmente el potasio y el
sodio) del interior al exterior de la membrana y recı́procamente.
En el exterior de la neurona también se encuentran distintos tipos de
iones, principalmente sodio, potasio y cloro, en distintas concentraciones.
Para el caso de la membrana neuronal en reposo del axón gigante del
calamar, (ver Tabla 1). Los tipos de iones y las concentraciones de éstos,
pueden variar, dependiendo del tipo de neurona que se considere.
Tabla 1
Lı́quido
Potencial de
Especie
Citoplasma extracelular equilibrio∗
iónica
(mM)
(mM)
(mV)
+
K
400
20
−75
N a+
50
440
+55
Cl−
52
560
−60
A− (iones 385
–
–
orgánicos)
La distribución de los principales iones a través de la membrana
neuronal en reposo del axón gigante del calamar.
∗
Se denomina ası́ al potencial de membrana en el que no hay flujo neto de iones a través
de la membrana
8
Actualmente se considera que el potencial de reposo de la membrana se
logra cuando los flujos pasivos y activos de N a+ , K + , Cl− están balanceados
(Fig. #3), lo cual se logra cuando las fuerzas de difusión, provocadas por las
distintas concentraciones de los iones, en el interior y exterior de la membrana
se equilibran con las fuerzas eléctricas, producidas por el hecho de que los
iones son partı́culas cargadas eléctricamente.
El potencial de acción se logra cuando la neurona se despolariza hasta
un nivel crı́tico; en ese momento tiene lugar un cambio súbito, donde las
barreras normales que se oponen al flujo de sodio y potasio desaparecen
momentáaneamente, estableciéndose un flujo local de iones, de suficiente
magnitud como para invertir el potencial de membrana, que alcanza unos
50 mV (positivo en el interior) y después se invierte de nuevo para establecer
el potencial de reposo. Pero la primera inversión (hacia el interior positivo)
ha producido una señal potente que se extiende y lleva la región adyacente de
la membrana a un nivel crı́tico, repitiéndose entonces el proceso (Fig. #3).
En la siguiente sección explicaremos más detalladamente esto y describiremos
los diversos modelos que existen al respecto.
Una ecuación importante en la fisiologı́a es la ecuación de Nernst, obtenida
en 1888 por W. Nernst. Esta ecuación relaciona el potencial de una membrana permeable a un ion con las concentraciones del ion en el interior y el
exterior de la membrana. La ecuación se aplica con la condición de que la
corriente total iónica sea igual a cero, es decir, cuando la corriente iónica
producida por el potencial es igual a la corriente producida por la diferencia
de concentraciones en el interior y el exterior de la neurona.
Supongamos que tenemos dos compartimientos A y B, donde se tiene una
solución que contiene un ion I, a diferentes concentraciones CA y CB y que
los compartimientos están separados por una membrana permeable al ion I
(Fig. #4)
Figura 4
A
ion I
CA +
ion I
−CB
V
9
B
Como el ion I se encuentra a diferentes concentraciones en A y B, entonces
se empezará a difundir del compartimiento de mayor concentración al de
menor concentración pero como el ion lleva una carga, el compartimiento
hacia el cual se está difundiendo empieza a aumentar su carga ası́ aumenta su
voltaje frenando la difusión del ion. El equilibrio dinámico se logrará cuando
el voltaje en el compartimiento que recibe los iones difundidos se equilibra
con la fuerza de difusión que surge debido a las distintas concentraciones del
ion. Es decir, cuando:
µ
¶
RT
CB
(1.1)
V =
ln
ZF
CA
donde V es el voltaje entre el compartimiento A y el compartimiento B, R
es la constante de los gases ideales, T la temperatura (en grados Kelvin), Z
la valencia del ion y F la constante de Faraday (la carga de una mole de
partı́culas univalentes), aproximadamente 9.65 × 104 C/mol. Esta ecuación
nos permite calcular el voltaje entre A y B en términos de las concentraciones
del ion en los compartimientos A y B, cuando el flujo neto de iones entre los
compartimientos es cero. Este voltaje entre A y B es llamado el potencial
de equilibrio del ion. En la Tabla 1 están los valores de los potenciales de
equilibrio de algunos iones.
En el caso que los compartimientos contengan más de un ion, con distintas concentraciones y que la membrana contenga canales independientes
para el paso de los iones, se puede calcular el potencial de equilibrio de un
solo ion, usando (1.1), con la condición de que la membrana tenga abierto,
únicamente, el canal de ese ion. Si todos los canales de la membrana están
abiertos simultáneamente, entonces el voltaje dependerá de las permeabilidades (facilidad con que un ion atraviesa la membrana (cm/seg)) relativas
de los iones y de las concentraciones de los iones. Goldman (1943) calculó el
potencial final promedio:
ÃP
P
RT
P + [I + ]A + Pi− [Ii− ]B
ln P i+ +i
V =
P
Zi F
Pi [Ii ]B + Pi− [Ci− ]A
!
(1.2)
donde Pi+ , Pi− representan las permeabilidades de las especies iónicas positivas y negativas respectivamente y [Ii+ ]∗ , [Ii− ]∗ son las concentraciones de estos
iones en el compartimiento ∗.
10
En el caso de la membrana de una neurona sabemos que en el interior
y el exterior de la misma se encuentran varios tipos de iones con distintas
concentraciones y permeabilidades. Entre los principales se encuentran el
sodio (N a+ ), el potasio (K + ) y el cloro (Cl− ). Aplicando (1.2) para estos
tres iones obtenemos:
"
KT
Pk+ [K + ]e + PN a+ [N a+ ]e + PCla− [Cl− ]i
V =
ln
F
Pk+ [K + ]i + PN a+ [N a+ ]i + PCla− [Cl− ]e
#
(1.3)
En la ecuación de Goldman se supone que el campo eléctrico que existe
a través de la membrana varı́a uniformemente, por lo que también se supone
que los movimientos de los iones a través de la membrana son independientes
unos de otros.
La importancia de este resultado es que nos permite hacer predicciones
cuantitativas del potencial de membrana cuando se tienen variaciones de las
concentraciones de los iones y de las permeabilidades de la membrana a los
distintos iones.
Usando la ecuación de Nernst y con ciertos experimentos de tipo eléctrico
se ha logrado demostrar que la membrana de una neurona en estado de reposo
no es permeable únicamente al potasio sino que también es permeable al
sodio, aunque la permeabilidad del sodio es muy pequeña con respecto a la
del potasio. Se considera que la permeabilidad con respecto a los iones de
sodio es aproximadamente el 1% de la correspondiente a los iones de potasio
y por ello debe de haber una “bomba” en la membrana que bombea iones
de sodio hacia el exterior de la membrana e iones de potasio al interior de la
misma.
También se logró establecer lo que ocurre cuando una neurona genera
un potencial de acción. El potencial de acción se genera cuando se logra
una despolarización de la membrana de aproximadamente 15 mV y en ese
momento ocurre que la permeabilidad de la membrana al sodio aumenta, por
la apertura de canales de N a+ (que dependen del voltaje), lo que hace que
el potencial se eleve. Esto a su vez permite que aumente la permeabilidad
de la membrana, por un nuevo aumento de los “canales” de sodio, lo que
hace que el potencial se eleve y ası́ sucesivamente, hasta que el potencial casi
alcanza el potencial de equilibrio del sodio (se establece en la membrana un
11
mecanismo de retroalimentación positiva entre el potencial de membrana y
la permeabilidad de ésta al sodio). Este potencial de equilibrio no se alcanza
ya que cuando se inicia el aumento de la permeabilidad de la membrana al
sodio, un corto tiempo después empieza a aumentar la permeabilidad de la
membrana al potasio y una entrada de cloro a la membrana. Además los
canales de sodio empiezan un proceso de inexcitabilidad y también se tiene
un flujo pasivo de iones de potasio.
En resumen, el potencial de acción se origina cuando la estimulación
eléctrica alcanza un cierto umbral. Su dinámica se encuentra determinada
por la relación que existe entre las distintas corrientes iónicas que se establecen en la membrana, tanto corrientes pasivas como activas.
En la sección 1.6 se mostrarán algunos modelos que describen este proceso
de manera cuantitativa.
Una manera equivalente de estudiar el potencial de membrana de una
neurona consiste en construir un circuito eléctrico que incorpore las tres
propiedades básicas que utiliza la membrana para generar señales eléctricas
- los canales iónicos, los gradientes de concentración de los iones, la capacidad de la membrana para almacenar carga eléctrica. Cada canal iónico se
puede representar como una resistencia (o como una conductancia que es
la inversa de la resistencia). La conductancia para un canal iónico serı́a en
cierta medida, equivalente a la permeabilidad de la membrana a ese ion. La
capacidad de la membrana para almacenar carga se puede representar como
una capacitancia y los gradientes de concentración como una baterı́a (Figs.
#5,6)
El modelo de circuito equivalente es muy importante desde el punto de
vista experimental ya que permite el uso de técnicas electrofisiológicas para
el estudio de las señales neurales, en lugar del uso de isótopos radiactivos que
se usan para medir los cambios de permeabilidad y difusión de iones.
1.5 TÉCNICAS PARA OBTENER INFORMACIÓN DEL
CEREBRO
En el proceso de obtener información sobre la estructura anatómica y
funcional del cerebro en sus diversos niveles –neurona, poblaciones de neu12
ronas y todo el cerebro– se han elaborado, con la ayuda de la tecnologı́a y
otras ciencias, un conjunto de técnicas e instrumentos que se han ido mejorando con el transcurso del tiempo y que han sido de gran utilidad para el
diagnóstico clı́nico, para los neuroanatomistas y los neurofisiólogos.
Para el estudio anatómico se han usado instrumentos como el microscopio
óptico, el microscopio electrónico, etc. y algunos métodos de tinción selectiva
como el método de Golgi, Nissl, Nauta, etc. El método de Golgi, inventado
por Camilo Golgi alrededor de 1875, es un método de tinción selectiva que
tiene la propiedad de que al teñir un grupo de neuronas, sólo algunas de
ellas, aproximadamente una de cada cien se tiñen por completo (dendritas,
cuerpo celular, axón) y las demás quedan intactas. Con este método, Santiago Ramón y Cajal logró establecer que el cerebro está constituido por
células separadas, que las interconexiones entre ellas no eran al azar, además
de describir la arquitectura de una gran cantidad de estructuras del cerebro.
Con el método de Nissl se logra destacar los cuerpos celulares de las neuronas,
pero no sus dendritas y axones. Por mucho tiempo los neuroanatomistas usaron estos dos métodos de tinción para estudiar anatómicamente el cerebro.
Alrededor de 1952 W. Nauta inventó un método de tinción selectiva, que
actualmente lleva su nombre. Según D. H. Hubel, este método es el primer
instrumento eficaz para cartografiar las conexiones entre las distintas estructuras del cerebro. El método de Nauta aprovecha el hecho de que, cuando
se destruye una neurona, la fibra nerviosa que procede de ella degenera y,
antes de desaparecer por completo, puede ser teñida de manera diferente de
sus vecinos normales.
Los avances tecnológicos y una mayor comprensión del funcionamiento
de las neuronas han permitido crear técnicas con las cuales, no sólo se logra
obtener información más precisa sobre las conexiones entre regiones cerebrales, identificar neuronas, sino también sobre las zonas activas del cerebro.
Existen técnicas que aprovechan el hecho de que las neuronas absorben
y sintetizan varias sustancias – proteı́nas, enzimas, etc. –. Algunas de éstas
sustancias son transportadas, a través del axón, desde el soma hasta los terminales axónicos (transporte anterógrado). También existe un transporte
de sustancias desde los terminales axónicos al cuerpo celular (transporte
retrógrado).
13
La autorradiografı́a de Cowan (1972) es una técnica que aprovecha el
transporte axónico anterógrado para determinar los terminales axónicos de
una neurona. Esta técnica consiste en inyectar, en alguna estructura cerebral, aminoácidos marcados radiactivamente, éstos penetran en el soma de
la neurona y se incorporan en las sustancias que son transportadas hasta los
terminales axónicos, donde se les puede detectar por un proceso de autorradiografı́a.
También se cuenta con técnicas que aprovechan el transporte axónico
retrógrado. Entre las más conocidas esta la técnica de la peroxidasa de
rábano, se inyecta en el sistema nervioso, la enzima peroxidasa de rábano,
ésta es absorbida por los terminales axónicos y transportada retrógradamente
hacia el cuerpo celular, donde puede ser detectada por diversas técnicas, y
ası́ conocer el origen del axón.
Para mayor información de otras sustancias que se usan para el estudio de
las conexiones entre regiones cerebrales, otros métodos para detectar sustancias (la inmunohistoquı́mica se basa en la detección de sustancias tisulares
por medio de anticuerpos), la hibridación in situ, etc. recurrir a [51] y [107].
En este último se hace una revisión detallada de la organización neuronal y
sináptica de la corteza cerebral.
Basándose en el hecho de que la glucosa es el combustible de las neuronas,
y éstas consumen más glucosa cuando están activas que cuando se encuentran
en reposo, Louis Sokoloff (1977) inventó una técnica que le permitió medir
el metabolismo de las neuronas.
La técnica consiste en inyectar en el torrente sanguı́neo un análogo
quı́mico de la glucosa (la 2–desoxiglucosa), que es absorbida por las neuronas,
mediante el mismo sistema que la glucosa, aunque una vez en el interior de
la neurona, ya no puede ser metabolizada. Si la 2–desoxiglucosa es marcada
radiactivamente e inyectada en el torrente sanguı́neo, la velocidad con que
se acumula en el interior de las neuronas, nos proporciona un ı́ndice de la
actividad metabólica de las neuronas.
Las técnicas que hemos mencionado hasta ahora requieren en algún momento de la destrucción del sujeto de estudio. Con el desarrollo de la tecnologı́a se crearon nuevas técnicas menos invasivas, que permitı́an obtener
información sobre la actividad del cerebro humano en vivo, por ejemplo la to14
mografı́a axial computarizada (TAC), la tomografı́a por emisión de positrones
(TEP) y la formación de imagenes por resonancia magnética (IRM), entre
otras. Las herramientas teóricas fundamentales para el análisis de los diferentes tipos de tomografı́as están vinculadas a la transformada de Radon y
los métodos matemáticos y computacionales de reconstrucción por proyecciones. Existe una vasta bibliografı́a de estos temas entre los que vale citar
a [53], [54] y [85].
La tomografı́a axial computarizada se aprovecha del hecho que los distintos tejidos absorben diferentes cantidades de energı́a de rayos X; mientras
más denso es el tejido, tanto más absorbe. Por ello un haz de rayos X muy
concentrado que atraviese el cuerpo emergerá del mismo a un nivel reducido
de intensidad que dependerá de los tejidos y órganos encontrados a su paso.
Si se dirigen a través del cuerpo haces planos de rayos X con diferentes
ángulos de inclinación, se recogerá información suficiente para reconstruir
una imagen de la sección del cuerpo, la cual debidamente procesada permite
reconstruir una imagen tridimensional de la intensidad de radiación absorbida
y por ende de la estructura interna del objeto estudiado. La tomografı́a por
emisión de positrones, permite cartografiar las estructuras cerebrales activas
ya que permite la medición de diversas actividades del cerebro, tales como
el metabolismo de la glucosa, el consumo de oxı́geno, el riego sanguı́neo y la
interacción con los fármacos. De todas estas variables, el flujo sanguı́neo ha
demostrado ser el indicador más fiable de la función cerebral en cada momento. Con la TEP se pueden localizar cambios de actividad con precisión
milimétrica.
La formación de imágenes por resonancia magnética (IRM) es un instrumento bastante común para diagnosticar lesiones en tejidos. Tuvo su
origen en una técnica de laboratorio denominada resonancia magnética nuclear (RMN), diseñada con el fin de explorar las caracterı́sticas quı́micas de
las moléculas. La (RMN) se basa en que muchos átomos se comportan como
diminutas agujas imantadas en presencia de un campo magnético. En estas condiciones, la aplicación de impulsos de una onda radioeléctrica a una
muestra, perturba los átomos de una manera muy precisa, y éstos en consecuencia emiten señales de radio detectables que identifican el número y
el estado de los átomos pertenecientes a la muestra. El ajuste del campo
magnético y de los impulsos de onda radioeléctrica proporciona información
15
sobre la muestra de estudio. En el caso del cerebro P. C. Lauterbur descubrió
que podı́an formarse imagenes por RMN mediante la detección de protones.
La utilidad de estas partı́culas, abundantes en el cuerpo humano, reside en
su gran sensibilidad a los campos magnéticos que los hace comportarse como
minúsculas agujas imantadas. De esta forma se obtuvieron excelentes imagenes anatómicas con detalle muy superior a las conseguidas mediante TAC.
El interés por aplicar la IRM a la obtención de imágenes cerebrales se
debe a que esta técnica es capaz de detectar una señal inaccesible a las exploraciones mediante TEP como por ejemplo, el aumento en la concentración
de oxı́geno que se produce en una zona donde se ha intensificado la actividad
neuronal. Otras ventajas sobre la TAC y otras técnicas son:
1. La señal procede directamente de cambios inducidos por causas funcionales en el tejido cerebral (el cambio en la concentración de oxı́geno
venenoso).
2. Para obtener una señal no hay que inyectar ninguna substancia, radiactiva o de otro tipo.
3. La IRM proporciona información anatómica y funcional sobre cada
sujeto, lo cual permite una identificación estructural exacta de las zonas
activas.
4. Adecuadamente equipada, puede observar en tiempo real el ritmo de
cambio de la señal de oxı́geno inducida en el flujo sanguı́neo.
Sin embargo a pesar de las ventajas que proporciona la combinación de
estas técnicas, ninguna de ellas puede determinar la evolución temporal del
intercambio de información en diferentes zonas del cerebro, lo cual es fundamental para estudiar cómo se coordinan en red diferentes zonas del cerebro
para la realización de una actividad dada. El obstáculo que enfrentan la
IRM y TEP para este objetivo es la gran diferencia de velocidades entre la
actividad neuronal y el ritmo de variación de los niveles de oxigenación en las
neuronas. Las señales viajan de una parte a otra del cerebro en 0.01 segundos
o menos, pero los cambios de flujo sanguı́neo y oxigenación de la sangre son
mucho más lentos y ocurren desde décimas de segundo hasta varios segundos más tarde. La IRM no puede seguir las “conversaciones” entre regiones
16
cerebrales. Los únicos métodos que responden con rápidez suficiente son las
técnicas de registro eléctrico, es decir, la electroencefalografı́a (EEG) que
detecta la actividad eléctrica del cerebro a partir de mediciones efectuadas
sobre el cuero cabelludo y la magnetoencefalografı́a (MEG), la cual mide el
comportamiento electromagnético que genera la actividad eléctrica en el interior del cerebro. El defecto de estas técnicas es que su resolución espacial
es menor que la de otras técnicas y la localización precisa de las fuentes de
actividad cerebral sigue siendo un problema muy complejo.
En 1929 el psiquiatra alemán Hans Berg demostró que era posible medir
–usando electrodos colocados en el cuero cabelludo– la actividad eléctrica
del cerebro humano. Demostró que bajo ciertas condiciones en el sujeto de
estudio (ojos cerrados, sin estı́mulos exteriores, relajado) el registro de la
actividad eléctrica seguı́a un patrón temporal rı́tmico al que llamó ritmo α.
Este patrón temporal tenı́a una frecuencia entre 9 y 10 c/s. Los orı́genes fisiológicos del ritmo α, siguen siendo motivo de controversia, algunos fisiólogos
consideran que el ritmo α se origina en el tálamo y que el tálamo le impone
este ritmo a la corteza, a través de las conexiones que existen entre ambos, mientras que otros consideran que este ritmo es una manifestación de la
actividad masiva de la corteza cerebral.
Posteriormente se han descubierto otros patrones temporales rı́tmicos que
se producen bajo ciertas condiciones corporales especiales. Por ejemplo, el
ritmo δ (0 a 5 Hz) se produce bajo sueño profundo, el ritmo β (15 a 30 Hz)
se produce bajo intensa actividad mental y otros más.
El EEG (medido en el cuero cabelludo o en otras regiones del cerebro)
se ha convertido en una herramienta clı́nica importante para el estudio de
ciertas patologı́as del cerebro, como la epilepsia, tumores cerebrales, etc.
Estas patologı́as aparecen en el EEG como patrones anormales de la actividad
eléctrica espontánea del cerebro. Los potenciales evocados (los potenciales
que se producen como repuesta a un estı́mulo) se han usado en gran medida
para el diagnóstico clı́nico de las afecciones sensoriales (auditivos, visuales,
etc.) y en la detección de problemas cognitivos [14], [92], [105].
Uno de los problemas fundamentales de la electroencefalografı́a consiste
en correlacionar los patrones espacio-temporales obtenidos mediante el EEG
con la actividad eléctrica del cerebro y con conductas observables. Si se quiere
17
elaborar modelos que describan la dinámica de la actividad eléctrica del cerebro y que se correlacionen con los patrones espacio-temporales obtenidos con
el EEG; se requiere tomar en cuenta las escalas espacio-temporales en que
están promediando las mediciones efectuadas ya que los datos electroencefalográficos son obtenidos promediando la actividad neuronal en distintas
escalas espacio-temporales. Existen EEG’s con distintas resoluciones tanto
espaciales como temporales. Por ello las variables a considerar en nuestros
modelos deben estar en congruencia con las escalas espacio-temporales usadas en el EEG. En el siguiente capı́tulo detallaremos los problemas que se
presentan en esta área de estudio.
Los modelos que más han sido estudiados son los que utilizan los principios de la teorı́a del campo electromagnético considerando al cerebro como
un medio continuo conductor, con una cierta geometrı́a (que depende de las
hipótesis del modelo) y cierta conductividad variable.
En el problema directo, se considera dada la ubicación de las fuentes
de corrientes impresas (zonas activas), su configuración, la geometrı́a de la
cabeza y las conductividades de los distintos comportamientos (cuero cabelludo, cráneo, corteza, lı́quido extracraneal, etc.) y se calcula la distribución
del potencial generado en estas condiciones, comparándose los resultados
teóricos con los obtenidos con técnicas electroencefalográficas. Mediante el
aporte de ambos datos se puede encontrar alguna información sobre la respuesta electrica de una zona activa dada.
El problema inverso consiste en obtener información sobre la ubicación de
las fuentes de corrientes impresas y las conductividades en diferentes partes
del cerebro usando los datos discretos que se obtienen con el EEG. Este
problema no tiene solución única, es decir, existe una multitud de distribuciones de fuentes de corriente en el cerebro que pueden producir la misma
distribución de potencial en el cuero cabelludo. Al no tener solución única
el problema inverso, usualmente se agregan restricciones, sugeridas por la
anatomı́a y fisiologı́a del cerebro, a las posibles fuentes que originan la distribución de potencial dada sobre el cuero cabelludo.
Usualmente estos modelos se restringen a la aproximación electrostática
por lo que son útiles para estudiar la distribución espacial de generadores de
los EEG pero no ası́ sus propiedades temporales [9], [105].
18
Para estudiar la dinámica espacio-temporal de la actividad eléctrica se
pueden construir ecuaciones plausibles, a partir de datos anatómicos y fisiológicos que involucren variables que sean congruentes con las escalas
espacio-temporales usadas en los distintos EEG, para estudiar sus soluciones
y su relación con las mediciones obtenidas a través del EEG.
Finalmente, otro método usado consiste en elaborar una jerarquı́a de
escalas espacio-temporales. En cada escala espacio-temporal se deducen las
ecuaciones que gobiernan las variables que se consideran idóneas en ese nivel
y se estudia la relación que existe entre un nivel y el nivel superior inmediato.
En la sección siguiente desarrollaremos con mayor detalle algunos de estos
modelos.
1.6 ALGUNOS MODELOS DE ESTRUCTURAS DEL
CEREBRO
En el estudio del cerebro se han elaborado diversos modelos para los
distintos niveles del mismo. Se cuenta con modelos a nivel de neurona, grupos
de neuronas, de subestructuras del cerebro y del cerebro completo. Cada uno
de estos modelos toma en cuenta ciertos aspectos del cerebro dependiendo
de los problemas que se quieren estudiar. Es importante recalcar que los
elementos y las relaciones que se toman en cuenta para la elaboración del
modelo están intimamente relacionados con el nivel del cerebro que se esté
considerando y con los fenómenos que se quieren explicar.
Uno de los primeros éxitos importantes en la modelación matemática en
el estudio del cerebro fue el logrado por A. Hodking y A. Huxley (H.H) en
1952 [57].
El modelo por ellos propuesto describe, por medio de un sistema de ecuaciones diferenciales no lineales, cómo varı́a temporalmente el potencial de
membrana en función de las corrientes iónicas, que fluyen a través de la membrana en el axón de un calamar. El modelo toma en cuenta tres corrientes
iónicas, la de potasio (IK ), la de sodio (IN a ) y una corriente de filtraje (IL ).
Considera que las corrientes de sodio y potasio fluyen a través de canales
cuya conductancia depende del potencial de membrana y del tiempo.
19
De acuerdo con las leyes de Kirchhoff, se tiene que:
Cm
dV
= IN a + IK + IL ,
dt
(1.4)
donde Cm es la capacitancia de la membrana. También se tiene:
IN a = gN a (V − VN a ) ,
(1.5)
IK = gK (V − VK ) ,
(1.6)
IL = gL (V − VL ) ,
(1.7)
donde g∗ es la conductancia del ion ∗ y V∗ el potencial de equilibrio de dicho
ion. Luego se tiene que
Cm
dV
= gN a (V − VN a ) + gK (V − VK ) + gL (V − VL ) .
dt
(1.8)
En este modelo se considera que gN a y gK son dependientes del potencial
V y gL es constante, es decir:
Cm
dV
= m3 hḡN a (V − VN a ) + n4 ḡK (V − VK ) + gL (V − VL ) ,
dt
(1.9)
donde ḡN a , ḡK son las conductancias máximas del sodio y el potasio y 0 <
m, n, h < 1 donde m y n son variables de activación y h una variable de
inactivación que obedecen las ecuaciones:
dm
= αm (1 − m) − βm m ,
dt
(1.10)
dn
= αn (1 − n) − βn n ,
dt
(1.11)
dh
= αh (1 − h) − βh h ,
dt
(1.12)
con coeficientes αm , βm , αn , βn , αh , βh que son funciones obtenidas empirı́camente
y que dependen del potencial V .
20
Hodking y Huxley describen un gran número de experimentos controlados que involucran corrientes, conductancias y potenciales, cuyos resultados
“coinciden” con las predicciones del modelo. El modelo (1.9), (1.10), (1.11)
y (1.12) se construyó bajo la hipótesis de que el potencial de membrana sea
uniforme a lo largo del axón, o se considere únicamente un pedazo pequeño
de membrana. Para estudiar situaciones en las que el potencial de membrana
varı́a temporal y espacialmente a lo largo del axón, este modelo no se puede
aplicar. Usando resultados de la teorı́a del cable se puede construir un sistema de ecuaciones diferenciales parciales no lineales que tomen en cuenta la
dependencia espacial y generalicen el sistema (1.9), (1.10), (1.11), (1.12).
En este caso el sistema de ecuaciones que se obtiene es
1 ∂ 2V
∂V
=
C
+ ḡN a m3 h(V − VN a )2πr
R ∂x2
∂t
+{ḡk n4 (V − Vk ) + ḡe (V − Ve )}2πr ,
(1.13)
∂m
= αm (1 − m) − βm m ,
∂t
(1.14)
∂h
= αh (1 − h) − βh h ,
∂t
(1.15)
∂n
= αn (1 − n) − βn n ,
(1.16)
∂t
donde r es el radio del axón, R es la resistencia por unidad de longitud en el
interior del axón y C la capacitancia de la membrana por unidad de longitud.
Para un resumen de los experimentos realizados por Hodking y Huxley y
de los resultados obtenidos consultar [21]. También se puede consultar [40].
El modelo es importante ya que describe la evolución temporal del potencial de membrana en función de las corrientes iónicas subyacentes. Sin
embargo, no toma en cuenta los eventos moleculares y quı́micos subyacentes
al proceso, y por lo tanto no es posible predecir cómo influyen en el proceso
los cambios en el nivel quı́mico y molecular. Además no es un modelo universal, es decir, se aplica únicamente al potencial de membrana de una neurona
del calamar gigante.
21
El modelo de H-H incorpora únicamente tres corrientes iónicas, sodio,
potasio y una corriente de filtración, y debido a esto, el modelo consiste de un
sistema de 4 ecuaciones diferenciales ordinarias. Existen otros modelos donde
se incorporan más corrientes iónicas, por ejemplo, Bucholtz, [16] elabora un
modelo de la L. P. neurona en el ganglio estomatogástrico del cangrejo. Este
modelo consiste de un sistema de 14 ecuaciones diferenciales ordinarias.
Otros neurofisiólogos elaboran modelos de neuronas que siguen esencialmente los principios del modelo de H-H, los cuales incorporan más corrientes
iónicas o toman en cuenta las variaciones de las concentraciones de ciertos
iones en el interior de la membrana. En [21] se presentan otros modelos de
neuronas, también se pueden ver otros modelos en [99], [100].
El estudio matemático de estos modelos es bastante complejo debido al
número de ecuaciones diferenciales que aparecen. Por lo general son “intratables” analı́ticamente y a veces sólo se pueden encontrar soluciones por
integración númerica. Se han realizado algunos estudios cualitativos de algunos de estos modelos [58], [111].
Una técnica que se usa frecuentemente es reducir, mediante algún mecanismo, el número de dimensiones del modelo manteniendo las principales
caracterı́sticas del mismo. Por ejemplo, Fitzhugh [28], [29], propone un modelo bidimensional para el modelo de H-H Abbot, Kepler [2] proponen un
método de reducción más sistemático que el de Fitzhugh, ya que reducen el
modelo de H-H a dos dimensiones, y Golomb, Guckenheimer, Gueron [39],
aplican este método para reducir el modelo de Buchholtz a 7 dimensiones.
Doya, Selverston [26] elaboran un método de reducción de variables usando
redes neuronales artificiales.
1.7 REDES NEURONALES
El estudio de la estructura y funcionamiento de una neurona es un paso
importante en el estudio del cerebro, pero el cerebro está compuesto de
aproximadamente 1011 neuronas densamente interconectadas, organizadas en
grupos de diversos tamaños. Las neuronas de un mismo grupo pueden interactuar entre sı́ y también pueden interactuar los grupos entre sı́ y frecuentemente se dan estos tipos de interacción de manera simultánea. La estructura
22
de estos grupos y sus mecanismos de interacción no son muy conocidos. Los
estudios experimentales de estos grupos de neuronas son difı́ciles de realizar.
Los modelos de redes neuronales (un conjunto de unidades elementales
interconectadas) son una herramienta importante para resolver muchos problemas en el área de la ingenierı́a de la computación, matemáticas, neurofisiologı́a, etc. En cada una de estas áreas se interesan en diferentes aspectos de
las redes neuronales, dependiendo de los problemas que se quieran resolver.
A los neurofisiólogos les interesa estudiar redes neuronales con la finalidad
de estudiar el cerebro humano. Ellos desean construir redes neuronales “realistas” que imiten ciertas capacidades del cerebro como las sensoriales, las de
aprendizaje y las relacionadas con la memoria. También desean construir redes que describan los mecanismos subyacentes a los patrones rı́tmicos que se
detectan en el cerebro. Para los ingenieros, computólogos y matemáticos, interesa construir redes neuronales con fines computacionales. A estas últimas
se les conoce como redes neuronales artificiales y frecuentemente tienen poca
relación con lo que ocurre en el cerebro.
Las redes neuronales se construyen siguiendo ciertas propiedades que tiene
el cerebro. Están constituidas por un conjunto de n unidades elementales
llamadas neuronas artificiales, unidades de procesamiento o nodos que están
conectados entre sı́.
Las propiedades de las redes dependen del número de neuronas de que
consta la red y de la forma en que se interconectan.
Cada nodo se modela considerando algunas caracterı́sticas de las neuronas reales. Como mencionamos anteriormente una neurona real recibe la
información por las dendritas, la procesa en el cuerpo celular y bajo ciertas
condiciones emite potenciales de acción.
Si la red tiene n nodos, a cada nodo i se le asignan tres cantidades:
la entrada neta Si al nodo i, Ui = Ui (t), el estado interno del nodo en el
instante t (potencial de membrana de una neurona real), la señal de salida
yi del nodo i (número de potenciales de acción que dispara la neurona por
unidad de tiempo) que es función del estado interno, es decir
yi = fi (Ui ) ,
23
donde fi transforma el estado interno Ui en una señal de salida yi .
La entrada Si al nodo i, usualmente se considera de la forma
Si =
n
X
Wij yi ,
j=1
donde Wij es el tercer tipo de magnitud asociada a un nodo y se llama peso
sináptico. Esta cantidad nos indica el grado de conexión entre el nodo i y el
nodo j. El peso Wij puede ser constante o dependiente del tiempo y puede ser
positivo o negativo según corresponda a una conexión excitatoria o inhibitoria
respectivamente. Se considera que el conocimiento se encuentra representado
en los pesos sinápticos y que los procesos de aprendizaje implican cambios
en los pesos sinápticos.
Los modelos de redes neuronales que se obtienen, dependen de cómo se
relacionan las cantidades asignadas a cada nodo.
El primer modelo que se construyó fue el modelo de Mc Culloch-Pitts [80]
que tiene las siguientes caracterı́sticas.
1. Wij constantes.

n
 1 , si P W y ≥ θ (la neurona “dispara”), θ umbral de la neurona i
i
ij i
i
2. yi = fi (Ui ) =
i=1

0,
en el otro caso .
En este caso se dice que el nodo o neurona artificial es binario, es decir
está activado o está desactivado dependiendo de si la magnitud Si rebasa un
umbral o no.
En otros modelos lo que se construye es una ecuación diferencial en Ui
que describe la evolución temporal de Ui :
dUi
= Gi (Ui , Si ), i = 1, 2, ..., n .
dt
Por ejemplo en el modelo de Hopfield [59], [60]:
τi
dUi
= −αi Ui + (Si + θi ) ,
dt
24
Si =
n
X
Wij gj (Uj ) con g, una función sigmoide .
j=1
En el modelo de Grossberg [43]:

τi
n
X

dUi
= −αUi + (γi − βi Ui )  Wij gi (Uj ) + θj  .
dt
j=1
Cuando se quiere estudiar redes neuronales que aprendan, es decir que
sean capaces de modificar sus pesos sinápticos, se debe agregar un sistema
de ecuaciones diferenciales que describan la evolución temporal de los pesos
sinápticos.
Las redes neuronales se forman con un conjunto de neuronas interconectadas. Dependiendo de cómo se agrupan las neuronas, las redes se clasifican
en:
1. Redes neuronales monocapa (1 capa). Una capa o nivel es un conjunto
de neuronas cuyas entradas provienen de la misma fuente (que puede
ser otra capa de neuronas) y cuyas salidas se dirigen al mismo destino.
2. Redes neuronales multicapa.
Con respecto a las conexiones en la red, si ninguna salida de las neuronas
es entrada de la misma capa o capas precedentes, se dice que la red es de
propagación hacia adelante. Cuando las salidas pueden ser conectadas como
entradas de neuronas de niveles previos o del mismo nivel, incluyéndose ellas
mismas además de conexiones con neuronas de las capas posteriores, se dice
que la red es de propagación hacia atrás.
Con respecto al aprendizaje de las redes neuronales, modificación de los
pesos sinápticos, las redes neuronales se pueden clasificar como: redes neuronales con aprendizaje supervisado y redes neuronales con aprendizaje no
supervisado. En las primeras se precisa de un agente externo (supervisor)
que controle el proceso de aprendizaje de la red, y que determina la respuesta
de la red a partir de una entrada determinada, en las segundas la propia red
controla su propio proceso de aprendizaje.
En ambos tipos de redes se requiere conocer las reglas de aprendizaje
–criterios que se siguen para modificar los pesos sinápticos–. El aprendizaje
25
supervisado se puede llevar a cabo de tres maneras: aprendizaje por corrección de error, por esfuerzo y estocástico. En el aprendizaje no supervisado
se puede llevar a cabo, entre otras maneras, el aprendizaje Hebbiano y el
aprendizaje competitivo y cooperativo. Para mayor información sobre este
tema recurrir a [11], [12], [31] y [55].
El primer modelo de red neuronal artificial, llamada Perceptron, fue desarrollado por Rosenblatt [98]. El Perceptron es una red neuronal con aprendizaje supervisado, con regla de aprendizaje de corrección de error, de propagación hacia adelante que consiste de dos capas, una de entrada - que tiene
n neuronas lineales - y la otra capa, que consta de una sola neurona, es la
capa de salida. El Perceptron es capaz de decidir cuándo una entrada presentada a la red pertenece a una de las dos clases que es capaz de reconocer.
La única neurona de salida del Perceptron realiza la suma ponderada de las
entradas, resta el umbral y pasa el resultado a una función de transferencia
de tipo escalón. La regla de decisión es responder +1 si el patrón presentado
pertenece a la clase A, o −1 si patrón pertenece a la clase B; A y B son las
dos clase que el Perceptron es capaz de reconocer, previo entrenamiento. El
Perceptron funciona cuando el conjunto de n−adas, cuyas componentes son
las entradas a las n−neuronas de la capa de entrada, pueden ser separadas
P
por un hiperplano n−dimensional - el correspondiente a N
i=1 Wi Xi − Θ = 0;
donde Wi es el peso correspondiente a la entrada Xi y Θ el umbral. En general
esto no siempre se puede hacer. Posteriormente se desarrolla el Perceptron
multicapa [11], [31].
El Adaline y el Madaline son redes neuronales desarrolladas por Widrow
y Hoff [118] en 1960.
La arquitectura del Adaline y del Madaline son esencialmente la misma
que la del Perceptron. Lo que cambia es la regla de aprendizaje ya que
Adaline y Madaline usan la llamada regla delta de Widrow - Hoff o regla del
mı́nimo error cuadrado (LMS), basada en la búsqueda de una expresión del
error entre la salida deseada y la salida lineal obtenida antes de aplicarle la
función de activación escalón (a diferencia de la salida binaria utilizada en
el caso del Perceptron). Sus principales aplicaciones están en el campo del
procesamiento de señales.
También existen redes con aprendizaje supervisado del tipo estocástico,
26
la máquina de Boltzman [56], etc., del tipo aprendizaje con refuerzo, Linear
Reward-Penalty o LR−P .
Con respecto a las redes con aprendizaje no supervisado se tienen: la
red de Hopfield, Learning Matrix, etc. que utilizan el aprendizaje Hebbiano
(ajuste de los pesos de las conexiones de acuerdo con la correlación de los
valores de activación de las dos neuronas conectadas).
Para las que utilizan el aprendizaje competitivo - cooperativo (las neuronas compiten o cooperan unas con otras con el fin de llevar a cabo una
tarea) ver [5], [7].
Para mayor información sobre la arquitectura, reglas de aprendizaje, aplicación y limitaciones de las redes que hemos mencionado y a otras más recurrir a: [31], [55].
Cuando se intenta construir modelos de redes neuronales para el estudio del cerebro, los nodos de las redes se pueden interpretar como una sola
neurona o como un grupo de neuronas.
Matemáticamente no existe ninguna distinción entre estos dos casos, pero
si se desea construir redes neuronales que describan ciertos fenómenos que
ocurren en el cerebro y que los resultados se correlacionen con mediciones
experimentales, la distinción sı́ es importante.
Se ha detectado actividad rı́tmica en algunas estructuras del cerebro:
oscilaciones en el tálamo, oscilaciones en el sistema bulbo olfativo, corteza,
etc. Algunos modelos de redes neuronales se elaboran tomando en cuenta
datos anatómicos y fisiológicos con la finalidad de describir estos fenómenos
temporales rı́tmicos. Por ejemplo:
Lopes da Silva [76], [77] elabora un modelo de red neuronal que intenta
explicar los orı́genes del ritmo α en el perro con redes neuronales creadas con
neuronas de la corteza y del tálamo.
Destexhe, Contreras, Sejnowski, y Steriade [24] elaboran un modelo de red
neuronal para neuronas del núcleo reticular talámico y obtienen resultados
teóricos de tipo oscilatorio con frecuencias de 7 a 14 Hz y otros en el rango
de 0.5 - 1 Hz dependiendo de los parámetros que se consideran.
27
Freeman [33], [34], [35], elabora un modelo que describe el modo en que
la información sensorial puede ser procesada en un modo distribuido para
una de las formas más simples de corteza. El modelo se elabora con base en
estudios anatómicos y fisiológicos sobre las interconexiones neuronales entre
el bulbo olfativo y la corteza en gatos y conejos e intenta relacionar los resultados teóricos con datos electroencefalográficos obtenidos en experimentos
con estos animales.
Menon [81], estudia el fenómeno de las oscilaciones temporales en el cerebro, considerando que éstas se deben a la interacción de neuronas inhibitorias
y neuronas excitatorias. Elabora un modelo en el cual considera dos grupos
de neuronas, una que tiene únicamente conexiones inhibitorias y el otro grupo
que tiene únicamente conexiones excitatorias. Estudia la dinámica de estos
dos grupos de neuronas interconectados elaborando un modelo de dos ecuaciones integro-diferenciales acopladas.
Nuñez comenta en [89], [92] que la mayorı́a de los modelos de redes neuronales no surgen para estudiar las bases fisiológicas de los EEG en el cuero
cabelludo sino para estudiar algunas funciones cerebrales, como por ejemplo
la memoria de corto plazo. También menciona que si se quiere “explicar”
por medio de redes neuronales los patrones rı́tmicos temporales que se miden con electrodos en el cuero cabelludo (EEG), no basta que se construyan
redes (aún basados en datos anatómicos y fisiológicos) que produzcan soluciones oscilatorias con la misma frecuencia que se obtienen en los EEG. Por
lo tanto, es necesario considerar la estructura espacial de la red para considerar los potenciales correspondientes y compararlos con los que se obtienen
en los EEG, ya que la magnitud de un potencial depende de la configuración
geométrica de sus fuentes.
Nuñez elabora un modelo de interacción neuronal en la corteza a la cual
considera como un medio continuo, donde la columna de corteza es la unidad
fundamental. El modelo consiste en un par de ecuaciones integrales acopladas
que relacionan la actividad sináptica excitatoria e inhibitoria en una columna
de corteza en un tiempo t con los potenciales de acción que se “disparan”
en las otras columnas en tiempos anteriores, tomando en cuenta las conexiones que existen entre las columnas y las velocidades de propagación de los
potenciales de acción, ver [87], [88], [89] y [92].
28
Nuñez elabora el modelo con la finalidad de estudiar las bases fisiológicas
de los patrones espacio temporales medidos con electrodos ubicados en el
cuero cabelludo. Considera que para el estudio de estos patrones espaciotemporales la columna de corteza es la unidad a considerar, ya que los EEG se
obtienen con electrodos que captan la actividad eléctrica promedio de varios
centı́metros cuadrados de corteza, además de la homogenidad que ofrece la
corteza a nivel de cada capa.
En el siguiente capı́tulo analizaremos con más detalle este modelo.
29
1.8 MODELOS DE MEDIO CONDUCTOR
En el funcionamiento del cerebro, la actividad eléctrica es un elemento
fundamental; la comunicación entre las neuronas y entre estructuras del
cerebro se realiza a través de ella. Otro elemento fundamental en el funcionamiento del cerebro es la actividad neuroquı́mica. Se conoce muy poco
de cómo interactúan ambos elementos en el funcionamiento del cerebro.
La actividad eléctrica del cerebro puede ser medida en muchas escalas
espacio-temporales y con diversos grados de resolución se puede medir la
actividad eléctrica de una sola neurona, de grupos de neuronas, y usando
electrodos colocados en el cuero cabelludo, la del cerebro en su totalidad.
Se conocen diversos patrones espacio-temporales de la actividad eléctrica
en distintos niveles, y se desea estudiar las bases fisiológicas subyacentes a
ellos.
Se han usado diversos modelos: redes neuronales y sistemas cooperativoscompetitivos de interacción neuronal, para describir los posibles mecanismos
que originan estos patrones espacio-temporales. Pero estos modelos, deben
ir acompañados de otros modelos que toman en cuenta otras caracterı́sticas
del cerebro como son las propiedades electromagnéticas derivadas de la activación de fuentes de corrientes impresas en el cerebro y por ello es necesario
estudiar cómo se propagan las corrientes eléctricas a través de los distintos
tejidos del cerebro.
No basta con determinar que ciertas redes neuronales correspondientes a
ciertas estructuras del cerebro producen oscilaciones teóricas que coinciden
con las oscilaciones detectadas con los EEG medidos en diversas partes del
cerebro. Los potenciales dependen de la distribución de las fuentes y de las
propiedades eléctricas del medio en que ellos se propagan.
En el estudio del cerebro, como medio conductor, una herramienta importante es la teorı́a del campo electromagnético que en el nivel microscópico
nos describe la relación entre cargas, corrientes eléctricas, campos eléctricos
y magnéticos y en el nivel macroscópico debe incluir ciertas relaciones experimentales para la permitividad (para aislantes), la conductividad (para
conductores) y la permeabilidad (para materiales magnéticos). Pero nuevamente, si se quieren elaborar modelos cuyos resultados se puedan relacionar
30
con datos experimentales, es necesario determinar en qué escala espacio temporal se está trabajando.
En cualquier escala espacio temporal que se use, las mediciones electroencefalográficas se obtienen como promedio de la actividad eléctrica de
un cierto volumen de tejido y en cierto perı́odo de tiempo y por tanto, las
conductividades, las permitividades y las permeabilidades magnéticas de los
tejidos deben considerarse en la misma escala. Además las variables de nuestro modelo deben estar en la misma escala espacio-temporal que la de las
mediciones.
El problema que se requiere estudiar es determinar la ubicación y
propiedades de las fuentes de corriente generadas por la actividad del cerebro a partir de las mediciones eléctricas y magnéticas efectuadas en el cuero
cabelludo o en otras estructuras del cerebro. Este problema no tiene solución
única; en el sentido de que la misma medición de potencial o campo magnético
puede ser generada por diferentes distribuciones de fuentes de corriente concentradas en la masa cerebral.
El problema inverso no tiene solución única, aunque se conozca la distribución de potencial (o campo magnético) en todos los puntos del cuero
cabelludo. Pero la situación real es todavı́a peor ya que las mediciones del
potencial se obtienen únicamente en un conjunto finito de puntos del cuero
cabelludo donde están ubicados los electrodos.
El problema es sumamente complejo ya que es necesario considerar la
geometrı́a y propiedades electromagnéticas de las distintas regiones que intervienen –corteza, cráneo, cuero cabelludo, etc.– con el objetivo de relacionar
los datos obtenidos en el cuero cabelludo con sus posibles fuentes.
Aún el problema directo, es decir, dada la distribución de fuentes de corriente, la conductividad de los diversos tejidos y la permeabilidad magnética,
determinar el campo electromagnético en cualquier punto, ya es bastante
complejo.
Usualmente lo que se hace, en ambos problemas, es introducir simplificaciones con base en información conocida sobre el cerebro y resolver el
problema con esas simplificaciones. Estas simplificaciónes dependen mucho
de la escala en que se está trabajando.
31
La primera simplificación importante que se introduce es el considerar
que se puede aplicar al cerebro la teorı́a electromagnética cuasi-estática [92].
Es conocido que en un medio conductor con conductividad σ las cargas se
distribuyen muy rápidamente, es decir, la densidad ρ de carga decrece como
σ
ρ ∼ e− εo t ,
donde εo es la constante dieléctrica. Para la conductividad del hueso
1
del orden de 200
y para la piel y la masa cerebral es mucho menor.
Por este motivo se tiene que
σ
εo
es
¯
¯
¯ ∂ρ ¯
σ − εσ t
¯
¯
¯
¯∼
e o
¯ ∂t ¯
εo
y como se conoce experimentalmente que el intervalo de tiempo en que la
densidad de corriente J~ en el cerebro varı́a apreciablemente, es mucho mayor
que εσo , entonces en la ecuación de continuidad
∂ρ
∇ · J~ +
=0
∂t
se puede despreciar
∂ρ
∂t
y queda
∇ · J~ = 0 .
Además se considera que las corrientes que circulan en el cerebro son de
dos tipos: Ohmicas, es decir, generadas por cargas en movimiento a través
del fluido extracelular y corrientes primarias o impresas, J~p , que se mueven a
través de las membranas celulares y aparecen en las zonas activas. Entonces
se tiene
~,
J~ = J~p − σ E
~ es el campo eléctrico.
donde E
Del desacoplamiento entre el campo eléctrico y magnético se tiene que
~ = 0 y, por lo tanto, E
~ = ∇u, donde u es el potencial del campo
rot E
electrostático.
32
De la ecuación de continuidad simplificada se obtiene que en la masa del
cerebro se cumple la ecuación
∇u =
1
∇ · J~p .
σ
Si además se considera un sistema de varias capas donde la región Ω1
ocupada por el cerebro está cubierta por otras regiones disjuntas Ωi (i ≥ 2)
que representan el fuido intracraneal, el cráneo, la piel, etc., entonces en
cada frontera de separación entre dos regiones Ωi y Ωi+1 , se deben satisfacer
condiciones del tipo
Ã
σi
∂u
∂ni
!−
Ã
= σi+1
∂u
∂ni
!+
+ f~i ,
u− = u+ ,
donde σi es la conductividad en cada Ωi , ni es el vector unitario normal
exterior a Ωi en la parte común de las fronteras de Ωi y Ωi+1 . El supraı́ndice
− o + indica que el valor lı́mite a la frontera de la cantidad indicada se toma
desde Ωi o desde Ωi+1 respectivamente.
Mediante f~i se denota la componente normal de la densidad de corriente
superficial, la cual en este modelo puede estar presente únicamente en la
frontera del cerebro, es decir, en la corteza cerebral.
Si se denota por Ωe la región cuya frontera exterior Se coincide con el
cuero cabelludo, entonces el problema único consiste en determinar J~p y f~1
a partir de las mediciones de u sobre Se .
Este problema está mal planteado y por ello se necesita introducir algoritmos de regularización en su solución.
Una de las simplificaciones más importantes en el estudio del problema
inverso consiste en representar f~1 y J~p en forma de un número finito de
dipolos puntuales y, en este caso, el problema se reduce a la localización de
sus posiciones y momentos, es decir, a encontrar 6N parámetros donde N es
el número de dipolos propuestos para la fuente activa [103].
33
1.9 CONCLUSIONES
Como hemos podido ver, existe un gran número de técnicas y métodos
para determinar anomalı́as de carácter anatómico o funcional en el cerebro,
ası́ como para el estudio de los patrones temporales ritmicos de la actividad
eléctrica del mismo. Sin embargo, ninguno de estos métodos puede describir
la interacción neuronal entre diferentes regiones del cerebro ni puede explicar el origen fisiológico de los patrones temporales rı́tmicos de la actividad
eléctrica.
El objetivo de este trabajo es estudiar modelos que son una generalización
de los introducidos por Paul Nuñez en [89], [92] que, en conjunto con los
modelos de medio conductor, permitan describir la evolución temporal de
la interacción neuronal masiva en la corteza cerebral y también dar alguna
explicación matemáticamente fundamentada de los orı́genes fisiológicos de
los patrones temporales rı́tmicos de la actividad eléctrica.
Esperamos además que este tipo de trabajo sirva como fundamento para
el desarrollo de software más apropiado para el análisis clı́nico de los EEG.
34
Capı́tulo 2
MODELACIÓN GLOBAL DE LA ACTIVIDAD SINÁPTICA DE LA CORTEZA CEREBRAL
2.1 INTRODUCCIÓN
Para el estudio de las correlaciones existentes entre ciertos estados fisiológicos y las mediciones electroencefalográficas –a nivel de cuero cabelludo
o de corteza– se han elaborado diversos modelos de interacción neuronal a
nivel de corteza-tálamo. Estos modelos se pueden dividir, grosso modo, en
modelos locales, modelos globales y modelos mixtos.
En los modelos locales [33], [34], [75], [76], [77], [95], se tratan de explicar
las propiedades espacio-temporales de los EEG con base en propiedades a
nivel de neurona, redes neuronales discretas o redes neuronales distribuidas
en el rango de milı́metros; las propiedades temporales de los EEG son determinadas por los retrasos sinápticos como son los tiempos de respuesta y de
decrecimiento de los potenciales postsinápticos.
En los modelos globales [87], [88], [89], [90], [92], se toman como elemento importante los sistemas de fibras corticorticales (cuyas longitudes se
encuentran en el rango de centı́metros) y el hecho de que la corteza es una
estructura cerrada. En estos modelos se considera que el EEG medido sobre el cuero cabelludo es fundamentalmente el resultado de la interacción de
neuronas corticales a través de los potenciales de acción. Se considera que
35
ciertos patrones espacio-temporales de los EEG se pueden explicar tomando
en cuenta los retrasos debido a la velocidad finita de propagación de los potenciales de acción que viajan por las fibras corticorticales y por el hecho de
que la corteza es una estructura cerrada.
Nuñez considera [92] que los modelos locales pueden ser importantes en
ciertos estados fisiológicos del cerebro y que en otros casos, como en aquellos
que interviene un mı́nimo de actividad cognitiva (sueño, coma, anestesia) se
pueden explicar con los modelos globales; pero que dada la extrema complejidad de los interacciones neuronales se requieren modelos que toman en
cuenta la interacción neuronal en las diversas escalas espacio-temporales [64],
[65], [66], [67], [68].
En este capı́tulo analizaremos el modelo global elaborado por Nuñez. Este
modelo será el punto de partida para la elaboración de nuestro modelo. El
modelo describe la evolución de la actividad sináptica en la corteza y consiste
en un par de ecuaciones integrales que relacionan la actividad sináptica en
cada columna de corteza en un instante de tiempo t con los potenciales de
acción que se generan en tiempos anteriores en las demás columnas de corteza,
tomando en cuenta las conexiones que existen entre las columnas de corteza
y las velocidades de propagación de los potenciales de acción. También se
pueden tomar en cuenta los retrasos sinápticos o incluir circuitos locales,
convirtiéndolo en un modelo local-global más completo.
Con el modelo global, sin tomar en cuenta los retrasos locales, se estudia
bajo qué condiciones se tienen soluciones de tipo ondulatorio y se comparan
estas soluciones con los datos obtenidos a través del EEG. Sus resultados
teóricos “coinciden” con los obtenidos experimentalmente para el ritmo α
pero no incluye otros estados fisiológicos.
2.2 BASES FISIOLÓGICAS PARA LA MODELACIÓN DE
LA ACTIVIDAD SINÁPTICA EN LA CORTEZA
CEREBRAL
La corteza cerebral es un sistema muy complejo que presenta interacciones
en distintos niveles: neurona, grupos de neuronas y entre grupos de neuronas.
Aun ası́ la corteza presenta ciertas caracterı́sticas anatómicas y fisiológicas
36
que permiten elaborar modelos de interacción neuronal cuyas variables y
parámetros se pueden correlacionar con las medidas electroencefalográficas
obtenidas en el cuero cabelludo.
El modelo que presentaremos toma como unidad fundamental la columna
de corteza y pretende describir las interacciones entre las columnas de corteza,
a través de sus interconexiones y la interacción de la corteza con otras estructuras del cerebro (principalmente el tálamo); más precisamente se describirá
la evolución espacio-temporal de la actividad sináptica de las columnas de
corteza.
Los elementos del modelo son:
a) h± (x, t) − La densidad de sinapsis activas o de actividad sináptica (+
para las sinapsis de tipo excitatorio, − para las de tipo inhibitorio), en
la columna x en el tiempo t. h± (x, t) representa el número de sinapsis
activas por cm3 en la columna x en el tiempo t (1/cm3 )
b) g(x0 , t) − La fracción del total de neuronas de la columna x0 que generan potenciales de acción en el tiempo t (adimensional)
c) R± (x0 , x, v) − El número de conexiones sinápticas por cm3 en la columna
x por unidad de área en x0 y por unidad de velocidad, es decir,
R± (x0 , x, v)dx0 dv representa el número de conexiones sinápticas por cm3
en la columna x producidas por las neuronas que se encuentran en las
columnas de x0 a x0 + dx0 y por cuyos axones los potenciales de acción
“viajan” con velocidad entre v y v + dv (+ para las conexiones de tipo
excitatorio, − para las de tipo inhibitorio).
d) S± (x) − La densidad de sinapsis en la columna x (+, para las sinapsis
excitatorias, − para las inhibitorias) (1/cm3 ). Se tiene que:
Z ∞Z
0
S
R± (x0 , x, v)dx0 dv = S± (x) ,
donde S es la superficie de la corteza.
37
e) F± (v) − La fracción del total de neuronas de la corteza por unidad de
velocidad (seg/cm), es decir, F± (v)dv representa la fracción del total de
neuronas de la corteza que generan su potencial de acción con velocidad
entre v y v + dv, (+ para las conexiones sinápticas excitatorias, − para
las inhibitorias).
f ) h0± (x, t) − La densidad de sinapsis activas o de actividad sináptica ( +
para las sinapsis de tipo excitatorio, − para las de tipo inhibitorio),
en la columna x en el tiempo t. Actividad sináptica producida por los
potenciales de acción que se generan en otras estructuras del cerebro,
fundamentalmente el tálamo y que “llegan” a la columna x en el tiempo
t. h0± (x, t) representa el efecto producido en la columna x en el tiempo
t por las estructuras internas del cerebro que tienen conexiones en la
columna x.
Con estas definiciones es fácil ver que si S es la superficie de la corteza
cerebral, entonces se tienen las ecuaciones:
h+ (x, t) =
ho+ (x, t)
h− (x, t) =
ho− (x, t)
+
+
Z ∞Z
0
Z ∞Z
0
d(x, x0 )
R+ (x , x, v)g(x , t −
) dx0 dv ,
v
S
(2.1)
d(x, x0 )
) dx0 dv ,
R− (x , x, v)g(x , t −
v
S
(2.2)
0
0
0
0
donde d(x, x0 ) es la distancia geodésica entre dos puntos de la corteza S.
Este sistema de ecuaciones integrales nos dice que la densidad de actividad
sináptica tiene dos componentes. La primera de ellas es debida al “input”
que reciben la corteza y las fibras corticorticales proveniente de estructuras
interiores durante un estado fisiológico fijo. La componente integral expresa
la idea de que los potenciales de acción que se disparan en la columna x0 se
trasladan a través de las fibras corticorticales causando actividad sináptica
en la columna x en un tiempo posterior que depende de la distancia que
separa x de x0 y de la velocidad de propagación de los potenciales de acción.
38
En general esta distancia no tiene que coincidir con la euclideana aunque
puede ser considerada como tal para el análisis cualitativo.
A continuación enumeraremos las hipótesis de trabajo y las principales
consideraciones anatómicas para el trabajo con este modelo.
1. No están considerados los retrasos debidos a la transmisión sináptica
–retrasos locales– sino, únicamente el debido a la velocidad finita de
propagación de los potenciales de acción a lo largo de los axones (fibras
corticorticales).
2. En el modelo (2.1), (2.2) se desprecia la existencia de circuitos de
retroalimentación entre la corteza y el tálamo ya que el número de fibras corticorticales es mucho mayor que las fibras talamocorticales. En
efecto las primeras representan el 99% del total [72]. No resulta dı́ficil
incluir en este modelo global la presencia de estos circuitos de retroalimentación. Un ejemplo importante de estos circuitos es el circuito
córtico-estriato-nigro-tálamocortical asociado a la actividad motora.
3. Las neuronas corticales están colocadas en columnas que presentan una
gran homogeneidad en su distribución espacial con respecto a la profundidad. Esto quiere decir que es imposible distinguir entre dos cortes
perpendiculares realizados en dos secciones diferentes de la corteza cerebral. Por este motivo podemos considerar que la corteza tiene una estructura estratificada uniforme compuesta por seis capas debido a las
diferencias en las densidades relativas de diferentes tipos de neuronas
localizadas a distintas profundidades corticales y que en una primera
aproximación las conexiones en el plano de la corteza son del mismo
tipo en todas las direcciones y por ello sólo dependen de la distancia
entre las columnas consideradas es decir:
R± (x0 , x, v) = R± (d(x, x0 ), v) .
(2.3)
Las neuronas corticales se pueden clasificar de acuerdo con la forma
del cuerpo celular, por la longitud de su axón y por la extensión y distribución espacial del árbol dendrı́tico. Además del tipo de transmisor
39
que libera, el blanco axonal, el tipo de sinapsis, etc. Para un estudio
más detallado sobre esta tema recurrir a [82].
En primera aproximación, las neuronas corticales las dividiremos en
dos grandes tipos: excitatorias (neuronas piramidales y estrelladas espinosas) e inhibitorias (todas las interneuronas). Las interneuronas
inhibitorias se dividen en varios tipos y se consideran que pueden tener
varias escalas espacio–temporales. Pero en nuestro caso tomaremos que
la inhibición ocurre rápidamente.
Las piramidales tienden a ocupar volúmenes cilı́ndricos con axones y
dendritas alineadas perpendicularmente a la superficie cortical y aproximadamente de 2/3 a 3/4 de las neuronas corticales son de este tipo.
Las neuronas estrelladas ocupan volumenes más esféricos.
El cuerpo celular de las neuronas ocupa alrededor de un 10% de la
superficie neuronal. La mayor parte del área superficial de la neurona
está compuesta por sus ramificaciones dendrı́ticas. El número de ramas
dendrı́ticas que parten de una neurona y que intersectan a una superficie esférica concéntrica, centrada en el cuerpo celular de la neurona,
decrece exponencialmente según una ley del tipo:
Número de ramas ≈ exp(−λ· distancia radial) donde λ−1 (cm) es la
escala caracterı́stica de longitud que determina el decrecimiento en e−1
veces del número total de conexiones corticales a partir de una neurona.
Por ejemplo, para las células estrelladas y las dendritas basales de las
células piramidales en las cortezas visual y motora del gato y el ratón
se ha comprobado experimentalmente que λ = 350 cm−1 .
La relación para el decrecimiento de las ramas generadas a partir de
una dendrita es de orden polinomial en lugar de exponencial.
4. Las interconexiones entre neuronas corticales pueden ser divididas en
dos grandes grupos: las conexiones intracorticales que son de rango
corto y que tienen una longitud promedio de menos de 1 mm en humanos y las conexiones corticorticales que tienen una longitud promedio
de hasta varios centı́metros en humanos.
La materia blanca subcortical está compuesta por las fibras de asociación corticorticales que están ubicadas paralelamente a la superficie
40
cortical (Fig. #7).
Mientras que las conexiones intracorticales pueden ser tanto excitatorias como inhibitorias, las fibras de asociación (axones de las células
piramidales) que constituyen la mayor parte de la materia blanca son
solamente conexiones de tipo excitatorio. Las neuronas piramidales que
envı́an axones a otras regiones interiores del cerebro como el tálamo,
terminan únicamente en sinapsis excitatorias.
En las neuronas piramidales existen otro tipo de neuronas intracorticales que también terminan en sinapsis excitatorias. Sin embargo la
cantidad de sinapsis excitatorias en humanos realizadas por las fibras
de asociación (corticorticales) es mucho mayor que las que son producidas por axones intracorticales, aproximadamente el 80% del total
de las sinapsis excitatorias son por las fibras de asociación. Aunque en
ratones es del orden del 50%.
Para mayor información sobre las fibras de asociación se puede consultar [72].
5. La mayorı́a de las fibras de asociación mielinadas tienen diámetros entre
1 y 2 µ y en ellas las velocidades de los potenciales de acción son
proporcionales al diámetro y fluctúan entre 6 y 12 m/seg.
6. El área de la superficie total de ambas mitades del cerebro, las cuales
están separadas por el cuerpo calloso es del orden de 1600 cm2 . Las 2/3
partes de esta superficie están en las fisuras incluyendo la fisura sagital
que separa las dos mitades del cerebro. De estas dos terceras partes
la mitad del área está contenida en las pequeñas fisuras incluyendo la
fisura sagital. Esto nos dice que el área efectiva del medio donde se
propagan las ondas es diferente de la superficie cortical.
Un problema importante en los modelos globales del EEG consiste en
saber cuándo ambas mitades del cerebro deben ser consideradas por
separado. La materia blanca está formada por fibras de aferentes y
eferentes que van y vienen del interior del cerebro y la médula espinal,
junto con fibras de asociación que conectan diferentes áreas corticales
41
y fibras comisurales que conectan las dos mitades del cerebro fundamentalmente a través del cuerpo calloso.
En humanos el número total de fibras comisurales que van de un hemisferio a otro del cerebro a través del cuerpo calloso es del orden de
2 × 108 , la cual es significativamente menor que el número de fibras de
asociación que entran a la corteza, que es del orden de 1010 . Esto hace
pensar que las fibras comisurales permiten una conexión muy débil entre los dos hemisferios cerebrales, por lo cual parece ser que el análisis
de los fenómenos ondulatorios debidos a la actividad cerebral debe restringirse a un sólo hemisferio que puede identificarse con la cuarta parte
de un elipsoide de revolución.
El cerebro humano tiene 45cm de circunferencia (sin contar las fisuras)
de la superficie: del polo occipital al polo frontal, a lo largo de la lı́nea
media. Esta idea resulta consistente con los resultados experimentales
obtenidos en pacientes a quienes se les ha medido la frecuencia del ritmo
α con los dos hemisferios acoplados y con cuerpo calloso cortado. Si
los dos hemisferios acoplados se comportaran como un sistema único
se esperarı́a que al efectuar el corte en el cuerpo calloso hubiera un
significativo incremento de la frecuencia del ritmo α, lo cual no ocurre.
7. Los potenciales del EEG medidos en el cuero cabelludo se deben a una
actividad sincrónica de neuronas corticales que cubren un área de al
menos varios centı́metros cuadrados de corteza. Por este motivo la
medición del EEG corresponde a un nivel jerárquico macroscópico de
modelación donde la homogeneidad e isotropı́a son dos consideraciones
importantes. En cambio, del hecho que las conexiones corticorticales
son inhomógeneas a niveles jerárquicos pequeños se concluye que los
modelos de actividad medida con microelectrodos requiere la inclusión
de parámetros de conexión preferencial.
Como las neuronas piramidales generan grandes momentos dipolares y
están orientadas perpendicularmente a la superficie de la corteza, ellas
son capaces de producir corrientes sincronizadas de suficiente magnitud
para ser medidas sobre el cuero cabelludo. Debido a esto y a la superioridad numérica de las neuronas piramidales junto con el hecho de que
42
son ellas las que producen la mayor parte de sinapsis excitatorias en la
corteza a través de sus axones (fibras corticorticales), es que son estas
neuronas las principales responsables de la medición del EEG sobre el
cuero cabelludo.
De aquı́ concluimos que para una primera modelación de la contribución
global de la corteza al EEG se debe considerar la actividad sináptica excitatoria producida por las neuronas piramidales ası́ como la actividad
sináptica inhibitoria producida sobre ellas por las neuronas estrelladas.
Como mencionamos al inicio de este capı́tulo, en esta tesis nos dedicaremos al estudio de modelos teóricos de la actividad global de la
corteza cerebral que dan lugar al EEG. Estos modelos tienen en cuenta
las consideraciones anatómicas que hemos descrito en los puntos 1-7
anteriores.
A diferencia de ellos, los modelos de circuitos locales deben incluir
parámetros de conexión preferenciales, en lugar de la homogeneidad
que se considera en los modelos de actividad global, lo cual se pone
de manifiesto en la selección de las funciones R± (x0 , x, v) para cada
vı́a del circuito. Si los circuitos son pequeños puede considerarse que
la velocidad de propagación v es infinita no siendo ası́ en los circuitos
en que aparecen interacciones de largo alcance, en los cuales el tiempo
que se demora al viajar un potencial de acción de un lugar a otro del
circuito puede ser del orden de las constantes de tiempo de las neuronas
piramidales.
Por otra parte, como ya mencionamos en el punto 1, los modelos locales deben considerar los retrasos temporales debido a la transmisión
sináptica, que son más pequeños que los retrasos debidos a la velocidad finita de propagación de los potenciales de acción a través de los
axones, y por ello no se consideran en los modelos globales.
43
2.3 DESCRIPCIÓN Y SIMPLIFICACIONES EN EL
MODELO GLOBAL DE ACTIVIDAD SINÁPTICA EN
UNA BANDA DE CORTEZA.
Para evitar repeticiones innecesarias, en este capı́tulo el subı́ndice + siempre se referirá a las sinapsis de tipo excitatorio y el subı́ndice − a las de tipo
inhibitorio
Aplicaremos el modelo general (2.1), (2.2) al estudio de la actividad del
EEG producida por una banda infinitamente delgada de la corteza que por
conveniencia, supondremos infinitamente larga.
Para el análisis cualitativo obviaremos la geometrı́a de la corteza y supondremos que la banda coincide con una recta.
En este caso, que llamaremos el caso unidimensional, los elementos del
modelo son los mismos, sólo cambian las unidades con que se trabajará. Ası́,
h± (x, t) representará el número de sinapsis activas por cm2 en la columna
x en el tiempo t (1/cm2 ). R± (x0 , x, v) representará el número de sinapsis,
por cm2 en la columna x en el tiempo t por unidad de longitud en x0 y por
unidad de velocidad (seg/cm4 ). Es decir, R± (x0 , x, v)dx0 dv representará el
número de sinapsis, por cm2 en la columna x en el tiempo t producidas por
las neuronas que están ubicadas desde la columna x0 hasta la columna x0 +dx0
y por cuyos axones, los potenciales de acción “viajan” con velocidad entre v
y v + dv, h0± (x, t) representará el número de sinapsis activas en la columna
x por cm2 , esta actividad sináptica es la producida por los potenciales de
acción que se generan en estructuras internas del cerebro, principalmente el
tálamo y que llegan en el tiempo t a la columna x, g(x0 , t) seguirá siendo una
cantidad sin dimensiones.
Todas las complicaciones del sistema de conexiones intracorticales e intercorticales están “escondidas” en la representación de las funciones R+ y
R− .
Introduciremos algunas hipótesis adicionales que nos permitirán obtener
una expresión para estas funciones que al menos cualitativamente reune las
caracterı́sticas principales de este sistema de conexiones.
44
Hipótesis 1
Comenzaremos suponiendo que cada columna está constituida por una
cantidad finita N (que no depende de la columna) de sistemas de fibras que
se superponen.
En primera aproximación se podrı́a considerar que N = 6 corresponde a
la distribución de las neuronas en las 6 capas en que usualmente se divide la
corteza.
Hipótesis 2
Como las sinapsis inhibitorias se producen sólo a nivel local, consideraremos que la densidad de actividad sináptica inhibitoria en una columna de
corteza x en un tiempo t se produce por los potenciales de acción que se generan en las columnas “cercanas” a la columna x en tiempos “anteriores”, es
decir, las columnas ubicadas fuera de una vecindad de la columna de corteza
x no contribuyen a la actividad sináptica inhibitoria en la columna x.
Es pertinente mencionar que Nuñez enfoca esta cuestión de manera diferente: considera que la velocidad de propagación de los potenciales de acción
que “viajan” por las fibras que tienen conexión sináptica de tipo inhibitorio
es “infinita”. Esta hipótesis la justifica por el hecho de que la distancia que
recorren los potenciales de acción por estas fibras es muy pequeña. Nosotros
consideramos que esta hipótesis no refleja que las sinapsis inhibitorias se producen sólo a nivel local, lo que refleja es que el “efecto” de los potenciales
de acción que “viajan” por las fibras que tienen conexión sináptica de tipo
inhibitorio es “sentida” instantáneamente en las columnas donde se tienen
estas sinapsis inhibitorias. En realidad, para tener en cuenta el hecho de
que las sinapsis inhibitorias se producen sólo a nivel local, debemos considerar que los potenciales de acción que “viajan” por las fibras que tienen
conexión sináptica de tipo inhibitorio, recorren esas distancias en tiempos
más pequeños que los usados por los potenciales de acción que “viajan” por
las fibras de largo alcance. Más adelante estudiaremos un modelo que tome
en cuenta esta hipótesis.
Desde el punto de vista matemático, nuestra hipótesis (el suponer que
la actividad sináptica inhibitoria en una columna de corteza es afectada
45
únicamente por los potenciales de acción que se generan en columnas cercanas) nos conducirá a estudiar, en el caso lı́mite, un proceso global de
propagación de ondas de actividad sináptica inhibitoria modulada por un
proceso local de propagación de ondas de actividad sináptica inhibitoria;
mientras que en el caso de Nuñez, lo que se estudia es un proceso global
de propagación de ondas de actividad sináptica excitatoria acoplada con un
proceso estacionario de actividad sináptica inhibitoria.
De esta forma definiremos, para cada n = 1, ..., N :
F±n (v) - La fracción del total de neuronas de la corteza por unidad de velocidad (seg/cm) cuyos potenciales “viajan” por las fibras del sistema
n-ésimo. Es decir, F±n (v)dv representará la fracción del total de neuronas de la corteza cuyos potenciales “viajan” por las fibras del n-ésimo
sistema con velocidad entre v y v + dv. Supondremos que esta cantidad
no depende de en cuál columna se ubiquen las neuronas.
Ahora la función R± (x0 , x, v) se puede representar en la forma
R± (x0 , x, v) =
N
X
n
R±
(x0 , x, v) .
(2.4)
n=1
n
R±
(x0 , x, v) (seg/cm2 ) representa la densidad de conexiones sinápticas en
la columna x por unidad de longitud en x0 y por unidad de velocidad ason
ciadas con las fibras del n-ésimo sistema. Ası́ R±
(x0 , x, v)dx0 dv representará
el número de sinapsis por cm2 en la columna x producidas por las neuronas
que están ubicadas desde la columna x0 hasta la columna x0 + dx0 y cuyos
potenciales de acción “viajan” por el n-ésimo sistema de fibras con velocidad
entre v y v + dv.
n
Supondremos que R±
se representa en la forma
n
n
R±
(x0 , x, v) = r±
(x0 , x)F±n (v) ,
(2.5)
n
(x0 , x) (1/cm3 ) será la densidad de conexiones sinápticas en la
donde r±
columna x, asociadas al n-ésimo sistema de fibras, por unidad de longitud en
x0 , ası́ rn (x0 , x)dx0 representará el número de sinapsis por cm2 en la columna
x producidas por neuronas que están ubicadas desde la columna x0 hasta la
columna x0 + dx0 y cuyos axones están ubicados en el n-ésimo sistema de
fibras.
46
n
Ahora r±
(x0 , x) se puede escribir como el producto de dos factores
n
r±
(x0 , x) = tn± (x0 , x)S±n (x) ,
(2.6)
donde tn (x0 , x) (1/cm) representa el número de neuronas por unidad de longitud en x0 , cuyos axones establecen conexión sináptica en la columna x y
dichos axones se encuentran en el n-ésimo sistema de fibras, S±n (x) representa
la densidad de sinapsis (1/cm2 ) en la columna x, producida por neuronas
cuyos axones se encuentran en el n-ésimo sistema de fibras. De (2.6) se tiene
n
que R±
se puede escribir como el producto de tres factores
n
R±
(x0 , x, v) = tn± (x0 , x)S±n (x)F±n (v) .
(2.7)
Ahora, haremos la siguiente suposición
Hipótesis 3
La densidad de fibras representadas por la función tn± y que conectan
neuronas de la columna x0 con la columna x es igual a la densidad de fibras
que conectan neuronas de la columna x con x0 . Esta hipótesis es fundamental
para la simplificación del modelo (2.1), (2.2) y significa que
tn± (x0 , x) = tn± (x, x0 ) .
(2.8)
Ya hemos comentado que la función tn± (x0 , x) decrece exponencialmente
con la distancia entre x0 , y x. Para que se cumpla la propiedad (2.8) supondremos además que:
Hipótesis 4
La función tn± (x0 , x) depende únicamente de la distancia entre los puntos
x y x en la forma:
0
λn± −λn± |x0 −x|
e
,
(2.9)
2
donde (λn± )−1 (cm) es la escala longitudinal que determina el decrecimiento
de conexiones entre dos puntos de la corteza.
tn± (x0 , x) =
47
De (2.4) a (2.9) obtenemos
R± (x0 , x, v) =
N
1X
n 0
S±n (x)λn± F±n (v)e−λ± |x −x| .
2 n=1
(2.10)
Hipótesis 5
En primera aproximación supondremos que cada sistema de fibras está
n
caracterizado por una velocidad v±
de propagación de potenciales de acción.
En este caso:
n
).
F±n (v) = δ(v − v±
(2.11)
N
1X
0
n −λn
S±n (x)λn± δ(v − v±
)e ± |x −x| .
2 n=1
(2.12)
Finalmente obtenemos
R± (x0 , x, v) =
Pasemos ahora a imponer las restricciones que exigiremos a la función
g(x0 , t).
Hipótesis 6
Los potenciales de acción que se generan en la columna x0 en el instante
t dependen de la densidad de actividad sináptica - potenciales postsinapticos - y de la intensidad de la actividad neuroquı́mica (concentraciones de
neurotransmisores) que se realice en dicha columna en ese instante.
En este trabajo sólo consideraremos la dependencia de g con respecto a
h+ y h− , es decir
g(x0 , t) = G(h+ (x0 , t), h− (x0 , t)) ,
(2.13)
donde la función G es suave y tiene un comportamiento sigmoidal cuya forma
exacta es muy dı́ficil de predecir aunque esto no es importante para el análisis
cualitativo que llevaremos a cabo posteriormente.
48
Es decir, la función G(h+ , h− ) es creciente con respecto a h+ para cada
h− fijo y decreciente con respecto a h− para cada h+ fijo.
Además existen curvas de nivel
γ(h+ , h− ) = α, γ(h+ , h− ) = β ,
(2.14)
tales que


0,






si γ(h+ , h− ) < α ,
G(h+ , h− ) =  0 < G < 1 , si α < γ(h+ , h− ) < β ,





 1,
(2.15)
si γ(h+ , h− ) > β .
Para cada h− fijo (respectivamente h+ fijo) las funciones G(·, h− ) (resp.
G(h+ , ·)) deben tener un comportamiento sigmoidal usual.
Más adelante propondremos una expresión explı́cita para las curvas de
nivel (2.14).
Las relaciones sigmoidales son muy frecuentes en los modelos del sistema
nervioso. Por ejemplo, también existe una relación de este tipo entre las
concentraciones de neurotransmisores y las frecuencias de disparo de las neuronas o también la frecuencia de actividad sináptica.
Todas las hipótesis anteriores nos conducen a una simplificación del modelo determinado por las ecuaciones (2.1) y (2.2)
h+ (x, t) = ho+ (x, t)
+ 12
PN
n=1
S+n (x)λn+
R∞
−∞
n
0
e−λ+ |x −x| G(h+ (x0 , t −
|x0 −x|
), h− (x0 , t
n
v±
−
|x0 −x|
))dx0
vn ±
,
(2.16)
h− (x, t) = ho− (x, t)
+ 12
PN
n=1
S−n (x)λn−
R x+² −λn |x0 −x|
0|
|x0 −x|
0
−
G(h+ (x0 , t − |x−x
),
h
(x
,
t
−
))dx0 ,
n
−
x−² e
v
vn
−
49
−
(2.17)
donde le parámetro ² > 0 corresponde al carácter local de la actividad inhibitoria.
Sustituyendo (2.13) en (2.16) y (2.17) obtenemos
h+ (x, t) = ho+ (x, t)
+
1
2
R
0
−λn
n ∞
n
+ |x −x| g(x0 , t
n=1 S+ (x)λ+ −∞ e
PN
−
|x0 −x|
)dx0
n
v+
(2.18)
,
h− (x, t) = ho− (x, t)
+
1
2
PN
n
n
n=1 S− (x)λ−
R x+² −λn |x0 −x|
0|
−
g(x0 , t − |x−x
)dx0 .
x−² e
vn
(2.19)
−
En el sistema (2.18), (2.19), existen un conjunto de simplificaciones no
mencionadas explı́citamente con anterioridad. Consideramos que es importante explicar estas simplificaciones, ya que además de provocarnos algunas
confusiones, el considerarlas nos permitirá la construcción de un modelo más
detallado, el cual también consistirá de un sistema de ecuaciones integrales.
Este sistema se podrá transformar en un sistema de ecuaciones diferenciales
parciales de tipo hiperbólico en donde podremos analizar con mayor detalle
las distintas “fuerzas” que actúan en el medio cortical. Es pertinente mencionar que al intentar transformar el sistema (2.18), (2.19), en un sistema de
ecuaciones diferenciales hiperbólicas nos surgieron algunas complicaciones
que no permitieron dicha transformación.
Algunas de las simplificaciones introducidas en el sistema (2.18), (2.19) y
no mencionadas por Nuñez son:
No distingue por capas la densidad de actividad sináptica ni la generación
de los potenciales de acción; este hecho nos obliga a considerar que el efecto
producido en la actividad sináptica en una columna x en un instante t, producido por los potenciales de acción que se generan en tiempos anteriores en
las demás columnas, se distribuye de manera uniforme sobre cada capa de
la columna x (en cada capa el efecto total se dividirá por N , donde N es
el número de capas). A nivel de columna tampoco distingue la contribución
de la actividad sináptica de cada capa en la columna x a la generación de
50
los potenciales de acción en esa columna. Otro problema de estos modelos
es que tampoco distinguen la actividad sináptica producida por neuronas de
diferente clase.
En la sección siguiente presentaremos otros modelos que subsanan algunas
de estas simplificaciones.
2.4 ALGUNOS MODELOS ALTERNATIVOS AL MODELO
DE NUÑEZ
En esta sección analizaremos algunos modelos alternativos al modelo de
Nuñez. Estos modelos subsanan algunas de las simplificaciones mencionadas
en la sección anterior. También analizaremos brevemente la relación de estos
modelos con el modelo de Nuñez y las relaciones que tienen entre sı́.
1. El modelo más general de capas serı́a considerar todas las interacciones
que se pueden establecer entre todas las capas de cada par de columnas,
ası́ obtendrı́amos el sistema de 2N 2 ecuaciones integrales siguiente para
un modelo de N capas.
0j
hij
± (x, t) = h± (x, t)
+
R∞R∞
0
ij
0
ij
0
−∞ R± (x , x, v)g (x , t −
|x−x0 |
)dx0 dv
vij
(2.20)
,
donde:
hij
± (x, t) − La densidad de actividad sináptica en la j-ésima capa de
la columna x en el tiempo t producida por potenciales de acción
que se generan en las otras columnas en tiempos anteriores y que
viajan por el i-ésimo sistema de fibras, junto con la provocada por
los potenciales de acción que se disparan en tiempos anteriores en
estructuras interiores del cerebro.
h0j
± (x, t) − La densidad de actividad sináptica en la j-ésima capa de
la columna x producida por potenciales de acción que se generan
51
en estructuras internas del cerebro y que tienen conexión sináptica
con la capa j-ésima de la columna x.
ij
R±
(x0 , x, v) − El número de conexiones sinápticas/cm2 en la j-ésima
capa de la columna x por unidad de longitud en la capa i-ésima
de la columna x0 y por unidad de velocidad.
g ij (x0 , t) − La fracción del total de neuronas de la capa i-ésima de
la columna x0 que generan potenciales de acción en el tiempo t y
que tienen contacto sináptico con neuronas de la capa j-ésima de
las otras columnas.
Con respecto al modelo de Nuñez se tienen las siguientes relaciones
R± (x, t) =
N X
N
X
ij
(x, t) ,
R±
(2.21)
i=1 j=1
gij (x, t) = αij (x)g(x, t) ,
h± (x, t) =
N X
N
X
ij
h± (x, t) ,
(2.22)
(2.23)
i=1 j=1
donde αij (x) es un coeficiente entre 0 y 1 que indica la parte de g(x, t)
que se disparan en la capa i-ésima de la columna x y que llegan a la
capa j-ésima de las otras columnas.
2. Otra alternativa serı́a considerar
hj± (x, t) = h0j
± (x, t)
+
1
N
R∞R∞
0
−∞
j
R±
(x0 , x, v)
donde:
52
PN
i=1
αij g i (x0 , t −
|x0 −x|
)dx0 dv
i
v±
,
(2.24)
hj± (x, t) − La densidad total de actividad sináptica en la j-ésima
capa de la columna x en el tiempo t.
g j (x0 , t) − La fracción del total de neuronas de la capa j-ésima de la
columna x0 que generan potenciales de acción en el tiempo t.
h0j
± (x, t) − La densidad de actividad sináptica en la capa j-ésima de
la columna x producida por potenciales de acción que se generan
en tiempos anteriores en estructuras internas del cerebro.
j
R±
(x0 , x, v) − El número de conexiones sinápticas/cm2 en la columna
x que llegan a la capa j-ésima de la columna x, por unidad de longitud en la columna x0 y por unidad de velocidad.
αij − Los pesos que corresponden a la contribución de la capa iésima de la columna x0 a la actividad sináptica en la capa j-ésima
de la columna x. Obviamente 0 ≤ αij ≤ 1 y podrı́a considerarse
P
αij = αij (x, x0 , t). También se tiene que N
j=1 αij ≤ 1.
Con respecto al modelo 1 se está eliminando la contribución particular
de cada capa de la columna x0 a la actividad sináptica en la capa
j-ésima de la columna x, además de que se tiene un problema similar
que en el sistema (2.18) (2.19) ya que no se distingue tampoco el efecto
de g i (x0 , t) en cada capa de la columna x, asumiendo que se distribuye
uniformemente sobre cada capa de x (la contribución a la actividad
0 0
sináptica en la capa j-ésima de la columna x serı́a g (xN ,t) ). La ventaja
que se tiene en este modelo es que g j (x0 , t) depende únicamente de
hj+ (x0 , t) y hj− (x0 , t); es decir de la actividad sináptica en esa capa y no
de la actividad sináptica en todas las capas de esa columna.
Se puede considerar que se cumple la relación
h± (x, t) =
PN
i=1
βi (x)hi± (x, t) ,
donde
βi (x) son ciertos pesos positivos tales que
PN
i=1
βi (x) = 1.
3. También se puede considerar el siguiente sistema
53
(2.25)
h̄j± (x, t) = h0± (x, t)
+
R∞R∞
0
j
0
j
0
−∞ R± (x , x, v)g (x , t
−
|x−x0 |
)dx0 dv
v
(2.26)
,
donde
h̄j± (x, t) − La densidad de actividad sináptica en la columna x en el
tiempo t producida por potenciales de acción que se generan en la
capa j-ésima de las otras columnas en tiempos anteriores más la
actividad sináptica en la columna x producida por potenciales de
acción que se generan en tiempos anteriores en otras estructuras
interiores del cerebro. En este caso no se distingue, como en el
caso 2, la actividad sináptica en la columna por capas, sino por
los potenciales de acción que provienen de cada capa en otras
columnas.
ḡ j (x0 , t) − La fracción del total de neuronas de la columna x0 que
generan potenciales de acción en el tiempo t y que viajan por las
fibras del j-ésimo sistema.
h0± (x, t) − La densidad de actividad sináptica en la columna x en el
tiempo t producida por los potenciales de acción que se generan
en tiempos anteriores en otras estructuras interiores del cerebro.
j
(x0 , x, v) − El número de conexiones sinápticas por cm2 en la
R±
columna x por unidad de longitud en la capa j-ésima de la columna
x0 y por unidad de velocidad.
El problema que se presenta en este modelo es que ḡ j (x, t) depende
de toda la actividad sináptica que ocurre en la columna x, es decir de
h̄1± (x, t), . . . , h̄N
± (x, t)(x, t).
Estos tres modelos y el modelo de Nuñez son las posibilidades más importantes que se pueden encontrar en un modelo de N capas. Estos 4 modelos
están relacionados. El modelo 1 es el modelo más detallado de los cuatro,
pero es el que incluye más ecuaciones y parámetros que nos dificultarı́an el
54
estudio del mismo. Estos modelos están relacionados entre sı́ de la siguiente
manera:
En el modelo 2 se está considerando la densidad de actividad sináptica
por capas en cada columna pero sin distinguir la contribución a esta actividad
de las distintas capas de las otras columnas, además tampoco se distingue
el efecto de una capa de una columna (g j (x, t)) sobre cada capa de las otras
columnas, se considera que el efecto (g j (x, t)) se distribuye de manera uniforme sobre las capas de cada columna afectada. La relación entre el modelo
2 y el 1 se da a través de considerar que
Ã
g ij
|x − x0 |
x0 , t −
vij
!
Ã
= αij g i
|x − x0 |
x0 , t −
vi
!
,
(2.27)
i.e. el efecto de la i-ésima capa de la columna x sobre la j-ésima capa de las
columnas afectadas es independiente de la capa que afecta y se distribuye de
manera “uniforme” sobre cada capa. También se cumple que
N
X
1 i 0
ij
R± (x , x, v) =
(x0 , x, v)
R±
N
j=1
(2.28)
y por último
hi± (x, t) =
N
X
ij
h± (x, t) .
(2.29)
hi± (x, t) ,
(2.30)
j=1
En el modelo 3 se tiene obviamente que
h± (x, t) =
N
X
i=1
donde h+ (x, t) es la densidad total de actividad sináptica que aparece en
el modelo de Nuñez. Nótese, además que g j (x0 , t) = ḡ j (x0 , t), pero g j sólo
depende de hj+ y h− con respecto a las variables h̄i+ , h− se tiene que ḡ j ;
depende de h̄1± , . . . , h̄N
± . Con más exactitud
j
j
j
ḡ = Ḡ (h+ , h− ) = Ḡ
ÃN
X
i=1
55
h̄i+ , h−
N
X
i=1
!
h̄i−
(2.31)
y como el modelo de Nuñez se relaciona con el modelo 2 a través de
h± (x, t) =
N
X
j
h± (x, t) ,
(2.32)
j=1
g j (x, t) =
g(x, t)
.
N
(2.33)
Entonces también el modelo (2.18) (2.19) se relaciona con el modelo 1.
El modelo que eligiremos para estudiar será el modelo 3 ya que es un poco
más detallado que el de Nuñez y resulta un modelo con menos ecuaciones y
parámetros que los otros dos modelos restantes.
En la siguiente sección desarrollaremos el modelo 3.
2.5 DESARROLLO Y ESTUDIO DE UN MODELO DE LA
ACTIVIDAD SINÁPTICA DE LA CORTEZA
CEREBRAL
En esta sección desarrollaremos el modelo 3, que como mencionamos en
la sección anterior es un poco más preciso que el del sistema (2.18), (2.19),
además este sistema se puede transformar en un sistema de ecuaciones diferenciales parciales. El modelo tiene también una conexión muy clara en el
sistema (2.18), (2.19). En este modelo distinguiremos la contribución de cada
capa de cada columna a la actividad sináptica de una columna de corteza
dada en un instante de tiempo y lo que estudiaremos es la evolución espaciotemporal de esa actividad sináptica tomando en cuenta la distribución de
conexiones sinápticas, las velocidades de propagación de los potenciales de
acción de las neuronas de la capa correspondiente. Para ser más precisos
introduciremos los elementos de que consiste el modelo.
Sean:
h̄j± (x, t) − La densidad de actividad sináptica (potenciales postsinápticos)
en la columna x en el tiempo t, número de sinapsis activas/cm2 , producida por potenciales de acción que se generan en las otras columnas y
que viajan por las fibras del sistema j-ésimo, más la actividad sináptica
56
en la columna x producida por potenciales de acción que se generan en
tiempos anteriores en otras estructuras interiores del cerebro.
ḡ j (x, t) − La fracción del total de neuronas de la columna x0 que generan
potenciales de acción en el tiempo t y que viajan por las fibras del
j-ésimo sistema.
h̄0± (x, t) − La densidad de actividad sináptica en la columna x en el tiempo
t producida por los potenciales de acción que se generan en estructuras interiores del cerebro en tiempos anteriores y que tienen conexión
sináptica en la columna x.
j
R±
(x0 , x, v) − El número de conexiones sinápticas/cm2 , del j-ésimo sistema en la columna x por unidad de longitud en x0 y por unidad de
velocidad.
La relación entre h̄j± (x, t), ḡ j (x, t) y h± (x, t), g(x, t) es la siguiente
h± (x, t) =
N
X
j
h̄± (x, t) .
(2.34)
j=1
La misma idea que se usó para obtener la ecuación (2.18) se puede aplicar
para obtener un sistema de ecuaciones integrales para h̄j± (x, t). Ası́ que para
j = 1, . . . , N
h̄j+ (x, t) = h̄0+ (x, t)
+
R ∞ −λj |x−x0 | j 0
1 j j
e +
ḡ (x , t
λ S (x) −∞
2 + +
−
|x0 −x|
)dx0
j
v+
(2.35)
.
Para h̄j− (x, t) se tiene
h̄j− (x, t)
=
h̄0− (x, t) +
Z x+²
j
1 j
|x − x0 |
0
j
S− (x)λ−
)dx0 . (2.36)
e−λ− |x−x | ḡ j (x0 , t −
j
2
x−²
v−
57
Ahora, sabemos que los potenciales de acción que se generan en cada capa
de cada columna (ḡ j (x, t)) dependen de la actividad sináptica que existe en
esa columna (no depende únicamente de la actividad sináptica en esa capa
por la forma en que se definió h̄j± (x, t)) y de la actividad neuroquı́mica en la
misma columna.
En este trabajo sólo consideraremos la dependencia de ḡ j (x, t) con respecto a (h̄j± (x, t), j = 1, . . . , N ) es decir consideraremos
N
X
ḡ j (x, t) = Ḡj (
h̄i+ (x, t),
N
X
hi− (x, t)) .
(2.37)
i=1
i=1
La forma exacta de esta dependencia es dı́ficil de predecir, pero cualitativamente consideramos que debe ser de tipo sigmoide en cada variable. Ası́
para cada k = 1, . . . , N y para cada h̄i+ (x, t), i 6= k y h̄j− (x, t) fijo debemos
k
N
1
N
tener que Ḡj (h̄1+ , ..., h̄k−1
+ , h̄+ , . . . , h+ , h̄− , . . . , h̄− ) es una función creciente de
tipo sigmoide.
Para cada 1 ≤ k ≤ N y para cada h̄1+ , . . . , h̄N
+ fijo debemos tener que
j
1
N
Ḡ (h̄+ , . . . , h̄+ , ·) es una función decreciente de tipo sigmoide.
Además existen curvas de nivel
j
N
1
γj (h̄1+ , ..., h̄j+ , . . . , hN
+ , h− ) = αj y γj (h̄+ , ..., h̄+ , . . . , h+ , h− ) = βj , (2.38)
tales que
Ḡj (h̄1+ , h̄j+ . . . h̄N
+ , h− ) =


0,













si γj (h̄1+ , h̄j+ . . . h̄N
+ , h− ) < αj ,
0 < Ḡj < 1 ,
si αj < γj (h̄1+ , h̄j+ . . . h̄N
+ , h − ) < βj ,
1,
si γj (h̄1+ , h̄j+ . . . h̄N
+ , h− ) > βj .
(2.39)
Más adelante propondremos una expresión explı́cita para las curvas de nivel
(2.38).
Con estas condiciones el sistema (2.35), (2.36) se transforma en
58
h̄j+ (x, t) = h̄0+ (x, t)
+ 12 λj+ S+j (x)
R∞
−∞
j
0
PN
e−λ+ |x −x| Ḡj (
i=1
h̄i+ (x0 , t −
|x0 −x| PN
), i=1
j
v+
hi− (x0 , t −
|x0 −x|
)dx0
i
v+
,
(2.40)
h̄j− (x, t) = h̄0− (x, t)
+ 12 λj− S−j (x)
R x+² −λj |x0 −x| j PN
0
P
|x0 −x|
0
i
0
−
Ḡ ( i=1 h̄i+ (x0 , t − |xv−x|
), N
j
i=1 h̄− (x , t − v i )dx ) .
x−² e
−
−
(2.41)
Como ya hemos mencionado anteriormente consideramos más adecuado
para el análisis del modelo, transformar el sistema (2.35), (2.36) en un sistema
de ecuaciones diferenciales parciales para h̄j+ (x, t), h̄j− (x, t). Esto último lo
haremos en la siguiente sección.
59
2.6 TRANSFORMACIÓN DEL SISTEMA DE ECUACIONES
INTEGRALES EN UN SISTEMA DE ECUACIONES
DIFERENCIALES PARCIALES
Usando la regla de Leibnitz para derivar bajo el signo de la integral,
∂ 2 h̄j
∂ h̄j
∂ h̄j
∂ 2 h̄j
podemos obtener ∂x± , ∂x2± , ∂t± , ∂t2± , no escribiremos estas derivadas ya que
además de constar de muchos términos, no son necesarias para los siguientes
desarrollos. Escribiremos solamente el sistema de ecuaciones diferenciales
parciales resultante, tanto para h̄j+ como para h̄j− .
Para h̄j+ se obtiene
∂ 2 h̄j+
∂t2
j
+ 2λj+ v+
·
j 2
)
(v+
∂ 2 h̄j+
∂x2
∂ h̄j+
∂t
j 2 j
+ (λj+ v+
) h̄+ =
j 0
) (x) ∂ h̄j+
2(S+
j
∂x
S (x)
−
µ
j 00
j
j 0
) (x)
(x)(S+
) (x))2 −S+
2((S+
+
j 2
(S+
) (x)
+
h
¶
j j
j j
ḡ (x, t) +
S+ (x) λj+ v+
+λj+ v+
i
∂ḡ j
(x, t)
∂t
¸
h̄j+
(2.42)
+ F+j (x, t) .
Para h̄j− se obtiene
∂ 2 h̄j−
∂t2
j
+ 2λj− v−
·
j 2
)
(v−
∂ 2 h̄j−
∂x2
−
∂ h̄j−
∂t
j 2 j
) h̄− =
+ (λj− v−
j 0
) (x) ∂ h̄j−
2(S−
j
∂x
S (x)
µ
j 00
j
j 0
) (x)
(x)(S−
) (x))2 −S−
2((S−
+
j 2
(S−
) (x)
−
h
j j
j j
ḡ (x, t) +
S− (x) λj+ v−
+λj− v−
i
∂ḡ j
(x, t)
∂t
¸
¶
h̄j−
(2.43)
+ F−j (x, t) + F j (x, t, ²) ,
donde
·
F+j (x, t)
=
∂ 2 h̄0+
∂t2
−
j 2
(v+
)
+
·
0
j ∂ h̄+
2λj+ v+
∂t
∂ 2 h̄0+
∂x2
−
¸
+
j 2 0
(λj+ v+
) h̄+
j 0
2(S+
) (x) ∂ h̄0+
j
∂x
S (x)
+
60
µ
+
j 0
j
j 00
2((S+
) (x))2 −S+
(x)(S+
) (x)
j 2
(S+
) (x)
¸
¶
h̄0+
,
(2.44)
·
F−j (x, t) =
∂ 2 h̄0−
∂t2
−
j 2
(v−
)
j
+ 2λj− v−
·
∂ 2 h̄0−
∂x2
−
∂ h̄0−
∂t
¸
j 2 0
+ (λj− v−
) h̄−
j 0
2(S−
) (x) ∂ h̄0−
j
∂x
S (x)
µ
+
j 0
j
j 00
2((S−
) (x))2 −S−
(x)(S−
) (x)
j 2
(S−
) (x)
−
·
j
j j
j j −λ− ² j
F j (x, t, ²) = − 21 λj− v−
S− (x) v−
λ− e
ḡ (x − ², t −
j
j
+ e−λ− ² ∂ḡ
(x − ², t −
∂t
²
j )
v−
j
j −λ− ² j
ḡ (x + ², t −
e
+ λj− v−
+
j
j
e−λ− ² ∂ḡ
(x
∂t
+ ², t −
²
j )
v−
j
¶
¸
h̄0−
,
(2.45)
²
j )
v−
j
j −λ− ² ∂ḡ
− v−
e
(x − ², t −
∂x
²
j )
v−
²
j )
v−
+
j
v−
j
j
e−λ− ² ∂ḡ
(x
∂x
¸
+ ², t −
²
j )
v−
.
(2.46)
Con la finalidad de presentar concisamente y más adecuadamente para
nuestros cálculos el sistema (2.42), (2.43), definiremos 4N operadores diferenciales parciales, actuando sobre las funciones de dos variables con derivadas
parciales hasta de segundo orden continuas. Denotaremos estos operadores
con los sı́mbolos Li,j
± , 1 ≤ i ≤ 2, 1 ≤ j ≤ N y estarán definidos como sigue
L1,j
± [u] =
·
L2,j
± [u]
=
j 2
)
(v±
∂2u
∂x2
−
∂2u
∂t2
j ∂u
+ 2λj± v±
+ (λ± v± )2 u ,
∂t
j 0
2(S±
) (x) ∂u
j
S (x) ∂x
±
µ
+
j 0
j
j 00
2((S±
) (x))2 −S±
(x)(S±
) (x)
j 2
(S±
) (x)
(2.47)
¶ ¸
u ,
(2.48)
donde u = u(x, t) es una función con derivadas parciales hasta de segundo
orden continuas.
61
Con esta nueva notación, el sistema (2.42), (2.43) se convierte en el sistema
j
2,j j
L1,j
+ [h̄+ ] = L+ [h̄+ ]
h
j j
j j
+ λj+ v+
S+ (x) λj+ v+
ḡ (x, t) +
i
∂ḡ j
(x, t)
∂t
(2.49)
+ F+j (x, t) ,
j
2,j j
L1,j
− [h̄− ] = L− [h̄− ]
h
j j
j j
+ λj− v−
S− (x) λj− v−
ḡ (x, t) +
i
∂ḡ j
(x, t)
∂t
(2.50)
+ F−j (x, t) + F (x, t, ²) .
Para reducir el sistema a un sistema equivalente más sencillo de estudiar
haremos uso de los siguientes resultados los cuales son muy fáciles de probar
j
Lema 1 Sean Li,j
± los operadores diferenciales definidos en (2.48) y sean L̄±
los operadores diferenciales parciales definidos por
L̄j± [u] =
∂2u
∂t2
2
j ∂ u
j ∂u
j 2
) u − v±
+ 2λj± v±
+ (λj± v±
.
∂t
∂x2
(2.51)
Entonces
L̄j± [ūj± ] =
1
j
S±
(x)
h
i
2,j j
j
L1,j
± [u± ] − L± [u± ] ,
(2.52)
donde
ūj± (x, t) =
uj± (x, t)
.
S±j (x)
Lema 2 Sea L̄j± como en el Lema 1, entonces
j
j
j
j
eλ+ v+ t L̄j+ [uj+ ] = ∂+2,j [ūj+ ] ,
eλ− v− t L̄j− [uj− ] = ∂−2,j [ūj− ] ,
62
(2.53)
(2.54)
donde
j
j
ūj± (x, t) = eλ± v± t uj± (x, t)
y
∂+2,j [u] =
2
∂ 2u
j 2∂ u
−
(v
)
,
+
∂t2
∂x2
(2.55)
∂−2,j [u] =
2
∂ 2u
j 2∂ u
−
(v
)
.
−
∂t2
∂x2
(2.56)
Aplicando el Lema 1 al sistema (2.49), (2.50) ese sistema es equivalente
al sistema
"
L̄j+ [H+j ]
"
L̄j− [H−j ]
#
∂ḡ j
1
j j
j
ḡ (x, t) +
λj+ v+
= j
F+j (x, t) + λj+ v+
(x, t) ,
∂t
S+ (x)
(2.57)
#
∂ḡ j
1
1
j j
j
ḡ (x, t) +
λj− v−
= j
F−j (x, t)+λj− v−
,
(x, t) +F j (x, t, ²) j
∂t
S− (x)
S− (x)
(2.58)
donde
H+j (x, t)
h̄j± (x, t)
= j
S± (x)
y aplicando el Lema 2 al sistema (2.57), (2.58), este sistema es equivalente a
"
∂+2,j [H̄+j ]
j
λj+ v+
t
=e
"
∂−2,j [H̄−j ]
j
λj− v−
t
=e
##
"
F+j (x, t)
∂ḡ j
j j
j j j
(x, t)
+
λ
v
λ
v
ḡ
(x,
t)
+
+
+
+
+
∂t
S+j (x)
#
"
,
(2.59)
#
F j (x, t, ²)
∂ḡ j
F−j (x, t)
j j
j j j
(x,
t)
+
+
λ
,
v
λ
v
ḡ
(x,
t)
+
−
−
−
−
∂t
S−j (x)
S−j (x)
(2.60)
63
donde
j
j
H̄±j (x, t) = eλ± v± t H±j (x, t) .
El sistema (2.57), (2.58) es un sistema de 4N ecuaciones diferenciales parP
i PN
i
ciales no lineales acopladas (recordemos que ḡ j (x, t) = Ḡj ( N
i=1 h̄+ ,
i=1 h− ),
j
donde lo único que le hemos pedido hasta ahora a Ḡ es que sea una función
sigmoide de cada una de las variables h̄i± ), además, la ecuación (2.59) depende de un pequeño parámetro ² > 0. Este sistema es sumamente complejo
de estudiar, ası́ que para hacer más tratable el sistema, estudiaremos primero
el caso N = 1, esto significa que se tiene únicamente una sola capa en cada
columna. En este caso el sistema a estudiar es
"
#
1
∂ḡ
L̄+ [H+ ] =
F+ (x, t) + λ+ v+ λ+ v+ ḡ(x, t) + (x, t) ,
S+ (x)
∂t
"
(2.61)
#
1
∂ḡ
1
L̄− [H− ] =
F− (x, t) + λ− v− λ− v− ḡ(x, t) + (x, t) + F (x, t, ²)
.
S− (x)
∂t
S− (x)
(2.62)
Es decir, tenemos dos ecuaciones diferenciales parciales con la ecuación
(2.62) dependiendo de un parámetro pequeño ² > 0, la función ḡ depende
únicamente de h̄± , de manera sigmoide en cada una de ellas
ḡ(x, t) = Ḡ(h+ , h− ) .
(2.63)
En la siguiente sección estudiaremos el sistema (2.61), (2.62) usando
técnicas de parámetro pequeño, es decir, estudiaremos primeramente el sistema degenerado, − el sistema (2.61), (2.62) con ² = 0 − , para después, con
base en este estudio, obtener información sobre el sistema no degenerado.
Además introduciremos una expresión para Ḡ, que dependerá de una sola
variable, la cual es una combinación lineal de h̄+ , y h̄− , a esta nueva variable
la llamaremos variable de activación, ya que el hecho de que se generen o no
potenciales de acción en una columna dependen de esta nueva variable.
64
2.7 LA ECUACIÓN DE ACTIVACIÓN PARA LA CORTEZA
CEREBRAL
Sabemos que ḡ(x, t) depende de una manera muy compleja de la actividad
sináptica excitatoria h̄+ , y de la actividad sináptica inhibitoria h̄− , i.e. ḡ(x, t)
es igual a Ḡ(h̄+ , h̄− ), cualitativamente consideramos que Ḡ(., h̄− ) debe ser
una función sigmoide creciente de h̄+ , para cada h̄− fijo y que debe ser una
función sigmoide decreciente de h̄− para cada h̄+ fijo. Este planteamiento
es sumamente general y difı́cil de estudiar, ası́ que estudiaremos un caso
particular de la dependencia de ḡ con respecto a h̄+ y h̄− . El caso que
estudiaremos está basado en lo siguiente:
1. Si V± representa el valor del potencial postsináptico en una sinapsis,
entonces V± h̄± (x, t) representará la densidad de potencial postsináptico
en la columna x en el tiempo t.
2. Habrá activación en la columna x (generación de potenciales de acción)
cuando V+ H+ (x, t) − V− H− (x, t) ≥ V0 , donde V0 es un cierto umbral
promediado de disparo de la columna x, en caso contrario no habrá
activación.
3. Si V+ H+ (x, t) − V− H− (x, t) es muy grande, habrá saturación de la
columna x en el tiempo t y entonces, no aumentará significativamente
la actividad de la columna en ese tiempo.
Basados en estos tres puntos consideramos que:
ḡ(x, t) = Ḡ(V+ H+ (x, t) − V− H− (x, t)) ,
(2.64)
donde Ḡ es una función sigmoide de una sola variable. Ḡ es una función
sigmoide de la densidad de potencial postsináptico de la columna.
En realidad, como veremos después, se pueden elegir varias expresiones
cuantitativas para Ḡ.
Llamaremos a la variable u(x, t) = V+ H+ (x, t) − V− H− (x, t) variable de
activación de la columna x.
65
La dinámica espacio-temporal de u(x, t) nos permitirá estudiar las zonas
de la corteza que están activas en un tiempo determinado. Para cada t, las
zonas activas de la corteza estarán caracterizadas por el conjunto
E(t) = {x en la corteza|u(x, t) ≥ V0 } .
Para obtener la ecuación diferencial parcial que determina la dinámica
de la variable de activación u(x, t) multiplicamos la ecuación L̄+ [H+ ] =
1
F (x, t) + λ+ v+ [λ+ v+ ḡ(x, t) + ∂ḡ
(x, t)] por V+ y le restamos V− L̄+ [H− ]
S+ (x) +
∂t
en ambos lados de la ecuación, ası́ se obtiene
"
L̄+ [u] = V+ λ+ v+
#
∂ḡ
F+ (x, t)
λ+ v+ ḡ(x, t) + (x, t) +V+
−V− L̄+ [H− ] . (2.65)
∂t
S+ (x)
Además tenemos la ecuación para las inhibiciones
"
L̄− [H− ] = λ− v−
#
∂ḡ
F− (x, t) F (x, t, ²)
λ− v− ḡ(x, t) + (x, t) +
+
.
∂t
S− (x)
S− (x)
(2.66)
Como mencionamos en la sección anterior, primeramente estudiaremos el
caso degenerado (² = 0) del sistema (2.65), (2.66)
"
L̄+ [u] = V+ λ+ v+
#
∂ḡ
F+ (x, t)
λ+ v+ ḡ(x, t) + (x, t) +V+
−V− L̄+ [H− ] . (2.67)
∂t
S+ (x)
Además tenemos la ecuación para las inhibiciones
L̄− [H− ] =
F− (x, t)
.
S− (x)
(2.68)
Usando el lema 2, se puede reducir el sistema (2.67), (2.68) al estudio del
sistema equivalente
"
L̄+ [u] = V+ λ+ v+
#
∂ḡ
∂u
F+ (x, t)
λ+ v+ ḡ(u) +
(u)
+ V+
− V− L̄+ [H− ] , (2.69)
∂u
∂t
S+ (x)
66
∂−2 (H̄− ) = e−λ− v− t
F− (x, t)
,
S− (x)
(2.70)
donde
H̄(x, t) = eλ− v− t H− (x, t) .
Utilizando la definición de L̄, ∂−2 y calculando L̄+ [e−λ− v− t H̄− (x, t)] y reordenando los términos, obtenemos el sistema
∂2u
∂t2
h
2
2 ∂ u
+ 2λ+ v+ ∂u
+ (λ+ v+ )2 u − v+
= V+ λ+ v+ λ+ v+ ḡ(u) + ḡ 0 (u) ∂u
∂t
∂x2
∂t
·
−V− e−λ− v− t (1 −
+ S+V+(x) F+ (x, t) +
2
v+
∂ 2 H̄−
2 ) ∂t2
v−
i
¸
−
+ 2(λ+ v+ − λ− v− ) ∂ H̄
+ (λ+ v+ − λ− v− )2 H̄−
∂t
1 F− (x,t) −λ− v− t
2 S (x) e
v−
−
,
(2.71)
2
∂ 2 H̄−
e−λ− v− t
2 ∂ H̄−
=
v
+
F− (x, t) .
−
∂t2
∂x2
S− (x)
(2.72)
Nótese que la ecuación (2.72) es independiente de ḡ, ası́ que este sistema
se puede interpretar como un proceso global de propagación de ondas de la
variable de activación, modulado por un proceso de propagación de ondas de
las inhibiciones.
67
2.8 CONSTRUCCIÓN DE UN MODELO
COMPARTIMENTADO DE ACTIVACIÓN DE LA
CORTEZA CEREBRAL.
En el modelo que hemos analizado en las secciones anteriores –con las
diversas formas propuestas para la distribución de conexiones– se eliminan
los efectos producidos en cada columna por columnas lejanas a ella. La eliminación de los “efectos lejanos”, se debe a la elección de una distribución
de conexiones que decrece exponencialmente con las distancias entre columnas con la misma escala de decrecimiento, en todas las direcciones sobre la
corteza. Una manera de tomar en cuenta estos “efectos lejanos”, es considerar
a la corteza cerebral dividida en n compartimentos disjuntos interconectados,
considerando la distribución de conexiones para cada par de compartimentos
de decrecimiento exponencial, pero con una escala de decrecimiento que depende del par de compartimentos considerados. Ası́, “los efectos” producidos
en una columna de cierto compartimento estarán clasificados dependiendo de
los compartimentos con que interactúa. En nuestro modelo consideramos que
las velocidades de propagación de los potenciales de acción no dependen de
los compartimentos que se están considerando. La división de la corteza
en n compartimentos, además de ser una mejor aproximación de la corteza
real que la propuesta en las secciones anteriores, nos permite correlacionar
nuestros resultados teóricos con las mediciones reales de actividad cortical.
Además nos permitirá estudiar el diagrama de conexiones funcional que se
establece en un estado fisiológico dado. Se podrı́a estudiar, también la interrelación que existe entre las distintas áreas de la corteza, siguiendo la división
propuesta por el anatomista K. Brodmann (Fig. # 9).
A continuación describiremos más detalladamente este modelo.
Como ya hemos mencionado, en el presente modelo consideramos a la
corteza cerebral dividida en n compartimentos ajenos Ωi , i = 1, . . . , n. Por
simplicidad y siguiendo las ideas de las secciones anteriores, consideraremos
una banda de corteza en la cual, las ondas que se propagan por ella, tienen
longitud de onda muy pequeña con respecto a la longitud de la misma. Con
estas simplificaciones puede pensarse en la corteza como una recta dividida
en n intervalos disjuntos, identificando los dos intervalos infinitos como uno
solo (considerando que los “efectos ” de los dos intervalos infinitos provocan
68
el mismo efecto sobre los otros intervalos).
Lo que deseamos es relacionar la densidad de actividad sináptica – excitatoria e inhibitoria – en una columna de corteza con los potenciales de acción
que se generan en las otras columnas en tiempos anteriores, pero como tenemos dividida a la corteza en n compartimentos, aumentará el número de
ecuaciones a considerar.
A continuación procederemos a definir las variables y parámetros que
intervienen en el modelo, deduciremos las ecuaciones en forma integral y al
igual que en las secciones anteriores transformaremos éstas en un sistema de
ecuaciones diferenciales parciales.
Como en las secciones anteriores, el subı́ndice + se referirá a las sinapsis
excitatorias y el − a las sinapsis inhibitorias.
Sean :
a) hij
± (x, t) − La densidad de sinapsis activas o de actividad sináptica por
cm2 en la columna x ∈ Ωj , en el tiempo t, producidos por los potenciales
de acción que se generan en el compartimento Ωi en tiempos anteriores
hij
± (x, t) representa el efecto de Ωi en la columna x ∈ Ωj en el instante
t (1/cm2 ).
b) f±j (x, t) − La densidad de sinapsis activas por cm2 en la columna x ∈ Ωj
en el instante t, producidas por los potenciales de acción que se generan
en tiempos “anteriores” en las estructuras internas del cerebro, principalmente el tálamo, y que tienen conexión sináptica con la columna
x ∈ Ωj . Supondremos que f−j (x, t) ≡ 0 para cada j.
c)
hj± (x, t) =
n
X
ij
h± (x, t) + f±j (x, t) .
(2.73)
i=1
Representará la densidad de actividad sináptica total en la columna x ∈
Ωj en el instante t, producida por los potenciales de acción generados
en todas las regiones de la corteza y las estructuras internas del cerebro,
principalmente el tálamo.
ij
d) R±
(x0 , x, v) − El número de conexiones sinápticas por cm2 en la columna
x ∈ Ωj por unidad de longitud en x0 ∈ Ωi y por unidad de velocidad, es
ij
decir, R±
(x0 , x, v)dx0 dv representa el número de conexiones sinápticas
69
por cm2 en la columna x ∈ Ωj que provienen de las neuronas que
están ubicadas desde la columna x0 hasta la columna x0 + dx0 y por
cuyos axones los potenciales de acción “viajan” con velocidad entre v
y v + dv.
e) g(x, t) − La fracción del total de neuronas en la columna x que generan
potenciales de acción en el instante t.
f ) g j (x, t) − La fracción del total de neuronas en la columna x ∈ Ωj que
generan potenciales de acción en el tiempo t. Obviamente
g(x, t) = g j (x, t) para x ∈ Ωj .
(2.74)
Lo que queremos conocer es la evolución espacio-temporal de la densidad
de actividad sináptica total en cada columna de cada compartimento de la
corteza (hij
± (x, t), x ∈ Ωj , j = 1, . . . , n). Para lograr esto, primero obtendremos la evolución espacio-temporal de la densidad de actividad sináptica
en cada columna de un compartimento dado de la corteza producida por
la interacción con otro compartimento de la misma (hij
± (x, t), x ∈ Ωj , i, j =
1, . . . , n) y le añadiremos la producida por las estructuras internas del cerebro
(f+j (x, t), j = 1, . . . , n).
Consideraremos que por cada par de compartimentos Ωi , Ωj se cumplen
las siguientes hipótesis :
1. Existe una única velocidad de propagación, v± , de los potenciales de
acción que viajan de Ωi a Ωj y que no depende de Ωi ni de Ωj .
ij
(x0 , x, v)) depende únicamente de la dis2. La densidad de conexiones (R±
tancia entre las columnas y tienen su decrecimiento exponencial, con
cada escala de decrecimiento (λij
± ) que sólo depende de las regiones
Ωi , Ωj .
R± (x0 , x, v) =
n
X
ij
R±
(x0 , x, v),
i=1
donde
70
x ∈ Ωj ,
(2.75)
ij
R±
(x0 , x, v)
=
 ij
ij
0 −λ |x−x0 |

Sj (x)δ(v − v± ) ,
 α± (x )e ±


0,
si x ∈ Ωj , x0 ∈ Ωi
en otro caso
(2.76)
0
R± (x , x, v) representará el número de conexiones sinápticas por cm2 en
x, por unidad de longitud en x0 y por unidad de velocidad. Sj (x) es la
ij
densidad de conexiones sinápticas en x ∈ Ωj (1/cm2 ) y α±
(x0 ) (1/cm)
es un factor de normalización, sujeto a la condición
n Z
X
i=1 R
ij
0
ij
α±
(x0 )e−λ± |x−x | dx0 = 1, para cada x ∈ Ωj .
(2.77)
En lo que sigue, identificaremos Ωi (i = 1, . . . , n) con el conjunto de
ij
ij
ij
,
(x0 ) = α±
(x0 ) constante; α±
los números reales y si suponemos α±
entonces de (2.77) obtenemos :
2
ij
n
X
α±
i=1
λij
±
= 1.
(2.78)
Debido al carácter local de las inhibiciones, podemos suponer que
ij
= 0 si i 6= j.
α−
Ası́, en lo que sigue consideraremos que :
ij
=
α+
λij
λjj
+
jj
= − .
, α−
2n
2
(2.79)
3. Los potenciales de acción que se generan en la columna x ∈ Ωj en el
instante t, dependen de la actividad sináptica total en x ∈ Ωi en el instante t y de la densidad de la actividad neuroquı́mica (concentración
de neurotransmisores) que se realice en dicha columna en ese instante.
Nosotros sólo consideraremos la dependencia de g i con respecto a la actividad total en la columna x en el instante t. Además, de la suposición
de que la actividad inhibitoria en la corteza sólo tiene un carácter local,
se concluye que hki
− (x, t) ≡ 0 si k 6= i.
71
Finalmente como las conexiones aferentes a la corteza son sólo excitatorias, tenemos que f−j (x, t) ≡ 0 y ası́ obtenemos
g i (x, t) = g i (
n
X
i
ii
hki
+ (x, t) + f+ (x, t), h− (x, t)) .
(2.80)
k=1
Siguiendo el mismo razonamiento expuesto en las secciones anteriores,
tenemos que :
hij
− (x, t)
Z ∞
jj
λij
|x − x0 |
0
+
e−λ+ |x−x | g(x0 , t −
)dx0 ,
=
Sj (x)
2n
v+
−∞
hjj
− (x, t) =
Z x+²
jj
λjj
|x − x0 |
0
−
Sj (x)
)dx0 ,
e−λ+ |x−x | g(x0 , t −
2
v−
x−²
(2.81)
(2.82)
donde el parámetro ² > 0 corresponde a la localidad de la actividad inhibitoria. De nuevo, debido al carácter local del proceso inhibitorio, consideraremos
que hij
− (x, t) ≡ 0 si i 6= j.
Introduciendo las nuevas funciones incógnitas
H+ij (x, t)
nhij
(x, t)
hjj
+ (x, t)
jj
=
, H− (x, t) = −
Sj (x)
Sj (x)
(2.83)
2
y realizando el mismo proceso de aplicar el operador ∂ 2 a cada lado de
∂x
(2.81), (2.82) y después pasar al lı́mite cuando ² → 0, se obtiene, debido a
(2.74)
2 ij
∂H+ij
∂ 2 H+ij
ij
ij
ij 2 2
2 ∂ H+
v
)
v
H
−
v
+
2λ
+
(λ
+
+
+
+ +
+
∂t2
∂t
∂x2
"
=
λij
+ v+
#
ij
λij
+ v+ g (x, t)
∂g ij
(x, t) ,
+
∂t
x ∈ Ωj ,
2 jj
∂ 2 H−jj
∂H−jj
jj
jj 2 2
jj
2 ∂ H−
+ 2λ− v−
+ (λ− ) v− H− − v−
= 0,
∂t2
∂t
∂x2
72
x ∈ Ωj ,
(2.84)
(2.85)
donde g ij (x, t) es la fracción del total de neuronas que generan potenciales
de acción en x ∈ Ωj debido a las sinapsis que provienen de neuronas en Ωi .
En lo que sigue supondremos que
g ij (x, t) = µij (x)g j (x, t) ,
(2.86)
donde g j (x, t) está definida en (2.74) y 0 ≤ µij (x) ≤ 1 es tal que
n
X
µi,j (x) = 1, para toda x ∈ Ωj .
i=1
Las funciones µij (x); i, j = 1, . . . , n determinan el diagrama de conexión
funcional en la corteza, correspondiente a un estado fisiológico, caracterizado
por una función de activación de neuronas g j (x, t) para cada compartimento
Ωj ; j = 1, . . . , n. El hecho de que µij (x) pueda ser diferente de µkj para
i 6= k nos permite considerar una distribución no homógenea de activación
en x ∈ Ωj correspondiente a los compartimentos Ωk y Ωi que no podı́amos
considerar en el modelo estudiado en las secciones anteriores.
Sea V± el valor absoluto del potencial postsináptico producido por una
sinapsis (+ excitatoria, − inhibitoria). Introduciremos las nuevas variables
h
uij (x, t) =
V+ (n − 1) H+ij (x, t) +
f j (x,t)
Sj (x)
i
− V− H−jj (x, t)
(V+ + V− )
,
(2.87)
i, j = 1, . . . , n , que llamaremos variable de activación parcial en x ∈ Ωj , en
el instante t, provocada por la acción de Ωi sobre Ωj .
De (2.84), (2.85), (2.86), (2.87) se obtiene el siguiente sistema de ecuaciones para las variables de activación :
2
∂ 2 uij
∂uij
ij
ij
2
2 ∂ uij
v
v
)
u
−
v
+
2λ
+
(λ
=
+
+
ij
+
+
+
∂t2
∂t
∂x2
"
λij
∂g j
+ v+ V+ (n − 1)
ij
j
µij (x) λ+ v+ g (x, t) +
(x, t)
(V+ + V− )
∂t
#
+ Fij (x, t) − Iij (x, t), x ∈ Ωj ; i, j = 1, . . . , n ,
73
(2.88)
donde

³
´
j
2 f (x,t)
∂
V+ (n − 1)  ∂ Sj (x)
ij
Fij (x, t) =
v
+
2λ

+
+
(V+ + V− )
∂t2
2 2
+(λij
+ ) v+
j
f (x, t)
2
− v+
Sj (x)
∂2
³
f j (x,t)
Sj (x)
∂x2
³
f j (x,t)
Sj (x)
´
∂t
´

 ,
"
(2.89)
#
2 jj
V−
∂ 2 H−ij
∂H−ij
ij
ij
2 jj
2 ∂ H−
Iij (x, t) =
+
2λ
.
v
+
(λ
v
)
H
−
v
+ +
+ +
−
+
V+ + V− ∂t2
∂t
∂x2
(2.90)
Nótese que, en (2.90), H−jj (x, t) sólo depende de las condiciones iniciales
de inhibición, debido a que satisface la ecuación (2.85).
Teniendo en cuenta (2.80) podemos suponer que en (2.88) se tiene
g j (x, t) = gj (uj ) ,
donde
uj (x, t) =
n
X
(2.91)
uij (x, t)
i=1
es la llamada variable de activación total en x ∈ Ωj .
De
Ã
!
n
X
f j (x, t)
ij
0≤
≤ 1, para cada x ∈ Ωj
H+ (x, t) +
Sj (x)
i=1
y 0 ≤ H−jj (x, t) ≤ 1, para cada x ∈ Ωj ,
se tiene que
−1 ≤ uj (x, t) ≤ 1, para cada x ∈ Ωj .
Entonces de (2.91) y el comportamiento cualitativo que, a partir del
conocimiento fisiológico, debe tener la función de activación gj , podemos
suponer que gj (uj ) es una función creciente continuamente diferenciable de
la variable uj , definida en todo R y tal que
74
gj =


 0,


creciente ,
1,
si uj ≤ −1 ,
si − 1 ≤ uj ≤ 1 ,
si uj ≥ 1 .
(2.92)
Para tener más arbitrariedad en la selección de gj , podemos suponer que
es una función sigmoide de la variable uj que está muy cercana al valor 0 en
uj = −1 y muy cercana al valor 1 en uj = 1.
Haremos ahora la suposición de que las conexiones excitatorias entre compartimentos diferentes son pocas en comparación con las que hay dentro de
cada compartimento, es decir,
λjj
+
¿ 1, si i 6= j; i, j = 1, . . . , n .
λij
+
(2.93)
Entonces, si ponemos
λjj
λjj
+ v+
+
= βj ,
= ²,
jj
λ− v−
λij
+
i 6= j
(2.94)
y pasamos a las nuevas coordenadas
jj
λjj
+ v+ t = τ, λ+ x = y
(2.95)
en (2.88) para todo i, j = 1, . . . , n, obtenemos un sistema con parámetro
pequeño ² para las variables de activación parcial en coordenadas adimensionales
²2
2
∂ 2 ũij
∂ ũij
2 ∂ ũij
+
2²
+
ũ
−
²
=
ij
∂τ 2
∂τ
∂y 2
"
#
V+ (n − 1)
∂g̃ j
µ̃ij (y) g̃ j (y, τ ) + ²
(y, τ ) +
(V+ + V− )
∂τ
(2.96)
+
h
i
1
˜ij (y, τ ) ;
F̃
(y,
τ
)
−
I
ij
2
(λij
+ v+ )
75
1
y ∈ Ωj , i, j = 1, . . . , n,
λjj
+
i 6= j ,
donde
ũij (y, τ ) = uij (x, t), µ̃ij (y) = µij (x),
g̃ j (y, τ ) = g j (x, t), F̃ij (y, τ ) = Fij (x, τ ),
I˜ij (y, τ ) = Iij (x, t).
Además se tiene
∂ 2 ũjj
∂ ũjj
∂ 2 ũjj
+
2
+
ũ
−
=
jj
∂τ 2
∂τ
∂y 2
"
#
V+ (n − 1)
∂g̃
µ̃jj (y) g̃ j (y, τ ) +
(y, τ ) +
(V+ + V− )
∂τ
(2.97)
i
h
1
˜jj (y, τ ) , 1 y ∈ Ωj .
F̃
(y,
τ
)
−
I
jj
(λ+ v+ )2
λjj
+
De (2.91) tenemos que
g̃ j (y, τ ) = g̃j (ũj ) ,
donde
ũj (y, τ ) =
n
X
ũij (y, τ )
i=1
es la variable de activación total de la región Ωj en coordenadas adimensionales.
76
En lo que sigue consideraremos que las variables de activación se obtienen
en forma de promedios espaciales en cada compartimento Ωj , es decir
ũij (y, τ ) = ũij (τ ), ũij (y) = ũij ,
g̃ j (y, τ ) = g̃ j (τ ), F̃ij (y, τ ) = F̃ij (τ ) ,
I˜ij (y, τ ) = I˜ij (τ ) .
De esta forma, el sistema (2.96), (2.97) se reduce al sistema
²2
"
d2 ũij
dũij
+
2²
+ ũij =
dτ 2
dτ
#
i
h
dg̃ j
1
V+ (n − 1) j
˜
g̃ (τ ) + ²
F̃
µ̃ij
(τ ) + ij
ij (τ ) − Iij (τ ) ,
(V+ + V− )
dτ
(λ+ v+ )2
i, j = 1, . . . , n,
(2.98)
i 6= j ,
"
#
∂g̃ j
d2 ũjj
dũjj
V+ (n − 1) j
g̃
(τ
)
+
+
2
+
ũ
=
µ̃
(τ )
jj
ij
dτ 2
dτ
(V+ + V− )
∂τ
(2.99)
+
i
h
1
˜jj (τ ) .
F̃
(τ
)
−
I
jj
2
(λjj
+ v+ )
Por simplicidad supondremos que la entrada aferente Fij (τ ) a la región
Ωj proveniente de estructuras internas del cerebro es pequeña, es decir :
F̃ij (τ ) = (λjj
+ v+ )G̃ij (τ ), i 6= j ,
donde G̃ij (τ ) está acotada.
77
(2.100)
Por otra parte, de (2.90) se tiene que
"
#
2 jj
1
V−
dH̃−jj (τ )
2 d H̃− (τ )
˜ij (τ ) =
I
²
+
2²
+ H−jj (τ ) , (2.101)
2
2
(V
+
V
)
dτ
dτ
(λij
v
)
+
−
+ +
donde H̃−jj (τ ) = H−jj (τ ) es, según (2.85), solución de la ecuación
2 jj
2 d H̃−
βj
2
dH̃−jj
+ 2βj
+ H̃−jj = 0 ,
dτ
dτ
donde βj viene definido en (2.94). De (2.102) se tiene que
(2.102)
H̃−jj (τ ) = (aτ + b)e−τ /βj ,
donde a y b son parámetros que dependen de las condiciones iniciales de
inhibición.
Finalmente, la ecuación (2.98) se transforma en
²2
d2 ũij
dũij
+ 2²
+ ũij =
2
dτ
dτ
#
"
V+ (n − 1)
dg̃ j
µ̃ij g̃ j (τ ) + ²
(τ ) + ²2 G̃ij (τ )
(V+ + V− )
dτ
−
(2.103)
V−
H̃ jj (τ ) para cada i, j = 1, . . . , n, i 6= j .
(V+ + V− ) −
En la aproximación de orden cero, para valores de τ no cercanos a 0, la
solución ũ0ij (τ ) de (2.103) satisface la ecuación
ũ0ij =
V+ (n − 1)
V−
µ̃ij g̃j (ũ0j ) −
H̃−jj (τ ), i, j = 1, . . . , n, i 6= j ,
(V+ + V− )
(V+ + V− )
(2.104)
donde
ũ0j (τ ) =
n
X
i=1
78
ũ0ij (τ ),
de (2.104) se concluye que
ũ0j =
n
X
ũ0ij (τ ) =
i=1
n
X
V+ (n − 1)
(n − 1)V− jj
µ̃ij + ũ0jj −
g̃j (ũ0j )
H̃− (τ ) ,
(V+ + V− )
(V
+
V
)
+
−
i=1,i6=j
(2.105)
de donde se tiene :
ũ0j =
V+ (n − 1)
(n − 1)V− jj
(1 − µ̃jj )g̃(ũ0j ) − ũ0jj +
H̃ (τ ) .
(V+ + V− )
(V+ + V− ) −
(2.106)
Si suponemos que gj es una función sigmoide con un solo punto de inflexión, entonces la ecuación (2.106) admite una solución única u0j , para
cualquier ũ0jj , si se cumple la condición
V+ (n − 1)
0
(1 − µ̃jj )g̃j (u) < 1, para cada j = 1, . . . , n .
(V+ + V− )
(2.107)
En este caso, aplicando al sistema (2.99), (2.103) los resultados de la
teorı́a de perturbaciones singulares para sistemas con parámetro pequeño y
de la ecuación (2.106), tendremos que para cualquier compacto contenido en
{τ > 0} la aproximación de orden cero de la variable de activación total,
satisface la ecuación
dũ0j
d2 ũ0j
+2
+ ũ0j =
2
dτ
dτ
"
#
dũ0
V+ (n − 1)
0
µ̃jj g̃j (ũ0j ) + g̃j (ũ0j ) j +
(V+ + V− )
dτ
#
"
dg̃j 0
V+ (n − 1)
d2 g̃j 0
(1 − µ̃jj )
(ũj ) + 2
(ũ ) + g̃j (ũ0j )
2
(V+ + V− )
dτ
dτ j
"
dH̃−jj
(n − 1)V− d2 H̃−jj
+
2
+ H̃−jj
−
(V+ + V− ) dτ 2
dτ
79
#
+
1
2
(λjj
+ v+ )
h
i
F̃jj (τ ) − I˜jj (τ ) .
Teniendo en cuenta la expresión de H̃−jj y que
1
V−
I˜jj (τ ) =
jj
2
(V+ + V− )
(λ+ v+ )
=
Ã
d2 H̃−jj
dH̃ jj
+ 2 − + H̃−jj
dτ
dτ
!
V−
e−τ /βj (βj − 1) [2aβj + (βj − 1)(aτ + b)] .
(V+ + V− )βj2
Obtenemos finalmente la ecuación
"
#
d2 ũ0j
V+ (n − 1)
0
(1 − µ̃jj )gj (ũ0j )
+
1−
(V+ + V− )
dτ 2
"
#
dũ0j
V+ (n − 1)
0
2−
(2 − µ̃jj )gj (ũ0j )
+
(V+ + V− )
dτ
Ã
dũ0j
V+ (n − 1)
00
−
(1 − µ̃jj )gj (ũ0j )
(V+ + V− )
dτ
ũ0j
1
−
+
(2.108)
#
"
F̃jj (τ )
2
(λjj
+ v+ )
!2
V+ (n − 1)
g̃j (ũ0j ) =
−
(V+ + V− )
nV−
(βj − 1)e−τ /βj [2aβj + (βj − 1)(aτ + b)] ,
(V+ + V− )βj2
donde el primer sumando del miembro derecho está asociado a la entrada
aferente a la región Ωj proveniente de estructuras interiores del cerebro y el
segundo sumando corresponde al control inhibitorio en Ωj .
La ecuación (2.108) describe, en la aproximación de orden cero con respecto a las potencias de ², la evolución temporal de la variable de activación
total en cada compartimento Ωj .
80
Nótese que el coeficiente (1 − µ̃jj ) corresponde a la contribución de los
compartimentos Ωj con i 6= j a la activación de Ωj .
Nótese además que bajo las suposiciones efectuadas se obtienen ecuaciones desacopladas para la variable de activación total en cada Ωj , aunque
como señalamos en el párrafo anterior el efecto de los restantes compartimentos está encerrado en el coeficiente (1 − µ̃jj ).
Finalmente hemos obtenido el teorema siguiente.
Teorema 1 Supongamos que la corteza cerebral está dividida en n compartimentos mutuamente ajenos Ωj , j = 1, . . . , n y que se cumplen las siguientes
hipótesis:
1. La función de distribución de densidad de conexiones entre dos localizaciones de la corteza está descrita por (2.75), (2.76), (2.79).
2. La función que describe la fracción del total de neuronas en una localización determinada que generan potenciales de acción debido a las
sinapsis que provienen de neuronas ubicadas en otros compartimentos
viene dada por las expresiones (2.86), (2.91) y (2.92).
3. Las conexiones excitatorias entre compartimentos diferentes son pocos
en comparación con los que hay dentro de cada compartimento, según
se expresa en (2.93).
4. Las entradas sinápticas aferentes a la corteza son pequeñas en el sentido
de (2.100).
5. Entonces para la aproximación de orden cero de los promedios espaciales en cada Ωj de la variable de activación total
u0j (τ )
=
n
X
u0ij (τ )
(ver 2.87),
i=1
se tiene la ecuación (2.108), si se cumple la condición (2.107).
81
Capı́tulo 3
ANÁLISIS MATEMÁTICO DEL MODELO
DE ACTIVACIÓN DE LA CORTEZA CEREBRAL.
3.1 UN MODELO MATEMÁTICO PARA UNA REBANADA
DE CORTEZA CEREBRAL
Un problema importante en fisiologı́a es determinar bajo qué estı́mulos
externos se puede lograr que una pequeña rebanada de corteza cerebral oscile en frecuencias determinadas de antemano. Con el modelo que hemos
obtenido es posible estudiar este problema. Como estarı́amos considerando
únicamente una pequeña rebanada de corteza cerebral, lo que tendrı́amos que
hacer en nuestro modelo, es considerar que no se tiene variación espacial; en
este caso, nuestro modelo serı́a un sistema de dos ecuaciones diferenciales
ordinarias no lineales. Lo que nos proponemos estudiar, es las condiciones
que debe cumplir el estı́mulo externo para que nuestro sistema de ecuaciones
tenga soluciones periódicas.
Ası́ que, el sistema de ecuaciones diferenciales ordinarias a estudiar es
82
d2 u
dt2
h
+ 2λ+ v+ du
+ (λ+ v+ )2 u = V+ λ+ v+ λ+ v+ ḡ(u) + ḡ 0 (u) du
dt
dt
·
−V− e−λ− v− t (1 −
2
v+
d2 H̄−
)
2
dt2
v−
i
¸
−
+ 2(λ+ v+ − λ− v− ) dH̄
+ (λ+ v+ − λ− v− )2 H̄−
dt
+f+ (t) ,
(3.1)
d2 H̄−
= f− (t) ,
(3.2)
dt2
donde f+ (t), f− (t) son funciones que están relacionadas con los estı́mulos
externos a la pequeña rebanada de corteza.
La ecuación (3.2) es muy fácil de resolver, sus soluciones son de la forma
H̄− (t) = R(t) + at + b ,
(3.3)
donde a, b son constantes y
R(t) =
Z t µZ s
0
0
¶
f− (z)dz ds .
(3.4)
Sustituyendo (3.3) en (3.1), obtenemos la ecuación diferencial ordinaria
no lineal que necesitamos estudiar
d2 u
dt2
h
+ 2λ+ v+ du
+ (λ+ v+ )2 u = V+ λ+ v+ λ+ v+ ḡ(u) + ḡ 0 (u) du
dt
dt
·
−V− e
−λ− v− t
(1 −
i
2
v+
d2 R
2 ) dt2
v−
i
(3.5)
+2(λ+ v+ − λ− v− )( dR
+ a) + (λ+ v+ − λ− v− )2 (R(t) + at + b)
dt
+f+ (t) .
A continuación realizaremos un cambio de variable que nos permitirá
reducir la ecuación a una ecuación diferencial equivalente más sencilla de
estudiar.
83
Lema 3 Sea τ = λ+ v+ t, entonces la ecuación (3.5) es equivalente a la
ecuación diferencial
d2 u
dτ 2
h
+ 2 du
+ u = V+ ḡ(u) + ḡ 0 (u) du
dτ
dτ
−V− e
−
λ− v− τ
λ+ v+
·
(1 −
2
v+
d2 R
2 ) dτ 2
v−
i
+ 2(λ+ v+ − λ− v− )( (λ+1v+ )2 dR
+
dτ
)
+(λ+ v+ − λ− v− )2 ( (λR(τ
2 +
+ v+ )
aτ
(λ+ v+ )3
+
b
)
(λ+ v+ )2
i
a
)
(λ+ v+ )2
(3.6)
+f+ (τ ) .
Demostración.
Dividimos la ecuación (3.5) por (λ+ v+ )2 y usamos la regla de la cadena
para obtener las derivadas necesarias de las funciones que aparecen en la
ecuación con respecto a τ .
En el siguiente resultado construiremos, a partir de la ecuación (3.6),
una familia de ecuaciones diferenciales de primer orden dependiente de un
parámetro. Esta familia de ecuaciones diferenciales de primer orden es equivalente a la ecuación diferencial de segundo orden (3.6), en el sentido que se
enuncia en la siguiente proposición.
Proposición 1 Si u es una solución de la ecuación diferencial (3.6) con
condiciones iniciales u(0) = u0 y u0 (0) = u00 , entonces u es solución de la
ecuación diferencial de primer orden
du
+ u = V+ ḡ(u) + e−τ (F+ (τ ) + K) ,
dτ
84
(3.7)
donde
F+ (τ ) =
Ã
Rτ
−
−V− e
0
λ− v− ν
λ+ v+
·
(1 −
2
v+
d2 R
2 ) dν 2
v−
+ 2(λ+ v+ − λ− v− )( (λ+1v+ )2 dR
+
dν
+
(λ+ v+ − λ− v− )2 ( (λR(ν)
2 +
+ v+ )
a
)
(λ+ v+ )2
aν
(λ+ v+ )3
+
i
b
)
(λ+ v+ )2
´
+ f+ (ν) eν dν
(3.8)
y
K = u00 + u0 − V+ g(u0 ) .
Recı́procamente, si u es una solución de (3.7) con u(0) = u0 , entonces, u es
una solución de (3.6) con condiciones iniciales u(0) = u0 y u0 (0) = u00 .
Demostración.
Si multiplicamos la ecuación (3.6) por eτ , obtenemos
eτ
h
d2 u
dτ 2
i
h
i
+ 2 du
+ u = eτ V+ ḡ(u) + ḡ 0 (u) du
+
dτ
dτ
dF+
dτ
Como
"
"
#
d2 u
d τ du
du
e ( + u) = eτ
+2 +u
2
dτ
dτ
dτ
dτ
#
y
V+
d τ
du
[e ḡ(u)] = V+ [ḡ(u) + ḡ 0 (u) ]
dτ
dτ
(3.6) se transforma en
"
#
d
d τ du
e ( + u) =
[V+ eτ ḡ(u) + F+ (τ )] .
dτ
dτ
dτ
Entonces
85
.
(3.9)
eτ (
du
+ u) = V+ eτ ḡ(u) + F+ (τ ) + K ,
dτ
donde
K = u00 + u0 + V+ ḡ(u0 ) ,
entonces
du
+ u = V+ ḡ(u) + e−τ (F+ (τ ) + K) .
dτ
Es decir, u es solución de la ecuación (3.7). El recı́proco se demuestra de
manera similar.
Ahora consideraremos que la entrada externa inhibitoria a la corteza cerebral es cero, es decir, tomaremos f− (t) ≡ 0.
En este caso la ecuación (3.6) se reduce a
"
d2 u
du
du
0
+
2
+
u
=
V
ḡ(u)
+
ḡ
(u)
+
dτ 2
dτ
dτ
−V− e
−
λ− v− τ
λ+ v+
"
#
Ã
a
aτ
b
2(λ+ v+ − λ− v− )
+ (λ+ v+ − λ− v− )2
+
2
3
(λ+ v+ )
(λ+ v+ )
(λ+ v+ )2
+f+ (τ )
y
F+ (τ ) =
Rτ
Ã
−
−V− e
0
Ã
2
+(λ+ v+ − λ− v− )
λ− v− ν
λ+ v+
"
2(λ+ v+ − λ− v− )
b
aν
+
3
(λ+ v+ )
(λ+ v+ )2
!#
a
(λ+ v+ )2
!
+ f+ (ν) eν dν
En el siguiente resultado, probaremos que si escogemos el estı́mulo externo
de tal manera que e−τ (F+ (τ ) + K) sea una función periódica en τ de perı́odo
86
!#
T , se pueden proponer condiciones sobre ḡ para que la ecuación diferencial
(3.7) tenga una única solución periódica, de perı́odo T .
Observación
Para lograr que e−τ (F+ (τ ) + K) sea una función periódica de perı́odo T ,
basta tomar
f+ (τ ) = p(τ ) + q(τ )
donde p(τ ) sea una función cotinua periódica de perı́odo T y limτ →∞ q(τ ) = 0.
Además, en nuestro caso se puede demostrar que la función q(τ ) decrece
a cero en forma exponencial cuando τ tiende a ∞.
En lo que sigue B = B[0, ∞) = {u : [0, ∞) → R| u es continua y acotada
en [0, ∞)}, para u ∈ B, ||u||∞ = supt∈[0 ∞] |u(t)| sabemos que B es un espacio
de Banach con ||u||∞ .
Sea M = {u ∈ B : u(0) = u0 }, M es un subespacio cerrado de B.
Teorema 2 Si e−τ (F+ (τ )+K) es una función continua en [0, ∞) y periódica
de perı́odo T , ḡ una función tal que V+ supu∈R |ḡ 0 (u)| < 1, entonces la
ecuación diferencial
du
+ u = V+ ḡ(u) + e−τ (F+ (τ ) + K)
dτ
tiene una única solución acotada en [0 ∞)
(3.10)
Demostración.
Por el método de variación de parámetros, u es una solución de (3.10) si
y sólo si
−τ
u(τ ) = e u0 +
Z τ
0
(e−τ es [V+ ḡ(u) + e−s (F+ (s) + K)])ds .
(3.11)
Vamos a construir un operador T : M → M que tiene un único punto fijo
y este punto fijo será una solución de la ecuación (3.11) y por lo tanto de la
ecuación (3.10).
87
Definamos T : M → M por
−τ
T u(τ ) = e u0 +
Z τ
0
(e−τ es [V+ ḡ(u) + e−s (F+ (s) + K)])ds .
(3.12)
1. Vamos a demostrar que T u es un elemento del espacio M
|T u(τ )| ≤ |e−τ u0 | +
≤ |u0 | + V+ e−τ
Rτ
0
|e−τ +s |V+ ḡ(u(s))|ds +
Rτ
0
(|e−τ (F+ (s) + K)|)ds
R
Rτ s
−τ τ
0 |F+ (s) + K|ds
0 e ds + e
≤ |u0 | + V+ e−τ (eτ − e0 ) + ||F+ (τ ) + K||∞ e−τ τ ,
por lo tanto
|T u(τ )| ≤ |u0 | + V+ (1 − e−τ ) + ||F+ (τ ) + K||∞ eτ τ
≤ |u0 | + V+ + ||F+ (τ ) + K||∞ e−τ τ .
Como limτ →∞ τ e−τ = 0 se tiene que |T u(τ )| ≤ m para cada τ ∈ [0 ∞).
Ası́ si u ∈ M entonces, T u ∈ M
2. Ahora vamos a probar que T es una contracción. Sean u1 , u2 ∈ M ,
entonces
|T u1 (τ ) − T u2 (τ )| ≤
Rτ
0
|e−τ +s [V+ (ḡ(u1 (s))) − (ḡ(u2 (s)))]|ds
≤ V+ e−τ supu∈R |ḡ 0 (u)|||u1 − u2 ||∞
R∞ s
0 e ds
= supu∈R |ḡ 0 (u)|||u1 − u2 ||∞ V+ [1 − e−τ ]
≤ Sup|ḡ 0 (u)|u∈R ||u1 − u2 ||∞ V+ , por lo tanto
||T u1 − T u2 ||∞ ≤ k1 ||u − 1 − u2 ||∞ , con k1 = V+ supu∈R |ḡ 0 (u)| < 1.
El operador T es una contracción y por el teorema del punto fijo de
Banach, T tiene un único punto fijo u, el cual es una solución de la ecuación
diferencial (3.10)
En resumen hemos probado que si e−τ (F+ (τ )+K) es una función periódica,
de perı́odo T , entonces la ecuación (3.10) tiene una única solución acotada
88
en [0, ∞). Por el Teorema 4.11 de la página 116 del libro [39], se tiene que
la ecuación diferencial (3.9) tiene una solución periódica de perı́odo T , por
la proposición 1, u es una solución periódica de la ecuación (3.6).
Para garantizar que la solución periódica u cumpla con
| u(τ ) |≤ max{V+ , V− } ,
es necesario considerar, en lugar del conjunto M , el conjunto cerrado M 0 ,
donde M 0 es
M 0 = {u ∈ B | u(0) = u0 y exista τ0 ≥ 0 tal que | u(τ ) |≤ max{V+ , V− }
para cada τ ≥ τ0 } .
Se puede demostrar que el operador T definido en (3.12) tiene la propiedad
de que
T (M 0 ) ⊂ M 0
y ası́ la prueba del teorema 2 se puede repetir.
3.2 CONTROLABILIDAD DE LOS ESTADOS
FISIOLÓGICOS DE LA CORTEZA CEREBRAL.
En esta sección estudiaremos el problema del control de los estados fisiológicos de la corteza cerebral. Utilizaremos como elemento de control a
las inhibiciones. Se considera que el papel de las inhibiciones es controlar el
proceso de excitación. Esta hipótesis se refleja en el sistema (3.1), (3.2), ya
que la función de densidad de actividad sináptica inhibitoria (H¯− ) aparece
en la ecuación para la función de activación (u), pero no al revés.
En lo que sigue consideraremos que
1. No existen entradas inhibitorias a la corteza cerebral, esto se traduce
en que f− (t) ≡ 0 en la ecuación diferencial (3.1).
2. Los tiempos caracterı́sticos de las inhibiciones son mucho menores que
los de las excitaciones. Por esto consideraremos que
λ + v+
¿ 1.
λ − v−
89
(3.13)
Ası́, la ecuación diferencial que estudiaremos es
d2 u
dt2
h
+ (λ+ v+ )2 u = V+ λ+ v+ λ+ v+ ḡ(u) + ḡ 0 (u) du
+ 2λ+ v+ du
dt
dt
i
−V− e−λ− v− t [2(λ+ v+ − λ− v− )α + (λ+ v+ − λ− v− )2 (αt + β)]
(3.14)
+f+ (t) .
ya que H̄− debe satisfacer la ecuación (3.2) y por lo tanto, debe ser de la
forma
H̄− (t) = αt + β .
Usando el cambio de variable τ = λ+ v+ t y el lema 3, la ecuación diferencial (3.14) es equivalente a la ecuación
d2 u
dτ 2
h
+ 2 du
+ u = V+ ḡ(u) + ḡ 0 (u) du
dτ
dτ
λ
−V− e
v
− λ− v− τ
+ +
i
h
2(λ+ v+ − λ− v− ) (λ+av+ )2 + (λ+ v+ − λ− v− )2 ( (λ+av+ )3 τ +
+f+ (τ ) .
(3.15)
Ası́, el problema de control que queremos estudiar es el siguiente: dados
u0 , u00 , T, uT y u0T , lo que queremos demostrar es que existen a, b números
reales y u solución de la ecuación (3.15) tal que
u(0) = u0 ,
u(T ) = uT
y
u0 (0) = u00 ,
u0 (T ) = u0T .
Es decir queremos ver si el sistema dinámico determinado por la ecuación
(3.15) es controlable.
Para estudiar este problema, nos auxiliaremos de la proposición 1, la cual
nos garantiza que u es una solución de (3.15) que satisface las condiciones
iniciales u(0) = u0 , u0 (0) = u00 si, y sólo si, u satisface la ecuación diferencial
de primer orden
90
i
b
)
(λ+ v+ )2 )
du
+ u − V+ ḡ(u) = (f0 (τ ) + ca,b (τ ) + K)e−τ ,
dτ
(3.16)
K = u0 (0) + u(0) − V+ ḡ(u(0)) ,
(3.17)
donde
f0 (t) es una función relacionada con las entradas excitatorias a la corteza
cerebral y ca,b (t) es la función
(" ³
ca,b (t) = V−
e
" ³
+
e
λ
λ
v
´ "Ã
1− λ− v− t
+ +
´
v
1− λ− v− t
+ +
#Ã
−1
!
#
#
λ − v−
t
1
1
−1
−
+
a
λ + v+
λ+ v+ λ+ v+
λ + v+
! )
λ− v−
−1 b .
λ+ v+
(3.18)
Las constantes a y b están relacionadas con las condiciones iniciales de la
inhibición.
El problema de control que estudiaremos es: dado u0 , u00 , T > 0, uT , u0T
¿existirán a, b números reales y u solución de la ecuación (3.16) tal que
u(0) = u0 ,
u(T ) = uT ,
?
0
u (0) =
u00
,
0
u (T ) =
u0T
.
Este problema es equivalente a probar que la función U : R2 −→ R2
definida por
à !
Ã
!
a
u(T, a, b)
U
= 0
b
u (T, a, b)
es sobreyectiva, donde u(t, a, b) es la solución de la ecuación diferencial (3.16)
que cumple la condición inicial u(0) = u0 .
Nota. Por la importancia que tienen las constantes a, b, es que ahora
denotaremos a las soluciones de (3.16) por u(t, a, b). La dependencia de las
condiciones iniciales se omitirá, ya que éstas condiciones estarán fijas.
91
Sabemos que U es una función diferenciable en R2 .
Para demostrar que U es suprayectiva, probaremos que U es un mapeo
abierto y también que U (R2 ) es un conjunto cerrado, y como R2 es conexo,
se debe tener que U (R2 ) = R2 .
Para simplificar la notación, introduciremos las siguientes funciones
v(t) = u(t, a + ², b) − u(t, a, b) ,
³
f1 (t) = e
λ
v
´ "Ã
1− λ− v− t
+ +
à ³
f2 (t) =
Como
λ+ v+
λ− v−
e
!
(3.19)
#
λ − v−
t
1
1
−1
−
+
, (3.20)
λ + v+
λ + v+ λ + v+
λ+ v+
λ
´
v
1− λ− v− t
+ +
!Ã
−1
!
λ − v−
−1 .
λ + v+
(3.21)
¿ 1 se tiene que f1 (t) > 0 y f2 (t) < 0 para cada t ≥ 0 .
Lema 4 Para cada 0 < τ ≤ T se cumple
f1 (τ )f2 (T ) − f2 (τ )f1 (T ) < 0 .
(3.22)
Demostración.
Sea
f3 (τ ) =
f1 (τ )
.
−f2 (τ )
(3.23)
Probaremos primero que f3 es una función decreciente en (0, T ], demostrando que f30 ≤ 0 en (0, T ].
Multiplicando y dividiendo (3.23) por
λ
e
v
−(1− λ− v− )τ
+ +
y sustituyendo (3.20) y (3.21) en (3.23), se tiene
92


1 
f3 (τ ) =
λ + v+
τ
λ
e
+
v
( λ− v− −1)τ
+ +
1
λ− v −
λ+ v +
−1
−1
.
(3.24)
Derivando (3.24), se tiene
³
f30 (τ )
e
λ− v−
−1
λ+ v+
´
³
"
τ
λ
v
´
1− λ− v− τ
=Ã ³
´
!2 1 − e
λ− v−
−1
τ
−1
e λ+ v+
+ +
Ã
! #
λ− v−
−
−1 τ .
λ+ v+
(3.25)
El primer factor de (3.25) es mayor que cero, ası́ que, para probar que
< 0 en (0, T ], necesitamos probar que el segundo factor de (3.25) es menor
que cero; para esto, consideremos la función
f30
³
f4 (τ ) = 1 − e
λ
´
v
1− λ− v− τ
+ +
Ã
!
λ − v−
−
−1 τ.
λ + v+
Como f4 (0) = 0, entonces, si probamos que f4 es decreciente en (0, T ], se
tendrı́a que f4 < 0 en (0, T ]. Pero
Ã
f40 (τ ) =
! ³
λ− v−
−1 e
λ+ v+
λ
v
´
1− λ− v− τ
+ +
Ã
!2 Ã ³
λ − v−
−
−1
λ + v+
e
λ
v
´
1− λ− v− τ
+ +
!
− 1 ≤ 0,
por lo tanto, f4 es una función decreciente en (0, T ]. Ası́, f4 (τ ) ≤ 0 para
τ ∈ (0, T ]. Esto último implica que f30 ≤ 0 en (0, T ], que es lo que querı́amos
probar.
Ahora, por ser f4 una función decreciente, para cada τ ∈ (0, T ],se tiene
f1 (t)
f1 (T )
≥
−f2 (t)
−f2 (T )
y como f2 (τ ) < 0 para cada τ , se tiene
f1 (τ )f2 (T ) − f2 (τ )f1 (T ) < 0 .
En lo que sigue, vamos a necesitar las siguientes derivadas
∂u0
∂u0
∂u
(T, a, b), ∂u
(T,
a,
b),
(T,
a,
b)
y
(T, a, b) .
∂b
∂a
∂a
∂b
Ası́ que vamos a calcular estas derivadas.
93
1.
∂u
(T, a, b) .
∂a
Como u(t, a, b), u(t, a + ², b) son soluciones de (3.16) se cumple que
du
(τ, a
dτ
+ ², b) + u(τ, a + ², b) − V+ ḡ(u(τ, a + ², b)) =
(3.26)
−τ
(f0 (τ ) + ca+²,b (τ ) + K)e
du
(τ, a, b)
dτ
,
+ u(τ, a, b) − V+ ḡ(u(τ, a, b)) =
(3.27)
−τ
(f0 (τ ) + ca,b (τ ) + K)e
,
entonces, restando (3.27) a (3.26) y usando (3.19) se tiene
dv
+ v − V+ (ḡ(u(τ, a + ², b)) − ḡ(u(τ, a, b))) = (ca+²,b (τ ) − ca,b (τ ))e−τ .
dτ
Ahora, tomando la parte lineal de ḡ(u(τ, a, b)), se tiene
dv
+ v − V+ (ḡ 0 (u(τ, a, b)))v = (ca+²,b (τ ) − ca,b (τ ))e−τ ,
dτ
pero como
ca+²,b (τ ) − ca,b (τ ) = V− ²f1 (τ )
y
v(T ) = V− ²
Z T
0
e−t e−
RT
t
(1−V+ ḡ 0 (u(τ,a,b)))dτ
f1 (t) dt ,
(3.28)
entonces, usando (3.19) y dividiendo (3.28) por ², se tiene
u(T, a + ², b) − u(T, a, b)
v(T )
=
²
²
= V−
94
Z T
0
e−t e−
RT
t
(1−V+ ḡ 0 (u(τ,a,b)))dτ
f1 (t) dt .
Tomando el lı́mite cuando ² −→ 0, se obtiene
Z T
RT
∂u
−t − t (1−V+ ḡ 0 (u(τ,a,b)))dτ
(T, a, b) = V−
e e
f1 (t) dt . (3.29)
∂a
0
De manera similar, se puede calcular
2.
∂u
(T, a, b)
∂b
.
∂u
(T, a, b).
∂b
Z T
RT
0
∂u
(T, a, b) = V−
e−t e− t (1−V+ ḡ (u(τ,a,b)))dτ f2 (t) dt.
∂b
0
3.
(3.30)
∂u0
(T, a, b).
∂a
Como u(t, a, b), u(t, a + ², b) son soluciones de (3.16) se cumple
u0 (T, a, b) = V+ ḡ(u(T, a, b)) − u(T, a, b) + (f (T ) + c(T ) + K)e−T (3.31)
y
u0 (T, a+², b) = V+ ḡ(u(T, a+², b))−u(T, a+², b)+(f (T )+c(T )+K)e−T ,
(3.32)
Restando (3.31) a (3.32), se obtiene
u0 (T, a + ², b) − u0 (T, a, b) = V+ [ḡ(u(T, a + ², b)) − ḡ(u(T, a, b))]
− [u(T, a + ², b) − u(T, a, b)] + [ca+²,b (T ) − ca,b (T )]
= [V+ ḡ 0 (u(T, a, b)) − 1] [u(T, a + ², b) − u(T, a, b)] + e−T ²V− f1 (T ) .
Entonces
u0 (T, a + ², b) − u0 (T, a, b)
=
²
95
"
#
u(T, a + ², b) − u(T, a, b)
[V+ ḡ (u(T, a, b)) − 1]
+ e−T V− f1 (T ) ,
²
0
Tomando el lı́mite cuando ² −→ 0, obtenemos
∂u0
∂u
(T, a, b) = [V+ ḡ 0 (u(T, a, b)) − 1] (T, a, b) + V− e−T f1 (T ) . (3.33)
∂a
∂a
De la misma manera, se obtiene
4.
∂u0
(T, a, b).
∂b
∂u0
∂u
(T, a, b) = [V+ ḡ 0 (u(T, a, b)) − 1] (T, a, b) + V− e−T f2 (T ) . (3.34)
∂b
∂b
Proposición 2 Sea U : R2 −→ R2 definida por
à !
Ã
!
a
u(T, a, b)
= 0
,
U
b
u (T, a, b)
donde u(t, a, b) es la solución de la ecuación diferencial (3.16) que cumple la
condición inicial u(0) = u0 . Entonces U es un mapeo abierto.
Demostración.
Basta probar que localmente U es un mapeo abierto. Para probar que U
es localmente un mapeo abierto, demostraremos que
 ∂u
(T, a, b)
 ∂a
DU (a, b) = det 
∂u0
(T, a, b)
∂a
para cada a, b elementos de R.
96
∂u
(T, a, b)
∂b
∂u0
(T, a, b)
∂b


 6= 0
(3.35)
DU (a, b) =
∂u
∂u0
∂u0
∂u
(T, a, b)
(T, a, b) −
(T, a, b) (T, a, b) .
∂a
∂b
∂a
∂b
(3.36)
Sustituyendo (3.33) y (3.34) en (3.36) se obtiene
"
#
∂u
∂u
DU (a, b) =
(T, a, b)
(T, a, b) [V+ ḡ 0 (u(T, a, b)) − 1] + V− e−T f2 (T )
∂a
∂b
"
#
∂u
∂u
∂u
− (T, a, b)
(T, a, b) [V+ ḡ 0 (u(T, a, b)) − 1] − V− (T, a, b)f1 (T ) ,
∂b
∂a
∂b
entonces
(
−T
DU (a, b) = V− e
)
∂u
∂u
(T, a, b)f2 (T ) −
(T, a, b)f1 (T )
∂a
∂b
.
(3.37)
Sustituyendo (3.29) y (3.30) en (3.37) se obtiene
"Z
DU (a, b) = V−2 e−T
T
0
−s −
e e
RT
s
0
(1−V+ ḡ (u))dθ
#
{f1 (s)f2 (T ) − f2 (s)f1 (T )} ds .
(3.38)
Por el lema 4, se tiene
f1 (t)f2 (T ) − f2 (t)f1 (T ) < 0 para 0 < t ≤ T
y por lo tanto el integrando en (3.38) es estrictamente menor que 0 en el
intervalo [0, T ], y por lo tanto
DU (a, b) < 0 para cada a, b ∈ R .
Ası́, hemos probado que U es localmente un mapeo abierto y por lo tanto
que U es un mapeo abierto en R2 .
Ahora demostraremos que U (R2 ) es un conjunto cerrado en R2 .
Primeramente probaremos el siguiente lema.
97
Lema 5 Para cada T > 0 sea BT = {(u(T, a, b), u0 (T, a, b)) : (a, b) ∈ A}. Si
BT es un conjunto acotado, entonces A es un conjunto acotado.
Demostración.
Primero demostraremos que los conjuntos
(
−T
α(T, a, b) = e
)
Z T
ca,b (s)ds : (a, b) ∈ A
0
(3.39)
y
{β(T, a, b) = ca,b (T ) : (a, b) ∈ A}
(3.40)
son acotados.
Como u es solución de la ecuación (3.16), es decir satisface la ecuación
du
+ u − V+ ḡ(u) = (f0 (τ ) + ca,b (τ ) + K)e−τ ,
dτ
entonces, por el método de variación de parámetros, u debe satisfacer
−τ
−τ
u(τ ) = e u(0)+e V+
Z τ
0
s
−τ
ḡ(u(s))e ds+e
Z τ
0
(f (s)+ca,b (s)+K)ds . (3.41)
Evaluando (3.41) en T se tiene
−T
u(T, a, b) = u(0)e
+ V+
Z T
0
ḡ(u(s))e
s−T
−T
ds + e
Z T
0
(f (s) + ca,b (s) + K)ds .
(3.42)
Despejando α(T, a, b) de (3.42) se tiene
−T
α(T, a, b) = u(T, a, b)−u(0)e
−V+
Z T
0
s−T
ḡ(u(s))e
ds−e
−T
Z T
0
(f0 (s)+K)ds ,
(3.43)
entonces
| α(T, a, b) |≤| u(T, a, b) | + | u(0)e−T | +
| V+
RT
0
s−T
ḡ(u(s))e
−T
ds | + | e
98
RT
0
(3.44)
(f0 (s) + K)ds | .
Ahora, por hipótesis existe N1 (T ) tal que
| u(T, a, b) |≤ N1 (T ) para cada (a, b) ∈ A .
(3.45)
Además, como ḡ es acotada por 1 en R y e−τ ≤ 1 para τ ≤ 0, se tiene
¯
¯
Z T
¯
¯
¯
¯
s−T
¯ ≤ V+ T
¯ V+
ḡ(u(s))e
ds
¯
¯
0
(3.46)
y como f0 (τ ) es continua en el intervalo [0, T ], se tiene que, existe M (T ) tal
que
¯
¯Z
¯ T
¯
¯
¯
¯
(f
0 (s) + K)ds ¯ ≤ M (T ) .
¯ 0
¯
(3.47)
Usando (3.45), (3.46), (3.47) en (3.44) se prueba que
| α(T, a, b) |≤ N1 (T ) + u0 + V+ (T ) + M (T ) para cada (a, b) ∈ A
(3.48)
que es lo que querı́amos probar.
Por (3.41) se tiene
u0 (T, a, b) + u(T, a, b) − V+ ḡ(u(T, a, b)) = (f (T ) + ca,b (T ) + K)e−T . (3.49)
Despejando β(T, a, b) de (3.49) se tiene que
β(T, a, b) = u0 (T, a, b) + u(T, a, b) − V+ ḡ(u(T, a, b)) − (f (T ) + K)e−T , (3.50)
entonces,
| β(T, a, b) |≤| u0 (T, a, b) | + | u(T, a, b) | +
(3.51)
| V+ ḡ(u(T, a, b)) | + | (f (T ) + K)e−T | .
Por la hipótesis existe N2 (T ) tal que
| u0 (T, a, b) |≤ N2 (T ), para cada (a, b) ∈ A .
(3.52)
Usando (3.45) y (3.52) se tiene que
| β(T, a, b) |≤ N2 (T ) + N1 (T ) + V+ + M1 (T ), para cada (a, b) ∈ A (3.53)
99
y por lo tanto se obtiene que
{β(T, a, b) = ca,b (T ) : (a, b) ∈ A}
es un conjunto acotado en R.
Ahora, demostraremos que a y b deben satisfacer un par de ecuaciones.
Z T
α(T, a, b) =
³
−V− e
λ
v
´
1− λ− v− T
+ +
ca,b (s)ds =
0
³
(
−
1
a + V− (−1)e
λ + v+
λ
v
´
1− λ− v− T
+ +
Ã
!
)
λ − v−
+1−
−1 T b
λ + v+
= A(T )a + B(T )b ,
donde
³
A(T ) = −V− e
³
(
λ
λ
v
+ +
v
´
− 1− λ− v− T
B(T ) = V− (−1)e
´
1− λ− v− T
+ +
1
,
λ + v+
Ã
!
λ − v−
−1 T
+1−
λ + v+
)
y
β(T, a, b) = V− f1 (T )a + f2 (T )b .
Por lo tanto, tenemos que se cumple
A(T )a + B(T )b = α(T, a, b) ,
(3.54)
V− f1 (T )a + V− f2 (T )b = β(T, a, b) .
Si probamos que



A(T )
B(T )
V− f1 (T ) V− f2 (T )
100



es invertible, se tendrı́a que



a
b




=
A(T )
−1
B(T )
Ã


V− f1 (T ) V− f2 (T )
α(T, a, b)
β(T, a, b)
!
para cada (a, b) ∈ A.
Pero como hemos demostrado que α(T, a, b), β(T, a, b) son conjuntos acotados para cada T > 0, se concluye que A es un conjunto acotado.
Probaremos que
A(T )f2 (T ) − B(T )f1 (T )
6= 0 .
V−2
³
A(T )f2 (T ) − B(T )f1 (T )
=e
V−2
³
Ã
λ
´
v
´
v
+ +
Ã
1
f2 (T )−
λ + v+
! !
λ− v−
+1−
− 1 t f1 (T ) .
λ+ v+
1− λ− v− T
−e
λ
1− λ− v− T
(3.55)
+ +
Como f2 (T ) < 0,
³
λ
´
v
1− λ− v− T
e
+ +
1
f2 (T ) < 0 .
λ+ v+
Sabemos también que f1 (T ) > 0, ası́ si probamos que
à ³
e
λ
´
v
1− λ− v− T
+ +
Ã
!
λ − v−
−1+
−1 T
λ + v+
!
< 0,
(3.56)
demostrarı́amos la desigualdad
A(T )V− f2 (T ) − B(T )V− f1 (T )
< 0.
V−2
Para probar (3.56), consideremos la función
à ³
f5 (t) = e
λ
v
´
1− λ− v− t
+ +
Ã
!
λ − v−
−1+
−1 T
λ + v+
101
!
.
Se tiene que
f5 (0) = 0
y su derivada es
! Ã ³
Ã
f50 (t)
Como
λ− v−
λ+ v+
λ − v−
t e
= 1−
λ + v+
λ
´
v
1− λ− v− t
+ +
!
−1
.
À 1, se concluye que
f50 (t) < 0 para cada t ≥ 0 .
Por lo tanto f5 es un función decreciente en [0, ∞) y como f5 (0) = 0, se
tiene
f5 (T ) ≤ f5 (0) = 0 .
Con esto último probamos (3.56).
Por lo tanto hemos demostrado que



A(T )
B(T )
V− f1 (T )



V− f2 (T )
es invertible. Que es lo que querı́amos demostrar.
En la siguiente proposición, probaremos que U (R2 ) es un conjunto cerrado
en R2 .
Proposición 3 Si U : R2 −→ R2 definida por


a
U




= 

u(T, a, b)


 ,
u0 (T, a, b)
b
donde u(t, a, b) es una solución de (3.16) con u(0) = u0 . Entonces U (R2 ) es
un conjunto cerrado en R2 .
102
Demostración.
Sea (uT , u0T ) ∈ U (R2 ), existe una sucesión {(an , bn ) : n ∈ N } tal que



lim 
U 
n→∞
an




 = lim 
n→∞
bn
u(T, an , bn )




=
u0 (T, an , bn )
uT


 .
u0T
³ ´
a
Por demostrar que existen b en R2 y u solución de (3.16) tal que

U

a




=
b
u(T, a, b)




=
u0 (T, a, b)
uT


 .
u0T
Como {u(T, an , bn ), u0 (T, an , bn ) : n ∈ N } es un conjunto acotado en R2 ,
por el lema 5 se tiene que
{(an , bn ) | n ∈ N }
es un conjunto acotado en R2 , por el teorema de Bolzano-Weierstrass, existe
una subsucesión convergente {(anj , bnj ) | j ∈ N } de {(an , bn ) | n ∈ N }. Sea
(a, b) el lı́mite de esta subsucesión.
Por la continuidad de las soluciones u(t, anj , bnj ) se tiene que
lim u(T, anj , bnj ) = u(T, a, b) ,
j→∞
ası́ u(T, a, b) = uT .
Además, como
u0 (T, anj , bnj ) = −u(T, anj , bnj ) + V+ ḡ(u(T, anj , bnj ))
+(f (T ) + canj ,bnj (T ) + K)e−T ,
lim (−u(T, anj , bnj )) = −u(T, a, b) ,
j→∞
103
lim (V+ ḡ(u(T, anj , bnj ))) = V+ ḡ(u(T, a, b)) ,
j→∞
y
lim ((f (T ) + canj ,bnj (T ) + K)e−T ) = (f (T ) + ca,b (T ) + K)e−T
j→∞
se tiene que
lim u0 (T, anj , bnj ) = u0 (T, a, b) .
j→∞
por lo tanto hemos demostrado que U (R2 ) es un conjunto cerrado en R2 .
En resumen, hemos probado que el mapeo
U : R2 −→ R2
definido por


a
U




=
u(T, a, b)


 ,
u0 (T, a, b)
b
donde u(T, a, b) es una solución de (3.16) con condición inicial u(0) = u0 es
un mapeo sobreyectivo.
Teorema 3 El sistema dinámico
d2 u
dτ 2
h
+ 2 du
+ u = V+ ḡ(u) + ḡ 0 (u) du
dτ
dτ
λ
−V− e
v
− λ− v− τ
+ +
i
h
2(λ+ v+ − λ− v− ) (λ+av+ )2 + (λ+ v+ − λ− v− )2 ( (λ+av+ )3 τ +
+f+ (τ )
(3.57)
es controlable en R2 , es decir, dados u0 , u00 y T, uT y u0T , existen a, b números
reales y u solución de la ecuación (3.57) tal que
u(0) = u0 ,
u(T ) = uT
y
u0 (T ) = u0T .
u0 (0) = u00 ,
104
i
b
)
(λ+ v+ )2 )
Demostración.
Sean u0 , u00 y T, uT y u0T dados. Por la proposición 2, U : R2 → R2
definido por


a
U




=
b
u(T, a, b)


 ,
u0 (T, a, b)
donde u(t, a, b) es una solución de
du
+ u − V+ ḡ(u) = (f0 (τ ) + ca,b (τ ) + K)e−τ ,
dτ
con
K = u0 (0) + u(0) − V+ ḡ(u(0)),
es un mapeo abierto en R2 . Por la proposición 3, U es un mapeo cerrado en
R2 , por lo tanto como R2 es conexo, se tiene que U es suprayectiva, es decir,
existen a, b ∈ R y u solución de
du
+ u − V+ ḡ(u) = (f0 (τ ) + ca,b (τ ) + K)e−τ ,
dτ
con
u(0) = u0 ,
u(T ) = uT ,
u0 (0) = u00 ,
u0 (T ) = u0T
y
K = u0 (0) + u(0) − V+ ḡ(u(0)) .
Pero por la proposición 1, u también es solución de (3.57), y ası́ hemos
demostrado el teorema.
105
3.3 EXISTENCIA DE ONDAS VIAJERAS EN LA CORTEZA
CEREBRAL
En esta sección estudiaremos la existencia de estados fisiológicos que se
corresponden con la existencia de ondas viajeras (ondas que se propagan por
el medio cortical sin perder su forma) en la ecuación (2.71).
En el caso de modelos matemáticos de neuronas se ha demostrado la
existencia de ondas viajeras [30], [48], [101]. También se ha estudiado este
problema en otro tipo de procesos de reacción difusión [102].
Buscaremos bajo qué condiciones la ecuación diferencial parcial
"
#
2
∂2u
∂u
∂u
2
2 ∂ u
0
+
2λ
v
+
(λ
v
)
u
−
v
=
V
λ
v
λ
v
ḡ(u)
+
ḡ
(u)
+
+
+
+
+
+
+
+
+
+
∂t2
∂t
∂x2
∂t
(3.58)
tenga soluciones de la forma: u(x, t) = φ(x − ct) para algunos valores de la
constante c.
Como en la sección anterior, haremos algunos cambios de variable, que no
alteran el resultado buscado, pero que nos permitirán construir una ecuación
diferencial parcial equivalente a (3.58) más sencilla de estudiar.
Si hacemos el cambio de variable τ = λ+ v+ t y y = λ+ x, la ecuación (3.58)
es equivalente a la ecuación
"
#
∂2u
∂u
∂ 2u
∂u
+2
+ u − 2 = V+ ḡ(u) + ḡ 0 (u)
.
2
∂τ
∂τ
∂y
∂t
(3.59)
Para esta ecuación buscaremos soluciones de la forma
u(y, τ ) = φ(y − cτ ) = φ(s) con s = y − cτ y tal que φ(−∞) = s0 = φ(∞)
donde s0 es constante.
Antes de estudiar este tema, enunciaremos algunos resultados importantes que usaremos más adelante en esta sección. Estos resultados no serán
demostrados. Para ver la demostración de los Teoremas 4 y 5, recurrir a [23].
Teorema 4 Sea A una matriz n × n y f : R → Rn una función continua
acotada Consideremos la ecuación diferencial
106
du
= Au + f (t) .
(3.60)
dt
La condición necesaria y suficiente para que la ecuación (3.60) tenga una
única solución acotada en R, es que σ(A) no intersecte al eje imaginario.
En este caso la solución acotada se expresa en la forma.
u(t) =
Z ∞
−∞
G(t − τ )f (τ )dτ ,
(3.61)
donde G(t) es la función de Green principal de A, la cual se define como
(
G(t) =
eAt P− ,
−eAt P+ ,
t>0
t < 0,
(3.62)
donde P− , P+ son, respectivamente los proyectores sobre los subespacios propios con parte real positiva, negativa.
Teorema 5 Consideremos la ecuación
du
= Au + f (t, u) .
dt
Supongamos que se cumplen las condiciones siguientes
(3.63)
1. A es una matriz n × n, tal que
σ(A) ∩ ({0} × R) = ∅ .
(3.64)
2. Sea ρ > 0. Existe M > 0 tal que
| f (t, u) |≤ M
para t ∈ R y | u |≤ ρ .
(3.65)
3. Existe q > 0 tal que
| f (t, u1 ) − f (t, u2 |≤ q | u1 − u2 | para | u1 |≤ ρ y | u2 |≤ ρ . (3.66)
107
4. Existe N > 0 y r > 0 tal que
| G(t) |≤ N e−rt para cada t ∈ R .
(3.67)
5. M y q cumplen con
M≤
ρr
r
y q<
.
2N
2N
(3.68)
Entonces, la ecuación (3.63) tiene una única solución acotada en R, tal
que
| u(t) |≤ ρ para cada t ∈ R .
(3.69)
Ahora, demostraremos un lema, que usaremos posteriormente.
Lema 6 Sea f : R2 → R2 continua. Supongamos que la ecuación diferencial
dx
= f (x),
x ∈ R2
(3.70)
dt
tiene un único punto estacionario x0 , el cual es un punto silla. Entonces
la condición necesaria y suficiente para que la ecuación (3.70) tenga una
solución homoclı́nica es que tenga una solución no trivial acotada en R.
Demostración.
La suficiencia es evidente.
Sea φ una solución no trivial acotada de (3.70). Como x0 es el único
punto estacionario de la ecuación (3.70), y como además es un punto silla, el
teorema de Poincaré-Bendixon, nos asegura que
1. limt→±∞ φ(t) = x0
o
2. ω(φ) = {x0 } = α(φ), donde ω(φ) y α(φ) son, respectivamente los
conjuntos ω-lı́mite y α-lı́mite de la órbita generada por la solución φ .
108
En ambos casos, φ serı́a la solución homoclı́nica no trivial buscada.
En lo que sigue, retomaremos el problema de la existencia de las ondas
viajeras para la ecuación (3.59)
Ã
!
∂u
∂2u
∂u
∂ 2u
+
2
+
u
−
= V+ ḡ 0 (u)
+ ḡ(u) .
2
2
∂τ
∂τ
∂y
∂τ
Es decir, buscaremos soluciones de la forma
u(y, τ ) = φ(y − cτ ) ,
(3.71)
donde φ es una función de una variable que cumple que el
lim φ(s) = φ0 ,
s→±∞
(3.72)
donde φ0 es una constante.
Ahora veremos que la existencia de soluciones de la ecuación (3.59) que
son de la forma (3.71) y que cumplen con la condición (3.72), es equivalente
a la existencia de soluciones de una ecuación diferencial ordinaria de segundo
orden, que cumplen con la condición (3.72).
Si u es una solución de la ecuación (3.59) de la forma (3.71), que cumple
con la condición (3.70), entonces usando la regla de la cadena, obtenemos
∂u = −c dφ ,
∂τ
ds
2
∂ 2 u = c2 d φ ,
∂τ 2
ds2
∂u = dφ ,
∂y
ds
2
∂ 2u = d φ .
2
∂y
ds2
Sustituyendo estas derivadas en la ecuación (3.59) y reordenando la ecuación,
obtenemos que φ debe ser solución de la ecuación
d2 φ
dφ
(c − 1) 2 + c(V+ ḡ 0 (φ) − 2)
+ (φ − V+ ḡ(φ)) = 0 .
(3.73)
ds
ds
Recı́procamente, si φ es solución de (3.73) que cumple con la condición
(3.72) entonces la función
2
u(y, τ ) = φ(y − cτ )
es una solución de la ecuación (3.59).
109
Observación.
En lo que sigue supondremos que 0 < c < 1, ya que esto corresponde a
una velocidad real de propagación, cV+ , la cual no puede ser superior a la
velocidad de propagación de los potenciales de acción.
También tenemos que la ecuación (3.73) es equivalente al sistema de ecuaciones de primer orden.
à !

d φ
=

dt ψ
0
1
1
1−c2
−2c
1−c2



à !
Ã
!
0
φ
V+
,
+
ψ
1 − c2 cḡ 0 (φ)ψ − ḡ(φ)
(3.74)
donde ψ(s) = dφ
ds
La existencia de ondas viajeras para la ecuación (3.59) es equivalente a
la existencia de soluciones homoclı́nicas para la ecuación (3.74) para algún
valor de c ∈ (0, 1).
De (3.73) es evidente que para c = 1 no puede haber soluciones homoclı́nicas para esta ecuación.
En lo que sigue supondremos que
ḡ 0 (φ) <
1
para cada φ ∈ R .
V+
(3.75)
Con esta condición veremos que el sistema (3.74) tiene un único punto
estacionario para cualquier valor de la constante c que cumple 0 < c < 1.
Además, este punto estacionario es independiente de c y es un punto
silla del sistema. Con este resultado, el problema de encontrar soluciones
homoclı́nicas para el sistema (3.74) se reduce al de buscar soluciones acotadas
en R, diferentes de la trivial de este sistema.
³ ´
φ
Como las soluciones estacionarias de (3.74) son de la forma 00 donde
φ0 es raı́z de la ecuación
φ − V+ ḡ(φ) = 0 .
De la condición V+ ḡ 0 (φ) < V1+ , se concluye que esta ecuación tiene una
única solución φ0 , que corresponde a una única solución estacionaria de
(3.74).
110
Para determinar si el punto estacionario es un punto silla, estudiaremos
la “parte lineal” del sistema en un entorno del punto estacionario.
Calculando la derivada del lado derecho de (3.74) en el punto estacionario
y haciendo algunos cálculos obtenemos que la “parte lineal” del sistema en
un entorno del punto estacionario es
à !

d φ

=
dt ψ
0

1
1
(V+ ḡ 0 (φ0 )
c2 −1
− 1)
c
(2
c2 −1


− V+ ḡ 0 (φ0 ))
à !
φ
.
ψ
(3.76)
Su ecuación caracterı́stica es :
·
¸
c
1
λ λ+ 2
(V+ ḡ 0 (φ0 ) − 2) + 2
(V+ ḡ 0 (φ0 ) − 1) = 0 .
c −1
c −1
(3.77)
Las raı́ces de (3.77) son
λ± =
c
(2
c2 −1
− V+ ḡ 0 (φ0 )) ±
q
( c2c−1 (2 − V+ ḡ 0 (φ0 ))2 +
1
(V+ ḡ 0 (φ0 )
c2 −1
2
− 1))
.
(3.78)
Ahora, como V+ ḡ (φ0 ) < 1 y 0 < c < 1, obtenemos 0 < λ+ y λ− < 0.
0
³ ´
φ
Por lo tanto 00
es un punto silla de (3.76) para cualquier valor de
0 < c < 1.
Según el lema 6, la existencia de soluciones homoclı́nicas de (3.74)
es
³ ´
φ0
equivalente a la existencia de soluciones acotadas diferentes de la trivial 0 .
Por el teorema 4, si se denota por

Ac = 

0
1
1
1−c2
−2c
1−c2



(3.79)
y G(t, c) la correspondiente función de Green principal, entonces las soluciones acotadas en R de (3.74) coinciden con las soluciones acotadas de la
ecuación integral
Ã
!
V+ Z ∞
0
φ
G(t − τ, c)
dτ .
=
1 − c2 −∞
cḡ 0 (φ(τ )ψ(τ ) − ḡ(τ )
ψ
à !
111
(3.80)
Ahora estudiaremos bajo qué condiciones la ecuación integral (3.80) tiene
soluciones acotadas.
Lema 7 La función de Green principal de Ac es




1 − c c2 − 1


t


1 c−1


e 

 ,

2




−1
1
+
c

si t > 0




1 + c 1 − c2


t




− 12 e c+1 

 ,




1
1−c
si t < 0 .
G(t, c) = 
(3.81)
Demostración.
Según (3.62), requerimos calcular eAc t y los proyectores P+ , P− .
Los valores propios de Ac son
1
1
,
.
(3.82)
1+c
c−1
Usando la forma canónica de Jordan de Ac , y realizando algunos cálculos,
que no escribiremos aquı́, obtenemos

eAc t
1 t 
= e 1+c 
2
1 + c 1 − c2
1
1−c


1 t 

 + e c−1 
2
1 − c c2 − 1
−1


 .
(3.83)
1+c
Los proyectores son
à !
Pc+
x
1
= 
2
y
à !
Pc−


1
x
= 
y
2
1 + c 1 − c2
1
112
à !
 x
,

y
1−c
1 − c c2 − 1
−1

1+c
(3.84)

à !
 x
.

y
(3.85)
Sustituyendo (3.83), (3.84) y (3.85) en (3.62), se obtiene




1 − c c2 − 1


t


1 c−1 


e 
 ,

2




−1
1
+
c

si t > 0



1+c


t 


1 c+1

−2e 




1
si t < 0 .
G(t, c) = 
1 − c2
(3.86)


 ,
1−c
Lema 8
−|t|
| G(t) |≤ (1 + c)e 1+c .
(3.87)
Demostración.
Se concluye del lema 7 y del hecho que la norma del operador


M (x, y) = 
a11 a12

à !
 x

a21 a22
y
cumple
| M |≤ 2K ,
donde
K=
max
1≤i≤2, 1≤j≤2
| aij | .
Para obtener estimados más simples –sin perder generalidad– del integrando en la ecuación integral, supondremos
1
con α > 0, β > 0 .
(3.88)
1 + e−αu+β
Obteniendo la primera y segunda derivada de ḡ, definida por (3.88) se
tiene
ḡ(u) =
ḡ 0 (u) = αḡ(1 − ḡ) ,
113
(3.89)
ḡ 00 (u) = αḡ 0 (1 − 2ḡ) ,
(3.90)
| ḡ 0 |≤ α ,
(3.91)
| ḡ 00 |≤ 2α2 .
(3.92)
Lema 9 Sea ρ > 0. Si
| φ |2 + | ψ |2 ≤ ρ2 ,
| φi |2 + | ψi |2 ≤ ρ2 para 1 ≤ i ≤ 2 .
Entonces, se tienen las siguientes desigualdades
| cḡ 0 (φ)ψ − ḡ(φ) |≤ cαρ + 1 ,
(3.93)
| (cḡ 0 (φ1 )ψ1 − ḡ(φ1 )) − (cḡ 0 (φ2 )ψ2 − ḡ(φ2 )) |≤ 2α(αρc + 1) .
(3.94)
Demostración.
Se sigue inmediatamente de (3.91) y (3.92).
Proposición 4 Supongamos que en el sistema (3.74) la función ḡ tiene la
forma
1
ḡ(φ) =
, α > 0, y β > 0 .
1 + e−αφ+β
1
}
Supongamos además que V+ α < 14 . Entonces para cada ρ > max{φ0 , 2α
existe 0 < c(ρ) < 1 tal que para cada 0 < c < ³c(ρ),
la única solución acotada
´
φ
de este sistema contenida en la bola cerrada { ψ ∈ R2 :| φ |2 + | ψ |2 ≤ ρ2 },
³
es la solución estacionaria
³
c(ρ) =
1+
1
αρ
+
1
4α2 ρV+
´
´
φ0
0 . Además c(ρ) viene dada por
+
r³
1+
1
αρ
+
2
114
1
4α2 ρV+
´2
−4
³
1
αρ
−
1
4α2 ρV+
´
. (3.95)
Observación.
³ ´
φ
La condición V+ α < 14 garantiza que 00 es el único punto estacionario
del sistema (3.74), donde φ0 es la única solución de la ecuación
φ − V+ ḡ(φ) = 0 .
Demostración.
Demostraremos que bajo las condiciones de esta proposición, el sistema
(3.74) cumple con las condiciones del teorema 4, y que por lo tanto este
sistema
en la bola cerrada
³ ´ debe tener una única solución acotada u
³ contenida
´
φ
φ0
2
2
2
2
{ ψ ∈ R :| φ | + | ψ | ≤ ρ }, pero como 0 es también una solución
acotada
del sistema contenida en la misma bola cerrada, se debe tener que
³ ´
φ0
0 es la única solución del sistema (3.74) en esta bola cerrada.
Pasemos a demostrar que el sistema (3.74) cumple con las hipótesis del
teorema 5.
Los valores propios de la matriz Ac son
1
,
1+c
1
,
c−1
por lo tanto Ac cumple con la condición (3.64).
Por el lema 9 se tiene
V+
1 − c2
¯Ã
¯Ã !¯
!¯
¯
¯ φ ¯
¯
V+
0
¯
¯
¯
¯
(cαρ
+
1)
para
¯
¯
¯≤
¯ ≤ ρ,
0
2
¯ cḡ (φ)ψ − ḡ(φ) ¯
¯
1−c
ψ ¯2
V+
1 − c2
¯Ã
!¯
¯
¯
0
¯
¯
¯
¯
0
0
¯ cḡ (φ1 ) − ḡ(φ1 ) − (cḡ (φ2 )ψ2 − ḡ(φ2 )) ¯
¯Ã
!
Ã
¯ φ
V+
φ2
1
¯
≤
2α(αρc
+
1)
¯
−
2
¯
1−c
ψ1
ψ2
para
¯Ã !¯
¯Ã !¯
¯ φ ¯
¯ φ ¯
1 ¯
2 ¯
¯
¯
¯
¯ ≤ρ y ¯
¯ ≤ ρ.
¯ ψ1 ¯
¯ ψ2 ¯
2
2
115
!¯
¯
¯
¯
¯
2
2
Ası́, también se cumplen las condiciones (3.66), (3.67), con
M = (cαρ + 1)
V+
,
1 − c2
q = 2α(αρc + 1)
V+
.
1 − c2
Por los lemas 7 y 8, se tiene que la función de Green principal, G(t, c),
de Ac cumple
−|t|
| G(t, c) |≤ (1 + c)e 1+c ,
por lo tanto, también se cumple (3.67) con
N =1+c y r =
1
.
1+c
Sólo nos falta probar la condición (3.68), es decir, nos falta probar que
bajo las condiciones de la proposición 4 se cumple
(cαρ + 1)V+
ρ
≤
,
2
1−c
2(1 + c)2
(3.96)
2α(αρc + 1)V+
1
<
.
2
1−c
2(1 + c)2
(3.97)
Multiplicando (3.96), (3.97) por 1−c2 y simplificando obtenemos que este
sistema de desigualdades es equivalente a
(cαρ + 1)V+ ≤
ρ(1 − c)
,
2(1 + c)
2α(αρc + 1)V+ <
(3.98)
1−c
.
2(1 + c)
(3.99)
Realizando algunas operaciones y reordenando obtenemos que el sistema
(3.98), (3.99) es equivalente al sistema
Ã
!
1
1
1
1
P1 (c) = c2 + 1 +
+
c+
−
≤ 0,
αρ 2V+ α
αρ 2V+ α
Ã
!
Ã
1
1
1
1
+ 2
c+
− 2
P2 (c) = c + 1 +
αρ 4α ρV+
αρ 4α ρV+
2
116
(3.100)
!
< 0.
(3.101)
Ası́ que necesitamos probar que
P1 (c) ≤ 0 ,
P2 (c) < 0 .
(3.102)
(3.103)
Para probar esto, obtengamos las raı́ces de los polinomios cuadráticos
P1 (c), P2 (c).
Las raı́ces de P1 son
³
λ1± (ρ) =
− 1+
1
αρ
+
1
2V+ α
´
±
r³
1
αρ
1+
+
1
2V+ α
´2
−4
³
1
αρ
−
1
2V+ α
2
´
.
(3.104)
Las raı́ces de P2 son
³
λ2± (ρ) =
− 1+
1
αρ
+
1
2V+ α
´
±
r³
1+
1
αρ
+
1
4ρV+ α2
´2
−4
2
Como por hipótesis V+ α <
1
4
yρ>
1
2α
³
1
αρ
−
´
1
4α2 ρV
+
.
(3.105)
se tiene
1
1
− 2
< 0,
αρ 4α ρV+
por lo tanto las raı́ces de P1 (c) son reales, una negativa, λ1− (ρ), y la otra
positiva λ1+ (ρ). Lo mismo ocurre con P2 (c), tiene una raı́z negativa, λ2− (ρ) y
otra raı́z positiva, λ2+ (ρ).
Entonces
P1 (c) ≤ 0 para λ1− (ρ) ≤ c ≤ λ1+ (ρ) ,
(3.106)
P2 (c) < 0 para λ1− (ρ) < c < λ2+ (ρ) .
(3.107)
Además
λ2+ (ρ) < λ1+ (ρ) ,
117
(3.108)
esto último se sigue de la expresión para cada uno de ellos y del hecho que
4α2 ρV+ = 2αρ(2V+ α) > (2V+ α) ya que ρ >
1
.
2α
También es fácil probar que λ1+ (ρ) < 1. Ası́, si tomamos 0 < c(ρ) =
< 1 se tiene que, para cualquier 0 < c < c(ρ)
λ1+ (ρ)
P1 (c) ≤ 0 ,
P2 (c) < 0 ,
que es lo que deseábamos demostrar.
Resumiendo, con las condiciones de la proposición 4, el sistema (3.74)
cumple con las condiciones del teorema 5, para cualquier 0 < c < c(ρ). Esto
es lo que querı́amos demostrar.
1 },
Con la proposición 4, hemos demostrado que si V+ α < 14 , ρ > max{φ0 , 2α
entonces el sistema (3.74) no puede tener soluciones homoclı́nicas para 0 <
c < c(ρ), por lo tanto, la no existencia de ondas viajeras para la ecuación
(3.59), cuando tomamos 0 < c < c(ρ).
1 }
Ahora pasaremos a analizar qué sucede cuando tomamos ρ > max{φ0 , 2α
y c(ρ) < c < 1.
1 } el espacio β ,
Para esto, consideremos, para ρ > max{φ0 , 2α
ρ
(Ã !
βρ =
à !
φ
φ
: R −→ R2 :
ψ
ψ
)
es continua, k φ
k2∞
+kψ
k2∞ ≤
ρ
2
.
Sabemos que βρ es un
de Banach. Consideremos la familia de
µ³ espacio
´ ¶
φ
operadores integrales F ψ , c definida por
ÃÃ !
F
!
à !
φ
φ
,c =
− F1
ψ
ψ
ÃÃ !
!
φ
,c ,
ψ
(3.109)
donde
ÃÃ !
F1
!
Ã
!
V+ Z ∞
φ
0
,c =
G(t − τ, c)
dτ .
cḡ 0 (φ(τ ))ψ(τ ) − ḡ(φ(τ ))
ψ
1 − c2 −∞
(3.110)
118
µ³ ´ ¶
φ
ψ , c actúe de βρ a βρ , c debe cumplir con la condición
Para que F
(3.100). De aquı́ concluı́mos que aquellos valores de c para los cuales es posible encontrar soluciones homoclı́nicas del sistema (3.74), son los que cumplen
λ2+ (ρ) < c < λ1+ (ρ) .
Ası́ consideremos
F : βρ × (λ2+ (ρ), λ1+ (ρ)) −→ βρ definido por (3.109) .
³
φ
´
Como 00 es un punto estacionario del sistema (3.74) para cada 0 < c <
1, tenemos que
ÃÃ !
F
!
φ
, c = 0 para cada λ2+ (ρ) < c < λ1+ (ρ) .
0
³³ ´
´
Entonces, si lograramos probar que DF φ00 , c es invertible, y que se
cumplen ³³
las otras
condiciones del teorema de la función implı́cita
³ en
´ un en´ ´
φ0
φ0
torno de
, c , se concluirı́a que existen U, V vecindades de 0 y c res0
pectivamente, tal que para cada c∗ ∈ V , existe una única L (c∗ ) =
tal que
ÃÃ !
³ ´
φ
ψ
∈U
!
φ
, c∗ = 0 ,
ψ
F
³ ´
φ
pero entonces, ψ serı́a una solución acotada en R, no trivial del sistema
(3.74) y por lo tanto una solución homoclı́nica del mismo sistema, a la cual
le corresponderı́a una onda viajera de la ecuación (3.59) con velocidad c∗ .
Ası́ que, en lo que sigue probaremos que el operador integral F con
λ1+ (ρ) < c < λ2+ (ρ), cumple con las condiciones del teorema de la función
implı́cita.
³³ ´ ´
Vamos ahora a calcular la derivada de Frechet de F en ψφ , c con
λ1+ (ρ) < c < λ2+ (ρ). Por la forma que tiene F, basta calcular la derivada
de Frechet de la función F1 . Además se tendrı́a que
ÃÃ !
DF
!
φ
, c = I − DF1
ψ
119
ÃÃ !
!
φ
,c .
ψ
(3.111)
³ ´
φ
Para hacer esto, calcularemos las derivadas de Gato de F1 en ψ con res³³ ´ ´
pecto a φ y ψ respectivamente. Denotaremos estas derivadas por δφ F1 ψφ , c ,
δψ F2
³³ ´
φ
ψ
´
, c respectivamente.
Lema 10 Bajo las condiciones de la proposición 4, el³³operador
´ ´ integral F
definido en (3.110), es de clase C 1 en un entorno de φ00 , c . Además la
derivada parcial de Frechet de F con respecto a
DF
³³ ´
φ
ψ
³ ´
φ
ψ
´
, la cual denotamos por
, c , viene dada por
ÃÃ !
DF
con
ÃÃ !
!
φ
, c = I − DF1
ψ
!
φ
, 0 = δφ F1
ψ
DF1
ÃÃ !
ÃÃ !
!
φ
,c ,
ψ
!
φ
, c + δψ F1
ψ
ÃÃ !
φ
,c
ψ
!
(3.112)
y
ÃÃ !
!
φ
,c ,
ψ
δφ F1
ÃÃ !
δψ F1
φ
,c
ψ
!
son las derivadas de Gato, con respecto a φ y ψ, respectivamente, las cuales
vienen dadas por
δφ F1
³³ ´
φ
ψ
´
, c (h) =

"
!
#
Ã
t
−τ
c+1
−V+  e c+1 Z ∞ c+1
e
[ cḡ 00 (φ(τ ))ψ(τ ) − ḡ 0 (φ(τ ) ]h(τ ) dτ
2 1 + c t
1
"
Ã
!
#


−τ
e c−1 Z t
c−1
+
e c−1
[ cḡ 00 (φ(τ ))ψ(τ ) − ḡ 0 (φ(τ )) ]h(τ ) dτ  ,
c + 1 −∞
1
t
ÃÃ !
δψ F1

!
Ã
(3.113)
!
−τ
−V+ c  e c+1 Z ∞ c+1
φ
c+1 0
, c (k) =
e
ḡ (φ(τ ))k(τ ) dτ
ψ
2 1 + c t
1
+
t
c−1
e
1+c
Z t
−∞
Ã
−τ
e c−1
t
!


c−1 0
ḡ (φ(τ ))k(τ ) dτ .

1
120
(3.114)
Demostración.
Observemos que
F1
³³ ´
φ
ψ
´
,c =
 t
"
Ã
!
#
−τ
−V+  e c+1 Z ∞ c+1
c+1
0
e
cḡ (φ(τ ))ψ(τ ) − ḡ(φ(τ )) dτ
2 1 + c t
1
"
Ã
!

#

−τ
e c−1 Z t
c−1
+
e c−1
(cḡ 0 (φ(τ ))ψ(τ ) − ḡ(φ(τ )) dτ  .
c + 1 −∞
1
t
Calculemos la derivada de Gato de F1 en el punto
³³ ´
φ
ψ
aφ
ÃÃ !
δφ F1



!
F1
φ
, c (h) = lim 
²→0
ψ

Ã
−τ
c+1
V+
e c+1 Z ∞ c+1

=−
lim
e

²→0
2
1+c t
1
t
³³
´
φ+²h
ψ
´
, c − F1
Ã
³³ ´
φ
ψ

dτ
(ḡ 0 (φ(τ ) + ²h(τ )) − ḡ(φ(τ )))ψ(τ )
c
²
ḡ(φ(τ )) + ²h(τ ) − ḡ(φ(τ ))
−
²
δφ F1
´
!#
!" Ã
Ã
Entonces
,c
(ḡ 0 (φ(τ ) + ²h(τ )) − ḡ 0 (φ(τ )))ψ(τ )
c
²
ḡ(φ(τ ) + ²h(τ ) − ḡ(φ(τ ))
−
²
t
φ
ψ
!" Ã
Ã
−τ
e c+1 Z t c+1
c−1
+
e
c + 1 −∞
1
´
, c con respecto
³³ ´
²
(3.115)
´
, c (h) =
121
!#
#)
dτ
.
!
!

Ã
!
#
"
t
−τ
V+  e c+1 Z ∞ c+1
c+1
00
0
−
[cḡ (φ(τ ))ψ(τ )h(τ ) − ḡ (φ(τ ))h(τ )] dτ
e
1
2 1 + c t
+
t
1+c
e
1+c
Z t
"
Ã
De manera similar, se obtiene la derivada de Gato de F1 en
respecto a ψ
ÃÃ !
δψ F1


#
c−1
[cḡ 00 (φ(τ )ψ(τ )h(τ ) − ḡ 0 (φ(τ ))h(τ )] dτ  .
1
−τ
e c−1
−∞
!

!
Ã
³³ ´
φ
ψ
´
, c con
!
−τ
φ
V+ c  e c+1 Z ∞ c+1
c+1 0
, c (k) =
e
ḡ (φ(τ ))k(τ )dτ
ψ
2 c + 1 t
1
t

!
Ã

−τ
c−1 0
e c−1 Z t c−1
ḡ (φ(τ ))k(τ )dτ .
+
e

1
1 + c −∞
t
Además, δφ F1
³³ ´
φ
φ
´
, c , δψ F1
³³ ´
φ
φ
´
, c son lineales como función de h y k
³ ´
φ
respectivamente y continuas como función de ψ , se tiene que la derivada de
³ ´
φ
Frechet de F1 con respecto a ψ existe y se puede calcular con la expresión
(3.112).
³³ ´ ´
También es fácil probar que DF ψφ , c es continua en un entorno de
(
³ ´
φ
ψ ,c
).
Lema 11 Bajo las condiciones de la proposición 4, DF
erador invertible.
Demostración.
Como
ÃÃ
DF
!
!
φ0
, c = I − DF1
0
122
ÃÃ
!
!
φ0
,c ,
0
³³ ´
φ0
0
´
, c es un op-
basta probar que
Evaluando (3.112) en
ÃÃ
DF1
!
φ0
,c
0
°
ÃÃ ! !°
°
°
φ0
°
°DF1
, c °° < 1 .
°
°
0
³³ ´ ´
φ0
0
!Ã !
h
= δφ F1
k
, c obtenemos
ÃÃ
!
!
φ0
, 0 (h) + δψ F1
0
ÃÃ
!
!
φ0
, c (k) , (3.116)
0
entonces
°
ÃÃ ! ! Ã !°
°
φ0
h °°
°
°DF1
,c
°
°
0
k °
∞
°
°
°
°
ÃÃ ! !
ÃÃ ! !
°
°
°
°
φ0
φ0
°
°
°
≤ °δφ F1
, c (h)° +°δψ F1
, c (k)°°
°
°
°
°
0
0
∞
.
∞
(3.117)
Obtengamos estimados para
°
ÃÃ ! !°
°
°
φ0
°
°δφ F1
, c °°
°
°
0
,
∞
Evaluando (3.113) en
δφ F1
³³ ´
φ0
0
´
³³ ´
φ0
0
°
ÃÃ ! !°
°
°
φ0
°
°δψ F1
, c °°
°
°
0
.
∞
´
, c obtenemos
, c (h) =

!
Ã
Ã
!


−τ
−τ
e c−1 Z t c−1
c+1
c−1
V+ ḡ(φ0 )  e c+1 Z ∞ c+1
h(τ )dτ +
h(τ )dτ
−
e
e

c + 1 t
1
1
2
1 + c −∞

t
Ã
t
!
Ã
!


−τ
−τ
c+1
e c−1 Z t c−1
c−1
V+ cḡ 0 (φ0 )  e c+1 Z ∞ c+1
k(τ )dτ −
k(τ )d(τ ) ,
+
e
e
c + 1 t

2
1
c + 1 −∞
1
t
t
entonces
°
°
ÃÃ ! !
°
°
φ
°
°
,
c
(h)
°
°δφ F1
°
°
0
∞
√

t
−τ
V+ ḡ (φ0 ) √ 2
e c+1 Z ∞ c+1
≤
c
+
2c
+
2
e
dτ k h k∞

2
c+1 t
0
+ c2 − 2c + 2
t
c−1
e
1+c
Z t
−∞
123
τ
e c−1 dτ k h k∞



.
Calculando las integrales impropias y simplificando
°
°
ÃÃ ! !
°
°
φ0
°
°δφ F1
, c (h)°°
°
°
0
≤
∞
´
√
V+ ḡ(φ0 ) ³√ 2
c + 2c + 2 + c2 − 2c + 2 k h k∞ .
2
Para 0 < c < 1 se tiene que
√
√
c2 + 2c + 2 es creciente,
c2 + 2c − 2 es decreciente.
Ası́ se tiene que
°
°
ÃÃ ! !
°
°
φ0
°
°δφ F1
, c (h)°°
°
°
0
√ ´
V+ ḡ 0 (φ0 ) ³√
5 + 2 k h k∞ .
2
(3.118)
√ ´
V+ ḡ 0 (φ0 )c ³√
5 + 2 k k k∞ .
2
(3.119)
≤
∞
De manera similar se obtiene
°
°
ÃÃ ! !
°
°
φ0
°
°
°δψ F1
, c (k)°
°
°
0
≤
∞
Sustituyendo (3.118), (3.119) en (3.117), obtenemos
°
ÃÃ ! ! Ã !°
¶
°
³√
√ ´ µq
φ0
h °° V+ ḡ 0 (φ0 )
°
2
2
(1 + c) 5 + 2
k h k∞ + k k k∞ .
°DF1
,c
°≤
°
2
k °
0
Como 0 < c < 1 y V+ α <
1
4
obtenemos
°
ÃÃ ! ! Ã !°
°
φ0
h °°
°
°DF
,c
° < 1,
°
k °
0
que es lo que querı́amos probar.
Finalmente, tenemos el siguiente teorema.
124
Teorema 6 Supongamos que en el sistema (3.74), la función ḡ tiene la
forma
1
ḡ(φ) =
, α, β > 0
1 + e−αφ+β
1
y que αV+ < 14 . Entonces para cada ρ > max{φ0 , 2α
}, existen c1 (ρ), c2 (ρ)
con 0 < c1 (ρ) < c2 (ρ) < 1, tales que, para cada c, que cumpla 0 < c < c1 (ρ),
la única solución acotada de³este
sistema, contenida en la bola cerrada β[0, ρ],
´
φ0
es la solución estacionaria 0 . Mientras que, para cada c, con c1 (ρ) ≤ c ≤
c2 (ρ), existe una solución acotada, no trivial, del sistema, la cual también es
una solución homoclı́nica, y a la que le corresponde una onda viajera de la
ecuación (3.59) con velocidad v. La velocidad real de traslación de esta onda
viajera es cV+ . Además c1 (ρ), c2 (ρ) son las raı́ces positivas de los polinomios
cuadráticos definidos en (3.100), (3.101).
Demostración.
La existencia de c1 (ρ), fue demostrada en la proposición 4. Por los lemas
10 y 11 el operador integral F definido en (3.101) cumple con
³³las´ condi´
ciones del teorema de la función implı́cita en un entorno de φ00 , c , es
decir, F
³³ ´
φ
ψ
´
, c es de clase C 0 en un entorno de
³³ ´
φ0
0
´
,0 , F
³³ ´
φ0
³0
´
,c = 0 y
´
φ0
0 ³ y´ c resφ
pectivamente tales que, para cada c∗ ∈ V , existe un único elemento ψ ∈ U
DF
³³ ´
φ0
0
´
, c es invertible. Entonces existen entornos U y V de
no trivial tal que
ÃÃ !
F
!
φ
, c∗ = 0 .
ψ
³ ´
φ
Entonces ψ es una solución no trivial, acotada de (3.80) y por el teorema
³ ´
φ
3, ψ es una solución acotada, no trivial del sistema (3.74), la cual, por el
lema 6, es una solución homoclı́nica del sistema (3.74).
3.4 ONDAS ESTACIONARIAS EN LA CORTEZA
CEREBRAL.
125
En esta sección estudiaremos las oscilaciones pequeñas –con la misma
frecuencia angular w para cada x ∈ R– de la ecuación
∂2u
∂t2
2
2 ∂ u
+ 2λ+ v+ ∂u
+ (λ+ v+ )2 u − v+
=
∂t
∂x2
h
V+ λ+ v+ λ+ v+ ḡ(u) + ḡ 0 (u) ∂u
∂t
i
(3.120)
alrededor de un estado de reposo de la misma (una solución u(x, t) ≡ u0 ).
Este problema corresponde a considerar una banda circular de ancho
infinitamente pequeño de corteza cerebral. Lo que queremos saber es si
en nuestro modelo existen ondas estacionarias que se correspondan con los
ritmos básicos del cerebro.
Para estudiar este problema consideraremos una solución u(x, t) ≡ u0 de
la ecuación (3.120). Tomaremos la parte lineal de esta ecuación alrededor de
u0 y en la ecuación lineal buscaremos soluciones de la forma
u(x, t) = u0 + L(x)eiwt ,
con L(x) una función periódica real, de periodo 2π y L(x) pequeño.
126
(3.121)
Observación.
En realidad, lo que buscamos son soluciones de la ecuación lineal de la
forma
u(x, t) = u0 + A(x) cos wt
donde A(x) es una función periódica real, de periodo 2π. Pero como se tiene
una ecuación lineal, basta encontrar soluciones de la forma (3.121) y después
tomar la parte real de estas soluciones.
Ahora, si u(x, t) ≡ u0 es una solución de (3.120) la “parte lineal” de esta
ecuación alrededor de u0 es
2
∂ 2 ū
∂ ū
0
2
0
2 ∂ ū
+
λ
v
(2
−
V
ḡ
(u
))
+
[(λ
v
)
(1
−
V
ḡ
(u
))]ū
−
v
=0
+
+
+
o
+
+
+
0
+
∂t2
∂t
∂x2
(3.122)
donde
ū = u − u0
Teorema 7 Sea u0 ∈ (0, V+ ) tal que, ḡ 0 (u0 ) = V2+ , supongamos además que
1 − λ2+ > 0. Entonces la ecuación (3.122) tiene soluciones de la forma
ū(x, t) = L(x)eiwt ,
con L(x) una función periódica real, de periodo 2π. Más concretamente, si
µq
w = v+
¶
1−
λ2+
,
(3.123)
ū(x, t) = A1 cos xeiwt
(3.124)
entonces
con A1 constante, es una solución de (3.122).
127
Observación
2
,
v+
Como ḡ 0 (u0 ) =
esto nos garantiza que
u(x, t) ≡ u0 es solución de la ecuación (3.120).
Demostración.
Sea w como en (3.123). Obtengamos las derivadas de (3.124)
∂ ū
∂t
= (iw)A1 cos xeiwt ,
∂ 2 ū
∂t2
= −w2 A1 cos xeiwt ,
∂ ū
∂x
= −A1 sen xeiwt ,
∂ 2 ū
∂x2
= −A1 cos xeiwt .
Sustituyendo estas derivadas en la ecuación (3.122) obtenemos
−w2 A1 cos xeiwt + λ+ v+ (2 − V+ ḡ 0 (u0 ))iwA1 cos xeiwt
2
+[(λ+ v+ )2 (1 − V+ ḡ 0 (u0 )]A1 cos xeiwt + v+
A1 cos xeiwt =
2
2
(λ2+ − 1)v+
A1 cos xeiwt + (λ+ v+ )2 A1 cos xeiwt + v+
A1 cos xeiwt = 0 ,
ya que ḡ 0 (u0 ) =
2
V+
q
y w = v+ 1 − λ2+ .
128
CONCLUSIONES.
El principal modelo que hemos estudiado, a pesar de ser un modelo que
tiene algunas simplificaciones importantes : modelo unidimensional sin curvatura, considerar una distribución de conexiones que depende únicamente de
la distancia entre columnas y con escalas de decrecimiento que no dependen
de las regiones corticales consideradas, las complicaciones intra e intercorticales están incluidas en las funciones R± , no tomar en cuenta explı́citamente
la actividad neuroquı́mica, etc., nos ha permitido obtener un conjunto de
resultados que pueden ser sujetos a ser estudiados experimentalmente y que
en principio coinciden cualitativamente con algunos datos experimentales ya
conocidos. Se debe destacar la relativa simplicidad del modelo obtenido y
la posibilidad inmediata de incluir, aunque sea parcialmente, algunos elementos no considerados en el modelo principal como por ejemplo, escalas de
decrecimiento que dependan de las regiones de corteza considerada, que la
función de generación de potenciales de acción ḡ dependa de la columna y
de la función de activación (ver el modelo de la sección 2.8) además de la
posibilidad de construir un modelo matemático bidimensional de la corteza
cerebral que además incluya, aunque sea parcialmente, la geometrı́a real de la
corteza cerebral. Obviamente el estudio matemático de este problema resulta
mucho más complejo. Consideramos que todo lo anteriormente enunciado,
son los elementos más valiosos de la tesis.
La inclusión en el modelo de la actividad neuroquı́mica nos parece que
es un problema sumamente complejo pero en primera aproximación, considerando que la liberación de neurotransmisores es proporcional a la actividad sináptica, podrı́a ser incluida en el mismo.
A continuación enlistamos las principales conclusiones de nuestro trabajo.
1. En este trabajo se han continuado desarrollando las ideas de Nuñez,
[86], [87], [89], [92], en al menos dos direcciones : se han introducido
modificaciones y agregado nuevos elementos al modelo. Además el
tratamiento matemático del modelo ha sido diferente y más completo
que el que le dió Nuñez.
2. En cuanto a resultados, se ha demostrado, que bajo ciertas condiciones,(teorema 6) existen ondas viajeras de actividad cortical que se
propagan por la corteza cerebral con una velocidad de propagación
menor que la velocidad de propagación de la actividad excitatoria.
129
También se ha demostrado,(teorema 2) que la actividad cortical de una
rebanada de corteza se puede acoplar con la de una entrada aferente
periódica sostenida. Esto ha sido demostrado experimentalmente [25],
[78]. Además en estos trabajos también se concluye, experimentalmente, el papel controlador de la actividad inhibitoria en la corteza.
Lo cual probamos en el teorema 3.
3. Los resultados obtenidos,(teorema 2, 6, 7) corresponden a estados fisiológicos muy especiales, ésto último se refleja en las restricciones
impuestas a la función de activación neuronal. Consideramos que es
necesario estudiar, posteriormente, la significación fisiológica de estas
restricciones.
4. Hemos obtenido un modelo (sección 2.8), que en conjunto con estudios
del problema inverso electroencefalográfico, nos permitirá estudiar, posteriormente, el diagrama de conexiones funcionales y de la función de
activación neuronal en un estado fisiológico determinado.
5. El modelo compartimentado (2.8) puede ser usado para el estudio de
otros sistemas del cerebro, como el sistema de los ganglios basales y el
sistema lı́mbico.
6. Posteriormente será necesario realizar simulaciones computacionales de
las soluciones del modelo con diversos valores de los parámetros involucrados y compararlos con mediciones experimentales.
130
Bibliografı́a
[1] Abbott, L. F. [1990] : A network of oscillators, J. Phys. A. Math.
Gen. 23 : 3835-3859.
[2] Abbot, L. F., F., Kepler, T. B. [1990] : Model neurons from
Hodgkin-Huxley to Hopfield In : Garrido (ed) Statistical mechanics of neural networks. Springer, Berlin, Heidelberg, New York.
[3] Abeles, M. [1982] : Local Cortical Circuits in Studies in Brain
Function 6, Springer Verlag, New York.
[4] Amick, C, J, [1991] : A problem in neural networks. Proceedings
of the Royal Society of Edinburgh. 118 : 225-236.
[5] Amari, S. [1972] : Competitive and cooperative aspects in dynamics of neural excitation and self-organization, IEEE Trans on
computers.
[6] Amari, S. [1982] : A mathematical theory of self-organizing nerve
sytems. In Proceedings of Biomathematics : Current Status and
Perspectives. Ricciardi, L. M., Scott, A., ed.
[7] Amari, S., Arbib, M. A (eds) [1982] : Competition and Cooperation in Neural Networks, Springer Verlag, New York.
[8] Amari, S., Takeuchi, A. [1978] : Mathematical theory on formation
of category detecting nerve cells. Biol. Cybernetics 29. 127-136.
[9] Amir, A. [1994] : Uniqueness of the generators of brain evoked
potential maps, IEEE Trans. on Biomedical-Engineering. Vol.
41,no.1-11.
131
[10] Andersen, P., Andersson, S.A. [1968] : Physiological basis of the
alpha rhytm. Appleton-Century-Crofts, New York.
[11] Anderson, J. A., Rosenfeld, E. [1988] : Neurocomputing Foundation of Research. Cambridge. The M.I.T. Press.
[12] Beale, R.,Jackson, T. [1994] : Neural Computing : An Introduction, Institute of Physics Publishing, Bristol and Philadelphia.
[13] Baumgärtel, H. [1985] : Analytic Perturbation Theory for Matrices
and Operators. Birkhäuser Verlag Basel. Boston-Stuttgart.
[14] Bogacz, J. y colaboradores. [1985] : Los potenciales evocados en
el hombre. Ateneo, Buenos Aires.
[15] Bower, J. M. [1992] : Modeling the nervous system. Trends in
Neuroscience 15 : 411-412.
[16] Buchholtz, F., Golowasch, J., Epstein, I. R., Marder, E. [1992]
: Mathematical model of an identified stomatogastric neuron.
J.Neurophysiol. 67 : 332-340.
[17] Cardinali, D. P. [1992] Manual de Neurofisiologı́a. Ediciones Dı́az
de Santos, S.A. Madrid, España.
[18] Chauvet, G. A. [1993a] : Non-locality in biological systems results
from hierarchy. Application to nervous system. Journal Math. Biol.
31. 475-486.
[19] Cichocki, A., Unbehauen, R. [1993] : Neural Networks for Optimazation and Signal Processing, John Wiley, New York.
[20] Coro, F. [1996] : Fisiologı́a Celular, Un enfoque biofı́sico, Textos
Cientı́ficos , BUAP, Puebla.
[21] Cronin, J. [1987] : Mathematical aspects of Hodgkin-Huxley neural
theory, Cambridge University Press.
[22] Cuffin, B.N. [1987] : The role of model and computacional experiments in the biomagnetic inverse problem.Phys.Med.Biol.
132
[23] Daleckii, Ju. L., Krein, M. G. [1974] : Stability of Solutions of
Differential Equations in Banach Space. American Mathematical
Society, Providence, Rhode Island.
[24] Destexhe, A., Contreras, D., Sejnowski, T. J., Steriade, M. [1994]
: A model of spindle rhytmicity in the isolated thalamic reticular
nucleus. Journal of Neurophysiology. Vol. 72, No. 2. 803-1818.
[25] Dickson, C. T., Alonso, A [1997] : Muscarinic induction of Synchronous population activity in the entorhinal cortex. Journal of
Neuroscience. 17 (17:6729-6744).
[26] Doya, K., Selverston, A.I. [1993] : Dimension reduction of biological neuron models by artificial neural networks, Neural Comp. Vol.
6 no. 4 : 697-717.
[27] Farkas, Miklos. [1994] : Periodics Motions, Springer Verlag, New
York.
[28] Fitzhugh, R. [1961] : Impulses and physiological stages in theoretical models of nerve membrane. Biophys, J. 1 : 445-466.
[29] Fitzhugh, R. [1969] : Mathematical models of excitation and propagation in nerve. In Biological Engineering H.P. Schwann.Mc Graw
Hill Publications, New York.
[30] Flores, G. [1991] : Stability Analysis for the Slow Traveling Pulse
of the Fitzhugh-Nagumo System, Siam J. Math. Anal. Vol. 22, no.
2, 392-399.
[31] Freeman J. A., Skapura, D.M. [1993] : Redes Neuronales. Algoritmos, aplicaciones y técnicas de programación. Addison-Wesley
Iberoamericana Wilmington, Delaware. U.S.A.
[32] Freeman, W.J.[1975] : Mass Action in the Nervous System, Academic Press, New York.
[33] Freeman, W.J. [1979] : Nonlinear dynamics of paleocortex manifested in olfactory EEG, Biological Cybernetics. 35 : 21-37.
133
[34] Freeman, W. J. [1992] : Predictions on neocortical dynamics derived from studies in paleocortex. In : Induced rhythms of the
brain, ed. E.Basar, T. H. Bullock. Birkhauser Boston Inc. 183-201.
[35] Freeman, W. J. [1994] : Neural mechanisms underlying destabilization of cortex by sensory input. Physica D 75 : 151 - 64.
[36] Fukushima, K. [1975] : Cognitron : a self organizing multilayered
neural network, Biolog. Cybernetics. Vol. 20 : 121-136.
[37] Fukushima, K. [1980] : Neocognitron : a self organizing neural
network model for a mechanism of pattern recognition unaffected
by shift in position, Biolog. Cybernetics. Vol. 36,193-202.
[38] Galusinski, C. [1998] : Existence and Continuity of Uniform Exponential Attractors of the Singularity Perturbed Hodgkin-Huxley
System. Journal of Differential Equations. 144 : 99-169.
[39] Golomb, D., Guckenheimer, J., Gueron, S. [1993] : Reduction of a
Channel Based Model for a Stomatogastric Ganglion LP Neuron.
Biol. Cybern. 68 : 129-137.
[40] Golomb, D., Rinzel, J. [1993] : Dynamics of globallly coupled
inhibitory neurons with heterogeneity., Physical Review. E. Vol.
48,no. 6.
[41] Golomb, D., Rinzel, J. [1994] : Clustering in globally coupled inhibitary neurons. Physics D. 72 : 259-282.
[42] Golomb, D., Wang, X. J., Rinzel, J. [1994] : Synchronization
Properties of Spindle Oscillations in a Thalamic Reticular Nucleus
Model. J. Neurophysiology. Vol. 72, no.3, 1109-1126.
[43] Grossberg, S. [1976] : Adaptative pattern classification and universal recording. Biolog. Cybernetics. Vol. 23 : 121-134 y 187-202.
[44] Guckenheimer, J., Holmes, P. [1983] : Nonlinear oscillations dynamical systems, and bifurcations of vector fields. Springer, Berlin,
Heidelberg, New York.
134
[45] Hale, J. K. [1963] : Oscillations in Nonlinear Systems. Dover, New
York.
[46] Hale, J. K. [1969] : Ordinary Differential Equations. Wiley Interscience. New York, London, Sydney, Toronto.
[47] Hale, J. ; Kocak, H. [1991] : Dynamics and Bifurcations, Springer
Verlag, New York.
[48] Hastings, S. [1975/1976] : On traveling wave solutions of the
Hodgkin-Huxley equations, Arch. Rational Mech. Anal. 60 : 229257.
[49] Hastings, S. [1976] : On existence of homoclinic and periodic orbits
for the Fitzhugh-Nagumo equations,Quart J. Math.27 : 123-134.
[50] Hebb, D. O. [1949] : The organization of behavior. Wiley, New
York.
[51] Heiner, L. ; Záborszky, L. [1989] : Neuroanatomical Tract–Tracing
Methods 2. Plenum Press. New York and London.
[52] Heller, L. [1990] : Return Current in encephalography Variationals
Principles, Biophysical Journal. Vol. 57 : 601-606.
[53] Herman G. T. [1980] : Image reconstruction from projections. The
fundamental of computerized tomography. Academic Press.
[54] Herman G. T., Tuy H. K., Langenberg K. J., Sabatier P.C. [1985]
: Basic methods of tomography and inverse problems. Malvern
Physics Series.
[55] Hilera, J. R., Martı́nez, V. J. [1995] : Redes neuronales artificiales.
Fundamentos, modelos y aplicaciones. Addison-Wesley Iberoamericana. Madrid, España.
[56] Hinton, G. Ackley, D. Sejnowski, T. [1984] : Boltzman machines
: Constraint satisfaction networks than learn. Carnegie-Mellon
University, Departament of Computer Science Technical Report.
(CMU-CS-84-119).
135
[57] Hodgkin, A. L. ; Huxley, A.F. [1952] : A quantitative description of
membrane current an its application to conductance and excitation
in nerve, J. Physiol. Lond.117 : 500-544
[58] Holden, A. V., Hyde, J., and Muhamad, M. A. [1991] : Equilibria,
Periodicity, Bursting and Chaos in Neural Activity, in Pasemann,
F., and Doebner, H. D. (eds), World Scientific Pub., Singapore :
97-187.
[59] Hopfield, J. J. [1982] : Neural networks and physical systems with
emergent collective computational abilities, Proc. Natl.Acad.Sci.
79 : 2554-58.
[60] Hopfield, J. J [1984] : Neurones with graded response have collective computational propierties like those of two states neurones.
Proccedings of National Academy of Science. 81 : 3088-92
[61] Hoppensteadt, F. [1986] : An Introduction to the Mathematics of
Neurons, Cambridge University Press.
[62] Hoppensteadt, F.C. [1993] : Analysis and Simulation of Chaotic
Systems, Springer Verlag, New York.
[63] Hoppensteadt, F., Izhikevich, E. M. [1997] : Weakly Connected
Neural Networks. Springer Verlag. New York, Inc.
[64] Ingber, L. ; Nuñez P. [1990] : Multiple scales of statistical physics
of Neocortex. Aplication to the Electroencephalography. Mathematical and Computer Modelling.
[65] Ingber, L. [1991] : Statistical mechanics of neurocortical interactions : A scaling paradigm applied to electroencephalography.,
Physical Review. A 44 : 4017-60.
[66] Ingber, L. [1994] : Statistical Mechanics of Multiple Scales of Neocortical Interactions.
[67] Ingber, L. [1994] : Statistical Mechanics of Neocortical Interactions. : Path-integral evolution of Short-term memory. Physical
Review E 49 : 4652-64.
136
[68] Ingber, L. [1995] : Statistical mechanics of multiple scales of neocortical interactions. In : Neocortical dynamics and human EEG.
rhythms, ed. P.L. Nuñez. Oxford University Press.
[69] Jahnsen, H. ; Llinas, R. [1984a] : Electrophysiological properties
of guinea-pig thalamic neurons : An in vitro studio. J. Physiol.
(London). 349 : 205-226.
[70] Jahnsen, H. ; Llinas, R. [1984b] : Ionic basis for the electroresponsiveness and oscillatory properties of guinea-pig thalamic neurons
in vitro. J. Physiol. (London) 349 : 227-2447.
[71] Kandel, E. R., Schwartz, J. H. ; Jessell, T. M. [1997] : Neurociencia
y conducta. Prentice Hall. Madrid.
[72] Katznelson, R. D. [1981] : Normal modes of the brain : Neuroanatomical basis and a physiological theoretical model. In : Electric fields of the brain : The Neurophysics of EEG. P.L. Nuñez.
Oxford University Press. New York.
[73] Kepler, T. B., Abbot, L. F., Marder, E. [1992] : Reduction of
conductance based neuron models. Biol. Cybern. 66 : 381-387.
[74] Kevorkian, J., Cole, J. D. [1996] : Multiple Scale and Singular
Perturbation Methods, Springer-Verlag, New York, Inc.
[75] Lopes da Silva, F. H. [1991] : Neural mechanisms underlying brain
waves : From neural membranes to networks electroencephagr.
Clin. Neurophysiol. 79 : 81-93.
[76] Lopes da Silva, F. H., van Rotterdam, A., Barts, P., van Heusden,
E., Burr, B. [1976] : Models of neurones populations the basic
mechanism of rhytmicity. Progress in Brain Research. 45 : 282308.
[77] Lopes da Silva, F. A. ; Hoeks, A., Smiths, H. and Zettergerg L. H.
[1974] : Model of Brain Rhytmic Activity, Kybernetik. 15 : 27-37.
[78] Lukatch, H. S.; Mac Iver, B. [1997] : Physiology, Pharmacology and Topography of cholinergic neocortical oscillations in vitro.
Journal of Neuroscience. 2427-2444.
137
[79] Mac. Gregor, R. J. [1987] : Neural and Brain Modeling. Academic
Press, Inc. San Diego. Cal.
[80] Mc Culloch, W. S. ; Pitts, W. H. [1943] : A logical calculus of ideas
immanent in nervous activity, Bull. Math. Biophys. 5 : 115-33.
[81] Menon, V. [1991] : Population oscillations in neuronal groups.
International Journal of Neural Systems. Vol. 2, no. 3, 237-62.
[82] Mischenko, E. F.; Rozov, N. Kh. [1980] : Differential Equations
with Small Parameters and Relaxation Oscillations, Plenum, New
York.
[83] Mountcastle, V. B. [1978] : An organizing principle for cerebral
function : the unit module and the distributed system. In : The
Mindful Brain : Cortical Organization and the Group Selective
theory of Higher Function, ed. G. M. Edelman and V. B. Mountcastle. MIT Press, Cambridge.
[84] Mountcastle, V. B. [1998] : Perceptual Neuroscience the Cerebral
Cortex. Harvard University Press. Cambridge, Massachusetts, and
London, England.
[85] Natterer, F. [1986] : The mathematics of Computerized Tomography. Teubner, Stuttgart.
[86] Nuñez, P. L. [1981] : A study of origins of the time dependencies
of scalp EEG : I. Theoretical basis, IEEE Trans. Bio-Med Eng. 28
: 271-280, 1981b.
[87] Nuñez, P. L. [1981] : A study of origins of the time dependencies of
scalp EEG: II. Experimental support of theory,IEEE Trans. BioMed Eng. 28 : 281-288, 1981c.
[88] Nuñez, P. L. [1988] : Global Contribution to cortical dynamics :
theoretical and experimental evidence for standing waves phenomena. In : Dynamics of Sensory and Cognitive Processing by the
Brain. (E. Baser, Ed.). Springer Verlag. New York.
138
[89] Nuñez, Paul L; Katznelson, R. D. [1981] : Electric Fields on the
Brain, the Neurophysics of EEG, Oxford University Press, New
York.
[90] Nuñez, P. [1989] : Generation of Human EEG by a combination of
long and short range neurocortical interactions, Brain topography.
1 : 199-215.
[91] Nuñez, P., Srinivasan, R. [1993] : Implications of recording
strategy for estimates of neurocortical dynamics with electroencephalography chaos. 3 : 257-66.
[92] Nuñez, P. L. [1995] : Neocortical Dynamics and Human EEG Rythms, Oxford University Press.
[93] Plonsey, R. [1969] : Bioelectric Phenomena, Mc Graw Hill, New
York.
[94] Petsche, H., Sterc, J. [1968] : The significance of the cortex for
the traveling phenomenon of brain waves. Electroencephalography
and Clinical Neurophysiology. 25 : 1222.
[95] Rall, W., Shepherd, G. M. [1968] : Theoretical reconstruction of
field potentials and dendrodendritic synaptic interactions in olfactory bulb. J. Neurophysiol. 31 : 884-915.
[96] Ramacher, U. [1993] : Hamiltonian Dynamics of Neural Networks.
Neural Networks. 6 : 547-557.
[97] Rinzel, J. ; Ermentrout, G. B. [1989] : Analysis of neural excitability and oscillations. In Methods of neural Modelling : From
Synapses to networks (C. Koch and I. Segev. Eds), Cambridge,
MA : MIT Press.
[98] Rosenblatt, F. [1957] : The perception a perceiving and recognizing
automation, , Cornell Aeronautical Laboratory Report.
[99] Rose, R. M. ; Hindmarsh, J. L. [1985] : A model of a thalamic
neuron, Proc. R. Soc. Lond.
139
[100] Rush, M. E., Rinzel, J. [1994] : Analysis of bursting in a thalamic
neuron model. Biol. Cybern. 71, 281-291.
[101] Sánchez, F., Maini, P. K. [1997] : Travelling wave phenomena
in non-linear diffusion degenerate Nagumo equations. Journal of
mathematical biology, 35 : 713-728.
[102] Sánchez, F., Maini, P. K. [1997] : Travelling wave phenomena in
some degenerate reaction diffusion equations. Journal of differential equations. Vol. 117, no. 2 : 281-319.
[103] Sarvas, J. [1987] : Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem. Phys. Med. Biol. Vol.32,
no.1 : 11-22.
[104] Silberstein, R. B. [1994] : Neuromodulation of neocortical dynamics. In : Neocortical dynamics and human EEG rhythms. ed. P. L.
Nuñez. Oxford University Press.
[105] Silberstein, R. B. [1994] : Steady-State Visually Evoked Potentials, Brain Resonances, and Cognitive Processes In : Neocortical
dynamics and human EEG rhythms. ed. P. L. Nuñez. Oxford University Press.
[106] Smith, C. U. M. [1982] : El cerebro. Sexta ed. Alianza Universidad
Madrid.
[107] Somogyi, P. ; Tamás, G. ; Lujan, R., Buhl, E. [1998] : Salient features of synaptic organisation in the cerebral cortex. Brain research
reviews, 26 : 113-135.
[108] Temam, R. [1998] : Infinite-Dimensional Dynamical Systems. 2 ed.
Springer Verlag, New York, Berlin, Heidelberg.
[109] Tijonov A.N., Goncharsky A.V. [1987] : Problemas mal planteados
de las ciencias naturales. (Editorial de la Universidad Estatal de
Moscú).
[110] Tikhonov, V. S. [1985] : Differential Equations. Springer Verlag.
Berlin , Heidelberg.
140
[111] Troy, W. C. [1978] : The bifurcation of periodic solutions in the
Hodgkin-Huxley equations. Quart. Appl. Math., Vol. 36, 73-83.
[112] Tyson, J. J. ; Keener, J. P. [1988] : Singular Perturbation Theory
of Traveling Waves in Excitable Medic (a review). Physica.D 32 :
327-361.
[113] van Rotterdam, A., Lópes da Silva, F. H., van den Ende, J.,
Viergeuer, M. A., Hermans, A. J. [1982] : A model of the spaticaltemporal characteristics of the alpha rhytm. Bull. Math. Biol. 4e :
283-305.
[114] Wang, X. J. ; Rinzel, J.[1993] : Spindle rhythmicity in the reticularis thalami nucleus : synchronization among mutually inhibitory
neurons. Neuroscience. 53 : 899-904.
[115] Wang, X. J. [1994] : Multiple dynamical modes of thalamic relay neurons : Rhythmic bursting and intermittent phase-locking.
Neuroscience. Vol. 59, no.1, 21-31.
[116] Weiss, T. F. [1996] : Celular biophysics, Vol. 1, Transport. MIT
Press, Cambridge, MA.
[117] Weiss, T. F. [1996] : Celular biophysics, Vol. 2, Electrical properties. MIT Press, Cambridge, MA.
[118] Widrow, B., Hoff, M. [1960] : Adaptative Switching Circuits
IREWESCON Convention Record, Part 4, 96-104.
[119] Wilson, H. R., Cowan, J. D. [1972] : Excitatory and inhibitory
interaction in localized populations of model neurons. Biophys.
Journal 12. 1-24.
[120] Wilson, H. R., Cowan, J. D. [1973] : A mathematical theory of
the functional dynamics of cortical and thalamic nervous tissue.
Kybernetik 13. 13-80.
[121] Wright, J.J. ; Liley, T. J. [1996] : Dynamics of the Brain at Global
and Microscopic scales : Neural Networks and the EEG, Behavioral and the brain sciences. 19 : 285-320.
141