Estimation of parameters in biophysically detailed neural models with simulation based inference

Abstract

A method, system, and apparatus for estimating parameters for a biological neural network model is disclosed. The method includes identifying model parameters for the biological neural network model and parameter ranges corresponding to the model parameters. Additionally, the method includes sampling from the parameter ranges to generate a dataset of simulated neural patterns and selecting summary statistics for the dataset of simulated neural patterns. Further, the method includes training an artificial neural network based on the summary statistics, and applying the artificial neural network to a neural activity pattern to estimate parameter distributions for simulated neural patterns within the dataset that match the neural activity pattern.

Description

(0312021-00054)ESTIMATION OF PARAMETERS IN BIOPHY SICALLY DETAILED NEURAL MODELS WITH SIMULATION BASED INFERENCECROSS REFERENCE TO RELATED APPLICATION

[0001] The application claims priority to U.S. Provisional Patent Application No. 63 / 567,108, filed March 19, 2024, entitled “ESTIMATION OF PARAMETERS IN BIOPHYSICALLY DETAILED NEURAL MODELS WITH SIMULATION BASED INFERENCE,” the contents of which are incorporated by reference herein in their entirety.GOVERNMENT LICENSE RIGHTS

[0002] This patent application was made with government support under grant numbers R01 AG076227, RF1 MH130415, 1U24NS129945, and R01 EB022889 awarded by the National Institutes of Health. The government has certain rights in this application.BACKGROUND

[0003] Transient low frequency oscillations in the 15-29 Hz beta-frequency band are referred to as Beta Events or beta bursts. Many studies have shown that Beta Events occur throughout the brain and that they often have a stereotypical waveform shape that resembles an inverted Ricker wavelet lasting -150 ms. Changes in transient beta activity have been associated with sensory processing and motor action, and the Beta Event waveform shape has been specifically associated with motor pathology in Parkinson’s disease and with the aging process.

[0004] Local field potential (LFP) motor cortical Beta Events are robustly present during reach and grasp in nonhuman primate (NHP) models. LFP Beta Events exhibit amplitude modulation during reaching tasks. For example, in response to a task cue, Beta Event amplitude is decreased. However, it is unknown what circuit mechanisms underlie Beta Event generation and modulation in the NHP motor cortex. The Human Neocortical Neurosolver (HNN) modeling framework connects macroscale neural signals to underlying cell and circuit-level activity. Previous studies have applied the HNN modeling framework to propose novel mechanisms for the cellular and circuit level generation of Beta Events. Simulation parameters used in the HNN modeling framework define proximal / distal dendritic drives, but predicting parameters fromspecific observations remains a challenge. Parameter distributions are essential for applying real data to biophysical models. The solution space for parameters is non-unique. That is, there exists multiple simulations / model configurations that replicate observed data.

[0005] Understanding the relationship between network connectivity and emergent neural dynamics is a fundamental and unsolved problem in neuroscience. Detailed biophysical models can simulate highly realistic neural circuits with biologically interpretable parameters, providing a powerful way to study neural dynamics. However, the use of these biophysical models is challenged by a large number of model parameters, computationally expensive simulations, and complex mappings from model parameters to simulation outputs. Previous work has demonstrated that deep neural networks (e.g., surrogate models) can be trained to approximate compartmental neuron models, offering simulation speeds that are orders of magnitude faster.SUMMARY

[0006] Example systems, methods, and apparatus are disclosed herein for estimation of parameters in biophysically detailed neural models with simulation based inference (SBI). SBI provides a method to produce posteriori distributions of HNN modeling parameters likely to reproduce observations (e.g., observed EEG data). Examples disclosed herein show how using a combination of SBI and HNN, parameter distributions can be predicted for early proximal followed by late distal dendritic inputs. Modulation of the parameters demonstrates that increased variance and / or strength of proximal dendritic inputs produces larger amplitude Beta Events.

[0007] In examples disclosed herein, surrogate models of individual multicompartment biophysical cells are connected using a spiking neural network (SNN) architecture to connect them into a biophysical network model. This construction permits efficient gradient-based optimization of cell-cell connectivity parameters, and the ability to optimize simulations to complex neural activity patterns. The effectiveness of this approach is demonstrated by using surrogate models to approximate a detailed model of the neocortex, the HNN. HNN is a large-scale detailed model of a cortical column designed to simulate electrical currents (i.e., current dipoles) underlying M / EEG signals with interpretability at the cell and circuit level.

[0008] Disclosed herein, a trained CNN-LSTM (Convolutional Neural Network, Long Short-Term Memory) machine-learning model was used to predict the response (e.g., voltage overtime) of individual neurons to a train of input spikes. In other examples, the machine-learning model used for time-series generation may be one or more of a CNN, an LSTM, a recurrent neural network, a gated recurrent unit network, transformers, or a hybrid architecture. The trained surrogate model disclosed herein can accurately predict compartmental voltages and current dipoles. The single-cell surrogate models disclosed herein may be connected in a cortical microcircuit based on the cell-cell connectivity structure defined by HNN. Techniques from SNN methods were adopted into the surrogate network model. Additionally, the spike threshold function was equipped with an approximate gradient (e.g., the fast sigmoid function) which is critical for performing backpropagation through the surrogate network model.

[0009] Aspects of the subject matter described herein may be useful alone or in combination with one or more other aspects described herein. Without limiting the foregoing description, in a first aspect of the present disclosure, a method for estimating parameters for a biological neural network model includes identifying model parameters for the biological neural network model and parameter ranges corresponding to the model parameters, sampling from the parameter ranges to generate a dataset of simulated neural patterns, selecting summary statistics for the dataset of simulated neural patterns, training an artificial neural network based on the summary statistics, and applying the artificial neural network to a neural activity pattern to estimate parameter distributions for simulated neural patterns within the dataset that match the neural activity pattern.

[0010] In accordance with a second aspect of the present disclosure, which may be used in combination with the first aspect, the method further comprises determining diagnostics for the estimated parameter distributions for the neural activity pattern.

[0011] In accordance with a third aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the diagnostics include one or more of a parameter recovery error (PRE), a posterior predictive check (PPC), and an overlap coefficient (OVL).

[0012] In accordance with a fourth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the biological neural model is one or more of, a biophysically detailed network model, the Human Neocortical Neurosolver (HNN), a spiking neural network, or a neural dynamics model.

[0013] In accordance with a fifth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the artificial neural network estimates the parameter distributions using one or more of normalizing flows, variational autoencoders, autoregressive models, and energy-based models.

[0014] In accordance with a sixth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the neural activity pattern comprises electrophysiological measurements, including one or more of local field potentials (LFP), electroencephalography (EEG), magnetoencephalography (MEG), single-unit spiking, multi-unit- spiking, and current source density.

[0015] In accordance with a seventh aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the neural activity pattern comprises bloodoxygenation-level-dependent (BOLD) signals.

[0016] In accordance with an eighth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the neural activity pattern comprises optical signals including one or more of calcium imaging and voltage sensors.

[0017] In accordance with a ninth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, selecting the summary statistics includes mathematically transforming the simulated neural patterns into a latent-space representation.

[0018] In accordance with a tenth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the neural activity pattern is one of a real neural activity pattern or a simulated neural activity pattern.

[0019] In accordance with an eleventh aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, a method for approximating simulations of a biological neural network model includes training an artificial neural network to approximate a surrogate model, the surrogate model corresponding to an individual cell of the biological neural network model, connecting one or more duplicates of the surrogate model into a surrogate network model, and optimizing the surrogate network model to fit simulations to target activity patterns.

[0020] In accordance with a twelfth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the biological neural network model is one ormore of, a biophysically detailed network model, the Human Neocortical Neurosolver (HNN), a spiking neural network, or a neural dynamics model.

[0021] In accordance with a thirteenth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the artificial neural network uses one or more of recurrent neural networks, convolutional neural networks, long short-term memory networks, gated recurrent unit networks, transformers, or hybrid architectures for time series generation.

[0022] In accordance with a fourteenth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the optimizing of the surrogate network model includes one or more of gradient-based methods, backpropagation, evolutionary algorithms, and Bayesian optimization.

[0023] In accordance with a fifteenth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the optimizing of the surrogate network model includes identifying model parameters for the biological neural network model and parameter ranges corresponding to the model parameters, sampling from the parameter ranges to generate a dataset of simulated neural patterns, selecting summary statistics for the dataset of simulated neural patterns, training an artificial neural network based on the summary statistics, and applying the artificial neural network to a neural activity pattern to estimate parameter distributions for simulated neural patterns within the dataset that match the neural activity pattern.

[0024] In accordance with a sixteenth aspect of the present disclosure, which may be used in combination with any other aspect disclosed herein, the method further includes combining the surrogate model and an output of the optimization into a single artificial neural network model.

[0025] In accordance with a seventeenth aspect of the present disclosure, any of the structure and functionality illustrated and described in connection with FIGS. 1 to 20 may be used in combination with any of the structure and functionality illustrated and described in connection with any of the other of FIGS. 1 to 20 and with any one or more of the preceding aspects.

[0026] In light of the present disclosure and the above aspects, an advantage of the present disclosure is the ability to define complex neural activity patterns as the optimization objective, enabling rapid identification of connectivity parameters that produce desired outputs.

[0027] Additional benefits of the combination of HNN and SBI disclosed herein include determining parameter optimization for pharmacokinetics. For example, methods disclosed herein can be used to estimate how a patient responds to a certain neuro-stimulant drug using an accurate neural model of the patient created from optimized parameters. Another benefit of the examples disclosed herein includes determining parameters to target as biomarkers. For example, the models disclosed herein can reveal which are the parameters key to having the model approximate recorded EEG data. Further benefits of the examples disclosed herein include the ability to create a simplified patient model that produces results almost identical to recorded EEG data.

[0028] Additional features and advantages are described in, and will be apparent from, the following Detailed Description and the Figures. The features and advantages described herein are not all-inclusive and, in particular, many additional features and advantages will be apparent to one of ordinary skill in the art in view of the figures and description. Also, any particular embodiment does not have to have all of the advantages listed herein and it is expressly contemplated to claim individual advantageous embodiments separately. Moreover, it should be noted that the language used in the specification has been selected principally for readability and instructional purposes, and not to limit the scope of the inventive subject matter.BRIEF DESCRIPTION OF THE FIGURES

[0029] Fig. 1 is a diagram of a simulation based inference (SBI) workflow to estimate parameters and possible parameter combinations, according to an example embodiment of the present disclosure.

[0030] Fig. 2 is a diagram showing a schematic of a Human Neocortical Neurosolver (HNN) model, according to an example embodiment of the present disclosure.

[0031] Figs. 3A-3C are diagrams showing resistor-capacitor (RC) circuit simulations, according to an example embodiment of the present disclosure.

[0032] Figs. 4A-4D are diagrams showing diagnostics to compare summary statistics on an RC circuit model, according to an example embodiment of the present disclosure.

[0033] Figs. 5A-5C are diagrams showing HNN simulations that mimic an RC circuit, according to an example embodiment of the present disclosure.

[0034] Figs. 6A-6D are diagrams showing SBI diagnostics of summary statistics in HNN, according to an example embodiment of the present disclosure.

