Ketosis regulates K+ ion channels, strengthening brain-wide signaling disrupted by age

Endogenous ketone bodies, including D-β-hydroxybuty-rate (D-βHB), are produced by the liver and can be utilized by cells as fuel when glucose is not readily available (Krebs, 1960). Accumulating evidence suggests that ketone metabolism may confer important neurological benefits (Barañano & Hartman, 2008; Henderson et al., 2009; Maalouf et al., 2007; Vanltallie et al., 2005). D-βHB has been found to increase ATP production and cardiac efficiency from 10% to 24% when added to a perfusion of glucose (Sato et al., 1995). Moreover, ketone uptake bypasses the insulin-dependent glucose transporter GLUT4, and thus can be metabolized as fuel even after neurons become insulin resistant (Sapolsky, 1986).

Insulin-independent metabolic pathways may be critical in the context of brain aging, as Type 2 diabetes mellitus and its associated decrease in GLUT4-dependent neuronal glucose utilization are linked to age-related brain hypometabolism (Baker et al., 2011; Soares et al., 2019) and cognitive decline (Antal et al., 2022; Beeri et al., 2004). Even in the face of impaired glucose metabolism, aging brains can still metabolize ketone bodies (Cunnane et al., 2016). Furthermore, ketones also influence cerebral metabolism by modulating glucose uptake in astrocytes (Valdebenito et al., 2016) and directly affecting neuronal excitability (Cauli et al., 2023; Karagiannis et al., 2021). Thus, supplementing ketone bodies as an alternative fuel source may have the potential to slow or arrest neurodegeneration (Zilberter & Zilberter, 2017).

Recent findings provide evidence supporting the ability of exogenous D-βHB to ameliorate mechanisms and biomarkers associated with brain aging. At the spiking-neuron scale, insulin and ketone bodies modulate neuronal excitability through the regulation of K⁺ ion gradients (Ma et al., 2007; Sweeney & Klip, 1998). At the circuit scale, direct application of D-βHB to the mouse hippocampal CA3-CA1 circuit reverses insulin resistance-induced deficits in neuronal excitability and axon conduction velocity, showing improvements that exceeded even baseline (control) values (Kula et al., 2024). Finally, at the whole-brain scale, aging has been linked to the destabilization (Mujica-Parodi et al., 2020) and desynchronization (Weistuch et al., 2021) of brain networks identified in resting-state functional magnetic resonance imaging (rsfMRI). Ketosis, induced by either diet or exogenous D-βHB administration, improved these measures of brain-wide coordination, even in young, healthy adults. Despite the theoretical relationships suggested between these independent results, the direct linking across scales—of mechanism with its emergent effects—has yet to be tested due to the technical constraints inherent in each experimental approach.

To address these challenges, we employ a novel experimental-computational approach for evaluating the contribution of individual biophysical components on brain-wide coordinated behaviors. To link fMRI-derived brain network instability and synchrony specifically to the postsynaptic potentials of pyramidal neurons, we established an EEG-derived biomarker for brain aging by analyzing a publicly available Lifespan resting-state EEG (rsEEG) dataset (ages 20-80 years, $N=201$⁠). To test how this biomarker changes with metabolism, we conducted a new Metabolic rsEEG experiment (⁠$μage=26.9$$±$$11.2$ years, $N=36$⁠) (Table S1) (Fig. 1A). Administering individually weight-dosed D-βHB versus calorically-matched glucose orally, we measured the effects of fuel type on synchrony (Fig. 1B) and brain network instability (Fig. 1C). To identify the set of mechanistic parameters consistent with our data, we then used our EEG results as emergent constraints on the Larter-Breakspear model, a conductance-based multiscale neural mass model (Breakspear et al., 2003). By simulating this model across a full range of key mechanistic parameters and simultaneously testing $38$ hypotheses in parallel (three types of trends, eight parameters), we were able to test whether measured differences in EEG-derived neural synchrony at the macroscale were consistent with effects expected following modulation of K⁺ at the neuronal micro-scale (Fig. 1D). Equally important, however, by testing alternative plausible hypotheses in parallel, we were able to assess not only the sensitivity but also the specificity of our results.

