The goal of our experimental study is to investigate the close link between neuronal brain activity and subsequent metabolic responses of the human organism at a systemic level. In 25 , the methods and results of the experimental study are described in detail. The study design is depicted in Figure 1.

In a randomized sham-controlled crossover design, a homogeneous group of 15 healthy young male volunteers with a body mass index of 23.2 ± 0.38 kg/m² is examined. Neuronal brain activity is stimulated by transcranial direct current stimulation (tDCS) during the time interval t ∈ [125, 145] (minutes). Transcranial stimulation of the brain causes transient effects on motor cortical excitability outlasting the stimulation period 26 . Sham stimulation serves as control condition. For sham stimulation, electrodes are placed at the same site without current stimulation. TDCS-induced effects on cerebral energy metabolism and systemic glucose regulation are measured. The study is carried out in accordance with the Declaration of Helsinki (2000) of the World Medical Association and has been approved by the ethics committee of the University of Lübeck. Each participant gave written informed consent.

³¹Phosphorus magnetic resonance spectroscopy (³¹P-MRS) allows performing non-invasive in vivo measurements of brain metabolites that are centrally involved in the energy metabolism. Phosphate metabolites such as adenosinetriphosphate (ATP), i.e., the sum of α-, β-, and γ-ATP, as well as phosphocreatine are measured in the cortex reflecting the overall high-energy phosphate turnover 27 . Here, the ratio of ATP and inorganic phosphate (Pi) is evaluated as an indicator of the intracellular energy status 25 28 . ³¹P-MR spectra are measured at times t = 115, 160, 205, 215, 225, 235, 245, 385 (minutes).

During euglycemic-hyperinsulinemic clamping, an insulin infusion at the predetermined fixed dosage of 1.5 mU (kg min)⁻¹ and a variable glucose infusion are administered in order to reach stable plasma glucose concentrations between 4.5 and 5.5 mmol/l. Under steady-state conditions of euglycemia, the glucose infusion rate equals glucose uptake by all tissues in the body 29 and is therefore a measure of glucose tolerance. To monitor the peripheral glucose metabolism, blood samples of glucose and insulin are regularly taken at times t = 0, 110, 155, 210, 220, 230, 240, 300, 340, 380 (minutes).

While blood glucose and insulin concentrations do not differ between conditions, overall cerebral high-energy phosphate measurements display a biphasic course upon tDCS as compared with the sham condition, see 25 and Figure 2. An initial energetic drop in the ATP/Pi ratio upon tDCS is observed. Subsequent ³¹P-MR spectra reveal a rapid increase above the control condition followed by a return of the ATP/Pi ratio to baseline levels. Glucose infusion rates show the same biphasic response to tDCS indicating that transcranial stimulation improves systemic glucose tolerance 29 30 31 . Measurements of the hypothalamus-pituitary-adrenal (HPA) hormonal system reveal decreasing concentrations of circulating stress hormones such as cortisol upon tDCS (compare 25 ).

The experimental data demonstrate that transcranial brain stimulation not only evokes alterations in local neuronal processes but also clearly influences brain energy homeostasis and peripheral metabolic systems regulated by the brain 25 . Hence, manipulating brain activity by tDCS affects peripheral metabolic regulation such as the glucose metabolism and related neuroendocrine mediators. Effects of tDCS on cerebral ATP/Pi, blood glucose, and insulin as compared with the sham stimulation are shown in Figure 2.

Nevertheless, the mechanisms underlying these experimental observations remain unknown. Concerning the specific mechanisms by which neuronal excitation, via a drop in high-energy phosphate content, improves glucose tolerance, one can only speculate at this point.

The objective of the following mathematical analysis is to gain insight into physiological mechanisms underlying the effects of brain stimulation on cerebral and peripheral energy metabolism. In order to clarify the underlying mechanisms in this context, we combine experimental data with the mathematical model introduced in the following subsection.

Physiological processes may be described by systems of ordinary differential equations,

dy dt t , p = f t , y , p with y 0 = y 0 ,

