Multi-timescale Modeling of Activity-Dependent Metabolic Coupling in the Neuron-Glia-Vasculature Ensemble

Glucose is the main energy substrate in the adult brain under normal conditions. Accumulating evidence, however, indicates that lactate produced in astrocytes (a type of glial cell) can also fuel neuronal activity. The quantitative aspects of this so-called astrocyte-neuron lactate shuttle (ANLS) are still debated. To address this question, we developed a detailed biophysical model of the brain’s metabolic interactions. Our model integrates three modeling approaches, the Buxton-Wang model of vascular dynamics, the Hodgkin-Huxley formulation of neuronal membrane excitability and a biophysical model of metabolic pathways. This approach provides a template for large-scale simulations of the neuron-glia-vasculature (NGV) ensemble, and for the first time integrates the respective timescales at which energy metabolism and neuronal excitability occur. The model is constrained by relative neuronal and astrocytic oxygen and glucose utilization, by the concentration of metabolites at rest and by the temporal dynamics of NADH upon activation. These constraints produced four observations. First, a transfer of lactate from astrocytes to neurons emerged in response to activity. Second, constrained by activity-dependent NADH transients, neuronal oxidative metabolism increased first upon activation with a subsequent delayed astrocytic glycolysis increase. Third, the model correctly predicted the dynamics of extracellular lactate and oxygen as observed in vivo in rats. Fourth, the model correctly predicted the temporal dynamics of tissue lactate, of tissue glucose and oxygen consumption, and of the BOLD signal as reported in human studies. These findings not only support the ANLS hypothesis but also provide a quantitative mathematical description of the metabolic activation in neurons and glial cells, as well as of the macroscopic measurements obtained during brain imaging.


Introduction
The mammalian brain exhibits remarkable processing power. It is at the same time energy efficient. The design features that allow such efficient computation are mapped in cellular and molecular components and their roles in information processing. Concurrently, these features are anchored in, and constrained by, the universal metabolic chains that provide energy to cells. Deciphering the metabolic code and the neural code are thus tandem requirements for a comprehensive understanding of brain function. Understanding the metabolic underpinnings of information processing is also of added value to understanding the etiology and progression of neuropsychiatric and neurodegenerative disorders [1,2]. The picture that emerges from this dynamical system will reflect the cooperative function of neurons, glia and the vascular system.
Glutamate, the brain's major neurotransmitter, effects numerous cascades and processes in brain cells [3,4]. Among them, astrocytes couple synaptic activity to energy metabolism via a sodium-dependent uptake of glutamate [5]. The ensuing cascade of molecular events leads to the glycolytic processing of glucose and the release of lactate by astrocytes. A comprehensive model of brain energy metabolism must consider oxidative and non-oxidative glucose consumption, intracellular and extracellular compartmentalization and transport of choke-point metabolic intermediates such as lactate and pyruvate, as well as feedback mechanisms that report local synaptic and intrinsic neuronal activity [6,7]. These pathways are in turn complicit in the molecular and cellular mechanisms that contribute to the still poorly understood readout of functional brain imaging [8].
The role of astrocytes and how they metabolically interact with neurons is well supported experimentally; some mostly theoretical considerations, however, have challenged this view. Magistretti and colleagues proposed that clearance of glutamate from the synaptic cleft by astrocytes could be coupled to glycolysis and subsequent lactate production [5]. Lactate produced in this way would then be transported to the extracellular space. Controversy remains surrounding directionality and timing of lactate flow in the brain; while a neuron-to-astrocyte lactate system (NALS) is proposed by some [9,10], an astrocyte-to-neuron direction (ANLS) is supported by a large set of experimental evidence [11]. Biophysical models also weigh-in on the conditions and sequences of events required for lactate production and consumption [12][13][14][15].
The existence of an extracellular pool of lactate likely used as an energy reservoir at the onset of stimulation has been observed in rats and humans [16,17]. The distribution of monocarboxylate transporters at the membrane of neurons and astrocytes supports the hypothesis of a net transfer of lactate from astrocytes to neurons through the extracellular space [18]. Glial cells have been observed to take up most glucose [19,20], while neurons are responsible for the largest part of brain oxygen consumption [21,22]. Additional evidence comes from the direct measure of NADH transients in brain slices, showing that neurons display early oxidative metabolism following presynaptic activity, while astrocytes display a delayed activation of glycolysis but no detectable oxidative response [23]. Nevertheless, the ANLS-hypothesis is still debated and challenged with arguments focusing now on the exact interpretation of the above observations [24,25].
Nicotine adenine dinucleotide, either oxidized or reduced (NAD + , or NADH), is a workhorse cofactor that acts as a central electron broker for metabolic redox cycles including glycolysis, the citric acid cycle (Krebs, TCA) and oxidative phosphorylation. Owing to its high UV wavelength absorption, it is also responsible for cellular auto-florescence. This coincidence makes it a useful indicator of metabolic activity. NADH is an important metabolic signal because it is produced or used during both mitochondrial activity and activation of the glycolytic pathway, and because it cannot diffuse freely through the mitochondrial membrane but needs to be transported by appropriate shuttles. Fluctuations of the NADH concentration measured in the appropriate cellular compartments can then indicate increased or decreased oxidative and glycolytic metabolism. A critical previous finding in this regard was the observation of early and late activity-dependent phases of metabolic activity with the early phase taking the form of a NADH "dip" and the late phase appearing as a NADH "overshoot" with a longer time constant of decay [23]. Interestingly, these phases also correlate with the fluctuations of the extracellular lactate concentration as determined in animals [16] and humans [26,27].
The emerging consensus is that the early phase represents NADH depletion in the dendrites of active neurons and that the overshoot represents glycolytic activity that results in the accumulation of NADH. This activity results in the high production of lactate in astrocytes as rapid glycolysis overtakes the subsequent consumption by oxidative pathways [23,26,27]. The accumulation and transportation of lactate between glial cells and neurons may in turn serve as an activity-dependent buffer that is informed by the neuronal release and glial uptake of glutamate [5]. It might also act as a signaling molecule to the vasculature [28] or to brain cells via binding to the G-protein coupled receptor GPR81 [29].
The preference for lactate over glucose as an energy substrate in neurons has been demonstrated in vivo as well as in vitro [30,31], as has a neuroprotective role for lactate in the case of insulin-dependent hypoglycemia [32] and other conditions [33,34]. The role of the ANLS in homeostatic maintenance involves the regulation of blood glucose [35] and sodium [36]. While the ANLS hypothesis, since its initial formulation [5], does not preclude the use of glucose by neurons as an energy substrate, it has been challenged by some studies defending the view that glucose, rather than lactate, is the sole energy substrate for oxidative metabolism in neurons [37,38].
Previous modeling efforts have advanced our knowledge of this functional metabolic network by demonstrating that lactate consumption by neurons occurs early in the stimulus regimen and that the early and late lactate transients correspond to the activity of two distinct populations of cells, neurons and glia. The current study builds on and complements those and other models [9,10,[12][13][14][15][39][40][41][42] and addresses unresolved mechanisms of neuron-glial metabolic and vascular coupling. The model is based on several previous studies [13,40] with five significant improvements: 1) the compartmentalization of NADH between cytosolic and mitochondrial compartments; 2) the linking of metabolic and Hodgkin-Huxley formalisms; 3) the input to the neuronal and astrocytic compartments formulated as a presynaptic glutamatergic stimulation; 4) the model explicitly and continuously updates reversal potentials; and 5) the model was constrained using in vitro data and correctly predicts in vivo results without the need of invoking glycogen (which is deliberately excluded from the model).
In this paper, we will show that a biophysical model of astrocyte-neuron metabolic interactions designed following these principles leads to the presence of an activity-dependent lactate shuttle from astrocytes to neurons and that this model can reproduce the evoked response of NADH in its various compartments as reported by Kasischke and colleagues [23]. We will subsequently show that our biophysical model correctly predicts qualitatively-and to some extent quantitatively-the evoked responses of tissue lactate and tissue oxygen as observed in the rat brain in vivo [16]. Finally, our biophysical model predicts the evoked responses of tissue lactate, of the BOLD signal and the glucose and oxygen consumption as observed in the human brain in vivo.

