Modelos cerebrales de gran escala: conectar estructura y dinámica

Un conectoma describe relaciones entre regiones, pero por sí solo no dice cómo cambia su actividad en el tiempo. Para obtener una trayectoria hay que asignar a cada región una dinámica local, definir cómo se transmite la influencia entre ellas y precisar qué registra el instrumento de medida. Una misma matriz anatómica puede sostener comportamientos muy distintos cuando cambian las ganancias, los retardos o las entradas. Construir un modelo cerebral consiste, en buena parte, en articular esas decisiones de manera explícita.

Imaginemos dos modelos que reproducen igualmente bien la correlación entre señales regionales en reposo. Uno emplea un acoplamiento fuerte con ruido independiente en cada región; el otro combina un acoplamiento débil con una entrada común. La calidad parecida del ajuste deja abierta la pregunta sobre el mecanismo que generó los datos. Una perturbación podría distinguirlos si altera una región y permite seguir cómo llega su efecto a las demás. El modelo generativo sirve para calcular esa predicción, siempre que la observación esté representada con fidelidad suficiente. La relación entre ambos planos tiene respaldo empírico y también límites claros. En un trabajo que midió conectividad estructural y funcional en las mismas personas, la estructura anatómica restringía la intensidad y la organización espacial de la conectividad funcional en reposo, aunque aparecían conexiones funcionales intensas entre regiones sin enlace estructural directo (Honey et al., 2009).

Plataformas como The Virtual Brain organizaron la simulación de redes cerebrales completas combinando conectividad estructural, modelos de dinámica regional y modelos de observación para distintas modalidades de imagen (Sanz-Leon et al., 2013). La formulación matemática de esa plataforma se describió después con una notación que hace corresponder cada componente de la ecuación con un componente del programa, incluidos los modelos directos que traducen la actividad simulada a EEG, MEG, EEG estereotáctico y resonancia funcional (Sanz-Leon et al., 2015). La utilidad de esta arquitectura depende menos del detalle de cada pieza que de la claridad de las interfaces entre ellas y de la posibilidad de averiguar qué parámetros quedan efectivamente respaldados por las observaciones. Las revisiones del campo presentan este enfoque como un marco capaz de integrar datos de varias modalidades y de producir predicciones que pueden ponerse a prueba o refutarse, sin desconocer los desafíos que siguen abiertos (Breakspear, 2017).

Construiremos un modelo pequeño con la misma arquitectura conceptual que uno de gran escala. Tres regiones bastan para separar estados, transmisión y medición. El ejemplo usará datos sintéticos y un filtro de observación deliberadamente simple. Su propósito es estudiar identificabilidad y respuesta a perturbaciones; en ningún momento pretende ser una simulación calibrada de un cerebro individual.

405.1. Las capas de una explicación generativa

Una formulación general puede escribirse como

Aquí $x_i$ es el estado de la región $i$, que puede reunir varias variables. La función $f_i$ describe su dinámica local. Un ejemplo clásico de esa dinámica es el par de ecuaciones de Wilson y Cowan para poblaciones localizadas de neuronas excitatorias e inhibitorias, cuyo análisis mostró histéresis y oscilaciones de ciclo límite según el tipo de estímulo (Wilson y Cowan, 1972). La matriz $C_{ij}$ organiza las conexiones y $H_{ij}$ especifica cómo la actividad del origen influye sobre el destino. $G$ es una escala global de acoplamiento, $d_{ij}$ un retardo de transmisión, $u$ una entrada conocida y $\eta$ una perturbación interna incluida en el modelo. Todos los términos deben expresarse en unidades compatibles.

La propia matriz $C$ depende de decisiones previas a la simulación. Un estudio sistemático de redes anatómicas obtenidas por difusión y tractografía comparó parcelaciones anatómicas y aleatorias en un amplio rango de escalas. Las decisiones gruesas sobre la organización de la red, como la presencia de propiedades de mundo pequeño, no dependieron de la escala. Los valores de parámetros como el agrupamiento, la longitud de camino o la eficiencia fueron consistentes entre parcelaciones con el mismo número de nodos y variaron considerablemente con ese número (Zalesky et al., 2010). Por analogía, un parámetro de acoplamiento estimado sobre una parcelación queda ligado a ella y su traslado a otra requiere justificación.

La observación exige una relación adicional,

El operador $\mathcal O$ puede depender de la historia de los estados, por ejemplo a través de un filtro. El vector $\phi$ reúne los parámetros de observación y $\varepsilon$ representa el error de medición. Mantener separadas ambas ecuaciones permite distinguir la variabilidad del proceso de la variabilidad del instrumento. Cuando las dos se introducen como un único término de ruido sin justificación, se pierde la posibilidad de anticipar cómo debería cambiar cada una al mejorar la adquisición.

Falta además una distribución o un valor para el estado inicial. En presencia de retardos no basta con fijar un vector en $t=0$, porque hay que especificar la historia completa durante un intervalo que cubra el mayor de los retardos. Dos simulaciones que comparten el valor final de esa historia pueden divergir, ya que los términos de acoplamiento consultarán más adelante instantes anteriores. La manera de inicializar forma parte del diseño del experimento numérico.