[0035] Figs. 7A-7D are diagrams showing a HNN-SBI recovery circuit mechanism of Beta Event amplitude, according to an example embodiment of the present disclosure.

[0036] Figs. 8A-8E are diagrams showing inference of local connectivity parameters from event related potential (ERP) waveforms, according to an example embodiment of the present disclosure.

[0037] Fig. 9 is a table showing simulation parameters and SBI training, according to an example embodiment of the present disclosure.

[0038] Figs. 10A-10C illustrate a surrogate model which accurately predicts compartmental voltage responses and net current dipole of individual neurons in HNN, according to an example embodiment of the present disclosure.

[0039] Figs. 11A-11F illustrate connectivity parameters optimized by the model of Fig. 10, according to an example embodiment of the present disclosure.

[0040] Fig. 12 is a diagram illustrating the functional significance of Beta Events, according to an example embodiment of the present disclosure.

[0041] Fig. 13 is a diagram illustrating local field potential (LFP) motor cortical Beta Events robustly present during reach and grasp in nonhuman primates (NHP), according to an example embodiment of the present disclosure.

[0042] Fig. 14 is a diagram illustrating LFP Beta Events exhibiting amplitude modulation during a reaching task, according to an example embodiment of the present disclosure.

[0043] Fig. 15 is a diagram illustrating how HNN connects macroscale signals to underlying neural activity, according to an example embodiment of the present disclosure.

[0044] Fig. 16 is a diagram illustrating early proximal and late distal dendritic inputs, according to an example embodiment of the present disclosure.

[0045] Fig. 17 is a diagram illustrating how SBI may be used to determine parameter distributions where solution space is non-unique, according to an example embodiment of the present disclosure.

[0046] Fig. 18 is a diagram illustrating observed LFP Beta Events, according to an example embodiment of the present disclosure.

[0047] Fig. 19 is a diagram illustrating predicted parameter distributions generated by SBI and HNN and modulation of Beta Event amplitude based on proximal dendritic inputs, according to an example embodiment of the present disclosure.

[0048] Fig. 20 is a system diagram illustrating a computer system capable of performing the methods disclosed herein, according to an example embodiment of the present disclosure.DETAILED DESCRIPTION

[0049] The terms “a,” “an,” “the” and similar referents used in the context of describing the invention (especially in the context of the following claims) are to be construed to cover both the singular and the plural, unless otherwise indicated herein or clearly contradicted by context.

[0050] All methods described herein can be performed in any suitable order unless otherwise indicated herein or otherwise clearly contradicted by context.

[0051] The use of any and all examples, or exemplary language (e.g., “such as”) provided herein is intended merely to better illuminate the invention and does not pose a limitation on the scope of the invention otherwise claimed. No language in the specification should be construed as indicating any non-claimed element essential to the practice of the invention.

[0052] The terms "comprise", "comprises", "comprised" or "comprising", "including" or "having" and the like in the present specification and claims are used in an inclusive sense, that is to specify the presence of the stated features but not preclude the presence of additional or further features.Introduction

[0053] Methods, systems, and apparatus are disclosed herein for estimation of parameters in biophysically detailed neural models with simulation based inference (SBI). Biophysically detailed neural modeling is a known framework used to study fast time scale neural dynamics. While challenging to construct, advances in computational resources have enabled the proliferation of detailed models from principled models of single neurons to large scale biophysically detailed neural networks that enable multi-scale interpretation from cell spiking tolocal field potentials to macroscale magneto- and electroencephalographic (MEG / EEG) signals. Numerous detailed models are now openly distributed to encourage their use and expansion.

[0054] A common goal of detailed neural modeling is to infer biophysical parameters in individual cells and / or network connections that can account for observed changes in neural activity over time. Given the large-scale nature of any detailed model, parameter inference is an inherently challenging problem with no solution by known methods. The difficulty of parameter inference is closely tied to the level of biophysical detail in the model. For example, as the number of parameters increases with more realistic models, the difficulty of parameter inference increases. In practice, parameters cannot all be estimated simultaneously, but model elements may be estimated in a serial fashion (e.g. cell dynamics followed by network connectivity) and fixed. Then, a limited set of parameters are chosen as the target for estimation. Even with this limited set, the parameter estimation process is complex. The problem is confounded by the fact that there may be many parameter configurations that produce an equally good representation of the data. Identifying unique biophysically constrained parameter sets that can account for observed neural dynamics, and differences across experimental conditions, is essential to the meaningful use of large scale biophysically detailed models for neuroscience research. For example, if a biophysically detailed model is to be used to infer circuit level mechanisms generating an EEG waveform that is a biomarker of a healthy compared to neuropathological condition, an approach not only to estimate the parameter distributions that can generate the waveforms but also to assess if the distributions are distinguishable is needed.

[0055] A powerful approach for estimating parameters in neural models is Bayesian inference. There is an extensive history of research applying Bayesian inference, and specifically the algorithm of variational inference, to estimate parameters in reduced models of neural activity, for example in the Dynamic Causal Modeling (DCM) framework that relies on reduced “neural mass models.” However, while compatible with reduced models that are mathematically tractable, the algorithm of variational inference is not compatible with detailed biophysical models due to its computational complexity, and specifically lack of access to a known likelihood function (see also Discussion). Additionally, mean field variational inference assumes that model parameters are independent, preventing the ability to capture biologically meaningful parameter interactions.

[0056] Simulation based inference (SBI) has been proposed as an alternative Bayesian inference framework to estimate parameters in detailed neural models. SBI overcomes the challenge of the lack of a likelihood function by leveraging advances in deep learning to perform density estimation. An advantage of SBI is that it only relies on a dataset of simulations from the model being investigated, rather than requiring knowledge of a parameter likelihood function, which is typically not accessible in large-scale biophysically detailed models. From the dataset of simulations, a neural density estimator (i.e., an artificial neural network specifically made to approximate probability distribution functions) can be trained to learn a mapping of observed neural dynamics (e.g., time series waveforms) to corresponding model parameter distributions. Another advantage of this approach is that SBI estimates a full distribution over model parameters that may account for the data and provides information about parameter interactions. Determining information about parameter interactions is not possible with optimization techniques that have historically been used in large-scale biophysical models such as COBYLA and genetic algorithms, because these techniques estimate only a single parameter configuration that best fits the data.

[0057] Fig. 1 outlines an overall SBI workflow 100 to estimate parameters and possible parameter combinations that can reproduce recorded neural activity, according to an example embodiment of the present disclosure. In a first step 102, a prior distribution of relevant model parameters and ranges is constructed. In a second step 104, a dataset of simulated neural activity patterns is generated with parameters sampled from the prior distribution. For example, parameter ranges are randomly sampled to simulate real world neural waveforms. In a third step 106, user defined summary statistics are chosen to describe waveform features of interest. In a fourth step 108, an artificial neural network (e.g., a specialized deep learning architecture) is trained to learn the mapping from neural activity constrained by summary statistics to underlying model parameters. In a fifth step 110, specific neural activity patterns of interest are fed into the trained artificial neural network, which subsequently estimates the parameter distributions for waveforms (e.g., neural activity patterns) of interest. For example, the artificial neural network outputs a distribution over the potential underlying model parameters for a neural activity pattern of interest. In a sixth step 112, parameter estimates for different waveforms can be compared. For example, diagnostics like parameter recovery (if the ground truth is known) or posterior predictive checks can be run.

[0058] Disclosed herein is an approach for applying SBI to estimate parameters underlying time series waveforms generated by biophysically detailed neural models. Important considerations for this approach include identifying the parameters and ranges over which the inference process will be performed (i.e., prior distribution), which necessarily depends on user- defined hypotheses, and selecting informed summary statistics of the waveform activity. Examples disclosed herein describe how diagnostics can be used to assess the uniqueness and quality of the posterior estimates and to assess the overlap of distributions from estimation applied to two different waveforms. These evaluation steps are important to resolve the uniqueness of distributions for two or more waveforms.

[0059] Examples disclosed herein are first described using a simplified example of a nonlinear resistor-capacitor circuit, and then extended to an example large-scale biophysically detailed modeling framework. The example large-scale biophysically detailed modeling framework, the Human Neocortical Neurosolver (HNN), was developed to study the multi-scale neural origin of human MEG / EEG signals. The foundation of HNN is a biophysically detailed model of a neocortical column, with layer specific synaptic activation representing thalamocortical and corti co-corti cal drive (see Fig. 2 described below). HNN has been applied to study the cell and circuit origin of commonly measured MEG / EEG signals, including low frequency oscillations (e.g. Beta Events and event related potentials (ERPs)), along with differences across experimental conditions. Examples disclosed herein demonstrate application of SBI to estimate parameter distributions that can account for variation in example hypothetical Beta Event and in ERP waveforms with selected parameter priors (see the first step 102 of Fig. 1) based on previous studies. Examples disclosed herein show that due to the model complexity some parameters can be inferred uniquely while others are indeterminate. The methods described herein may be used to guide future applications of SBI in a wide variety of applications that use detailed models to study neural dynamics.Materials and MethodsResistor-Capacitor Circuit Simulations

[0060] Various resistor-capacitor (RC) circuit examples are described herein. RC circuits are a building block of HNN-type models and offer an SBI setup with time series inputs and the presence of indeterminacies. RC circuit simulations may be performed using the odeint ordinarydifferential equation (ODE) solver of the SciPy Python package. Described below in the results section is a description of how RC circuit simulations were performed in examples disclosed herein.Human Neocortical Neurosolver

[0061] Biophysical modeling of cortical activity underlying MEG / EEG signals was performed using the EINN. The example HNN model described herein comprises 100 pyramidal neurons and 33 inhibitory neurons in both layers 2 / 3 and 5 for a total of 266 neurons and represents a localized small patch of neocortex. To accurately reproduce macroscale electrical signals, HNN utilizes multi-compartment pyramidal neuron models, and synchronous intracellular current flow of aligned layer 2 and 5 pyramidal neuron dendrites is assumed to be the generator of the primary electrical current sources underlying recorded extracranial MEG / EEG signals, due to their length and parallel alignment.

[0062] Figs. 2A-2D illustrate a schematic of a Human Neocortical Neurosolver (HNN) model. Fig. 2A shows a diagram 200 of local network connections of the HNN model including excitatory AMPA and NMDA synaptic connections 202 originating from pyramidal neurons 204 and inhibitory GABAa and GABAb synaptic connections 206 originating from inhibitory interneurons 208. The inhibitory basket cells 208 are modeled as single compartment point neurons given their negligible impact in producing the recorded electrical currents, but are none-the-less crucial to the local network dynamics. The pyramidal cells 204 in the model connect to other cells with both AMPA and NMDA synapses, while basket cells produce GABAa and GABAb synaptic connections.

