How can conditional flow matching capture multimodal single-neuron electrophysiology responses?

Why one voltage trace per stimulus is the wrong target

The same repeated stimulus can elicit very different voltage responses from the same cell, visible in trial-to-trial variability in spike timing, firing rate and burstiness. Biophysically detailed multi-compartment models capture this variability by keeping an ensemble of deterministic 'hall-of-fame' parametrizations, each tuned by evolutionary optimization to match multimodal recordings. That ensemble is an imperfect substitute: each member is itself deterministic, so the spread reflects uncertainty in the model parametrization rather than genuine trial-to-trial variability, and obtaining it costs on the order of hundreds of thousands of CPU hours.

Existing machine-learning surrogates inherit the same limitation. Neural networks, physics-informed neural networks and neural operators map a stimulus to a single voltage response, and they are evaluated predominantly in the regular firing regime, leaving the highly nonlinear transition around firing threshold unclear. The paper's framing is that single-neuron electrophysiology should instead be posed as conditional generation of a membrane potential V given an applied current I, directly targeting p(V|I). That reframing is what makes multimodality in spike timing something the model must represent rather than average away.

How flow matching turns current into a voltage distribution

The model is a conditional generative model trained with flow matching, using a velocity field conditioned on the input current. Rather than regressing one trace, it transports samples toward the conditional voltage distribution, so a single conditioning input can yield an ensemble of responses. Conditioning on the full time-aligned stimulus waveform, not just a scalar amplitude, is what later allows the same architecture to handle more general stimulation protocols.

Integration choices were checked before the main results. Euler, midpoint and fourth-order Runge-Kutta were compared at matched numbers of function evaluations in terms of sample quality; all three converge to a common limit by 200 NFE, set by error in the learned velocity field rather than by the integrator. Midpoint converges fastest, reaching that common limit in 10 to 20 NFE, and the reported results use 20 midpoint steps (40 NFE) so they are not confounded by integration error, although fewer steps can still produce high-quality samples.

What the training data and comparison baseline are

Evaluation focuses on one parvalbumin-positive (PVALB) neuron and one vasoactive intestinal peptide (VIP) neuron from an inhibitory interneuron Patch-seq dataset of the human middle temporal gyrus. The two types were chosen because they exhibit qualitatively distinct firing regimes: PVALB neurons are fast-spiking interneurons that orchestrate high-frequency gamma oscillations (30-80 Hz), while VIP neurons operate on a slower timescale, which lets the authors test generalization across electrophysiological phenotypes.

Voltage response datasets were generated from a single neuron morphology and 60 HoF conductance parametrizations from multi-objective evolutionary optimization, of which 50 were used for training. For each model, 1,680 stimuli were drawn, giving 100,800 current-voltage pairs, with every stimulus paired with a single model so a conditioner is almost never observed more than once. Responses were simulated in NEURON at a 0.02 ms timestep over a 400 ms activation window, cropped to 45 ms before onset, the activation window and 115 ms after it, then low-pass filtered and decimated by a factor of three, giving 9,334 samples at delta t = 0.06 ms. The 50 training models contribute 84,000 traces, split 90/10 into 75,600 training and 8,400 validation examples. DC step amplitudes were drawn from a truncated skew-normal distribution on [-0.11, 0.28] nA, the range over which the HoF models are experimentally validated, with the mode placed in the peri-threshold window; the same distribution was used for the NOBLE training data, making the comparisons fair. Comparisons are against the HoF models, NOBLE as trained by the original authors, and real experimental human cortex recordings.

Generated ensembles match feature distributions at unseen currents

At unseen current amplitudes, the generated ensembles for both PVALB and VIP closely resemble those produced by the HoF models, capturing both response dynamics and their variability. Spike-triggered averages show the median and interquartile range of generated responses closely matching the HoFs, with experimental curves lying within the same range, and F-I curves show generated firing-rate distributions tracking the HoFs across stimulation amplitudes while experimental measurements align closely with the median response.

Because pointwise voltage error is sensitive to small temporal shifts, the paper also computes electrophysiological features per response. Generated distributions closely match the HoFs across features and amplitudes, with explained variance above 95% for several features across amplitudes. Performance is lower for inverse inter-spike interval (ISI) values. Separating agreement into location and spread, the normalized mean bias stays within about 0.25 sigma for every feature and amplitude; the largest discrepancy is in dispersion, with generated inter-spike intervals under-dispersed at 0.10 nA, consistent with the lower explained variance for inverse ISI.

Wasserstein-2 distances between reference and generated feature distributions remain small relative to the underlying variance. Averaged across all amplitudes for the PVALB neuron, distances are 25.3 microseconds for AP width, 0.595 mV for AHP depth, 1.14 mV for mean AP amplitude, 0.533 Hz for inverse ISI and 0.334 Hz for mean frequency. Because averaging across all amplitudes could hide regime-specific failure modes, the analysis is repeated across three regimes: threshold amplitudes, where at least one HoF model produces a steady-state response and at least one produces a firing response; low-firing amplitudes, with mean spike counts between 5 and 17; and high-firing amplitudes, with mean spike counts between 18 and 37. Performance remains consistent with no regime-specific failure mode, though the threshold regime is noticeably more challenging for several features.

Threshold and depolarization block separate the two approaches

