Computational singular perturbation analysis of brain lactate metabolism

Summary

Lactate in the brain is considered an important fuel and signalling molecule for neuronal activity, especially during neuronal activation. Whether lactate is shuttled from astrocytes to neurons or from neurons to astrocytes leads to the contradictory Astrocyte to Neuron Lactate Shuttle (ANLS) or Neuron to Astrocyte Lactate Shuttle (NALS) hypotheses, both of which are supported by extensive, but indirect, experimental evidence. This work explores the conditions favouring development of ANLS or NALS phenomenon on the basis of a model that can simulate both by employing the two parameter sets proposed by Simpson et al. (J Cereb. Blood Flow Metab., 27:1766, 2007) and Mangia et al. (J of Neurochemistry, 109:55, 2009). As most mathematical models governing brain metabolism processes, this model is multi-scale in character due to the wide range of time scales characterizing its dynamics. Therefore, we utilize the Computational Singular Perturbation (CSP) algorithm, which has been used extensively in multi-scale systems of reactive flows and biological systems, to identify components of the system that (i) generate the characteristic time scale and the fast/slow dynamics, (ii) participate to the expressions that approximate the surfaces of equilibria that develop in phase space and (iii) control the evolution of the process within the established surfaces of equilibria. It is shown that a decisive factor on whether the ANLS or NALS configuration will develop during neuronal activation is whether the lactate transport between astrocytes and interstitium contributes to the fast dynamics or not. When it does, lactate is mainly generated in astrocytes and the ANLS hypothesis is realised, while when it doesn’t, lactate is mainly generated in neurons and the NALS hypothesis is realised. This scenario was tested in exercise conditions.

Keywords:

Glucose metabolism – Algorithms – Evolutionary rate – Neurons – Glucose – Reaction dynamics – Reactants – Astrocytes

Introduction

The human brain uses glucose as its main source of energy. Even though the adult brain accounts for ∼2% of the body weight, it consumes ∼20% of glucose-derived energy. Brain neurons have the highest energy demand [1] and the brain becomes dependent on glucose for proper functioning.

Glucose is metabolized to pyruvate in the presence oxygen, in both astrocytes and neurons. Pyruvate can follow two different pathways: oxidative or aerobic (producing energy in the mitochondria, in the presence of oxygen) and non-oxidative or anaerobic (producing energy and lactate in the cytosol, before the involvement of oxygen). Even though brain metabolism has been the subject of interest and research for decades, there is still a lot of discussion on to what extent energy derived from glucose is consumed by the astrocytes and the neurons themselves, and how much is shuttled between astrocytes and neurons in the form of lactate [2]. Two different hypotheses were introduced for the shuttling of lactate; one advocating the flow of lactate from astrocytes to neurons and another one advocating the flow of lactate from neurons to astrocytes. These two hypotheses are supported by experimental data, which are indirect “because current technology does not have adequate spatiotemporal resolution to quantify metabolic activity in single cells in vivo” [3].

The astrocyte-to-neuron lactate shuttle (ANLS) was proposed on the basis of glutamate-evoked increases in glucose utilization and lactate release by cultured astrocytes [4, 5]. According to ANLS, the neuronal activity triggers a glucose uptake in astrocytes, which leads to large increases in production of lactate. Lactate is then released to the extracellular space and transports to nearby neurons, where it is used as substrate for energy production (for detailed review see Ref. [5]). The validity of the ANLS hypothesis has been supported by experimental evidence over the years; many of them can be found in Refs. [6, 7].

On the other hand, the neuron-to-astrocyte lactate shuttle (NALS) states that lactate is generated by neurons and taken up by astrocytes during activation. Due to the high energy demands in neuronal activation, neurons take up glucose and transfer lactate to astrocytes, since the neuronal glucose transporter, GLUT3, appears to have higher transport rate compared with the astrocytic glucose transporter, GLUT1 [8–12]. Similarly to the ANLS hypothesis, there is experimental evidence that supports the NALS hypothesis [13–15].

In the context of the ANLS hypothesis, neurons are unable to increase their glycolytic activity in response to neuronal activation [16–18]. Therefore, enhanced glycolysis due to neuronal activation can take place only in astrocytes. On the other hand, in the context of the NALS hypothesis, neuronal activation enhances glycolysis in neurons [19–21].

A large portion of the scientific community is currently split between the validity of either the ANLS or NALS hypothesis, while another portion advocates that both hypotheses can be valid depending on the operating conditions. There is a large body of literature advocating the validity of ANLS [6, 22–32], as this is supported by indirect experimental findings [4, 18, 32–44]. The main criticism of the ANLS hypothesis stems from the NALS supporters [8, 9, 45–54], and the related indirect experimental evidence [10, 13–15, 55–64]. On the other hand, recent findings suggest that perhaps brain responds differently to different environmental challenges, depending on the availability of the fuel source. For example, it was shown that under inadequate supply of glucose, alternative substrates contribute to brain energetics [7, 11], such as lactate and/or ketones during exercise, hypoxia, hypoglycemia or intense brain activity [65, 66]. Therefore, under different brain environmental conditions either the ANLS or NALS hypothesis can be valid. The identification of the conditions under which ANLS or NALS manifests is one of the objectives of this manuscript.

The very beginning of the brain chemistry can be traced back to 1884 [67] but early evidence on metabolism in brain tissue can be found in Himwich’s early work [68] dated in 1951. Since then, the brain energy metabolism network has been investigated extensively mainly through in vitro and in vivo experiments and measurements on rats (e.g., Refs. [34, 69–76]) and humans (e.g., Refs. [55, 70, 77–79]).

On the other hand, the development of mathematical models that aim to simulate the metabolic interactions between metabolites into and among the various cerebral compartments has been flourishing very recently. Aubert et al. were among the first to develop mathematical brain metabolism models in order to investigate the connections between electrical activity, energy metabolism and hemodynamics [80, 81]. These models served as baselines for the development of more detailed ones with various investigation objectives like: the ANLS hypothesis [82], the initial rapid decrement in the extracellular lactate concentration [83], the blood oxygenation level-dependent (BOLD) signal in functional MRI (fMRI) [84], the inclusion of glycogen dynamics in astrocytes [85]. The model proposed by Simpson et al. [9] was developed with the purpose of addressing the role of transport in cerebral metabolism and the limitations of the previously developed models related to the concentrations and kinetic properties of the glucose transporters (GLUTs) and monocarboxylate transporters (MCTs) in the brain. With this work, Simpson et al. provided the setup for both the ANLS and the NALS hypotheses. Although, this model was criticized for its simplisity, its increased accuracy compared to in vivo rodent data lead to its utilization (with small parameters’ adjustments) for studying in vivo data from humans at high field functional magnetic resonance spectroscopy (fMRS) by Mangia et al. [10].