where y : R × R m → R n and f : R × R n × R m → R n are time-depending functions with t∈[0,T]. Here, y ₀ denotes the initial conditions. We collect unknown model parameters to be estimated in the vector p = (p ₁, …, p _m )^⊤. In the following, given parameters regarded as constant are assembled in a vector c.

A collection of mathematical models of this kind describing interactions of main components of the human glucose metabolism can found in the book of Chee and Fernando 17 . Most of these models are based on the glucostatic or lipostatic theory. Exemplarily, one could mention the well-known Minimal Model 32 , the Ackerman Model 14 , or more recently models published in 19 21 . Additionally, mathematical models developed in the context of the Selfish Brain Theory can be found in the review paper by Chung and Göbel 33 .

Here, we regard a mathematical model of the human whole body energy metabolism considering the brain not only as energy consumer but more importantly as a superordinate controller, compare 12 . The model includes energy metabolites in separated compartments, energy fluxes between these compartments, and signals directing energy fluxes within the organism, see Figure 3.

The brain-centered model of the energy metabolism is given by the system of four ordinary differential equations,

dA dt = p 1 G A − p 4 , dG dt = − p 1 G A − c 1 GI + p 2 R G + G ext , dI dt = p 3 A − c 2 I + I ext , dR dt = c 1 GI − p 2 R G − p 5 .

Cerebral energy content, i.e., high-energy phosphates in the brain, is denoted by A. Furthermore, G is the blood glucose concentration, and R specifies energy resources in the body, which comprise available energy reserves foremost in liver, muscle, and fat tissue. Note, various types of energy, such as fat, glycogen, and glucose, are combined in the energy resources compartment. In addition to energy metabolites, the model contains the control signal, identified as blood insulin concentration.

Conceptually, our model bases on conservation of energy. In general, stimulatory influences are modeled as proportional relations and inversely proportional relations describe prohibitive influences. The glucose flux from the blood into the brain crossing the blood brain barrier is proportional to intensified by a factor p ₁ [M s⁻¹]. This factor quantifies the glucose flux rate from the blood into the brain across the blood brain barrier. Suppose the cerebral energy content A is low, the facilitated diffusion process from the blood compartment into the brain via the passive glucose transporter GLUT1 is accelerated. Contrary, high levels of A inhibit this flux. This can be seen as “energy on demand” of the brain. This concept has been published in 12 . Glucose G needs to be available in the blood compartment to reach the brain.

Our model combines energy resources and metabolites, such as glycogen, glucose and lactate, in the compartment G. The energy flux from the resources into the blood compartment, which is composed of several sub-mechanisms, is proportional to the energy resources R and anti-proportional to the actual blood glucose concentration G with a flux rate with a factor p ₂ [M s⁻¹], i.e., p ₂ R / G. This flux includes endogenous glucose production by the liver amongst others.

The hormone insulin acts not only as local response to the blood glucose concentration. Moreover, it is regarded as central feedback signal of the brain with an insulin secretion factor p ₃ [s⁻¹]. Notice that with low cerebral energy, ventromedial hypothalamic centers inhibit pancreatic insulin secretion 34 35 36 . Being an anorexigenic hormone, peripherally secreted insulin is a key feedback signal to the brain reducing food intake and systemic glucose uptake 13 15 . Energy consumption of the brain and energy consumption by the periphery are denoted by p ₄ [M s⁻¹] and by p ₅ [M s⁻¹].

Insulin-dependent glucose uptake from the blood into the energy resources compartment is modeled by c ₁ GI with a parameter c ₁ [(M s)⁻¹]. This flux mainly comprises glucose uptake into the peripheral stores, i.e., muscle and fat tissue. To accelerate this flux, glucose and insulin need to be available in the blood at the same time in order to activate glucose uptake via the insulin-dependent glucose transporter GLUT4.

Degradation of insulin is supposed to be of first order with the insulin degradation rate c ₂ [s⁻¹]. External glucose infusion is denoted by G _ext (t) [M s⁻¹], insulin infusion is I _ext (t) [M s⁻¹].

