¿Qué demostraría que una dinámica cerebral es crítica?

482.1. Qué afirma una hipótesis de criticalidad

Un registro de actividad cerebral contiene episodios de tamaños muy dispares y, al representar su frecuencia en ejes logarítmicos, aparece un tramo aproximadamente recto. Ese hallazgo puede motivar una hipótesis de criticalidad, aunque la hipótesis dice bastante más que la existencia de eventos grandes y pequeños. Sostiene que el sistema opera cerca de una transición dinámica con propiedades particulares de propagación, respuesta y escalamiento. Evaluarla exige precisar qué transición se postula y qué observaciones permitirían distinguirla de otros mecanismos capaces de generar distribuciones parecidas.

La palabra criticalidad recibe definiciones precisas dentro de cada familia de modelos. En un proceso de ramificación ideal, el punto crítico separa la extinción segura de la actividad de su posible crecimiento sin límite. En otros sistemas designa una transición entre fases con un parámetro de orden y relaciones de escalamiento propias. Trasladar propiedades de una familia a otra requiere justificar las equivalencias pertinentes. Una señal compatible con ramificación crítica no garantiza que el circuito comparta todas las propiedades que la física estadística atribuye a los sistemas críticos, y menos aún que esas propiedades optimicen alguna función cerebral.

Beggs y Plenz (2003) registraron potenciales de campo local con una matriz de 60 electrodos en cultivos organotípicos de corteza de rata y en cortes agudos de corteza de rata adulta. Con una ventana temporal igual al intervalo medio entre eventos, encontraron una distribución de tamaños de avalancha con exponente cercano a −3/2, duraciones con pendiente inicial próxima a −2 y un parámetro de ramificación estimado cercano a uno. El trabajo incluyó manipulaciones farmacológicas y comparaciones entre subconjuntos de electrodos que se comentan más adelante. Se trató de preparaciones in vitro, de modo que extender sus conclusiones a la dinámica cerebral en general exige evidencia adicional sobre escalas, modo de observación y condiciones de funcionamiento.

El capítulo reconstruye una prueba más exigente que reconocer una recta en un gráfico, a partir de un modelo de ramificación y de sus rivales. Los números proceden de ejemplos sintéticos construidos para aislar un argumento y no constituyen evidencia sobre tejido nervioso ni evidencia clínica.

482.2. Definir el evento que se analiza

Una avalancha observada suele definirse como una secuencia de intervalos temporales con actividad, delimitada antes y después por intervalos sin actividad. Esa definición convierte un registro continuo, o una sucesión de eventos, en unidades analizables. El tamaño puede contar eventos, canales activos o amplitud acumulada; la duración puede contar intervalos o expresarse en tiempo físico. Estas opciones no son intercambiables, y antes de comparar exponentes hay que saber qué magnitud se midió y cómo el procedimiento de detección transformó la señal en episodios.

El ancho del intervalo temporal puede unir eventos de procesos distintos o fragmentar uno prolongado. Beggs y Plenz informaron que el exponente de tamaños variaba con ese ancho y que el valor cercano a −3/2 aparecía al usar, en cada preparación, el intervalo medio entre eventos, unos 4 ms en promedio. Ese criterio se fija sin mirar el exponente, lo que lo separa del procedimiento inaceptable de ensayar anchos hasta obtener el valor esperado. Aun así, el exponente informado queda ligado a una regla de discretización. Dos estudios solo son comparables si sus reglas son equivalentes, y la elección debe acompañarse de un análisis de sensibilidad.

El umbral de detección introduce otra transformación. Una señal de campo local resume la actividad sináptica de poblaciones y depende de la geometría del registro; contar sus deflexiones no equivale a contar las descargas de todas las neuronas. Cambiar el umbral modifica el número de eventos y la estructura de los silencios. Que un patrón resista elecciones razonables es informativo, aunque no aporta especificidad causal, porque un mecanismo alternativo también puede ser robusto frente a esos cambios.

El tamaño del sistema observado impone otros límites. Una matriz de electrodos no cubre todas las conexiones, de modo que un intervalo sin señal puede indicar ausencia de actividad detectada en la región cubierta sin que la actividad del circuito se haya extinguido. Un corte en la distribución puede reflejar el tamaño físico de la red, la ventana de observación o la sensibilidad instrumental. Beggs y Plenz compararon subconjuntos de 60, 30 y 15 electrodos con la misma separación y observaron que el corte aparecía en el número máximo de electrodos de cada subconjunto, lo que vincula el corte con la cobertura del instrumento.

482.3. Un proceso de ramificación como hipótesis concreta

Consideremos un evento inicial que produce un número aleatorio de descendientes, cada uno de los cuales repite la misma regla de manera independiente. Si la media de descendientes es menor que uno, la actividad se extingue con probabilidad uno y su tamaño total esperado es finito. Si la supera, existe una probabilidad positiva de crecimiento indefinido en el modelo ideal. En el valor uno la extinción sigue siendo segura cuando la descendencia tiene varianza positiva, pero el tamaño total esperado diverge. Esa media, el parámetro de ramificación $m$ , resume una regla generadora concreta y no debe confundirse con cualquier cociente entre conteos observados sucesivos.