Combining the models proposed in [9, 10, 81, 82], Di Nuzzo et al. [8] introduced a new model to investigate the effects of cellular transport, metabolic capacity, specialization for energy use and metabolite trafficking during normal and enhanced cerebral activity. Aiming on the steady state conditions of the metabolic pathways and their relation to neuronal activation, Occhipinti et al. initially proposed a model [86] and used statistical tools to examine the reaction fluxes and the transport rates. Later this model was revised in [87] in order to investigate the preferred energy substrate (i.e., glucose or lactate) in neurons at high activity. They further investigated the reasons of production of γ-aminobutyric acid (GABA), by revising accordingly their latest model in [88]. By modifying and extending the models of Occhipinti et al. [87, 88] for glutamatergic neurons and astrocyte-GABAergic neuron cellular complex, respectively, Calvetti at al. [89] proposed a very detailed steady state model comprising of three cells, each one subdivided into cytosol and mitochondria, where astrocytes, glutamatergic, and GABAergic neurons interact through a common extracellular space (ECS). Finally, the models proposed in [82, 84] were used recently as baseline for the development of a more sophisticated model by Jolivet et al. [22], in order to investigate the ANLS hypothesis and its connection to cerebral blood flow.

Sensitivity analysis and flux analysis are the standard and widely used techniques for the investigation of such complex models. The former technique quantifies the effect of a perturbation in the model’s parameters to the evolution of the system. It is often utilized to obtain insights into the metabolic properties of the system by examining the derivative of fluxes with respect to each model parameters [8, 85], but also to quantify the relative contribution of input parameters on the features of an output attribute, such as BOLD signal [90, 91]. On the other hand, flux balance analysis targets on finding the physiological feasible configuration of reaction fluxes and transport rates capable on maintaining a steady state. It is often utilized to verify the hypotheses on the role of different pathways related to a brain function that have been proposed and to predict the activation levels of pathways at different inputs according to the initial hypotheses [86, 92]. Flux balance analysis can be applied in a mathematical model even when only the stoichiometry is known, and therefore it has been widely used and extended in both deterministic [93, 94] and stochastic frameworks [87, 95, 96]. Flux balance analysis is indeed useful in the interrogation of biochemical systems, but its value is limited since it can lead to conclusions that are only global in character. In addition, sensitivity coefficients (determined through conventional local sensitivity analysis) can be misleading, since a parameter can be related to insignificant sensitivity coefficient and yet be important [97]. Moreover, even in the early 80s works it was acknowledged that there are functional connections between the sensitivity coefficients that they cannot be identified by their mere inspection, thus, additional treatment of the local sensitivity coefficients was proposed [98]. As such, none of these methods can provide details on the temporal evolution of the set of biochemical reactions and metabolites that are most influential on the dynamics and progress of the system [94]. Moreover, they cannot identify the chemical reactions that control the evolution of a metabolite or a variable in general. Such identifications are of paramount importance since they can facilitate the associated physical understanding and eventually can lead to the development of efficient drugs.

Here, the Computational Singular Perturbation (CSP) algorithm will be employed for the analysis of a brain lactate metabolism model. CSP is an algorithmic method for asymptotic analysis developed in the late 1980s by Lam and Goussis [99–101]. It exploits the fact that systems describing the time evolution of biological kinetics models, are driven by processes characterized by a wide range of time scales. In such systems, the processes characterized by the fastest time scales become quickly equilibrated (exhausted). The system is thus confined to evolve within the related surfaces of equilibria that emerge in phase space and is driven by the processes that are characterized by the slow time scales. CSP can provide both mathematical expressions that approximate the generated surfaces of equilibria and the simplified system that approximates the slow evolution along these surfaces. CSP was initially developed to treat large and complex chemical kinetic mechanisms in the context of reacting flows [100, 102], but was later employed for the analysis of other physical problems, like computational mechanics [103], atmospheric science [104], etc., including biological systems [105–108] and pharmacokinetics [109, 110]. CSP exhibits many advantages compared to other similar approaches. Firstly, it is fully algorithmic, requiring no input from the investigator, apart from the detailed model and the accuracy that the simplified model is required to provide. In addition, the CSP algorithm can identify the physical processes in the model that contribute to the emergence of the surfaces of equilibria, those that control the slow (long-term) evolution of the system and those that are responsible for the development of the fast and slow time scales. Finally, CSP can identify the variables whose evolution is characterised by the fast or slow time scales. This understanding allows the investigator to identify the processes that can control not only the slow evolution of the whole system, but also the evolution of each individual metabolite; e.g., see the applications in pharmacokinetics [109, 110].

The CSP algorithm refers to the Geometric Singular Perturbation (GSP) theory, which introduced slow invariant manifolds (low dimensional surfaces on which trajectories evolve according to the slow time scales) and the fast fibers (along which trajectories approach the slow manifold according to the fast time scales) [111–113]. At each point on the manifold, CSP provides basis vectors that span the fast and slow subdomains of the tangent space. These vectors are provided by two iterative procedures; one that improves the accuracy of the basis vectors that are tangent to the manifold and one that improves the accuracy of the basis vectors aligned with the fast fibers. Considering the small parameter ϵ, equal to the ratio of the characteristic fast and slow time scales, each of these iterations increases the accuracy of the CSP basis vectors by O ( ϵ ) [114–117]. With this set of basis vectors being available, the slow manifold is approximated by the relations that result from the negligible projection of the vector field along the fast basis vectors and the system that governs the flow on the manifold is approximated by the projection of the vector field along the slow basis vectors [100]. Due to the multi-scale character and the complexity of the models that are currently of interest, approaches based on the GSP theory are becoming popular in the field of biology [118–122]. The major advantage of CSP is that is algorithmic and does not require the mathematical model to be cast in a non-dimensional form [105–110].

In the current study, CSP and its algorithmic tools are employed in order to analyse brain lactate metabolism, using the model proposed by Simpson et al. [9], as modified by Mangia et al. [10]. This model was selected for the following reasons. Firstly, it is a fairly simple model, thus, an ideal candidate for the demonstration purpose of the proposed method (i.e., CSP and its tools). Secondly, the model of Mangia et al. [10] provides a basis for both the ANLS and the NALS hypotheses. Therefore, both hypotheses can be examined, in view of the same mechanism, which can lead to a fair comparison of the results. Thirdly, it is among the very few models that was fitted against human in vivo data [10, 28, 84], in contrast to the majority of the models that are calibrated against in vitro human data [8, 80–84, 88, 89] or in vivo rat data [9, 28, 85].

In the following, results obtained from the CSP analysis of a simple brain lactate metabolism model will be reported. By investigating the fast and slow dynamics of the model, the couplings among the glucose and lactate pathways are revealed in the context of both the ANLS and the NALS hypotheses. In addition, the rate limiting steps in each of the two pathways are identified. These findings are crucial in order to understand the mechanism for lactate shuttling in the brain, and for the identification of the most conducive point in the pathways for the control of the process. These findings will be demonstrated in normal and exercise conditions of neuronal activation.

Methods

The computational model

The brain lactate metabolic network is analyzed by utilizing the model introduced in Simpson et al. [9], as modified by Mangia et al. [10], which simulates the kinetic behavior of glucose and lactate among the cerebral compartments. As shown in Fig 1 the model consists of 5 compartments: the endothelium (e), the basal lamina (bl), the astrocyte (a), the interstitium (int) and the neuron (n). Glucose (Glc) and lactate (Lac), which are the only chemical species that are accounted for in the model, are oxidized in astrocytes and neurons in order to produce energy. The transport or diffusion of these two substances from one compartment to another is denoted as a reaction. As shown in Fig 1, there are 12 reversible and 6 irreversible reactions in the model; see S1 Text for a detailed description of the model. At this point it suffices to state that the interactions of Glc and Lac through the various compartments is governed by a 10-dimensional system of Ordinary Differential Equations (ODEs):