To meet the simplified notation from Equation (1) we collect the state variables y = (A, G, I, R)^⊤ and the parameters p = (p ₁, …, p ₅)^⊤. The constants c = (c ₁, c ₂)^⊤ and time-dependent external infusions G _ext , I _ext are not explicitly shown in (1). Notice that all parameters, constants, and states are non-negative. Model properties he been investigated in detail and it has been shown that the model realistically reproduces qualitative and quantitative behavior of the whole body energy metabolism even for a large class of physiological interventions (see 11 12 for details).

To accommodate the characteristics of the experimental study, the dynamical system (2) slightly differs from the model introduced in 12 . First, glucose and insulin infusions are administered. Secondly, ingestion of food is neglected since no food intake occurs during the examinations.

In the following, we introduce the general technique to estimate the model parameters p. The solution y of the dynamical system (2) varies with respect to the model parameters p. In order to validate and implement model predictions, the mathematical model needs to be compared to experimental data. Here, we analyze the model behavior in identifying the unknown model parameters p = (p ₁, …, p ₅)^⊤ of our dynamical system in Brain-centered energy metabolism model section using experimental data presented in Experimental study section.

In general, parameter identification problems for ordinary differential equations can be stated as follows

min p , y 0 1 2 ∥ σ − 2 y τ , p − d ∥ 2 2 + S λ , p subject to dy dt = f t , y , p and y τ 1 = y 0 .

Equation (3) states a classical constrained optimization problem, where constraints are given by an initial value problem 37 . With ∥·∥₂ denoting the Euclidian norm, we minimize the distance between the model solution at the times τ = (τ ₁, …, τ _k )^⊤, where the data are measured and given data d = d 1 , … , d k ⊤ ∈ R k ; the minimization is constrained by the validity of the mathematical model. With σ 2 ∈ R denote the variance of the data assumed to be independent and identically distributed. The term S(λ,p) is a convex regularization inducing prior knowledge on the parameters for which we will provide details later. For convenience, we assume the model function f to be continuously differentiable. Note that the optimization problem (3) corresponds to the maximum likelihood estimator including a prior and is a standard formulation of a parameter identification problem 38 . Note that estimator (3) infers that the errors are independently and identically normal distributed.

We seek to find p ^ ∈ R m minimizing Equation (3). The optimization problem (3) can only be targeted by numerical optimization methods. Problem (3) is a typical inverse problem since data and model are given and we aim to identify the model parameters. This inverse problem is well known to experimentalists and various methods have been established to solve such type of parameter estimation problems. Most commonly used are single and multiple shooting methods. Both methods face certain advantages and disadvantages. Single shooting methods are easy to implement but are not robust to initial guesses of the model parameters and optimization algorithms are likely to fail or to find “non-optimal” local minimizer. Multiple shooting methods are by far more robust and are shown to face better convergence. However, multiple shooting methods solve constrained optimization problems, which dramatically increase the algorithmic complexity 37 39 .

Here, we follow an approach similar to methods proposed by Ramsay et al. 40 , Chung and Haber 41 , or Poyton et al. 39 . Equation (3) can be restated as

min p , y 0 1 2 ∥ σ − 2 y τ , p − d ∥ 2 2 + S λ , p subject to dy dt − f t , y , p 2 = 0 and y τ 1 = y 0 .

The equivalence of optimization problem (3) and (4) stays true for any appropriate integral norm ∙ (here we choose the L ² –norm). Constrained optimization problems such as (4) are commonly approximated by performing a Lagrangian relaxation 38 . We get:

min p , y 0 1 2 ∥ σ − 2 y τ , p − d ∥ 2 2 + α 2 dy dt − f t , y , p 2 + S λ , p subject to y τ 1 = y 0 ,

with the Lagrangian multiplier a ≥ 0. Notice that for increasing a unconstrained optimization problem (5) becomes optimization problem (4). Next, we choose a standard discretize-then-optimize approach to solve problem (5), numerically. We let u = u 1 , … , u k ⊤ ∈ R k be a approximation of y at the data points τ ₁, …, τ _k and D _t be a finite differences operator approximating dy / dt at τ ₁, …, τ _k , then we can restate the optimization problem (5) as the discretized and unconstrained optimization problem.