El tamaño total esperado se obtiene sumando generaciones: la inicial contiene un evento; la siguiente, $m$ eventos esperados; la posterior, $m^2$ . El resultado depende solo de la media y vale para cualquier distribución de descendencia, no únicamente para la de Poisson. La suma geométrica converge cuando $m<1$ . A medida que $m$ se acerca a uno, una perturbación pequeña genera una respuesta acumulada mayor y más persistente.

E[Z_t\mid Z_0=1]=m^t,\qquad E[S]=\sum_{t=0}^{\infty}m^t=\frac{1}{1-m},\qquad 0\leq m<1.

Con $m$ igual a 0,6, 0,9 y 0,98, los tamaños esperados son 2,5, 10 y 50 eventos. La varianza sí depende de la descendencia: si esta tiene varianza $\sigma^2$ , la de $S$ es $\sigma^2/(1-m)^3$ . Para descendencia de Poisson, $\sigma^2=m$ y las varianzas son exactamente 9,375, 900 y 122 500. Esa dispersión creciente tiene consecuencias de diseño, porque pocas perturbaciones entregan una estimación muy inestable de la respuesta media, y un promedio pequeño no descarta una respuesta esperada grande si los episodios raros que la sostienen no aparecieron durante el registro.

\operatorname{Var}(S)=\frac{m}{(1-m)^3}\ \ (\text{descendencia de Poisson}),\qquad m^t=e^{-t/\tau},\qquad \tau=-\frac{1}{\ln m},\quad 0<m<1.

El tiempo característico de decaimiento de la respuesta media está definido para $0<m<1$ ; tiende a cero cuando $m$ tiende a cero y diverge cuando $m$ tiende a uno. En unidades de generación, los tres valores del ejemplo dan aproximadamente 1,96, 9,49 y 49,50. Si una generación corresponde a un intervalo $\Delta t$ , el tiempo físico es $-\Delta t/\ln m$ ; con $\Delta t=4$ ms y $m=0{,}98$ resultan unos 198 ms. Una estimación de $m$ que no declara el intervalo de discretización no permite comparar registros adquiridos con resoluciones distintas.

El punto crítico tiene además una firma distribucional. Con descendencia de Poisson de media uno, el tamaño total sigue la distribución de Borel, cuya probabilidad para un tamaño $n$ se aproxima a $n^{-3/2}/\sqrt{2\pi}$ . Si la descendencia crítica tiene varianza finita y positiva $\sigma^2$ , la probabilidad de que la actividad sobreviva más allá de la generación $t$ se comporta asintóticamente como $2/(\sigma^2 t)$ , lo que corresponde a una densidad de duraciones con exponente −2. De ahí provienen los valores −3/2 y −2 con que suelen compararse los datos de avalanchas, y su aparición en un registro no basta para identificar el proceso.

482.4. La observación finita modifica la prueba

Todo experimento registra durante un horizonte finito. Si se observan las generaciones cero a cien, la respuesta acumulada esperada con $m=0{,}98$ es aproximadamente 43,50, inferior a 50; la diferencia corresponde a actividad posterior a la ventana. En el punto crítico la respuesta media por generación permanece en uno y la suma hasta ese horizonte vale 101. Una ventana finita puede volver comparables valores cuyas propiedades asintóticas difieren mucho, y el análisis debe tratar la censura de manera explícita.

E[S_H]=\sum_{t=0}^{H}m^t= \begin{cases}(1-m^{H+1})/(1-m),&m\neq1,\\H+1,&m=1.\end{cases}

La criticalidad exacta es una propiedad límite del modelo. En un circuito finito, la saturación, la refractariedad y los límites de recursos modifican la propagación, y no cabe esperar una divergencia literal. La hipótesis tiene que traducirse a predicciones sobre cómo cambian respuestas y cortes al modificar el tamaño del sistema o la distancia a la transición, en vez de exigir una singularidad inobservable o aceptar cualquier distribución amplia como evidencia.

La distribución de Borel muestra hasta qué punto un proceso subcrítico puede pasar por crítico. Para $m<1$ conserva el factor $n^{-3/2}$ , multiplicado por un corte exponencial de escala $1/(m-1-\ln m)$ , que con $m=0{,}98$ ronda los 4900 eventos. El cociente entre las probabilidades subcrítica y crítica vale prácticamente uno para cien eventos, cerca de 0,83 para mil y de 0,37 para cinco mil. Con avalanchas de algunos cientos de eventos, la forma de la distribución no separa ambos valores; la diferencia se manifiesta en el corte, justamente donde también actúan el tamaño del sistema y la ventana de observación.

También hay que distinguir una dinámica cercana al punto crítico de un mecanismo que se autoorganiza hacia él. La segunda afirmación propone procesos que alcanzan y mantienen ese régimen sin ajuste externo preciso, y un parámetro cercano a uno no los identifica; sostenerla exige examinar adaptación, recuperación tras cambios impuestos y la relación entre los mecanismos de regulación y el régimen alcanzado.