Computational singular perturbation (CSP) method and algorithmic tools

The Computational Singular Perturbation (CSP) method relies on the multi-scale character of the system under investigation and leads to an algorithm that delivers everything the traditional singular perturbation technique does. As the system evolves, CSP identifies the subspaces in phase space, in which the fast and slow dynamics act. The components of the model that tend to drive the system along the fast directions equilibrate, generating thus constraints in which the system is bound to evolve. This evolution is driven by components that act along the slow directions. With the fast and slow subspaces provided by CSP at each point in time, the components of the model that are mainly responsible for the generated constraints and for driving the system are easily identified.

According to the CSP approach, the vector field g(y) is resolved in N modes, so that Eq (1) is cast in the form [100, 125]:

The M constraints in Eq (3) f^r = (β^r ⋅ S₁)R¹ + … + (β^r ⋅ S_K)R^K ≈ 0 (r = 1, …, M) are the result of significant cancellations among some of the additive terms (β^r ⋅ S_k)R^k (k = 1, …, K). The reactions that contribute significantly to the formation of each of the M constraints are identified by the Amplitude Participation Index (API):

The formation of the M constraints and the dynamics of the slow system in Eq (3) are characterized by the M fastest time scales and by the fastest of the N − M slow ones, respectively. These time scales are approximated by the inverse of the eigenvalues of the Jacobian J, τ_n = |λ_n|⁻¹ (n = 1, …, N). The reactions that contribute significantly to the generation of these time scales are identified by the Time scale Participation Index (TPI):

Each chemical species associates differently to each exhausted CSP mode, e.g., a chemical species can relate mostly to the m-th CSP mode (m = 1, …, M) and much less to the rest. The relation of the m-th CSP mode to the various chemical species is assessed by the CSP Pointer (Po):

The reactions that contribute the most to the evolution of the system within the M constraints, according to the N-dim. system of ODEs in Eq (3), are identified by the slow Importance Index (II):

Results

The dynamics of the system

The evolution of glucose and lactate concentrations in the various cerebral compartments is simulated during neuronal activation by the model described briefly in the previous section and in detail in S1 Text. The variations of the concentrations of metabolites during neuronal activation (electrical and chemical excitability of neurons) have been recorded experimentally by fMRS [10]. It was shown that neuronal activation begins (ends) with a rapid increase (decrease) in lactate and a decrease (increase) in glucose in the human cortex. In order for the model employed to reproduce this response of the concentrations, the parameters in the model were properly adjusted in the periods of neural activation. In particular, the system was allowed to rest on its baseline steady-state conditions for the first 1000 s and then the parameters associated with lactate (glucose) utilization were adjusted in order to generate increased (decreased) levels of lactate (glucose) during the first period of neuronal activation P₁ (1020 s < t < 1240 s). These adjustments were removed at the end of P₁ and the system was let to evolve towards the next period, P₂ (1320 s < t < 1540 s). As in P₁, adjustments in the parameters were introduced throughout P₂ in order to produce decreased (increased) lactate (glucose) levels. After the completion of P₂, the adjustments were removed and the system was allowed to reach steady state conditions. For a detailed explanation of the parameter adjustments during P₁ and P₂ see Ref. [10].

Multi-scale character and modal decomposition

The multi-scale character of the model is established by wide gaps among the fast and slow time scales τ_n and by the ensuing negligible amplitude of the fast modes fⁿ.

There are N = 10 variables in the model employed, so that according to Eq (2) there are 10 amplitudes, which are related to 10 time scales. Fig 2 displays the evolution of the N time scales τ_n that characterise the dynamics of the metabolic network under consideration and the amplitudes fⁿ of the N CSP modes, for both ANLS (top) and NALS (bottom) cases. The period considered includes the periods P₁ and P₂ of neuronal activation, as indicated in the figure. The rapid “step” deviations either in time scale or in amplitude evolution relate to the parameter adjustments during P₁ and P₂. These adjustments have a negligible influence in the time scales, but a significant one in the amplitudes.

In particular, in the left part of Fig 2 it is shown that in both the ANLS and NALS cases the time scale spectrum extends from O(1) to O(10³) s. These time scales are all dissipative (the related eigenvalues are all real and negative), so they are generated by processes that tend to drive the system towards a fixed point. There is a general agreement in the evolution of the time scales in the ANLS and NALS cases. The largest deviations are recorded in the 5th and 8th time scales: τ₅ is somewhat slower in the NALS case, while τ₈ is much slower. The most notable difference in the two cases is that the largest time scale gap is encountered between the 8th and the 9th scales in the ANLS case and between the 7th and the 8th scales in the NALS case. This finding suggests that there is one more constraint established in the ANLS case.

The corresponding amplitudes fⁱ of the CSP modes are displayed in the right part of Fig 2, where the 5 first modes are denoted with dotted lines and the last 5 with solid ones. Before the neuronal activation (t < 1020 s), where the system has already attained its steady state conditions, the two largest amplitudes are O(10⁻²³) and O(10⁻²²) and relate to the slowest 9th and 10th CSP modes respectively; indicating that all CSP modes have been exhausted (i.e., all 10 modes are inactive). As the neuronal activation manifests itself at t = 1020 s, all amplitudes increase. The largest amplitudes during the neuronal activation periods P₁ and P₂ in both ANLS and NALS cases are those of the two slowest 9th and 10th CSP modes, reaching O(10⁻¹⁶), followed by the amplitude of the 8th mode. In the ANLS case, the amplitude of the 8th CSP mode relaxes very quickly after its rapid increase at the start of neuronal activation. In contrast, in the NALS case the amplitude of the 8th CSP mode attains a relatively large value throughout the period of neuronal activation. This is in agreement to the finding that the largest fast/slow time scale gap is encountered between the 8th and the 9th time scales in the ANLS case and between the 7th and the 8th time scales in the ANLS case; i.e., Fig 2 shows that during neuronal activation only the amplitudes of the slow modes are large. After the end of P₂, all amplitudes decay according to the related time scale as the system returns to its steady state; the slowest decay exhibited by the amplitude of the slowest 10th mode.

Our study will focus on the periods P₁ and P₂ of neuronal activation, which are displayed with dashed red arrows in both the left and right panels of Fig 2. Period P₁ relates mostly to the glycolytic uptake in astrocytes, while period P₂ relates to the Lac oxidation in both neurons and astrocytes [10]. Both periods display quite similar dynamical behavior. In particular, in both the ANLS and NALS cases the amplitudes of all modes grow at the start of the two activation periods, but the majority of them relax very quickly. In the ANLS case the number of exhausted time scales is M = 8 for the largest part of P₁ and P₂. Therefore, the 9th and 10th CSP modes are considered to be the slow ones (i.e., the ones that drive the system in P₁ and P₂), while the 1st to the 8th are the exhausted fast ones (i.e., their negligible amplitudes define the confines in which the system evolves). In contrast, in the NALS case, the number of exhausted time scales is M = 7 for the same periods. Now, the 8th, 9th and 10th CSP modes are considered slow, while the 1st to the 7th modes are fast.