min p , u Φ p , u = 1 2 ∥ σ − 2 u − d ∥ 2 2 + α 2 ∥ D t u − f τ , u , p ∥ 2 2 + S λ , p .

Notice that we can neglect the remaining constraint in (5) since y ₀ = u ₁ is already included in the search parameters and is therefore always fulfilled. Equation (6) describes the general parameter estimation framework. This optimization problem has the advantages that it is robust and the unconstrained nature allows to use fast gradient-based methods, for details see 40 41 . Notice that if the data points τ ₁, …, τ _k are t dense, one may want to utilize a spline function s with knots τ and coefficients q at dense points ξ = ξ 1 , … , ξ l ⊤ to capture the dynamic of the differential equation. Then Equation (6) reads

min p , q Φ p , q = 1 2 ∥ σ − 2 s τ , q − d ∥ 2 2 + α 2 ∥ s ' ξ , q − f ξ , s ξ , q , p ∥ 2 2 + S λ , p .

One further choice to make is choosing the regularization term S(λ,p). A most common choice for the regularization term is

S λ , p = λ 2 ∑ p i − p ¯ i p ¯ i 2 ,

where p ¯ ∈ R m is a given parameter representing physiological parameter values and λ ≥ 0. With this we have established a parameter estimation method to tackle the model of Brain-centered energy metabolism model section.

Next, we present the parameter estimation setup for the model in Section 2.2 with the given data from Section 2.1. As derived in Section 2.2, all parameters of our model (2) have a physiological interpretation. Since insulin is at hyperphysiological levels in our experimental examination we will consider the dependent parameters c ₁ and c ₂ to be constant. Flakoll PJ, Wentzel LS, Rice DE, Hill JO and Abumrad NN 42 quantify the insulin-dependent whole body glucose uptake corresponding to c ₁ ≈ 0.06 (pM min)⁻¹. Information about insulin clearance is given by 43 resulting in c ₂ ≈ 1.4 min⁻¹. This yields the fixed parameter values c = (c _1, c ₂) = (0.06, 1.4)^⊤.

For the parameter p ¯ = p ¯ 1 , … , p ¯ 5 ⊤ , we choose physiologically relevant parameter values gathered from the literature. We use p ¯ 1 ≈ 0.15 mM/min for glucose transport rate across the blood brain barrier 44 . Baron AD and Clark MG 45 specify the glucose flux between peripheral stores and blood with p ¯ 2 ≈ 0.6 mM/min. For the insulin secretion rate, we choose p ¯ 3 ≈ 20 min⁻¹, see 46 . The work by Flakoll PJ, Wentzel LS, Rice DE, Hill JO and Abumrad NN 42 provides the maximal rate of glucose utilization. This grants an insight into the peripheral energy consumption p ¯ 5 ≈ 0.7 mM/min and cerebral energy consumption p ¯ 4 ≈ 0.175 mM/min since the brain uses up to 20% of total body glucose 10 . Hence, we choose the regularization value p ¯ = 0.15 , 0.6 , 20 , 0.175 , 0.7 ⊤ .

In order to establish global convergence, we choose a Monte Carlo sampling technique of the initial guess p ₀. We pick 500 randomly chosen normally distributed samples with mean p ¯ . For each sample we calculate the minimizer and choose the overall minimizer p ¯ to be the minimizer of the objective function Φ, see Equation (4).

To solve the optimization problem (4) numerically we use a Gauss-Newton method with Armijo line search, see 38 . The regularization parameters a and λ are set to a = 5 10⁻³ and λ = 10⁻⁷ for our numerical investigations. By empirical observations these values lead to a good balance between under- and overfitting.