El primer bloque comprueba las cantidades de las secciones 482.3 y 482.4. El segundo verifica la asíntota crítica de Borel, el corte para $m=0{,}98$ y la ley de supervivencia de la duración.

from math import log, isclose

def response(m, horizon=100):
    """Respuesta media acumulada hasta el horizonte, iniciada por un evento."""
    return sum(m**t for t in range(horizon + 1))

for m, mean, var in [(.6, 2.5, 9.375), (.9, 10, 900), (.98, 50, 122500)]:
    assert isclose(1/(1 - m), mean)          # válido para cualquier descendencia con media m
    assert isclose(m/(1 - m)**3, var)        # requiere descendencia de Poisson (varianza = m)
    tau = -1/log(m)                          # dominio 0 < m < 1
    print(m, mean, var, round(tau, 2), round(response(m), 2))
assert isclose(response(.98), 43.50164176115707)
assert response(1) == 101

from math import exp, lgamma, log

def borel_logpmf(n, m):
    """log P(S=n) para el tamaño total con descendencia Poisson(m) (distribución de Borel)."""
    return -n*m + (n - 1)*log(n*m) - lgamma(n + 1)

# Punto crítico: P(S=n) se aproxima a n**(-3/2)/sqrt(2*pi)
for n in (10, 100, 1000, 10000):
    ratio = exp(borel_logpmf(n, 1.0)) * n**1.5 * (2*3.141592653589793)**.5
    print('n =', n, 'cociente con la asíntota:', round(ratio, 4))
assert abs(ratio - 1) < 1e-4

# Cerca del punto crítico: corte exponencial de escala 1/(m - 1 - log m)
m = .98
scale = 1/(m - 1 - log(m))
print('escala de corte para m = 0.98:', round(scale))
for n in (100, 1000, 5000, 20000):
    print(n, round(exp(borel_logpmf(n, m) - borel_logpmf(n, 1.0)), 4))

# Duración crítica: P(T > t) * t tiende a 2 (Poisson(1), varianza 1)
q = 0.0
for t in range(1, 10001):
    q = exp(q - 1)       # probabilidad de extinción antes de la generación t
print('t * P(T > t) con t = 10000:', round(10000*(1 - q), 4))
assert abs(10000*(1 - q) - 2) < .01

482.5. El submuestreo puede ocultar persistencia

Para representar actividad sostenida se añade una entrada externa. El número de eventos del paso siguiente sigue una distribución de Poisson cuya media es la actividad actual multiplicada por $m$ , más una entrada constante $h$ . Con $0\leq m<1$ y $h>0$ existe una media estacionaria finita, $\mu=h/(1-m)$ . Si solo se registra una fracción de los eventos, la serie observada suma a la variación de la actividad otra introducida por el muestreo, y esta segunda fuente debilita la relación aparente entre pasos consecutivos.

Z_{t+1}\mid Z_t\sim\operatorname{Poisson}(mZ_t+h), \qquad Y_t\mid Z_t\sim\operatorname{Binomial}(Z_t,q), \qquad \mu=\frac{h}{1-m}.

Con $m=0{,}9$ y $h=5$ , la media estacionaria es 50 y la varianza, aproximadamente 263,16. Al observar una décima parte de los eventos, la media registrada es 5 y la varianza, aproximadamente 7,13. La pendiente de la regresión de un conteo observado sobre el anterior es aproximadamente 0,3321, muy por debajo del parámetro verdadero. Leerla como medida directa de propagación haría parecer mucho menos persistente una dinámica cuyo tiempo característico es de unos 9,5 pasos.

El sesgo tiene una explicación algebraica que descansa en un supuesto que conviene enunciar completo. No basta con que cada conteo observado sea binomial dado el conteo real del mismo instante. Se requiere que, dado el recorrido completo de la actividad, los conteos observados en instantes distintos sean condicionalmente independientes, que cada uno dependa solo de la actividad de su instante y que el muestreo no influya sobre la dinámica. Bajo esos supuestos, la ley de la covarianza total anula la covarianza condicional entre instantes distintos, y la covarianza observada para rezagos positivos es $q^2$ veces la de la actividad real, mientras la varianza observada incorpora además ruido binomial. Como la covarianza de la actividad decae como $m^k$ , el cociente entre las covarianzas de segundo y primer rezago recupera $m$ en la población. Con datos finitos ese cociente es inestable, y ajustar la caída a lo largo de varios rezagos suele ser más informativo.

V_Z=\frac{\mu}{1-m^2},\quad V_Y=q^2V_Z+q(1-q)\mu,\quad \operatorname{Cov}(Y_t,Y_{t+k})=q^2m^kV_Z\quad(k\geq1).

Si la fracción observada cambia de forma persistente, como cuando un electrodo pierde contacto durante un período o la detectabilidad depende del estado del tejido, aparece un término de covarianza ajeno a la dinámica y la fórmula deja de valer.

