# Achieving Occam’s Razor: Deep Learning for Optimal Model Reduction

Botond B. Antal  Department of Biomedical Engineering, State University of New
York at Stony Brook, Stony Brook, NY, USA  Anthony G. Chesebro  Department of
Biomedical Engineering, State University of New York at Stony Brook, Stony
Brook, NY, USA  Helmut H. Strey  Department of Biomedical Engineering, State
University of New York at Stony Brook, Stony Brook, NY, USA  Laufer Center for
Physical and Quantitative Biology, State University of New York at Stony
Brook, Stony Brook, NY, USA  Lilianne R. Mujica-Parodi  Department of
Biomedical Engineering, State University of New York at Stony Brook, Stony
Brook, NY, USA  Laufer Center for Physical and Quantitative Biology, Stony
Brook University, Stony Brook, NY, USA  Corey Weistuch<EMAIL_ADDRESS>Department of Medical Physics, Memorial Sloan Kettering Cancer Center, New
York, NY, USA

###### Abstract

All fields of science depend on mathematical models. Occam’s razor refers to
the principle that good models should exclude parameters beyond those
minimally required to describe the systems they represent. This is because
redundancy can lead to incorrect estimates of model parameters from data, and
thus inaccurate or ambiguous conclusions. Here, we show how deep learning can
be powerfully leveraged to address Occam’s razor. FixFit uses a feedforward
deep neural network with a bottleneck layer to characterize and predict the
behavior of a given model from its input parameters. FixFit has three major
benefits. First, it provides a metric to quantify the original model’s degree
of complexity. Second, it allows for the unique fitting of data. Third, it
provides an unbiased way to discriminate between experimental hypotheses that
add value versus those that do not. In two use cases, we demonstrate the broad
applicability of this method across scientific domains. To validate the method
using a known system, we apply FixFit to recover known composite parameters
for the Kepler orbit model. To illustrate how the method can be applied to
less well-established fields, we use it to identify parameters for a multi-
scale brain model and reduce the search space for viable candidate mechanisms.

##  I Introduction

Mathematical models are commonly used to describe the dynamical behavior of
physical systems. Yet model parameters are not mere mathematical descriptions.
The accurate estimation of parameter values can often yield deep mechanistic
insight, whether with respect to the properties of particles [1], interactions
among genetic networks [2], or the generation of neuronal signaling [3].

A fundamental challenge in parameter fitting and the construction of models
stems from parameter redundancies[4, 5]. Parameter degeneracy is a
particularly problematic in multi-scale models, where emergent measured values
exist many layers above their mechanistic parameters. This can lead to many
different combinations of parameters fitting the observed data equally well, a
phenomenon described as overdetermined or “sloppy models” [6, 7, 8]. Finding
these non-unique solutions is also difficult, as parameter interactions can
introduce numerous local minima that hinder algorithms for data fitting [9].
Consequently, there is often a trade-off between using an easier-to-interpret
model with fewer parameters that fails to describe the system accurately and
using a highly-detailed model that risks redundancy of its parameters [10].

A solution to this trade-off is to identify and account for parameter
dependencies. One of the most widely-used tools to quantify interactions among
parameters is the Fisher Information Matrix, which can be acquired through a
nonlinear least-square Levenberg-Marquardt algorithm. This procedure, for
given data, finds the locally best-fitting parameter values and their
covariance matrix [11, 8]. In overdetermined models, these covariances are
strong for pairs of parameters that are not separable. While the Fisher
Information Matrix only characterizes linear interactions, other methods can
also uncover nonlinear interactions [7, 12]. For example, one can determine
parameter interactions by approximating the model’s behavior as a function of
its parameters using the Taylor series expansion [13]. Nevertheless, a common
limitation of these methods is that they do not provide a functional form for
the redundancies, and therefore cannot guide the parameter fitting process. A
sensible strategy is to combine the redundant parameters into composite
parameters with a unique best fit. Although suitable methods exist, they rely
on studying the model analytically or through its local behaviors and
derivatives [14], a strategy that would not be feasible for highly complex
models. Thus, there remains a need for a more general approach.

Here we show how one may identify and estimate the largest set of lower-
dimensional latent parameters uniquely resolved by the data (Figure 1). The
framework, FixFit, builds on previously established methods and consists of
three steps. First, we find these latent parameters using a neural network
with a bottleneck [15]. We optimize the same architecture at variable
bottleneck widths to identify the optimal latent dimension [16, 17, 18]. The
ability of neural networks to approximate any function allows our method to be
applied to models of arbitrary complexity [19, 20, 21]. Next, FixFit
establishes the relationship between latent and original parameters using
global sensitivity analysis [22]. Finally, we fit latent parameters to data
using global optimization techniques [23].

Here we provide two use-cases for FixFit. To establish validity against a
known system, we first demonstrate its ability to recover nonlinear parameter
combinations for the Kepler orbit model [24]. To demonstrate its potential for
scientific discovery, we then use FixFit to identify previously undiscovered
parameter redundancies in the multi-scale Larter-Breakspear neural mass model
[25, 26].

Figure 1: FixFit compresses interacting parameters into a latent
representation that can be uniquely inferred from data. A: A schematic
representing the goodness of fit landscape of a model with two interacting
parameters. These interactions cause multiple parameter combinations to fit
experimental data equally (the redundant maxima in red). B: The same landscape
but with the two interacting parameters first combined into a single latent
variable. In contrast to the native parameters, numerical fitting over latent
variables will converge to a unique solution. C: FixFit generates such unique
latent representations using a neural network with an encoder, bottleneck
layer, and a decoder. After determining the optimal number ($k$) of latent
(bottleneck) nodes (see Methods), the neural network is trained on pairs of
parameters and their corresponding outputs from computational simulations of a
model of interest. Following training, the bottleneck layer will include a
representation of input parameters that is uniquely inferable from output
data. With the latent representation established, the decoder section of the
neural network can be combined with an optimizer to infer unique parameters in
latent space from previously unseen output data. In addition, the encoder part
can be combined with sensitivity analysis to determine the influence of input
parameters on the latent representation. This enables us to characterize
changes in different output samples in terms of underlying parameters that
would not otherwise be accessible.

##  II Methods

###  II.1 Models

####  II.1.1 Kepler orbit model