[0063] In addition to the local circuitry, HNN models extrinsic inputs to the column via layer specific synaptic excitatory drives generated by predefined patterns of action potentials presumed to come from exogenous brain areas (e.g., thalamus and other cortical regions). In general these are referred to as proximal and distal drives, reflecting the effective location of the synapses on the pyramidal neuron proximal and distal dendrites. Fig. 2B shows a diagram 210 of a proximal exogenous input connection pattern. Fig. 2C shows a diagram 220 of distal exogenous input connection pattern. The proximal drive reflects so called feedforward inputs from lemniscal thalamus, and the distal drive reflecting inputs from either non-lemniscal thalamus, or “feedback”connections from other cortical regions. These proximal and distal drives induce excitatory post- synaptic currents that drive current flow up and down the aligned pyramidal cell dendrites (see red and green arrows in Figs 2B-2C) that can generate positive and negative deflections in the primary electric current dipole of source localized MEG / EEG signals. Intracellular current flow due to these extrinsic inputs, as well as induced local spiking dynamics, combine to produce the recorded MEG / EEG signal. Fig. 2D shows a diagram 230 of a 3D rendering of the full neocortical column model.Posterior Diagnostics

[0064] When working with a posterior approximation (), it may be useful to characterize its behavior in different regions of the parameter space. For a given parameter configuration (0o) and simulated output (xo ~ p(x | 0o)), how concentrated 0(0 | xo) is around 0o can be quantified. This quantity of is referred to as the parameter recovery error (PRE) is defined as the Wasserstein distance between the Uth marginal of (0 | xo), and a Dirac delta centered at 0o[k], This can be empirically estimated as shown in Equation 1 below.

[0065] (Equation 1)

[0066] Using Equation 1, N samples {& , . . &N} <D(i / x0) are generated from the posterior approximation conditioned at x0- p(x I o), which is an observation from the simulator at ground truth i%. Note that &ais composed of k distinct values for each individual parameter, therefore there are k distinct PRE values. Additionally, each parameter i?0[k] is mapped from its range defined in the prior distribution, to the range (0,1). Therefore, the maximum PRE is 1.0, indicating the worst possible recovery, whereas a PRE of 0.0 indicates perfect recovery.

[0067] The previous diagnostic quantified how well represents the relationship between x0and &0in terms of the parameter space. Alternatively, the relationship in the observation space can be assessed. Specifically, given samples from the approximate posterior &, ~ 07 / x0), simulations x, ~ p(x i?,) are generated and assessed for how close x, is to the conditioningobservation0. This is known as a posterior predictive check (PPC) and is quantified as the root mean squared error between x0and the x, as shown in Equation 2 below.

[0068] (Equation 2)

[0069] Further, in many applications it is useful to characterize how well two distributions can be distinguished from one another. Accordingly, a distribution overlap coefficient (OVL) may be introduced. Given the approximate posterior , OVL is defined as shown in Equation 3 below. In Equation 3, xo and xi are two different observations whose posterior distributions may be compared on the marginal distribution for the k-th parameter. To calculate the OVLk numerically, an evenly spaced grid of 50d for a prior distribution over d E N+ parameters may be used.

[0070] (Equation 3)ResultsApproach to applying SBI in biophysically detailed models that simulate time series data

[0071] SBI has been established as an approach to estimate parameters in detailed biophysical models that simulate time series data. This approach overcomes the challenges of applying Bayesian inference in highly detailed non-linear models, namely estimation of complex posterior distributions that exhibit parameter interactions, by leveraging advances in likelihood- free inference and deep learning. Described below is a review of the SBI process with the mathematical description of each of the steps outlined in Fig. 1 provided herein.

[0072] SBI is used as disclosed herein to estimate parameters and possible parameter distributions that can account for an observed neural dynamic (e.g. time series waveform). In mathematical terms, this goal is stated as follows. Given an observation x and model with parameters 0, SBI seeks to create an approximation 0(0 | x) of a posterior parameter distribution p(9 | x) such that Equation 4 below is satisfied. As shown in Equation 4, Bayes’ rule specifies aclosed form for the desired posterior distribution p(0 | x) as being proportional to the likelihood p(x | 0) multiplied with the prior p(0).

[0073] (Equation 4)

[0074] In detailed biophysical models the likelihood function p(x | 0), which encodes the relationship between model parameters 9 and outputs x, is often analytically intractable but can be approximated from a large number of simulations. The novelty of SBI is that it circumvents likelihood evaluations altogether, and instead approximates the posterior distribution directly from a simulated dataset of model outputs xi ~ p(x | 0i) with parameter values sampled from a user defined prior distribution 0i ~ p(0). There are numerous known approaches to achieve this goal, each with their own unique considerations, benefits, and challenges. Disclosed herein are examples that use a deep learning architecture known as a conditional neural density estimator (0), which is a function that takes an observation x as input, and returns a probability density function defined over the parameter space. More specifically, the conditional neural density estimator utilizes normalizing flows following standard practices described herein in order to estimate the probability density function. In other examples, the conditional neural density estimator may utilize one or more of variational autoencoders, autoregressive models and energy-based models to perform the probability density estimation. The detailed steps in applying SBI (Fig. 1) are as follows.Steps 1-2: Define prior and generate training data

[0075] SBI begins with the user choosing the parameters of interest to be estimated, and the range and statistical distribution of values (e.g. uniform distribution) those parameters can take. This constitutes the prior distribution p(0) over model parameters 0 to be inferred. A well- constructed prior distribution is important because it encodes the assumptions and hypotheses of the inference problem considered, and strongly impacts any resulting predictions. Creating a good prior distribution requires domain expertise to choose meaningful parameters, and a biologically realistic range of values. Even so, predictions are not predetermined by the prior, as uncertainty can be encoded using flat and / or uninformative priors where the probability mass is evenly spread over the desired parameter range. Nevertheless, this aspect is highly important for detailed neural models where the inferred parameters represent a small subset of the total set of parameters.

[0076] With the prior constructed, a simulated dataset of observations xl :N (e.g., simulated neural patterns) is generated by simulating a large number of N time series using model parameters 01 :N drawn from the prior parameter distribution.Steps 3-4: Training neural network based on chosen summary statistics to describe the time series waveform

[0077] With the simulated dataset, a specialized deep learning architecture known as a conditional neural density estimator is trained to approximate the posterior distribution p(0|x) that can account for chosen summary statistics for any observation x. The output of the conditional neural density estimator is a distribution of parameters which can generate simulations close to the summary statistics of the conditioning observation xO. The neural density estimator <I» can be trained directly from the entire time series x, or from summary statistics s = S(x) which constitute a lower dimensional vector of values.

[0078] The choice of the summary statistics (s) should aim to obtain (9 | s) ~ p(9 | s) ~ p(9 | x), meaning that the posterior estimates are well enough approximated from s alone. A more in depth description of the choice and role of summary statistics is given after describing Steps 5- 6.Steps 5-6: Estimate and compare parameter distributions for distinct waveforms

[0079] With a trained neural density estimator, users can feed in new waveforms and assess the predicted parameter distributions underlying their generation. The new waveforms represent neural activity patterns which may be an electrophysiological measurement, including one or more of local field potentials (LFP), electroencephalography (EEG), magnetoencephalography (MEG), single-unit spiking, multi-unit-spiking, and current source density. In other examples, the waveform represents a neural activity pattern of blood-oxygenation-level-dependent (BOLD) signals. Additionally, in other examples, the neural activity pattern may be an optical signal such as calcium imaging and / or voltage sensors. The neural activity patterns fed into the trained neural density estimator may be real neural data of one of the types listed above or may be simulated neural data for one of these types. The simulated neural patterns of the simulated dataset may correspond to the type of neural activity pattern of the new waveform.

[0080] Diagnostics that assess the quality of the distribution can then be performed (see Materials and methods for details on the calculation of each diagnostic). If the ground truthparameters underlying the waveform of interest are known, users can calculate the dispersion of the posterior around this ground truth using parameter recovery error (PRE). If the ground truth is unknown, users can use posterior predictive checks (PPC) to assess if the parameter distribution consistently produces waveforms close to the conditioning observation. Finally, uniqueness of posterior estimates for two different waveforms can be assessed by directly comparing the overlap coefficient of the distributions (OVL).The important role of summary statistics in parameter inference

[0081] Using the inference framework described above, examples disclosed herein emphasize the role of summary statistics and how their selection directly impacts predicted parameter distributions. Inference on models with full time series outputs is challenging because the observations are high dimensional. This challenge comes in the form of interpretability and model misspecification due the simulator not capturing finer characteristics of the real data generating process.

[0082] Summary statistics can either be hand-crafted, leveraging domain expertise and hypotheses regarding the data, or they can be automatically extracted. A summary statistic s = S(x) is sufficient if all the relevant information for mapping the full observation to its underlying parameters is retained. More specifically, sufficiency is satisfied if p(0 | s) = p(0 | x). In practice, truly sufficient summary statistics are rare, but they can still provide a close approximation to inferences achieved when using the full observations x. In the case of hand-crafted features, the hand-crafted features can be readily interpretable, and come associated with hypotheses on their physiologic significance depending on previous research. Alternatively, full time series informed approaches like principal component analysis (PC A) may do a better job at retaining more complex relationships between observations and parameters. PCA is a known mathematical technique which involves transforming data into a latent-space representation. As an alternative to PCA, automatic extraction of summary statistics through embedding networks may also be used.

[0083] Given the current ambiguity around summary statistic selection, principled approaches to compare alternate approximations of the posterior distribution are essential. To this end, disclosed herein are examples which allow for an intuitive understanding of the role of different summary statistics. Additionally, diagnostic analyses that can be used to quantitatively compare desirable properties of posterior approximations produced using different summarystatistics, such as PPC, PRE, and OVL are introduced. These analyses are detailed in the examples below.RC circuit: a simple example to describe indeterminacies that can occur with time series inference

[0084] A difficulty with likelihood free inference is that the models studied typically do not permit access to a ground truth posterior distribution, over which inferences can be validated. To highlight this challenge and better understand how decisions in the SBI pipeline impact the resulting approximate posterior, the SBI pipeline is first applied on a model where the ground truth posterior is known; namely an RC circuit model.

[0085] The equation for the RC circuit simulations is shown in Equation 5 below.

[0087] In Equation 5, V(t) represents the voltage response of an RC circuit to a current injection Ie(t). In this example, capacitance C = 6 F, resistance R = 1 , and constant voltage source E = 0 V are used.

[0088] Figs. 3A-3C illustrate RC circuit simulations. In Fig. 3A, simulated voltage and current dipole waveforms are shown for two exemplar parameter configurations with latency between positive 1+ and negative I- square current injections, At = 0 (waveform 312, 322), and At 0 (waveform 314, 324). An example chart 310 of the original ’’Raw” waveforms, as well as an example chart 320 of the PCA transformed waveforms with the first 30 components (PCA30) are shown to demonstrate that this summary statistic retains almost identical information. Fig. 3B illustrates posterior distributions 330 showing the inferred values that can generate the waveforms from Fig. 3A (PCA30 used to generate distributions) and demonstrates that when the latency At between the inputs is zero, their amplitudes are indeterminate as visible with a positive correlation between the parameters 1+ and I-. Fig. 3C shows a schematic 340 of how the RC circuit is driven by positive 1+ and negative I- square current injections, where the amplitude and latency At between pulses serve as parameters. This simulation parallels HNN simulations described below(Fig. 5) in which a single excitatory proximal / distal input with variable synaptic conductances and latencies produce positive / negative deflections in the current dipole.