Las regiones pueden compartir la ecuación local y diferir en sus parámetros. Esa heterogeneidad puede representar diferencias conocidas de antemano o funcionar simplemente como flexibilidad adicional para el ajuste. Conviene no confundir ambas situaciones. Si una distribución regional de parámetros se estima libremente a partir del mismo registro que el modelo intenta explicar, su correspondencia con la biología necesita una validación independiente.

405.2. Una red lineal que permite ver todas las piezas

Para que la complejidad local no oculte el argumento, usaremos una sola variable por región:

Los estados representan desviaciones adimensionales respecto de una referencia, de modo que en general pueden tomar valores negativos. Elegimos constantes temporales de 10, 20 y 30 ms y la matriz

Cada fila suma uno. Esta normalización fija la influencia entrante total en la representación elegida, aunque al mismo tiempo borra las diferencias de fuerza total entre destinos. En un modelo empírico esa consecuencia tendría que justificarse. En nuestro caso facilita comparar el efecto de una ganancia global $G=0.6$ con el de constantes locales distintas.

La matriz normalizada no es simétrica, a pesar de que podría provenir de un camino anatómico simétrico antes de la normalización por filas. La asimetría de sus coeficientes, por lo tanto, no demostraría una asimetría de las conexiones biológicas, porque es una propiedad de la transformación utilizada. Esta observación prolonga la distinción, planteada en el capítulo 401, entre matriz estructural y operador dinámico.

Fijamos retardos de 5 ms para la transmisión entre las regiones 1 y 2 y de 10 ms entre las regiones 2 y 3. La historia inicial será nula. Aplicaremos a la primera región un pulso de amplitud uno entre los 20 y los 40 ms. La región 1 responde directamente al pulso, mientras que las otras dos solo reciben influencias transmitidas. El recorrido de la señal combina tiempos locales y retardos, de manera que el máximo de una respuesta no aparece simplemente después de sumar los tiempos de conducción.

Sin retardos, el sistema tiene matriz dinámica $A=T^{-1}(-I+GC)$, donde $T$ es la matriz diagonal de constantes temporales. Los discos de Gershgorin de $A$ tienen centros en −0.1, −0.05 y −0.033 ms$^{-1}$ y radios de 0.06, 0.03 y 0.02 ms$^{-1}$, respectivamente. Los tres quedan en el semiplano izquierdo porque, en cada fila, la magnitud del término de fuga supera la suma de los acoplamientos. Este argumento garantiza la estabilidad del sistema instantáneo. Para el sistema con retardos hace falta un análisis propio, que la sección siguiente aborda.

405.3. Cómo entra un retardo en la estabilidad

Busquemos soluciones proporcionales a $e^{st}$ en el sistema lineal sin entrada. La condición característica toma la forma

Los términos exponenciales hacen que el problema deje de ser un polinomio de grado igual al número de regiones, porque la ecuación pasa a tener infinitas raíces asociadas a la historia. En circuitos con retroalimentación negativa, un retardo puede desfasar la corrección y favorecer oscilaciones. Su efecto depende de los signos y de las ganancias, y no todo retardo desestabiliza una red.

Nuestro ejemplo pertenece a una clase en la que los retardos no pueden desestabilizar. Si se divide cada ecuación por $\tau_i$, el sistema sin entrada toma la forma $\dot x(t)=A_0x(t)+A_1x(t-d_1)+A_2x(t-d_2)$, con $d_1=5$ ms y $d_2=10$ ms. La matriz $A_0=-T^{-1}$ es diagonal. Las matrices $A_1$ y $A_2$ tienen elementos no negativos porque reúnen los acoplamientos $GC_{ij}/\tau_i$ asociados a cada retardo. Su suma es $T^{-1}(-I+GC)$, cuyos valores propios están en el semiplano izquierdo. Para sistemas lineales positivos de este tipo, con matriz instantánea de Metzler, matrices retardadas no negativas y uno o varios retardos constantes, la estabilidad asintótica equivale a la del sistema sin retardos cuya matriz es la suma de todas ellas, de modo que no depende de los valores de los retardos (Kaczorek, 2009). Para el caso de un solo retardo, Haddad y Chellaboina (2004) dieron una condición necesaria y suficiente de estabilidad asintótica válida para todo valor del retardo, formulada también sobre la suma de la matriz instantánea y la retardada. En consecuencia, la red de tres regiones es estable para cualquier par de retardos constantes. Una simulación con los retardos multiplicados por 10 y por 50 lo ilustra, aunque la recuperación se vuelve más lenta a medida que los retardos crecen.

El alcance de esta garantía depende de la positividad del sistema. Con acoplamientos inhibitorios o de signo mixto, alguna matriz retardada deja de ser no negativa y el teorema ya no es aplicable. Eso no implica que la estabilidad pase a depender necesariamente del retardo. Implica que, si la independencia existe, debe establecerse por otra vía, por ejemplo analizando la ecuación característica. Con dinámicas locales no lineales, los resultados lineales pueden aplicarse a la linealización en torno a un equilibrio cuando esa linealización conserva la estructura positiva. En ese caso informan sobre la estabilidad local de ese equilibrio. La estabilidad lejos de él requiere herramientas específicas para sistemas no lineales con retardo.