In the following, the CSP diagnostic tools discussed previously are utilized in order to acquire physical understanding and derive conclusions about the evolution of the system during neuronal activation. Interested in the long-term evolution of the system, the analysis will focus on the slow dynamics. However, in order to understand the influence of the surface of equilibria, which is generated by the M constraints, on the behavior of the system, a brief presentation of the established constraints is presented first. These constraints involve the most intense components of the system, the equilibration of which allows the less intense ones to drive the system.

Fast dynamics and established constraints

The fast dynamics of the system is responsible for the constraints that are established when the fast CSP modes become exhausted (f^r ≈ 0, r = 1, M). As previously stated, during neuronal activation, the number of exhausted modes is identified to be M = 8 and M = 7 in the ANLS and NALS cases, respectively. These fastest CSP modes were analysed with the CSP tools and the resulting diagnostics are displayed in Table 1 for the ANLS case and in Table 2 for the NALS case (the diagnostics of the non-exhausted 8th mode in the NALS case are included in Table 2 for comparison with the exhausted 8th mode in the ANLS case). The results displayed in the two Tables were computed at t = 1100 s, which belongs to P₁. These diagnostics are indicative of those throughout P₁ and P₂; only 1-2% differences in the CSP diagnostics are recorded. Along with the Pointer Po, the Amplitude Participation Index API and the Time scale Participation Index TPI, the product b^m⋅S_k is displayed in Tables 1 and 2. When b m · S k = O ( 1 ), the stoichiometric vector of the k-th reaction has a significant component along a_m. Since a_m is a fast CSP basis vector, b m · S k = O ( 1 ) indicates that the k-th reaction is a fast one and can potentially contribute to the fast dynamics.

As stated previously, the m-th exhausted CSP mode is distinguished by (i) the reactions that participate in the constraint expressed by the relation f^m ≈ 0 (identified by the API index, Eq (4)), (ii) the reactions that generate the time scale τ_m that characterizes the formation of this constraint (identified by the TPI index, Eq (5)) and (iii) the variables (concentrations) that relate the most to this mode (identified by the Po index, Eq (6)). These CSP diagnostics for the exhausted modes are displayed in Tables 1 and 2.

The main observation from Tables 1 and 2 is that among the eight fastest modes 4 relate to the Glucose path and 4 relate to the Lactate path; i.e., the 1st, 4th, 5th and 8th modes relate to variables and reactions along the Glucose path (Glc-path modes) and the 2nd, 3rd, 6th and 7th ones relate to variables along the Lactate path (Lac-path modes). In particular,

• The pointed by Po variables of the Glc-path modes (Lac-path modes), relate only to Glc (Lac).

• The reactants of the reactions identified by TPI to contribute to the development of the four time scales τ₁, τ₄, τ₅, τ₈ that relate to the Glc-path modes (τ₂, τ₃, τ₆, τ₇ that relate to the Lac-path modes), involve only Glc (Lac).

• The reactants of the reactions identified by API to contribute to the constraints f^m ≈ 0, associated to the four Glc-path modes, involve only Glc. To a large extent, a similar situation arises in the four constraints that relate to the Lac-path modes; i.e., the reactants of almost all contributing reactions involve Lac. The only exception in this case is that in the constraint related to the 6th Lac-path mode the Glc-consuming reaction 7 provides a small contribution in the ANLS case, while the Glc-consuming reaction 9 provides a small contribution in the NALS case.

Therefore, considering the couplings of the two paths in the expressions that approximate the manifold M:

• the Lactate path couples to the Glucose path through the 6th mode; this coupling takes place inside the astrocyte compartment in the ANLS case (reaction 7) and in the neuron compartment in the NALS case (reaction 9) and

• the Glucose path is decoupled from the Lactate path in both the ANLS and NALS cases.

The conclusion is that a perturbation in the reactions on the Lactate path will not alter the Glucose-related constraints. On the other hand, a perturbation in the reactions on the Glucose path will alter the Lactate constraints; as stated previously, this influence is exercised inside the astrocyte compartment in the ANLS case and in the neuron compartment in the NALS case.

Let us now examine the four constraints that develop along the Glucose path, on the basis of the CSP diagnostics displayed in Tables 1 and 2; i.e., those related to the 1st, 4th, 5th and 8th modes in the ANLS case and those related to the 1st, 4th and 5th modes in the NALS case. The 1st (fastest) mode denotes in both ANLS and NALS cases an equilibration between the 2nd and the 4th reactions:

A schematic representation of these constraints that develop along the Glucose path in the ANLS and NALS cases is shown in Fig 3. Clearly, the Glc-related constraints involve all five compartments along the Glucose path and all means of transport and diffusion among them, with the contribution of the 3rd reaction being negligible, since its rate is sufficiently small; see S1 Text.

The four constraints that develop along the Lactate path relate to the 2nd, 3rd, 6th and 7th modes, for both the ANLS and NALS cases. On the basis of the CSP diagnostics displayed in Tables 1 and 2, the following conclusions can be reached. The 2nd mode denotes mainly the equilibration between the forward and backward directions of the 14th reaction, corrected by the net contribution of the 16th and 15th reactions:

A schematic representation of these constraints that develop along the Lactate path is shown in Fig 3. As with the Glc-related constraints, the four Lac-related constraints involve all five compartments along the Lactate path and all means of transport and diffusion among them, with the contribution of the 13th reaction being negligible, since its influence is overshadowed by that of the 14th and 16th reactions; see S1 Text.

The system evolves on a surface of equilibria, approximated by these constraints, on which the slow dynamics become the characteristic ones. Having established, from the data displayed in Tables 1 and 2, that the species related the most to the slow dynamics are Glc^a, Glcⁿ, Lac^a and Lacⁿ (since they relate the least to the fast dynamics), we proceed forward by analyzing the metabolic profiles of neurons and astrocytes.

Slow dynamics

Interested in examining the dynamics of neurons and astrocytes interaction, the slow evolution of the system during neuronal activation is analyzed with the CSP tools. The capability to algorithmically determine the underlying physics of the slowly evolving system, within the constraints imposed by the exhausted modes, is demonstrated in this section.

Characteristic dynamics and driving reactions during the slow evolution

As stated in the previous section, during neuronal activation (periods P₁ and P₂), only the 9th and 10th modes are active in the ANLS case and only the 8th, 9th and 10th modes are active in the NALS case. Given that τ₉ < τ₁₀ in the ANLS case, the slow evolution of the system is characterized by τ₉. However, since τ₈ ≈ τ₉ < τ₁₀ in the NALS case, the slow evolution of the system is characterized by both τ₈ and τ₉. Tables 3 and 4 display CSP diagnostics (Po, TPI) of the slow modes in the ANLS and NALS cases, respectively, computed at t = 1100 s.