[0089] The example RC circuit simulations described here are parameterized in a way that enables comparison to similar simulations in the biophysically detailed simulations described below for HNN (Fig. 3C). For example, the RC circuit is driven with two square wave pulse injections, with positive and negative amplitudes lasting 20 ms each. Two parameters (1+ and I ) controlled the amplitude of positive and negative square pulse current injections, such that the sum le determined the final injected current for each time step. A third parameter, latency At, controlled the time delay between the two pulses. Specifically, At shifts the negative current pulse in time, while the positive current pulse remains fixed. These inputs play a similar role as excitatory proximal and distal inputs in HNN (Figs 2B-2C). Table 900 of Fig. 9 details the prior distribution over the parameters used to generate training examples for SBI.

[0090] The relationship between current injection amplitude and the RC circuit voltage response strongly depends on if At is zero or non-zero. When At is zero, the total injected current sums to a single square pulse since the negative and positive pulses perfectly overlap. Waveform 312 of Fig. 3A shows an example of this where a -0.2 mA square pulse current injection is delivered at 80 ms producing a voltage response with an exponential rise and decay with one peak. Since any combination of 1+ and I- which sum to -0.2 mA will produce an identical voltage, there will be an indeterminacy when attempting to infer these parameters from the voltage response. In contrast, waveform 314 of Fig. 3A shows an example with a non-zero latency At 0. Specifically, the first current injection with 1+ = 0.3 mA is delivered at 80 ms, and the second current injection with I- = 0.5 mA at 117.5 ms, therefore At = 37.5 ms. The voltage response exhibits two unique peaks due to the offset between the square pulses. Since there is only one combination of 1+ and I- that can produce this waveform, their values can be inferred exactly from the waveform. In other words, the amplitude parameters underlying the voltage response can be inferred exactly only when At 0. Note that since we are approximating the posterior distribution, even values close to At ~ 0 will still produce an indeterminacy.

[0091] To visualize and interpret posterior distributions produced by SBI, a sufficient number of random samples from the distribution must first be drawn. Since posterior distributions are often multidimensional, it is useful to plot the samples using a “pair-plot” 330. The diagonalof the pair-plot 330 is used to visualize the univariate (e.g., marginal) distribution of each parameter, here using a kernel density estimate of n = 1000 posterior samples. The squares below the diagonal visualize the bivariate distribution for pairs of parameters by plotting the samples explicitly on a scatter plot. This example highlights that even with simple simulators, indeterminacies can easily arise necessitating the use of flexible posterior approximators (e.g., masked autoregressive flows) compared to mean field variational inference which cannot handle parameter interactions.

[0092] Note that PC A with 30 components (PCA30) was used as the summary statistic in this example (chart 320) rather than the full time series to avoid the potential computational issues of conditioning posteriors on high dimensional data, while still retaining the majority of the waveform variance (explained variance=0.883). Chart 320 plots the inverse transformed PCA30 waveform to highlight that the summary statistic retains almost identical information.

[0093] The results in Fig. 3B show that the expected posterior distribution described above can be recovered with SBI. As shown in Fig 3B(a-e), when conditioned on the voltage response with At = 0, any distribution involving 1+ or I- (data 332, 334, 336, 338, 340, 342) will exhibit an indeterminacy (i.e., multiple recovered values along one dimension). The high correlation between 1+ and I- of 0.998 (p ~ 0) demonstrates that the indeterminacy is characterized by a strong linear interaction between these parameters. Specifically, the line 332 in Fig 3B corresponds to all values in which the amplitudes sum to a constant value of V = -0.2, and the resulting voltage waveform is identical. In contrast, the voltage response with At 0 (data 344, 346, 348, 350, 352, 354) produces a posterior distribution concentrated on a single point 346 around the ground truth parameters, correlation between 1+ and I-: 0.375; p < le-33).Diagnostics enable comparison of posterior estimates using different summary statistics

[0094] In the previous example, PCA30 was utilized as a summary statistic to learn a low dimensional representation of the voltage time series. PCA is a common choice for dimensionality reduction, and has been used in historical MEG / EEG inference work with only the first 3-4 principal components. However, it is not guaranteed that PCA, which aims to only capture variance, will retain the features that best allow SBI to map waveforms to simulation parameters. An alternative approach is to leverage domain-expertise to construct hand-crafted summarystatistics specific to the model and inference problem. Unfortunately, it cannot be known a priori which summary statistic will allow SBI to perform best, necessitating quantitative diagnostics that allow a systematic comparison. Disclosed herein are two simple hand-crafted summary statistics, as well as the posterior diagnostics PRE and PPC, which can build an intuitive understanding of how emphasizing different summary statistics can impact the final estimates produced by SBI.

[0095] Figs. 4A-4D illustrate diagnostics to compare summary statistics on the RC circuit model. In Fig. 4A, summary statistics applied to the simulated time series included: PCA30 402, PCA4 404, Peak 406 (amplitude and timing of max / min) 406, and BandPower 408 (four bands between dotted lines). PCA30 402 and PCA4 404 plots show the associated inverse transformed signal. Two exemplar simulations with pulse latencies At = 0 (waveforms 402a, 404a, 406a, 408a) and At 0 (waveforms 402b, 404b, 406b, 408b) are shown. Fig. 4B illustrates conditioning the approximate posterior distribution on the At = 0 time series produces indeterminacies for all summary statistics, such that the ground truth (dotted lines) current injection amplitudes (1+ and I-) cannot be uniquely recovered. PCA30 412, PCA4 414, and Peak 416 features exhibit a linear interaction between parameters for the At = 0 time series (data 412a, 414a, 416a, 418a), whereas the ground truth is recovered for the At 0 (data 412b, 414b, 416b, 418b) time series. BandPower 418 produces non-linear interactions for both time series.

[0096] The first hand-crafted summary statistic, peak 406, is defined as a four dimensional vector including the amplitude and timing of the maximum and minimum peaks of the simulated voltage response. It can be readily seen that these features reflect the underlying simulation parameters. Upon visual inspection of the voltage response with At 0 (waveform 406b), the height of the maximum and minimum peak directly correspond to the parameters 1+ and I-, and the distance between these peaks correspond to the latency parameter At. As discussed below, Peak 406 features permit inference that is close to that achieved with PCA30 402 and also PCA4 404. The second hand-crafted summary statistic is defined as a four dimensional vector including the band power (BandPower 408) of common frequency ranges used to study neural oscillations. Specifically, the beta (13-30 Hz), low-gamma (30-50 Hz), and high-gamma (50-80 Hz) ranges were considered, as well as the aggregate band power of the alpha and lower frequency ranges (0- 13 Hz). As discussed below, this feature was intentionally selected as a cautionary example of asummary statistic that is ill-suited for the inference problem, but has a basis in previous neuroscience and Bayesian inference literature.

[0097] Described next are two diagnostics that allow comparison of desirable properties of the posterior distribution for different summary statistics, namely PRE and PPC (see the sixth step 112 of Fig. 1). PRE values were calculated over a grid, with 10 points for every parameter dimension, spanning the range of the prior distribution. Brighter colors indicate higher dispersion of the posterior around the ground truth parameters defined by each square. Errors tend to be concentrated around At = 0 for PCA30, PCA4, and Peak features. In Fig. 4C, local parameter recovery error (PRE) heatmaps are shown over the parameter grid for PCA30 422, PCA4 424, Peak 426 and BandPower 428 and in Fig. 4D, local posterior predictive check (PPC) heatmaps are shown for PCA30 432, PCA4 434, Peak 436 and BandPower 438.

[0098] One of the advantages of inspecting local posterior diagnostics (“local” as in specific to the pair (xO, 00)) is the ability to identify patterns. Fig. 4C plots the PRE for the I- parameter, with respect to different ground truth values of 1+ itself, and the At parameter. A first square 422a and a second square 422b mark the ground truth values used to generate the waveforms in Fig. 4A. It was observed that the summary statistics PCA30 422, PCA4 424, and Peak 426 all exhibit a pattern where At values near zero produce a larger PRE compared to the rest of the parameter grid. This is due to the indeterminacy in 1+ as seen in the posterior samples for At = 0 (data 412a, 414a, 416a, 418a) of Fig 4B. It is apparent, however, that PCA30 values produce the lowest PRE values, even near a At of zero. In contrast, the BandPower 428 summary statistic produces a posterior distribution with complex indeterminacies for both observations. This results in high PRE values across the entire parameter grid as seen in BandPower 428, indicating that this summary statistic is not effective at recovering the ground truth parameters.

[0099] Defining gPRE as the global average value over this grid, the gPRE values for the parameter 1+ were low for PCA30 (0.01 ± 0.02), PCA4 (0.03 ± 0.04), and Peak (0.04 ± 0.04) features. In contrast, the larger gPRE for Band Power features (0.23 ± 0.14) demonstrates that the ground truth for 1+ is recovered much less accurately. The gPRE for I- followed a similar pattern, whereas the parameter At was generally well recovered for all ground truths (see Data and code availability for results of diagnostics for all parameters and conditions).

[0100] The PPC is a method to describe how well samples from the posterior match the conditioning observation. Given a well-estimated posterior distribution p(0 | x) and a conditioning observation xO ~ p(x | 00), one would expect simulations xi ~ p(x0 | 0i) to be close to the original conditioning observation. Unlike the PRE heatmaps, the PPC plots shown in Fig 4D do not exhibit obvious patterns with respect to the underlying parameter grid, and instead the summary statistics are well characterized by the global average PPC (gPPC). Similar to the PRE analysis, the gPPC values were relatively low for PCA30 (0.008 ± 0.004 mV), PCA4 (0.012 ± 0.006 mV), and Peak (0.027 ± 0.015 mV) summary statistics, whereas BandPower (0.133 ± 0.053 mV) was substantially larger.

[0101] These diagnostics demonstrate PCA30 performs the best for this inference problem. Additionally, the local PRE analysis revealed differences between summary statistics that were not apparent with the global diagnostics. It is important to note that neither of these diagnostics quantify the closeness of the approximation to the true posterior. For instance, if there is an interaction between parameters of the model causing a parameter indeterminacy, then the PRE will always be non-zero, since the posterior distribution will be spread in the parameter space. This is the case for the posterior of the RC circuit with At = 0 in Fig 3B. Similarly, if model simulations are stochastic, then a given set of parameters may map to multiple equally valid outputs, producing a non-zero PPC. Nevertheless, both diagnostics provide useful information to compare desirable properties of the posterior approximation.Subthreshold HNN simulations mimicking the RC circuit shows that inference with summary statistics that account for the full time series waveform perform best

