Simulation plays a central role in scientific discovery. In many applications, the bottleneck is no longer running a simulator; it is choosing among large families of plausible simulators, each corresponding to different forward models/hypotheses consistent with observations. Over large model families, classical Bayesian workflows for model selection are impractical. Furthermore, amortized model selection methods typically hard-code a fixed model prior or complexity penalty at training time, requiring users to commit to a particular parsimony assumption before seeing the data. We introduce PRISM, a simulation-based encoder-decoder that infers a joint posterior over both discrete model structures and associated continuous parameters, while enabling test-time control of model complexity via a tunable model prior that the network is conditioned on. We show that PRISM scales to families with combinatorially many (up to billions) of model instantiations on a synthetic symbolic regression task. As a scientific application, we evaluate PRISM on biophysical modeling for diffusion MRI data, showing the ability to perform model selection across several multi-compartment models, on both synthetic and in vivo neuroimaging data.

Objective.Brain-computer interfaces (BCIs) suffer from accuracy degradation as neural signals drift over time and vary across users, requiring frequent recalibration that limits practical deployment. Our goal is to develop a framework that eliminates the need for separate calibration phases by enabling continual, real-time model adaptation to new users and changing signal characteristics.Approach.We propose EDAPT, a task- and model-agnostic framework for continual online learning. EDAPT first establishes a robust baseline decoder through population-level pretraining (PRE) on data from multiple users. It then personalizes this model during deployment using supervised continual finetuning on a trial-by-trial basis. Due to its modular design, EDAPT can be composed with unsupervised domain adaptation (UDA) techniques to further address distribution shifts.Main results.We validate EDAPT across nine datasets, three BCI paradigms, and four deep learning architectures. EDAPT consistently improves decoding accuracy over static models for nearly all subjects and datasets, raising mean balanced accuracy from 0.80 to 0.87 on representative datasets (table3). Ablation studies confirm that the combination of PRE and online finetuning is the primary driver of this performance gain, with further improvements on some datasets when using UDA techniques. We demonstrate real-time feasibility of the framework, with adaptation latencies under 200 ms on consumer-grade hardware. Our scaling analysis further reveals that decoding accuracy is primarily determined by the total pretraining data budget, rather than its specific allocation between subjects and trials.Significance.These findings demonstrate that continual online learning is a practical and effective strategy for creating high-performance, user-adaptive BCIs. By systematically addressing the bottleneck of model recalibration, EDAPT reduces a major barrier to the widespread adoption of BCI technology and helps advance neurotechnology toward robust, user-friendly, real-world applications.

Neural activity fluctuates over a wide range of timescales within and across brain areas. Experimental observations suggest that diverse neural timescales reflect information in dynamic environments. However, how timescales are defined and measured from brain recordings vary across the literature. Moreover, these observations do not specify the mechanisms underlying timescale variations, nor whether specific timescales are necessary for neural computation and brain function. Here, we synthesize three directions where computational approaches can distill the broad set of empirical observations into quantitative and testable theories: We review (i) how different data analysis methods quantify timescales across distinct behavioral states and recording modalities, (ii) how biophysical models provide mechanistic explanations for the emergence of diverse timescales, and (iii) how task-performing networks and machine learning models uncover the functional relevance of neural timescales. This integrative computational perspective thus complements experimental investigations, providing a holistic view on how neural timescales reflect the relationship between brain structure, dynamics, and behavior.

Automated methods for discovering mechanistic simulator models from observational data offer a promising path toward accelerating scientific progress. Such methods often take the form of agentic-style iterative workflows that repeatedly propose and revise candidate models by imitating human discovery processes. However, existing LLM-based approaches typically implement such workflows via hand-crafted heuristic procedures, without an explicit probabilistic formulation. We recast model discovery as probabilistic inference, i.e., as sampling from an unknown distribution over mechanistic models capable of explaining the data. This perspective provides a unified way to reason about model proposal, refinement, and selection within a single inference framework. As a concrete instantiation of this view, we introduce ModelSMC, an algorithm based on Sequential Monte Carlo sampling that represents candidate models as particles which are iteratively proposed and refined by an LLM, and weighted using likelihood-based criteria. Experiments on real-world scientific systems illustrate that this formulation discovers models with interpretable mechanisms and improves posterior predictive checks. More broadly, this perspective provides a probabilistic lens for understanding and developing LLM-based approaches to model discovery.