It is shown in Table 3 that in the ANLS case the pointed variables and the reactions contributing to the related τ_i of the two slow modes refer to either the Lactate path (9th mode) or to the Glucose path (10th mode). As it will be shown next, this feature is due to the fact that these two modes are practically decoupled. Naturaly, the pointed variables in the two modes are the ones related the least with the fast dynamics, as it is evident by the results displayed in Table 1.

Regarding the 9th mode, the Lacⁿ-consuming reaction 18 and the Lac^a-consuming reaction 17, which represent the Lac oxidation in neurons and astrocytes, respectively, account for more than 50% of the total contribution to the generation of τ₉. On the basis of their negative TPI, it is concluded that these reactions tend to drive the system to a fixed point. On the other hand, the contributions to τ₉ of the two directions of the 16th reaction, which represent the transport of Lac from interstitium to neuron and vice versa, are of opposing nature; the Lac^int-consuming forward direction promoting the movement towards a fixed point and the backward opposing it. Their contributions cancel each other, so that the net influence of the 16th reaction is negligible. The dominance of the 17th and 18th reactions in the generation of τ₉ is reflected in the large Po value of their reactants Lac^a and Lac^int. This dominance implies that the duration of the neuronal activation is mainly controlled by these two reactions. For example, a decrease of the Lac oxidation rate via reaction 18 will lead to the increase of τ₉; i.e., it will lead to a slower evolution of the process. This influence of the 18th reaction will be demonstrated later by properly perturbing R¹⁸.

The results displayed in Table 3 also show that in the context of the 10th mode in the ANLS case, significant contributions to the generation of τ₁₀ are provided by the Glc^int-consuming reaction 4b, the Glcⁿ-consuming reaction 6b and the Glc^a-consuming reaction 5f. Other reactions provide additional contributions, mainly via the already established constraints. Counting the net contribution of each reaction (by adding the contributions of the forward and backward direction), the most significant influence originates from the 1st reaction, followed by the 2nd and then by the 4th.

Table 4 shows that, as in the ANLS case, the pointed variables and the reactions contributing to the related τ_i of the three slow modes in the NALS case refer to either the Lactate path (9th mode) or to the Glucose path (8th and 10th mode). The dynamics of the two slowest modes (9th and 10th) are similar to the equivalent modes in the ANLS case with small variations. However, in contrast to the ANLS case, in the NALS case the 8th mode is accounted in the slow component of the vector field and its dynamics characterise the evolution of the system, along with the dynamics of the 9th mode (since τ₈ ≈ τ₉). It is shown in Table 4 that the largest net contribution in generating τ₈ originates from the Glc^a-consuming 5th reaction, while the net contributions of the 6th and 4th reactions are much smaller. In comparison to the ANLS case, the dynamics of the slowest 10th mode in the NALS case mainly differs in that Glc^a is more pointed and Glcⁿ is less, while the differences in the 9th mode are minor.

Comparing the slow dynamics of the ANLS and NALS cases, Tables 3 and 4 show that the slow dynamics along the Glc-path shifts from neurons towards astrocytes, as it is evident by the sum of the CSP Pointer of Glcⁿ and Glc^a due to the slow modes (D G l c n 9 + 10 = 0.526 and D G l c a 9 + 10 = 0 . 280 in the ANLS case and D G l c n 8 + 9 + 10 = 0.970 and D G l c a 8 + 9 + 10 = 0 . 772 in the NALS case; i.e., a 84.5% and a 175.7% increase of the sum of slow Pointers of Glcⁿ and Glc^a, respectively). This shift of the slow dynamics along the Glc-path is accompanied by an opposite shift of the fast dynamics; i.e. from astrocytes towards neurons. This feature is clearly evident in Fig 3, by the significantly reduced contribution of the Glc^a-related 5th reaction in the established constraints in the NALS case and the ensuing reduced influence of Glc^a in the fast dynamics. In turn, as shown in Tables 1 and 2, this last feature results in a reduced involvement in the fast dynamics of the Glc^a-consuming 7th reaction, that couples the Glc-path to the Lac-path.

The decoupling of the Po and TPI diagnostics along the Glc- and Lac-paths in the context of the slow modes in the ANLS and NALS cases, as indicated by Tables 3 and 4, is reasonable for the Glc-path, since the concentration of Glc in the various compartments is not dependent on that of Lac, as shown in Fig 1. However, this is not that obvious for the Lac-path, since it is influenced by the Glc-path via the 7th and 9th reactions. In order to investigate this issue further, the reactions contributing the most to the evolution of Glc^a, Glcⁿ, Lac^a and Lacⁿ (thereafter, Glc^a,n and Lac^a,n) will be identified, with the use of the Importance Index II. The results displayed in Table 5 identify (i) the reactions that influence most the evolution of Glc^a,n and Lac^a,n and (ii) the corresponding values of II. The positive (negative) sign of the contribution of a reaction displayed in Table 5, denotes its effect in producing (consuming) the related species.

Table 5 shows that in the ANLS case the evolution of Glc^a,n is mainly determined by the R^1f (Glc transport from serum to endothelium), R^1b (Glc transport from endothelium to serum) and R^2f (Glc transport from endothelium to basal lamina). Smaller influences of the rates R⁷ and R^2b towards depleting Glc^a,n are also reported, while the even smaller influences of the two rates R^4f and R^4b cancel each other. Regarding Lac^a,n, Table 5 shows that its evolution is mainly determined by the rates R⁷ (glycolytic rate in astrocytes) and R¹⁸ and R¹⁷ (Lac oxidative rate in neurons and astrocytes, respectively); the first tending to increase Lac^a,n, while the latter two to decrease it. Smaller influences are reported by R^6f and R^6b, which cancel each other. Table 5 shows that R⁹ has a small influence in driving both Glc^a,n and Lac^a,n.

Regarding the reactions that drive Glc^a,n and Lac^a,n in the NALS case, the differences in the results displayed in Table 5 relative to the ANLS case are mainly due to the decreased influence of fast dynamics on Glc^a. As it is shown in Table 5, the evolution of Glc^a is now mainly determined by reactions that either consume it or produce it (R^5f, R^5b, R⁷), while the recorded net contributions by the 1st, 2nd and 6th reactions are smaller. The reactions that mainly determine the evolution of Glcⁿ in the NALS case are similar to those identified in the ANLS case. The main difference is the decreased influence of the 7th reaction and the increased one of the 9th. This is also the main difference in the reactions that drive Lac^a,n; apart from R¹⁸ and R¹⁷ that play a major role in the ANLS case, now R⁹ exhibits the strongest influence, while that of R⁷ is significantly diminished.

From the ANLS to the NALS configuration