[0102] Building from the RC circuit, described below is a nearly identical inference problem in our large-scale biophysically detailed model constructed to study the neural mechanisms of human MEG / EEG, HNN. Figs. 5A-5C illustrate HNN simulations that mimic an RC circuit. The HNN simulations of Figs. 5A-5C reflect the nearly identical parameter configuration as the RC circuit in Fig. 3. In Fig. 5 A, simulated current dipole waveforms are shown for two exemplar parameter configurations with At = 0 (waveforms 514, 524) and At0 (waveforms 512, 522). The original “Raw” simulated waveform 510 is plotted in comparison with the PCA inverse transformed waveform 520 with 30 components. In Fig. 5B, posteriordistributions 530 showing the inferred values that can generate the waveforms from panel A demonstrate that when the latency between the inputs is zero (data 531, 532, 533, 534, 535, 536), their amplitudes are indeterminate as visible with a positive correlation between the parameters P and D. Fig. 5C illustrates a schematic 540 of HNN simulations in which a single excitatory proximal / distal input with variable synaptic conductances and latencies produce positive (red) / negative (green) deflections in the current dipole.

[0103] As described further in Materials and methods, HNN is a neocortical column model with multiple layers and biophysically detailed local cell types. The local network receives exogenous excitatory synaptic input through layer specific pathways that effectively synapse on the proximal and distal dendrites of the pyramidal neurons, representing “feedforward” and “feedback” input from thalamus and higher order cortical areas. These inputs are simulated with external “spikes” that activate layer specific excitatory synapses (see Figs 2B-2C, and reduced schematic in Fig. 5C). This synaptic activity induces current flow within the pyramidal neuron dendrites, which is summed across the population to simulate a net current dipole that is directly comparable to that measured with MEG / EEG. Several previous studies have shown that patterns of activation of the local network through these pathways can generate commonly measured MEG / EEG current dipole signals such as event related potentials and low frequency brain rhythms.

[0104] To set up an inference problem that is comparable to the RC circuit example, patterns of drive to the network that create subthreshold current flow in the pyramidal neurons are considered, effectively “disconnecting” the network, because local synaptic interactions depend on local cell firing. Simulations with spiking dynamics will be demonstrated in the following section. For simplicity, the subthreshold current flow in the L5 pyramidal cells is first described only (Fig. 5C). Specifically, HNN simulations were run with single exogenous spikes that activate excitatory synapses on the proximal and distal dendrites of L5 pyramidal cells. Synaptic excitation of distal synapses generates current flow down the dendrites (e.g. see green arrow 542 of Fig 5C), and excitation of proximal dendrites generates current flow up the dendrites (e.g. see red arrow 544 of Fig 5C). A delay between these the time of the two driving spikes can create a net current dipole signal that is analogous to that observed in the RC circuit for a non-zero time delay between the applied currents (see Fig 5A, curves 512, 522). Further, when the delay between the spikes iszero (see Fig 5A, curves 514, 524) an indeterminacy in parameter estimation can occur, as described below.

[0105] With this set up, SBI was applied to infer parameters that mimic those of the RC circuit example, using PCA30 as the chosen summary statistic to constrain the inference problem. Namely, the strength of proximal and distal excitatory inputs, referred to as P and D, parameterized as the maximal conductance at their respective synapses and the latency At between the two inputs. The proximal input time was fixed, and the distal input time varied with At. The prior distribution over P and D was set to ensure that all simulations remained subthreshold (see Table 900 of Fig. 9 for prior distribution ranges).

[0106] Similar to the RC circuit example, simulations with At = 0 produce a current dipole with a reduced amplitude (Fig 5A, curves 514, 524), whereas At 0 produces a clear positive and negative peak (Fig 5A, curves 512, 522). As shown in Fig 5B(data 531, 532, 533, 534, 535, 536), when conditioned with At = 0, any posterior distribution involving P and D will exhibit an indeterminacy. This indicates that the proximal and distal inputs can compensate within a small range to produce similar current dipole waveforms. Unlike the RC circuit, this interaction does not span the full range of input strengths, and instead is more tightly concentrated around the ground truth (Fig. 5B(stars 537)).

[0107] The example shown in Figs. 6A-6B demonstrates that the choice of summary statistics impact the learned posterior distribution approximation, and that diagnostics can be used to evaluate the quality of the parameter estimation. Figs. 6A-6D illustrate SBI diagnostics of summary statistics in HNN. In Figs. 6A-6D, the analysis shown in Figs. 4A-4D is repeated on a simplified HNN simulations for comparison. Fig. 6A illustrates summary statistics included PCA30 602, PCA4 604, Peak 606, and BandPower 608. Two exemplar simulations with input latencies At = 0 (waveforms 602a, 604a, 606a, 608a) and At 0 (waveforms 602b, 604b, 606b, 608b) are shown. Fig. 6B shows the approximate posterior when conditioned on both the exemplar waveforms for PCA30 612, PCA4 614, Peak 616, and BandPower 618. The At = 0 time series produces a positive correlation between P and D for all summary statistics.

[0108] When At 0, both PCA30 612 and PCA4 614 produced a posterior that is localized around the ground truth (data 612a, 614a). When At = 0, the posterior was still concentrated but with a slight indeterminacy (data 612b, 614b) that was less prominent than the analogoussimulation in the RC circuit (data 412b, 414b). The Peak summary statistic produced a posterior that is well clustered around the ground truth for At 0 (data 616a), but for At = 0 exhibited a much more striking indeterminacy (data 616b). In contrast, BandPower 618 was insufficient for ground truth recovery for both the At 0 and At = 0 as was the case for the RC circuit simulations. Interestingly, the indeterminacy for the At = 0 waveform with BandPower 618 features is distinct from the RC circuit (Bandpower 418) and exhibits a clear linear interaction.

[0109] In Fig. 6C, local parameter recovery error (PRE) is shown for PCA30 622, PCA4 624, Peak 626, and BandPower 628. Unlike the RC circuit, PCA30 622 and PCA4 624 permit better ground truth recovery even when At is near zero. In contrast, Peak 626 features have poor parameter recovery similar to the RC example. In Fig. 6D, local posterior predictive checks (PPC) are shown for PCA30 632, PCA4 634, Peak 636, and BandPower 638. PCA30 632 and PCA4 634 produce the values across the parameter grid.

[0110] Desirable properties of the posterior estimates were quantified, with the local PRE and PPC analysis described above. At At « 0, PRE values were large when using BandPower 628 features, relatively smaller for PCA4 624 and Peak 626, and almost completely disappears for PCA30 622. The PPC values for the BandPower 638 show that inference using this summary statistic produced results that were highly dissimilar to the conditioning observations, while PCA30 632, PCA4 634, and Peak 636 exhibit a much lower PPC values. Local PRE and PPC analysis largely agrees with the summary statistic performance captured by the RC circuit PPC heatmaps (Figs. 4C-4D).

[0111] Global PRE and PPC values confirm a similar ranking with PCA30 performing best (P - gPRE: 0.005; gPPC: 8.675e-6 nAm) and BandPower exhibiting higher values for both diagnostics (P - gPRE: 0.098; gPPC: 4.250e-5 nAm) (see Data and code availability for results of diagnostics for all parameters and conditions).

[0112] In summary, considerately chosen summary statistics like Peak features can perform well, but leveraging information from the entire time series using PCA30 produced consistently lower PRE and PPC values. The highly effective parameter recovery across the entire parameter grid for PCA30 suggests that SBI permits a near unique mapping from dipole waveform to parameters when the summary statistic accounts for the full time series waveform and the parameters are kept in a subthreshold regime. In the next example, it is shown that this is not truein general. Even when using information from the full dipole waveform with PC A, inference in HNN simulations that include stochasticity can produce substantial parameter indeterminacies.SBI expands previously proposed mechanisms of subthreshold Beta Event simulations in HNN and shows stochastic simulation can lead to indeterminacies

[0113] Figs. 7A-7D illustrate HNN-SBI recovery circuit mechanism of Beta Event amplitude described in previous studies. The HNN modeling framework has been applied to propose novel mechanisms for the cellular and circuit level generation of transient low frequency oscillation in the 15-29 Hz beta-frequency band, referred to as Beta Events or beta bursts. Many studies have shown that Beta Events occur throughout the brain and that they often have a stereotypical waveform shape that resembles an inverted Ricker wavelet lasting ~ 150 ms. Fig. 7A shows a schematic 700 of Beta Event simulations in HNN. Beta events are generated by a simultaneous burst of subthreshold proximal excitatory inputs 702 and distal excitatory inputs 704 to L5 pyramidal neurons 706. Changes in transient beta activity have been associated with sensory processing and motor action, and the Beta Event waveform shape has been specifically associated with motor pathology in Parkinson’s disease and with the aging process. HNN provides potential mechanistic explanations for how changes in waveform shape may emerge. A novel Beta Event mechanism derived from HNN showed that Beta Events can arise from the dendritic integration of coincident bursts of subthreshold proximal and distal dendritic excitatory synaptic inputs to cortical pyramidal neurons, such that the distal drive is effectively stronger and lasts one beta period (~50 ms); a prediction that was supported by invasive laminar recordings in mice and monkeys and high-resolution MEG in humans. More specifically, HNN reproduced Beta Events with the observed stereotypical shape when stochastic trains of action potentials were simulated to drive the proximal and distal dendrites of the pyramidal neurons, nearly simultaneously. Bursts of input whose mean timing and standard deviation where chosen from Gaussian distributions activated excitatory synapses in a proximal and distal connection pattern, as shown in Fig 7A. The inverted Ricker waveform shape depended on the standard deviation of the proximal burst being broader than the distal burst, with the proximal burst occurring over ~ 150 ms and the distal burst occurring over ~50 ms, and the mean time of each burst being the same, i.e. reminiscent of At = 0above. The proximal drive pushes current flow up the pyramidal neuron dendrites, while the distal drive pushes it back down (see Fig. 7A).

[0114] Previous work also showed that the amplitude of the prominent middle trough depended on the variance of the distal drive, such that a parametric lowering of the variance pushed more current flow down the dendrites generating an increased amplitude and sharper peak. This prior study did not perform an automated parameter inference, but rather the results were based on hand tuning and parameter sweeps. Extending these prior results, the SBI methods can be applied to FINN to infer distributions of proximal and distal drive variance that can account for different waveform shapes, and the same diagnostics used above can be used to assess the quality of the estimates. Note that the waveforms analyzed below are simulated Beta Events using parameters from previous HNN studies, and not real recorded data.

[0115] A prior distribution was defined over proximal and distal input variance (Po2, Do2) using the parameters for Beta Event generation as shown in Fig. 9. All other parameters, including the number of spikes, mean input time, and synaptic conductances were all held fixed. The choice emphasizes the need for an a priori hypothesis based on domain expertise to constrain the prior distribution to a tractable subset of all the possible model parameters. PCA30 was used as the summary statistic motivated from the results above. The HNN-SBI workflow was run to obtain a posterior distribution approximation that allows us to infer PG2 and Dc2 for a given waveform (Fig 7C). Example waveforms of real neural activity consisting of large (720) and small (722) amplitude Beta Events, generated with different values of Do2, are shown in Fig. 7B together with the proximal (bottom) and distal (top) spike histograms that generated each waveform. The corresponding estimated posteriors 740 for these examples are shown in Fig 7C. Fig. 7C shows posterior distributions 740 conditioned on large (data 741, 743, 745) and small (data 742, 744, 746) amplitude Beta Events demonstrate that lower distal variance produces a larger amplitude Beta Event. Overlap coefficients (OVE) quantifying the separability of the marginal posterior distributions conditioned on each waveform are shown on the diagonal for the corresponding parameters. The posterior density over the large amplitude Beta Events occupies low values for DG2 in the range of 0-5 ms, while the small amplitude Beta Event produces a posterior density with higher distal variance in the range of 5-10 ms2.

