Neurochemistry-enriched dynamic causal models of magnetoencephalography, using magnetic resonance spectroscopy

We present a hierarchical and empirical Bayesian framework for testing hypotheses about synaptic neurotransmission, based on the integration of ultra-high field magnetic resonance spectroscopy (7T-MRS) and magnetoencephalography data (MEG). A first level dynamic causal modelling of cortical microcircuits is used to infer the connectivity parameters of a generative model of individuals’ neurophysiological observations. At the second level, individuals’ 7T-MRS estimates of regional neurotransmitter concentration supply empirical priors on synaptic connectivity. We compare the group-wise evidence for alternative empirical priors, defined by monotonic functions of spectroscopic estimates, on subsets of synaptic connections. For efficiency and reproducibility, we used Bayesian model reduction (BMR), parametric empirical Bayes and variational Bayesian inversion. In particular, we used Bayesian model reduction to compare models of how spectroscopic neurotransmitter measures inform estimates of synaptic connectivity. This identifies the subset of synaptic connections that are influenced by neurotransmitter levels, as measured by 7T-MRS. We demonstrate the method using resting-state MEG (i.e., task-free recording) and 7T-MRS data from healthy adults. We perform cross-validation using split-sampling of the MEG dataset. Our results confirm the hypotheses that GABA concentration influences local recurrent inhibitory intrinsic connectivity in deep and superficial cortical layers, while glutamate influences the excitatory connections between superficial and deep layers and connections from superficial to inhibitory interneurons. The method is suitable for applications with magnetoencephalography or electroencephalography, and is well-suited to reveal the mechanisms of neurological and psychiatric disorders, including responses to psychopharmacological interventions.


