Estimating a brain network predictive of stress and genotype with supervised autoencoders

Abstract Targeted brain stimulation has the potential to treat mental illnesses. We develop an approach to help design protocols by identifying relevant multi-region electrical dynamics. Our approach models these dynamics as a superposition of latent networks, where the latent variables predict a relevant outcome. We use supervised autoencoders (SAEs) to improve predictive performance in this context, describe the conditions where SAEs improve predictions, and provide modelling constraints to ensure biological relevance. We experimentally validate our approach by finding a network associated with stress that aligns with a previous stimulation protocol and characterizing a genotype associated with bipolar disorder.


Introduction
Mental illnesses are among the most debilitating medical conditions due to their prevalence (Lewinsohn et al., 1998) and the vast toll they take on individuals (Monkul et al., 2007) and society (Leon et al., 1998).Yet the brain and the malfunctioning neural dynamics responsible for mental illness remain largely a mystery.Because of this, the National Institutes of Health is actively encouraging research (Insel et al., 2013) to understand the brain at all levels of organization.Of the many types of electrophysiology commonly studied, our work focuses on understanding the dynamics of region-specific voltages in the brain, known as local field potentials (LFPs).Field potentials have proven to be good candidates for understanding dynamics because the learned features are predictive of neuropsychiatric conditions (Veerakumar et al., 2019) and can generalize to new individuals (Y.Li et al., 2017), allowing us to characterize and predict diseases using features that apply to the general population rather than a specific individual.
Recent developments have allowed neuroscientists to test hypotheses on neural dynamics by directly modifying the brain dynamics underlying actions and behaviours.For example, optogenetics (Boyden et al., 2005) is a set of techniques that use a virus to modify the cells in specific brain regions to respond to light by firing action potentials.The neural dynamics can then be modified by shining light through implanted fibre optic cables into these modified regions, inducing aggregate cell behaviour that in turn affects the LFPs (Buzsáki et al., 2012).The benefits of this development cannot be overstated.Neuroscientists can now use experimental manipulation to argue for the causality of cell types in behaviour, rather than mere correlation.But more importantly, it opens the door for targeted neurostimulation treatments of mental illnesses.Many of the current treatment methods for psychiatric disorders rely on drugs, which can have debilitating side effects (Short et al., 2018) due to their broad distribution in the brain and the low temporal resolution at which they impact brain cells.Direct manipulation of the faster dynamics responsible for the illness has the potential to cure these diseases without the negative consequences of medication.
However, we must move beyond mere predictive models to obtain the necessary insight into causal mechanisms for such interventions to succeed.For stimulation to be effective we must find a neural circuit responsible for the observed dynamics (Kravitz et al., 2010) rather than just finding predictors that correlate with the behaviour.Such collections of neural circuits, heretofore referred to as networks, inferred by a statistical method must be biologically relevant and provide enough information about brain function to find targets for stimulation.This is analogous to creating a harmony in an orchestra: one must be able to identify the out-of-tune instrument rather than merely detecting the errant notes produced.
This analogy of an orchestra is particularly apt for neuroscience where it is commonly assumed that a small number of these unobserved networks give rise to the high-dimensional dynamics (Medaglia et al., 2015).This viewpoint has often led researchers to use factor analysis as a method for understanding brain activity networks.Many factor models have been created to incorporate varying assumptions about the latent networks and resulting dynamics, including independent component analysis (ICA) (Du et al., 2016), non-negative matrix factorizations (NMF) (Lee & Seung, 1999), and principal component analysis (PCA) (Feng et al., 2015).Some factor models have been developed to specifically identify possible locations where stimulation would be effective (Gallagher et al., 2017).
Our objective is to find a brain network in a mouse model associated with two traits: being in a stressful state and having a genotype that has been linked to bipolar disorder, Clock-Δ19, compared to wild-type littermate controls using the data described by Carlson et al. (2017).Mice with the Clock-Δ19 genotype have been shown to be resilient to stressful situations (Dzirasa et al., 2011;Roybal et al., 2007;Sidor et al., 2015), so understanding how mice with and without the genotype differ can assist in understanding of stress response networks.The data in this setting are LFP recordings from both Clock-Δ19 and wild-type mice in a variety of situations ranging from not stressful to highly stressful.It is often a scientific goal to find a single network that differentiates between mice in different groups.This goal necessitates that the predictive information is located in a single factor, and this factor must be interpretable to neuroscientists in order to find candidate locations for stimulation.Manipulation of the dynamics will not be possible unless all these conditions are met.
Factor models that focus exclusively on modelling the electrophysiology data often fail to find networks that are predictive of traits of interest (i.e., stress condition, genotype).This suggests using joint factor models that assume common latent variables underlying both electrophysiology and traits.Unfortunately, joint factor models place a very strong prior on the predictive coefficients that can substantially degrade predictive performance (Hahn et al., 2013), particularly if the number of estimated factors is smaller than the number of true factors.The true number of networks responsible for neural dynamics are numerous, certainly more than the number of factors we can reliably estimate.In addition, most of these networks are unrelated to the traits of interest.These irrelevant networks can be particularly dominant, for example those related to motion (Khorasani et al., 2019) or blinking (Joyce et al., 2004).One potential solution is to increase the weight on the predictive component when fitting the joint model.However, we find that although such an approach can do well in-sample, it fails to accurately infer predictive networks for test subjects based on information on electrophysiology data alone and hence cannot address our goals.
Supervised autoencoders (SAEs) (Ranzato & Szummer, 2008) have arisen as an alternative to classical joint factor models.An autoencoder is a method for dimensionality reduction that maps inputs to a lower dimensional space using a neural network known as an encoder.The latent space is then used to reconstruct the original observations by a neural network called a decoder (Goodfellow et al., 2016).Both the encoder and the decoder are learned simultaneously to minimize a reconstructive loss.Supervised autoencoders add another transformation from the latent factors to predict the outcome of interest, thus encouraging the latent factors to represent the outcome well (Xu et al., 2017).These have been used with great success in image recognition, especially in the context of semi-supervised learning (Le et al., 2018;Pezeshki et al., 2016;Pu et al., 2016).In this work, we show how SAEs can be adapted to solve a difficult problem in neuroscience: finding a single network suitable for stimulation to modify a trait.Finding this predictive network has been a difficult problem in neuroscience, as experimental predictive ability has failed to match theoretical expectations with traditional inference techniques.We show empirically on synthetic data that model misspecification in generative models is a substantial contributor to these difficulties.We then demonstrate on both synthetic data and our experimental data set that our SAE-based approach is able to successfully identify such a predictive network, even under substantial misspecification.Interpreting this network leads to natural conclusions on potential stimulation targets.Previous studies have experimentally modified this target to successfully modulate animal behaviour (Carlson et al., 2017).Together, these contributions provide substantial evidence for the promise of our method as a useful tool for experimentalists to develop stimulation methods as potential treatments for mental illness.Notably, the proposed methodology has been used to successfully design a targeted neurostimulation protocol (Block et al., 2022;Mague et al., 2022), providing further evidence of this claim.
This paper is organized as follows: Section 2 contains a description of the data and motivation, along with defining joint factor analysis and demonstrating drawbacks under model misspecification.Section 3 derives our SAE approach while considering some basic properties and issues to motivate modifications of previous SAE frameworks.Section 4 provides two synthetic examples to illustrate the benefits of SAEs, one with a standard NMF model and another using synthetic LFPs.Section 5 shows that our approach learns predictive networks of genotype and stress, describes how the inferred 'stress network' could be modified through stimulation, and relates the network to previous literature.In Section 6, we provide concluding remarks and disucss further directions for research.Code for reproducing these results can be found at https://github.com/carlson-laband the LFP data is publicly hosted at the Duke University Research Data Repository (Carlson et al., 2023).