[0116] For both small and large amplitude Beta Events, there is a clear indeterminacy in Po2 . Unlike Do2, the distribution of Po2 is widely spread over the range of the prior distribution, indicating that Po2 cannot be accurately recovered from the waveform. To quantify the separability of these distributions, we can employ the distribution overlap coefficient (OVL) which varies on a scale of [0,1] such that 1 indicates complete overlap, and 0 indicates no overlap. Unsurprisingly, the Po2 distributions for the large and small Beta Events produce an OVL of 0.799 due to the clearly visible high degree of overlap, whereas when comparing the Dc2 distributions they exhibit almost no overlap with an OVL of 2e-10.

[0117] It should be understood that unlike the simulations from the previous sections, where all parameters were deterministic, the exogenous input times considered here are stochastic. As a result, there is not a 1 : 1 mapping from parameters to simulation output. To highlight this, Fig 7B shows simulations with parameters drawn from their respective posterior distributions (black traces). This visually represents the PPC diagnostic, where simulations from each posterior are close to the conditioning observation, but not a perfect match. In this setting, the stochasticity in the simulations make it so that measures such as PPC’s are not guaranteed to be zero even with a perfect approximation of the ground truth posterior distribution. In Fig. 7D, PRE heatmap 760 of DG2 shows accurate parameter recovery when the ground truth parameters of DG2 is small (dark colors at top of heatmap). The pattern observed in the local PRE heatmap 760 of DG2 indicates that this parameter is accurately recovered from simulations generated with a small Do2.

[0118] In summary, the Beta Event example demonstrates how SBI can be used to estimate parameter distributions for given time series waveforms, and to compare potential mechanisms underlying different waveforms. For the hypothetical example comparison shown, the variance of the distal drive could be uniquely inferred while the variance of proximal drive could not. Further, parameter estimates are more accurately recovered for simulations with small distal variance.SBI reveals parameter interdependencies for suprathreshold Event Related Potential simulations in HNN

[0119] In the Beta Event example above, the effective strength of the proximal and distal input were maintained in a range where the activity of the cells remained subthreshold, which naturally limits the dynamic range of the simulation. Next, a more complex example is considered,in which the cells are driven to a suprathreshold spiking regime to show that this additional complexity can lead to parameter estimation indeterminacy that indicates a compensatory interaction between parameters. Importantly, such parameter interactions would not be revealed with other estimation methods such as mean field variational inference given their assumption of Gaussian and independent parameter distributions (i.e., the Laplace approximation).

[0120] Figs. 8A-8E illustrate inferring local connectivity parameters from ERP waveforms. Fig. 8A shows a schematic 800 of ERP simulations in HNN. Evoked activity is driven by a fixed sequence of proximal-distal-proximal exogenous inputs. SBI is used to infer the maximal conductance strength (g) of local excitatory / inhibitory connections to the proximal / distal dendrites of L5 pyramidal neurons for example waveforms. In Fig. 8B, chart 820 shows exemplar simulated ERPs with differing local connectivity strengths chosen from a defined prior distribution (described in the text), along with the fixed timing of the sequence of exogenous inputs for each simulation.