Wilting y Priesemann (2018) desarrollaron un estimador de regresiones a múltiples pasos para este problema. Para cada rezago $k$ calculan la pendiente $r_k$ de la regresión de la actividad observada en $t+k$ sobre la observada en $t$ y ajustan las pendientes al modelo $r_k=b\,m^k$ , donde $b$ absorbe la atenuación por muestreo y es común a todos los rezagos. Suponen que la actividad completa admite una representación autorregresiva de primer orden y que la observada es, en esperanza, una función lineal de la completa. Advierten que el método requiere series largas y procesos bien aproximados por esa representación, e incluyen pruebas de aplicabilidad y de no estacionariedad. Su aporte consiste en mostrar que una inferencia aparentemente directa cambia cuando se representa de manera explícita cómo se observa el sistema.

482.6. Una simulación con observación parcial

La simulación genera doscientos mil pasos tras descartar mil de estabilización y observa, con muestreo binomial independiente en cada paso, una décima parte de cada conteo. Con la semilla fija, la media observada es aproximadamente 4,9947, la pendiente de un paso 0,3321 y el cociente entre covarianzas de segundo y primer rezago 0,9074, próximo al parámetro 0,9. El ajuste log-lineal de las diez primeras pendientes entrega $m\approx0{,}8998$ y $b\approx0{,}3696$ , frente a un valor teórico de 0,3690 para $b$ .

import numpy as np

rng = np.random.default_rng(482)
m = .9; h = 5.; q = .1; n = 200000; burn = 1000
z = np.empty(n + burn, dtype=int); z[0] = 50
for t in range(1, len(z)):
    z[t] = rng.poisson(m*z[t - 1] + h)
z = z[burn:]
# Muestreo binomial independiente en cada paso, dado el recorrido completo de Z
y = rng.binomial(z, q).astype(float)

def acov(x, k):
    c = x - x.mean()
    return np.mean(c*c) if k == 0 else np.mean(c[:-k]*c[k:])

c0, c1, c2 = acov(y, 0), acov(y, 1), acov(y, 2)
mu = h/(1 - m); vz = mu/(1 - m*m); vy = q*q*vz + q*(1 - q)*mu
naive = q*q*m*vz/vy
assert np.isclose(naive, .33210332103321044)
assert abs(c1/c0 - naive) < .01
assert abs(c2/c1 - m) < .03
print('Media observada:', y.mean())
print('Pendiente un paso:', c1/c0, '; razón de covarianzas:', c2/c1)

# Regresión a múltiples pasos: r_k = b * m**k, k = 1..10, ajuste log-lineal
ks = np.arange(1, 11)
r = np.array([acov(y, k)/c0 for k in ks])
slope, intercept = np.polyfit(ks, np.log(r), 1)
print('m estimado:', np.exp(slope), '; b estimado:', np.exp(intercept), '; b teórico:', q*q*vz/vy)
assert abs(np.exp(slope) - m) < .01

El segundo bloque conserva la dinámica, pero hace que la fracción observada alterne entre 0,05 y 0,15 según una cadena de dos estados que cambia, en promedio, cada cincuenta pasos. Cada conteo es binomial condicionado al conteo real y a la fracción observada de ese instante. Condicionado solo a la actividad real, es una mezcla de binomiales y el muestreo tiene dependencia temporal. El cociente entre covarianzas de segundo y primer rezago sube a 0,934, y el de quinto y cuarto, a 0,940, porque la covarianza de la fracción observada decae como $0{,}96^k$ , más despacio que $0{,}9^k$ , y se suma a la de la dinámica. Un estimador que supone muestreo independiente atribuiría esa persistencia instrumental al circuito y sobrestimaría su cercanía al punto crítico.

import numpy as np

# Contraejemplo: la fracción observada cambia de forma persistente en el tiempo.
# Y_t | Z_t, q_t es binomial. Dado solo Z, hay mezcla binomial y dependencia temporal.
rng = np.random.default_rng(4820)
m = .9; h = 5.; n = 200000; burn = 1000
z = np.empty(n + burn, dtype=int); z[0] = 50
for t in range(1, len(z)):
    z[t] = rng.poisson(m*z[t - 1] + h)
z = z[burn:]
switch = rng.random(n) < .02
state = np.cumsum(switch) % 2                  # cadena de dos estados, permanencia media de 50 pasos
q_t = np.where(state == 0, .05, .15)
y = rng.binomial(z, q_t).astype(float)

def acov(x, k):
    c = x - x.mean()
    return np.mean(c[:-k]*c[k:])

for k in (1, 4):
    print('razón de covarianzas', k + 1, '/', k, ':', acov(y, k + 1)/acov(y, k))
assert acov(y, 2)/acov(y, 1) > .92            # sesgo hacia arriba respecto de m = 0,9

482.7. Una cola de potencia sin propagación entre unidades

Una distribución amplia también puede surgir de la heterogeneidad de tasas. Para cada observación se extrae una tasa positiva de una distribución gamma y luego un tiempo de espera exponencial con esa tasa. Las tasas pequeñas producen tiempos muy largos. No hay red de descendencia ni transición de propagación y, sin embargo, la mezcla de exponenciales genera una cola de potencia, aunque cada distribución condicionada tenga una escala característica.

