Poblaciones excitatorias e inhibitorias: estabilidad y respuesta colectiva

Una región cortical puede mantener una actividad global casi constante mientras sus neuronas disparan de forma irregular. También puede oscilar sin que ninguna de sus unidades, aislada, posea el ritmo del conjunto. Para estudiar estos comportamientos conviene abandonar el voltaje de una célula particular y trabajar con variables que describan poblaciones. Ese paso abre el análisis de los circuitos de retroalimentación, aunque obliga a precisar qué actividad se promedia y a qué escala temporal.

El modelo de Wilson–Cowan describe la interacción entre una población excitatoria y otra inhibitoria mediante ecuaciones diferenciales no lineales acopladas (Wilson y Cowan, 1972, 1973). Su interés reside en que vincula arquitectura, respuesta y estabilidad con las herramientas del plano de fases. Cualquier interpretación depende de los supuestos de agregación y de la forma elegida para la función de respuesta. Por esa razón, un parámetro ajustado a una señal macroscópica no equivale de manera automática a una cantidad sináptica medida.

El desarrollo que sigue usa una variante de tasas con dos variables y parámetros escogidos para que el equilibrio pueda verificarse a mano. A partir de ese equilibrio se obtienen el jacobiano, la respuesta a perturbaciones pequeñas y la sensibilidad a entradas externas. El ejemplo muestra además que un aumento de la actividad inhibitoria observada no implica que la excitación haya disminuido. En un circuito cerrado, cada respuesta debe leerse a la luz de los recorridos de realimentación que la producen.

403.1. Qué representan las dos poblaciones

Sean $E(t)$ e $I(t)$ las actividades normalizadas de dos poblaciones. En la variante que usaremos toman valores entre cero y uno, de modo que no representan milivoltios ni números absolutos de neuronas. Pueden interpretarse como tasas escaladas o como una medida agregada de activación, siempre que la observación empleada para estimarlas sea coherente con esa elección. La letra $I$ , que en capítulos anteriores designaba una corriente, nombra aquí a la población inhibitoria; el contexto y las unidades bastan para distinguir ambos usos. Para evitar una ambigüedad adicional, la matriz identidad se escribirá $\mathbb{1}$ .

Las ecuaciones del modelo son

\tau_E\dot E=-E+F(w_{EE}E-w_{EI}I+P),

\tau_I\dot I=-I+F(w_{IE}E-w_{II}I+Q), \qquad F(u)=\frac{1}{1+e^{-u}}.

Los pesos $w$ son no negativos y los signos de la excitación y la inhibición aparecen de forma explícita. $P$ y $Q$ son entradas efectivas que pueden absorber desplazamientos de umbral. Las constantes $\tau_E$ y $\tau_I$ fijan la rapidez con que cada población relaja hacia su respuesta. Los argumentos de la sigmoide son adimensionales; si se parte de tasas físicas, los factores de escala correspondientes quedan incorporados en los pesos. Los supuestos de promedio temporal y poblacional en que descansan los modelos de tasas se tratan con detalle en los textos de neurociencia teórica (Dayan y Abbott, 2001).

Se trata de una variante de tasas dentro de la familia Wilson–Cowan. La formulación original describe la proporción de células activas en cada población, descuenta mediante un factor multiplicativo la fracción que permanece refractaria y llega a ecuaciones diferenciales tras un promedio temporal grueso (Wilson y Cowan, 1972; Chow y Karimipanah, 2020). El trabajo de 1973 extendió el planteamiento a un continuo espacial de tejido cortical (Wilson y Cowan, 1973; Destexhe y Sejnowski, 2009). La versión simplificada facilita el cálculo, a cambio de impedir que resultados o parámetros obtenidos con otra versión se trasladen sin revisión. El nombre del modelo designa una familia de ecuaciones; el análisis lo determinan las ecuaciones concretas que se escriben.

En el borde $E=0$ , la derivada de $E$ vale $F(\cdot)/\tau_E$ , cantidad estrictamente positiva porque la logística nunca alcanza el cero. En $E=1$ vale $[F(\cdot)-1]/\tau_E$ , estrictamente negativa por la razón simétrica. Lo mismo ocurre con $I$ en sus bordes. El campo vectorial apunta entonces hacia el interior en todo el contorno del cuadrado unidad, que resulta positivamente invariante para el flujo continuo, de manera que las tasas normalizadas permanecen acotadas. Un integrador explícito con un paso excesivo puede violar esta propiedad. Una tasa negativa en la salida numérica revelaría en ese caso una falla del cálculo y carecería de valor como predicción del modelo.

La sigmoide introduce una ganancia que depende del estado. Cerca de su centro, una variación de la entrada produce una variación relativamente grande de la salida, mientras que en las zonas saturadas el mismo cambio apenas modifica la actividad. Un circuito con pesos idénticos puede, por tanto, tener una estabilidad distinta según su punto de operación. Hablar de «más excitación» sin especificar el estado y la ganancia deja indeterminada la consecuencia dinámica.

403.2. Un equilibrio verificable

Elijamos $\tau_E=\tau_I=10$ ms y

w_{EE}=8,\quad w_{EI}=10,\quad w_{IE}=10,\quad w_{II}=2,\quad P=1,\quad Q=-4.

En $E^*=I^*=1/2$ , la entrada efectiva de la población excitatoria es $8/2-10/2+1=0$ y la de la inhibitoria es $10/2-2/2-4=0$ . Como $F(0)=1/2$ , ambas derivadas se anulan. El valor negativo de $Q$ representa un desplazamiento efectivo respecto del centro de la sigmoide. Nada en el modelo permite identificarlo con una corriente negativa medida en las neuronas inhibitorias.

Encontrar un equilibrio equivale a resolver a la vez las ecuaciones de ambas nulclinas. La nulclina excitatoria reúne los puntos donde $E=F(w_{EE}E-w_{EI}I+P)$ ; la inhibitoria se define de manera análoga. Solo en su intersección cesa la evolución de las dos variables. Sobre una nulclina aislada, la otra variable sigue cambiando y desplaza al sistema fuera de ella.