The differences between the ANLS and NALS configurations lies in the different values of the parameters of reactions 3, 5, 7, 8, 9, 10, 17 and 18. In order to investigate the influence of these different values, each of the changes in the value of the parameters, leading from the ANLS to the NALS configuration, was implemented separately. It was shown that only the changes in the parameters of the reactions 5, 7, 8, 9 and 10 produce significant deviations from the solution obtained with the ANLS model. The influence of these changes is demonstrated in Fig 4, where the evolution of the net rates R¹⁵ = R^15,f − R^15,b and R¹⁶ = R^16,f − R^16,b is displayed. These two net rates relate to the Lac flow from astrocytes to interstitium and from there to neurons, in the ANLS (A) and NALS (N) cases. As shown in the figure, in the ANLS case R¹⁵ and R¹⁶ are positive, while in the NALS case are negative. In addition, Fig 4 displays the evolution of these two rates in the case where each of the changes in the parameters of the reactions that lead from the ANLS to the NALS configuration is implemented separately; e.g., “R¹⁵ A + 7, 8” indicates the profile of R¹⁵ when the ANLS configuration is considered with the exception of the parameter in the rate of reactions 7 and 8 for which the change that leads to the NALS configuration is accounted for (see S1 Text for a list of all changes in the parameters that lead from the ANLS to the NALS configuration).

It is demonstrated in Fig 4 that the most influential change in the shift from the ANLS to the NALS configuration is produced by the parameter that modifies R⁵, followed by the one modifying R⁹ and R¹⁰. In contrast, the change in the parameter that modifies R⁷ and R⁸ shifts the ANLS configuration away from the NALS one. This finding is in agreement with the substantially different influence of R⁵ in the fast dynamics of the ANLS and NALS cases, discussed previously; i.e., in the ANLS case R⁵ is fast enough that allows the coupling of the Glc and Lac pathways via the 7th reaction (so that the flow of Lac is directed from astrocytes to neurons), while in the NALS case R⁵ is slow enough, so that the coupling between the two pathways is mainly established via the 9th reaction (so that the flow of Lac is directed from neurons to astrocytes).

Assessment of CSP diagnostics

The conclusions reached on the basis of the various CSP diagnostic tools will be assessed here, by perturbing the rates of selected reaction rates. For simplicity, only the ANLS case will be considered; similar conclusions are reached in the NALS case. In particular, consider the rate R¹⁸, which does not contribute to the generation of the eight constraints, as shown in Table 1. In addition, the TPI results in Table 3 suggested that the rate R¹⁸ has a significant contribution to the generation of the characteristic time scale τ₉. Finally, the II results in Table 5 suggested that R¹⁸ tends to decrease Lac^a,n. Therefore, a decrease of the rate constant of R¹⁸ is expected to increase τ₉ and Lac^a,n. Indeed, as it is shown in the left part of Fig 5, a 20% decrease of R¹⁸ causes Lacⁿ, Lac^a and Lac^int to reach a constant value at t ≈ 1805 s, while in the unperturbed case this state is achieved faster, at t ≈ 1235 s. In addition, the levels of Lac^a,n are higher in the perturbed case, as expected. On the other hand, the increase of the rate R¹⁰ by 20% has negligible influence on the evolution of Lac^a,n, as shown at the right part of Fig 5. This is an expected result, since R¹⁰ has a negligible contribution to the eight constraints, to τ₉ and to the slow evolution of the system; see Tables 1, 3 and 5.

A further demonstration of the usefulness of the CSP diagnostic tools emerges from the physical interpretation of the II results displayed in Table 5. In particular, consider the response of Lac^a,n and Glc^a,n when subjected to perturbations of R¹⁸, R⁹, R^1f and R^3f during neuronal activation. Such perturbations might affect the M constraints or the reduced model; see Eq (3). According to the API values in Table 1, R^1f is the only rate from this set that influences the constraints, via the 4th and the 5th modes that relate to Glc^e and Glc^int. According to the II values in Table 5 that assesses the influence of the reactions in the sow model, (i) R¹⁸ has the largest contribution towards decreasing Lac^a,n levels, while has negligible influence in the Glc^a,n, (ii) R⁹ exhibits a small influence in increasing Lac^a,n levels and a smaller influence in decreasing Glc^a,n levels, (iii) R^1f is the major contributor towards increasing Glc^a,n, but has negligible influence in the evolution of Lac^a,n and (iv) R^3f has a negligible contribution to the evolution of all variables.

The predicted influence of perturbations of R¹⁸, R⁹, R^1f and R^3f is demonstrated in Figs 6 and 7. Specifically, in the left part of Fig 6 it is shown that the decrease of R¹⁸ leads to increased Lac^a,n levels, but leaves Glc^a,n unaffected. On the right part of the same figure, it is shown that the increase of R⁹ leads to increased Lac^a,n levels and decreased Glc^a,n levels. Note that the response of Lac^a,n to a perturbation of R¹⁸ is more intense than that of Glc^a,n and Lac^a,n to a similar perturbation of R⁹, in accordance to the related II results displayed in Table 5. Also, in the left part of Fig 7 it is shown that the decrease of R^1f leads to increased Glc^a,n levels, but leaves Lac^a,n levels unaffected. Finally, in the right part of the same figure it is shown that an increase of R^3f has a negligible influence to both Glc^a,n and Lac^a,n, as predicted.

Investigating an exercise scenario

It was shown previously that a decisive factor on whether the ANLS or NALS configuration will develop is whether the 5th reaction contributes to the fast dynamics or not. In particular, it was shown that in the context of the ANLS (NALS) configuration the 5th reaction contributes (does not contribute) to the fast dynamics, allowing thus the major coupling between the Glc and the Lac paths to be established in astrocytes via the 7th reaction (in neurons via the 9th reaction). As a result, the Lac flow is directed from astrocytes to neurons in the ANLS case and in the opposite direction in the NALS case. These conclusions were reached for normal conditions (NC); among others, low Lac concentration in serum (Lac^s = 1 mM).

However, under exercise conditions (EC) the Lac concentration in the serum reaches large values, about 15 − 25 mM (“all-out” maximal exertion) [123, 124]. A high Lac concentration in serum will certainly influence the concentrations along the Lac-path. Nonetheless, this factor will leave the concentrations of Glc unaffected, since the Glc-path is not coupled to the Lac-path. As a result, it is not expected to influence the degree to which the 5th reaction contributes to the fast dynamics; i.e., the conditions under which the Glc-path dynamics favor the development of ANLS or NALS. The only feature that might be affected is the flow direction of Lac between astrocytes and neurons, given the changes in the concentrations along the Lac-path that result by the high Lac concentration in the serum. These issues will be examined next.

In the exercise scenario considered here, the serum Lac concentration was increased to Lac = 20 mM (from Lac = 1 mM in NC). As expected, high serum Lac conditions, introduce significant changes to all Lacⁱ, but leave unaffected those of all Glcⁱ in both ANLS and NALS cases. This is displayed in the left panels of Fig 8, where the ratio of the concentrations computed under exercise conditions (EC) and under normal conditions (NC), y^i,EC/y^i,NC, is displayed. This ratio equals unity for Glcⁱ in all the 5 compartments, while it is slightly larger than unity for the Lacⁱ in all compartments, with the exception of the ratio for Lac^e which is around 6.