El papel funcional de los retardos en modelos de gran escala es una cuestión abierta. En simulaciones con conectividad de primate, la combinación de ruido y retardos de propagación fue necesaria para que emergieran las fluctuaciones coherentes de la red por defecto (Ghosh et al., 2008). Una revisión posterior de tres modelos de neocorteza de primate atribuyó a la dinámica local, a los retardos de transmisión y al ruido contribuciones clave en la aparición de esos patrones (Deco et al., 2011). Se trata de resultados obtenidos con modelos. Indican qué mecanismos producen el fenómeno dentro de cada formulación, aunque dejan pendiente la demostración de que esos sean los mecanismos efectivamente operantes en el cerebro.

En modelos anatómicos suele aproximarse $d_{ij}=\ell_{ij}/v_{ij}$, donde $\ell$ representa una longitud y $v$ una velocidad de conducción. La longitud obtenida por un procedimiento de reconstrucción de tractos y la velocidad efectiva son cantidades de naturaleza distinta. Estimar una velocidad global supone que las diferencias no modeladas pueden resumirse en un único parámetro común. Si en cambio se deja variar libremente cada velocidad, la dimensión del problema inverso aumenta mucho.

El paso temporal limita la representación de los retardos. Redondear $d/h$ al entero más cercano altera el tiempo de transmisión. Cuando los retardos no son múltiplos del paso puede interpolarse la historia, comprobando la precisión de la interpolación. Un resultado que cambia con variaciones del retardo menores que esa discretización necesita revisión numérica antes de interpretarse como una diferencia fisiológica.

La duración de la simulación también importa. El tramo inicial conserva efectos de la historia impuesta. En estudios de actividad espontánea es habitual descartar un transitorio, cuya longitud debería justificarse comprobando que los observables se han estabilizado. Usar siempre la misma duración puede dejar residuos en modelos lentos o desechar datos sin necesidad en modelos rápidos.

405.4. La actividad simulada todavía no es el registro

Añadamos un filtro de observación

Usaremos $\tau_h=50$ ms como valor pedagógico. Este filtro no reproduce una respuesta hemodinámica humana calibrada. Su función es mostrar cómo una medición puede suavizar y retrasar los estados. Si el objetivo fuera la señal BOLD, haría falta una transformación neurovascular con su propia escala temporal, como el modelo del globo, que puede integrarse en un modelo hemodinámico de entrada, estado y salida dependiente de la actividad sináptica subyacente (Friston et al., 2000), además del muestreo y del ruido correspondientes a la adquisición.

El filtro atenúa las variaciones rápidas y desplaza la latencia del máximo. Por eso la comparación de los picos de $y$ entre regiones no identifica por sí sola los retardos de transmisión neuronal. Las diferencias pueden deberse a $\tau_i$, a $d_{ij}$, a $\tau_h$ o a combinaciones de ellos. Escribir la ecuación de observación hace visibles todas esas alternativas.

Los factores $a_i$ representan ganancias de medición. Si multiplicamos un estado por una constante y dividimos su ganancia observacional por la misma constante, la salida puede quedar intacta mediante una reparametrización compatible de la dinámica. Esta simetría de escala es un ejemplo de indeterminación estructural. El algoritmo de ajuste puede devolver un número concreto para cada factor sin que los datos los hayan fijado por separado.

La correlación añade otra pérdida de información sobre la escala. La correlación de Pearson entre dos señales no cambia si a cada una se le suma una constante o se la multiplica por un factor positivo. Ajustar únicamente correlaciones deja sin determinar las amplitudes absolutas y las ganancias observacionales. Un modelo que reproduce bien una matriz funcional puede fallar al predecir la amplitud de la respuesta a una perturbación. Ese fracaso informa sobre lo que la función objetivo había dejado libre.

En modelos de EEG o MEG, una matriz de proyección puede mezclar las fuentes antes de añadir ruido. Esa mezcla crea asociaciones entre sensores incluso cuando ciertas fuentes eran independientes. En otros contextos, el filtrado o la regresión de componentes compartidos modifican las dependencias entre señales. La observación y el procesamiento deben incluirse en la comparación entre datos simulados y adquiridos, aplicando a ambos los mismos procedimientos siempre que sea posible.

405.5. Calibración e identificabilidad

Calibrar consiste en seleccionar parámetros que aproximen las observaciones o resúmenes de ellas. La función objetivo puede incluir espectros, autocorrelaciones, matrices de conectividad o respuestas evocadas. Cada elección privilegia ciertas propiedades del sistema. Un error pequeño en una función concreta no demuestra que el modelo sea adecuado respecto de todo lo que es capaz de generar.

La identificabilidad pregunta si observaciones suficientes permiten distinguir los parámetros entre sí. Se habla de identificabilidad estructural cuando la pregunta se plantea para datos ideales, continuos y sin ruido. En nuestro modelo aparece una indeterminación exacta cuando la dinámica solo depende del producto $GC$. Reemplazar $C$ por $aC$ y $G$ por $G/a$ conserva ese producto, de manera que si se estima la escala de ambas cantidades sin una restricción externa el problema carece de solución única. Normalizar $C$ fija una convención a la que la interpretación de $G$ debe referirse siempre.

La identificabilidad práctica puede ser deficiente aunque no exista una simetría exacta. Un registro corto, ruidoso o muy filtrado puede volver casi indistinguibles dos constantes temporales. En ese caso los intervalos de incertidumbre amplios describen una limitación de la información disponible. La regularización y las distribuciones previas ayudan a estabilizar el ajuste, aunque introducen información o preferencias que conviene separar de lo aprendido a partir del registro.