[0121] The example considered describes simulations of a sensory evoked response or event related potential (ERP). HNN has been applied to study source localized ERPs in primary somatosensory cortex from tactile stimuli and in primary auditory cortex from auditory stimuli. In both cases, ERPs were simulated with a sequence of exogenous input that represented an evoked volley of drive to the local circuit that occurs after a sensory stimulus. The drive sequence consisted of a proximal drive representing the initial feedforward input from the thalamus, followed by a distal drive representing feedback input from higher order cortex, followed by a second proximal drive representing a loop of re-emergent thalamocortical drive (see schematic red arrows 802 and green arrows 804 in Fig 8A. These drives were strong enough to generate spiking interactions in the local network and induced current flow up and down the pyramidal neuron dendrites to generate a current dipole ERP waveform analogous to those experimentally recorded (Fig 8A). Note, here we are not examining recorded data, but only example simulations. The specific timing of this exogenous drive sequence for example simulations is shown with arrows 822 and arrows 824 in Fig 8B. The parameters regulating the timing and the strength of these drives were fixed to the same values for the different conditions considered.

[0122] HNN has also been applied to infer neural mechanisms underlying differences in ERP waveform shapes recorded across different experimental conditions. For example, in knownmethods, HNN was applied to infer the circuit mechanisms underlying differences in ERP waveforms emerging from perceptual threshold tactile stimulation (namely brief finger taps that were detected 50% of the time) that were reported as detected (felt) or not-detected (not felt). In that study, earlier and larger amplitude ERP peaks on detected trials could be reproduced with decreases in the timing and increases in the strength of the exogenous inputs. Here, instead of trying to reproduce any empirical findings, SBI is applied to examine the influence of changes in local network connectivity on the ERP waveform as a proof of concept example that examines a small subset of parameters distinct from our prior investigation. The results lay the foundation for future application using empirical data.

[0123] The example begins by simulating example ERPs with different peak amplitudes as show in in Fig. 8B as follows. ERPs were generated using a sequence of exogenous proximal- distal-proximal inputs (Fig. 8A and as described above). The parameters representing the timing and strength of the sequence of exogenous proximal and distal input to the local circuit were chosen to be those distributed with the HNN software representing an example tactile evoked response simulation and fixed to those values (see Table 900 of Fig. 9). A prior distribution was then defined over parameters that define a small subset of the local network excitatory and inhibitory connectivity (see Fig. 2A). These parameters included the maximum conductance (g) of layers 2 and 5 excitatory AMPA (EL2 / EL5) and inhibitory GABAa (IL2 / IL5) connections to the L5 pyramidal cell. Specifically, EL2 and IL2 pertains to synapses on the distal dendrites, EL5 on proximal dendrites, and IL5 on the soma. Note that there exist more local network connections than those varied here as shown in Fig 2A and this chosen prior distribution was not based on any hypothesis or prior knowledge of the impact of local parameters on the ERP, but rather as a tractable example.

[0124] Two example ERPs produced by networks with different local connectivity are shown in Fig. 8B. The ground truth parameters that created these waveforms are shown with stars in Fig 8D. Despite being activated by an identical exogenous input sequence, it is apparent that the local network connectivity differences lead to dramatically distinct current dipole ERP waveforms and corresponding spiking activities. Fig. 8C shows spike raster plots 830, 835 of cell specific firing for the two ERP simulation conditions from Fig. 8B. In Fig. 8C the spiking activity associated with each waveform is largely distinguished by the activity of L5 pyramidal neurons(data 836), with more firing in Condition 2 (waveform 826). For Condition 1 (waveform 828), the first proximal input leads to the beginning of a sustained negative deflection in the current dipole, which persists during the subsequent distal input due to prolonged activation of the L5 basket cells which inhibit the L5 pyramidal neuron soma to pull current flow down the dendrites. Once this inhibition ends, the L5 pyramidal neuron is able to spike, pushing current flow back up the dendrites and the subsequent volley of proximal drive continues to push current flow up and down the pyramidal neuron dendrites due to a similar spiking dynamic. This is in contrast to Condition 2 (waveform 826) in which L5 pyramidal neuron spiking starts almost immediately after the first proximal input and persists, pushing current flow up the dendrites to create a sustained positive deflection in the current dipole that persists through the entire simulation.

[0125] Fig. 8D shows posterior distributions 840 over local connectivity parameters alongside ground truth parameters (stars on diagonal) for conditioning observations. A strong interaction between excitatory / inhibitory distal inputs (EL2 and IL2) is visible in the lower square. Overlap coefficients (OVL) quantifying the separability of the marginal posterior distributions conditioned on each waveform are shown on the diagonal for the corresponding parameters. Fig. 8D shows the results of applying the HNN-SBI framework with the PCA30 summary statistic to estimate the ground truth parameters that generated the ERP waveforms described above. It is apparent that the posterior distributions conditioned on each waveform place high probability mass around the corresponding ground truth parameters (Fig. 8D, stars on diagonal), but also exhibit strong indeterminacies. For example, for each condition, there is a clear interaction between EL2 and IL2, such that as one parameter increases the other also increases, suggesting that these two parameters can compensate one another in a limited range to maintain a constant waveform. We can also observe that between the two conditions there are apparent differences in L5 pyramidal neuron spiking, as well as the sustained negativity (blue) and positivity (orange) observed in the current dipole due to the complex dynamics that each parameter configuration creates.

[0126] To quantify the separability of these distributions, we calculated the OVL coefficient for the marginal distributions of all parameters (Fig. 8D diagonal). Estimated parameter distributions of the synapses on the layer 5 distal dendrites, EL2 and IL2, exhibit a small amount of overlap across conditions with OVL values of 0.190 and 0.011 respectively. The parameter distribution of the synapses on layer 5 somas however were much more distinguishable for the twoconditions, exhibiting OVL values of 2.28e-5 for EL5, and 1.59e-13 for IL5. We additionally performed diagnostic analysis of the local PRE to see how recovery of IL2 changes as a function of IL2 itself, and EL2. Fig. 8E shows local parameter recovery error (PRE) for distal inhibition IL2 indicates errors are higher for observations generated with strong excitatory EL2 and weak inhibitory IL2 distal connections. Fig 8E exhibits a clear pattern indicating that the recovery of IL2 is worse when the network exhibits low IL2 and large EL2.

[0127] In summary, the ERP example provides another demonstration of how SBI can be used to estimate parameter distributions for given time series waveforms and to compare potential mechanisms underlying different waveforms. For the hypothetical example comparison shown, EL5 and IL5 were uniquely inferred with distributions for the waveforms compared exhibiting very little overlap. Distributions for EL2 and IL2 were also separable, albeit with slightly more overlap, and a marked interaction between these two parameters.Discussion

[0128] Recent developments in likelihood-free inference techniques have enabled predictions of parameter distributions for detailed biophysical models like HNN with a level of detail and complexity that was simply not feasible without example methods disclosed herein. Disclosed herein are step-by-step methods to employ SBI in detailed models of neural time series data and to assess the quality and uniqueness of the estimated parameter distributions. First, a simplified RC-circuit example which exemplified the possibility of parameter interactions was described and limitations with chosen summary statistics that do not consider the full time series waveform were highlighted. How distributions of biophysical parameters that can account for a given time series waveform can be inferred using an integrated HNN-SBI workflow applied to two common MEG / EEG signal motifs (Beta Events and ERPs) was demonstrated. Additionally, this disclosure demonstrates how to assess overlap of the distributions from two different waveforms. Examples disclosed herein provide useful examples and methods to highlight critical decisions in the inference pipeline. There are several major takeaways from this disclosure. First, highly nonlinear biophysically detailed neural models like HNN are not suitable for Bayesian estimation methods that require access to a likelihood function (e.g., variational inference) or that approximate posterior distributions with independent Gaussians (i.e., Laplace approximation).Rather they necessitate a method that can estimate complex posterior distributions and parameter interactions from a simulated dataset (e.g. SBI with masked autoregressive flows).

[0129] Second, an important initial step in the SBI process is to identify a limited set of parameters and a range of values for those parameters that are assumed to be able to generate the waveform of interest and variation around it (i.e., the prior distribution). Due to the large-scale nature of biophysically detailed network models, it is not possible to perform SBI on all parameters at once. The choice of the prior distribution represents a hypothesis about parameters that are assumed to contribute to variation in the waveform. This hypothesis can be informed by domain knowledge of the question of interest. In the HNN examples shown, the hypothesized parameters of interest for estimation were the strength of the proximal and distal excitatory synaptic drive for the Beta event simulation, and local excitatory and inhibitory connectivity for the ERP simulation; these parameters were chosen only for didactic purposes. All other parameters were fixed based on previous studies.

[0130] Third, posterior diagnostics like PRE and PPC, are valuable tools to guide decisions in the inference pipeline, e.g. optimal summary statistic selection, and OVL can be used to assess the uniqueness of estimated distribution for two different waveforms. Fourth, when estimating parameters that account for time series waveforms, summary statistics informed by the full time series such as PCA are the most effective at retaining essential information for mapping recorded signals to underlying parameters. While hand-crafted summary statistics, such as peak latency or amplitude, can permit an accurate mapping for certain waveform features, in some cases, their selection may be insufficient to identify unique parameters distributions.Comparison with inference in MEG / EEG neural modeling frameworks that rely on dynamic causal modeling

[0131] While there are several modeling frameworks for simulating MEG / EEG signals, the other frameworks that use likelihood-based Bayesian inference to estimate parameters fall in the category of Dynamic Causal Modeling (DCM). It is important to emphasize that while the HNN-SBI framework conceptually overlaps with DCM, they are two fundamentally distinct techniques which address different questions. At its base, DCM combines variational inference, a computationally efficient Bayesian inference algorithm, with neural mass models, as well as an observation model which translates simulated activity to experimental measures (e.g., MEG / EEG).Neural mass models refer to a specific class of neural models where a single variable represents the aggregate activity of large neural populations (e.g. population spike rates). The inferred parameters in the DCM framework most often represent the coupling strength between distinct neural masses (e.g. population nodes).

[0132] By making simplifying anatomical and physiologic assumptions, neural mass models in the DCM framework can be employed in a large variety of inference problems due to their computational efficiency. However, their ability to make precise biophysical predictions on cellular and local circuit level processes is limited as the parameters are an abstraction representing population level activity. For example, DCM employs the mean-field assumption, meaning that biologically important compensatory interactions such as excitatory / inhibitory (E / I) balance cannot be directly characterized in terms of synaptic conductance. Further, this means that DCM is not capable of representing the parameter indeterminacies with interactions that we showed can occur in the HNN-SBI framework. There are, however, advantages of DCM over the HNN-SBI framework that access to a known likelihood function and other simplifying assumptions allow. For example, one critical question that the HNN-SBI framework is currently not suited to address is inference with multiple spatially separated cortical sources. While theoretically possible, the high computational demands of HNN-SBI severely limit the ability to explore multi-area interactions, and highlight the importance of using neural mass models and DCM in the analysis of whole-brain neuroimaging data. Recent work has shown that neural mass models can also be integrated with the SBI-framework for whole-brain studies, highlighting the adaptability of the SBI methodology to a wide variety of neural models.Comparison to other biophysically detailed neural modeling studies and estimation techniques

[0133] The simulation process disclosed herein extends prior work using SBI to estimate parameters in detailed neural models. Prior work applying SBI to neural models has included a single compartment Hodgkin-Huxley model, and the stomatogastric ganglion model, both of which include an extensive parameter set, but contain significantly less detail and are smaller scale models than HNN. Additionally, non-amortized inference was performed in these models using sequential neural posterior estimation, allowing a much larger parameter set to be inferred, but only for a specific single observation. In contrast, the use of sequential methods was omitted toperform inference on multiple observations using the same trained neural density estimator, but at the expense of the number of parameters that can be inferred simultaneously. The SBI workflow applied here used the neural density estimation technique known as masked autoregressive flows. There are currently a large number of neural density estimation techniques beyond this choice, each offering distinct advantages such as sample efficiency, expressivity, and likelihood evaluation.

[0134] Concerns have been raised about the limits of such tools in the domain of Bayesian inference for scientific discovery. Unfortunately, there currently exist very few techniques for the validation of posterior distributions learned through neural density estimation beyond PPC and PRE diagnostics shown here. One promising work is simulation-based calibration, which plays a similar role as PPC’s by measuring properties the posterior approximation should satisfy if it is close to the ground truth. It is important to note, however, that this technique assesses the quality of the posterior approximation for the marginals of each parameter separately. More research in the domain of multidimensional calibration in the context of likelihood-free inference will be crucial to better represent complex parameter interactions like local E / I balance, as shown in Fig. 8. Nearly all parameter estimation techniques in high-dimensional biophy sically detailed neural modeling will be computationally expensive. Indeed, SBI with HNN has a high computational load (for each example shown here, 100,000 simulations were run on a computing cluster in parallel over 512 CPU cores).

[0135] While the upfront computational costs are high, there are advantages to SBI over other estimation techniques, for example COB YEA estimation, which has also been applied in HNN. The main distinguishing factor is that SBI makes use of every simulation to build an accurate approximation of the posterior distribution for many waveforms. In contrast, COBYLA uses simulations to iteratively search for an optimized parameter set for a single waveform. Once trained, the neural density estimator in the SBI framework can be applied again on new time series waveforms (that fall within the prior distribution) without retraining. As shown in the results, the posterior distribution is an object with several utilities. The mapping between observations and ground truth parameters is emphasized herein, but there are alternative uses such as parameter optimization via non-amortized inference, as well as building a more basic understanding of the model itself. Further, significant research efforts currently underway have the potential to improvethe computational cost of parameter estimation in biophysical neuron models enhancing the accessibility.Other important future directions

[0136] Disclosed herein it is shown that PCA is the appropriate choice when compared to simple hand-crafted summary statistics when performing SBI on time series waveforms. However, PCA is constrained to preserve high-variance features, when in fact low-variance features may also be critical for identifying certain parameters. An important line of future work is the improvement of methods to learn summary statistics from neurophysiological signals that can help identify features of the signal that are essential for accurate parameter estimation. A promising development in this domain is the use of embedding networks that are trained to estimate summary statistics simultaneously with the neural density estimator used to approximate the posterior distribution that can account for those summary statistics. Currently, it is unclear if existing methods to train these embedding networks coupled to neural density estimators are sufficient and require further analysis. The HNN-SB1 example disclosed herein focused on making inferences by constraining only to one output of the model; namely simulated current dipole waveforms. However, due to the multi-scale nature of the HNN model there are many other model outputs that could help constrain the inference problem, such as cell spiking profiles, and / or local field potential signals. A major advantage of the Bayesian framework is the ability to flexibly integrate multiple features into the parameter estimation. If additional multiscale data is known, properties of this data can provide further summary statistics over which the inference problem can be constrained.Discussion Conclusion

[0137] Using detailed neural models in a Bayesian framework is the product of significant developments in machine learning, biophysical modeling, and high-performance computing that have evolved largely independently. Results disclosed herein demonstrate that large-scale biophysically detailed models, like HNN, are now amenable to Bayesian methods via the SBI framework, an approach that has not been feasible in the past. However, this novel combination produces new conceptual and technical challenges that must be addressed to effectively use these techniques. It is apparent that the combination of HNN with SBI is a step forward for making mechanistic inferences underlying MEG / EEG biomarkers, with the potential to provide novelcircuit-level predictions on disease and neural function. These results lay the foundation for similar integration of SBI into the growing number of biophysically detailed neural modeling frameworks to advance neuroscience discovery.Supporting InformationPrior distribution setup and sampling

[0138] Prior samples Oi E Rd were generated by using the PyTorch Uniform distribution on the interval [0,1). The values in each dimension were then linearly mapped to the range of their corresponding parameter values. For parameters specifying the maximum conductance g (nanosiemens, nS) of synaptic connections in HNN, the values were additionally exponentiated in base 10 after being mapped to the appropriate range.

[0139] The prior support and transform function for the parameters of examples shown in the text.Simulation and SBI training

[0140] Prior samples and simulations were all generated and stored in the form of NumPy binary arrays before neural density estimator training. The SBI Python package was used for all neural density estimator training and posterior evaluation. A masked autoregressive flow architecture was utilized for approximation of the posterior distribution. Posteriors for all examples were trained using a dataset of 100,000 samples from the prior distribution. Gaussian white noise was added to training observations xi. The variance of the Gaussian noise added to observations was 0.01 for RC circuit simulations, and le-5 for HNN simulations.

[0141] All analysis was performed on the Expanse supercomputing cluster managed by XSEDE and the Neuroscience Gateway. HNN simulations were generated using the Dask distributed scheduler configured for the SLURM workload manager.

[0142] Diagnostic heatmaps were constructed by defining a grid over the support of the prior with a range of [0.05, 0.95]d, with a resolution of 10 samples in each dimension d.Additional Example

[0143] In the following embodiment, surrogate models of individual cells are implemented using a spiking neural network (SNN) architecture to connect them into a biophysical networkmodel. This construction permits efficient gradient-based optimization of cell-cell connectivity parameters, and the ability to optimize to complex neural activity patterns. The effectiveness of this approach is demonstrated by using surrogate models to approximate a detailed model of the neocortex, the HNN.

[0144] Figs. 10A-10C illustrates a surrogate model 1000 of an individual cell which can accurately predict compartmental voltage responses and net current dipole of individual neurons in HNN. Fig. 10A shows the HNN cortical column model is composed of 4 distinct cell types: L2 / 3 and L5 excitatory pyramidal (e) and inhibitory (i) neurons. Fig. 10B shows a diagram 1030 illustrating how the surrogate model is trained to take synapse activation times as input, and outputs compartmental voltages and current dipole signals. Fig. 10C shows predictions 1040 of the surrogate model (teal) compared to L5e neuron simulations (black).

[0145] Figs. 11A-11E illustrate connectivity parameters optimized by the model of Fig. 10. Fig 11A shows power spectral density 1100 of dipole produced by initial 1102 and final 1104 connectivity parameters. Fig. 11B shows a chart 1100 illustrating the optimizer efficiently maximized high frequency band power. Fig. 11C shows a chart 1120 displaying simulated dipole shows high frequency oscillations 1122 with final connectivity parameters 1124. Fig. 11D shows spike raster 1130 of initial parameters. Fig. 1 IE shows spike raster 1140 of final parameters. Fig. 1 IF shows optimization 1150 of cell-cell connection weights on each epoch.

[0146] As a proof-of-concept, the strength of connectivity among neurons that gives rise to 15-60 Hz oscillations from noisy background drive was inferred (Figs. 11A-11C), with corresponding predictions on cell spiking activity (Figs. 11C-11F). These methods open detailed biophysical modeling to questions that have been previously restricted to more abstract mathematically-tractable models with limited biological interpretability. Important future applications include the study of multi-area network models and cortical traveling waves.

[0147] HNN is a large-scale detailed model of a cortical column designed to simulate electrical currents (i.e., current dipoles) underlying M / EEG signals with interpretability at the cell and circuit level. The 4 principle cell types in the model are L2 / 3 and L5 excitatory pyramidal (e) and inhibitory (i) neurons. A CNN-LSTM (Convolutional Neural Network, Long Short-Term Memory) implemented in PyTorch was trained to predict the response (e.g., voltage) of individual neurons to a train of input spikes. The response of the individual neurons is returned as a timeseries that represents the voltage of the neuron. The training and validation sets were produced by activating all excitatory AMPA synapses with random spike trains of 10 Hz Poisson noise (independently across compartments). The final trained surrogate model accurately predicted compartmental voltages and current dipoles. For L5e, the correlation between the HNN simulation and surrogate model predictions on a held-out validation set was 0.843 and 0.861 for the membrane potentials at the soma and apical tuft, and 0.772 for the current dipole.

[0148] The single-cell surrogate models were then connected in a cortical microcircuit based on the cell-cell connectivity structure defined by HNN, with 50 excitatory and 16 inhibitory cells in each layer. A potential benefit of surrogate modeling is the ability to use optimizers such as gradient-based optimizers to tune the connectivity structure (Adam implemented in PyTorch used below). In other examples, the optimizer may additionally or alternatively include one or more of backpropagation, evolutionary algorithms, and Bayesian optimization. A naive implementation of a surrogate network model cannot be optimized because thresholding functions used to determine spiking are not differentiable. Therefore, techniques were adopted from SNN methods and the spike threshold function was equipped with an approximate gradient (the fast sigmoid function) which is critical for performing backpropagation through the surrogate network model. The major benefit of this approach is the ability to define complex neural activity patterns as the optimization objective, enabling rapid identification of connectivity parameters that produce desired outputs.

[0149] The network’s ability to optimize connectivity parameters was tested by tasking it to maximize spontaneous high frequency oscillations (15-60 Hz) in a network driven by 10 Hz Poisson noise. The loss function was defined as L = -log( R 60 15 SX(co) dco) and effectively maximizes 15-60 Hz power in the dipole signal X(t). The optimizer rapidly increased high frequency power from an initial connectivity configuration that produced low frequency activity in the alpha (8-12 Hz) range to a final configuration that produced substantially more activity in the beta (15-30 Hz) range after 10 epochs. Spike rasters of the initial and final configurations show that both oscillations are due to synchronous firing of L5e cells. Inspection of connectivity modifications reveals changes were made exclusively to L5e targeting projections, reflecting L5 pyramidal cells as the main contributor to M / EEG dipoles. L2i— >L5e connections werestrengthened, whereas L5i^L5e connections were weakened, suggesting that L2 and L5 inhibitory neurons play opposing roles in determining oscillation frequency in this proof-of-concept example.Additional Figures

[0150] Fig. 12 is a diagram 1200 illustrating the functional significance of Beta Events. Beta Events, which are observed in MEG / EEG and LFP, consist of a triphasic waveform in a 50 ms period. Beta events perform a potential role in the somatosensory cortex, where reducing timing and / or frequency of Beta Events can lead to decreased tactile detection. Beta Events are widely observed throughout the cortex, but their role in the motor cortex is not entirely known.

[0151] Fig. 13 is a diagram 1300 illustrating local field potential (LFP) motor cortical Beta Events which are robustly present during reach and grasp in nonhuman primates (NHP). Peak aligned LFP Beta Events were present during a cued grasp instructed delay (CGID) experiment. However, it is not entirely known how Beta Events present during phases of the task.

[0152] Fig. 14 is a diagram 1400 illustrating LFP Beta Events exhibiting amplitude modulation during a reaching task. Fig. 14 shows that in response to a task cue at t = 0 ms, Beta Event amplitude decreased (p < 0.001). However, it is unknown what circuit mechanisms underlie Beta Event generation and modulation in the NHP motor cortex.

[0153] Fig. 15 shows diagrams illustrating how HNN connects macroscale signals to underlying neural activity. Schematic 1502 illustrates an HNN simulations including excitatory proximal / distal input. The simulation parameters of the HNN define proximal / distal dendritic drives. Charts 1504 show results of a simulated Beta Event using HNN. Examples disclosed herein describe methods to predict simulation parameters for HNN or another biological neural model from specific observations.

[0154] Fig. 16 is a diagram 1600 illustrating early proximal and late distal dendritic inputs for the HNN simulation schematic 1502 shown in Fig. 15.

[0155] Fig. 17 is a diagram illustrating how SBI may be used to determine parameter distributions for a biological neural model such as HNN where solution space is non-unique (e.g., there exists multiple simulations / model configurations that replicate observed data). SBI can produce posteriori distributions of parameters likely to reproduce observations.

[0156] Fig. 18 is a diagram illustrating observed LFP Beta Events in voltage as a function of time. Curve 1802 shows small observed Beta Events whereas curve 1804 shows large observed Beta Events. Fig. 19 is a diagram illustrating predicted parameter distributions generated by a combination of SBI and HNN and modulation of Beta Event amplitude based on proximal dendritic inputs. The predicted parameters predict early proximal followed by late distal dendritic inputs. Additionally, regarding modulation of Beta Event amplitude, the predicted parameters show increased variance / strength of proximal dendritic inputs produces larger amplitudes.

[0157] Fig. 20 illustrates a computer system 2000, as may be used for parameter prediction for biological neural network models, according to embodiments of the present disclosure. The computing device 2000 may include at least one processor 2002, a memory 2004, and a communications interface 2010.

[0158] The processor 2002 may be any processing unit capable of performing the operations and procedures described in the present disclosure. In various embodiments, the processor 2002 can represent a single processor, multiple processors, a processor with multiple cores, and combinations thereof.

[0159] The memory 2004 is an apparatus that may be either volatile or non-volatile memory and may include RAM, flash, cache, disk drives, and other computer readable memory storage devices. Although shown as a single entity, the memory 2004 may be divided into different memory storage elements such as RAM and one or more hard disk drives. As used herein, the memory 2004 is an example of a device that includes computer-readable storage media, and is not to be interpreted as transmission media or signals per se.

[0160] As shown, the memory 2004 includes various instructions 2006 that are executable by the processor 2002 to provide an operating system 2012 to manage various features of the computer system 2000. Additionally, the memory 2004 includes an artificial neural network 2008 to provide various functionalities to users of the computer system 2000, which include one or more of the features and functionalities described in the present disclosure.

[0161] The communications interface 2010 facilitates communications between the computer system 2000 and other devices, which may also be computing devices as described in relation to Fig. 20. In various embodiments, the communications interface 2010 includes antennas for wireless communications and various wired communication ports. The computer system 2000may also include or be in communication, via the communications interface 2010, one or more input devices (e g., a keyboard, mouse, pen, touch input device, etc.) and one or more output devices (e.g., a display, speakers, a printer, etc.).

[0162] Although not explicitly shown in Fig. 20, it should be recognized that the computer system 2000 may be connected to one or more public and / or private networks via appropriate network connections via the communications interface 2010. It will also be recognized that software instructions may also be loaded into a non-transitory computer readable medium, such as the memory 2004, from an appropriate storage medium or via wired or wireless means.

[0163] Accordingly, the computer system 2000 is an example of a system that includes a processor 2002 and a memory 2004 that includes instructions that (when executed by the processor 2002) perform various embodiments of the present disclosure. Similarly, the memory 2004 is an apparatus that includes instructions that, when executed by a processor 2002, perform various embodiments of the present disclosure.Conclusion

[0164] It should be understood that various changes and modifications to the presently preferred embodiments described herein will be apparent to those skilled in the art. Such changes and modifications can be made without departing from the spirit and scope of the present subject matter and without diminishing its intended advantages. It is therefore intended that such changes and modifications be covered by the appended claims.

Claims

CLAIMSThe invention is claimed as follows:

1. A method for estimating parameters for a biological neural network model, the method comprising: identifying model parameters for the biological neural network model and parameter ranges corresponding to the model parameters; sampling from the parameter ranges to generate a dataset of simulated neural patterns; selecting summary statistics for the dataset of simulated neural patterns; training an artificial neural network based on the summary statistics; and applying the artificial neural network to a neural activity pattern to estimate parameter distributions for simulated neural patterns within the dataset that match the neural activity pattern.

2. The method of claim 1, further comprising determining diagnostics for the estimated parameter distributions for the neural activity pattern.

3. The method of claim 2, wherein the diagnostics include one or more of a parameter recovery error (PRE), a posterior predictive check (PPC), and an overlap coefficient (OVL).

4. The method of claim 1, wherein the biological neural model is one or more of, a biophysically detailed network model, the Human Neocortical Neurosolver (HNN), a spiking neural network, or a neural dynamics model.

5. The method of claim 1, wherein the artificial neural network estimates the parameter distributions using one or more of normalizing flows, variational autoencoders, autoregressive models, and energy-based models.

6. The method of claim 1, wherein the neural activity pattern comprises electrophysiological measurements, including one or more of local field potentials (LFP),electroencephalography (EEG), magnetoencephalography (MEG), single-unit spiking, multi-unit- spiking, and current source density.

7. The method of claim 1, wherein the neural activity pattern comprises blood- oxygenation-level-dependent (BOLD) signals.

8. The method of claim 1, wherein the neural activity pattern comprises optical signals including one or more of calcium imaging and voltage sensors.

9. The method of claim 1, wherein selecting the summary statistics includes mathematically transforming the simulated neural patterns into a latent-space representation.

10. The method of claim 1, wherein the neural activity pattern is one of a real neural activity pattern or a simulated neural activity pattern.

11. A method for approximating simulations of a biological neural network model, the method comprising: training an artificial neural network to approximate a surrogate model, the surrogate model corresponding to an individual cell of the biological neural network model; connecting one or more duplicates of the surrogate model into a surrogate network model; and optimizing the surrogate network model to fit simulations to target activity patterns.

12. The method of claim 11, wherein the biological neural network model is one or more of, a biophysically detailed network model, the Human Neocorti cal Neurosolver (HNN), a spiking neural network, or a neural dynamics model.

13. The method of claim 11, wherein the artificial neural network uses one or more of recurrent neural networks, convolutional neural networks, long short-term memory networks, gated recurrent unit networks, transformers, or hybrid architectures for time series generation.

14. The method of claim 11, wherein the optimizing of the surrogate network model includes one or more of gradient-based methods, backpropagation, evolutionary algorithms, and Bayesian optimization.

15. The method of claim 11, wherein the optimizing of the surrogate network model includes: identifying model parameters for the biological neural network model and parameter ranges corresponding to the model parameters; sampling from the parameter ranges to generate a dataset of simulated neural patterns; selecting summary statistics for the dataset of simulated neural patterns; training an artificial neural network based on the summary statistics; and applying the artificial neural network to a neural activity pattern to estimate parameter distributions for simulated neural patterns within the dataset that match the neural activity pattern.

16. The method of claim 11, further including combining the surrogate model and an output of the optimization into a single artificial neural network model.

Images