The clearest distinction from NOBLE appears in the firing threshold regime. At 0.02 nA the HoF data contain both firing and non-firing responses, and the generative model captures both, while NOBLE produces no firing responses at all. Firing probability as a function of stimulation amplitude shows the generative model closely following the HoF distribution, whereas NOBLE shifts toward higher firing thresholds. The authors' explanation is that NOBLE minimizes a pointwise relative L4 objective over a multimodal distribution; near threshold, where firing and non-firing responses coexist at the same amplitude, that objective can prefer an intermediate response rather than either mode, suppressing spiking and shifting the predicted threshold toward higher currents. The firing threshold is precisely where the underlying dynamics sit near a bifurcation, which is why it discriminates between the two approaches more sharply than either the sub- or supra-threshold regime.

A second separation appears at high amplitudes, where the VIP HoF ensemble approaches depolarization block. Differences across HoF parametrizations place models on either side of this transition, with some continuing to fire while others enter a low-firing state, producing a multimodal conditional response distribution. At 0.28 nA the HoF spike-count distribution separates into distinct low- and high-firing modes; the generative model recovers both modes and their approximate probability mass, whereas NOBLE spreads mass across the range and resolves neither mode. The paper notes that this difference would be obscured by the mean response alone.

Matching feature distributions alone does not guarantee that temporal structure is preserved, since spike timing determines when information reaches downstream neurons. Spike-time densities across representative amplitudes for both PVALB and VIP show the generative model closely following the HoF spike-time densities across regimes, while NOBLE shows larger discrepancies and, as before, does not capture spikes at amplitudes within the threshold regime.

Time-varying stimuli need the matching training distribution

Conditioning on the full time-aligned stimulus waveform allows more general protocols, so a separate model was trained on a new dataset of time-varying current waveforms using the same model, optimizer and training hyperparameters as the DC-step model. Time-varying stimuli were built by summing a random number of sine, step, ramp and chirp basis components with independently sampled amplitude, onset and duration, with at least one ramp always included so no pure square pulse occurs in that dataset; evaluation on step currents is therefore a test of transfer.

Performance on square pulse currents remains comparable across threshold, low-firing and high-firing regimes with no new regime-specific failure mode, so expanding the training distribution does not substantially degrade reproduction of DC-step response distributions. On unseen time-varying stimuli, the DC-trained model generalizes to some non-DC inputs by producing spikes in the correct regions, but blur in the generated ensembles indicates poor ISI consistency. The paper's explanation is architectural: because the same kernels apply everywhere, any input that locally resembles a DC step within the kernel's receptive field is indistinguishable from a step, so the model stitches locally step-like predictions across the trace, giving good local behavior but poor global consistency. Performance on the ISI distribution suffers more once input frequency exceeds 40 Hz, at which point kernel weights learned for locally step-like behavior are no longer suitable. Training on time-varying currents requires no architectural change and lets the model capture time-varying voltage dynamics and generalize across stimulus families.

What to trust and what remains open

The evidence for this framework is bounded in ways that matter for reading the results. Evaluation rests on two human cortical interneuron types, one PVALB and one VIP neuron, each generated from a single neuron morphology and a set of 60 HoF conductance parametrizations, with 50 used for training. The reference distribution is the HoF ensemble, which the paper itself describes as reflecting uncertainty in model parametrization rather than genuine trial-to-trial variability, so agreement with HoFs is agreement with that surrogate target, not direct validation against biological variability. Agreement with experimental recordings is reported for feature distributions, spike-triggered averages and F-I curves, not as a full distributional match.

The reported results use 20 midpoint steps (40 NFE) to avoid confounding by integration error, and the paper notes that fewer steps can still produce high-quality samples, so the sampling budget is a design choice rather than a hard requirement. The DC-trained model's failure on time-varying inputs above roughly 40 Hz is a property of that training distribution, not of the architecture, since the time-varying-trained model uses the same architecture. Extending the framework from single neurons to interacting populations, and using distillation or few-step generation to approach the speed of deterministic surrogates such as NOBLE, are stated as future directions rather than demonstrated results.

The decision-relevant split is regime, not average error

The practical priority this paper establishes is that the choice between a deterministic surrogate and a conditional generative one should be made on the regimes where the conditional response distribution is multimodal, not on average error across amplitudes. Aggregated Wasserstein-2 distances stay small for both approaches, and the paper explicitly notes that averaging across all amplitudes could hide regime-specific failure modes. The failures that matter cluster at the firing threshold, where firing and non-firing responses coexist at the same stimulus amplitude, and near depolarization block, where the VIP ensemble splits into low- and high-firing modes.

That is why the threshold regime is the sharpest test in this study: it is where the underlying dynamics sit near a bifurcation, and it is where NOBLE's deterministic relative L4 objective suppresses spiking and shifts the predicted firing threshold toward higher currents, and where it blurs the gap between modes at depolarization block. A reader comparing surrogates for single-neuron electrophysiology should therefore ask whether the model reproduces the modes and their approximate probability mass at these specific amplitudes, and whether the training distribution covers the stimulus family being queried, since the DC-trained model's ISI consistency degrades above roughly 40 Hz input frequency while the time-varying-trained model handles those inputs without architectural change.

Sources used in this answer