The increased concentrations of Lacⁱ modify the flow of lactate between astrocytes and neurons. As shown in the right panels of Fig 8, the high serum Lac concentration causes the rate of Lac-transport from astrocytes to interstitium (R¹⁵ = R^15,f − R^15,b > 0) to decrease and from interstitium to neurons (R¹⁶ = R^16,f − R^16,b > 0) to increase in the ANLS case. In contrast, in the NALS case the Lac-transport from neurons to interstitium (-R¹⁶ > 0) decreases and from interstitium to astrocytes (-R¹⁵ > 0) increases. Fig 8 shows that in both the ANLS and NALS cases the flow of lactate from basal lamina to intertitium (R¹⁴) increases. While in the ANLS case Lac flows from basal lamina to intertitium (R¹⁴ > 0) in both NC and EC, in the case of NALS the minor flow from intertitium to basal lamina (R¹⁴ < 0) in NC is significantly reversed in EC. A schematic presentation of the changes in the net rates R¹⁴, R¹⁵ and R¹⁶ is displayed in Fig 9. Clearly, the increased values of R¹⁶ in the ANLS case and of −R¹⁵ in the NALS case are in agreement with the nature of ANLS and NALS assumptions, respectively; since Lac is transported towards the neurons in the ANLS case and towards the astrocytes in the NALS case increases. However, the decreased values of R¹⁵ in the ANLS case and of −R¹⁶ in the NALS case are in disagreement; since the transport of Lac from astrocytes in the ANLS case and that from neurons in the NALS case decreases. In fact, as shown by the profile of R¹⁶ in the NALS case, at about t = 1300 s, Lac is transported to neurons (R¹⁶ > 0) and not from them, as the NALS hypothesis dictates.

In summary, it was shown that an increased serum Lac concentration resulted in an increased transport of Lac to neurons in the ANLS case, mainly due to the increased diffusion from basal lamina, accompanied by a decreased transport out of astrocytes. In contrast, an increased transport of Lac to astrocytes was recorded in the ANLS case, also due to the increased diffusion from basal lamina, accompanied by a decreased transport out of neurons. It can be shown that these features are in full agreement with the equilibrations that develop in the context of the 2nd and 3rd modes, as they are stated in Tables 1 and 2, which mainly involve the rates R¹⁶, R¹⁵ and R¹⁴ and to a lesser degree R¹³ in normal and exercise conditions and the rate R¹² only in exercise conditions.

Discussion

Fast and slow dynamics influence the response of the system in both the ANLS and NALS frameworks in two distinct manners. The fast dynamics are responsible for the constraints that are established in phase space, within which the system evolves driven by processes responsible for the slow dynamics.

Constraints develop along both the Glucose and Lactate paths. Only one of these constraints couples the two paths; i.e., the Lactate to the Glucose path, while the latter is uncoupled to the former, see Fig 1. Along the Glucose path in both the ANLS and NALS cases the following constraints related to Glc-transport develop:

• from endothelium to basal lamina (R^2f − R^2b) and from there to interstitium (R^4f − R^4b),

• from serum to endothelium (R^1f − R^1b) and from there to basal lamina (R^2f − R^2b), with the contribution of a higher order correction from the Glc-transport in and out of the interstitium (R^4f − R^4b, R^5f − R^5b, R^6f − R^6b); the latter contribution based on the previous establishment of the constraint in (i),

• from interstitium to neurons (R^6f − R^6b), from serum to endothelium (R^1f − R^1b) and only in the NALS case from astrocyte to interstitium (R^5f − R^5b), with the contribution of a higher order correction from the Glc-transport from basal lamina to interstitium (R^4f − R^4b) due to the constraints established previously in (i) and (ii), and

• only in the ANLS case, from astrocytes to interstitium (R^5f − R^5b) and from there to neurons (R^6f − R^6b), with the glycolytic rate in astrocytes R⁷ contributing as a higher order correction.

As Fig 3 shows, these constraints involve all five compartments along the Glucose path. Along the Lactate path in both the ANLS and NALS cases the following constraints related to Lac-transport develop:

• from basal lamina to interstitium (R^14f − R^14b), with the Lac-transport from interstitium via reactions 15 (R^15f − R^15b) and 16 (R^16f − R^16b) providing a higher order correction,

• from astrocytes to interstitium (R^15f − R^15b) and from there to neurons (R^16f − R^16b),

• from astrocytes to interstitium (R^15f − R^15b) and from there neurons (R^16f − R^16b), with a small contribution from reactions 7 and 9 from the Glucose path in the ANLS and NALS cases, respectively, and

• from serum to endothelium (R^11f − R^11b) and from there to basal lamina (R^12f − R^12b), with only a small contribution from reaction 16 which is possible after the equilibria in (i), (ii) and (iii) are established.

As with the constraints related to the Glucose path, Fig 3 shows that these constraints involve all five compartments along the Lactate path.

Clearly, the major difference in the ANLS and NALS cases is the contribution to the fast dynamics of Glc-transport from astrocytes to interstitium (R^5f − R^5b). In the ANLS case this Glc-transport participates in both the fast dynamics and in the ensuing related constraints (large TPI and API in the 4th, 5th and 8th modes; see Table 1). However, in the NALS case this participation is greatly diminished (small TPI and API in the 4th and 5th modes, while the 8th mode is now a slow one; see Table 2). These features suggest that the strong coupling ot the Lac- to the Glc-path that is established in the ANLS case (mainly via the 7th reaction inside the asrtocytes, in the context of the 6th and 8th modes) diminishes significantly in the NALS case (where this coupling is now established mainly via the 9th reaction inside the neurons, in the context of only the 6th mode).

The finding that—due to the behavior of R⁵—the Lac- to the Glc-path coupling is mainly established inside the astrocytes in the ANLS case and inside the neurons in the NALS case determines the conditions promoting the ANLS or NALS development. In the ANLS case R⁵ participates in the equilibria of the 4th and 5th modes, which allow the Glc concentrations in both the astrocytes and neurons to be coupled directly with that in the serum. The formation of Lac from Glc is more intense in astrocytes (R⁷) than in neurons (R⁹), so that the typical for the ANLS case flow of Lac from astrocytes to neurons is generated. In the NALS case the participation of R⁵ to the equilibria of the 4th and 5th modes is significantly reduced, so that only the Glc concentration in neurons is coupled directly with that in the serum. Therefore, the formation of Lac from Glc in neurons (R⁹) dominates, so that the typical for the NALS case flow of Lac from neurons to astrocytes is generated.

The influence of the fast dynamics in determining the development of ANLS or NALS is fully reflected in the slow dynamics that characterize the evolution of the system within the established constraints. In particular, as shown in Table 5, in the ANLS case the Glc^a,n level is mainly determined by the rate of reactions in the first half of the Glc-path (R¹, R²) and the glycolytic rates in astrocytes (R⁷), while the Lac^a,n level is mainly determined by rates in the second half of the Lac-path (R¹⁷ and R¹⁸) and by the rate R⁷. In essence, the glycolytic rate in astrocytes (R⁷) couples the Lac oxidative rates in astrocytes and neurons (R¹⁷ and R¹⁸) with the Lac inflow rate from the serum, via reactions at the start of the Glc-path (R¹, R²). As a result, the flow along the Lac pathway in the ANLS case is directed from astrocytes to interstitium (R¹⁵ = R^15,f − R^15,b > 0) and from there to neurons (R¹⁶ = R^16,f − R^16,b > 0). In contrast, Table 5 shows that in the NALS case the same coupling is now established mainly via the glycolytic rate in neurons (R⁹), while the glycolytic rate in astrocytes (R⁷) exhibits a smaller influence. In this case the flow along the Lac pathway is directed from neurons to interstitium (R¹⁶ = R^16,f − R^16,b < 0) and from there to astrocytes (R¹⁵ = R^15,f − R^15,b < 0). In both cases, the direct coupling of the Lac oxidative rates in astrocytes and neurons (R¹⁷ and R¹⁸) with the Lac inflow rate from the serum via the reactions in the first of the Lac-path (R¹² and R¹³) has no significant influence. Finally, in agreement to the model depicted schematically in Fig 1, Table 5 shows that in both the ANLS and NALS cases Glc^a,n is independent from Lac^a,n.