Methods
This in silico model represents a dynamic and integrative analysis of compartmentalized metabolism and its relation to neuronal signaling in the central nervous system. The model was designed based on knowledge of the underlying biophysics and required input parameters and equations from multiple species, time and spatial scales. The model is inspired by previous work from Aubert and colleagues [13,40] and consists of four compartments: neuron, astrocyte, capillary and extracellular space (see Fig. 1). These compartments are referred to by the subscripts n, g (for glia/astrocyte), c and e respectively. In addition, the neuronal and astrocytic compartments are further divided between cytosolic and mitochondrial sub-compartments to account for the compartmentalization of nicotinamide adenine dinucleotide (NADH). These are referred to by the superscripts cyto and mito. Transport between compartments is noted with the subscripts of both compartments; for instance, transport from the neuronal compartment to the extracellular space is labeled with ne or, conversely, en.
The model is formulated as a series of 33 differential equations adapted from previous work [13,40] with the following improvements: the compartmentalization of (mostly) NADH between the cytosolic and mitochondrial compartments; the model was joined to a Hodgkin-Huxley-type model [43,44]; the input to the neuronal and astrocytic compartments is formulated as a glutamatergic input and not anymore as an abstract stimulation; the model explicitly models sodium entry and extrusion in both the neuronal and astrocytic compartments, continuously updating the corresponding reversal potentials. The model, which for the first time bridges mathematical descriptions of energy metabolism and Hodgkin-Huxley equations, was constrained on in vitro data and correctly predicts in vivo results.
Parameters that were difficult to determine experimentally such as transport constants were left free to vary [10,42]. Free parameters were then optimized so that the model reproduces the experimental results presented in [23]. The number of free parameters was maintained as small as possible by enforcing constraints on the value of metabolites at steady state. Simulations with randomized fluctuations in parameter values (up to plus or minus 10% of reported values) did not reveal significant changes in behavior of the model. This approach successfully predicts qualitatively and quantitatively in vivo measurements in rodents and in humans (see Results, Fig. 5 and 6).
The pooling together of equations from various sources is unfortunately necessary for construction of such a broad and multi-dimensional model. There is no single source of equations that can be tapped for this model. Note however that most equations come from only two published sources. The variables, their steady-state values and the corresponding governing equations are given in Table 1 [40] and were originally introduced in ref. [45] for the equations describing metabolism and by Buxton and colleagues [46] for the equations describing the vascular dynamics. Equations (A.9)-(A.12), (A.14) and (A.15) are original. Equations (A.23)-(A.26) describe the neuronal membrane excitability following the Hodgkin-Huxley formalism and the neuronal calcium dynamics. They are adapted from [44].
All the fluxes and currents appearing in Table 1, as well the equations describing the dynamics of the gating variables, are given in Table 2. Like for Table 1, these equations are taken from [40,[44][45][46] except equations (A.35)-(A.37) which are original. All rates and state variables are given per unit cell volume (neuron or astrocyte) or per unit capillary volume to the exception Figure 1. Model structure. The model is divided in four main compartments: a neuronal compartment, an astrocytic compartment, the extracellular space and a vascular compartment. Neurons and astrocytes are further divided between a cytosolic sub-compartment and a mitochondrial sub-compartment to account for the compartmentalization of oxidative and glycolytic metabolisms. Both neurons and astrocytes contain a metabolic network including glycolytic enzymes, lactate dehydrogenase, glucose and lactate transporters, NADH shuttles, oxidative metabolism, phosphocreatine and the Na, K-ATPase electrogenic pump. Additionally, the neuronal compartment contains voltage-and calcium-gated ion channels following the Hodgkin-Huxley formalism. The system as a whole is driven by two independent inputs (red). First, a glutamatergic presynaptic population activates AMPA receptors on the neuronal membrane and excitatory amino-acid transporters (EAATs) on the astrocytic membrane. The activation of both AMPA receptors and EAATs leads to an increase of the intracellular sodium concentration which activates the energy consuming Na, K-ATPase pump and subsequent metabolic processes. Activation of AMPA receptors also depolarizes neurons and might lead to the generation of action potentials, which will also lead to an increase in intracellular sodium in the neuronal compartment via opening of voltage-gated sodium channels. Second, the cerebral blood flow (CBF) is modulated as a separate input. For comparison to in vitro experiments using acute brain slices, the only input is the presynaptic neuronal population while the supply of oxygen and glucose that normally comes from the CBF is held constant as a proxy for the laminar flow of a controlled perfusing solution. Finally, the extracellular space is the place where cells and capillaries exchange metabolites Additional equations ADP x is given as a function of the ATP concentration (x stands for n or g). It reads: with A = AMP x +ADP x +ATP x = 2.212 mM the total adenine nucleotide concentration and  q AK = 0.92 the adenylate kinase equilibrium constant [40,45]. As a consequence: with u x = q 2 AK +4Áq AK Á(A/ATP x -1).
Oxygen concentration at the end of the capillary Calcium-dependent potassium current Flow out of the venous balloon Input to the model The model receives input from a presynaptic excitatory population. Glutamate released by excitatory presynaptic neurons drives the intracellular sodium concentration in neurons and astrocytes and activates AMPA receptors on neurons, thus inducing a synaptic current I syn . The presynaptic population contains N exc excitatory neurons discharging at frequency f exc (t). This presynaptic population thus generates an excitatory conductance g exc (t) given by: with g = 7.8Á10 -6 mSÁcm -2 Ásec the total surface under the conductance evoked by one excitatory event [47,48]. The corresponding synaptic current is then given by: with ψ n the neuronal membrane voltage and E AMPA = 0 mV the reversal potential of AMPA ionotropic receptors. It is estimated that about two thirds of the current generated at AMPA receptors is due to a flow of sodium ions [49]. Sodium also flows through voltage-dependent sodium channels when the neuron is active (I Na ). As a consequence, the sodium drive to the neuron is approximated by: with S m ÁV n = 2.5Á10 4 cm -1 the ratio between neuronal membrane surface and neuronal volume and F = 9.64853Á10 4 CÁmol -1 the Faraday constant. Finally, the glutamate is cleared from the synaptic cleft by excitatory amino acid transporters located on the astrocyte membrane. Those transporters use the electrochemical sodium gradient to transport glutamate with a stoichiometry of three sodium ions for one glutamate molecule. We thus write the sodium drive to the astrocyte as follows: with Δ glut = 2.25Á10 -5 mM a constant, which corresponds to the total amount of glutamate released in the synaptic cleft by each presynaptic action potential multiplied by the ratio between synaptic and astrocytic fractional volumes. For the sake of simplicity, we assume that f exc (t) always follows the same temporal dynamics exponentially decaying from f 0 = 3.2 Hz to f 1 = 0.5 Hz with a time constant t f = 2.5 sec and N exc = 1500.

Cerebral blood flow
Following in vivo measurements in rodents [50,51], the cerebral blood flow is modeled as a piecewise double exponential function delayed in time by t 1 relatively to the onset of stimulation t 0 . It reads: with F 0 = 0.012 sec -1 [46]. Typical values are t 0 = 0 sec and t 0 = 1 sec, t end being the time at which stimulation ends. Two distinct simulation scenarios are considered to mimic in vitro and in vivo conditions. In the in vivo scenario, Equation (7) is used while in the in vitro scenario, the capillary state variables remain constant at their steady-state while the rest of the variables are left free to vary. Our simulations have shown that this is almost equivalent to taking a constant blood flow F(t) = F 0 .

BOLD signal
The blood-oxygen-level-dependent (BOLD) signal is computed following [52]. It is written as a function of the deoxyhemoglobin concentration (dHb) and of the venous volume (V v ): with dimensionless parameters k 1 = 2.22, k 2 = 0.46 and k 3 = 0.43 [40]. The steady-state values of deoxyhemoglobin (dHb 0 ) and venous volume (V v,0 ) are given in Table A1.

Optimization procedure
As noted already by Aubert and colleagues [53], most models of energy metabolism concentrate on erythrocytes, muscles or other organs such as the liver. It is also not clear whether or not parameters drawn from experiments could be directly injected as such into a model without spatial dimensions and without diffusion processes like ours. To circumvent this problem, we proceeded as follows: First, we chose target steady-state values for the concentration of metabolites following measures reported in the literature. Specifically, we chose the concentration of intracellular sodium following [54], the concentration of intracellular glucose, phosphoenolpyruvate, pyruvate, adenosine triphosphate and phosphocreatine following [55], and glyceraldehyde-3-phosphate following [56]. Finally, the NADH concentration in all four compartments where it appears in the model was chosen following [56] and calculations based on results by Kasischke and colleagues [23].
We then optimized a subset of model parameters (see Table 3) by fitting its predictions to the temporal dynamics of NADH fluorescence as measured by Kasischke et al. [23]. Namely, the dynamics of the NADH concentration in various compartments was extracted empirically from Fig. 4D in [23]. Data points were then fitted with sums of exponentials in order to obtain continuous curves. We then optimized the model by minimizing the distance between the temporal dynamics of NADH in the model and the one in the smoothed curve obtained from [23] using least-square distance as the error measure and using the downhill simplex algorithm. After the optimization converged, we rounded the value of the optimal parameter set and recomputed the steady-state value. All along optimization, we checked that the steady-state was stable by computing its Jacobian matrix (first order approximation) [45]. The parameter set in Table 3 is the set resulting from this procedure.

Numerics
The simulations were run in MATLAB (The Mathworks, Natick MA, USA). The model was integrated with the ordinary differential equation solver with fixed and optimized parameters (ode15s) that is adapted to stiff systems. We used a time step @t = 10 -4 sec when the neuron is spiking and @t = 1 sec starting one second after the end of presynaptic stimulation. @t = 10 -4 sec is smaller than the fastest time constant appearing in the Hodgkin-Huxley equations [tau h (-80 mV) = 6.4·10 -4 sec]. The second time step (@t = 1 sec) is small enough for the slow metabolic processes and maintains simulation time and memory usage to reasonable values for an average desktop PC. Simulations take a couple of minutes to execute on a recent laptop.

Results
We developed a model of the coupling between neuronal activity and metabolic response in neurons and astrocytes. The model employed to simulate the neural-glial-vascular (NGV) functional system is composed of four distinct computational units representing a neuron, an astrocyte, a capillary and the extracellular space (Fig. 1). The core of our model is composed of the compartmentalized model of brain energy metabolism recently proposed by Aubert and Costalat [13,40]. This model connects a model of erythrocyte glycolytic metabolism [45,56] together with the so-called "Balloon model" of blood flow dynamics [46]. From this starting point, we added a precise description of neuronal membrane excitability formulated within the Hodgkin-Huxley framework [57]. Channels dynamics is drawn from a model proposed by Wang [44]. It includes all the standard Hodgkin-Huxley currents plus a high-threshold calcium current and a calcium-gated potassium current inducing spike-frequency adaptation. The Hodgkin-Huxley model is connected to the metabolic pathways through the electrogenic Na, K-ATPase pump which is responsible for a net outward current and concomitant ATP consumption. We modified the metabolic pathways to include compartmentalization of NADH between the cytosol and mitochondria. To do so, we developed a very simple model of mitochondrial respiration and added NADH malate-aspartate shuttles between the cytosol and mitochondria, drawing inspiration from a model by [58]. Finally, the model is driven by external input modeled as a global excitatory presynaptic activity and coordinated increase of the cerebral blood flow. The presynaptic population is coarsely described through a time-dependent excitatory conductance. This conductance drives sodium flow in neurons through AMPA receptors and action potential-generating voltage-gated sodium channels, and in astrocytes through excitatory amino acid transporters which co-transport glutamate using the sodium gradient. The model is illustrated in Fig. 1 and extensively described in the Methods section.
We first tested the model for its responsiveness to an excitatory stimulus ( Fig. 2A) and recorded its voltage response (Fig. 2B). In response to this stimulus, the model generated action potentials within the initial 7 sec of the stimulation. Because of spike-frequency adaptation and of the time course of the excitatory stimulus, the frequency of elicited action potentials quickly decreased until the neuron eventually ceased to fire (Fig. 2B inset). Response trajectories of intracellular sodium, in both the astrocyte (red) and neuron (blue), showed significant differences in both amplitude and duration, with astrocytes exhibiting a smaller but more sustained response and a delayed recovery (Fig. 2C).
We then examined the time course of critical intermediates in energy metabolism in response to the same excitatory stimulus. We first focused on the concentration changes of adenosine triphosphate (ATP) and phosphocreatine (PCr) in the glial and neuronal compartments (Fig. 3). In both cases, the response to the excitatory stimulation, evidenced as a consumption of these energy rich metabolites, was slower in the astrocytic (red) than in the neuronal compartment (blue) (Fig. 3A). And while the decrease in glial ATP surpassed that seen in the neuron (Fig. 3C), the consumption of PCr predominated in the neuron (Fig. 3B). In both compartments, the resulting decrease in ATP concentration was very limited.

NADH responses evoked in vitro are consistent with the ANLS hypothesis
Next, we examined the trajectories of nicotinamide adenine dinucleotide (NADH) in response to the stimulation in the astrocytic, neuronal and mitochondrial compartments (Fig. 4A). The dashed lines represent experimental data from Kasischke et al. [23]. Fig. 4A shows the temporal evolution of the concentration of NADH in the astrocytic cytosol, in the neuronal mitochondria and averaged over the whole tissue as evoked by a 20 sec stimulation episode (see Methods). All three curves are in excellent quantitative agreement with the results reported by Kasischke et al. [23]. In particular, the NADH concentration in the neuronal mitochondria displays an initial dip of about-10% indicating a strong increase of the oxidative metabolism in neurons (Fig. 4B). It then returns towards its baseline before the presynaptic bombardment has finished and finally slightly overshoots in the poststimulus period. On the contrary, the NADH concentration in the astrocytic cytosol increases significantly only about 10 sec after the onset of the stimulation and displays a long-lasting monophasic behavior. This corresponds to a strong and sustained increase of the glycolysis in this compartment (Fig. 4B). The initial dip in the neuronal mitochondria is the result of consumption of NADH to produce ATP. A recovery and rebound results when NADH is produced from the consumption of lactate imported into the neuronal cytosol from the extracellular space (Fig. 4C). In both Fig. 4B and C, it can be seen that the astrocytic response is slower than the neuronal response. In particular, both oxygen and glucose consumption increase immediately at the beginning of the stimulation in the neuronal compartment. Partially supporting this metabolic activity, neurons immediately start to import lactate from the extracellular space (Fig. 4C). On the contrary, the increase in glucose consumption by astrocytes (Fig. 4B) is more gradual and the increase in lactate export by astrocytes to the extracellular space is slightly delayed (Fig. 4C). The initial release of presynaptic glutamate with subsequent neuronal activity and reuptake into astrocytes lead to the increase in intracellular sodium concentration and activation of the Na, K-ATPase imposing, along with the conversion of glutamate to glutamine in astrocytes, an increased metabolic demand. However, as can be seen in Fig. 2C, the increase in intracellular sodium is slower and more gradual in the astrocytic compartment leading to the 10 second delay in the glial cell metabolic response to the stimulation. Finally, the dynamics of tissue NADH (Fig. 4A) is mirrored in the predicted tissue and extracellular lactate concentrations (Fig. 4D). The utilization of glucose and oxygen by neurons and astrocytes during this 20 sec stimulation episode is shown in Fig. 4B. For model optimization, we imposed that the largest fraction of glucose goes to astrocytes while the largest fraction of oxygen goes to neurons [42]. We observed that this bias is further increased during stimulation. The neuronal oxygen utilization immediately increased at the onset of stimulation in register with the initial dip of the NADH in the neuronal mitochondria. This is consistent with reports that the astrocytic fraction of glucose utilization increases during stimulation [19].  Fig.  2 and 3). The dotted lines indicate corresponding in vitro data reproduced from Kasischke et al. [23]. B. Oxygen utilization (thick lines) and glucose utilization (thin lines) by neurons (blue) and astrocytes (red) as evoked by the same 20 sec stimulation episode as in A. Glucose utilization is multiplied by 6 to allow a direct comparison to oxygen utilization. C. Net transport of lactate between the four compartments of the model during the same 20 sec stimulation episode as in A. Lactate is exported by the astrocyte to the extracellular space (thick red line) and imported by the neuron from the extracellular space (thick blue line; net import is negative by convention in this model). The thin red line denotes the activity of the lactate dehydrogenase converting pyruvate into lactate in the astrocytic cytosol while the thin blue line denotes the activity of the lactate dehydrogenase converting lactate into pyruvate in the neuronal cytosol (again, the negative sign is a convention). Note how increase in lactate-to-pyruvate conversion precedes the increase in net lactate import by neurons, while pyruvate-to-lactate conversion follows the increase in net lactate export by astrocytes. D. Relative fluctuations of tissue lactate (pink) and of extracellular lactate (black) during the same 20 sec stimulation episode as in A to C. Both lines are superimposed and barely distinguishable doi:10.1371/journal.pcbi.1004036.g004 The biophysical model correctly predicts rodent in vivo oxygen and lactate evoked responses One of the hallmarks of a successful model is its ability to reproduce and explain empirical observations. We thus now turn to an in vivo situation and compare the predictions of the mathematical model we designed in the precedent sections to two experiments carried out in rats. Tissue oxygen and lactate during stimulus ( Fig. 5B and C), and lactate transfer between compartments were compared (Fig. 5D).
Upon stimulation, CBF increases after a delay of~1 sec (see Equation 7), quickly peaks before relaxing to an elevated plateau. This pattern matches neurovascular responses observed in Figure 5. Predicted evoked responses of tissue lactate and of tissue oxygen in vivo in rodents during a 60 sec stimulation episode. A. Temporal evolution of the cerebral blood flow chosen as an input to the model during a simulated 60 sec stimulation episode in vivo (grey area). This specific time course closely matches measurements in rodents during functional forepaw or whisker stimulation [50,51]. Note that the blood flow only significantly increases approximately 1 sec after the onset of activation (inset). B. Relative fluctuations of the intra-parenchymal oxygen concentration as evoked by a 60 sec stimulation episode with the blood flow as in A. These results closely match the experimental results of Ances et al. [50] (their Fig. 1). C. Relative fluctuations of tissue lactate (pink line) and of extracellular lactate (black line) during the same 60 sec stimulation episode as in A to B. These results closely match the experimental results of Hu and Wilson [16] (their Fig. 1). D. Net transport of lactate between the four compartments of the model during the same 60 sec stimulation episode as in A to C. Like in the in vitro case (Fig. 4), lactate is exported by the astrocyte to the extracellular space (thick red line) and imported by the neuron from the extracellular space (thick blue line; net import is negative by convention in this model). A small amount of lactate is exported from the extracellular space to the capillary at baseline and this export increases by 69% after the end of the stimulation (pink line). The thin red line denotes the activity of the lactate dehydrogenase converting pyruvate into lactate in the astrocytic cytosol while the thin blue line denotes the activity of the lactate dehydrogenase converting lactate into pyruvate in the neuronal cytosol (again, the negative sign is a convention) doi:10.1371/journal.pcbi.1004036.g005 Figure 6. Predicted evoked responses of tissue lactate, CMR glc , CMR O2 , oxygen-glucose index (OGI) and BOLD in vivo in humans. A. Temporal evolution of the cerebral blood flow chosen as an input to the model during a simulated 900 sec stimulation episode in vivo (grey area). This specific time course closely matches in vivo measurements in humans during imaging experiments [52,60]. B. Relative fluctuations of tissue lactate concentration during the same stimulation episode as in A. The model predicts an initial lactate dip followed by a 60% increase sustained till the end of the stimulation. The presence of a dip matches experimental data from Mangia et al. [17]. C and D. Cerebral metabolic rate of glucose consumption (CMR glc ) and cerebral metabolic rate of oxygen consumption (CMR O2 ) during the same 900 sec stimulation episode as in A to B. In both cases, the light area corresponds to the contribution of the astrocytic compartment towards the total tissue consumption, while the dark area corresponds to the contribution of the neuronal compartment. While glucose consumption increases by about 40%, the increase is mostly due to the astrocytic compartment (light red) with the neuronal glucose utilization even slightly decreasing at the onset of activation (dark red). On the contrary, while oxygen utilization increases by about 10%, most of this increase is due to the neuronal compartment (dark blue) with the astrocytic oxygen utilization being almost constant (light blue). E. The predicted ratio of CMR O2 to CMR glc or oxygen-glucose index (OGI) during the same 900 sec stimulation episode as in A to D. rodents in response to sustained sensory stimulation (for instance by mechanical activation of the whiskers [50,51] (Fig. 5A)). The 1 sec delay seen in panel A is hardcoded in the model (see Equation 7). This matches our own observations that CBF only starts to increase above its baseline~0.5-1 sec after the onset of stimulation [51,59].
Because neuronal activity slightly precedes functional hyperemia, the oxygen concentration initially dips below its resting value before rising as cerebral blood flow finally increases after the 1 sec delay (Fig. 5B inset). The dip in oxygen below baseline levels after stimulation has ceased reflects the rapid decrease in the replenishment rate by blood at a time when oxygen is still consumed to replenish the ATP that is used to fuel the Na, K-ATPase pump.
Extracellular lactate is consumed throughout the stimulation. Its concentration initially dips until the cerebral blood flow increases and leads to a sustained overproduction of lactate (Fig. 5C). At rest, the astrocytic compartment exports lactate to the extracellular space, part of which is taken up by the neuronal compartment for energy production (Fig. 5D), the rest being exported to the circulation. Upon stimulation, export of extracellular lactate to the circulation is reduced while import into neurons is increased. Export of lactate to the extracellular space by the astrocytic compartment is also increased but with a delay and explains the initial dip in concentration. In the recovery period, all transports slowly return back to their baseline values explaining the long lasting overshoot of extracellular and tissue lactate. Export of extracellular lactate to the circulation is durably increased in the recovery period ( Fig. 5D; pink line). Fig. 5 shows results of simulations independent from the simulations that yielded Fig. 2 to 4 where the model was constrained to reproduce the experimental results from ref. [23]. In these new simulations, not only is the temporal course of tissue oxygen qualitatively predicted, but the amplitude of fluctuations is, to some extent, quantitatively predicted as well. Our model predicts that the oxygen pressure first drops by −1.7% (inset), then overshoots at +17.7% before stabilizing 2.4% above its baseline in the last 20 sec of the stimulation. Finally, it undershoots to −2.6% in the post-stimulus period. These figures are to be compared with the values reported by [50], namely, an initial drop at−1.8%, an overshoot at +19.9%, a stabilization 1.7% above the baseline and a final undershoot at −3.1%.
The consumption of lactate closely tracks the stimulation dependent oxygen consumption but lacks the inflections corresponding to CBF changes as it is less affected by the blood flow (Fig. 5D). Lactate is exported by the astrocyte to the extracellular space (thick red line) and imported by the neuron from the extracellular space (thick blue line; net import is negative by convention in this model). A small amount of lactate is exported from the extracellular space to the capillary at baseline and this export increases by 69% after the end of the stimulation (pink line). The thin red line denotes the activity of the lactate dehydrogenase converting pyruvate into lactate in the astrocytic cytosol while the thin blue line denotes the activity of the lactate dehydrogenase converting lactate into pyruvate in the neuronal cytosol (again, the negative sign is a convention). These net transfers all contribute to the evolution of the tissue and extracellular lactate concentrations (Fig. 5C). This prediction closely matches the experimental results of Hu and Wilson [16] (their Fig. 1).