While age-specific fertility rates (ASFRs) provide the most extensive record of reproductive change, their aggregate nature masks the underlying behavioral mechanisms that ultimately drive fertility trends. To recover these mechanisms, we develop a likelihood-free Bayesian framework that couples an individual-level model of the reproductive process with Sequential Neural Posterior Estimation (SNPE). This allows us to infer eight behavioral and biological parameters from just two aggregate series: ASFRs and the age-profile of planned versus unplanned births. Applied to U.S. National Survey of Family Growth cohorts and to Demographic and Health Survey cohorts from Colombia, the Dominican Republic, and Peru, the method reproduces observed fertility schedules and, critically, predicts out-of-sample micro-level distributions of age at first sex, inter-birth intervals, and family-size ideals, none of which inform the estimation step. Because the fitted model yields complete synthetic life histories, it enables behaviorally explicit population forecasts and supports the construction of demographic digital twins.

The temporal order of a sequence of events has been thought to be reflected in the ordered firing of neurons at different phases of theta oscillations. Here we assess this by measuring single neuron activity (1,420 neurons) and local field potentials (921 channels) in the medial temporal lobe of 16 patients with epilepsy performing a working-memory task for temporal order. During memory maintenance, we observe theta oscillations, preferential firing of single neurons to theta phase and a close relationship between phase of firing and item position. However, the firing order did not match item order. Training recurrent neural networks to perform an analogous task, we also show the generation of theta oscillations, theta phase-dependent firing related to item position and, again, no match between firing and item order. Rather, our results suggest a mechanistic link between phase order, stimulus timing and oscillation frequency. In both biological and artificial neural networks, we provide evidence supporting the role of phase of firing in working-memory processing.

Mergers of binary neutron stars emit signals in both the gravitational-wave (GW) and electromagnetic spectra. Famously, the 2017 multi-messenger observation of GW170817 (refs. 1,2) led to scientific discoveries across cosmology3, nuclear physics4-6 and gravity7. Central to these results were the sky localization and distance obtained from the GW data, which, in the case of GW170817, helped to identify the associated electromagnetic transient, AT 2017gfo (ref. 8), 11 h after the GW signal. Fast analysis of GW data is critical for directing time-sensitive electromagnetic observations. However, owing to challenges arising from the length and complexity of signals, it is often necessary to make approximations that sacrifice accuracy. Here we present a machine-learning framework that performs complete binary neutron star inference in just 1 s without making any such approximations. Our approach enhances multi-messenger observations by providing: (1) accurate localization even before the merger; (2) improved localization precision by around 30% compared to approximate low-latency methods; and (3) detailed information on luminosity distance, inclination and masses, which can be used to prioritize expensive telescope time. Additionally, the flexibility and reduced cost of our method open new opportunities for equation-of-state studies. Finally, we demonstrate that our method scales to long signals, up to an hour in length, thus serving as a blueprint for data analysis for next-generation ground- and space-based detectors.

Simulating aging in 3D brain MRI scans can reveal disease progression patterns in neurological disorders such as Alzheimer's disease. Current deep learning-based generative models typically approach this problem by predicting future scans from a single observed scan. We investigate modeling brain aging via linear models in the latent space of convolutional autoencoders (MRExtrap). Our approach, MRExtrap, is based on our observation that autoencoders trained on brain MRIs create latent spaces where aging trajectories appear approximately linear. We train autoencoders on brain MRIs to create latent spaces, and investigate how these latent spaces allow predicting future MRIs through linear extrapolation based on age, using an estimated latent progression rate . For single-scan prediction, we propose using population-averaged and subject-specific priors on linear progression rates. We also demonstrate that predictions in the presence of additional scans can be flexibly updated using Bayesian posterior sampling, providing a mechanism for subject-specific refinement. On the ADNI dataset, MRExtrap predicts aging patterns accurately and beats a GAN-based baseline for single-volume prediction of brain aging. We also demonstrate and analyze multi-scan conditioning to incorporate subject-specific progression rates. Finally, we show that the latent progression rates in MRExtrap's linear framework correlate with disease and age-based aging patterns from previously studied structural atrophy rates. MRExtrap offers a simple and robust method for the age-based generation of 3D brain MRIs, particularly valuable in scenarios with multiple longitudinal observations.