The classical Kepler orbit model captures the gravitational motion of two
orbiting planetary bodies through four input parameters (Figure 2A) [24]: the
mass of the orbiting body ($m_{1}$), the mass of the body fixed at the focus
point ($m_{2}$), the closest distance between the two bodies ($r_{0}$) and the
corresponding initial angular velocity ($\omega_{0}$). To facilitate later
training of a neural network, the units of each parameter were chosen such
that the parameter values had comparable magnitudes. Thus, $m_{1}$ and $m_{2}$
were expressed in $kg$, $r_{0}$ in $m$, and $w_{0}$ in $day^{-1}$. The center
of mass of the orbiting body evolves according to:

$m_{1}\frac{d^{2}r}{dt^{2}}-m_{1}r\omega^{2}=\frac{Gm_{1}m_{2}}{r^{2}}$ (1)

where $G$ is the universal constant of gravitation (approximately
$0.5m^{3}kg^{-1}day^{-2}$), $r$ is the distance between the two bodies, and
$\omega$ is the angular velocity $\frac{d\theta}{dt}$. As is apparent, the
mass of the orbiting body ($m_{1}$) cancels in the above equation and thus
cannot be resolved from data. More broadly, the dynamics are highly
constrained by the conservation of energy and angular momentum. As a result,
the orbit $r$ can be readily solved as a function of the angle
$\theta\in[0,2\pi)$ swept by the orbiting object (starting at $\theta=0$)
[24]. These solved orbits are simple ellipses and thus can be described using
only two pieces of information, the eccentricity $e$ (the shape) and the semi-
latus rectum $l$ (the size), as follows:

$r(\theta)=\frac{l}{1+e\cos(\theta)}$ (2)

with

$e=\bigg{|}\frac{r_{0}^{3}\omega_{0}^{2}}{Gm_{2}}-1\bigg{|}$ (3)

$l=\frac{r_{0}^{4}\omega_{0}^{2}}{Gm_{2}}$ (4)

The above two equations establish a ground truth for the degeneracy that we
aimed to characterize and overcome during parameter inference. Specifically,
parameter $m_{1}$ is canceled out and is therefore a completely redundant
parameter. The remaining three parameters ($m_{2}$, $r_{0}$, $\omega_{0}$) can
be compressed into two terms and still uniquely map to outputs [24].

####  II.1.2 Multi-scale brain model

The Larter-Breakspear model (Figure 3) describes the dynamics of a group of
coupled brain regions, each modeled as a population of neurons [25, 26, 27,
28]. Each brain region captures the averaged synaptic processes and voltage-
dependent ion transport ($\text{K}^{+}$, $\text{Na}^{+}$, and
$\text{Ca}^{2+}$) of its constituent neurons. These effective dynamics are
described through three state variables: mean excitatory membrane voltage
V(t), mean inhibitory membrane voltage Z(t), and the proportion of open
potassium channels W(t). Note that although $V$, $Z$, and $W$ depend on time,
this notation is omitted in the following equations for clarity. These states
of each $i$th brain region evolve according to:

$\frac{dV_{i}}{dt}=\\\
-(g_{\text{Ca}}+r_{\text{NMDA}}a_{ee}[(1-c)Q_{V}+cQ_{i}^{\text{network}}])m_{\text{Ca}}(V_{i}-V_{\text{Ca}})\\\
-(g_{\text{Na}}m_{\text{Na}}+a_{ee}[(1-c)Q_{V}+cQ_{i}^{\text{network}}])(V_{i}-V_{\text{Na}})\\\
-g_{\text{K}}W_{i}(V_{i}-V_{\text{K}})-g_{L}(V_{i}-V_{L})-a_{ie}Z_{i}Q_{z}+a_{ne}I_{0}$
(5)

$\frac{dZ_{i}}{dt}=b(a_{ni}I_{0}+a_{ei}V_{i}Q_{V})$ (6)

$\frac{dW_{i}}{dt}=\phi\frac{m_{\text{K}}-W_{i}}{\tau_{\text{K}}}$ (7)

Where $Q_{V}$/$Q_{Z}$ are excitatory/inhibitory mean firing rates, and
$m_{\text{Na}}$, $m_{\text{K}}$, $m_{\text{Ca}}$ are ion channel gating
functions. These are computed as:

$Q_{V}=0.5Q_{V_{\text{max}}}(1+\tanh(\frac{V-V_{T}}{\delta}))$ (8)

$Q_{Z}=0.5Q_{Z_{\text{max}}}(1+\tanh(\frac{V-V_{Z}}{\delta}))$ (9)

$m_{\text{ion}}=0.5(1+\tanh(\frac{V-T_{\text{ion}}}{\delta_{\text{ion}}}))$
(10)

Brain regions are connected through the coupling term
$Q_{i}^{\text{network}}$, which is scaled with a global coupling constant $c$,
to produce whole-brain-scale dynamics. $Q_{i}^{\text{network}}$ is given by:

$Q_{i}^{\text{network}}=\frac{\sum_{j}u_{i,j}Q_{V_{j}}}{\sum u_{i,j}}$ (11)

We considered 78 brain regions selected from the Desikan-Killiany atlas
included in FreeSurfer, with inter-regional structural connectivity
($u_{i,j}$) that was determined from diffusion tensor imaging (DTI) data from
13 healthy human adults [29, 27]. All parameters of the Larter-Breakspear
model are described in Table 1.

Table 1: Parameters of the Larter-Breakspear multi-scale brain model Parameter |  Description
---|---
$V_{\text{Na}}$ |  $\text{Na}^{+}$ reversal potential
$V_{\text{K}}$ |  $\text{K}^{+}$ reversal potential
$V_{\text{Ca}}$ |  $\text{Ca}^{2+}$ reversal potential
$V_{\text{L}}$ |  Leak channels reversal potential
$g_{\text{Na}}$ |  $\text{Na}^{+}$ conductance
$g_{\text{K}}$ |  $\text{K}^{+}$ conductance
$g_{\text{Ca}}$ |  $\text{Ca}^{2+}$ conductance
$g_{\text{L}}$ |  Leak channels conductance
$T_{\text{Na}}$ |  $\text{Na}^{+}$ channel threshold
$T_{\text{K}}$ |  $\text{K}^{+}$ channel threshold
$T_{\text{Ca}}$ |  $\text{Ca}^{2+}$ channel threshold
$\delta_{\text{Na}}$ |  $\text{Na}^{+}$ channel threshold variance
$\delta_{\text{K}}$ |  $\text{K}^{+}$ channel threshold variance
$\delta_{\text{Ca}}$ |  $\text{Ca}^{2+}$ channel threshold variance
$V_{T}$ |  Excitatory neuron threshold voltage
$Z_{T}$ |  Inhibitory neuron threshold voltage
$\delta$ |  Variance of thresholds
$Q_{V_{\text{max}}}$ |  Excitatory population maximum firing rate
$Q_{Z_{\text{max}}}$ |  Inhibitory population maximum firing rate
$a_{ee}$ |  Excitatory-to-excitatory synaptic strength
$a_{ei}$ |  Excitatory-to-inhibitory synaptic strength
$a_{ie}$ |  Inhibitory-to-excitatory synaptic strength
$a_{ne}$ |  Non-specific-to-excitatory synaptic strength
$a_{ni}$ |  Non-specific-to-inhibitory synaptic strength
$I_{0}$ |  Subcortical excitatory input
$b$ |  Time scaling factor
$\phi$ |  Temperature scaling factor
$\tau_{\text{K}}$ |  $\text{K}^{+}$ relaxation time
$r_{\text{NMDA}}$ |  NMDA/AMPA receptor ratio
$c$ |  Global region-to-region coupling constant