The biophysical model correctly predicts human glucose and oxygen utilization during brain activation
We then compared new simulations to measurements from human subjects [26,52,60]. As in Fig. 5, the neurovascular response was adapted from experimental measurements (Fig. 6A). Simulated levels of lactate (Fig. 6B), the cerebral metabolic rate for glucose (Fig. 6C), the cerebral metabolic rate for oxygen (Fig. 6D), the ratio of cerebral uptake of O 2 to cerebral uptake of glucose (Oxygen-Glucose Index or OGI) (Fig. 6E), as well as the blood-oxygen-level dependent signal (Fig. 6F) all supported a validation of the model.
The model qualitatively and to some extent quantitatively predicted well-known features of these various macroscopic observables, including the BOLD signal [52,60], despite extrapolation from multiple sources. For instance, after an initial dip [17], the tissue lactate concentration reached a plateau that last until the end of the stimulation [61]. The initial dip in extracellular lactate (Fig. 6B inset) is representative of a surge in lactate consumption at the beginning of neural activity, and a symptom of the intrinsic latency in the start-up of the ANLS and of functional hyperaemia. The time lag between the onset and offset of neuronal activation and the onset and offset of CBF also explain other transients such as the slow recovery of the lactate concentration after the cessation of neural activity [61]. Relative increase in glucose and oxygen consumption also matched experimental results and led to a decreased OGI during stimulation [26].