Amortized simulation-based inference (SBI) methods train neural networks on simulated data to perform Bayesian inference. While this strategy avoids the need for tractable likelihoods, it often requires a large number of simulations and has been challenging to scale to time series data. Scientific simulators frequently emulate real-world dynamics through thousands of single-state transitions over time. We propose an SBI approach that can exploit such Markovian simulators by locally identifying parameters consistent with individual state transitions. We then compose these local results to obtain a posterior over parameters that align with the entire time series observation. We focus on applying this approach to neural posterior score estimation but also show how it can be applied, e.g., to neural likelihood (ratio) estimation. We demonstrate that our approach is more simulation-efficient than directly estimating the global posterior on several synthetic benchmark tasks and simulators used in ecology and epidemiology. Finally, we validate scalability and simulation efficiency of our approach by applying it to a high-dimensional Kolmogorov flow simulator with around one million data dimensions.

Biophysical neuron models provide insights into cellular mechanisms underlying neural computations. A central challenge has been to identify parameters of detailed biophysical models such that they match physiological measurements or perform computational tasks. Here we describe a framework for simulating biophysical models in neuroscience-JAXLEY-which addresses this challenge. By making use of automatic differentiation and GPU acceleration, JAXLEY enables optimizing large-scale biophysical models with gradient descent. JAXLEY can learn biophysical neuron models to match voltage or two-photon calcium recordings, sometimes orders of magnitude more efficiently than previous methods. JAXLEY also makes it possible to train biophysical neuron models to perform computational tasks. We train a recurrent neural network to perform working memory tasks, and a network of morphologically detailed neurons with 100,000 parameters to solve a computer vision task. JAXLEY improves the ability to build large-scale data- or task-constrained biophysical models, creating opportunities for investigating the mechanisms underlying neural computations across multiple scales.

Synaptic plasticity is the basis of learning and memory, but the link between synaptic changes and neural function remains elusive. Here, we used automated search algorithms to obtain thousands of strikingly diverse quadruplets of excitatory(E)-to-E, E-to-inhibitory(I), IE, and II plasticity rules, cooperating to stabilize recurrent spiking networks. Despite the fact that quadruplets were selected for homeostasis, more than 90% of them performed well in simple and more difficult memory tasks such as novelty detection, contextual novelty and sequence replay. Co-activity was crucial, i.e., most rules failed in isolation. Our purely local, unsupervised plasticity rules could also help solve computer games such as pong. Our work showcases automated discovery augmenting human intuition to find en masse solutions for high dimensional problems. ### Competing Interest Statement The authors have declared no competing interest.

Extracting the relationship between high-dimensional neural recordings and complex behavior is a ubiquitous problem in neuroscience. Encoding and decoding models target the conditional distribution of neural activity given behavior and vice versa, while dimensionality reduction techniques extract low-dimensional representations thereof. Variational autoencoders (VAEs) are flexible tools for inferring such low-dimensional embeddings but struggle to accurately model arbitrary conditional distributions such as those arising in neural encoding and decoding, let alone simultaneously. Here, we present a VAE-based approach for calculating such conditional distributions. We first validate our approach on a task with known ground truth. Next, we retrieve conditional distributions over masked body parts of walking flies. Finally, we decode motor trajectories from neural activity in a monkey-reach task and query the same VAE for the encoding distribution. Our approach unifies dimensionality reduction and learning conditional distributions, allowing the scaling of common analyses in neuroscience to today's high-dimensional multi-modal datasets.