Una herramienta útil son los perfiles de verosimilitud o de error. Se fija un parámetro, se reajustan los demás y se observa cuánto empeora el objetivo. Una región extensa de valores casi equivalentes revela incertidumbre. Este procedimiento permite detectar tanto indeterminaciones estructurales, que se manifiestan como relaciones funcionales entre parámetros, como indeterminaciones prácticas debidas a la cantidad y calidad limitadas de los datos. Sirve además para planificar experimentos (Raue et al., 2009). También puede examinarse la matriz de sensibilidad de las observaciones respecto de los parámetros. Si dos columnas se parecen mucho, las modificaciones correspondientes producen efectos difíciles de distinguir con el protocolo empleado.

Las pruebas de recuperación son especialmente informativas. Se generan observaciones con parámetros conocidos, se añade el ruido previsto y se ejecuta el procedimiento completo de estimación. Si en ese entorno controlado no se recuperan los valores o las predicciones con precisión suficiente, la interpretación de datos reales debe revisarse. Este paso se recomienda como parte de cualquier análisis basado en modelos, junto con su equivalente para la comparación entre modelos, que consiste en verificar si los datos simulados con un modelo son atribuidos a ese mismo modelo por el procedimiento de selección (Wilson y Collins, 2019). La prueba debe abarcar varias parametrizaciones y niveles de ruido, porque la recuperación de un único caso favorable aporta poca evidencia.

La selección entre mecanismos exige incluir alternativas. Si solo se ajusta una familia de modelos, cualquier desacuerdo tenderá a absorberse en sus parámetros. Un cambio en la observación que el modelo no representa puede aparecer, por ejemplo, como un cambio de acoplamiento. Comparar modelos con estructuras diferentes permite examinar esa posibilidad, aunque tampoco garantiza que el mecanismo verdadero esté entre los candidatos.

405.6. Una simulación con historia explícita

El programa siguiente integra el ejemplo por el método de Euler con un paso de 0.1 ms. Conserva el historial completo de estados y usa retardos que son múltiplos exactos del paso. Incluye el filtro de observación y verifica la convergencia comparando el resultado con el que se obtiene al dividir el paso por dos.

import numpy as np

C = np.array([[0,1,0],[.5,0,.5],[0,1,0]], dtype=float)
delay = np.array([[0,5,0],[5,0,10],[0,10,0]], dtype=float)
tau = np.array([10.,20.,30.])

def simular(dt):
    t = np.arange(round(400/dt)+1)*dt
    x = np.zeros((len(t),3))
    z = np.zeros_like(x)
    pasos = np.rint(delay/dt).astype(int)
    for k in range(len(t)-1):
        entrada = np.zeros(3)
        if 20 <= t[k] < 40:
            entrada[0] = 1.0
        acople = np.zeros(3)
        for i in range(3):
            for j in range(3):
                q = k - pasos[i,j]
                pasado = x[q,j] if q >= 0 else 0.0
                acople[i] += C[i,j]*pasado
        x[k+1] = x[k] + dt*(-x[k] + .6*acople + entrada)/tau
        z[k+1] = z[k] + dt*(x[k]-z[k])/50.0
    return t,x,z

t,x,z = simular(.1)
t2,x2,z2 = simular(.05)
print('Error máximo al dividir el paso:', np.max(abs(x-x2[::2])))
print('Picos de estados:', x.max(axis=0))
print('Picos de observación filtrada:', z.max(axis=0))

Para cada enlace, el programa consulta el estado de la región de origen en el instante retrasado, con la convención de índices de la matriz declarada. El filtro usa el estado al comienzo del paso, como corresponde al método de Euler explícito. Emplear por descuido el estado recién actualizado cambiaría el esquema y su error. Esa variante podría ser válida como otro integrador, siempre que se analice y se documente.

Al ejecutar el programa, la diferencia máxima entre las trayectorias calculadas con pasos de 0.1 y 0.05 ms es de aproximadamente $9\times10^{-4}$, pequeña frente a la amplitud de la respuesta. Los máximos de los estados son cercanos a 0.874, 0.141 y 0.051 y se alcanzan hacia los 40.0, 51.7 y 87.1 ms. Después del filtro, los máximos bajan a 0.266, 0.066 y 0.033 y se desplazan a 53.9, 87.4 y 132.7 ms. La observación añade así entre 14 y 46 ms de latencia según la región. Comparar los picos filtrados como si midieran tiempos de transmisión llevaría a sobrestimar los retardos.

Para estudiar el efecto del muestreo podemos tomar uno de cada cincuenta puntos del resultado con paso de 0.1 ms, lo que equivale a una observación cada 5 ms. La simulación interna sigue usando el paso fino, porque integración y adquisición son operaciones diferentes. Un registro poco frecuente puede perder transitorios aunque la trayectoria interna se haya calculado con precisión.

Si añadimos ruido observacional, conviene fijar una semilla para reproducir el ensayo y repetirlo con otras semillas para evaluar la variabilidad. Si añadimos ruido de proceso continuo, su discretización requiere el factor adecuado, proporcional a $\sqrt{dt}$ en el caso de incrementos de Wiener. Sumar el mismo ruido en cada paso con independencia de $dt$ altera la intensidad efectiva del proceso e invalida la comparación de convergencia.