Con forma 1/2 y tasa unitaria, la probabilidad de que el tiempo supere $t$ es $(1+t)^{-1/2}$ y la densidad decae asintóticamente con exponente −3/2. La coincidencia con el exponente discutido para tamaños de avalancha no establece una identidad de fenómenos, porque aquí se modelan tiempos de espera. El ejemplo demuestra algo más acotado: la forma de una cola no identifica por sí sola su mecanismo. Priesemann y Shriki (2018) observan algo análogo en unidades de Poisson sin interacción con una tasa común que varía lentamente, cuyos intervalos entre eventos forman una mezcla de exponenciales que puede aproximarse a una ley de potencia.

\Lambda\sim\operatorname{Gamma}(1/2,\text{tasa}=1),\qquad T\mid\Lambda\sim\operatorname{Exponencial}(\Lambda),

\Pr(T>t)=\int_0^\infty e^{-\lambda t}f_\Lambda(\lambda)\,d\lambda =(1+t)^{-1/2},\qquad f_T(t)=\tfrac12(1+t)^{-3/2}.

La variable $1+T$ sigue una distribución de Pareto exacta sobre los valores mayores o iguales a uno, con exponente de densidad 3/2. El estimador continuo de máxima verosimilitud de la sección 482.9, con umbral igual a uno, entrega un valor cercano a 1,5, lo que solo confirma que el ajuste reconoce la distribución generada. Como ese exponente es menor que dos, $T$ tiene además media infinita: la media divergente que en el proceso de ramificación aparece en el punto crítico surge aquí de la heterogeneidad de tasas.

En una muestra sintética de cien mil observaciones, las proporciones de tiempos mayores que 3 y que 99 son 0,50322 y 0,10034, frente a valores poblacionales de 0,5 y 0,1, y el exponente estimado es aproximadamente 1,4979, con un error estándar cercano a 0,0016. Una señal estadística puede ser muy nítida y seguir siendo inespecífica respecto de la hipótesis física que interesa.

import numpy as np

rng = np.random.default_rng(1482)
rates = rng.gamma(.5, 1, 100000)                # forma 1/2, escala 1 (tasa 1)
times = rng.exponential(1/rates)                # escala = 1/tasa
alpha_hat = 1 + len(times)/np.log1p(times).sum()
for threshold, probability in [(3, .5), (99, .1)]:
    estimate = np.mean(times > threshold)
    se = np.sqrt(probability*(1 - probability)/len(times))
    assert abs(estimate - probability) < 6*se
    print('Supervivencia', threshold, estimate, 'teoría', probability)
assert abs(alpha_hat - 1.5) < .02
print('Exponente de densidad estimado:', alpha_hat, '; error estándar aproximado:', (alpha_hat - 1)/np.sqrt(len(times)))

482.8. Alternativas más próximas a la observación neuronal

Para evaluar registros neuronales se necesitan rivales biológicamente plausibles que reproduzcan también el procedimiento de observación. Touboul y Destexhe (2017) estudiaron una red de neuronas de integración y disparo con excitación e inhibición, el modelo de Brunel, en un régimen sincrónico irregular alejado de cualquier transición. Sus avalanchas de descargas mostraron exponentes de tamaño y duración cercanos a −1,42 y −2,11, además de colapso de las formas promedio. Como control construyeron un sustituto de unidades independientes: procesos de Poisson sin interacción directa cuya tasa común era la parte positiva de un proceso de Ornstein-Uhlenbeck. El sustituto también produjo leyes de potencia aparentes, con exponentes cercanos a −1,47 y −1,9, y colapso de formas.

La conclusión de ese trabajo se refiere a la suficiencia de los indicadores: si un mecanismo no crítico reproduce una señal, la señal aislada no identifica criticalidad. De ahí no se sigue que todos los registros cerebrales sean no críticos ni que una hipótesis concreta quede descartada sin comparación específica. Priesemann y Shriki (2018) analizaron unidades de Poisson sin interacción sometidas a una tasa común variable y encontraron que la tasa debe variar lentamente respecto de la duración de las avalanchas y recorrer un intervalo amplio, y que aun así solo se obtienen leyes de potencia aproximadas, con cortes comparables a los experimentales. Propusieron como diferencias medibles una dependencia del ancho de ventana mucho mayor que la de un proceso de ramificación cercano al punto crítico y una relación prácticamente lineal entre tamaño medio y duración.

Un control adecuado respeta las propiedades que no se pretende poner a prueba. Aleatorizar por completo una serie elimina a la vez las tasas variables, la dependencia temporal y la organización espacial, y cualquier diferencia con ese control admite demasiadas explicaciones. Resultan más informativos los controles sucesivos que conservan unas propiedades y alteran otras. Mantener la tasa compartida y eliminar la interacción directa, como en el sustituto de Touboul y Destexhe, responde a una pregunta distinta de la que se examina al conservar las interacciones y cambiar la conectividad. Esos autores advierten, además, que usar como tasa del sustituto la extraída de la propia red impediría saber si las leyes de potencia quedaron codificadas en esa función.