It was shown that the shift in the Glc to Lac activity along the Glc-path from the astrocytes to neurons, in the NALS case relative to the ANLS case is reflected in both the fast and slow dynamics. This shift is accompanied by a weakening of the Glc^a-consuming/producing reactions 3f, 3b, 5f, 5b, 7 and 8 and a strengthening of the Glcⁿ-consuming reactions 9 and 10, as shown in Fig 10. In the same figure it is shown that in the NALS case the glycolytic rate in astrocytes (R⁷) decreases and the glycolytic rate in neurons (R⁹) increases, relative to the ANLS case, supporting the findings reported previously on the influence of these reactions in coupling the oxidative rates in astrocytes and neurons (R¹⁷ and R¹⁸) with the Lac inflow rate in the serum.

The large Lac concentrations in the serum considered in the context of the exercise conditions, were shown to affect the flow of Lac, but left unaffected the conditions under which the Glc concentrations in astrocytes and neurons are coupled directly with that in the serum, since the Glc-path is decoupled from the Lac-path in the model considered. In particular, it was shown that under exercise conditions the increased rate of Lac diffusion from basal lamina to interstitium modified the rate of Lac transport between astrocytes and neurons; i.e., in the ANLS case contributed to the increased transport to neurons and the decreased transport from the astrocytes and in the NALS case contributed to the increased transport to astrocytes and the decreased transport from neurons. It is noted here that although there were significant changes in the initial conditions and the values of some parameters in the cases of normal and exercise conditions, the basic dynamical features were not altered. Specifically, there were not changes recorded in the reactions contributing to the time scales and the emerging constraints, the variables related to the modes and the number of exhausted modes. Only moderate changes were noticed on the degree to which some reactions contributed to the constraints and to the evolution of the system within these constraints.

The findings reported here are supported by experimental results in the literature. The high glycolytic rate of astrocytes (R⁷) in the ANLS case has been detected experimentally [4, 16, 69, 135, 136], since neurons cannot afford to sustain high glycolytic rate (R⁹) [6, 16]. In addition, the dominance of the glycolytic uptake in neurons (R⁹) has also been detected in the NALS case during neuronal activation [13–15]. Here, the dominance of R⁷ over R⁹ in the ANLS case and of R⁹ over R⁷ in the NALS case, was identified both in the fast dynamics (i.e., in the generated constraints) and the slow dynamics (i.e., in the reactions that drive the system within these constraints).

The results reported in Tables 3, 4 and 5, regarding the TPI and II indices, show that the Lac oxidative metabolism in neurons (R¹⁸) has a larger impact than the oxidative metabolism in astrocytes (R¹⁷) in both the ANLS and NALS cases. This is consistent with the experimental finding that the Lac oxidative metabolism in neurons is more important than the respective metabolism in astrocytes [77, 135]. In addition, the identifications reported here in Table 5, show that the concentration of Lac^a,n is, among others, determined by the rates R⁷ and R⁹, which act towards increasing the rate of change of Lac^a,n. These rates have been reported to increase during neuronal activation in both ANLS (only R⁷) [16–18] and NALS (both R⁷ and R⁹) [19–21] cases, resulting in increasing Lac^a,n concentrations. This result is agreement with the elevated lactate levels that have been reported during neuronal activation [35, 137]. This feature has also been characterised as a moderate to slow process [12], which is also in agreement with our results, due to the aforementioned identifications in the slow modes of the system.

At this point it must be emphasized that the model employed here neglects metabolic processes which are currently considered important. Among them are (i) formation of pyruvate from glucose, (ii) oxidative phosphorylation in mitochondria, (iii) TCA cycle and (iv) consumption/production of ATP and NADH [8, 22, 88]. Clearly, it is possible that consideration of these processes might modify (i) the processes that are mainly responsible for the generation of the fast and slow dynamics and (ii) the conditions for which ANLS or NALS manifest. The degree to which these additional processes influence the response of the system, as it was examined here by analyzing algorithmically the fast and slow dynamics, is currently under investigation.

Finally, the ease by which the CSP tools lead to the determination of the factors that promote the development of ANLS or NALS must be highlighted. CSP allowed for these identifications by identifying the components of the model that are responsible for the fast and the slow dynamics, which are responsible for the development of the constraints in phase space (the fast) and for driving the system within these constraints (the slow).

Conclusions

A computational model describing the evolution of the brain lactate metabolism under neuronal activation was investigated using algorithmic tools of asymptotic analysis. The model introduced in Simpson et al. [9] and modified by Mangia et al. [10] was considered, which conforms to both the ANLS and NALS hypotheses. The analysis was carried out using the CSP method, which provides algorithmic tools for the analysis of multi-scale systems, so it is not hindered by the complexity and/or the size of the system under investigation and its tools can deliver systems-level understanding. CSP was employed in order to: (i) investigate the fast/slow dynamics of the system and the metabolic profiles of neurons and astrocytes during neuronal activation and (ii) compare the dynamics of the ANLS and NALS hypotheses, under normal and exercise conditions.

It was shown that the ANLS or NALS configuration develops depending on whether the Glc rate of transport from astrocytes to interstitium (R⁵) contributes or not to the generated constraints and the fast dynamics. In the ANLS case R⁵ participates in the equilibria of the 4th and 5th modes, which allow Glc in the astrocytes and neurons to be coupled directly with that in the serum. Lac formation from Glc is more intense in astrocytes (R⁷) than in neurons (R⁹), so that the typical for the ANLS case flow of Lac from astrocytes to neurons is generated. In the NALS case the participation of R⁵ to the equilibria of the 4th and 5th modes is significantly diminished, so that only Glc in neurons is coupled directly with that in the serum. Therefore, Lac formation from Glc is mainly taking place in neurons (R⁹), so that the typical for the NALS case flow of Lac from neurons to astrocytes is generated.

As demonstrated, CSP is not hindered by the size or the complexity of the mechanism under investigation, therefore more complicated and larger models of brain metabolism can be considered. In many Central Nervous System pathogenies (like Parkinson’s and Alzheimer’s diseases) or in metabolism-related diseases (e.g., diabetes), where the degeneration of the brain structures causes abnormalities in brain function, glucose and lactate metabolism is altered; see for example [7, 138, 139]. CSP could be effectively employed in order to identify optimal strategies for the control of the system towards desired outcomes, contributing thus to the development of therapeutic treatments.

Supporting information

S1 Text [pdf]
The brain lactate metabolism model.

S2 Text [pdf]
The computational singular perturbation methodology.