Next, to produce a signal compatible with functional MRI, a commonly utilized
neuroimaging modality, the simulated region-specific excitatory signals were
transformed into blood-oxygen-level-dependent (BOLD) signal via the Balloon-
Windkessel model with standard parameters [27, 30]. Finally, we derived
functional connectivity (FC) from a simulated BOLD signal to yield a phase-
invariant signal. We quantified FC with all-to-all Pearson correlation
coefficients among the 78 regions, resulting in a 78-by-78 correlation matrix
for each simulation.

###  II.2 Data generation

####  II.2.1 Kepler orbit model

Training and validation data for the neural network were generated using
equations 2 \- 4. The four input parameters, $m_{1}$, $m_{2}$, $r_{0}$, and
$\omega_{0}$, were sampled with a four dimensional Sobol sequence (SciPy
v1.7.1 [23]). All four parameters were drawn from a range of $[0.1,1]$ (Figure
S1 [31]). Next, a subset of samples was rejected based on an eccentricity
criterion. A raw parameter set was discarded if $e>0.95$ or $e<0.7$. These two
conditions ensured all resultant orbits were ellipse-shaped with moderate
eccentricity. The final sample sizes were 2,276 for training and 253 for
validation. Output space consisted of $r(\theta)$ values computed at 100
evenly distributed $\theta$ values within the range of $[0,2\pi]$ for each
input parameter set. The computed $r(\theta)$ values were log-transformed to
narrow down their range, thus ensuring a favorable output space for the neural
network (Figure S2 [31]).

####  II.2.2 Multi-scale brain model