La comparación debe incluir también cómo se producen las señales agregadas. Touboul y Destexhe calcularon señales semejantes a potenciales de campo a partir de las corrientes sinápticas de otra red simulada y obtuvieron estadísticas de picos con escalamiento parecido en un régimen sincrónico, en uno asincrónico y en un sustituto de Poisson de tasa constante, con exponentes cercanos a −1,4. Si un experimento solo dispone de esa señal macroscópica, su capacidad de distinguir mecanismos queda muy limitada, y la salida puede requerir observaciones complementarias, perturbaciones o una afirmación más restringida. Hacer más complejo el modelo sin mejorar su vínculo con el instrumento añade parámetros sin recuperar información discriminante.

482.9. Ajustar una ley de potencia de manera defendible

Clauset, Shalizi y Newman (2009) explican por qué la inspección visual en ejes logarítmicos no basta. Muchas distribuciones parecen rectas en un rango limitado, y la regresión por mínimos cuadrados sobre histogramas transformados puede producir errores sistemáticos grandes en el exponente que los errores de la propia regresión no revelan. El análisis debe declarar si el modelo se propone para toda la distribución o solo para la cola, porque esa decisión determina cuántos datos participan y qué conclusión cabe extraer. Un buen ajuste de los eventos más grandes no describe necesariamente cómo se generan los pequeños.

Para una cola continua con umbral conocido $x_{\min}$ , el estimador de máxima verosimilitud usa los logaritmos de las observaciones divididas por el umbral, con error estándar aproximado $(\hat\alpha-1)/\sqrt{n}$ . En datos discretos, como los conteos de electrodos, el estimador exacto exige resolver numéricamente una ecuación con la función zeta de Hurwitz; la aproximación habitual reemplaza el umbral por $x_{\min}-1/2$ y, según los autores, tiene un error cercano al uno por ciento o menor cuando $x_{\min}$ es del orden de seis o mayor. Aplicar la fórmula continua a tamaños enteros pequeños puede sesgar el resultado. Cuando el umbral se elige con los datos, los autores recomiendan minimizar la distancia de Kolmogorov-Smirnov entre las distribuciones empírica y ajustada, y su incertidumbre forma parte del resultado.

\hat\alpha=1+n\left[\sum_{i=1}^{n}\ln\frac{x_i}{x_{\min}}\right]^{-1}\ \ (\text{continuo}),\qquad \hat\alpha\simeq1+n\left[\sum_{i=1}^{n}\ln\frac{x_i}{x_{\min}-\tfrac12}\right]^{-1}\ \ (\text{discreto, aproximado}).

La bondad de ajuste se contrasta con un remuestreo semiparamétrico en el que cada muestra sintética repite la estimación completa, incluida la del umbral; usar como referencia una distribución que supone parámetros conocidos, cuando se ajustaron con la misma muestra, es inadecuado. En el criterio de los autores, un valor $p$ menor o igual que 0,1 descarta la ley de potencia, y un valor mayor solo indica que no pudo descartarse con la potencia disponible. El procedimiento supone, además, observaciones independientes e idénticamente distribuidas. Touboul y Destexhe señalan que estos métodos, eficaces con datos estacionarios, fallan al evaluar las colas de datos no estacionarios, situación frecuente cuando miles de avalanchas proceden de un mismo registro prolongado.

Después deben compararse alternativas pertinentes. Clauset y colaboradores emplean cocientes de verosimilitud, con la prueba de Vuong para modelos no anidados, frente a las distribuciones lognormal, exponencial y exponencial estirada, y una prueba adaptada para la ley de potencia con corte exponencial, que es un modelo anidado. Cuando el cociente no es significativo, el resultado se declara no concluyente. Una preferencia estadística es siempre relativa al conjunto comparado, y elegir una distribución por costumbre o por afinidad con una teoría atractiva no resuelve la comparación. En registros limitados, la conclusión más precisa puede ser que existe una cola amplia cuya forma exacta y cuyo mecanismo aún no se han identificado.

482.10. Perturbar para observar respuesta

El modelo de ramificación predice cómo evoluciona la actividad después de un evento adicional, lo que ofrece un contraste distinto de la distribución espontánea de avalanchas. Priesemann y Shriki sostienen que separar las propiedades internas del sistema de las impuestas desde fuera requiere, en último término, intervenciones causales, porque una entrada externa suficientemente elaborada podría imitar actividad crítica ante cualquier análisis correlacional. La perturbación debe relacionarse con el parámetro de control propuesto. Si modifica muchas propiedades a la vez, no puede atribuirse cualquier cambio de respuesta a la distancia respecto del punto crítico.

Las manipulaciones de Beggs y Plenz ilustran esa dificultad. En los cultivos, la picrotoxina, antagonista de receptores GABA_A, redujo la inhibición, destruyó la forma de potencia y produjo una distribución bimodal con una joroba de eventos grandes que los autores describieron como descargas epileptiformes; tras el lavado, el exponente inicial volvió a valores cercanos a −3/2. El resultado es compatible con un desplazamiento hacia un régimen supercrítico, pero la intervención alteró el balance entre excitación e inhibición y no solamente un parámetro de ramificación. En los cortes agudos, por otra parte, la actividad se indujo con NMDA y un agonista dopaminérgico D1, de modo que el régimen observado dependía de esa inducción.