Que esta intersección exista no demuestra que sea única. Wilson y Cowan (1972) describieron en su modelo regímenes con histéresis y varios estados estables, además de actividad de ciclo límite, de modo que la multiplicidad de equilibrios es una posibilidad real dentro de la familia. Para buscar otras intersecciones se emplean integraciones desde condiciones iniciales diversas, continuación de ramas y argumentos de monotonía. Una solución hallada por un algoritmo local identifica solo un candidato. El espacio de estados acotado facilita una exploración sistemática, aunque una malla finita no sustituye una demostración global.

El análisis que sigue se restringe a la respuesta local alrededor del equilibrio explícito. Esa restricción permite calcular la estabilidad sin pretender describir todas las trayectorias posibles. Si una perturbación aleja mucho al sistema, las pendientes de la sigmoide cambian y la aproximación lineal pierde validez. El punto de operación forma parte del resultado.

403.3. Del jacobiano a una oscilación amortiguada

Denotemos por $g_E=F'(u_E^*)$ y $g_I=F'(u_I^*)$ las ganancias locales. Para la logística se cumple $F'(u)=F(u)[1-F(u)]$ , de modo que en el equilibrio ambas valen $1/4$ . El jacobiano es

J=\begin{pmatrix} (-1+g_Ew_{EE})/\tau_E&-g_Ew_{EI}/\tau_E\\ g_Iw_{IE}/\tau_I&(-1-g_Iw_{II})/\tau_I \end{pmatrix} =\begin{pmatrix} 0.10&-0.25\\ 0.25&-0.15 \end{pmatrix}\ \mathrm{ms}^{-1}.

La traza vale −0.05 ms $^{-1}$ y el determinante 0.0475 ms $^{-2}$ . Los valores propios son

\lambda_{\pm}=-0.025\pm0.216506\,i\quad\mathrm{ms}^{-1}.

La parte real negativa indica decaimiento local. La parte imaginaria introduce una rotación en el plano, por lo que una perturbación pequeña genera una oscilación amortiguada. La envolvente lineal decae con un tiempo característico de 40 ms y el período es $2\pi/0.216506\approx29.02$ ms, que corresponde a unos 34.46 ciclos por segundo en la aproximación local.

La frecuencia obtenida pertenece a los parámetros elegidos. Que caiga dentro del intervalo que la electrofisiología suele llamar gamma no identifica ningún mecanismo en un registro real. Conviene distinguir además una oscilación amortiguada tras un pulso de una oscilación sostenida autónoma. En la primera, la amplitud se extingue cuando cesa la perturbación; en la segunda, el sistema mantiene un ciclo estable sin necesitar una secuencia periódica de entradas.

La matriz contiene una contribución autoexcitatoria local positiva en la población $E$ . Si esa ecuación se aislara, el término diagonal 0.10 amplificaría cualquier desviación. El acoplamiento con $I$ aporta una retroalimentación que, con estos parámetros, estabiliza el conjunto. Examinar cada diagonal por separado haría perder el mecanismo de regulación recíproca, porque la estabilidad pertenece al sistema acoplado.

Para un sistema real de dos dimensiones, una traza negativa junto con un determinante positivo equivale a que ambos valores propios tengan parte real negativa. Un determinante negativo implica un punto de silla. Si el discriminante $\operatorname{tr}(J)^2-4\det(J)$ es negativo, los valores propios son complejos conjugados; en el ejemplo vale −0.1875 ms $^{-2}$ . Estas condiciones son locales y se aplican al jacobiano de la variante especificada. Una red con más variables o con retardos de transmisión requiere otros criterios.

403.4. Entrada externa y respuesta estacionaria

Introduzcamos un incremento pequeño y constante $\delta P$ . Cerca del equilibrio,

\frac{d}{dt}\begin{pmatrix}\delta E\\\delta I\end{pmatrix} =J\begin{pmatrix}\delta E\\\delta I\end{pmatrix} +\begin{pmatrix}0.025\\0\end{pmatrix}\delta P.

El factor 0.025 procede de $g_E/\tau_E$ . Si el sistema permanece estable, el desplazamiento estacionario lineal satisface

\begin{pmatrix}\delta E_\infty\\\delta I_\infty\end{pmatrix} =-J^{-1}\begin{pmatrix}0.025\\0\end{pmatrix}\delta P =\begin{pmatrix}0.078947\\0.131579\end{pmatrix}\delta P.

Una entrada dirigida a la población excitatoria eleva ambas actividades. En este ejemplo, la respuesta inhibitoria normalizada resulta incluso mayor que la excitatoria. Interpretar un incremento observado de $I$ como prueba de que la causa inicial fue una mayor entrada a la población inhibitoria sería un error, ya que la respuesta puede surgir del reclutamiento de la inhibición por la excitación. Atribuir el cambio exige considerar la estructura del circuito y las demás entradas compatibles con la observación.

Si $\delta P=0.01$ , la aproximación predice incrementos de aproximadamente 0.000789 y 0.001316. Al resolver el equilibrio de las ecuaciones no lineales con ese mismo incremento se obtienen valores que difieren de los lineales en menos de $10^{-7}$ . Para incrementos mayores, la ganancia cambia y los desplazamientos dejan de escalar en proporción exacta. Comparar varias amplitudes delimita el rango en que la respuesta lineal ofrece una descripción suficiente.

La respuesta transitoria contiene información adicional. Un pulso breve puede elevar primero $E$ y reclutar después $I$ , con un retorno al equilibrio que incluye sobrepaso. El desplazamiento estacionario no registra ese orden temporal. Dos parametrizaciones pueden compartir respuestas finales y diferir en sus latencias o en su amortiguamiento, de modo que una evaluación que mida solo el estado final no podrá distinguirlas.