Data and model
In this section, we introduce the electrophysiological data analysed in this paper and the scientific motivation for a latent variable model.We highlight the need for a single latent variable to predict the experimental outcomes; here, these outcomes relate to an animal model of stress and a genotype associated with bipolar disorder.As part of this section, we highlight the importance of robustness to certain types of model misspecification (here, primarily latent dimensionality), which are nearly ubiquitous in neuroscience.Unfortunately, typical inference strategies for latent variable models are not robust to such misspecification, as we show in simulations.We show misspecification has a strong deleterious effect on previous approaches to learn predictive factor models.

Electrophysiology: the tail suspension test
The LFPs analysed in the tail suspension test (TST) came from 26 mice of two genetic backgrounds (14 wild type and 12 Clock-Δ19).The Clock-Δ19 genotype has been used to model bipolar disorder (van Enkhuizen et al., 2013).Each mouse was recorded for 20 min across 3 behavioural contexts: 5 min in its home cage (non-stressful), 50 min in an open field (arousing, mildly stressful), and 100 min suspended by its tail (highly stressful).Data were recorded from 11 biologically relevant brain regions with 32 electrodes (multiple electrodes per region) at 10,000 Hz.These redundant electrodes were implanted to allow for the removal of faulty electrodes and electrodes that missed the target brain region.We chose to average the signals per region to yield an 11-dimensional time series per mouse.This was done because the brain region represents the smallest resolvable location when modelling multiple animals; multiple electrodes function as repeated measurements and averaging allows us to reduce the variance of the measured signal.A visualization of these data and the experimental design can be seen in Figure 1.
We want to determine a single brain network that predicts stressful activity, so we consider all data from the home cage as the negative label (stress-free condition) and all data from the other two conditions as the positive label (stressed condition).A second prediction task is to determine brain differences between the genetic conditions (i.e., what underlying differences are there between the wild type and bipolar mouse model?).There is strong evidence to support the belief that the behaviourally relevant aspects of electrophysiology are frequency-based power within brain regions, as well as frequency-based coherence between brain regions (Hultman et al., 2016;Uhlhaas et al., 2008).We discretized the time series into 1-s windows to model how these spectral characteristics change over time.Windows containing saturated signals were removed (typically due to motion artefacts).While there are methods that can characterize the spectral features on a continuous time scale (Prado & West, 2010), the behaviour we deal with changes over a longer scale than the observed dynamics.Consequently, it is more effective to discretize and obtain sharper spectral feature estimates (Cohen, 1995) that are more amenable to factor modelling.
We chose to extract the relevant quantities from the recorded data prior to modelling rather than extracting spectral features in the modelling framework for simplicity; the extra modelling step would substantially increase the number of parameters in the model.The features related to power were computed from 1 to 56 Hz in 1 Hz bands using Welch's method (Welch, 1967), which is recommended in the neuroscience literature (Kass et al., 2014).We chose 56 Hz as a threshold to avoid the substantial noise induced at 60 Hz from the recording equipment, as prior literature has demonstrated that much of the meaningful information is contained in the lower frequencies.We calculated mean squared coherence at the same frequency bands for every pair of brain regions (Rabiner & Gold, 1975).This procedure converted each window of data into 3,696 non-negative power and coherence features.Figure 2 shows two LFPs and the associated features calculated from the recordings.