La comparación debería incluir intensidades pequeñas, para examinar la respuesta local, y condiciones en que aparezca saturación. Una respuesta grande puede reflejar cercanía a una transición, pero también excitabilidad elevada o amplificación no crítica. La forma temporal de la respuesta, su dependencia del estado inicial y su cambio al modificar parámetros restringen más, siempre que el modelo prediga un conjunto conectado de resultados; acumular indicadores favorables sin una teoría cuantitativa que los vincule crea una apariencia de convergencia sin identificar un mecanismo común.

La recuperación necesita una referencia clara. Volver al promedio no equivale siempre a recuperar la misma dinámica, porque pueden haber cambiado las correlaciones, la variabilidad o la composición de la actividad. Una recuperación lenta tampoco es específica de criticalidad y puede proceder de entradas lentas o de adaptación. Una prueba útil distinguiría si la escala temporal aumenta al acercarse a la frontera postulada o si permanece fijada por una entrada externa que varía con independencia del circuito.

La repetición de perturbaciones puede, además, modificar el sistema por adaptación, fatiga o cambios de conectividad, de modo que el orden y los intervalos entre ensayos forman parte del diseño y una respuesta promedio que mezcla etapas distintas puede ocultar ese cambio.

482.11. Escalamiento y tamaño finito

Una hipótesis de escalamiento relaciona observables obtenidos con distintos tamaños o resoluciones y es, por eso, más exigente que ajustar una pendiente en cada registro. El modelo debe especificar cómo cambian los cortes, las amplitudes o las formas al modificar la escala. Si se ajustan libremente muchos factores de reescalamiento, la superposición de curvas aporta poca información. Conviene declarar qué exponentes se estimaron y cuáles estaban predichos, y conservar grados de libertad suficientes para que una discrepancia pueda contradecir la propuesta.

Cambiar la cantidad de electrodos observados no equivale necesariamente a cambiar el tamaño físico del sistema. La primera operación altera el muestreo; la segunda puede modificar la dinámica y las condiciones de borde. En la comparación de Beggs y Plenz con 60, 30 y 15 electrodos la red registrada era la misma, de modo que el desplazamiento del corte informa sobre la observación antes que sobre el tamaño del circuito. El submuestreo puede introducir cortes y sesgos aunque el circuito permanezca idéntico, y una evaluación sólida compara las predicciones del modelo dinámico con las del modelo de observación.

Las relaciones entre exponentes descansan también en supuestos. Si el tamaño de una avalancha de duración $T$ se concentra alrededor de una media que crece como $T^{\gamma}$ , un cambio de variables liga los tres exponentes. Con −3/2 y −2, la relación predice $\gamma=2$ .

P(S)\propto S^{-\alpha_S},\qquad P(T)\propto T^{-\alpha_T},\qquad E[S\mid T]\propto T^{\gamma}\ \Longrightarrow\ \gamma=\frac{\alpha_T-1}{\alpha_S-1}.

Comprobar esa coherencia aporta más que una pendiente aislada. Touboul y Destexhe encontraron que en la red de Brunel y en su sustituto el tamaño medio crecía con exponentes cercanos a 1,5 y 1,4, pese a que las formas promedio colapsaban; en el proceso con tasa variable de Priesemann y Shriki el exponente es prácticamente uno. En esos casos la relación discrimina donde la forma de las distribuciones no lo hace. Touboul y Destexhe indican, sin embargo, que el exponente que obtienen es consistente con el que citan para datos de cultivos in vitro. Si así fuera, la desviación respecto de la relación no afectaría solo a los modelos no críticos y obligaría a precisar qué clase de transición se postula para los datos.

La incertidumbre debe considerar además la dependencia entre eventos y entre registros. Miles de avalanchas de una sola preparación no equivalen a miles de sistemas biológicos independientes, y la variación entre episodios, sesiones y preparaciones condiciona el diseño del remuestreo y la lectura de los intervalos.

482.12. Qué afirmación permite cada resultado

Una distribución compatible con una ley de potencia sostiene una afirmación estadística sobre un rango observado. Una respuesta temporal compatible con un modelo de ramificación añade información dinámica, y una comparación que descarta las entradas compartidas como explicación suficiente restringe los rivales. El informe debe indicar qué pasos se completaron, sin condensarlos en una etiqueta que sugiera una identificación más fuerte que la permitida por la evidencia.

La cercanía a un valor crítico debe acompañarse de su incertidumbre y de la resolución temporal con que se estimó. Una estimación de 0,98 no equivale a uno, aunque ambos valores sean difíciles de distinguir con datos limitados. Wilting y Priesemann aplicaron su estimador a registros de descargas en corteza prefrontal de mono, corteza visual de gato e hipocampo de rata y obtuvieron valores entre 0,963 y 0,998, con mediana de 0,984 y tiempos de autocorrelación de cien a dos mil milisegundos. Describieron ese régimen como reverberante, distinto tanto del asincrónico irregular como del crítico. Afirmar criticalidad exacta, autoorganización o ventaja funcional exige argumentos que una buena estimación del decaimiento temporal no contiene.