Formalmente, la respuesta a una entrada que varía en el tiempo se obtiene mediante la exponencial matricial y una convolución. Para $\delta x(0)=0$ ,

\delta x(t)=\int_0^t e^{J(t-s)}b\,\delta P(s)\,ds, \qquad b=(0.025,0)^{\mathsf T}.

La expresión indica que la entrada se transforma según todos los modos del sistema. Si solo se observa una combinación $y=c^{\mathsf T}\delta x$ , algunos modos pueden aparecer con muy poco peso. La ausencia de una oscilación evidente en la salida no demuestra, por tanto, que el modelo interno carezca de modos complejos.

403.5. Cambiar excitación manteniendo el punto de operación

Para estudiar la influencia de $w_{EE}$ sin desplazar el equilibrio, definamos $P=5-w_{EE}/2$ y conservemos los demás parámetros. Con esa elección, $E^*=I^*=1/2$ sigue siendo solución y las ganancias locales permanecen en $1/4$ . La manipulación es un diseño matemático controlado. En un tejido real, el cambio de una conductancia rara vez llega acompañado del reajuste exacto de otra entrada, de modo que el experimento formal debe describirse como lo que es.

La traza pasa a ser $0.025w_{EE}-0.25$ y el determinante $0.0775-0.00375w_{EE}$ . Para $w_{EE}=10$ la traza se anula y el determinante vale 0.04, con valores propios puramente imaginarios $\pm0.2i$ ms $^{-1}$ . La traza cruza el cero con derivada 0.025 respecto del parámetro, distinta de cero. Con $w_{EE}=12$ los valores propios son $0.025\pm0.1785i$ ms $^{-1}$ y el equilibrio se ha convertido en un foco inestable.

Estas son las condiciones lineales de una bifurcación de Hopf. Para caracterizarla falta verificar una condición no lineal de no degeneración, que suele expresarse mediante el primer coeficiente de Lyapunov. Su signo decide si el ciclo que emerge es estable o inestable y a qué lado del valor crítico aparece (Kuznetsov, 2004). Una serie temporal oscilatoria obtenida en un barrido finito no sustituye esa caracterización, porque puede corresponder a un transitorio largo o a un artefacto de la integración. El análisis de bifurcaciones del modelo de dos poblaciones en función de dos parámetros divide el plano paramétrico en regiones de comportamiento cualitativamente equivalente. En algunas de ellas el sistema llega a un estado estacionario o a un ciclo límite según la condición inicial (Borisyuk y Kirillov, 1992; Cowan et al., 2016).

La afirmación «aumentar la excitación produce oscilaciones» resulta demasiado amplia. Su consecuencia depende del desplazamiento del equilibrio (que arrastra consigo las ganancias) y de la intensidad de la realimentación inhibitoria. Si la sigmoide entra en saturación, su pendiente disminuye y un peso mayor deja de implicar una ganancia efectiva mayor. El circuito debe analizarse en el estado al que efectivamente llega. En el trabajo original, la frecuencia del ciclo límite variaba de forma monótona con la intensidad del estímulo (Wilson y Cowan, 1972), lo que ilustra que el ritmo de una oscilación autónoma tampoco es una constante del circuito.

Modificar $\tau_I$ introduce otra diferencia. La velocidad de reclutamiento inhibitorio altera el desfase entre poblaciones y el amortiguamiento. Dos circuitos con los mismos pesos y el mismo equilibrio pueden responder de forma distinta por sus constantes temporales. En modelos con retrasos de transmisión, el estado relevante incluye la historia reciente y el análisis de estabilidad deja de reducirse a los dos valores propios de una matriz ordinaria.

403.6. Balance y significado de la observación

La expresión «balance entre excitación e inhibición» puede referirse a corrientes, tasas, pesos o efectos netos. Sohal y Rubenstein (2019) subrayan que tanto la excitación como la inhibición son entidades multidimensionales, compuestas por subtipos celulares con dinámicas en varias escalas temporales, por lo que un cociente único simplifica en exceso lo que ocurre en un circuito. En nuestro equilibrio, $E=I$ porque los parámetros se eligieron para producir actividades normalizadas iguales. Esa igualdad no establece equivalencia de corrientes ni de número de sinapsis. Los pesos son distintos y las variables carecen de las unidades que harían falta para sostener una equivalencia fisiológica.

Un cociente $E/I$ tampoco resume la estabilidad. En la familia de la sección 403.5 el cociente se mantiene en uno mientras la traza cruza el cero. Dos sistemas con el mismo cociente pueden recuperarse con rapidez, oscilar con un amortiguamiento lento o ser inestables. La dinámica depende de la relación entre variaciones, que el jacobiano expresa localmente, además de los niveles medios.

La señal adquirida necesita su propia ecuación de observación. Con $y=E+I$ , una trayectoria produce un registro; con $y=E-I$ , la misma trayectoria puede mostrar cancelaciones distintas. Ninguna de estas fórmulas constituye por sí sola un modelo de EEG o de señal BOLD. Su función aquí es mostrar que una salida agregada depende del modo en que las variables internas se convierten en medición.

Una medición electrofisiológica macroscópica depende de corrientes y de la geometría del tejido, mientras que una señal hemodinámica introduce otras transformaciones y otras escalas temporales. Las plataformas de simulación de redes cerebrales completas, como The Virtual Brain, combinan modelos de dinámica local con la generación de señales de fMRI, EEG o MEG con el propósito de inferir mecanismos a partir de esas señales (Sanz-Leon et al., 2013). Breakspear (2017) revisa los supuestos de este enfoque, la evidencia que lo respalda y los desafíos que siguen abiertos. Ajustar una oscilación de tasas directamente a una señal filtrada puede hacer que los parámetros neuronales absorban características del proceso de observación. Mantener separadas ambas capas ayuda a reconocer esa ambigüedad. El capítulo 405 desarrollará esta arquitectura de forma explícita.