405.7. Diseñar una perturbación que reduzca ambigüedad

Supongamos que dos parametrizaciones reproducen por igual la correlación basal. Una combina transmisión rápida con filtrado lento; la otra, transmisión lenta con filtrado rápido. Una estimulación breve acompañada de un registro con mejor resolución temporal podría separar sus predicciones. Si solo se observa la salida muy filtrada, la diferencia puede seguir siendo invisible. El diseño tiene que considerar la observación disponible además de los estados internos simulados.

La idea de estudiar el acoplamiento entre regiones mediante entradas diseñadas y conocidas tiene una formulación consolidada en el modelado causal dinámico, que identifica en un marco bayesiano los parámetros de influencia de las entradas, de acoplamiento entre estados y de modulación de ese acoplamiento a partir de respuestas observadas (Friston et al., 2003). En modelos de cerebro virtual, un trabajo basado en inferencia mediante simulaciones informó que la estimulación sistemática permitió estimar de manera fiable la degradación de un subconjunto de conexiones que sin estimulación no era identificable (Hashemi et al., 2024). Ese resultado procede de datos sintéticos generados por el propio modelo y muestra la utilidad potencial de perturbar, sin establecer todavía su rendimiento con registros reales.

Una evaluación prospectiva útil define primero la perturbación, calcula las predicciones de cada candidato y selecciona un contraste en el que difieran de manera robusta frente al ruido. Después se realizan las observaciones y se comparan con esas predicciones. Reajustar todos los parámetros a la luz del resultado puede ser necesario para mejorar el modelo, pero convierte esa fase en aprendizaje. La validación requerirá entonces una observación independiente.

La perturbación tiene límites de factibilidad. Un impulso matemático sobre una variable latente no siempre corresponde a una intervención física realizable. Un dispositivo puede afectar varias regiones o componentes, con una intensidad y una duración conocidas solo de manera imperfecta. Esa incertidumbre debe incorporarse al operador de entrada $B$ y al análisis de sensibilidad. La capacidad formal para controlar un estado no demuestra que exista una técnica capaz de hacerlo de ese modo.

También conviene comparar con referencias sencillas. Un modelo generativo detallado debería aportar algo relevante frente a una predicción lineal, a la simple persistencia temporal o a una regresión con variables clínicas, cuando ese sea el uso propuesto. La complejidad puede justificarse por su valor interpretativo, por su capacidad de extrapolación o por su rendimiento predictivo, siempre que esas ventajas se evalúen de forma concreta. Un ajuste visualmente atractivo no basta para establecer ninguna de ellas.

405.8. Un control analítico antes de introducir heterogeneidad

Parte de la implementación puede verificarse eliminando los retardos y asignando a las tres regiones una constante común de 20 ms. La matriz $C$ tiene valores propios 1, 0 y −1. Con $G=0.6$, los valores propios del operador $(-I+GC)/20$ son −0.02, −0.05 y −0.08 ms$^{-1}$, que corresponden a tiempos de decaimiento de 50, 20 y 12.5 ms.

Los vectores propios asociados tienen una interpretación directa. El modo más lento, proporcional a $(1,1,1)^{\mathsf T}$, es una desviación común de las tres regiones; la recurrencia positiva compensa parte de la fuga y prolonga su recuperación. El modo intermedio, proporcional a $(1,0,-1)^{\mathsf T}$, opone entre sí los dos extremos de la red y deja quieta la región central. El más rápido, proporcional a $(1,-1,1)^{\mathsf T}$, opone la región central a los extremos. Una trayectoria cualquiera combina estos modos según la entrada y el estado inicial, por lo que una única constante temporal efectiva medida sobre ella puede mezclar contribuciones de varios.

Esta versión permite comparar la simulación con la solución exacta $x(t)=e^{At}x(0)$ en ausencia de entrada. Con $x(0)=(1,0,0)^{\mathsf T}$ y un paso de 0.1 ms, la diferencia máxima entre el resultado de Euler y la solución exacta durante 400 ms es del orden de $9\times10^{-4}$. Si el código fallara en este caso, habría que corregirlo antes de introducir retardos y observación. El control tiene un alcance limitado y a la vez valioso, porque comprueba índices, signos y factores temporales en una situación que dispone de una referencia independiente.

Para el método de Euler explícito, un modo real negativo $\lambda$ se convierte en el multiplicador $1+h\lambda$. La estabilidad numérica exige que su módulo sea menor que uno. En el modo más rápido del ejemplo, esto impone $h<25$ ms. Cumplir esa desigualdad evita la divergencia numérica de ese modo, aunque no asegura una buena precisión. Un paso cercano al límite distorsionaría fuertemente la trayectoria. La convergencia se evalúa comparando los observables al reducir el paso.

La heterogeneidad posterior modifica los modos. Con constantes de 10, 20 y 30 ms, los valores propios ya no se obtienen dividiendo los de $C$ por una misma cifra. Hay que formar el operador completo $T^{-1}(-I+GC)$. En este caso sus valores propios son aproximadamente −0.115, −0.05 y −0.019 ms$^{-1}$, con tiempos de decaimiento cercanos a 8.7, 20 y 53.8 ms. El orden de las multiplicaciones importa porque $T$ y $C$ no conmutan. Este detalle se pierde con facilidad al trasladar a código regional una derivación hecha para el caso homogéneo.