Las afirmaciones funcionales necesitan medir la función propuesta. Beggs y Plenz complementaron sus registros con simulaciones de redes de propagación hacia adelante en las que la transmisión de información era máxima con un parámetro de ramificación cercano a uno. Un resultado así muestra que cierto régimen maximiza una medida bajo las reglas del modelo, pero no que el tejido observado la optimice ni que esa optimización se traduzca en mejor funcionamiento clínico. Para eso habría que definir entradas, salidas y costos y contrastarlos en el sistema de interés. Mantener separadas la descripción dinámica y la valoración funcional evita usar una propiedad física como sinónimo de salud.

482.12.1. Un resultado negativo también puede discriminar

Supóngase que el análisis encuentra una cola amplia, pero la respuesta provocada decae con una escala temporal que apenas cambia cuando se modifica el parámetro propuesto como control. Si el modelo había predicho una variación apreciable en ese intervalo experimental, el resultado debilita esa formulación, aunque la distribución espontánea siga pareciendo atractiva. Antes de descartarla conviene examinar la potencia de la prueba y la fidelidad de la manipulación. Si ambas fueron suficientes, mantener la explicación exige una modificación explícita con nuevas predicciones; invocar indefinidamente factores no medidos convertiría la hipótesis en una descripción capaz de acomodarse a cualquier observación.

Un resultado nulo significa otra cosa si la manipulación no alteró la variable relevante. Modificar una concentración externa no garantiza un cambio conocido y uniforme en la ganancia efectiva de propagación, porque el circuito puede compensarlo. El experimento debe distinguir la intervención aplicada del cambio fisiológico conseguido, y medir un indicador intermedio ayuda a evaluar ese vínculo, aunque introduce su propio error. Explicar un resultado negativo exige precisar qué etapa se observó y cuál sigue siendo supuesta.

La comparación entre laboratorios requiere registrar las decisiones de procesamiento que alteran el objeto estadístico, porque dos grupos pueden llamar avalanchas a objetos construidos con umbrales, referencias eléctricas o reglas de concatenación distintas. La reproducción debe incluir la adquisición y la transformación de la señal, además del programa de ajuste, y las simulaciones del instrumento no sustituyen una validación independiente en nuevas preparaciones.

También hay que separar el fracaso de una predicción del fracaso de una familia teórica. Un modelo de ramificación homogéneo puede fallar porque ignora la inhibición o la heterogeneidad espacial, mientras otra teoría de transición sigue siendo plausible. Esa posibilidad justifica formular alternativas concretas, pero no autoriza a leer cualquier desviación como favorable. Una ampliación del modelo es aceptable cuando explica discrepancias identificadas y genera consecuencias nuevas que pueden ponerse a prueba. Un informe útil termina indicando qué propiedad se detectó, qué modelos se compararon y qué resultado futuro cambiaría la interpretación.

482.13. Alcance de la conclusión

Demostrar una hipótesis concreta de criticalidad exige un modelo lo bastante definido para producir predicciones discriminantes y un diseño que observe sus consecuencias pertinentes. La distribución de avalanchas participa en esa evaluación, pero debe relacionarse con la propagación, la respuesta a perturbaciones y los efectos de escala. Los ejemplos mostraron que la observación parcial puede ocultar persistencia y que un muestreo con memoria propia puede fabricarla; que un proceso subcrítico puede imitar la distribución crítica en un rango amplio; y que una cola de potencia puede surgir sin ramificación. También mostraron bajo qué supuestos varios rezagos permiten recuperar información dinámica.

Una conclusión defendible podría establecer que una preparación, en condiciones determinadas, es compatible con una dinámica cercana a una transición y que los rivales examinados explican peor ciertas respuestas, dejando registrado qué alternativas siguen abiertas y qué observaciones cambiarían el juicio. Si la evidencia se reduce a una recta en un gráfico logarítmico, corresponde una conclusión estadística más modesta. Lo que separa ambas situaciones es el número de mecanismos que el diseño permite distinguir y la calidad de las mediciones que conectan el modelo con la actividad registrada.

Referencias

Beggs, J. M., y Plenz, D. (2003). Neuronal avalanches in neocortical circuits. The Journal of Neuroscience, 23(35), 11167–11177. https://doi.org/10.1523/JNEUROSCI.23-35-11167.2003

Clauset, A., Shalizi, C. R., y Newman, M. E. J. (2009). Power-law distributions in empirical data. SIAM Review, 51(4), 661–703. https://doi.org/10.1137/070710111. Lectura asistida de la versión de autor arXiv:0706.1062.

Priesemann, V., y Shriki, O. (2018). Can a time varying external drive give rise to apparent criticality in neural systems? PLOS Computational Biology, 14(5), e1006081. https://doi.org/10.1371/journal.pcbi.1006081

Touboul, J., y Destexhe, A. (2017). Power-law statistics and universal scaling in the absence of criticality. Physical Review E, 95(1), 012413. https://doi.org/10.1103/PhysRevE.95.012413. Lectura asistida de la versión de autor arXiv:1503.08033.

Wilting, J., y Priesemann, V. (2018). Inferring collective dynamical states from widely unobserved systems. Nature Communications, 9(1), 2325. https://doi.org/10.1038/s41467-018-04725-4