La estimación del cociente entre excitación e inhibición a partir de potenciales de campo ilustra la misma dependencia. Gao et al. (2017) derivaron de un modelo computacional una relación entre ese cociente y la pendiente del espectro de potencia, que contrastaron con registros de ratas y macacos, incluida la administración de un anestésico general en macacos. El valor así obtenido hereda los supuestos del modelo que conecta corrientes sinápticas con espectro. Su uso en otra población o con otra técnica de registro requiere comprobar que esos supuestos se mantienen antes de leerlo como medida directa.

En psiquiatría, la alteración del balance entre excitación e inhibición se ha propuesto como marco organizador para investigar mecanismos de trastornos del neurodesarrollo como el autismo, con evidencia proveniente en gran medida de modelos murinos (Sohal y Rubenstein, 2019). En ese contexto, una diferencia de amplitud o de frecuencia puede ser compatible con cambios en varios pesos, entradas o constantes temporales. El modelo sirve para construir hipótesis y diseñar perturbaciones que las pongan a prueba. Atribuir un mecanismo particular requiere evidencia capaz de discriminar entre alternativas. Un parámetro denominado inhibición sigue siendo una cantidad inferida dentro de una representación, salvo que exista una calibración independiente que conecte su escala con una medida biológica.

403.7. Simulación, controles y perturbaciones

El código siguiente verifica el equilibrio y los valores propios. Después compara la evolución no lineal tras una perturbación pequeña del estado con la predicción lineal. Usa NumPy y SciPy; el paso de salida fijado por t_eval no coincide con el paso interno del integrador adaptativo.

import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import expm
from scipy.special import expit

def rhs(t, x):
    E, I = x
    return np.array([
        (-E + expit(8*E - 10*I + 1))/10,
        (-I + expit(10*E - 2*I - 4))/10
    ])

xeq = np.array([0.5, 0.5])
J = np.array([[0.10, -0.25], [0.25, -0.15]])
assert np.allclose(rhs(0, xeq), 0)
print('Valores propios:', np.linalg.eigvals(J))

ts = np.linspace(0, 200, 1001)  # milisegundos
delta = np.array([0.001, 0.0])
sol = solve_ivp(rhs, (0, 200), xeq+delta, t_eval=ts,
                rtol=1e-9, atol=1e-11)
assert sol.success
lineal = np.array([xeq + expm(J*t) @ delta for t in ts])
print('Error máximo lineal:', np.max(np.abs(sol.y.T-lineal)))
print('Ganancia estacionaria:', -np.linalg.solve(J, [0.025, 0]))

Con estos valores, el error máximo entre la trayectoria no lineal y la lineal es de unos $1.4\times10^{-8}$ y la ganancia estacionaria reproduce los coeficientes 0.078947 y 0.131579 de la sección 403.4.

La comprobación de tolerancias debe acompañar el uso del resultado. Puede repetirse el cálculo con tolerancias más estrictas y verificar que los observables cambien mucho menos que la diferencia investigada. Reducir la amplitud inicial permite examinar si el error entre el sistema no lineal y su linealización disminuye como se espera cerca del equilibrio. En este ejemplo el error cae aproximadamente con el cubo de la amplitud, desde unos $1.4\times10^{-5}$ para una desviación inicial de 0.01 hasta unos $1.4\times10^{-8}$ para 0.001. El exponente tres se explica porque la segunda derivada de la logística se anula en su centro y ambas poblaciones operan justamente allí, de modo que desaparecen los términos cuadráticos del desarrollo de Taylor. Con amplitudes todavía menores, el error deja de bajar a ese ritmo porque alcanza el piso fijado por las tolerancias del integrador. Ese control evalúa la aproximación local, además de la precisión numérica.

Una segunda comprobación calcula el jacobiano por diferencias finitas y lo compara con la derivada analítica, lo que permite detectar signos o factores temporales equivocados. El incremento debe ser lo bastante pequeño para aproximar la derivada y lo bastante grande para evitar la cancelación numérica. Con diferencias centradas, el error máximo es de unos $2\times10^{-4}$ con un incremento de $10^{-2}$ , baja a unos $10^{-11}$ con $10^{-6}$ y vuelve a crecer hasta unos $4\times10^{-9}$ con $10^{-9}$ . Comparar varios incrementos permite identificar así una zona estable de estimación. El bloque siguiente, que se ejecuta a continuación del anterior, reúne estos controles y calcula también la amplificación transitoria que se discute en la sección 403.12.

def jac_fd(f, x, h):
    cols = [(f(0, x + h*e) - f(0, x - h*e))/(2*h) for e in np.eye(len(x))]
    return np.column_stack(cols)

for h in (1e-2, 1e-4, 1e-6, 1e-9):
    print('h =', h, 'error:', np.max(np.abs(jac_fd(rhs, xeq, h) - J)))

for a in (1e-2, 1e-3):
    d = np.array([a, 0.0])
    s = solve_ivp(rhs, (0, 200), xeq + d, t_eval=ts,
                  rtol=1e-12, atol=1e-14)
    lin = np.array([xeq + expm(J*t) @ d for t in ts])
    print('amplitud =', a, 'error:', np.max(np.abs(s.y.T - lin)))

tt = np.linspace(0, 200, 20001)
sv = np.array([np.linalg.norm(expm(J*t), 2) for t in tt])
k = int(np.argmax(sv))
print('Máxima amplificación:', sv[k], 'en t =', tt[k], 'ms')

Para estudiar la respuesta a entradas se modifica $P(t)$ en lugar del estado inicial. Un pulso rectangular debe especificar amplitud y duración. Un impulso ideal pertenece a otra formulación matemática, de modo que un pulso finito solo lo aproxima si conserva su área. Si se comparan regiones con distintas constantes temporales, mantener la amplitud o mantener el área conduce a experimentos diferentes.

403.8. Una perturbación inhibitoria con respuesta paradójica