The goal of this section is to show that parameter values of our dynamical system (2) can be chosen such that the resulting mathematical model is able to approximate the experimental measurements of the study described in Section 2.1. To validate the mathematical model, we first identify the model parameters by mean values of the experimental data for cerebral ATP/Pi, blood glucose, and insulin concentration in the sham stimulation condition. Note, for the energy resources compartment no experimental data are available, and we arbitrarily set R = 50 [mmol/l]. We let the corresponding variance σ² = 1000 be fairly high, which allows the optimization method to pick adequately different values for R without contributing much to the function Φ. Thereby, we obtain the minimizer p ^ = 0.0049 , 0.0683 , 11.3768 , 0.0025 , 0.8331 ⊤ with objective function value Φ = 6.7611 10⁻⁵ for the placebo group.

In order to test the accuracy of the estimation, we calculate the distance between data and calculated model by the relative mean squared error

θ p = 1 K ∑ k = 1 K y k − d k d k 2 ,

where K is the number of data points. The vector d contains the experimental data of cerebral high-energy phosphates, blood glucose, and insulin concentrations, (experimental data and model prediction are shown in Figure 4). The vector y specifies the respective solution of the dynamical system (2) evaluated at the data time points in the interval [110, 385] min. We obtain the relative mean squared error θ p ^ = 7.6012 · 10⁻³.

The parameter estimates for the mean data from the tDCS intervention read p ^ = 0.0092 , 0.0474 , 10.8068 , 0.0006 , 0.8507 ⊤ , with Φ = 2.0410 10⁻⁴ and θ p ^ = 9.6891 · 10⁻³ and are illustrated in Figure 5.

Since we conduct parameter identification with 500 random initial values p ₀ as well as the values Φ and θ are small, it is likely that the estimated values p ^ are the global minimum. Therefore, the model parameters can be chosen such that the mathematical curves closely approximate the data. Note that the estimated parameter values p ^ stay not too far from the initial guesses p ¯ , taken from the literature 4 17 39 40 . However, uncertainties in the data may amplify uncertainties in the estimated parameters.

With the estimated parameter values p ^ we solve the forward problem. Figures 4 and 5 show the model curves y = (A, G, I, R)^⊤ compared with experimental data for sham and tDCS interventions. The model predictions closely approximate the measurements. Ideally a glycemic-hyperinsulinemic clamp generates constant blood glucose and plasma insulin levels. Constant glucose and insulin levels by themselves do not contribute much insight to the parameter estimation. However, the glucose and insulin fluxes G _ext and I _ext are recorded and are included in the parameter estimation. Note that blood glucose predictions are slightly unsteady due to the external glucose infusion.

To validate the mathematical model, we perform a sensitivity analysis. The results of the sensitivity analysis for the two central parameters p ₁ and p ₄ are illustrated in Figures 6 and 7, respectively. We find that model predictions of cerebral ATP A, blood glucose G, and insulin I are sensitive to changes in these parameters. Modifying model parameters p ₂, p _3, p ₅ and initial values y ₀ also affects the predicted profiles, indicating their sensitivity to changes in these parameters as well (analysis not shown). Energy resources R are sensitively affected mainly by variations in peripheral energy consumption p ₅.

While Section 3.2 considers the population of placebo and tDCS data, we next investigate each individual data set using the same settings as before. We identify the parameters p ^ 1 , … , p ^ 5 for each volunteer participating in the experimental study separately. Thereby, we are able to identify reasonable parameter values for each subject. Box-and-whisker diagrams of the estimated parameter values are shown in Figure 8.

We statistically investigate the effects of tDCS on estimated parameter values. Each model parameter has a physiological interpretation. Thereby, our examinations provide insight into underlying physiological mechanisms that are not measurable in the experimental study.

We identify outliers p ^ _ i below quartile Q ₁ and above quartile Q ₃ by p ^ i < Q 1 − 1.5 Q 3 − Q 1 and p ^ i > Q 1 + 1.5 Q 3 − Q 1 , respectively (see red plus in Figure 8). Individuals with extreme outliers are plotted individually. Mean values and standard deviations of the estimated parameters p ^ are listed in Table 1.