Binary black holes (BBHs) in eccentric orbits produce distinct modulations in the emitted gravitational waves (GWs). The measurement of orbital eccentricity can provide robust evidence for dynamical binary formation channels. We analyze 57 GW events from the first, second, and third observing runs of the Laser Interferometer Gravitational-Wave Observatory (LIGO)-Virgo-Kamioka Gravitational Wave Detector (KAGRA) (LVK) Collaboration using a multipolar aligned-spin inspiral-merger-ringdown waveform model with two eccentric parameters: eccentricity and relativistic anomaly (assuming a quasicircular merger-ringdown). This is made computationally feasible with the machine-learning code Dingo, which accelerates inference by 2-3 orders of magnitude compared to traditional inference techniques. First, when using a uniform prior on the eccentricity, we find eccentric aligned-spin against quasicircular aligned-spin log10 Bayes factors of 1.84 to 4.75 (depending on the glitch mitigation) for GW200129, 3.0 for GW190701 and 1.77 for GW200208_22. We infer egw;10 Hz (egw;20Hz) to be 0.27 & thorn;0.10-0.12 (0.16 & thorn;0.04-0.05) to 0.17 & thorn;0.14-0.13 (0.1 & thorn;0.05-0.04) for GW200129, 0.54 & thorn;0.12-0.30 (0.31 & thorn;0.12-0.13) for GW190701 and 0.39 & thorn;0.23-0.23 (0.21 & thorn;0.08-0.08) for GW200208_22. Second, we find log10 Bayes factors between the eccentric aligned-spin versus quasicircular precessing-spin hypothesis between 1.43 and 4.92 for GW200129, 2.61 for GW190701 and 1.23 for GW200208_22. Third, our analysis does not show evidence for eccentricity in GW190521, which has an eccentric aligned-spin against quasicircular aligned-spin log10 Bayes factor of 0.04. Fourth, we estimate that if we neglect the spin-precession and use an astrophysically motivated prior on the rate of eccentric BBHs, the probability of one out of the 57 events being eccentric is greater than 99.5% or (100-8.4 x 10-4)% (depending on the glitch mitigation). Fifth, we study the impact on parameter estimation when neglecting either eccentricity in quasicircular models or higher modes in eccentric models for GW events. These results underscore the importance of including eccentric parameters in the characterization of BBHs for the upcoming observing runs of the LVK Collaboration and for future detectors on the ground and in space, which will probe a more diverse BBH population.

Simulation-based inference (SBI) is an established approach for performing Bayesian inference on scientific simulators. SBI so far works best on low-dimensional parametric models. However, it is difficult to infer function-valued parameters, which frequently occur in disciplines that model spatiotemporal processes such as the climate and earth sciences. Here, we introduce an approach for efficient posterior estimation, using a Fourier Neural Operator (FNO) architecture with a flow matching objective. We show that our approach, FNOPE, can perform inference of function-valued parameters at a fraction of the simulation budget of state of the art methods. In addition, FNOPE supports posterior evaluation at arbitrary discretizations of the domain, as well as simultaneous estimation of vector-valued parameters. We demonstrate the effectiveness of our approach on several benchmark tasks and a challenging spatial inference task from glaciology. FNOPE extends the applicability of SBI methods to new scientific domains by enabling the inference of function-valued parameters.

A central challenge in many areas of science and engineering is to identify model parameters that are consistent with prior knowledge and empirical data. Bayesian inference offers a principled framework for this task, but can be computationally prohibitive when models are defined by stochastic simulators. Simulation-based Inference (SBI) is a suite of methods developed to overcome this limitation, which has enabled scientific discoveries in fields such as particle physics, astrophysics, and neuroscience. The core idea of SBI is to train neural networks on data generated by a simulator, without requiring access to likelihood evaluations. Once trained, inference is amortized: The neural network can rapidly perform Bayesian inference on empirical observations without requiring additional training or simulations. In this tutorial, we provide a practical guide for practitioners aiming to apply SBI methods. We outline a structured SBI workflow and offer practical guidelines and diagnostic tools for every stage of the process – from setting up the simulator and prior, choosing and training inference networks, to performing inference and validating the results. We illustrate these steps through examples from astrophysics, psychophysics, and neuroscience. This tutorial empowers researchers to apply state-of-the-art SBI methods, facilitating efficient parameter inference for scientific discovery.