El mismo ejemplo permite calcular una respuesta menos intuitiva. Aumentemos ligeramente $Q$ , la entrada directa a la población inhibitoria. La aproximación estacionaria es

\begin{pmatrix}\delta E_\infty\\\delta I_\infty\end{pmatrix} =-J^{-1}\begin{pmatrix}0\\0.025\end{pmatrix}\delta Q =\begin{pmatrix}-0.131579\\-0.052632\end{pmatrix}\delta Q.

Ambas actividades disminuyen en el nuevo equilibrio lineal, de modo que la población que recibió una entrada positiva termina con menor actividad. No se trata de un error de signo. El aumento inicial de $I$ reduce $E$ . Esa reducción retira suficiente excitación de la población inhibitoria como para que su actividad final caiga por debajo del valor de partida. El cálculo con las ecuaciones no lineales y $\delta Q=0.01$ reproduce ambos descensos con la misma precisión que en la sección 403.4.

La respuesta inicial y la final tienen signos opuestos. Justo después de incrementar $Q$ , con los estados todavía en sus valores originales, $\dot I$ es positiva. La retroalimentación modifica luego el resto del sistema. Un registro que observe solo los primeros milisegundos y otro que promedie el estado tardío pueden describir efectos contrarios de la misma perturbación. La ventana temporal forma parte de la interpretación del experimento.

El signo puede derivarse sin invertir numéricamente toda la matriz. Las variaciones estacionarias satisfacen

(1-g_Ew_{EE})\delta E+g_Ew_{EI}\delta I=0,

-g_Iw_{IE}\delta E+(1+g_Iw_{II})\delta I=g_I\delta Q.

El determinante del sistema es

D=(1-g_Ew_{EE})(1+g_Iw_{II})+g_Eg_Iw_{EI}w_{IE}.

En el ejemplo vale 4.75 y la sensibilidad inhibitoria es $(1-g_Ew_{EE})g_I/D$ . Como $D=\tau_E\tau_I\det(J)$ , el determinante es positivo siempre que el equilibrio sea estable. El signo de la sensibilidad depende entonces solo de $1-g_Ew_{EE}$ . Dentro de esta linealización, la respuesta paradójica aparece exactamente cuando la recurrencia excitatoria, por sí sola, superaría a la fuga ( $g_Ew_{EE}>1$ ). Esa situación recibe el nombre de red estabilizada por inhibición. En el ejemplo, $g_Ew_{EE}=2$ , de modo que la población inhibitoria sostiene la estabilidad de un régimen que sin ella sería inestable.

Tsodyks et al. (1997) describieron este efecto en un modelo de células principales con conexiones recurrentes acopladas a interneuronas. Mostraron que aumentar la entrada inhibitoria externa dirigida a las interneuronas puede elevar su tasa de disparo, con la condición principal de que la red excitatoria sea inestable en ausencia de retroalimentación inhibitoria. El caso tratado aquí es la imagen especular de ese resultado. Ozeki et al. (2009) registraron intracelularmente en la corteza visual primaria del gato y encontraron que un estímulo en la periferia del campo receptivo reducía la inhibición que recibían las células, contra la expectativa de un aumento de la inhibición lateral. Interpretaron la supresión como un retiro de excitación propio de una red estabilizada por inhibición. Mediante estimulación optogenética de interneuronas positivas para parvalbúmina, Sanzeni et al. (2020) observaron supresión paradójica de la actividad inhibitoria en las cortezas visual, somatosensorial y motora del ratón. En promedio, el efecto apareció solo con expresión transgénica de la opsina, mientras que la expresión viral no lo produjo. Sadeh y Clopath (2021) revisan este conjunto de pruebas y sus consecuencias computacionales.

La inferencia inversa exige más cautela. Mahrach et al. (2020) analizaron respuestas a la fotoestimulación de interneuronas positivas para parvalbúmina en la corteza motora y en la corteza de barriles del ratón. Encontraron que los modelos de dos poblaciones no explicaban sus datos y que, en modelos de cuatro poblaciones con varias clases de interneuronas, la respuesta paradójica indicaba solo que la excitación recurrente no era demasiado intensa, sin constituir evidencia de estabilización por inhibición. La equivalencia entre efecto paradójico y red estabilizada por inhibición vale para el modelo de dos poblaciones. Extenderla a arquitecturas con más tipos celulares requiere un análisis propio.

El cálculo ofrece, en suma, una hipótesis experimental específica que no basta para identificar un mecanismo en una señal cerebral. Hay que conocer qué población se perturba, qué componente se observa y si la entrada actúa de forma selectiva. Una estimulación que alcance simultáneamente otros elementos puede producir otra respuesta. La magnitud y la duración de la perturbación también deben mantenerse dentro del rango en que la linealización es adecuada.

El ejemplo obliga a separar cantidades que suelen confundirse. La entrada que reciben las neuronas inhibitorias, la actividad que esas neuronas muestran y el efecto inhibitorio que ejercen sobre otra población pueden cambiar en direcciones diferentes durante un mismo ensayo. Tomar una de ellas como representación completa de las demás eliminaría precisamente la dinámica que el modelo explica.

403.9. La respuesta a una entrada periódica

Una perturbación sinusoidal permite estudiar cómo el circuito selecciona frecuencias. En la linealización, la transferencia desde $\delta P$ hacia las dos poblaciones se obtiene de $(s\mathbb{1}-J)^{-1}b$ . En nuestro caso,

H_E(s)=\frac{0.025(s+0.15)}{s^2+0.05s+0.0475}, \qquad H_I(s)=\frac{0.00625}{s^2+0.05s+0.0475}.

Al evaluar en $s=i\omega$ se obtienen la amplitud y la fase para cada frecuencia angular. En $\omega=0$ , las ganancias coinciden con las sensibilidades estacionarias calculadas antes, 0.078947 y 0.131579. Esta coincidencia proporciona una verificación interna entre dos derivaciones del mismo modelo.