We first applied domain knowledge to reduce the 30 parameters of the Larter-
Breakspear model to a smaller subset of parameters that is more tractable for
parameter inference. We fixed 19 parameters that we determined to be unlikely
to be sensitive to different biological conditions to default values. All
default values were taken from [27]. The remaining eleven input parameters of
the model were drawn from biologically relevant ranges (Table 2) using a Sobol
sequence. The model was simulated with input parameters still on their
original scales. Later, all eleven input parameters were scaled to within a
range of $[0,1]$ for training the neural network (Figure S3 [31]). Simulations
were performed using Neuroblox.jl, a Julia library optimized for high-
performance computing of dynamical brain circuit models
(http://www.neuroblox.org); while the Larter-Breakspear model is technically
unitless, the parameters are scaled so that a single timestep is 1 ms. [26].
Simulated time series were converted to BOLD using the Balloon-Windkessel
model [30] $(T_{R}=0.8\text{s})$ [32] and then bandpass filtered
($0.01<f<0.1\>Hz$) [33] to quantify functional connectivity [34]. We computed
functional connectivity among the 78 brain regions from the processed time
series and retained values above the diagonal to discard duplicates. The
resultant 3,003 Pearson correlation values per sample constituted the output
space for the neural network (Figure S4 [31]). Simulated data were subject to
two exclusion criteria to ensure biologically realistic behavior: time series
with non-oscillatory behavior or with a mean FC larger than $0.3$ were
discarded. As a result, there were 4,730 retained samples for training and 526
for validation.

Table 2: Investigated parameter subset of the Larter-Breakspear multi-scale brain model Parameter |  Description |  Range
---|---|---
c |  Global region-to-region coupling constant |  $[0.2,0.5]$
$\delta$ |  Variance of thresholds |  $[0.64,0.7]$
$g_{\text{Ca}}$ |  $\text{Ca}^{2+}$ conductance |  $[0.95,1.05]$
$V_{\text{Ca}}$ |  $\text{Ca}^{2+}$ reversal potential |  $[0.95,1.01]$
$g_{\text{K}}$ |  $\text{K}^{+}$ conductance |  $[1.95,2.05]$
$V_{\text{K}}$ |  $\text{K}^{+}$ reversal potential |  $[-0.75,-0.65]$
$g_{\text{Na}}$ |  $\text{Na}^{+}$ conductance |  $[6.6,6.8]$
$V_{\text{Na}}$ |  $\text{Na}^{+}$ reversal potential |  $[0.48,0.58]$
$a_{ee}$ |  Excitatory-to-excitatory synaptic strength |  $[0.33,0.39]$
$a_{ei}$ |  Excitatory-to-inhibitory synaptic strength |  $[1.95,2.05]$
$r_{\text{NMDA}}$ |  NMDA/AMPA receptor ratio |  $[0.20,0.30]$

###  II.3 Neural network

####  II.3.1 Kepler orbit model

We utilized a network with a fully connected architecture and a bottleneck
layer in the middle to enforce compression. The network structure was
determined empirically based on a trade-off between model accuracy and
training speed. There were two hidden layers before and two after the
bottleneck layer, all of which had a $\tanh$ activation function. The final
output layer, by contrast, was given a linear activation function to map to
the output space. By the universal approximation theory of neural networks,
hidden layers before and after the bottleneck layer had 14 and 110 nodes [35,
36]. The neural network was implemented with TensorFlow (version 2.6.0 [37]).
For exact details of the applied architecture, see Figure S5 [31].

####  II.3.2 Multi-scale brain model

A similar architecture was used for the brain model with minor adaptations to
a significantly wider output space (3,003 values per sample). The encoder
included two fully connected hidden layers with 21 nodes in each, whereas the
decoder had a single hidden layer with 3,013 nodes. All hidden layers were
given $ReLu$ activation, whereas the bottleneck and output layers were both
implemented with a linear activation function. Details of the network are
described in Figure S6 [31].

###  II.4 Bottleneck analysis

To determine the number of uniquely resolvable latent parameters for a given
model, we trained our neural network architecture on the previously described
training data at varying bottleneck layer dimensions ($k\in\\{1,2,3,4,5\\}$).
At each $k$ increment, ten replicate training runs were performed
independently. We employed 5,000 epochs during training to ensure model
accuracy was not limited by training length. Both analyses involved batches of
256 samples. Model weights were updated by Adam optimizer [38] using mean
squared error as the metric to optimize. We employed early stopping during
optimization; we stopped training if accuracy had not improved for 200
subsequent epochs. From each replicate, we extracted the minimal validation
error and the corresponding model weights. These were the outputs considered
during subsequent analyses. For $k$ dimensions that were equal to or greater
than the underlying complexity, validation error was expected to not decrease
further with $k$. To distinguish the ideal $k$ from overparameterized
solutions, we chose the smallest $k$ for which the error was not statistically
significant from the minimum error across all $k$ values.

###  II.5 Global sensitivity analysis

Following the acquisition of a latent representation, we determined the
influence of the original parameters on the latent parameters using global
sensitivity analysis. We used the encoder (all layers before the bottleneck of
the optimized neural network) to compute data pairs of input parameters and
corresponding latent parameters. Since these pairs were unevenly sampled due
to our filtering steps, we employed structural and correlative sensitivity
analysis (SCSA), a method that handles non-uniform sampling and accounts for
correlations among inputs, to compute global sensitivities [22]. SCSA
partitions sensitivity into uncorrelated and correlated contributions. Of
these, we used the uncorrelated component $S^{unc}_{ij}$, reflecting the
exclusive contribution of each input parameter $x_{i}$, to more sparsely
identify drivers of each latent variable $y_{j}$:

$S_{i,j}^{unc}=\sum_{s=1}^{N}(f_{p_{i,j}}(x_{i}^{(s)}))^{2}/\sum_{s=1}^{N}(y_{j}^{(s)}-\bar{y_{j}})^{2}$
(12)

Where $s$ is the sample index, $N$ is the total number of samples, and
$f_{p_{i,j}}$ is a data-driven sensitivity function for $x_{i}$ with respect
to $y_{j}$ [22]. We used SALib’s (v1.4.5) implementation of SCSA [39, 40] to
determine sensitivities. The above procedure was applied to every $j$-th
latent parameter separately to derive each row of the $I\times J$ sensitivity
matrix for $I$ input and $J$ latent parameters.

###  II.6 Global fitting

We employed global optimization to infer native and latent parameters from
output data. In both cases, parameters were first normalized to a hypercube
($p_{i}\in[0,1]$). For latent parameters, we determined bounds for the
hypercube based on ranges that we observed for the latent parameters. Next, we
used a basin-hopping algorithm (SciPy v1.7.1 [23]) to find the best-fitting
parameters. Basin-hopping combines a global step-taking routine and a local
optimizer to find a global minimum across the parameter space. The global
step-taking routine was initialized at 0.5 within the hypercube, and each step
involved a random displacement of coordinates with a step size of 0.2. Local
minima were found at each step, including the initial point, using the
Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. The goodness of fit was
evaluated based on the residual sum of squares. Steps were accepted with a
probability $P$ determined by the change in the value of the cost function
$f(p)$.

$P=exp(-\frac{f(x_{new})-f(x_{old})}{T})$ (13)

Notice that the acceptance rate was 100% for steps that improved the objective
but still nonzero for steps that yielded worse objectives. This acceptance
rate is tuned through a temperature parameter $T$ individually adjusted for
each model. This allows the algorithm to explore the landscape within the
hypercube and thus greatly increases the likelihood of finding the global
optimum.

##  III Results

Figure 2: Recovering the known parameter redundancies of the Kepler orbit
model. A: Diagram of an example Keplerian orbit, the corresponding input
parameters (red), and the model outputs ($\theta$, $r$ pairs). B: Validation
error of the compressed representation found by FixFit as a function of the
bottleneck dimension ($k$). Multiple replicates were performed at each $k$
using a stochastic optimizer. Shown for each $k$ are the individual data
points, mean, and standard error. FixFit correctly identified the underlying
redundancy of the Kepler orbit model by indicating saturating error at $k=2$
(see Methods). In subsequent panels, we utilized one of the fitted neural
network replicates at $k=2$. C: Objective landscape of two of the original
four dimensions. Each point corresponds to an optimizer evaluation and the
corresponding objective value (log-transformed sum of squares error) with
respect to two parameters, $m_{2}$ and $\omega_{0}$ (scaled between $0$ and
$1$). The broad distribution of low error evaluations suggests a parameter
redundancy. Consequently, standard fitting could not reliably identify the
underlying ground truth parameters (marked by dashed lines). D: Objective
landscape of the two latent dimensions. In contrast to the previous case, the
same optimizer procedure converged to the correct, unique minimum in the
latent parameter space ($L_{1}$,$L_{2}$) identified by FixFit. E: SCSA global
sensitivities of the latent parameters to the four original parameters. Higher
values of sensitivity (green) indicate a stronger influence. Considering the
closed-form solution as a reference, SCSA correctly identified that parameter
$m_{1}$ had no influence on outputs and that the remaining three parameters
$m_{2}$, $r_{0}$, and $\omega_{0}$ together determined the two latent
parameters.

###  III.1 Recovering the known parameter redundancies of the Kepler orbit
model

The Kepler model [24] describes the elliptical orbit (pairs of angles
($\theta$) and radii ($r$), see Figure 2A) of two gravitationally-attracting
bodies as a function of four input parameters ($m_{1}$, $m_{2}$, $r_{0}$,
$\omega_{0}$) (Methods). Ellipses, however, can be entirely described by two
composite shape parameters, eccentricity ($e$) and the semi-latus rectum
($l$). As the dependencies between the four input and two composite parameters
are known analytically, our results can be compared against a ground truth.

After generating simulated data using the original model, we trained a neural
network with a variable-width bottleneck layer ($k=1...4$, (Figure S7 [31]))
to determine, from data, the underlying complexity of the Kepler model (see
Methods). The validation error saturated at $k=2$, suggesting that the
underlying minimal representation, consistent with the ground truth, is two-
dimensional (Figure 2B). In subsequent analyses, we thus utilized a trained
neural network from one of the replicates at $k=2$ to encode a minimal
representation of the four input parameters.

Figure 3: The Larter-Breakspear brain model connects microscopic properties of
neurons with patterns of brain-wide activity. This model produces brain
dynamics by considering neuronal properties on multiple spatial scales
(bottom). On the bottom scale, it takes into account transmembrane ion
transport (bottom left), which is the basis of neuronal dynamics. Ion
transport is facilitated by voltage-gated ion channels that are separately
specified for $\text{K}^{+}$, $\text{Na}^{+}$, and $\text{Ca}^{2+}$ ions.
Next, a mean-field approximation aggregates single-cell behavior into the
population level of neural masses. An interacting pair of an excitatory (+)
and an inhibitory (-) population is integratively modeled to form individual
brain regions. Finally, on the whole-brain scale, a network of 78 connected
brain regions (connected based on diffusion tensor imaging) is considered to
complete the spatial span of the model (bottom right). The simulated region-
specific brain dynamics are then transformed to achieve compatibility with
functional MRI, a non-invasive neuroimaging modality (top). Simulated time
series (top right) are first converted into blood-oxygen-level dependent
(BOLD) signals using a hemodynamic response function. Then, to remove phase-
specific information, all-to-all Pearson correlations are computed among the
time series resulting in a 78-by-78 functional connectivity (FC) matrix (top
left).

To demonstrate the practical issue of inferring degenerate model parameters,
we attempted to uniquely infer all four original parameters from simulated
Kepler orbits. As expected, the global optimizer found a wide set of equally
optimal solutions across the parameter space (Figure 2C). By contrast, global
optimization in latent parameter space (Figure S8 [31]) (through the decoder
section of the neural network, see Methods) recovered a unique solution that
coincided with the expected ground truth (Figure 2D).

Next, we demonstrate that, by employing Structural and Correlative Sensitivity
Analysis (SCSA), FixFit can quantify relationships between the original and
latent model parameters (see Methods). Applied to the Kepler model (Figure
2E), the approach provides global sensitivities consistent with closed-form
solutions for the two intermediate terms $e$ and $l$. Firstly, as expected, no
sensitivity was detected with respect to $m_{1}$. As such, $m_{1}$ was a
completely redundant parameter. The remaining three parameters, $r_{0}$,
$m_{2}$ and $\omega_{0}$ together determined $L_{1}$. Finally, $r_{0}$ alone
influenced $L_{2}$. Both ground truth terms, $e$ and $l$, by contrast,
included all three of these parameters (equations 3 and 4). However, these
analytical expressions share a common term with respect to the three input
parameters and differ only by an additional power of $r_{0}$ in $l$. By
exploiting this relationship, FixFit separates the influence of $r_{0}$ from
that of $m_{2}$ and $\omega_{0}$, thus providing a sparser latent
representation compared to the ground truth.

###  III.2 Identifying novel parameter relationships in a multi-scale brain
model

The Larter-Breakspear model is commonly used to connect microscopic neuronal
properties with emergent brain activity, such as that measured using
functional magnetic resonance imaging (fMRI) [41] (Figure 3). This is achieved
by modeling the voltage-gated ion dynamics (here, $\text{Na}^{+}$,
$\text{K}^{+}$, and $\text{Ca}^{2+}$) of a population of neurons called a
“neural mass” [25, 26, 28]. Previous studies have further shown that, by
coupling multiple neural masses according to the inter-regional connectivities
measured through Diffusion Tensor Imaging (DTI), one can simulate semi-
realistic fMRI dynamics [27]. As a challenge, this model has many parameters
such that even after assigning identical parameter values to all 78 of the
resulting brain regions and fixing biologically inert parameters, the model
still has eleven remaining parameters corresponding to coupled biological
processes (see Methods). Furthermore, these simulated fMRI activities are
typically studied using functional connectivity (FC) or the matrix of region-
to-region Pearson correlations [34], thus resulting in a reduction of
information. To address this issue, we applied FixFit to resolve potential
redundancies and, more broadly, characterize new relationships among the
parameters. Microscopic model parameters, such as ion channel conductivity,
are often implicated in disease mechanisms and potential therapeutic targets,
which makes their determination valuable [42].

As with the previous example, we first trained neural networks with various
bottleneck widths (Methods) to quantify the information content of our fMRI
simulations (Figure S9 [31]). The minimum validation error occurred at $k=4$
(Figure 4A), suggesting that four latent parameters ($L_{1}$, $L_{2}$,
$L_{3}$, $L_{4}$) capture the composite effect of the original eleven
parameters.

Using a representation acquired at $k=4$, we next give an example of how to
infer latent parameters from model outputs and interpret the results in input
space. For this purpose, we first perturbed the $\text{Na}^{+}$ reversal
potential ($V_{\text{Na}}$) and considered two conditions
($V_{\text{Na}}=\\{0.48,0.54\\})$ while keeping all other parameters fixed
(Figure 4B). We observed that increasing $V_{\text{Na}}$, on average, reduced
all correlations and introduced more anti-correlations across the brain. We
then performed parameter inference in latent space in both conditions (Figure
S10 [31]), and in each case, the global optimizer converged to a single unique
solution. The two solutions differed only along the $L_{1}$ and $L_{2}$ axes,
whereas detected changes with respect to $L_{3}$ and $L_{4}$ were minimal.

To interpret the observed shifts in the fitted latent parameters, we again
applied global sensitivity analysis (Figure 4C). We found three parameters,
$g_{\text{Na}}$, $a_{ee}$, and $r_{\text{NMDA}}$ had no effect on the latent
parameters and therefore did not affect FC model outputs. Since no shifts in
this example were observed in $L_{3}$ and $L_{4}$, 5/8 of the remaining
parameters can be eliminated. This reduced the candidate parameters from
eleven down to three, with $V_{\text{Na}}$ contributing to the shift in
$L_{1}$ and $L_{2}$ (Figure 4C) (consistent with Figure 4B). As demonstrated
by this example, FixFit can be used to narrow down potential mechanisms for
the observed changes in latent parameters. This reduction makes the remaining
uncertainty tractable to resolve in specialized experiments that target
individual parameters on the single neuron scale.

Figure 4:  Identifying novel parameter relationships in a multi-scale brain
model. A: Validation error of the Larter-Breakspear model as a function of the
bottleneck dimension $(k)$. Shown for each $k$ are the individual data points,
mean, and standard error. The minimal error occurs at $k=4$, implying that the
eleven model parameters can be summarized by four latent parameters while
still uniquely mapping to output space. B: Simulated functional connectivity
matrices (left) and corresponding fitting results in latent space (right) for
two different values of $V_{\text{Na}}$. Fitting converged to global minima in
both conditions and indicated shifts in $L_{1}$ and $L_{2}$ (right, orange).
C: Global sensitivities, as computed through SCSA, of the four latent
parameters to the 11 original model parameters. Latent parameter sensitivities
to input parameters are color-coded. Three input parameters, $g_{\text{Na}}$,
$a_{ee}$, and $r_{\text{NMDA}}$, do not appear to influence model outputs.
$V_{\text{Na}}$ (marked with a red bracket), which we investigated in our
perturbation experiment, influenced latent parameters $L_{1}$ and $L_{2}$ but
not $L_{3}$ and $L_{4}$. These sensitivities are consistent with our fitting
results in panel B, where we only detected shifts in $L_{1}$ and $L_{2}$.

##  IV Discussion

While many methods exist for fitting model parameters to data, most of them
are limited to situations where there is a need for a single, optimal
solution. However, many different parameter sets often fit the data equally
well due to complex parameter dependencies and the limited information content
of experimental measurements [14]. Therefore, to fully assess what model
features data can and cannot resolve, we need tools that describe the
interaction between these solutions. In this paper, we have presented FixFit,
a data-driven approach that uniquely identifies a minimal set of latent
features from the set of interacting model parameters. Our approach correctly
identified that, as previously established, two latent parameters were
sufficient to characterize the Kepler orbit model. Furthermore, we
demonstrated in the Larter-Breakspear brain network model that we can use
changes in such composite features to narrow down parameter candidates that
shift in response to experimental perturbations.

To validate FixFit, we used two synthetic use cases; however, our framework
can be adapted to infer parameters from new experimental data. Crucially,
experimental data can include a significant amount of measurement noise which
can introduce additional redundancies among model parameters. Our approach can
be used to gauge how resolvable complexity depends on the amount of noise.
Various methods are available to supplement simulations with noise and yield
realistic model outputs [43, 44, 45, 46].

The bottleneck layer within FixFit finds a minimal and sufficient latent
representation without any additional constraints [18], but, as with other
dimensionality reduction methods, [47, 48], the representation is not unique.
This, however, provides FixFit with an additional layer of flexibility that
can be exploited. For example, variance maximization (Varimax) or sparsity
constraints may be placed during training to select, among equally accurate
alternatives, representations where latent variables are functions of the
fewest possible input parameters [47, 49].

The idea of achieving a minimal and sufficient representation of a physical
system by implementing compression on dynamical systems has been previously
explored [50]. Our approach builds upon this notion by linking the compressed
representation to model parameters, which enables us to address redundancies
in the original parameter space. However, to ensure the utility of this
representation discovered by our tool, it must be interpreted in relation to
the input parameters. This can be a challenge due to the complexity of the
encoder that maps between these two elements. Various methods are available to
interpret the components of an already-trained neural network. These include
sensitivity analysis [51], feature importance scoring [52], activation
maximization [53], and gradient-based saliency mapping [54]. Thanks to the
modular structure of our approach, it has the advantage of adaptability and
can be interpreted using any one of these methods.

For the demonstration of our framework, we used SCSA, a sensitivity analysis
variant, to quantify the impact of each model parameter on the latent
parameters [22]. However, as with the other mentioned approaches, sensitivity
analysis does not provide explicit formulas for the nonlinear parameter
relationships identified by the neural network. Besides revealing more
information about relationships between input parameters, such formulas would
allow latent parameters to be understood as composite parameters that may be
easier to interpret (i.e., eccentricity in the Kepler orbit model). This, in
turn, could allow FixFit to develop intuitive, simplified models. Symbolic
regression offers a possible avenue to achieve this by approximating the
encoder [55, 56] or by itself being incorporated into the neural network [57,
58].

Bayesian techniques are also a viable means of inferring parameters, with the
benefit of providing uncertainties for the estimates, even in the presence of
moderate degeneracies [59, 60, 61]. Nevertheless, they are prone to fail in
cases of pronounced structural non-identifiability, where due to bad mixing of
the sampling algorithm, the true values are not captured [12]. In addition,
assessing interplay among parameters from posteriors becomes challenging in
high dimensional cases. In contrast, our approach adeptly handles structural
non-identifiability as the latent representation ensures interpretability even
for high dimensional models, such as the eleven-parameter brain network model
we show here. Although distinct in their approaches, Bayesian techniques and
our approach could be combined synergistically to harness their individual
strengths. Once the low dimensional latent representation is established,
Bayesian methods could be applied efficiently to consider priors and quantify
uncertainties of the latent parameters in the presence of noise.

Future studies could utilize our framework to rank different experimental
designs based on how well they resolve specific parameters of interest in
simulations. In neuroscience, for example, brain-wide activities can be
measured in a multitude of ways, including by fMRI, electroencephalography
(EEG), and magnetoencephalography (MEG). FixFit could use simulations of these
modalities from models (e.g., Larter-Breakspear) to determine which parameters
each can resolve and, as a result, when to use them. Further, our framework
can suggest targeted follow-up experiments to augment pre-existing data and
inferred latent parameters. Related methods of optimal design aim to select
experimental conditions such as temperatures, concentrations, and sample sizes
to most easily resolve the parameters of a fixed model [62, 63, 64]. Thus, our
tool could work synergistically with these prior approaches to optimize future
data collection.

##  V Conclusion

In conclusion, we present FixFit, a method to identify unique representations
of redundant model parameters of computational models. Our tool enables
scientists to select the most informative measurement modalities for their
experiments, particularly in fields such as neuroscience, where it is still
ambiguous what information experiments can resolve. Additionally, using our
approach, researchers can utilize an extensive library of existing
computational models for parameter inferences. These efforts can be further
enhanced by specialized experiments that address the remaining parameter
uncertainties, leading to greater synergy between modeling and
experimentation.

## Code availability

All relevant code and data can be accessed at www.lcneuro.org/analytic/fixfit.

## Acknowledgements

Research presented here was funded by the following: National Science
Foundation/White House Brain Research Through Advancing Innovative
Technologies (BRAIN) Initiative, United States (NSFNCS-FR 1926781, LRMP),
Baszucki Brain Research Fund, United States (LRMP), NIHGM MSTP Training Award,
United States (T32-GM008444, AGC), and the Marie-Josée Kravis Fellowship in
Quantitative Biology (CW). In addition, we would like to thank David Hofmann
for valuable discussions.

## References

  * Group  _et al._ [2012] G. Group, M. Baak, M. Goebel, J. Haller, A. Hoecker, D. Kennedy, K. Mönig, M. Schott, and J. Stelzer, Updated status of the global electroweak fit and constraints on new physics, The European Physical Journal C 72, 2003 (2012).
  * D’haeseleer  _et al._ [1999] P. D’haeseleer, X. Wen, S. Fuhrman, and R. Somogyi, Linear modeling of mrna expression levels during cns development and injury, in  _Biocomputing’99_ (World Scientific, 1999) pp. 41–52.
  * Hodgkin and Huxley [1952] A. L. Hodgkin and A. F. Huxley, Currents carried by sodium and potassium ions through the membrane of the giant axon of loligo, The Journal of physiology 116, 449 (1952).
  * Prinz  _et al._ [2004] A. A. Prinz, D. Bucher, and E. Marder, Similar network activity from disparate circuit parameters, Nature neuroscience 7, 1345 (2004).
  * Kacprzak and Fluri [2022] T. Kacprzak and J. Fluri, Deeplss: Breaking parameter degeneracies in large-scale structure with deep-learning analysis of combined probes, Physical Review X 12, 031029 (2022).
  * Gutenkunst  _et al._ [2007] R. N. Gutenkunst, J. J. Waterfall, F. P. Casey, K. S. Brown, C. R. Myers, and J. P. Sethna, Universally sloppy parameter sensitivities in systems biology models, PLoS computational biology 3, e189 (2007).
  * Chis  _et al._ [2011] O.-T. Chis, J. R. Banga, and E. Balsa-Canto, Structural identifiability of systems biology models: a critical comparison of methods, PloS one 6, e27755 (2011).
  * Transtrum  _et al._ [2015] M. K. Transtrum, B. B. Machta, K. S. Brown, B. C. Daniels, C. R. Myers, and J. P. Sethna, Perspective: Sloppiness and emergent theories in physics, biology, and beyond, The Journal of chemical physics 143 (2015).
  * Weise [2009] T. Weise, Global optimization algorithms-theory and application, Self-Published Thomas Weise 361 (2009).
  * Ryom and Treves [2023] K. I. Ryom and A. Treves, Speed inversion in a potts glass model of cortical dynamics, P R X Life 1, 013005 (2023).
  * Moré [2006] J. J. Moré, The levenberg-marquardt algorithm: implementation and theory, in  _Numerical Analysis: Proceedings of the Biennial Conference Held at Dundee, June 28–July 1, 1977_ (Springer, 2006) pp. 105–116.
  * Wieland  _et al._ [2021] F.-G. Wieland, A. L. Hauber, M. Rosenblatt, C. Tönsing, and J. Timmer, On structural and practical identifiability, Current Opinion in Systems Biology 25, 60 (2021).
  * Pohjanpalo [1978] H. Pohjanpalo, System identifiability based on the power series expansion of the solution, Mathematical biosciences 41, 21 (1978).
  * Cole [2020] D. Cole,  _Parameter redundancy and identifiability_ (CRC Press, 2020).
  * Kramer [1991] M. A. Kramer, Nonlinear principal component analysis using autoassociative neural networks, AIChE journal 37, 233 (1991).
  * Tishby  _et al._ [2000] N. Tishby, F. C. Pereira, and W. Bialek, The information bottleneck method, arXiv preprint physics/0004057 (2000).
  * Tishby and Zaslavsky [2015] N. Tishby and N. Zaslavsky, Deep learning and the information bottleneck principle, in  _2015 ieee information theory workshop (itw)_ (IEEE, 2015) pp. 1–5.
  * Achille and Soatto [2017] A. Achille and S. Soatto, On the emergence of invariance and disentangling in deep representations, arXiv preprint arXiv:1706.01350 125, 126 (2017).
  * Hornik  _et al._ [1989] K. Hornik, M. Stinchcombe, and H. White, Multilayer feedforward networks are universal approximators, Neural networks 2, 359 (1989).
  * Hornik [1991] K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural networks 4, 251 (1991).
  * Csáji  _et al._ [2001] B. C. Csáji _et al._ , Approximation with artificial neural networks, Faculty of Sciences, Etvs Lornd University, Hungary 24, 7 (2001).
  * Li  _et al._ [2010] G. Li, H. Rabitz, P. E. Yelvington, O. O. Oluwole, F. Bacon, C. E. Kolb, and J. Schoendorf, Global sensitivity analysis for systems with independent and/or correlated inputs, The journal of physical chemistry A 114, 6022 (2010).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Scipy 1.0: fundamental algorithms for scientific computing in python, Nature methods 17, 261 (2020).
  * Bate  _et al._ [2020] R. R. Bate, D. D. Mueller, J. E. White, and W. W. Saylor,  _Fundamentals of astrodynamics_ (Courier Dover Publications, 2020).
  * Larter  _et al._ [1999] R. Larter, B. Speelman, and R. M. Worth, A coupled ordinary differential equation lattice model for the simulation of epileptic seizures, Chaos: An Interdisciplinary Journal of Nonlinear Science 9, 795 (1999).
  * Breakspear  _et al._ [2003] M. Breakspear, J. R. Terry, and K. J. Friston, Modulation of excitatory synaptic coupling facilitates synchronization and complex dynamics in a biophysical model of neuronal dynamics, Network: Computation in Neural Systems 14, 703 (2003).
  * Endo  _et al._ [2020] H. Endo, N. Hiroe, and O. Yamashita, Evaluation of resting spatio-temporal dynamics of a neural mass model using resting fmri connectivity and eeg microstates, Frontiers in computational neuroscience 13, 91 (2020).
  * Chesebro  _et al._ [2023] A. G. Chesebro, L. R. Mujica-Parodi, and C. Weistuch, Ion gradient-driven bifurcations of a multi-scale neuronal model, Chaos, Solitons & Fractals 167, 113120 (2023).
  * Desikan  _et al._ [2006] R. S. Desikan, F. Ségonne, B. Fischl, B. T. Quinn, B. C. Dickerson, D. Blacker, R. L. Buckner, A. M. Dale, R. P. Maguire, B. T. Hyman, _et al._ , An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest, Neuroimage 31, 968 (2006).
  * Friston  _et al._ [2003] K. J. Friston, L. Harrison, and W. Penny, Dynamic causal modelling, Neuroimage 19, 1273 (2003).
  * [31] See supplemental material at ¡url¿.
  * Mujica-Parodi  _et al._ [2020] L. R. Mujica-Parodi, A. Amgalan, S. F. Sultan, B. Antal, X. Sun, S. Skiena, A. Lithen, N. Adra, E.-M. Ratai, C. Weistuch, _et al._ , Diet modulates brain network stability, a biomarker for brain aging, in young adults, Proceedings of the National Academy of Sciences 117, 6170 (2020).
  * Biswal  _et al._ [1996] B. Biswal, E. A. Deyoe, and J. S. Hyde, Reduction of physiological fluctuations in fmri using digital filters, Magnetic resonance in medicine 35, 107 (1996).
  * Biswal  _et al._ [1995] B. Biswal, F. Zerrin Yetkin, V. M. Haughton, and J. S. Hyde, Functional connectivity in the motor cortex of resting human brain using echo-planar mri, Magnetic resonance in medicine 34, 537 (1995).
  * Lu  _et al._ [2017] Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang, The expressive power of neural networks: A view from the width, Advances in neural information processing systems 30 (2017).
  * Kidger and Lyons [2020] P. Kidger and T. Lyons, Universal approximation with deep narrow networks, in  _Conference on learning theory_ (PMLR, 2020) pp. 2306–2327.
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, _et al._ , Tensorflow: Large-scale machine learning on heterogeneous systems (2015).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Herman and Usher [2017] J. Herman and W. Usher, SALib: An open-source python library for sensitivity analysis, The Journal of Open Source Software 2, 10.21105/joss.00097 (2017).
  * Iwanaga  _et al._ [2022] T. Iwanaga, W. Usher, and J. Herman, Toward SALib 2.0: Advancing the accessibility and interpretability of global sensitivity analyses, Socio-Environmental Systems Modelling 4, 18155 (2022).
  * Logothetis [2008] N. K. Logothetis, What we can do and what we cannot do with fmri, Nature 453, 869 (2008).
  * de Lores Arnaiz and Ordieres [2014] G. R. de Lores Arnaiz and M. G. L. Ordieres, Brain na+, k+-atpase activity in aging and disease, International journal of biomedical science: IJBS 10, 85 (2014).
  * Scales and Snieder [1998] J. A. Scales and R. Snieder, What is noise?, Geophysics 63, 1122 (1998).
  * Fuller [2009] W. A. Fuller,  _Introduction to statistical time series_ (John Wiley & Sons, 2009).
  * Gravel  _et al._ [2004] P. Gravel, G. Beaudoin, and J. A. De Guise, A method for modeling noise in medical images, IEEE Transactions on medical imaging 23, 1221 (2004).
  * Goldwyn and Shea-Brown [2011] J. H. Goldwyn and E. Shea-Brown, The what and where of adding channel noise to the hodgkin-huxley equations, PLoS computational biology 7, e1002247 (2011).
  * Kaiser [1958] H. F. Kaiser, The varimax criterion for analytic rotation in factor analysis, Psychometrika 23, 187 (1958).
  * Xu  _et al._ [2003] W. Xu, X. Liu, and Y. Gong, Document clustering based on non-negative matrix factorization, in  _Proceedings of the 26th annual international ACM SIGIR conference on Research and development in informaion retrieval_ (2003) pp. 267–273.
  * Ghasedi Dizaji  _et al._ [2017] K. Ghasedi Dizaji, A. Herandi, C. Deng, W. Cai, and H. Huang, Deep clustering via joint convolutional autoencoder embedding and relative entropy minimization, in  _Proceedings of the IEEE international conference on computer vision_ (2017) pp. 5736–5745.
  * Chen  _et al._ [2022] B. Chen, K. Huang, S. Raghupathi, I. Chandratreya, Q. Du, and H. Lipson, Automated discovery of fundamental variables hidden in experimental data, Nature Computational Science 2, 433 (2022).
  * Zhang [2019] P. Zhang, A novel feature selection method based on global sensitivity analysis with application in machine learning-based prediction model, Applied Soft Computing 85, 105859 (2019).
  * Murdoch  _et al._ [2019] W. J. Murdoch, C. Singh, K. Kumbier, R. Abbasi-Asl, and B. Yu, Definitions, methods, and applications in interpretable machine learning, Proceedings of the National Academy of Sciences 116, 22071 (2019).
  * Montavon  _et al._ [2018] G. Montavon, W. Samek, and K.-R. Müller, Methods for interpreting and understanding deep neural networks, Digital signal processing 73, 1 (2018).
  * Selvaraju  _et al._ [2017] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, Grad-cam: Visual explanations from deep networks via gradient-based localization, in  _Proceedings of the IEEE international conference on computer vision_ (2017) pp. 618–626.
  * Schmidt and Lipson [2009] M. Schmidt and H. Lipson, Distilling free-form natural laws from experimental data, science 324, 81 (2009).
  * La Cava  _et al._ [2021] W. La Cava, P. Orzechowski, B. Burlacu, F. O. de França, M. Virgolin, Y. Jin, M. Kommenda, and J. H. Moore, Contemporary symbolic regression methods and their relative performance, arXiv preprint arXiv:2107.14351 (2021).
  * Petersen  _et al._ [2019] B. K. Petersen, M. L. Larma, T. N. Mundhenk, C. P. Santiago, S. K. Kim, and J. T. Kim, Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients, arXiv preprint arXiv:1912.04871 (2019).
  * Kim  _et al._ [2020] S. Kim, P. Y. Lu, S. Mukherjee, M. Gilbert, L. Jing, V. Čeperić, and M. Soljačić, Integration of neural network-based symbolic regression in deep learning for scientific discovery, IEEE transactions on neural networks and learning systems 32, 4166 (2020).
  * Carlin and Chib [1995] B. P. Carlin and S. Chib, Bayesian model choice via markov chain monte carlo methods, Journal of the Royal Statistical Society: Series B (Methodological) 57, 473 (1995).
  * Jaakkola and Jordan [2000] T. S. Jaakkola and M. I. Jordan, Bayesian parameter estimation via variational methods, Statistics and Computing 10, 25 (2000).
  * Gonçalves  _et al._ [2020] P. J. Gonçalves, J.-M. Lueckmann, M. Deistler, M. Nonnenmacher, K. Öcal, G. Bassetto, C. Chintaluri, W. F. Podlaski, S. A. Haddad, T. P. Vogels, _et al._ , Training deep neural density estimators to identify mechanistic models of neural dynamics, Elife 9, e56261 (2020).
  * Chaloner and Verdinelli [1995] K. Chaloner and I. Verdinelli, Bayesian experimental design: A review, Statistical science , 273 (1995).
  * Liepe  _et al._ [2013] J. Liepe, S. Filippi, M. Komorowski, and M. P. Stumpf, Maximizing the information content of experiments in systems biology, PLoS computational biology 9, e1002888 (2013).
  * Smucker  _et al._ [2018] B. Smucker, M. Krzywinski, and N. Altman, Optimal experimental design, Nat. Methods 15, 559 (2018).