Difficulties with joint factor models
After extracting the relevant quantities described above, we have high-dimensional data (which scales by the square of the brain regions and number of frequencies considered).A natural approach for characterizing the data is to use latent variable models.These models posit that an unobserved latent space (factors or in the neuroscience field 'networks') generates the observed dynamics, with the mapping between the two (loadings) explaining the correlations in the observed covariates.As the dimensionality of the latent space is often far lower than the observed space, this allows for efficient representation of the correlations in the data (Bishop, 2006).The covariates in neuroscience often have strong correlations, so it is unsurprising that latent variable models have been extensively used (Bassett & Sporns, 2017).The latent representation must also contain information relevant to the outcomes relating to stress and genotype.
There is a vast literature of methods that reduce dimensionality while retaining predictive information.One widely studied approach is sufficient dimensionality reduction (SDR), which exploits the assumption that the response is conditionally independent of the predictors given a projection of the predictors into a lower dimensional subspace.Examples of this technique include sliced inverse regression (K.-C.Li, 1991) and extensions (Fukumizu et al., 2009;Wu et al., 2009).However, these approaches require stringent assumptions on the data and the outcome, which limits their utility in neuroscience.In particular, these models induce a strong prior on the predictive coefficients, which leads results to be sensitive to model misspecification.This difficulty can be seen particularly clearly with misspecification of the latent dimensionality (Hahn et al., 2013).
Joint factor models are appealing modelling choices to incorporate information about both the covariates and the outcomes, as demonstrated by the popularity and use of supervised probabilistic PCA (Yu et al., 2006), adding supervised variables in highly related topic models (McAuliffe & Blei, 2007), and the use of the concept in deriving supervised dictionary learning algorithms (Mairal et al., 2009).We now formally introduce these models.To define our notation, let {x i } i=1,...,N ∈ X be the measured dynamics (predictors) and {y i } i=1,...,N ∈ Y be the corresponding trait or outcome.The latent factors are denoted {s i } i=1,...,N ∈ R L , Θ are the parameters relating the factors to x, and Φ are the parameters relating the factors to y.To handle intractable integrals involved in marginalizing over the latent factor distributions (Gallagher et al., 2017), two common approaches are to optimize a variational lower bound similar to a variational autoencoder (Kingma & Welling, 2013) or alternatively to estimate the factors jointly with the model parameters (Mairal et al., 2009).We choose to proceed with the latter, due to its previous usage in similar applications.
Given this strategy, a generic objective function for a joint factor model is max (1) The tuning parameter μ controls the relative weight between reconstruction and supervision.
Including μ is somewhat atypical, and we note that setting μ = 1 recovers the log-likelihood of a standard joint factor model.In practice, it may be important to modify μ to upweight the importance of prediction by setting μ ≪ 1.This term can be thought to correspond to a fractional posterior p(x | Θ) μ p(y | x), which has been previously used in robust statistics (Bhattacharya et al., 2019), and is highly related to the common practice of modifying the noise variances on x and y in joint factor models (Yu et al., 2006).This tuning parameter is frequently used in machine learning (Mairal et al., 2009).Given this objective function, we find the optimal values of the latent factors, Θ, and predictive coefficients Φ.This is commonly done with stochastic methods due to their computational advantages in large data sets.However, joint models are not without their own drawbacks, as we now demonstrate using analytic solutions with an L 2 loss on a misspecified model.Furthermore, as motivation for our approach, we show that SAEs are not as affected by this type of misspecification.We generated a synthetic data set corresponding to supervised probabilistic PCA (Yu et al., 2006), a simple joint linear model.We set the number of predictors to 20 with a single outcome and a latent dimensionality of 3. The largest factor was unpredictive of the outcome.Additionally, we limited the inferred dimensionality to 2. This mimics our statistical models that may under-specify the true number of brain networks.We fit a linear joint model and an SAE, which we will fully describe in Section 3, with two components over a range of supervision strengths using the analytic solutions and show the results in Figure 3.At low supervision strengths (large values of μ), both the joint model and SAE focus on reconstructing the data.At higher strengths (small values of μ), they sacrifice some of the reconstruction for better predictions of the outcome.However, the joint model simply overfits the factors by making them overly dependent on y, leading to an effect we refer to as 'factor dragging'.We define this as A visualization of the features used for cross-spectral factor analysis (CSFA)-NMF.On the left, we show a 1-s window for the LFPs in two brain regions.The middle shows the power spectra for each LFP from 1 to 56 Hz estimated using Welch's method.The right shows the mean squared coherence, which we used to quantify the coherence between the two LFPs.We show two channels for clarity; in practice we have many channels and calculate all pairwise coherences.
where g(•) and f (•) represent the mapping to s when both x and y are observed and when just x is observed, respectively.This represents the difference between the factor estimates when the outcome is known (training) and when the outcome is unobserved (testing).For a joint factor model in this set-up, this corresponds to g(x i , When the latent space is dedicated to reconstructing x, this difference will be small.However, as μ is increasingly small, the values of s will be largely influenced by the predictive loss on y when it is observed, and this difference will become increasingly substantial when y is not observed.This is particularly relevant to our applications, as the dimensionality of x often corresponds to between 3,000 and 10,000 features.In order to influence the latent space, we often must roughly 'balance' the reconstructive loss and predictive loss, corresponding to a value of μ ≈ 1/3,000.The rightmost plot of Figure 3 shows this pattern, which is increasingly large with strong supervision.Large discrepancies will lead to poor predictions, as it indicates that the absence of knowledge of the outcome y (inherent to prediction) dramatically affects the latent representation, and by extension the reconstruction/prediction. In contrast, our SAE approach, which uses the same map from x regardless of whether y is observed, circumvents these issues and is broadly applicable as we fully explain in Section 3. Some efforts have been made to overcome these known issues, particularly task-driven dictionary learning (Mairal et al., 2011).While task-driven dictionary learning addresses this specific issue, it is developed specifically as a matrix decomposition method and lacks the flexibility as to be incorporated in many neuroscience models, such as Gaussian process-based models (Gallagher et al., 2017).While these models are not explored in this work, our approach naturally generalizes to these models.
There have been efforts to fix this problem, such as 'cutting-the-feedback' approaches (Liu et al., 2009;McCandless et al., 2010;Plummer, 2014).In Markov chain Monte-based methods with this approach, the samples of the latent variables are drawn using only the electrophysiology (hence the influence from the outcome is cut).The corresponding maximum likelihood technique fits a factor model exclusively on the electrophysiology and then fits a predictive model using the estimated factors.The problem of cutting-the-feedback is immediately apparent; dominant networks are often associated with motion, maintaining homeostasis, or other irrelevant behaviours.The latent space will likely not contain information relevant to prediction when the outcome is ignored (Jolliffe, 1982).
An alternative solution for factor models is to substantially increase the number of latent components.This will capture the relevant information, but the information will be scattered across multiple networks.This is particularly harmful in our application, as the latent network itself is used to choose stimulation targets (Carlson et al., 2017).If the predictive information is located in a single network, neuroscientists can choose regions that are active and well connected in the predictive network to stimulate or run downstream tests on a single variable with higher statistical power.If the generative model is a good representation of the neural dynamics, this will result in a global stimulation of the predictive network resulting in altered behaviour.However, if the predictive information is scattered across multiple networks, it is impossible to identify which connected regions of the networks will result in altered behaviour.We will show that all these issues arise when predicting stress and genotype using electrophysiology on the TST data set, and that our SAE-based approach addresses these issues, resulting in substantially improved performance.

An SAE-based approach
We now introduce a novel SAE-based approach to learn a single predictive latent factor, and discuss reasons for superior performance in practice.We then derive analytic solutions for SAEs in a linear model.These solutions are used to contrast the behaviour of SAEs to traditional inference techniques in an illustrative example.Finally, we analyse appropriateness of SAEs in biological applications based on common assumptions.

Why do SAEs yield superior predictions compared to joint factor models?
Given a new sample x * , prediction of y * given x * is straightforward once Θ and Φ have been estimated; simply compute point estimates of the factors and outcome as ŝ * = arg max log p(s | x * ) and ŷ * = arg max log p(y | ŝ * ).However, these predictions are often quite poor.This issue does not stem from the inference method, as the example from Section 2 maximized Θ and Φ directly and yet still yielded poor predictions.Instead, the issue is due to the model attempting to explain all the variance of the outcome using the latent factors (Hahn et al., 2013).If there is a high amount of weight on predictive performance, a joint model will discover that the best predictor of y is y itself.This influence is most keenly felt when L is small, as it is difficult to force the model to prioritize prediction of y without increasing the dependence on y for proper estimation of s.
While there are a plethora of methods for selecting L, such as cross-validation or carefully designed priors (Bhattacharya & Dunson, 2011), these methods will be largely unhelpful in improving predictive performance.These methods value a parsimonious representation of the joint likelihood, corresponding to fewer latent factors.This conflicts with a prior focussed on predictive ability which would encourage more factors.While cross-validation of the predictive accuracy over the dimensionality seems like an appealing option, in practice it still yields subpar predictions as compared to a purely predictive model and has difficulty finding a single predictive latent factor.
SAEs are an alternative to joint models that modify a deep feedforward network to include an autoencoder in order to reconstruct the predictors.Let us define the predictive feedforward network as ŷ = f Φ (A(x)), where A : X → R L is the encoder, and the reconstructive autoencoder is a composition of the encoder and decoder x = g Θ (A(x)).For estimation, it is common to place losses on the reconstructive and predictive parameters as well as the latent factors as regularization.This yields an objective function, min (3) To maintain the desired interpretability of joint models, we choose the parameters and losses to correspond to the parameters and negative log likelihoods from the joint model.An SAE learns a mapping from x to s that is used to predict y due to the supervision loss, incorporating the relevant information during training.However, this mapping is only dependent on x in contrast to joint models which depend on both variables as inputs.Once the SAE is estimated, at test time x alone is used in prediction and reconstruction of the electrophysiology.Because of this, SAEs are limited to modelling the variance in y that can be predicted by x.From this point of view, the reconstruction loss can be considered a deep-learning version of regularization discussed in Hahn et al. (2013), and was explicitly motivated as such in Le et al. (2018).In these situations, the reconstruction loss is placed on the latent space functions as a method to reduce the complexity of the model, corresponding to a structural risk minimization (SRM) approach (Vapnik, 1992).This technique has been used successfully to regularize complex models (Guyon et al., 1992;Kohler et al., 2002).However, there is one substantial difference between our approach and a typical SRM approach.With an SRM, the complexity penalty shrinks with an increasing number of observations (corresponding to μ → 0).In SAEs, μ is a fixed constant, as it is critical that the factors maintain biological relevance even with large numbers of observations.
3.2 When are SAEs an appropriate predictive model?
The previous section makes it clear that, relative to joint models, SAEs improve predictive performance and reduce the consequences of misspecification of the true latent dimensionality.However, it is difficult to intuitively see what these differences are, as most formulations of Equations ( 1) and (3) do not have analytic solutions and are learned with stochastic methods (Bottou et al., 2018).We develop novel analytic solutions when the likelihood is replaced with an L 2 loss.This yields a form similar to PCA, which can be considered the limiting case of probabilistic PCA as the variance of the conditional distribution vanishes (Bishop, 2006).We assume that both the predictors and outcome are demeaned.For convenience, we will define matrix forms of the data and factors: The solution for the concatenation of W and D of the joint factor model, [W T , D T ] T , can be found as an eigendecomposition of the matrix When Y is known, the factor estimates can be computed using the fixed-point equation detailed in the Appendix.When the outcome is unknown, the factor estimates are simply the projection of the data onto the latent space, S = (W T W) −1 W T X.This solution form is well known from PCA of the predictors and outcome jointly, except for the minor extension with μ.
The solution for the concatenation of W and D under the SAE formulation as defined above is also an eigendecomposition with a slightly different matrix where P X = X T (X T X) −1 X is a projection matrix onto X.The estimate for A can be computed via a fixed-point equation.We provide the details of the fixed-point equations and the derivations in the Appendix.The important aspect to note is that the SAE only models the variance of y that it is linearly predictable by x due to the term P X .Thus, given that the latent space is only computed with x, this formulation forces the model to find predictive factors using only the predictors.We can use these analytic formulae to examine the distributional assumptions under which SAEs provide good predictive models.Previous work by Le et al. (2018) showed that a reconstruction loss added to a predictive model provided a bound on the generalization error in linear SAEs via sensitivity analysis (Bousquet & Elisseeff, 2002).This corresponded to predictive gains with more common deep networks.However, bounds provided by sensitivity analysis can be loose (Zhang et al., 2017) and ignore properties of the data generating process (Bu et al., 2020).This led to a claim that reconstruction penalties empirically never harm performance.While this is often true in practice, an intuitive counterexample to this claim is where the latent structure is not highly related to the outcome.The reconstruction loss will force the model to focus on high-variance predictors unrelated to the outcome, harming predictive performance.On the other hand, having the predictive information correlate with high variance latent factors yields substantial improvements, as estimates of high variance components can converge with small numbers of samples (Shen et al., 2016).Therefore, it is important to evaluate whether a data set is likely to benefit from an SAE-based approach before analysis.
We empirically explore the distributional assumptions of SAEs on random matrices using the analytic solutions previously derived.Specifically, we show that SAEs are beneficial when the predictive information lies on a low dimensional manifold, especially one that correlates strongly with high-variance components.To demonstrate the first claim, we generated 200 samples of predictors x ∈ R 300 with zero mean and unit variance with a single latent component s ∈ R explaining 70% of the variance.The outcome y ∈ R was generated as y i = (1 − λ)x i β + λαs i + ϵ i .By allowing λ to vary, we can control how important the low-dimensional manifold is in predicting the outcome, as shown on the left of Figure 4. SAEs are largely ineffective when the low-dimensional manifold is largely irrelevant for predictions (small values of λ).However, the performance dramatically improves as the manifold becomes increasingly influential (large values of λ).Models that place high value on the reconstruction loss are more effective at taking advantage of this low-dimensional structure (high values of μ).
To demonstrate the second claim, we generated data x ∈ R 300 with zero mean and unit variance such that x = (1 − λ)s i1 w 1 + λs i2 w 2 + η i .The outcome y ∈ R was generated as y i = s i2 + ϵ i , using exclusively the second latent factor.By varying λ, we can explore the importance of having predictive information align with high variance components in the predictors, as is shown on the right of Figure 4.The models that place high value on the reconstruction (larger values of μ) largely ignore the predictive information when it is correlated with little variance in the predictors, while the model that emphasizes prediction overfits (smaller values of μ).However, with slight increases in predictive variance, the model with small values of μ is quicker to incorporate information relevant to the outcome.Finally, when most variance in the predictors aligns with the predictive component, the models that place higher emphasis on reconstruction exploit the regularization most effectively and have the lowest generalization error.
The sparsity-in-dimensionality assumption of SAEs aligns better with network-based neuroscience than an alternative sparsity assumption in terms of the predictors (Tibshirani, 1996).Many of our predictors are highly correlated, such as power in adjacent frequency bins in the same brain region.It is likely that both predictors will be similarly predictive of the outcome.LASSO tends to shrink the coefficient for the slightly less useful predictor to 0. An SAE would ideally find the network responsible for the variance observed in both predictors and relate this latent network to the outcome.Given the strong evidence in favour of network-based neuroscience and previous experiments demonstrating the relevance of our analysed features, it is reasonable to expect that SAEs will yield improved predictive performance relative to purely predictive or joint models when the relevant networks are moderately-to-highly expressed, but are potentially detrimental when the network has comparatively low expression.

How does predictive sparsity impact SAEs?
The SAE approach as previously defined only ensures a predictive latent subspace.However, our applications require that a single predictive network be learned, while the remaining networks model less relevant dynamics.In related settings, it is common to place a hard sparsity constraint on the parameters Φ, limiting the predictive influence to a single (or small number) of factors (Gallagher et al., 2017).This can lead to undesirable consequences when minimizing the loss in Equation (3) when the decoder involves a difficult-to-optimize objective (Ulrich et al., 2015).In particular, one locally optimal solution is to use the predictive factors to highly overfit the outcome, while the associated loadings vanish.This corresponds to a latent space factorized into a predictive and generative model; the predictive factor has no biological significance as it explains none of the dynamics, while the unpredictive factors model the electrophysiology with one-fewer latent dimensions.This local optimum is not apparent with the linear models due to more stable training properties.However, with more complex models, such as the Gaussian process model of Gallagher et al. (2017), this is a common problem.We address this issue using an SAE approach by replacing the penalization terms on the loadings and the latent factors with a normalization constraint on the loadings for each factor.Let Θ l be the loadings associated with factor l.Then, our objective function corresponding to Equation ( 3) is replaced with an objective function min where the non-negativity constraints are applied elementwise.For convenience we chose c = �� p √ for each factor.However, the choice is arbitrary; latent variable models inherently lack multiplicative identifiability.If the latent variables are multiplied by a constant k and the loadings by 1/k, the new reconstruction kA(x)Θ/k = A(x)Θ will remain unchanged.This is not an issue for selecting a stimulation target, as we choose targets based on relative importance.Incorporating these constraints is straightforward with modern packages.Specifically, we define hidden unconstrained variables Ξ, and then set the loadings as a transformation of these hidden variables Θ = h(Ξ) such that h(Ξ) fulfils the constraint.In our case, we used the softplus function, softplus(x) = log (1 + exp (x)), which smoothly maps from the full real line to the positive reals.This unconstrained optimization with a transformation is straightforward to implement, particularly with automatic differentiation.Specifically, by placing the softplus rectification as the final step, the remaining parameters are able to be learned without constraints.The objective in Equation ( 6) can be learned with respect to Φ, Ξ, and A, rather than optimizing the loadings directly.Other options, such as learning an unconstrained A and projecting into the feasible (nonnegative) region, would require specialized implementations.
In SAEs with deep decoders, implementing this type of identifiability constraint is not straightforward.However, our application requires that the SAE correspond to a joint model.Thus, the decoder corresponds to a shallow network, and there are natural normalization constraints on the associated loadings.In this model, a natural constraint on the NMF loadings is to require that the norm of each factor be equal to a constant.We enforce this constraint by first learning unconstrained variables and then dividing by the norm to obtain the loadings after the model is learned.
This normalization constraint also possesses a scientific justification.As a single factor is used to make predictions, we can choose targetable features based on the importance of the supervised factor loadings relative to the other factor loadings (e.g., the features uniquely important to the supervised network).This relevance is evaluated by dividing each entry in the supervised factor loadings by the sum of that entry in all factor loadings.The features with the highest ratios in the supervised network are potential candidates for stimulation.Additionally, as the network associated with our target behaviour is often relatively small in variance (hence the initial need for SAEs), normalizing the loadings for each factor to a constant provides a convenient method to ensure that the supervised network is not deemed irrelevant due to smaller loadings, as would be the case with penalization-based identifiability constraints.

Validation of our approach with synthetic data
We now provide two synthetic examples that demonstrate that our SAE-based approach can recover a single predictive factor that is robust to misspecification.In the first example, we generate data from a known NMF model and show that our SAE-based approach can recover a predictive factor under misspecification in the latent dimensionality.In the second demonstration, we generate synthetic LFP dynamics using an alternative latent variable model (Gallagher et al., 2017) and extract a predictive factor using our NMF model.The first example validates two aspects: first, that the estimated factor is predictive and second, that the estimated loadings for the SAE match the true loadings.The second example validates that our NMF approach can robustly extract predictive factors in a realistic simulation of experimental conditions.In both examples, we show that alternative approaches (cutting-the-feedback and supervised dictionary learning) fail under our experimental conditions.

Recovering a synthetic NMF component
We generated synthetic data such that the latent factors and the conditional distribution of the observations given the factors had truncated normal distributions.This aligns with the likelihood assumed with an L 2 loss and non-negativity constraints.We chose a latent dimensionality of 10 and a 100-dimensional observation space to match the assumption that the number of networks is substantially lower than the observed dimensionality.The latent factors were independent with distinct variances.The outcome came from a Bernoulli distribution where the probability depended exclusively on the lowest-variance component.The data generation process can be written as We fit an NMF model with five components to match the assumption that the number of estimated factors is smaller than the true number of networks.We evaluated three methods for inferring the model, cutting-the-feedback, joint modelling, and our SAE-based approach.For each model we maintained identical identifiability constraints on the loadings corresponding to Equation ( 6).We measured the predictive performance of the learned factor using the area under the curve (AUC), which can be interpreted as the probability of the predictive network being higher in a positive sample as compared to a negative sample.This was evaluated on a test set.We also quantified the accuracy of the learned loadings as the cosine similarity between the true and estimated loadings.The cosine similarity between two vectors v 1 and v 2 is defined as This quantity measures how well the two vectors align with 1 indicating perfect alignment and 0 indicating orthogonality.The cosine similarity allows us to answer how similar different brain networks are.Highly similar networks will involve similar regions with similar covariates.As such, they will largely align with cosine similarities close to 1. On the other hand, networks that involve different covariates will have lower cosine similarities.The permutative non-identifiability is accounted for in the SAE and joint models by the fact that a single component is supervised, so similarity is measured between the supervised loadings and the true loadings.In the joint NMF, this similarity is computed between the loadings associated with the most predictive factor and the true loadings.Figure 5 and Table 1 show the results of the predictive AUC of the single estimated predictive factor, along with the cosine similarity between the estimated loadings and the true network characteristics.For comparison we have included the values for the true model and logistic regression.The latter provides an estimate of the predictive ability of a standard linear model.Cosine similarity ranges from 0 to 1, and for the specific generated data set the cosine similarity between two randomly learned loadings is 0.64.Values larger than this indicate that the model successfully incorporated information from the true model into the estimated network.
The factor estimated via cutting-the-feedback contains no predictive information and does not align with the true network.This is unsurprising, as none of the estimated factors using reconstructive loss alone will contain information related to a low-variance factor as defined.More interestingly, the joint factor model provides no benefits over a cut model.This seems surprising as incorporating information related to the outcome during estimation should improve predictive performance.However, the figure on the right shows why this is not the case.The AUC on the training set of the joint model is far higher than the SAE, and higher than the theoretical bound provided by the true model.However, predictive ability returns to random chance when the factors are estimated in the test set without knowledge of the outcome.This is the 'factor dragging' problem stated in Section 2. The model overfits the predictive factors without modifying the corresponding loadings, making prediction dependent on knowledge of the outcome.Our SAE approach avoids this overfitting and accurately characterizes the network.

Extracting predictive features in synthetic LFPs
We now validate our SAE-NMF approach using synthetic LFPs to match our experimental conditions as closely as possible.The data were generated using a previously developed factor model specifically designed for analysing LFPs, referred to as Cross-Spectral Factor Analysis (Gallagher et al., 2017).These models use Gaussian processes in a multiple kernel learning framework to represent the spectral features in a low-dimensional space of latent factors.This approach functions as a generative model given a prior on the latent factors.
We initialized a CSFA model with 30 latent components and 8 measured brain regions and generated 20,000 samples of synthetic LFPs using the sampling method described in Section A.6.The draws associated with a single latent factor were used to generate a synthetic binary outcome.The associated power spectrum of a subset of the regions is plotted in Figure 6.This data set reflects many of our assumptions of brain dynamics; the latent variables have a substantial amount of sparsity, a small number of latent variables are responsible for the outcome of interest, and the spectral features associated with the latent variable are sufficient to characterize the brain network.
We then calculated the power and coherence features of the synthetic data using the same procedure described in Section 2. This allows us to examine the performance of our approach under realistic circumstances.We compare our supervised approach to a logistic regression model using the observed features, a sequentially fit NMF model, and a joint model.Each NMF model was estimated with only 5 latent factors instead of the 30 factors used to generate the data.Therefore, our example demonstrates the effect of misspecification, not only in the latent dimensionality but also in the observational likelihood of the data under realistic assumptions.
The results are shown in Table 2.The cut model fails to capture any relevant information to the outcome, which we attribute to the substantial misspecification.Correspondingly, the joint model  : Applied Statistics, 2023, Vol. 72, No. 4 also fails to yield accurate predictions as it predicts nearly perfectly on the training data but does not incorporate that information into the associated loadings.However, our SAE approach is able to discover relevant features for prediction.In fact, it slightly improves on the pure logistic regression model aligning with the theoretical results of Le et al. (2018).
The results of these two synthetic examples suggest that our approach will be able to accurately recover the relevant dynamics associated with the outcome of interest and be robust to the types of misspecification that can be reasonably expected when analysing LFPs.

Estimating networks associated with stress and bipolar disorder
We have demonstrated why and how our SAE approach can yield improved predictive performance in synthetic LFP data.Furthermore, our approach can accurately estimate the associated network characteristics in a known model and handle misspecification of both the latent distribution and observational likelihoods.We now demonstrate the predictive abilities of SAEs on the TST data set.We use an NMF model with a single supervised factor to force all the predictive information into a single brain network.Our objective is to find a single network predictive of stress and a single network predictive of genotype.We fit a model with 10 factors and a set supervision strength of μ = 1.The complete details of the implementation and the features are provided in the Appendix.
We compared the results from our approach to the two most relevant competitors, a joint factor model with a single supervised factor and 'cutting-the-feedback' with different numbers of latent  factors.Since our objective is to find a single predictive network, we tested two measures of predictive ability when comparing to cutting-the-feedback.The first was the predictive ability using the entire latent subspace and the second was the predictive ability using the single most predictive factor on the training set.One method for improving the predictive performance with cutting-the-feedback approaches is to increase the number of latent factors, so we compare the performance of a 25-factor model in addition to the corresponding 10-factor model.We also compared the reconstructive losses to ensure biological relevance of the learned model.Table 3 shows the predictive and reconstructive metrics using all stated models for predicting stress, along with 95% confidence intervals.It is apparent that our SAE-based approach is substantially better than cutting-the-feedback or fitting a joint factor model in predictive ability.Our approach with a single factor achieves an AUC of 0.94 relative to 0.76 for both the cut model and joint model.Even with more factors, a cut model approach still yields worse predictions.In fact, increasing the latent dimensionality can make the prediction using a single factor worse as the relevant information becomes divided between an increasing number of factors.The reconstruction losses were better for the cut and joint models at 0.043 and 0.045, respectively, which is unsurprising as the encoder constrains how well the SAE can adapt the factors to each observation.However, the SAE still explains a substantial portion of the variance with a reconstruction loss of 0.051 compared to 0.09 loss from using the mean, indicating that the estimated factors are still biologically relevant.
Choosing the optimal number of factors is a challenging problem.We note that several approaches exist to select the dimensionality in factor models, but they do not naturally apply in our context.For example, there are a number of Bayesian hierarchical methods that define priors for helping choose the dimensionality (Bhattacharya & Dunson, 2011) or penalization-based schemes to balance generative capability with complexity (Minka, 2000).It is relatively common in machine learning to choose the number of components based on heuristics of explained variance (Ferré, 1995) or by choosing the number of components through a cross-validation procedure.However, it is non-trivial to extend these existing techniques to the needs of our application, as it is challenging to include the hierarchical priors in the SVAE formulation, and a focus exclusively on prediction does not lend itself to finding a system that reconstructs the data well and can be explained as networks.Instead, for this application we chose 10 factors to match previously learned models (Gallagher et al., 2017).This allowed for the generative factors to account for the baseline variance while not rendering the predictive factor irrelevant in the reconstruction of the electrophysiology.
The sign of the predictive model associated with the supervised factor in the SAE was consistently negative across all the data splits.This indicates that the learned network was less active in stressful conditions as compared to a non-stressful environment.Figure 7 shows one method for visualizing the network.Each segment along the edge represents power in the labelled brain region.The 'spokes' in the circle represent coherence between the two specified regions.Stimulation targets are chosen as region/frequency combinations that are particularly influential in the supervised network.These influential combinations will have large loadings in spectral power and coherences between that region and other brain regions.By stimulating this 'central' region, we can influence the entire network and its associated behaviour (Mague et al., 2022).
In this particular case, the network negatively associated with stress is characterized by power in the low frequencies of 1-4 Hz in the prefrontal cortex (IL_Cx and PrL_Cx), thalamus (Md_thal), amygdala (BLA), and ventral tegmental area (VTA) along with high levels of coherence between these regions.It also is characterized by coherence at these frequencies within the hippocampus (mDHipp and IDHipp).Since this network is lower in the stressful situation, stimulation in the prefrontal cortex or VTA at low frequencies should modify the dynamics in a stressful situation to more closely align with the dynamics observed in non-stressful situations.Such stimulation would answer whether the network is causal, and potentially ameliorate the impact of stressful behaviour.We note that a prior stimulation procedure was used that stimulated the thalamus (Md_thal) in a phase-locked manner to 3-7 Hz waves in the prefrontal cortex (Carlson et al., 2017).This stimulation procedure aligns with the network discovered by our approach, and reduced time immobile in the TST, which is consistent with a reduction in stress.Differentiating Clock-Δ19 vs wild type dynamics is a far more difficult task as compared to characterizing stress.This would unsurprisingly indicate that tail suspension and non-stressful situations are more distinctive than the baseline dynamics due to modification of a single gene.Yet also here our SAE approach outperforms prediction using cut or joint models as shown in Table 4.In this case, adding extra factors to the cut model yields a closer (but still inferior) predictive ability to the SAE model.However, the previous task showed that a strategy relying on increasing latent dimensionality to improve predictions with a single factor will give inconsistent results.Our SAE approach provides a more reliable method for finding a single network correlated with the outcomes of interest in neuroscience.
The network found in each of the splits was consistently negatively correlated with the Clock-Δ19 genotype and the network found in one of the splits is shown in Figure 8.That this network is negatively associated with Clock-Δ19 indicates that this stimulation in Clock-Δ19 mice would make their dynamics more closely align with the wild type population.This network is Figure 7. Visualization of the brain network found to be predictive of experimental condition.We first normalize the loadings of the specific factor relative to the loadings of all factors.These values evaluate how much influence the network had on a specific feature.We then thresholded for the largest coefficients; any covariate that exceeded this threshold was included.Each segment of the circle represents a different brain region, and the numbers (10, . . ., 50) represent frequencies in Hertz.Any power features that exceeded the threshold were included as a band on the rim.Any coherence features were plotted as a 'spoke'.As an example interpretation, this network substantially influenced the power in the prefrontal cortex (PrL_Cx and IL_Cx) at low frequencies, as well as the coherence between the prefrontal cortex and the acumbens, amygdala, thalmus, and ventral tegmental area (Acb_Sh, BLA, Md_Thal, and VTA).

926
Talbot et al.    characterized by power in the prefrontal cortex at low frequencies similar to the stress network.However, this network is weaker in the hippocampus (IDHip and mDHip).This would indicate that a reasonable location for stimulation is the prefrontal cortex at low frequencies.Our emphasis in this work has been finding networks that can be used for stimulation, and therefore we have largely refrained from interpreting the scientific conclusions derived from the estimated networks in the TST data set.Nevertheless, it is easy to see the similarities between the estimated stress network and the estimated Clock-Δ19 network.In fact, the cosine similarity between the networks in Figures 5 and 6 is 0.99, whereas if we pick a random network in the task model and a random network in the genotype model, the average cosine similarity is 0.38.Based on this, our models would suggest that the differences in genotype are similar to the observed differences induced by stressful and non-stressful conditions.It is encouraging to note that this observation has matched previous scientific experimentation (Murata et al., 2005).Additionally, this implies that the stress felt by the Clock-Δ19 animals is lower than wild type (WT) animals, which aligns with the finding that they are more active during the TST (Carlson et al., 2017).We note that these findings show an appealing aspect of our approach: by using an interpretable and biologically relevant generative model in our SAE, we have been able to scientifically analyse the differences between the experimental groups.
As a final comparison, we can also empirically verify our claim that SAEs predictive ability is robust to the latent dimensionality.Figure 9 shows the predictive ability of our SAE-based approach, a SDL approach, and a 'cutting-the-feedback' as a function of latent dimensionality, predicting both genotype and stress condition using the most predictive factor.These results were obtained evaluating the general population rather than a mouse-by-mouse basis to match standard statistical evaluation techniques.We can see on both tasks that the SAE maintains a consistently high predictive ability with all latent dimensionalities.However, the other two methods not only fail to match this predictive ability but also fluctuate rather than providing a clear monotonic trend.This illustrates a substantial problem that has traditionally plagued network-based neuroscience; even simple parameter choices can cause dramatic shifts in predictive results.

Conclusions
We have developed an approach for finding a single stimulatable brain network that correlates with a condition or behaviour.Our strategy is based on supervised autoencoders, which are more conducive to finding a single network predictive of an outcome of interest than joint factor models.We are able to find a brain network that is correlated with stress and a network that is correlated with a genotype linked to bipolar disorder.This approach naturally leads to candidate brain regions for stimulation to modify the network.
While the objective of this work is developing an approach for identifying relevant brain networks, we also provide broad conceptual contributions to understand why, how, and when SAEs are effective model choices.This is done by deriving SAEs in a manner that elucidates the source of superior predictive ability relative to joint factor models, and using analytic solutions in a simple case to illustrate why carefully designed SAEs are useful in our applications.

Figure 1 .
Figure 1.On the left, we show example LFP recordings in a subset of the regions measured.The LFPs were divided into 1-s windows to yield discrete observations.On the right, we provide a simple diagram of the TST experiment.The mice were recorded in the home cage as a non-stressful environment.They were then recorded for 5 min in a mildly stressful open field test (OFT) and then 10 min in a highly stressful TST.

Figure 3 .
Figure 3.There are three true latent factors but the model has two latent factors.We vary the supervision strengths (μ) over a range of values and report the MSE for three cases: a joint model where the factors are estimated with the trait known (Joint: Y known), a joint model where the trait is unknown (Joint: Y unknown), and an SAE.
, and S = [s 1 , . . ., s N ] ∈ R L×N .We let the model parameters Θ = W and Φ = D for matrices W ∈ R p×L , D ∈ R q×L .For the associated SAE model, we assume a linear encoder, A(x) = Ax.

Figure 4 .
Figure 4.The plot on the left shows effect of sparsity in dimensionality on the true risk (generalization error).The plot on the right shows how having predictive information align with high variance components of the covariates affects generalization error.

Figure 5 .
Figure 5.The ROC curves on the training and testing sets for the estimated joint model and SAE.

Figure 6 .
Figure 6.The relative importance of the predictive CSFA factor in three of the eight generated regions along with phase offset (indicates directionality).
Bold values mark highest performance.

Figure 8 .
Figure 8.The predictive network of the Clock-Δ19 genotype, with an identical interpretation to Figure 7.

Figure 9 .
Figure9.The impact of latent dimensionality on the SAE predictions of task.SAE refers to our supervised autoencoder approach, while CTF refers to the 'cutting-the-feedback' approach.

Table 1 .
Prediction and reconstruction

Table 2 .
Prediction using synthetic local field potentials

Table 3 .
Prediction of stress.