In order to investigate the effects of tDCS on underlying physiological mechanisms in the systemic energy metabolism, we compare mean estimated parameter values for tDCS and sham stimulation. Identified parameter ratios between tDCS and sham condition are shown and statistically significant values are highlighted in Table 2. Each ratio reflects the effect of tDCS on the respective physiological process described by the model parameter. To test if mean estimated parameter values upon tDCS and sham stimulation significantly differ, we conduct an analysis of variance, a two sample F-test. Hereby, a P-value < 0.1 was considered significant.

^aSignificance, P = 0.04, ^bSignificance, P = 0.08, and ^cSignificance, P = 0.03. Statistical analyses were performed by analysis of variance to identify main effects in experimental condition comparisons.

Glucose flux p ^ 1 across the blood brain barrier as well as cerebral energy consumption p ^ 4 significantly increase upon tDCS compared with sham stimulation. In contrast, the energy flux p ^ 2 from peripheral tissue into the blood significantly decreases in the stimulation condition. Insulin production rate p ^ 3 slightly decreases, and peripheral energy consumption p ^ 5 does not vary significantly between different experimental conditions.

The physiological interpretation of the obtained ratios reads as follows: TDCS excites cortical neurons. Our results verify that this leads to significantly increased cerebral energy consumption p ^ 4 . Accordingly, our results reflect an increase in cerebral glucose uptake p ^ 1 from the blood into the brain crossing the blood brain barrier in response to excitation-induced energy consumption. This mechanism assures an adequate energy supply of the brain.

Moreover, we analyze cerebral energy supply on demand. Therefore, we calculate the ratio of energy flux p ^ 1 across the blood brain barrier to cerebral energy consumption p ^ 4 . This yields p ^ 1 / p ^ 4 ≈ 14.37 for tDCS compared to p ^ 1 / p ^ 4 ≈ 1.95 for sham stimulation and the ratios significantly differ (P-value < 0.01). Our result reflects that the allocation mechanism of the brain is distinctly enhanced by transcranial stimulation. This reveals beneficial effects of neuronal activation on cerebral energy homeostasis and glucose metabolism.

Insulin production rate p ^ 3 is reduced upon tDCS. This reduction may be seen as enhanced allocation mechanism as well. With low blood insulin concentrations less energy is transported into the peripheral stores via the insulin-dependent glucose transporter GLUT4. Hence, available glucose can be transported into the brain via the insulin-independent transporter GLUT1.

Our results are in line with the hypothesis that the brain supplies itself with glucose in dependence of its own needs 3 9 . The allocation mechanism of the brain is activated with low high-energy phosphate content in the brain and the available energy is allocated to the brain.

Our findings explain the experimentally observed biphasic effect of tDCS on cerebral energy content, compare 25 . Neuronal excitation increases cerebral energy consumption p ^ 4 leading to a drop in high-energy phosphate content in the brain. Our results reveal a significant increase in glucose transport rate p ^ 1 across the blood brain barrier upon tDCS. Furthermore, cerebral energy supply on demand p ^ 1 / p ^ 4 is enhanced. These observations explain the subsequent rise of cerebral high-energy phosphates above basal levels. The increase in cerebral glucose uptake in response to excitation-induced energy consumption thereby provides a conceivable explanation for the experimental observations. The allocation mechanism of the brain ensures cerebral energy supply.

Additionally, we find a significant decrease in energy flux p ^ 2 from peripheral tissue into the blood. This result is in line with the experimental measurements revealing suppression of the HPA system. Stress hormones such as cortisol induce gluconeogenesis mainly in the liver and thereby enhance the glucose flux from the resources into the blood. Reduced peripheral glucose release in conjunction with increased cerebral glucose uptake provides explanations for improved glucose tolerance.

Peripheral energy consumption p ^ 5 remains almost identical since the subjects are not physically active during examinations and therefore, the result meets the expectations. The obtained results are physiologically reasonable. Our investigations provide evidence for the experimentally observed changes in cerebral and peripheral energy metabolism upon tDCS. The results shed light on underlying physiological processes that are not measurable within the scope of a human experimental approach.