El máximo de amplitud de una función de transferencia no tiene por qué coincidir con la parte imaginaria del valor propio, porque el numerador y el amortiguamiento influyen sobre el pico. La población observada también importa, ya que $H_E$ y $H_I$ tienen numeradores diferentes. En el ejemplo, $|H_E(i\omega)|$ alcanza su máximo cerca de $\omega\approx0.2170$ ms $^{-1}$ (unos 34.54 Hz) con una ganancia de 0.607, mientras que $|H_I(i\omega)|$ lo alcanza cerca de 0.2151 ms $^{-1}$ (unos 34.23 Hz) con una ganancia de 0.577. Ambos picos difieren de la parte imaginaria del valor propio, 0.2165 ms $^{-1}$ , además de diferir entre sí. Una frecuencia dominante en una salida puede reflejar la combinación de modos, entrada y observación, además de la dinámica interna.

Una entrada con espectro amplio puede excitar varias frecuencias a la vez. Si contiene más potencia cerca de un modo, la salida puede mostrar un pico sin que exista una oscilación autónoma. Para distinguir una respuesta forzada de un ciclo sostenido conviene examinar qué ocurre al retirar la entrada y cómo cambia la amplitud con su intensidad. El análisis debe considerar también el ruido capaz de mantener fluctuaciones alrededor del equilibrio.

En una simulación, el barrido puede hacerse aplicando perturbaciones pequeñas de distintas frecuencias y descartando el transitorio de cada ensayo. Después se estiman la amplitud y el desfase del régimen periódico. Si la perturbación es demasiado grande aparecen armónicos y la transferencia lineal deja de describir toda la salida. Ese desacuerdo puede señalar una no linealidad real del modelo, sin que haya un error del integrador.

403.10. Nulclinas y geometría alrededor del punto de trabajo

La sigmoide logística puede invertirse mediante $\operatorname{logit}(E)=\ln[E/(1-E)]$ . La nulclina excitatoria se escribe, para $0<E<1$ ,

I=\frac{w_{EE}E+P-\operatorname{logit}(E)}{w_{EI}}.

Su pendiente en el equilibrio es $(w_{EE}-1/[E(1-E)])/w_{EI}$ , que vale 0.4. La nulclina inhibitoria puede expresarse como $E=[\operatorname{logit}(I)+w_{II}I-Q]/w_{IE}$ ; su pendiente $dE/dI$ vale 0.6 y, localmente, $dI/dE$ vale $5/3\approx1.667$ .

Estas curvas ayudan a interpretar los signos del flujo. A un lado de la nulclina excitatoria $E$ aumenta y al otro disminuye. La nulclina inhibitoria organiza de manera semejante la segunda componente. Seguir ambas regiones de signos permite comprender el giro de una trayectoria, aunque la clasificación cuantitativa de la estabilidad requiere también las constantes temporales y las derivadas.

Una representación del plano de fases debe mostrar escalas comparables y distinguir las nulclinas de las trayectorias. Las flechas de un campo normalizado muestran la dirección del flujo y ocultan su velocidad. Con esa normalización, una zona lenta puede parecer tan rápida como otra. Para estudiar tiempos de recuperación conviene complementar la geometría con la serie temporal o con mapas de la magnitud de la velocidad.

Los bordes de la sigmoide añaden otra precaución. La transformación logit diverge al aproximarse a cero o a uno. Una gráfica que recorte el rango puede ocultar parte de la curva o sugerir una intersección inexistente. Los cálculos de equilibrio deben realizarse sobre las ecuaciones originales, con comprobación de residuos, especialmente cerca de la saturación.

403.11. Fluctuaciones alrededor del equilibrio

Una extensión lineal estocástica escribe $d\delta x=J\delta x\,dt+B\,dW_t$ , donde $W_t$ es un proceso de Wiener y $B$ fija cómo entra el ruido en cada población. Si $J$ es estable, la covarianza estacionaria $\Sigma$ satisface la ecuación de Lyapunov

J\Sigma+\Sigma J^{\mathsf T}+BB^{\mathsf T}=0.

La covarianza depende de la dinámica y del modo en que entra el ruido. Dos matrices $B$ distintas pueden producir asociaciones diferentes con el mismo jacobiano. A la inversa, cambios de $J$ y del ruido pueden compensarse y dejar covarianzas semejantes. Una matriz de correlación, por sí sola, no determina todos esos componentes.

Si el ruido afecta directamente a ambas poblaciones de manera correlacionada, puede aumentar su asociación sin que cambie ningún peso. Si entra solo en una, la otra recibe sus efectos a través del circuito. Comparar la covarianza con la respuesta a una entrada conocida aporta restricciones de distinta naturaleza. La primera describe fluctuaciones bajo la excitación espontánea modelada; la segunda examina cómo se transforma un estímulo controlado.

Un experimento de recuperación puede simular varios mecanismos de ruido y ajustar después un modelo que supone solo uno. Si los pesos estimados cambian para absorber esa diferencia, queda demostrado un sesgo posible del procedimiento. En datos reales, esa observación justifica analizar la sensibilidad a los supuestos sobre el ruido antes de interpretar cambios paramétricos como cambios sinápticos.

Las secciones anteriores no pretenden convertir dos ecuaciones en una descripción completa de una región. Muestran que un modelo pequeño genera predicciones temporales, frecuenciales y estadísticas susceptibles de contrastarse por separado. Cuantas más de esas predicciones se evalúen con información pertinente, mejor delimitado quedará el conjunto de mecanismos compatibles con la observación.

403.12. Estabilidad asintótica y amplificación transitoria

Que todos los valores propios tengan parte real negativa garantiza el retorno de perturbaciones suficientemente pequeñas hacia el equilibrio. No implica que cualquier medida de su amplitud disminuya desde el instante inicial. Nuestro jacobiano proporciona un ejemplo calculable. Escribamos $v=(\delta E,\delta I)^{\mathsf T}$ y consideremos la longitud euclidiana de esa desviación en las actividades normalizadas. Para la dinámica lineal,