Neurons have classically been characterized by their anatomy, electrophysiology, and molecular markers. More recently, single-cell transcriptomics has enabled an increasingly fine genetically defined taxonomy of cortical cell types, but the link between the gene expression of individual cell types and their physiological and anatomical properties remains poorly understood. Here, we develop a hybrid modeling approach to bridge this gap: our approach combines statistical and mechanistic models to predict cells' electrophysiological activity from gene expression patterns. To this end, we fit Hodgkin-Huxley-based models for a wide variety of cortical cell types by using simulation-based inference while overcoming the mismatch between model and data. Using multimodal Patch-seq data, we link the estimated model parameters to gene expression using an interpretable linear sparse regression model. Our approach identifies the expression of specific ion channel genes as predictive of biophysical model parameters including ion channel densities, implicating their mechanistic role in determining neural firing properties.

The meaning of a visual image depends on context-a mouse sees an expanding visual stimulus when approaching a dark refuge or when a cat approaches them, and distinguishing between the two is a matter of life and death. The superior colliculus (SC) is an evolutionarily ancient hub essential for survival behaviors like approach and avoidance of threats.1,2 We therefore combined virtual reality and neural recordings to ask whether matching visual stimuli to self-motion alters behavior and neural activity in SC. We first measured locomotion behavior and neural activity while animals approached an object in virtual reality or while the same object loomed at them. In both contexts, vision dominated activity in the superficial layers of SC (SCs), whereas locomotion had more influence on activity in the intermediate layers (SCim). In addition, animals instinctively slowed their locomotion when nearing the object or when the object neared them. To directly test animals' ability to distinguish self from object motion, we replayed the visual images generated during object approach. Locomotion behavior often changed during replay, showing animals can determine whether visual motion is matched to their self-movement. Further, decoders trained on locomotion behavior or on population activity in SC, particularly in SCim, were able to reliably discriminate self-movement and object movement contexts. We conclude that both mouse behavior and SC activity distinguish the context of visual motion and can thus discriminate motion arising from an animal's own movement and that of an external agent.

Multi-compartment Hodgkin-Huxley models are biophysical models of how electrical signals propagate throughout a neuron, and they form the basis of our knowledge of neural computation at the cellular level. However, these models have many free parameters that must be estimated for each cell, and existing fitting methods rely on intracellular voltage measurements that are highly challenging to obtain in vivo. Recent advances in neural recording technology with high-density probes and arrays enable dense sampling of extracellular voltage from many sites surrounding a neuron, allowing indirect measurement of many compartments of a cell simultaneously. Here, we propose a method for inferring the underlying membrane voltage, biophysical parameters, and the neuron's position relative to the probe, using extracellular measurements alone. We use an Extended Kalman Filter to infer membrane voltage and channel states using efficient, differentiable simulators. Then, we learn the model parameters by maximizing the marginal likelihood using gradient-based methods. We demonstrate the performance of this approach using simulated data and real neuron morphologies.

The ice shelves buttressing the Antarctic ice sheet determine the rate of ice-discharge into the surrounding oceans. The geometry of ice shelves, and hence their buttressing strength, is determined by ice flow as well as by the local surface accumulation and basal melt rates, governed by atmospheric and oceanic conditions. Contemporary methods resolve one of these rates, but typically not both. Moreover, there is little information of how they changed in time. We present a new method to simultaneously infer the surface accumulation and basal melt rates averaged over decadal and centennial timescales. We infer the spatial dependence of these rates along flow line transects using internal stratigraphy observed by radars, using a kinematic forward model of internal stratigraphy. We solve the inverse problem using simulation-based inference (SBI). SBI performs Bayesian inference by training neural networks on simulations of the forward model to approximate the posterior distribution, allowing us to also quantify uncertainties over the inferred parameters. We demonstrate the validity of our method on a synthetic example, and apply it to Ekström Ice Shelf, Antarctica, for which newly acquired radar measurements are available. We obtain posterior distributions of surface accumulation and basal melt averaging over 42, 84, 146, and 188 years before 2022. Our results suggest stable atmospheric and oceanographic conditions over this period in this catchment of Antarctica. Use of observed internal stratigraphy can separate the effects of surface accumulation and basal melt, allowing them to be interpreted in a historical context of the last centuries and beyond.