To determine whether stabilization of brain networks as modulated by fuel source seen in rsfMRI (Mujica-Parodi et al., 2020; Weistuch et al., 2021) implicates neuronal synchrony, we conducted resting-state EEG scans on a cohort of healthy adults (⁠$N=36$⁠, $μage=26.9±11.2$ years, 18 female) who were tested under both a ketotic (ketone burning) and glycolytic (glucose burning) condition. Inclusion/exclusion criteria were screened using a survey completed before enrollment in the study. Potential participants were excluded for any of the following reasons: chronic usage of alcohol (7 drinks/week for women, 14 drinks/week for men), recreational drug use, use of psychotropic medication within the past 30 days, use of medications that affect glucose and/or insulin utilization, a history of kidney disease, heart attack, stroke, myxedema, epilepsy, dementia, or other neurological disorders, difficulty when swallowing, current pregnancy or breastfeeding, or inability to provide informed consent. Scanning occurred at Stony Brook University’s Health Science Center. The study was approved by the Institutional Review Board of Stony Brook University (IRB2021-00018), and all participants provided informed consent.

All subjects were tested twice (1–14 days apart), both times following an overnight fast (subjects were instructed to eat no food for at least 8 hours before testing but were allowed unrestricted water). Following a baseline resting-state EEG scan, subjects drank either of two fuel sources. In the ketotic session, subjects drank a ketone sports drink, deltaG® Sports Supplement (TdeltaS Ltd, Thame, UK), dosed at 395 mg/kg. During the glycolytic session, the same subjects drank a bolus of glucose (Glucose Tolerance Test Beverages, Fisher Scientific, Inc.; Hampton NH) calorie-matched to the ketone bolus. The order of the bolus (whether a subject received a glucose bolus first or a ketone bolus first) was pseudo-randomized, with approximately half of all female subjects (⁠$N=9$⁠) and exactly half of all male subjects (⁠$N=8$⁠) drinking the glucose bolus during the first scanning day. The resting-state EEG scans were then repeated 30 minutes after the administration of the bolus, as prior experiments using magnetic resonance spectroscopy (MRS) showed peak glucose and ketone concentration in the brain approximately 30 minutes after consumption of the bolus (see Fig. S2 in Mujica-Parodi et al. (2020)). Blood glucose and ketone levels were measured three times throughout the experiment: at baseline, 10 minutes following the bolus, and 62 minutes following the bolus using a Precision Xtra Blood Glucose & Ketone Monitoring System (Abbott Laboratories). Our data analyses quantify network reorganization and neural phase synchrony changes in response to changing energy constraints (i.e., cognitive demand, fuel).

During the resting-state portion of the EEG scan, subjects underwent a total of 16 blocks, each lasting 60 seconds, with 8 EO blocks and 8 EC blocks. The blocks were interleaved. During the EO blocks, subjects fixated on a white orienting cross on a black background. Before the resting-state scan, subjects were instructed not to blink to minimize ocular artifacts and to keep motion to a minimum. The resting-state stimulus (a white cross in the center of a black background) was presented on a computer screen placed in front of the seated subject using PsychoPy v3.0 (Peirce, 2007). All data were collected in a shielded, dark, soundproofed Faraday cage using the ActiveTwo Biosemi™ electrode system from 65 (64 scalp, 1 ocular) electrodes arranged according to the international standard 10–20 system (Oostenveld & Praamstra, 2001) at a sampling frequency of 4096 Hz. The ocular (VOEG) electrode was placed below the left eye. Our experimental and pre-processing designs, especially those of the resting-state EEG scans, were modeled after the paradigm used by Leipzig’s LEMON rsEEG dataset group (Babayan et al., 2019) to be able to minimize confounding factors when directly comparing results between our experiment and this large-scale neuroimaging dataset.