\frac{d}{dt}\|v\|^2=v^{\mathsf T}(J+J^{\mathsf T})v =0.2(\delta E)^2-0.3(\delta I)^2.

Los términos cruzados se cancelan. Una desviación inicial exclusivamente excitatoria aumenta su longitud durante el primer instante, aunque las partes reales de ambos valores propios sean negativas. La trayectoria cambia después de dirección por el acoplamiento con la inhibición y termina decayendo. El circuito puede amplificar inicialmente una perturbación sin dejar de ser estable.

La solución exacta muestra cómo conviven ambos resultados. El jacobiano puede escribirse como $J=-0.025\,\mathbb{1}+K$ , con

K=\begin{pmatrix}0.125&-0.25\\0.25&-0.125\end{pmatrix}, \qquad K^2=-\omega^2\mathbb{1},\qquad \omega=\sqrt{0.046875}\approx0.216506.

Por tanto,

e^{Jt}=e^{-0.025t}\left[\mathbb{1}\cos(\omega t)+\frac{K}{\omega}\sin(\omega t)\right].

El factor exponencial controla el decaimiento de largo plazo, mientras el término matricial reorganiza la dirección y la amplitud durante la oscilación. Para una desviación inicial en $E$ , la componente inhibitoria aparece con un término proporcional al seno. La respuesta completa no puede reconstruirse usando únicamente el tiempo de decaimiento de cuarenta milisegundos.

Esta amplificación transitoria es posible porque el jacobiano no es normal, es decir, $JJ^{\mathsf T}\neq J^{\mathsf T}J$ y sus vectores propios no son ortogonales. Murphy y Miller (2009) mostraron que, cuando una retroalimentación inhibitoria estabiliza una excitación recurrente intensa, esa falta de normalidad equivale a una conectividad oculta de tipo anterógrado entre patrones de actividad. El resultado es una amplificación selectiva de algunos patrones sin prolongación de su duración, mecanismo que llamaron amplificación balanceada. Hennequin et al. (2014) extendieron la idea a redes con excitación recurrente intensa y aleatoria, estabilizada por una inhibición ajustada con criterios de control óptimo, capaces de amplificar de forma transitoria estados de actividad específicos.

Una forma de cuantificar la amplificación consiste en calcular, para cada tiempo, el mayor valor singular de $e^{Jt}$ . Su máximo describe la mayor amplificación posible de la longitud entre todas las direcciones iniciales de tamaño unitario, dentro de la aproximación lineal. En el ejemplo, ese máximo vale aproximadamente 1.46 y se alcanza hacia los 6.2 ms. La pregunta difiere de medir la respuesta a un pulso concreto, porque el pulso puede no excitar la dirección que alcanza ese máximo.

La medida de longitud también exige una elección. Aquí ambas actividades son adimensionales y comparten escala, por lo que la norma euclidiana resulta transparente. Si una variable expresara disparos por segundo y otra una concentración, sumar sus cuadrados sin normalización sería arbitrario. Una matriz de ponderaciones puede definir una magnitud relevante para la observación o para el costo de la perturbación. Cambiarla altera la amplificación medida aunque no modifique los valores propios del sistema.

En una simulación no lineal se puede aumentar gradualmente la amplitud de la perturbación y comparar la trayectoria con esta solución. Para desviaciones muy pequeñas, la discrepancia debe ser reducida. Al crecer la desviación cambian las pendientes de las sigmoides y la aproximación local pierde precisión. El umbral en que ocurre esa pérdida depende de la dirección y del observable, de modo que conviene informarlo mediante errores de trayectoria concretos.

Existe además una restricción global sencilla en nuestra variante. Como se vio en la sección 403.1, el campo vectorial apunta hacia el interior en todos los bordes del cuadrado de actividades admisibles, que es por ello invariante. Las trayectorias que comienzan allí permanecen dentro, incluso cuando una perturbación deja de ser pequeña.

La invariancia impide una divergencia ilimitada de la actividad, aunque no garantiza un único equilibrio ni excluye oscilaciones persistentes para ciertos parámetros. La estabilidad local, la amplificación transitoria y la organización global del plano de fases responden a preguntas diferentes. Estudiarlas por separado permite interpretar un aumento temporal de la señal sin confundirlo con inestabilidad y reconocer cuándo una predicción basada en el equilibrio necesita ampliarse.

El alcance de la aproximación puede comprobarse con dos perturbaciones de igual tamaño, una en cada población. Una desviación unitaria puramente excitatoria alcanza una norma máxima de 1.27 hacia los 4.0 ms. Una desviación unitaria puramente inhibitoria comienza disminuyendo, como indica el término $-0.3(\delta I)^2$ , aunque la rotación posterior la lleva a superar su tamaño inicial hasta 1.12 hacia los 8.9 ms. Sus máximos y latencias difieren pese a pertenecer al mismo circuito. Si se midiera únicamente la mayor amplitud de respuesta, esa diferencia podría atribuirse a un cambio de parámetros que nunca ocurrió. Conservar el lugar y la forma de la estimulación hace pertinente la comparación. Repetir después ambas entradas con otro parámetro sí ofrece un contraste entre circuitos, porque mantiene controlada la dirección inicial. La simulación separa así los efectos del sistema de los efectos del protocolo que lo perturba.

403.13. Una inferencia que el modelo permite y otra que deja abierta

Supongamos que un registro muestra una recuperación más lenta después de una perturbación. En el modelo, eso podría corresponder a valores propios con parte real más cercana a cero, interpretación que sugiere examinar la amplificación y la sensibilidad. Sin embargo, el mismo registro podría resultar de un filtro de observación más lento, de una entrada persistente o de una combinación de estados ocultos. El resultado observable debe confrontarse con cada una de esas explicaciones.