Introduction
This paper introduces an empirical Bayesian methodology for inferring synaptic physiology from human, in vivo neurophysiological recordings and in vivo measurements of neurochemicals, including neurotransmitters. It is an example of a broader class of enriched dynamic causal models (DCMs), informed by multimodal data that can incorporate neurochemistry, molecular pathology or selective loss of neuronal sub-populations and their synapses. Magnetic resonance spectroscopy (MRS) can be used to estimate individual differences in key neurochemical concentrations in local brain regions (Blüml and Panigrahy, 2012). It has been used to quantify neurochemical changes in neurological and neuropsychiatric disorders, such as Alzheimer's disease (Wang et al., 2015), frontotemporal lobar degeneration-related syndromes (Murley et al., 2020, Murley et al., 2022, and schizophrenia (Jelen et al., 2018, Cohen et al., 2015, Duarte and Xin, 2019, as well as normal aging (Boumezbeur et al., 2010). Such spectroscopy can be predictive of task performance and the response to pharmacological interventions (Chowdhury et al., 2012, Stagg et al., 2011, McColgan et al., 2020, Gawne et al., 2020, Schmitz et al., 2017, Adams et al., 2021. Differences in neurotransmitter levels influence synaptic transmission, which in turn affects the generators of magneto-and electro-encephalographic signals (M/EEG).
Here, we characterise the relationship between physiology and neurotransmitter levels, using a hierarchical Bayesian approach . We present a method for specifying and comparing the evidence for hierarchical models (using parametric empirical Bayes) that encodes the effect of neurotransmitter levels (as measured with MRS) on different combinations of synaptic parameters in a neural mass model of MEG recordings. Our approach complements previous studies, where correlational analysis has been used to investigate relationships between spectroscopy measures and electrophysiology (Rideaux, 2021, Steel et al., 2020, Kober et al., 2001. We envisage that our method may be useful for characterising the synaptic deficits in several neuropsychiatric conditions (e.g., (Adams et al., 2021, Limongi et al., 2021).
MRS is a non-invasive neuroimaging modality used to estimate biochemical concentrations, including the neurotransmitters glutamate and gamma-aminobutyric acid (GABA) (McColgan et al., 2020, Stagg andRothman, 2013). Glutamatergic neurotransmission mechanisms include release, reuptake into astrocytes, conversion to glutamine then glutamate, and vesicular repacking (Gruetter et al., 2001, Pellerin et al., 2007, van der Graaf, 2010. In contrast, GABA cycling is predominantly neuronal where, following a release phase, there is presynaptic GABA reuptake and vesicular repackaging (Gruetter et al., 2001). The majority of GABA and glutamate is intracellular (Myers et al., 2014), but total GABA and glutamate concentrations correlate with neurophysiological features such as gamma oscillatory power and corticospinal excitability (Lally et al., 2014, Greenhouse et al., 2017.
In what follows, we show that neurotransmitter levels can be used as prior constraints on the estimation of effective synaptic connectivity, in biophysically-informed models of cortical function. For example, glutamate and GABA concentrations control the dynamics of the excitatory amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors and inhibitory GABAergic receptors respectively (Rideaux, 2020, Gruetter et al., 2001. The aim of the current work was to compare the evidence for alternate hypotheses regarding the relationships between glutamate/GABA concentrations and synaptic connectivity, quantified using non-invasive MRS recordings and M/EEG data, respectively. The first step in establishing this relationship is to infer synaptic parameters from neurophysiological observations. Since Hodgkin and Huxley (Hodgkin and Huxley, 1952), there have been many approaches to examine micro-, meso-and macroscale brain dynamics with many options for model identification (Nelson and Rinzel, 1998, Steyn-Ross and Steyn-Ross, 2010, Robinson et al., 2003, Terry et al., 2022, Jirsa et al., 2014, Deco et al., 2008. DCMs build directly on the modelling framework established by Hodgkin and Huxley, who proposed a linear dynamical system model to explain the relation between conductance dynamics and ion current dynamics. Later inclusion of nonlinear dynamics improved performance, balancing complexity with accuracy in a way that remains highly relevant to DCM (e.g., (Nelson and Rinzel, 1998)). The hypothesis testing machinery in DCM-enabling inferences on the posterior distributions of neuronal model parameters-balances model complexity and accuracy in generative (i.e., forward) models of neuroimaging data. The open source platform of DCM allows researchers from diverse disciplines to formulate hypotheses and test them within a fairly standard framework. However, innovations are still required to improve group studies, hierarchical modelling of individual differences, and clinical applications.
We use dynamic causal modelling to invert a canonical microcircuit model of resting state MEG data in healthy adults (Friston et al., 2012, Friston et al., 2003. This entails the variational Bayesian inversion of biologically-informed forward models of neurophysiological observations, under the Laplace assumption; i.e., assuming Gaussian posterior probabilities over unknown parameters (Zeidman et al., 2022. In DCM, gradient optimisation of variational free energy is used for approximating the posterior probability density over unknown model parameters, and the model evidence (i.e. marginal likelihood) . The free energy provides a lower bound on logmodel evidence, which represents the model accuracy adjusted for complexity. Model evidence associated with alternate hypotheses-about the underlying generators of dataare compared using Bayesian model reduction (BMR). This identifies the most likely explanation for the empirical data (Kass andRaftery, 1995, Friston andPenny, 2011). Crucially, BMR enables the computationally efficient evaluation of posteriors and model evidence under a reduced prior; i.e., a model specified in terms of new prior constraints Penny, 2011). We estimate the parameters and evidence for a full model of a given dataset and then use BMR to evaluate posteriors and model evidences under alternative priors , Friston and Penny, 2011. At the second (e.g., group) level, hierarchical or parametric empirical Bayes (PEB) allows one to include empirical priors of interest . Usually, second level models apply priors that are conserved over multiple participants. This means the first level corresponds to a within-subject analysis, while the second level is a between-subject or group analysis (using subject specific posterior estimates as inferred by the first level). The combined hierarchical model can then be assessed in terms of its evidence, and subjected to BMR to test different hypothesis at the within or between-subject level.
In DCM, the prior mean and covariance of unknown parameters are specified to accommodate physiological interpretability and model stability, respectively. Effectively, (informative) priors provide constraints that enable the inversion of otherwise overparameterised models (Friston et al., 2003). Using DCM, differences among individuals (or groups) can be characterised in terms of post hoc associations between connectivity parameters and clinical, pathological or cognitive measure of interest, or by differential model evidence, or conditional densities over models' probabilities (Adams et al., 2021, Stephan et al., 2009, Rae et al., 2016, Passamonti et al., 2012. However, a more principled Bayesian methodology is to incorporate between-subject variables as priors on the generative model of brain physiology. Bayesian model reduction and group-DCM can then be used to test whether such (empirical) priors increase model evidence compared to conventional (or 'flat') priors. The source of empirical priors could be demographics (e.g. age), the burden of neuropathology (e.g. from PET scanning), or multi-modal measure (e.g. MRS). In this study, the DCM (at the first and group levels) leverages the spectroscopy information to impose constraints on synaptic physiology and assess their 'goodness' in terms of model evidence. If the embedding of empirical priors into a DCM leads to higher model evidence, one can infer the importance of the measured process in the generation of observed neurophysiology. We use MRS data, noting that GABA and glutamate changes are associated with many neurological conditions.
In this paper, we compare the evidence for different hypotheses about the association between MRS estimates of neurotransmitter concentration and synaptic function. See Figure   1 for an overview of the method. For a given set of connections, we use a hierarchical optimisation (a greedy search followed by a constrained optimisation) to identify the (linear or nonlinear) function of the MRS estimates-imposing prior constraints-that maximises model evidence at the group level. We then use model reduction to determine which subsets of synaptic parameters are sensitive to individual differences in neurotransmitter concentration. This approach complements previous work by (Stephan et al., 2009, Sokolov et al., 2019 who used diffusion-weighted imaging tractography as empirical priors on connectivity parameters that are inferred from functional magnetic resonance imaging time series. However, unlike previous work, we do not assume a one-to-one relation between MRS measures and connections in the generative model. Rather, we seek evidence to find which synaptic connections are informed by the empirical measures.   (Murley et al., 2020) with permission.
In the following sections, we describe the multi-modal dataset (MEG and 7T-MRS) and the first-level DCM used to fit cross-spectral density data features. We then describe the second (group) level PEB model of the mapping between neurotransmitter concentration and synaptic physiology. We present the results of first-level and group-level DCMs, using MRS data as empirical priors. Finally, we discuss the potential applications and limitations of the method. A glossary and definitions of acronyms and variables used in this paper are provided in Tables 1 to 4.

Participants
Eleven healthy adults (age 63-77 years, five women) participated after providing written informed consent. The study had ethical approval from the Cambridge 2 Research Ethics Committee. They were free from neurological or major psychiatric disorders and took no regular medication. Each participant underwent 7T MRI and MRS (Murley et al., 2020) and task free resting-state MEG data, on separate days.

7T MRS data
Ultrahigh field magnetic resonance data were acquired using a Siemens Terra scanner at the Wolfson Brain Imaging Centre using short-echo semi-LASER sequence (repetition time/echo time = 5000/26 ms, 64 repetitions (Deelchand et al., 2015)), and VAPOR water suppression calibration (Gruetter and Tkáč, 2000). A 2cm isotropic cubic voxel was located over the right inferior frontal gyrus (RIFG). The MRS voxel was placed using conventional anatomical landmarks by the same operator (AGM) prior to each scan. Anatomical variants (such as absent diagonal branches of the Sylvian fissure) could in principle affect voxel placement accuracy. However, there was a very close overlap of the prefrontal MRS voxels across individuals ( Figure 1 of (Murley et al., 2020)). Details of the voxel placement and MRS quality assurance metrics are provided with open access in Murley et al. (2020).
The MRS voxel size was constant across participants (2cm isotropic). This voxel will contain varying partial tissue volumes, and derive from marginally different brain regions among individuals. However, we used a MRS sequence with outer voxel suppression pulses to optimise localisation and used regression analysis to correct for variation in age, sex and partial volume before modelling (see (Murley et al., 2020) for details). We used rigorous localisation standards, together with a generalized linear model to remove the effect of age, sex and correct for partial volume effects. Residual glutamate and GABA values were used for subsequent analysis (Murley et al., 2020). Note that ultrahigh field MRS (7T-MRS) can better distinguish between glutamate from within the Glx "peak", and GABA, compared to high field MRS (3T-MRS).

MEG data:
Resting state MEG data were collected during two recordings, each of five minute duration, on a different day to the MRS, with eyes closed using an Elekta Vector View system with 204 planar gradiometers and 102 magnetometers. MEG data were recorded continuously with 1000 Hz sampling rate. Participants' horizontal and vertical eye movements were recorded using bipolar electro-oculogram and electro-cardiogram electrodes. Five head position indicator coils were placed on an electroencephalography cap to track the head position.
Three fiducial points (nasion, left and right pre-auricular) and a minimum of 100 head shape points were digitised using Polhemus digitization.
The data were pre-processed using the Elekta Neuromag toolbox (Elekta Oy), with MaxFilter v2.2.12 for detection and interpolation of bad sensors, and signal space separation to remove external noise from the data and correct for head movement correction. The data were then high-pass filtered at 1 Hz, stop-band filtered around [22 to 24] Hz and [47 to 51] Hz and divided into epochs of one second duration. We used the Field trip Toolbox  for detection and removal of eye movement artefacts and discontinuities.
We applied empirical Bayesian inversion in SPM12 for source inversion and extraction of the right IFG source time series for subsequent analyses (Litvak et al., 2011). We concatenated the data of the two recordings and then divided the source data into two sets, one comprising the odd numbered epochs ('odd data') and one with the even numbered epochs ('even data') for each participant to assess the reliability of the results . We separately averaged power spectral responses of the odd and even datasets. The two power spectral densities (referred to as Odd PSDs and Even PSDs) are used as the data features in the DCM of cross-spectral density.

First level analysis using dynamic causal modelling of resting states MEG data
To infer the neurophysiological parameters generating observed resting state MEG data, we used "DCM for cross spectral density" in SPM12 (Friston et al., 2012, Moran et al., 2007, Moran et al., 2011. Spectral features in the resting state MEG were used to infer synpatic parameters of a biophysically-informed neuronal mass model, together with model evidence.
DCM for cross spectral density assumes that recorded electrophysiological oscillations are due to finite responses of neuronal dynamics, under endogenous random fluctuations (Basar et al., 2012, Haken, 1977. We used a conductance based neuronal mass model as shown in Figure    Each population is represented by a Morris-Lecar model (Moran et al., 2013). The dynamic of each population is govened by stochastic differential equations that emulate the dynamics of pre/post-synaptic potentials, firing rates and membrane conductances. In a typical neruonal population, the dynamics of membrain potentials, , and conductances of an ion channel, * , are govered by the following equations: In equation 1, is the membrane capacitance, is random endogenous fluctuation, is passive leak current where its conductance is constant, * is conductance associated with ion channels/receptors with time constant * , . The term in the equation 1 is the non-negative afferent presynaptic firing from population , which is scaled by afferent intrinsic connectvity .
The scaled presynaptic firing rates are a proxy for neurotransmitter levels measured by MRS.
Therefore, we hypothesise a relation between the MRS glutamate and GABA measures and excitatory and inhibitory connections, which scale presynaptic actvity. In other words, the MRS glutamate and GABA can be used as empirical priors on excitatory/inhibitory connections in the generative model.
Mathematically the temporal dynamics of the conductance based model can be written in the canonical form of a dynamical system as follows (Friston et al., 2012):  Table 2.
One can approximate the dynamics of equation (2) with the (first order) linearised model ̇= ( * ) + ( * denotes the Jacobian at * ) with the spectral response: Given the spectral response of the system, the generative model of neuronal activity, ( ), can be calculated as follows: In equation 4, ( , ) represents the spectrum of the observation noise, which is a sum of common and source specific noise. The spectral response in sensor space can also be generated by inclusion of a forward electromagnetic model into equation 4, which is denoted by . ( ) ( is the gain and is the head model), as follows: In equation 5, ( ) is the cross spectra of the MEG data and ~(0, ) is a random effect (with unkown covariance). Because we perform the DCM analysis in source space, the forward electromagnetic model in this case just applies a scaling parameter. In the DCM, we do not pre-define which neuronal population(s) contribute to MEG data. Instead, we estimated the degree to which each population (e.g. superficial, deep pyramidal or inhibitory population) contributed to the generation of the MEG data (see Pereira et al 2021 for details specific to conductance based models). Recently clinical and interventional applications support this approach (Adams et al., 2021, Gilbert et al., 2016, Symmonds et al., 2018. The unkown parameters in the DCM are specified as log-scale values. This means that the parameter vector in DCM is a random variable, which is given by = 0 exp(̂). Here, 0 is a biologically informed scaling for the parameter and ̂ is a random variable, with prior normal density ̂0 ~(0, ∑ ) of zero mean and covariance ∑ . This expresses the belief about the range over which parameters can vary. Note that scale parameters of this sort cannot be less than zero; for example, distances and lengths or rate and time constants or precisions and variances.

Second level analysis: group DCMs informed by empirical priors
In many translational neuroscience studies, the goal is to test hypotheses at the group level, where subject-specific information is given by summary statistics. In PEB, this subject-specific information is the posterior over key model parameters at the first level (Zeidman et al., 2019a, Zeidman et al., 2019b. At the second level the objective is to test whether synaptic connections depend on between subject variables, such as age, disease-severity, genetics, diffusion tensor imaging and MRS . In other words, PEB seeks to explain intersubject variability on one or more first level model parameters. Mathematically, we denote a vector of model parameters at the first level DCM, over all participants, by a column vector (1) (superscript '1' denotes the first level analysis) with dimension × 1 ( number of participants and is the number of parameters for each participant). Then the generative model at the second level (i.e., random effects on the parameters) is given by : Where ∈ × is the design matrix with ≥ 1 covariant (the first regressor of X is always equal to one and reflects the group mean) and a column vector (2) (superscript '2' denotes the second level) contains the second level parameters. The symbol ⊗ is the Kronecker product and is the × identity matrix. The random effects have a Gaussian distribution (2)~( 0, (2) ) (where (2) is precision matrix or inverse of the covariance). The precision matrix is parameterised with a single (hyper-precision) parameter, γ, as follows ): (2) = ⊗ ( 0 + − 1 ) In equation 7, 0 ∈ × is the lower bound on the precision, defined with a small positive value). The (hyper)parameter, , scales a precision matrix 1 ∈ × , which is (by default) 16 times the prior precision of the group mean (Zeidman et al., 2019b): this hyperprior ensures that random effects are small in relation to prior uncertainty about the parameter in question.
The critical question is whether an alternative prior increases or decreases model evidence . This is an optimisation problem, where the objective function is the free energy of the PEB model. Inversion of the hierarchical model is computationally expensive if all parameters of the first level need to be inferred every time the second level parameters change (Raman et al., 2016). However, BMR can be used to re-evaluate first level posteriors under updated second level parameters . This significantly improves the efficiency of the inversion for PEB models.

Investigating relation between MRS and synaptic connections:
2.4.1 Problem setting: We test whether excitatory and inhibitory synaptic connections depend on (i.e., are functions of) glutamate and GABA measures, respectively. We used PEB to specify and compare the evidence for different functions of MRS measures, as illustrated in Figure 3.  To identify which synaptic parameters were sensitive to neurotransmitter levels, we grouped the inferred synaptic connections into excitatory and inhibitory subsets. In detail, for the To cross-validate the results, we separately tested the relationship between inferred synaptic parameters and MRS by splitting the MEG data into odd and even numbered epochs (odd-Data and even-Data). The ensuing procedure is shown in Figure 4.

Functional form of the empirical MRS priors:
The functional form of the MRS mapping is not known a priori. We therefore limit the search space over the MRS transformations to continuous and monotonic polynomial maps and sigmoid nonlinear functions, which cover a wide range of linear and nonlinear forms. The class of polynomials provides an approximation (Bishop, 2006) to any nonlinear monotonic form (Spivak, 2020): Here, are unknown hyperparameters that need to be inferred by finding the values that Animal experiments suggest that changes in GABA and glutamate concentrations are related to neuronal responses via a sigmoid relationship (Benardete and Kriegstein, 2002, Dyke et al., 2017, Chebib et al., 2009). The general form of a sigmoid nonlinearity can be parameterised as follows: Here, the unknown hyperparameters 1 , 2 are the slope and threshold parameters of the sigmoid transformation, respectively.

Hyperparameter estimation:
To estimate unknown hyperparameters, we first sorted the MRS measures from small to large values, which defined an interval of the real line (called the domain of the MRS data). We Each parameter in the map, , varies with resolution (or step-size) Δ . Note that γ and Γ have the same functional forms but with continuous and discretised domains. Having established the appropriate ranges (see Table 4), we used a hierarchical scheme to estimate the parameters of the map, Γ ( 3 Results: 3.1 DCM for MEG: The spectral activity from both data sets per participant (even and odd data) are used as the data features for a conventional "DCM cross spectral density" to infer the parameters of the conductance based canonical microcircuit. The predicted and observed data are provided in the supplementary material Figure 13. Over all subjects, the mean variance of observed CSDs explained by the predicted CSD was 98%. The comparison of predicted and observed spectral data shows that the synaptic parameter estimates can replicate the spectral patterns of both even and odd data; with an example shown in Figure 5 (see the supplementary material for the remaining comparisons). The correlation between the free energies associated with DCM inversions of odd/even PSD data is shown in Figure 5.  Figure 13.

Panel (b) shows the correlation between free energies of odd and even datasets over participants (points in the plot). This correlation shows a moderately good level of agreement between the odd and even data across all participants. Panel (c) illustrates correlations between inferred synaptic connection between odd and even epochs. The inhibitory connections are shown in red and excitatory connections are shown in black.
To confirm the predictive validity of DCM, we calculated the correlation between the inferred synaptic parameters from odd and even PSDs data sets: see Figure 5. Most parameters were estimated reliably (figure 5). There were exceptions, such as the excitatory connections of the deep pyramidal and inhibitory connections and superficial pyramidal cells and inhibitory cells.
The low correlation between estimates of some parameters is not surprising for complex nonlinear models with conditional dependencies among parameter estimates (Litvak et al., 2019). It implies that different combinations of synaptic parameters may generate similar physiological data. This means hypotheses may therefore be better tested in terms of model comparison, rather than focussing on maximum, a posteriori parameter estimates of single parameters (Rowe et al., 2010). In the following section, we use the reliable connectionswith a correlation of greater than one half-to examine the effect of neurotransmitter concentrations as measured by MRS.

Associations between synaptic parameters and MRS measures:
The most likely relationship between the synaptic parameters and neurotransmitter levels was identified through the hyperparameter optimisation. We grouped the inhibitory connections into four self-connections and three inter-lamina connections. By grouping the parameters in this way, we ask whether MRS GABA influences recurrent intra-laminar connectivity, or inhibition between layered populations. We assess the evidence for GABA effects, as mediated through second order polynomial and sigmoid transformations of the MRS measures. We systematically varied the hyperparameters of the sigmoid function for each subset of self-inhibitory synaptic connections and compared the resulting model evidence. Figure 6 shows the consistency of the resulting free energy over odd and even epoch datasets.
The correlation plot between the free energies of odd and even data models shows that the ranking of model evidence is consistent, where each 'model' is a hypothesis about how MRS data supplies prior constraints on synaptic parameters. Such agreement is an important validation of hypothesis testing, based on selecting models with the greatest evidence. Using PEB to make inferences at the between-subject level inherits this validity. The high reliability suggests that both group inversions converge to the same (global) minima. For each subset of self-inhibitory connections, the sigmoid MRS GABA hyperpriors maximised the model evidence. The winning models for each subset were then examined to determine which MRS mapping is most likely for each subset of self-inhibitory connections. As shown in Figure 7, priors using a sigmoid transformation of GABA concentration provided the most likely account of intersubject variation in synaptic connectivity; specifically, the recurrent connections or self-inhibition of superficial and deep pyramidal cells.    We repeated this procedure using a second order polynomial (the first order polynomial is contained within this model space) as the functional form (see supplementary material). We then compared the evidence of the winning models in each analysis; namely, the evidence for the mapping between (sigmoid MRS, self-inhibition), (sigmoid MRS, inter regions), (polynomial MRS, self-inhibition) and (polynomial MRS, inter regions), as shown in Figure 10.
The results suggest that the sigmoid transformation is the most likely functional form to explain intersubject variability in the inhibitory recurrent (self) connections in superficial and deep layers. 3.2.2 MRS glutamate data as an empirical prior for excitatory synaptic parameters: We performed the proposed analysis to identify the relationship between glutamate concentration and excitatory synaptic connections. We considered different combinations of the excitatory connections over sigmoid and second order polynomial functional forms of the relationship between MRS and synaptic parameters. The correlations of models' free energy over odd and even data is shown in Figure 11. The winning model over the complete search space confirms that MRS glutamate is linked to excitatory connections from superficial to deep layers and superficial to inhibitory interneurons (model 3).

Discussion:
We present a framework for specifying and comparing hypotheses of the relationship between synaptic physiology and neurotransmitter levels, based on a combination of restingstate MEG and MRS. This enables one to specify and compare the effect of biomarkers of inter-subject variation as informative priors on synaptic parameters in a time efficient and Bayes-optimal manner. The method is applicable to generative models of evoked and restingstate time series, recorded by magnetoencephalography, electroencephalography or functional magnetic resonance imaging techniques. In principle, the individual differences specified here by MRS could be replaced by other imaging modalities, or regional neuropathology assays.
We illustrate how the method can be used to test, non-invasively, hypotheses about the relationship between synaptic physiology and neurotransmitter concentrations in humans.
This goes beyond the correlation between M/EEG data features and MRS estimates of neurotransmitter concentration (McColgan et al., 2020), in using formal, generative models of the neurophysiological observations. It is computationally efficient, using first level DCM once and re-evaluating the DCM parameters and model evidence analytically using Bayesian model reduction. The model evidence is thereby compared between models with and without the MRS priors, in contrast to the reinversion of the DCMs for alternative priors (Stephan et al., 2009).
We provide evidence that a non-linear transform of the MRS GABA measures offers the best explanation for inter-subject variability in inhibitory recurrent (i.e., self) connectivity of superficial and deep neuronal populations. In addition, a sigmoid transform of the MRS glutamate provides the best explanation of inter-subject variations in excitatory connections from superficial to deep layers. The identification of a sigmoid form for the MRS transform is interesting, because experimental changes of GABA and glutamate concentrations exhibit a sigmoid relationship to neuronal responses, though postsynaptic gains (Benardete and Kriegstein, 2002, Dyke et al., 2017, Chebib et al., 2009. The relationship between superficial and deep layer connectivity is especially relevant to studies of neuropsychiatric and neurological disorders with abnormalities in GABA and glutamate, e.g., schizophrenia, movement disorders and dementia (Adams et al., 2021, Murley et al., 2020, Murley et al., 2022. Using MRS data as an empirical prior on the synaptic connections may provide valuable information about the impact of neuropathology on synaptic function and the response to treatment. There are several limitations to this study. We focus on the relationship between neurochemical concentrations and synaptic physiology in one cortical region. The choice of the right inferior frontal gyrus was motivated by our interest in frontal lobe function and its impairment in frontotemporal lobar degeneration (Eliasova et al., 2014, Hughes et al., 2015, Murley et al., 2020, Murley et al., 2022 and other neuropsychiatric disorders. However, the method can be applied to multiple regions. For instance, one could combine neurophysiology and spectroscopy from principal nodes in the default mode or salience networks; and test whether MRS data at different regions are associated with extrinsic (between-source) and intrinsic (within-source) connections. Such an approach could test whether local neurotransmitter concentrations in one source influence the activity and connectivity of other sources. The spatial resolution of MRS data is limited. The average neurotransmitter concentrations are captured over multiple cortical columns. Here, the GABA estimate is used as a marker of between-subject differences in neurotransmitters for the lateral frontal cortex.
We assume that the voxel-wise estimate approximates the neurochemical concentrations in the neurophysiological source used to extract time series for DCM inversion, noting that the source lies within the MRS voxel (Murley et al., 2020).
Our study size was modest, although similar to (Stephan et al., 2009)  Vector of inferred parameters at first level DCM over all subjects with dimensions × 1 ( number of subjects and number of inferred parameters in each DCM (2) Second level parameters with dimensions × 1 ( number of covariates) Γ Scalar function of MRS data which is parametrised with a vector of parameters 1 Vector of inferred parameters at first level DCM from all subjects that are influenced by MRS data.
(2) Precision matrix at the second level with dimensions × ( number of parameters at the first level)

Variable Description
Exogenous input (scalar) to membrane equation and conductance equation within each population. Mean depolarisation of a neuronal population. It is a scalar variable per each population.

( )
The neuronal firing rate (scalar function) -a sigmoid squashing function of depolarisation Lead field vector mapping from (neuronal) states to measured (electrophysiological) responses ( ), ( ), ( ) Spectral density of (neuronal) state fluctuations, observation noise and measurement, respectively. These are vector valued functions. System Jacobian or derivative of system flow with respect to (neuronal) states (matrix valued function).

( ) = [ ( )]
First order kernel mapping from inputs to responses; c.f., an impulse response function of time. This is the Fourier transform of the transfer function.

( ) = [ ( )]
Transfer function of frequency, modulating the power of endogenous neuronal fluctuations to produce a (cross spectral density) response. This is the Fourier transform of the kernel.     Figure 14 illustrates model comparison of second level PEB associated with self-inhibitory connections (top) and inter-regional connections (bottom plot) given the class of first and second order polynomial maps for MRS data as their regressors. Figure 15 illustrates the definition of models in Figures 7, 9 and 11. Subject 1 Subject 2 Subject 4 Subject 3 Figure 13. First level DCM predicted vs observed data for each subject even and odd PSD data. The variance explained across all data (even/odd) is 98%.