The EEG pre-processing was performed using EEGLab (version 2020.0) (Delorme & Makeig, 2004). Full resting-state data were downsampled from 4096 Hz to 512 Hz and bandpass filtered between 0.1 and 40 Hz using a Hamming-windowed FIR filter. The data were then separated into the eyes open (EO) and eyes closed (EC) conditions. These two conditions were pre-processed separately from this point on due to the differences in ocular artifacts in each condition; however, the pre-processing steps performed were identical. Bad channels were removed, and noisy portions of data were identified and removed using EEGLab’s Artifact Subspace Reconstruction (ASR) algorithm. Independent component analysis (ICA) was performed on the data using the infomax algorithm in EEGLab (runica), and non-neural components of the time-series were identified and removed using ICLabel (Pion-Tonachini et al., 2019). The reference was then set to average. The data were separated into frequency bands and time-series extracted using MNE-Python (version 0.21.1) (Gramfort et al., 2013).

Lifespan rsEEG data from the Leipzig LEMON dataset were obtained in already pre-processed form. Raw data were downsampled from 2500 Hz to 250 Hz and bandpass filtered between 1 and 45 Hz using an 8th order Butterworth filter. The data were separated into EO and EC conditions for subsequent pre-processing. Outlier channels were rejected, and data intervals containing high peak-to-peak fluctuations or high-frequency noise were identified and removed by visual inspection. Data dimensionality was reduced using principal component analysis before the use of independent component analysis (ICA) to identify and remove components reflecting eye- or heartbeat-related artifacts. Further pre-processing details may be found in the dataset documentation (Babayan et al., 2019).

Global Field Synchronization (GFS) was first introduced by König et al. (2001) to estimate differences in functional connectivity of brain processes in EEG frequency bands between a population of neuroleptic-naive schizophrenic patients and healthy controls. In contrast to other measures of phase synchrony such as phase-locking value, which can only measure synchrony between two time-series, GFS is a global measure of neural phase synchrony that does not rely on the a priori selection of brain regions to be studied. When applied to EEG data, GFS has the added benefit of being reference-independent and more easily interpretable without the use of source models (Michels et al., 2012). Further research using the measure has found changes in GFS in those with Alzheimer’s disease and mild cognitive impairment (König et al., 2005), during REM sleep (Achermann et al., 2016), and during general anesthesia (Nicolaou & Georgiou, 2014).

To calculate GFS values for EEG, the data were first pre-processed. Following pre-processing, sensor-level EEG time-series were divided into non-overlapping, consecutive 2-second epochs. Each epoch was frequency transformed using a fast Fourier transform (FFT, Tukey window, taper size $α=0.2$⁠), which yields the real and complex component of the signal of each electrode for each frequency value (1–40 Hz, step size = 0.1 Hz). These components are then plotted as vectors on the complex plane, with the magnitude of the vector representing the power of the signal at that frequency, the angle of the vector as measured from the real axis representing the phase, and the vector origin representing the reference of all EEG signals. Subsequently, a scatter plot of vector endpoints in the complex plane is generated for each frequency value, a representative of which can be seen in Figure 1B. The more these endpoints approximate a straight line, the more the signals are in phase or anti-phase. The more scattered the endpoints, the less the signals are in phase or anti-phase. To quantify the shape of the cloud formed by the vector endpoints, the points are entered into a two-dimensional principal component analysis, as principal components are orthogonal by definition. The resulting GFS value per epoch for a particular frequency $f$ is then determined by calculating the ratio of the eigenvalues (⁠$λ1(f)$ and $λ2(f)$⁠) of these two principal components, as expressed in Equation 1.

Finally, the overall GFS value for a particular frequency is obtained by taking the mean of the GFS values of all consecutive epochs for that frequency, creating a spectrum as shown in Figure 1B. To examine differences in GFS across age (using the Leipzig LEMON rsEEG dataset) and condition (ketotic vs. glycolytic), overall GFS values were categorized into their corresponding, classically-accepted frequency bands (da Silva, 2013), with endpoints adjusted to accommodate the filtering performed in pre-processing: $δ$ (1–3.5 Hz), $θ$ (4–7.5 Hz), $α$ (8–13 Hz), $β$ (14–30 Hz), low $γ$ (30–40 Hz), and broad (1–40 Hz).