Una estrategia experimental compara respuestas a perturbaciones dirigidas a componentes distintos. Si solo se estimula $E$ y se observa una mezcla de ambas poblaciones, pueden quedar combinaciones de parámetros indistinguibles. Estimular otra entrada o registrar una segunda salida puede mejorar la identificación. La utilidad de la medición adicional depende de que genere predicciones diferentes entre los modelos compatibles con la primera.

En una simulación es posible conocer el mecanismo verdadero y construir una prueba de recuperación. Se generan datos con una parametrización, se añade ruido y se ajustan varias alternativas. Si el procedimiento confunde sistemáticamente los pesos, no conviene interpretar con precisión esos mismos parámetros en un estudio empírico más incierto. Esta prueba debe realizarse antes de otorgar significado clínico a una estimación.

También puede compararse un modelo de dos poblaciones con otro que incorpora una variable lenta adicional o más tipos celulares. El estudio de Mahrach et al. (2020) citado en la sección 403.8 es un ejemplo, ya que sus datos optogenéticos exigieron pasar de dos a cuatro poblaciones. Un mejor ajuste del modelo ampliado obliga a considerar su flexibilidad y su capacidad predictiva fuera de la muestra. Si la variable añadida permite anticipar respuestas a protocolos nuevos, aporta evidencia útil sobre la necesidad de esa estructura. Si solo mejora la reproducción del registro usado para ajustarla, la ganancia explicativa queda menos clara.

El recorrido del capítulo deja varios resultados delimitados. Se verificó un equilibrio, se calculó un foco estable y se estimó su respuesta a entradas pequeñas, constantes o periódicas. Se construyó una familia en la que el equilibrio permanece fijo mientras cambia la estabilidad local. También se obtuvo la condición exacta de la respuesta paradójica dentro de la linealización. Estas operaciones permiten formular preguntas precisas sobre circuitos colectivos. Atribuir una señal cerebral a una alteración específica de la excitación o de la inhibición requiere, además, resolver el problema de la observación y contrastar mecanismos rivales con perturbaciones capaces de distinguirlos.

Referencias

Borisyuk, R. M. y Kirillov, A. B. (1992). Bifurcation analysis of a neural network model. Biological Cybernetics, 66(4), 319–325. doi:10.1007/BF00203668
Breakspear, M. (2017). Dynamic models of large-scale brain activity. Nature Neuroscience, 20(3), 340–352. doi:10.1038/nn.4497
Chow, C. C. y Karimipanah, Y. (2020). Before and beyond the Wilson–Cowan equations. Journal of Neurophysiology, 123(5), 1645–1656. doi:10.1152/jn.00404.2019
Cowan, J. D., Neuman, J. y van Drongelen, W. (2016). Wilson–Cowan equations for neocortical dynamics. The Journal of Mathematical Neuroscience, 6(1). doi:10.1186/s13408-015-0034-5
Dayan, P. y Abbott, L. F. (2001). Theoretical neuroscience: Computational and mathematical modeling of neural systems. MIT Press.
Destexhe, A. y Sejnowski, T. J. (2009). The Wilson–Cowan model, 36 years later. Biological Cybernetics, 101(1), 1–2. doi:10.1007/s00422-009-0328-3
Gao, R., Peterson, E. J. y Voytek, B. (2017). Inferring synaptic excitation/inhibition balance from field potentials. NeuroImage, 158, 70–78. doi:10.1016/j.neuroimage.2017.06.078
Hennequin, G., Vogels, T. P. y Gerstner, W. (2014). Optimal control of transient dynamics in balanced networks supports generation of complex movements. Neuron, 82(6), 1394–1406. doi:10.1016/j.neuron.2014.04.045
Kuznetsov, Y. A. (2004). Elements of applied bifurcation theory. Springer. doi:10.1007/978-1-4757-3978-7
Mahrach, A., Chen, G., Li, N., van Vreeswijk, C. y Hansel, D. (2020). Mechanisms underlying the response of mouse cortical networks to optogenetic manipulation. eLife, 9, e49967. doi:10.7554/eLife.49967
Murphy, B. K. y Miller, K. D. (2009). Balanced amplification: A new mechanism of selective amplification of neural activity patterns. Neuron, 61(4), 635–648. doi:10.1016/j.neuron.2009.02.005
Ozeki, H., Finn, I. M., Schaffer, E., et al. (2009). Inhibitory stabilization of the cortical network underlies visual surround suppression. Neuron, 62(4), 578–592. doi:10.1016/j.neuron.2009.03.028
Sadeh, S. y Clopath, C. (2021). Inhibitory stabilization and cortical computation. Nature Reviews Neuroscience, 22(1), 21–37. doi:10.1038/s41583-020-00390-z
Sanz-Leon, P., et al. (2013). The Virtual Brain: A simulator of primate brain network dynamics. Frontiers in Neuroinformatics, 7. doi:10.3389/fninf.2013.00010
Sanzeni, A., Akitake, B., Goldbach, H. C., Leedy, C. E., Brunel, N. y Histed, M. (2020). Inhibition stabilization is a widespread property of cortical networks. eLife, 9, e54875. doi:10.7554/eLife.54875
Sohal, V. S. y Rubenstein, J. L. R. (2019). Excitation-inhibition balance as a framework for investigating mechanisms in neuropsychiatric disorders. Molecular Psychiatry, 24(9), 1248–1257. doi:10.1038/s41380-019-0426-0
Tsodyks, M. V., Skaggs, W. E., Sejnowski, T. J. y McNaughton, B. L. (1997). Paradoxical effects of external modulation of inhibitory interneurons. The Journal of Neuroscience, 17(11), 4382–4388. doi:10.1523/JNEUROSCI.17-11-04382.1997
Wilson, H. R. y Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal, 12(1), 1–24. doi:10.1016/S0006-3495(72)86068-5
Wilson, H. R. y Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik, 13(2), 55–80. doi:10.1007/BF00288786