405.9. Una región observada puede ocultar un modo completo

En la versión homogénea y sin retardos, el vector $v=(1,0,-1)^{\mathsf T}$ cumple $Cv=0$. Su evolución es un decaimiento con tasa −0.05 ms$^{-1}$. Si el único sensor mide la región central, $y=x_2$, ese modo produce una salida nula durante toda su trayectoria. Los extremos pueden desviarse en sentidos opuestos sin dejar huella en la observación elegida. En términos formales, la matriz de observabilidad construida con este sensor tiene rango 2 en lugar de 3.

Dos estados iniciales que difieren solo en una componente a lo largo de $v$ generan la misma serie en la región central. El problema no se resuelve registrando durante más tiempo ni reduciendo el ruido a cero, porque la información falta por la estructura misma de la observación. Añadir una medida de uno de los extremos puede romper la indeterminación, siempre que su calidad permita distinguir esa contribución.

El ejemplo conecta los modelos cerebrales con la teoría de la observabilidad. Una gran cantidad de puntos temporales no equivale a observar todas las direcciones relevantes del estado. Por eso la elección de sensores o modalidades debe considerar qué combinaciones quedan ocultas. En modelos más complejos, un modo puede ser débilmente visible en lugar de exactamente invisible. Nuestro propio ejemplo lo muestra al reintroducir las constantes heterogéneas. Con 10, 20 y 30 ms, el mismo sensor central vuelve formalmente observable el sistema, pero el menor valor singular de la matriz de observabilidad es cercano a $7\times10^{-4}$, frente a un valor máximo cercano a 1 cuando el tiempo se mide en milisegundos. La reconstrucción es posible en principio y muy frágil frente al ruido.

Los filtros imponen otros límites. Un modo rápido puede alcanzar el sensor y quedar muy atenuado por la observación y el muestreo. Recuperarlo exige amplificar componentes en las que quizá domina el ruido. La invertibilidad formal de un filtro no garantiza una reconstrucción estable. Antes de atribuir una trayectoria latente precisa, interesa estudiar su sensibilidad al ruido y a los parámetros de observación.

Una simulación ayuda a diseñar esa evaluación. Se generan estados con componentes conocidas, se proyectan a los sensores y se intenta reconstruirlos. Los errores se comparan por modo o por variable, además del promedio global. Un error pequeño concentrado en direcciones irrelevantes para la pregunta puede ser tolerable, mientras que una dirección mal reconstruida que determine la decisión puede invalidar el análisis.

405.10. Qué resúmenes se utilizan para ajustar

Supongamos que calibramos el modelo para reproducir una matriz de correlaciones. Dos parametrizaciones pueden igualarla y diferir en el espectro o en la respuesta a pulsos. Añadir esos resúmenes impone restricciones adicionales, aunque también modifica el problema de optimización. Conviene escogerlos en función de la hipótesis que se quiere evaluar en vez de acumular medidas hasta lograr una apariencia de exhaustividad.

Las escalas de error deben ser comparables. Una diferencia espectral calculada sobre cientos de frecuencias puede dominar a una diferencia de correlación calculada sobre unos pocos pares por el solo número de términos. Normalizar o ponderar cada componente cambia lo que se considera un buen ajuste. Los pesos requieren justificación y análisis de sensibilidad. El óptimo obtenido depende de esa función objetivo.

Cuando se conocen las incertidumbres de los resúmenes, el error puede ponderarse respecto de ellas, de modo que una medida muy ruidosa influya menos que una estimada con precisión. Estimar esas incertidumbres a partir de pocos datos introduce, a su vez, error. Una verosimilitud explícita facilita la interpretación del ajuste cuando es apropiada, mientras que una distancia entre resúmenes debe aclarar qué aproximación representa.

Un procedimiento basado en simulaciones puede explorar miles de combinaciones de parámetros y conservar aquellas que producen resúmenes cercanos a los observados. El conjunto aceptado describe la compatibilidad con la tolerancia y los criterios elegidos. No constituye automáticamente una distribución posterior exacta. Los métodos de computación bayesiana aproximada formalizan esta idea y evitan evaluar la verosimilitud, a costa de supuestos y aproximaciones cuyo impacto debe examinarse en cada aplicación (Sunnåker et al., 2013). La interpretación depende del muestreo de parámetros, de la función de distancia, de la tolerancia y de la distribución previa utilizada.

La evaluación fuera de muestra debe reservar información que el ajuste no haya utilizado, sea otro intervalo, otra sesión o una perturbación. Separar puntos vecinos de una misma serie con fuerte autocorrelación produce una prueba demasiado fácil. Roberts et al. (2017) revisaron estrategias de validación cruzada para datos con estructura temporal, espacial o jerárquica y proponen separar bloques que respeten esa dependencia, porque la partición aleatoria tiende a producir estimaciones optimistas del error. Si el objetivo es anticipar una condición nueva, la validación debe aproximarse a esa transferencia.

405.11. Recuperación bajo un mundo ligeramente equivocado

Una primera prueba genera y ajusta los datos con exactamente la misma familia de modelos. Sirve para detectar errores de implementación y para evaluar la identificabilidad dentro de ese mundo idealizado. Si el procedimiento falla allí, necesita revisión. Si funciona, todavía puede fallar cuando la realidad incorpora estructuras que el modelo no representa. Superar esta prueba es necesario, aunque no suficiente.