We used a conductance-based neural mass model (Breakspear et al., 2003) to test the effect of ion gradient dynamics and excitatory-inhibitory subpopulation coupling on GFS. This neural mass model has been previously validated to capture properties of fMRI and EEG resting-state activity (Endo et al., 2020) while incorporating both local ion dynamics and interregional connectivity. Following the mathematics of prior work (Endo et al., 2020), we built a 78-region whole-brain simulation using Neuroblox, a Julia library optimized for high-performance computing of dynamical brain circuit models (https://github.com/Neuroblox/Neuroblox.jl). Within the model, we varied ion (Na⁺, K⁺, Ca²⁺) gradients and conductances, intraregional connectivity (excitatory-excitatory/inhibitory), and interregional connectivity across a range of biophysically plausible values (Chesebro et al., 2023; Endo et al., 2020; Roberts et al., 2019). Ion parameters were changed because metabolic changes have been shown to alter ion dynamics in microscale experiments (Baeza-Lehnert et al., 2019). Excitatory-excitatory/inhibitory and global connectivity parameters were varied as a control to ensure that the effects in the model do not arise from just any changes to inter- or intra-regional connectivity. The ranges of the values changed can be found in Table 1. The whole-brain simulation for each parameter set was repeated with 10 different sets of initial conditions to sample across the distribution of simulation outcomes. The simulated membrane potentials of excitatory neurons were averaged within each region to generate 78 neural mass signals, which were transformed into EEG signals through multiplication with an average lead field (63 EEG channels $×$ 78 model ROIs) generated by Endo et al. (2020). Each 10-minute simulation was sampled at a rate of 1000 Hz and took $∼$2 minutes to generate. The first 1 minute of each simulated EEG signal was discarded before calculating GFS to allow the simulation to equilibrate.

Brain network instability (described for use with fMRI data in Mujica-Parodi et al. (2020)) is a measure used to describe the persistence of brain networks over time. It can be considered a measure of dynamic functional connectivity (Hutchison et al., 2013), or a quantification of the frequency of switching between metastable brain states over various temporal scales (Roberts et al., 2019). To calculate brain network instability values for EEG, the data were first pre-processed as described above. Following pre-processing, sensor-level EEG time-series were divided into non-overlapping, consecutive epochs, or windows, with a 10 second window chosen for the results shown in Figures 2C and 2D. Within each window $Wi$ an all-to-all, signed correlation matrix between all time-series is calculated, meaning these resulting matrices have a size of $n×n$⁠, with $n$ being the number of clean channels remaining after pre-processing. Distance (in units of window size) between window pairs is decided by a value of $τ$ chosen for the instability calculation. Instability is then calculated for each possible value of $τ$ by taking the Frobenius, or L2, norm of the difference in correlation matrices of window pairs, and then taking the average of all norms. In the form of an equation, instability is written as

where $τmax$ is the number of windows available, defined as the length of the time-series divided by the window size rounded down to the nearest multiple of the window size. For example, an EEG recording with a length of 55 seconds divided into window sizes of 10 seconds means $τmax=5$⁠. Computing instability for $τ=2$ for this example recording would yield

When calculating instability, the window sizes within which correlation matrices are calculated and eventually subtracted from one another set a natural limit to the timescale at which changes in global network connectivity can be seen. As the temporal resolution of EEG is notably greater than that of fMRI, the window sizes used in the calculation of instability for the majority of this work are smaller than the windows used for our previous work (10 seconds as compared to 24 seconds) as more time points are available for calculation. By varying the window sizes used to calculate network instability, we determined that brain networks destabilize as a function of age only for networks whose detectable correlation differences persist for $≥$10 seconds (Fig. S2), and as such made the choice of window size to be 10 seconds for subsequent analyses. Furthermore, to ensure the increase in instability with age seen in the Leipzig LEMON dataset was not due to motion (as subject motion tends to increase with age (Seto et al., 2001)), we examined the relationship between the mean frame displacements (in mm) measured during the subjects’ corresponding rsfMRI scan and their brain network instability (⁠$τ=1$⁠, window size = 10 seconds) within the broad (1-40 Hz) frequency band of the subjects’ EC, rsEEG scan. We found no correlation between these two variables.

Two-sided independent t-tests were used to determine whether mean GFS and instability differed significantly between the younger and older cohorts of the Leipzig LEMON rsEEG dataset within frequency bands (for GFS, Fig. 2A) and within $τ$ values (for instability, Fig. 2C). Two-sided paired t-tests were used to determine whether mean GFS and instability differed significantly between the glucose and ketone conditions within frequency bands (for GFS, Fig. 2B) and within $τ$ values (for instability, Fig. 2D). As instability values per $τ$ (Fig. 2C, D) are by definition not independent, no multiple-comparisons correction was applied to these results. All other statistics (Fig. 2A, B) were corrected for multiple comparisons using the Benjamini–Hochberg procedure with a significance level set at $α$ = 0.05.

Using the physiological and mathematical boundary conditions discussed in Chesebro et al. (2023) to inform the model, the model equations are constructed as a system of three variables: mean excitatory membrane voltage $V(t)$⁠, mean inhibitory membrane voltage $Z(t)$⁠, and the proportion of potassium channels open at a given time $W(t)$⁠. Note that while $V(t)$⁠, $Z(t)$⁠, and $W(t)$ are all time-dependent, we omit this dependence in the following equations for ease of reading. Given this understanding, the neural mass model is defined as:

In these equations, $QV$ and $QZ$ are the mean firing rates for excitatory and inhibitory cell populations, respectively. These are computed as

The individual ion channel gating functions (⁠$mNa$⁠, $mK$ and $mCa$⁠) take the form

where $mion$ is the fraction of voltage-dependent channels open at any given time. Default values and descriptions for all constants in these equations are given in Table 1. Note that parameter values are unit-less to scale to a reasonable modeling range (i.e., $V,Q∈(−1,1)$ and $W∈(0,1)$⁠), and the integration time step $dt$ is in milliseconds.

We note that the three ions of interest are modeled in three different manners. Sodium serves as the dominant shape determinant of the neural mass spiking activity as it has the highest net positive conductance coupled to its ion channel gating function. Calcium serves as a secondary support of the spiking activity, providing some of the amplitude to the spiking activity. However, because the calcium gating function is also coupled to the excitatory population firing rate and the ratio of NMDA/AMPA receptors, calcium more importantly provides feedback to the neural mass firing rate. Finally, due to its more detailed modeling as a separate differential equation, potassium plays a unique role in determining the frequency of spiking activity (at larger reversal potentials) and the duration of the refractory period (at smaller reversal potentials). As a consequence of this extra modeling step, potassium also has a different ion gradient landscape than sodium and calcium.

While the excitatory (pyramidal) cell population $V(t)$ is modeled using the ion dynamics described above, the inhibitory population $Z(t)$ is a purely phenomenological model, receiving only the excitatory input via $aei$ and a background current via $ani$⁠. Although this serves to model the relationship between inhibitory interneurons and the excitatory pyramidal cells (as in Larter et al. (1999)), it does imply a caveat when interpreting the effects of altered ion gradients.

Since the model is a hybrid of a biophysically detailed excitatory neuron population and a phenomenological inhibitory population, claims regarding how closely this model resembles true biological neurons are necessarily constrained. However, the advantage of this neural mass model is that it produces physiologically interesting dynamics (e.g., burst-spiking) that are more common in next-generation neural mass models (Taher et al., 2022) than in traditional (e.g., Wilson-Cowan) oscillatory models.

Equations 5–7 describe a single neural mass comprising two subnetworks. Coupling between pairs of neural masses (⁠$i$ and $j$⁠) can also be achieved through connection terms:

Here, $c$ is the coupling constant, $QVj$ is the mean excitatory firing rate of region $j$⁠, and $ui,j$ is the strength of connection between regions $i$ and $j$⁠. To ensure that overall input current is approximately constant, the balancing between interregional and self-coupling takes the form of competitive agonism, where $c$ is the weight of interregional coupling and $(1−c)$ is self-coupling. The associated multi-regional neural mass model equations are then given by:

In this work, we use Equations 12–14, varying the microscale parameters (Table 1) therein to observe the effects on macroscale neural synchrony. Following prior work (Endo et al., 2020), we computed a 78 region whole-brain model. To determine connectivity between regions, we used a DTI-derived structural connectivity map averaged across healthy individuals. This connectivity has previously been shown to produce simulations with the conductance-based neural mass model used in this work that are consistent with both EEG and fMRI scale measurements (Endo et al., 2020). This map is the same as in Endo et al. (2020), having been provided by the authors as supplemental data. This is the standard approach employed by toolboxes such as The Virtual Brain (Sanz Leon et al., 2013).

Using the Leipzig LEMON rsEEG dataset, a comparison of Global Field Synchronization (GFS, see Methods) between the younger (aged 20 to 40 years, N = 138) and older (aged 55 to 80 years, N = 63) cohorts within classically defined frequency bands (⁠$δ$ (1–3.5 Hz), $θ$ (4–7.5 Hz), $α$ (8–13 Hz), $β$ (14–30 Hz), low $γ$ (30–40 Hz), and broad (1–40 Hz), (da Silva, 2013)) showed age to decrease neural synchrony in the broad, $θ$⁠, $α$⁠, and low $γ$ bands (Fig. 2A, Table 2). Likewise, older brains showed destabilization of their networks (Fig. 2C).

Changes in GFS values (⁠$Δ$GFS) for the metabolic study cohort (N = 36, $μage=$ 26.9 $±$ 11.2 years) were calculated by subtracting the average fasting condition GFS from the average experimental (glycolytic or ketotic) condition GFS within each band. Comparing the pre- and post-bolus synchrony values showed decreased synchrony after administration of the glucose bolus and increased synchrony after administration of the ketone bolus for all frequency bands in the rsEEG eyes-open (EO) condition (Fig. 2B, Table 2). Following glucose challenge, brain networks destabilized; in contrast, following ketone challenge, brain networks stabilized (Fig. 2D) in rsEEG. Both results were observed in the broad (1–40 Hz) frequency band of eyes-closed (EC) rsEEG.

We next sought to identify the mechanistic basis for age and metabolism-related changes observed in rsEEG GFS. To do so, we used a bottom-up approach by varying microscale parameters in the neural mass model to generate 2082 simulated rsEEG signals, from which we computed GFS and instability. We observed high sensitivity of GFS to changes in all (Ca²⁺, K⁺, and Na⁺) Nernst potentials, Ca²⁺ and K⁺ channel conductances, and excitatory $→$ inhibitory synaptic strength (Fig. 3B). Only variation of the K⁺ Nernst potential within this model explained both the magnitude of the changes seen with older age (2.5% decrease) and following D- $β$HB consumption (5.6% increase) in broadband GFS (Fig. 3C). Sensitivity of brain network instability to the parameter changes was binary: either it was not sensitive to changes in parameters, or network instability changed in non-consistent ways and with large magnitudes that did not reflect the smaller magnitudes seen in human data (Fig. S1). Furthermore, brain network instability was found to be uncorrelated with GFS within the broad (1-40 Hz) frequency band of Leipzig’s LEMON rsEEG dataset (p = 0.090, $r2$ = 0.014, N = 201).

It has been hypothesized that brain aging results from an “energy crisis” in the brain, in which decreased glucose oxidation capability leads to constrained ATP availability for neurons (Jensen et al., 2020; Hoyer, 1982), disrupting neuronal signaling and thus dysregulating the neural circuits that underlie cognitive processing. The fact that metabolism modulated both our EEG- and fMRI-derived biomarkers for brain aging provides robust support for this hypothesis. Our findings, based on a conductance-based neural mass model, suggest that the dysregulation of K⁺ ion gradients, crucially governed by ATP, is the primary driver of the observed changes in brain-wide EEG synchronization.

Dysfunction of Na⁺/K⁺-ATPase causes depolarization of the membrane potential, and thus desynchronization between brain regions (Chesebro et al., 2023). Ion pumps are a sink of metabolic resources in the brain, and thus Na⁺/K⁺-ATP-ase dysregulation is also expected under metabolic constraints, which would further exacerbate ion gradient dysregulation (Baeza-Lehnert et al., 2019). Changes in potassium reversal potential have been implicated in a number of different age-related processes (Scott et al., 1988; Sesti et al., 2010). The depletion of the potassium gradient can stem from damage to potassium channels by reactive oxygen species (Sesti et al., 2010) or from calcium signaling dysregulation (Farajnia et al., 2015; Power et al., 2002), or by a combination of accumulated insults, leading to impairment.

This conductance-based neural mass model provides a biophysically detailed simulation of ion dynamics while allowing for a whole-brain simulation by abstracting away population details. While this derivation preserves local ion dynamics (Larter et al., 1999), it does not fully capture emergent properties of neural populations (e.g., the model preserves cortical wave dynamics (Roberts et al., 2019), but not power spectra). We found that GFS was sensitive to changes in the biophysical model parameters (Fig. 3B), reflecting changes that were on the same scale of magnitude as the changes in GFS seen in the human (non-simulated) data. However, the absolute magnitude of the spectra produced using our modeling approach did not match that of spectra produced from our human rsEEG data. Furthermore, brain network instability was either not sensitive at all to changes in parameters, or for those parameters that did alter network instability, changed it in non-consistent ways and with large magnitudes of change that did not accurately reflect the smaller magnitudes seen in human data (Fig. S1). These limitations may be due to the inability of this neural mass model to capture certain emergent properties of neural populations and have been partially addressed in next-generation neural mass models (Byrne et al., 2020). Brain network instability and neural synchrony were found to be uncorrelated, implying the two metrics quantify distinct neural features, of which the features of the latter are better captured by the conductance-based neural mass model used in this work. There are several potential advances in the field which may lead to better capturing of both neural features. Some examples are the incorporation of networks of firing neurons that model ion and metabolic kinetics to existing models (Dutta et al., 2023), including additional metabolic variables when considering the BOLD signal fluctuations (as is becoming common in systems pharmacology models (Spiros et al., 2014)), and the development of next-generation neural mass models that follow a more “human-like” power spectrum (Byrne et al., 2020). As of yet, however, these models do not fully incorporate the same biophysical detail coupled with scalability as this conductance-based neural mass model does. Fusing these approaches is a topic of current interest in the field of computational neuroscience.

Prior literature has established ion gradient regulation to be the most significant energy sink within the brain (Baeza-Lehnert et al., 2019; Meyer et al., 2022). This is consistent with our results linking emergent whole-brain network effects to ion gradient regulation. Our further isolation of the consistent candidate mechanisms to K⁺ Nernst potentials is also in agreement with in vitro experimental evidence that ketosis directly impacts neuronal K⁺ regulation (Ma et al., 2007). In conclusion, our findings indicate promising avenues for future research, directly examining the association between age-related hypometabolism and disrupted neuronal K⁺ gradients, and their reversal by ketosis. Furthermore, the experimental validation of our findings highlights the role of generative models, guided by clinical neuroimaging, in connecting molecularly targeted therapies to patient outcomes while simultaneously showcasing their utility in mitigating confirmation bias by considering multiple hypotheses simultaneously.

Data and code are available at our website: https://www.lcneuro.org/data/eeg. Leipzig LEMON data are publicly available at http://fcon_1000.projects.nitrc.org/indi/retro/MPI_LEMON.html.

L.R.M.-P., C.W., and K.C. designed research. H.v.N. and C.P. performed research and collected data. H.v.N. analyzed data. H.v.N., A.G.C., and H.H.S. performed computational modeling. K.C. contributed reagents. H.v.N., L.R.M.-P., and A.G.C. wrote the paper. H.v.N., A.G.C., K.C., H.H.S., C.W., and L.R.M.-P. edited the paper.

The research described in this paper was funded by the White House Brain Research Through Advancing Innovative Technologies (BRAIN) Initiative (grants NSFNCS-FR 1926781 to L.R.M.-P.) and the Baszucki Brain Research Fund, United States (L.R.M.-P.). CW acknowledges support from the Marie-Josée Kravis Fellowship in Quantitative Biology. A.G.C. was also supported by the NIHGM MSTP Training Award, United States T32-GM008444.

The intellectual property covering the uses of ketone bodies and ketone esters is owned by BTG Plc., Oxford University Innovation Ltd., and the NIH. K.C., as inventor, will receive a share of the royalties under the terms prescribed by each institution. K.C. is a director of TΔS Ltd., a company spun out of the University of Oxford to develop products based on the science of ketone bodies in human nutrition.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00163.