Discussion
There is mounting evidence in support of a metabolic link between neurons and glial cells. The most prominent experimentally-based conceptual model of neuron-glia metabolic coupling, the astrocyte-neuron lactate shuttle (ANLS), has raised controversies, not so much for the proposed energy-dependent link between the two cell types, but because certain details and mechanisms are debated [11,62]. Although experimental evidence gathered over the last two decades largely supports it [11], experimentally untangling this system is challenging.
Our detailed biophysical model of the NGV ensemble expands on previous models in four distinct ways. First, a shuttling of lactate from astrocytes to neurons emerged in response to activation. Second, the model is consistent with increased neuronal oxidative metabolism and delayed increased astrocytic glycolysis for generating the activity-dependent NADH transients. Third, the model correctly predicts the dynamics of tissue lactate and oxygen as observed in vivo in rats. Fourth, the model correctly predicts with good quantitative precision the temporal dynamics of tissue lactate, CMR glc , CMR O2 and of the BOLD signal as reported in human studies. These findings not only support the ANLS hypothesis but also provide a quantitative mathematical description of the metabolic activation in neurons and astrocytes, as well as the macroscopic measurements obtained with brain imaging techniques.
The blood oxygen-level dependent (BOLD) signal which forms the basis of the functional magnetic resonance imaging (fMRI) technology reports fluctuations in brain activity, the molecular and cellular mechanisms of which are still incompletely understood [8,63]. Although somewhat limited in representing detailed processes, our use of the Buxton model was sufficient to correctly predict known features of the BOLD signal (Fig. 6F). Ultimately, modeling efforts that build on our work will have to include more detailed descriptions of blood flow regulation. For instance, blood flow regulation in the brain was recently suggested to happen in the microvasculature at the capillary level by active dilation of pericytes [64]. Subsequent efforts will need to focus on the daunting task of modeling the numerous pathways that relate neuronal activity to functional hyperaemia [65].
Recent modeling of the neuron-astrocyte cross-talk during oscillations linked to blood oxygenation levels verified the possibility that the slow fMRI BOLD signals might reflect the spontaneous ongoing activity of neuroglial networks [66]. Our results support this view by accurately modeling results from human imaging experiments [52,60]. A course for future modeling will be to examine and model data from multi-modal imaging experiments [67].
As noted in the introduction, controversy surrounds the directionality of lactate flow in the brain with a neuron-to-astrocyte direction proposed by some. Here, following the arguments delineated in Jolivet et al. [42], we imposed that the largest fraction of glucose should go to astrocytes (there is no controversy that neurons are responsible for the vast majority of oxygen consumption). We derive confidence in our model from the fact that it correctly predicts a vast array of in vivo experimental findings while being only loosely constrained by in vitro experimental findings and by the imposed compartmentalization of glucose uptake between neurons and astrocytes. As argued in Jolivet et al. [42], metabolic shuttling between the astrocytic and neuronal compartments originates from the imbalance between the high oxygen consumption of neurons and their limited glucose utilization for ATP production purposes. Unlike astrocytes, neurons are unable to up-regulate their rate of glycolysis in response to increased activity due to constitutive inhibition of the rate limiting glycolytic enzyme Pfkfb3 (6-phosphofructo-2-kinase/fructose-2, 6-bisphosphatase 3) [68]. This mechanism is crucial to neuronal defenses against reactive oxygen species. Additional in vitro mechanisms support that compartmentalization includes the existence of glutamate-induced glycolysis in astrocytes [5,69], glutamate-induced inhibition of glucose transport in neurons [70] and the involvement of extracellular increases in potassium inducing an Na, K-ATPase-dependent activation of glycolysis in astrocytes [71]. However, it is to be noted that if the proportion of glucose directly taken up by neurons was to be increased, the astrocyte-to-neuron lactate shuttle would be reduced in amplitude, and its direction eventually reversed if neurons were consuming glucose in excess of what they oxidize (see [42] for further discussion of this question).
Metabolic phenotypes have been suggested that convey a metabolic identity depending on how they utilize the various oxidative and non-oxidative pathways [14]. The re-equilibration of the constituents of extracellular space, a process involving physical flushing mediated by astrocytes, is also now thought to be one of the key housekeeping functions of sleep [72] and the disruption of normal metabolic processes is suggested to underlie the progression of neurodegenerative diseases such as Alzheimer's [73]. Lactate, whether sourced from glia or plasma, is associated with neuroprotection [32][33][34]. Assuming the ANLS is indeed taking place, we are left to question: What is it good for?
Neurons do not appear to suffer functional consequences as a result of their metabolic peculiarities. Lactate can sustain prolonged firing in neurons more efficiently than glucose in culture and can preferentially support activity in both resting and active states in vitro and in vivo [30,31,74]. Recently, it was shown in the subfornical organ that this pathway can also affect the dynamics of the local neural network by modulating the excitability of GABAergic neurons through the regulation of ATP-dependent potassium channels [36]. It is not necessary perhaps to preclude the use of glucose by both neurons and glia under certain circumstances as has been suggested by a computational model studying the ATP supply to neurons under hypoxic conditions [75], and as is indeed suggested by our own results (see Fig. 4B, 5D and 6C). Finally, extracellular lactate might also act as a messenger to the vasculature [28] and it is thus possible that the ANLS plays a role as one of the pathways regulating functional hyperaemia.
Astrocytes also support memory formation by supplying neurons with lactate [76]. So central is the ANLS to the normal function of the brain that learning, as measured by LTP and long term memory formation in the hippocampus of rats, is abolished by interfering with the transport of lactate from astrocytes to neurons [76]. Consistent with those findings, lactate, but not glucose, has been show to induce the expression of plasticity genes such as Arc, Zif 268 and BDNF in vitro in neurons and in vivo [77]. Further, the mechanism by which glucose enhances memory storage has been shown to involve the neuronal consumption of lactate [78].

Summary
We present here the first temporal multi-scale model of the NGV that accurately reflects experimental observations in multiple settings and organisms. These findings not only support the ANLS hypothesis but also provide a quantitative mathematical description of the metabolic activation in neurons and astrocytes, as well as of the macroscopic measurements obtained with functional brain imaging techniques.