Podemos añadir una variable lenta al generador y ajustar después el modelo de una variable por región. También podemos generar retardos heterogéneos y ajustar una velocidad única, o introducir ruido compartido cuando el ajuste supone independencia. Estas alteraciones del mundo simulado muestran hacia qué parámetros se desplaza el error. Un aumento aparente del acoplamiento, por ejemplo, puede ser la compensación de una entrada omitida.

La comparación debe medir tanto la recuperación de parámetros como la calidad predictiva. Es posible que los parámetros pierdan su interpretación mecanística y que ciertas predicciones sigan siendo útiles. En ese caso el modelo puede servir como aproximación predictiva dentro del dominio evaluado. Lo que debe revisarse entonces es la afirmación sobre mecanismos, mientras que la utilidad del modelo puede mantenerse en otros aspectos.

Otra prueba introduce una intervención que modifica el proceso de observación sin alterar la dinámica neuronal. Si el procedimiento estima después un cambio de conectividad, queda identificada una vía de confusión. Una modalidad adicional, una calibración del instrumento o un modelo observacional más flexible pueden reducirla. La mejora debe comprobarse con una nueva evaluación, ya que añadir parámetros no resuelve por sí mismo la identificación.

Por último, el modelo detallado puede compararse con uno reducido bajo el mismo presupuesto de datos. Si el reducido predice igual de bien y ofrece una incertidumbre más estable, quizá resulte preferible para esa aplicación. Si el detallado anticipa correctamente una perturbación que el reducido no puede representar, esa ventaja justifica su mayor complejidad. La decisión se toma a partir de resultados concretos y del propósito de uso, conservando la arquitectura necesaria para explicarlos.

405.12. Elegir una entrada que distinga parámetros

Una vez formulado el modelo, cabe preguntar qué estímulo ayudaría a estimar un parámetro determinado. Consideremos primero la versión sin retardos, $\dot x=A(\theta)x+Bu$, con $\theta$ escalar. La sensibilidad del estado, $s=\partial x/\partial\theta$, satisface

siempre que la entrada y su matriz $B$ no dependan de $\theta$. Si la condición inicial es conocida e independiente del parámetro, el cálculo parte de $s(0)=0$. Esta ecuación puede integrarse junto con la trayectoria. Su solución indica cuánto cambiaría el estado ante una variación pequeña del parámetro en torno al valor utilizado.

Para la ganancia global de nuestra red, $\partial A/\partial G=T^{-1}C$. Como este término aparece multiplicado por $x(t)$, la sensibilidad depende de los modos que la entrada haya excitado. Un estímulo que apenas activa las rutas de acoplamiento puede aportar poca información sobre $G$ aunque produzca una señal local intensa. La amplitud visible de una respuesta y la información que contiene sobre un parámetro son cosas distintas.

La observación transforma esa sensibilidad. Con un sensor lineal instantáneo $y=Hx+\varepsilon$, la derivada de la media observada es $Hs$. Si el sensor anula la dirección en que cambia el estado, la sensibilidad útil desaparece. Cuando hay un filtro de observación, también debe propagarse su dependencia temporal. Diseñar el estímulo sobre los estados ocultos sin pasar por el sensor puede sobrestimar mucho su capacidad para distinguir modelos.

Supongamos, solo para ilustrar el cálculo, errores gaussianos independientes con varianza conocida $\sigma^2$ en los instantes de medición. La información de Fisher para el parámetro, que actúa como medida local de información, es

donde $\mu$ es la media predicha de la señal. El criterio favorece los diseños que producen cambios observables grandes en relación con el ruido. Si los errores están correlacionados, la suma simple debe sustituirse por una forma cuadrática que incorpore su covarianza. Muchos puntos vecinos correlacionados aportan menos información que el mismo número de mediciones independientes.

Con varios parámetros, las sensibilidades forman las columnas de una matriz. Dos columnas casi proporcionales indican que los cambios en los parámetros correspondientes producen efectos parecidos y pueden compensarse entre sí. Un experimento puede ser muy sensible a una combinación de parámetros y apenas distinguir sus componentes. En ese caso, aumentar la amplitud de la misma entrada mejora la precisión en la dirección ya visible sin resolver necesariamente la ambigüedad.

En nuestra arquitectura podemos comparar pulsos aplicados en una región extrema y en la región central, manteniendo controladas la duración y la energía. También podemos observar ambas regiones durante la recuperación. La comparación se resuelve calculando las sensibilidades después de aplicar el mismo filtro y el mismo muestreo previstos para la adquisición. Si un diseño separa mejor la ganancia global de una constante local, ofrece una razón concreta para preferirlo.

El criterio es local porque depende de los valores de parámetros supuestos inicialmente. Un diseño excelente en torno a un valor puede ser pobre en otra región plausible del espacio de parámetros. Por eso conviene evaluarlo sobre un conjunto de parámetros compatibles con datos previos. Puede buscarse un buen desempeño promedio o protegerse frente al caso menos informativo, explicitando el criterio adoptado, que refleja el propósito de la adquisición.

