Track 4 Workshop — FAIR Brain Data Science Bootcamp
May 12, 2026
virtual-twin.github.io/Brain-Simulation-Workshop
To predict outcome, model the whole network, not a catalogue of regions.
A model has to reproduce the resting repertoire, not only stimulus responses.
Direct comparison with data, no translation gap.
| Ingredient | Symbol | Source |
|---|---|---|
| Connectome (weights, delays) | \(A_{ij},\ \tau_{ij}\) | DTI / tractography |
| Local dynamics | \(f\) | Neural mass / mean-field model |
| Coupling | \(g\) | Network propagation, with delays |
| Observation | \(h\) | BOLD / EEG / MEG forward model |
The next slides unpack each building block.
Diffusion MRI → tractography → parcellation → connectivity matrix. The output is the structural backbone \(A_{ij}\) (weights) and \(\tau_{ij}\) (delays).
A mass \(m\) on a spring with stiffness \(k\).
Perpetual oscillation at \(\omega = \sqrt{k/m}\).
State Equations: \[ \dot{x} = v \] \[ \dot{v} = - \frac{k}{m}*x \]
State Variables
| Variable | Initial Value | Unit | Description |
|---|---|---|---|
| \(x\) | 2.0 | \(\mathrm{m}\) | Displacement from equilibrium |
| \(v\) | 0.0 | \(\frac{\mathrm{m}}{\mathrm{s}}\) | Velocity |
Parameters
| Parameter | Value | Unit | Description |
|---|---|---|---|
| \(k\) | 0.0001 | \(\frac{\mathrm{N}}{\mathrm{m}}\) | Spring stiffness constant |
| \(m\) | 1.0 | \(\mathrm{kg}\) | Mass of the oscillating body |
Same spring-mass system as before.
What changes?
What do we see?
Same spring-mass system. Only the representation changes.
What changes?
What do we see?
Same equations and same initial condition. Only the mass changes.
What changes?
What do we see?
Now we modify the toy model toward a more realistic oscillator.
What changes?
\[ \dot{v} = -\frac{k}{m}x - \frac{c}{m}v + g_{\mathrm{eff}} \]
What do we see?
Spring Oscillator
Damped Spring 
JansenRit1995 Lumped-parameter cortical column model with three neural populations: pyramidal cells, excitatory interneurons, and inhibitory interneurons. Each population is modeled by a PSP block (2nd-order linear filter) and a sigmoidal potential-to-firing-rate transformation.
Derived Variables \[ v_{pyr} = y_{1} - y_{2} \]
State Equations \[ \dot{y_{0}} = y_{3} \] \[ \dot{y_{3}} = - y_{0}*a^{2} - 2*a*y_{3} + A*a*\operatorname{Sigm}{\left(y_{1} - y_{2} \right)} \] \[ \dot{y_{1}} = y_{4} \] \[ \dot{y_{4}} = - y_{1}*a^{2} - 2*a*y_{4} + A*a*\left(p + C_{2}*\operatorname{Sigm}{\left(C_{1}*y_{0} \right)} + K*c_{intercolumn}\right) \] \[ \dot{y_{2}} = y_{5} \] \[ \dot{y_{5}} = - y_{2}*b^{2} - 2*b*y_{5} + B*C_{4}*b*\operatorname{Sigm}{\left(C_{3}*y_{0} \right)} \]
with \[ \operatorname{Sigm}{\left(v \right)} = \frac{2*e_{0}}{1 + e^{r*\left(v_{0} - v\right)}} \]
State Variables
| Variable | Initial Value | Unit | Description |
|---|---|---|---|
| \(y_{0}\) | 0.1 | \(\mathrm{mV}\) | |
| \(y_{3}\) | 0.1 | \(\frac{\mathrm{mV}}{\mathrm{s}}\) | |
| \(y_{1}\) | 0.1 | \(\mathrm{mV}\) | |
| \(y_{4}\) | 0.1 | \(\frac{\mathrm{mV}}{\mathrm{s}}\) | |
| \(y_{2}\) | 0.1 | \(\mathrm{mV}\) | |
| \(y_{5}\) | 0.1 | \(\frac{\mathrm{mV}}{\mathrm{s}}\) |
Parameters
| Parameter | Value | Unit | Description |
|---|---|---|---|
| \(A\) | 3.25 | \(\mathrm{mV}\) | Maximum amplitude of excitatory PSP (EPSP) |
| \(B\) | 22 | \(\mathrm{mV}\) | Maximum amplitude of inhibitory PSP (IPSP) |
| \(C\) | 135 | — | Lumped connectivity constant (C = C1) |
| \(a\) | 100 | \(\mathrm{s}^{-1}\) | Reciprocal of EPSP time constant |
| \(b\) | 50 | \(\mathrm{s}^{-1}\) | Reciprocal of IPSP time constant |
| \(v_{0}\) | 6 | \(\mathrm{mV}\) | PSP for 50% firing rate (sigmoid midpoint) |
| \(e_{0}\) | 2.5 | \(\mathrm{s}^{-1}\) | Half of maximum firing rate |
| \(r\) | 0.56 | \(\mathrm{mV}^{-1}\) | Steepness of sigmoid transformation |
| \(p\) | 220 | \(\mathrm{s}^{-1}\) | Mean external input (afferent pulse density) |
| \(K\) | 0 | — | Per-node coupling strength scaling (K[i] scales coupling received by node i) |
Derived Parameters
\[ C_{1} = C \] \[ C_{2} = 0.8*C \] \[ C_{3} = 0.25*C \] \[ C_{4} = 0.25*C \]
Kuramoto (1975)
Phase oscillator model. Each node evolves with its own natural frequency \(\omega\); coupling \(I\) is added by the network.
State Equations \[\dot{\theta} = \omega\]
State Variables
| Variable | Initial value | Unit |
|---|---|---|
| \(\theta\) | \(\sim\mathcal{U}(0,2\pi)\) | — |
Derived Variables \[activity = \sin(\theta)\]
Parameters
| Parameter | Value | Unit | Description |
|---|---|---|---|
| \(\omega\) | 1.0 | — | Natural frequency |
Simulation settings
| Duration | 100 s |
| Step size | 0.01 s |
| Integrator | Heun |
| Nodes | 3 (V1, V2, MT) |
| Coupling \(K\) | 0.0 (uncoupled) |
Kuramoto (1975) — same model as before, now with coupling enabled.
State Equations \[\dot{\theta} = I + \omega\]
Coupling: \(I\)
\[ I = \frac{K \sum_{j=0}^{-1 + N} - {g}_{i,j} \sin{\left(\theta_{i} - \theta_{j} \right)}}{3} \]
Receives \(\theta\) from connected regions.
Pre-synaptic: \(c_{\text{pre}} = \sin(\theta_j - \theta_i)\)
Post-synaptic: \(c_{\text{post}} = \dfrac{K \cdot gx}{3}\)
Coupling parameters
| Parameter | Value | Description |
|---|---|---|
| \(K\) | 0.4 | Global coupling strength |
Network
| Edge | Weight |
|---|---|
| V1 → V2 | 0.7 |
| V2 → MT | 0.4 |
| V1 → MT | 0.2 |
Simulation settings
| Duration | 100 s |
| Step size | 0.01 s |
| Integrator | Heun |
\[dS_i = \left[f_d(S_i, \theta^d, C_i, I_i) \right]dt + g(S_i, \theta^g)\, dW_i\]
\[ C_i = f_c^{\text{post}}\Big(\sum_j A_{ij}\, f_c^{\text{pre}}(S_i, S_j(t-\tau_{ij}), \theta^c), S_i, \theta^c\Big) \]
Coupling
Time series \(x(t)\) at fixed \((a, x_0)\)
Phase plane flow geometry at fixed \(a\)
Bifurcation diagram regimes as \(a\) varies
\[ \dot x = a\,x \quad\Rightarrow\quad x(t) = x_0\,e^{a t} \]
Solution: \[ x(t) = x_0\,e^{a t} \]
Equilibrium: \[ \dot{x}(t) = ax(t) = 0? \] \[ \quad\Rightarrow\quad x_{eq} = 0 \]
\[ \dot z = (a + i\omega)z - |z|^2z, \qquad r_{\mathrm{cycle}} = \sqrt a \]
fixed point loses stability \(\Rightarrow\) stable oscillation appears
\[ F(x, a)=0, \qquad (x,a)=(x(s),a(s)) \]
TVBO + tvboptim simulation of Generic2dOscillator on the Desikan-Killiany network. Nodes are colored by their instantaneous \(V(t)\) and projected onto axial (top) and sagittal views; the regime map places each node at its current \((I_{\mathrm{eff}}, V)\) on the bifurcation diagram.
Inference: which \(\theta\) are consistent with the data?
What do we need to do?
Find global coupling \(G\) and recurrence \(w\) so simulated FC matches empirical FC in a Reduced Wong-Wang model on the Desikan-Killiany parcellation.
Optimization over a parameter landscape:
\(\mathcal{L}(\theta) = \mathrm{RMSE}(\hat{FC}(\theta),\, FC_{\text{obs}})\).
jax.vmap and jax.pmap under the hood
In more than 2–3 dimensions, random search beats grid for the same budget.
Why? Real loss surfaces have low effective dimensionality, most parameters barely matter. Grid spends its budget evenly across all axes; random samples every axis at every trial (Bergstra and Bengio 2012).
problem = {"num_vars": 3,
"names": ["G", "w", "sigma"],
"bounds": [[0.001, 0.7], [0.001, 0.7], [0.001, 0.01]]}
samples = sobol_sample.sample(problem, 256) # Saltelli design
state.coupling.G = DataAxis(samples[:, 0])
state.dynamics.w = DataAxis(samples[:, 1])
state.noise.sigma = DataAxis(samples[:, 2])
losses = ParallelExecution(loss, Space(state, mode="zip")).run()
Si = sobol_analyze.analyze(problem, np.array(losses))
Take-away: the FC loss lives on an effective 2D subspace.
DataAxis evaluation
Everything so far - grid, random / Sobol, Genetic algorithms - searches by sampling the simulator.
But a brain network model with regional parameters has \(d \sim 10^2\)–\(10^4\) knobs (e.g. one \(w_i\) per Desikan-Killiany region, plus per-region noise, plus …). No sampling budget reaches that regime.
Run forward once, run backward once. Reverse-mode AD sweeps the computation graph in reverse, accumulating \(\partial\mathcal{L}/\partial\theta_i\) for every \(i\) at the same time → \(\text{cost}(\nabla_\theta \mathcal{L}) \approx 3\text{–}30 \times \text{cost}(\mathcal{L})\).
Gradient cost in tvboptim, normalized by one simulation. Reverse-mode (blue) is flat at \(\sim 20\times\) across \(d \in [1, 10^4]\), independent of dimension.
In JAX: jax.grad(loss). One line, gradients with respect to every knob — what made per-region fits tractable.
\[ \theta_{t+1} \;=\; \theta_t - \eta \,\nabla_\theta\, \mathcal{L}(\theta_t) \]
A point estimate \(\theta^\star\) tells you one good fit. A posterior tells you all fits consistent with the data.
\[ \underbrace{p(\theta \mid y)}_{\text{posterior}} \;\propto\; \underbrace{p(y \mid \theta)}_{\text{likelihood}} \;\cdot\; \underbrace{p(\theta)}_{\text{prior}} \]
A prior \(p(\theta)\) is a probabilistic statement about parameters before you see this dataset. Three regimes:
Flat / uninformative: Uniform over a physiological range. “I know the bounds, nothing else.” Posterior is dominated by data.
Weakly informative: Normal around a literature value with generous spread. Regularizes, prevents pathological fits.Same simulator, same data, different priors → different posteriors. The prior is where domain knowledge enters the inference, cleanly and auditably.
Same model, same wide priors (\(G, w \sim \mathcal{U}(0.001, 0.7)\)). Different inference algorithm.
HMC (left) and SVI (right) on the same model. Background: grid-scan loss landscape. Red points: posterior samples; red contours: KDE; histograms: marginals. HMC bends along the valley; the SVI Gaussian aligns with it locally but cannot curve.
HMC / NUTS — gold standard.
SVI — fast approximation.
Not a bug in automatic differentiation, the dynamics are telling you the loss is non-smooth at fine scales.
Mitigations
→ each step uses the output of the previous as a warm start.
Original
Specification
Strongly informative: derived from independent measurements: anatomy, tau / amyloid PET maps, receptor density, EEG power. “This patient’s region has elevated tau” → prior on local excitation per region.- id: 6
label: "VEP: different columns, with delay
(Fig. 11)"
dynamics:
name: JansenRit1995_Delayed_VEP
parameters:
A: {value: 3.25, unit: mV}
B: {value: 22, unit: mV}
C: {value: 135}
a: {value: 100, unit: s^-1}
b: {value: 50, unit: s^-1}
v0: {value: 6, unit: mV}
e0: {value: 2.5, unit: s^-1}
r: {value: 0.56, unit: mV^-1}
p:
value: 220
unit: s^-1
shape: "(n_nodes,)"
distribution: *dist_p_uniform
a_d: {value: 30, unit: s^-1}
K: {value: 0, shape: "(n_nodes,)"}
derived_parameters: *jr_derived
functions: *jr_functions
coupling_terms: *jr_coupling
state_variables:
y0: {equation: {rhs: y3}}
y3:
equation:
rhs: >-
A*a*Sigm(y1 - y2)
- 2*a*y3 - a**2*y0
y1: {equation: {rhs: y4}}
y4:
equation:
rhs: >-
A*a*(p + P
+ C2*Sigm(C1*y0) + z0)
- 2*a*y4 - a**2*y1
y2: {equation: {rhs: y5}}
y5:
equation:
rhs: >-
B*b*(C4*Sigm(C3*y0))
- 2*b*y5 - b**2*y2
z0: {equation: {rhs: z1}}
z1:
equation:
rhs: >-
A*a_d*K*c_intercolumn
- 2*a_d*z1 - a_d**2*z0
derived_variables:
v_pyr: {equation: {rhs: y1 - y2}}
output: [v_pyr]
network:
number_of_regions: 2
number_of_nodes: 2
conduction_speed: {value: 0}
nodes:
- {id: 0}
- id: 1
parameters:
B: {value: 17.6}
C: {value: 108}
edges:
- source: 0
target: 1
directed: true
- source: 1
target: 0
directed: true
coupling: *sigmoidal_coupling
events: *flash_events
integration: *integration_vep
execution:
random_seed: 42
explorations:
VEP_delayed_fig11:
mode: zip
n_trials: 40
average: trials
parameters:
K:
element_domains:
- element: 0
explored_values: [300, 100]
- element: 1
explored_values: [1000, 4000]Reproduced
bsplot
Christoph Huettl
Rico Schmitt
Lia Domide
Romina Baila
Jil Meier
Halgurd Taher
Konstantin Buelau
Dionysios Perdikis
Michael Schirner
Leon Stefanovski
Claude Code
Petra Ritter
Marius Pille
marius.pille [at] bih-charite.de
Leon Martin
leon.martin [at] bih-charite.de
Leon Stefanovski
leon.stefanovski [at] charite.de
 |
 |
 |
|---|