Ninguna ganancia de información matemática autoriza por sí misma una intervención sobre personas. La entrada del ejemplo pertenece a una simulación. En una aplicación experimental, las opciones disponibles y sus restricciones éticas y técnicas definen primero el conjunto de diseños admisibles. Dentro de ese conjunto, las ecuaciones ayudan a escoger mediciones que reduzcan una incertidumbre específica y a anticipar cuándo una nueva observación aportará poco a la pregunta mecanística.

La predicción que guíe el diseño debe quedar formulada antes de obtener el registro, indicando qué señal cambiará, en qué intervalo y con qué diferencia esperada entre los candidatos. Esa especificación previa vuelve evaluable la adquisición propuesta y permite interpretar también un resultado que mantenga la ambigüedad.

405.13. Qué debe quedar al terminar un modelo

Un modelo reproducible necesita conservar la matriz utilizada y sus transformaciones, las ecuaciones con sus parámetros, la historia inicial y las reglas de integración. También debe registrar el operador de observación y el procesamiento aplicado a la salida. Sin estos elementos, dos equipos pueden compartir el nombre de una plataforma y estar ejecutando, en realidad, experimentos diferentes.

La evaluación final debe distinguir la verificación numérica, la calibración y el contraste empírico. La verificación comprueba que el programa resuelve lo que se especificó. La calibración encuentra parámetros compatibles con un conjunto de datos. El contraste empírico pregunta si las predicciones se sostienen frente a observaciones que podrían haberlas contradicho. Estas fases pueden repetirse durante el desarrollo, pero la evidencia que produce cada una es de naturaleza distinta.

En nuestro ejemplo hemos verificado una arquitectura con estados regionales, transmisión retardada y observación filtrada. También hemos visto cómo cambian sus respuestas al modificar un componente sin confundirlo con otro. El capítulo siguiente estudiará medidas de coordinación entre esas trayectorias. Su interpretación seguirá dependiendo de la arquitectura que las generó y de la parte de ella que llega realmente a la señal observada.

Referencias

• Breakspear, M. (2017). Dynamic models of large-scale brain activity. Nature Neuroscience, 20(3), 340–352. doi:10.1038/nn.4497
• Deco, G., Jirsa, V. K. y McIntosh, A. R. (2011). Emerging concepts for the dynamical organization of resting-state activity in the brain. Nature Reviews Neuroscience, 12(1), 43–56. doi:10.1038/nrn2961
• Friston, K. J., Harrison, L. y Penny, W. (2003). Dynamic causal modelling. NeuroImage, 19(4), 1273–1302. doi:10.1016/S1053-8119(03)00202-7
• Friston, K. J., Mechelli, A., Turner, R. y Price, C. J. (2000). Nonlinear responses in fMRI: The Balloon model, Volterra kernels, and other hemodynamics. NeuroImage, 12(4), 466–477. doi:10.1006/nimg.2000.0630
• Ghosh, A., Rho, Y., McIntosh, A. R., et al. (2008). Noise during rest enables the exploration of the brain's dynamic repertoire. PLoS Computational Biology, 4(10), e1000196. doi:10.1371/journal.pcbi.1000196
• Haddad, W. M. y Chellaboina, V. (2004). Stability theory for nonnegative and compartmental dynamical systems with time delay. En Proceedings of the 2004 American Control Conference (vol. 2, pp. 1422–1427). IEEE. doi:10.23919/ACC.2004.1386775
• Hashemi, M., Ziaeemehr, A., Woodman, M., et al. (2024). Simulation-based inference on virtual brain models of disorders. Machine Learning: Science and Technology, 5(3), 035019. doi:10.1088/2632-2153/ad6230
• Honey, C. J., Sporns, O., Cammoun, L., et al. (2009). Predicting human resting-state functional connectivity from structural connectivity. Proceedings of the National Academy of Sciences, 106(6), 2035–2040. doi:10.1073/pnas.0811168106
• Kaczorek, T. (2009). Stability of positive continuous-time linear systems with delays. En 2009 European Control Conference (ECC) (pp. 1610–1613). IEEE. doi:10.23919/ECC.2009.7074637
• Raue, A., Kreutz, C., Maiwald, T., et al. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923–1929. doi:10.1093/bioinformatics/btp358
• Roberts, D. R., Bahn, V., Ciuti, S., et al. (2017). Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography, 40(8), 913–929. doi:10.1111/ecog.02881
• Sanz-Leon, P., Knock, S. A., Spiegler, A., et al. (2015). Mathematical framework for large-scale brain network modeling in The Virtual Brain. NeuroImage, 111, 385–430. doi:10.1016/j.neuroimage.2015.01.002
• Sanz-Leon, P., Knock, S. A., Woodman, M., et al. (2013). The Virtual Brain: A simulator of primate brain network dynamics. Frontiers in Neuroinformatics, 7, 10. doi:10.3389/fninf.2013.00010
• Sunnåker, M., Busetto, A. G., Numminen, E., et al. (2013). Approximate Bayesian computation. PLoS Computational Biology, 9(1), e1002803. doi:10.1371/journal.pcbi.1002803
• 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, R. C. y Collins, A. G. E. (2019). Ten simple rules for the computational modeling of behavioral data. eLife, 8, e49547. doi:10.7554/eLife.49547
• Zalesky, A., Fornito, A., Harding, I. H., et al. (2010). Whole-brain anatomical networks: Does the choice of nodes matter? NeuroImage, 50(3), 970–983. doi:10.1016/j.neuroimage.2009.12.027