Analysis of a mathematical model for brain lactate kinetics.

The aim of this article is to study the well-posedness and properties of a fast-slow system which is related with brain lactate kinetics. In particular, we prove the existence and uniqueness of nonnegative solutions and obtain linear stability results. We also give numerical simulations with different values of the small parameter ε and compare them with experimental data.

1. Introduction. The brain is an organ with high energy needs. While it represents only 2% of the body weight it grabs at least 20% of its total energy needs [9]. The consumed energy can come from many forms such as glutamate, glucose, oxygen and also lactate [3]. Energy is necessary to support neural activity. Gliomas are the most frequent primary brain tumors (more than 50% of brain cancer cases according to the ICM institute). Like other cancers, they lead to alterations of cells' energy management. In particular, lactate creation, consumption, import and export of a glioma cell seem to play a key role in the cancer development [5]- [11]. Today, neuroimaging techniques allow an indirect and noninvasive measure of cerebral activity. It also enables measurement of various metabolic concentrations such as lactate and measurement of important biological parameters such as the relative cerebral blood volume (allowing relative cerebral blood flow calculations). But because energy management in healthy and tumoral cells and glioma growth can be difficult to observe and explain experimentally, we propose to use mathematical modeling to help to describe and understand cells energy changes.
To the best of our knowledge, only a few mathematical models have been proposed to study lactate fluxes in the brain and the interconnections with energy, see [3] for example. We aim herein at analyzing a model first described in [2].
Our paper is organized as follows. We first present the mathematical model proposed to describe the mechanisms of interest. We then investigate its wellposedness and derive bounds on the solutions. Indeed, such an analysis is necessary to justify how a mathematical model is well-adapted to a biological problem. We also analyze the limit model and study its steady state. Bounds on the solutions are important as they are related with the viability domain of the cell. Furthermore, as mentioned in [8], a therapeutic perspective is to have the steady state outside the viability domain where cell necrosis occurs. Additionnally, we present numerical simulations with different values of the small parameter ε and compare them with experimental data. We finally discuss our results.
2. Mathematical modeling. The present model is reduced in order to follow in a simpler way lacate kinetics between a cell and the capillary network in its neighborhood. It is built in vivo which means that we have to consider loss and input terms for both intracellular and capillary lactate concentrations. We denote by u ε the intracellular lactate concentration and v ε the capillary lactate concentration where ε stands for the volume separating the compartments. They are given in mM. One of the main parameters in the model is ε. Indeed to manage the blood flow, vessels dilate and modify their volume. We cannot model this phenomenon in a simple mathematical way. It is thus important to know how variations of their volume correlated with variations of ε impact the whole dynamics.
First, there is a lactate cotransport through the brain blood. It is taken into account by a simplified version of an equation for carrier-mediated symport. This nonlinear term depends on the maximum transport rate between the blood and the cell T > 0 and the modified Michaelis-Menten positive constant for both intracellular and capillary lactate concentrations (k and k respectively).
Then a cell can equally produce and consume lactate, but also export surplus lactate to neighboring cells. We denote by J the balance sheet of the whole phenomenon. The function J is a nonnegative function depending on t and u ε seen as a regulatory term. It is assumed to be bounded by a constant B J and Lipschitz continuous.
Next there is a blood flow contribution to capillary lactate depending on both arterial and venous lactates. We denote by L > 0 the arterial lactate concentration. We also define the blood flow F . The function F is a positive bounded continuous function (F 1 < F < F 2 ) seen as a forcing term.
Finally, we have the following ODE's, for t ∈ R + : ). (2. 2) The initial condition is given by : The model is biologically described in [2], to which we refer the interested readers for a better understanding of this process. Figure 1. Schematic representation of lactate exchanges in a local brain part. There is a cotransport through the brain-blood barrier, a blood flow, cell creation and consumption and interactions between a cell and its neighborhood. Interactions are described in the main text.
Since we have for nonnegative u 1 and u 2 : then we can rewrite (2.1)-(2.2), setting to have ∀t ∈ R + : where H is globally Lipschitz continuous with respect to the second variable. We finally conclude, thanks to the Cauchy-Lipschitz theorem, that we have existence and uniqueness of the solution to the system ∀t ∈ R + .
Bounds on the solution. By means of (2.2), we have ∀t which implies, using Gronwall's lemma, that One can see, using the above formula, that we have Theorem 3.1. We can exhibit a sufficient, but not necessary, condition to ensure a bound on u ε . When it exists, let B J ∈ R be such that : and : In that case, we have, setting z = Bv k +Bv + B J T and ∀t ∈ R + : for which a positive solution exists if and only if B J < T (1− Bv k +Bv ), i.e. (C 2.1 ) holds. From a biological point of view, this condition means that at each time the lactate uptake by a cell (from itself or its neighborhood) cannot be larger than the lactate it can purge through the blood. Otherwise, the cell lactate concentration increase may not be limited.
We set z = Bv k +Bv + B J T and have, thanks to (C 2.1 ) : Let t ∈ R + be such that: Then, which yields We finally deduce that

Remark 2.
This condition on an upper bound on J is sufficient but not necessary. We can actually also consider functions J which do not satisfy (C 2.1 ) but for which, when u ε is large, we can exhibit a better upper bound satisfying (C 2.1 ). It is thus sufficient to write that u ε is bounded from above with a mild (C 2.1 ) condition. For example, we can take J = J test such that, ∀(x, s) ∈ R + × R + : for positive constants G J , L J , C j and ε J such that G J > L J and G J < L J + T k k +Bv . The second condition means that the lactate creation of the cell is smaller than its consumption and purge through the blood, so that the cell is able to manage lactate excess. Then, for x Even though we cannot assert that such a function J is biologically relevant, it is improbable to find relevant functions J leading to a fatal lactate increase in the cell. In fact, it is logical to expect functions decreasing in x, since, when a cell has more substrate than necessary for it to live on its own, it does not have to import or create more.
We have already proved that u ε and v ε are nonnegative functions. We can exhibit lower bounds by using the same method. Indeed, Then, if Similarly, ∀t ∈ R + : Then u ε (t) M v k k leads to u (t) 0. This shows that The upper bound B v on v ε can be derived even when the initial datā v 0 andū 0 depend on ε and are bounded with respect to this parameter. To do so, we use the same method, adapting the final step so that The lower bounds M v and M u can be obtained in the same way for initial data depending on ε.
Stability of the equilibrium. For constant J and F , an equilibrium for (2.1)-(2.2) has been found in [4] : It has also been proven that this unique stationary point is a node, hence a locally stable equilibrium. However, this equilibrium does not always exist. For existence, the parameters need to satisfy : Remark 4. We have already shown that v ε (t) L + T F = B v , ∀t ∈ R + . We verify that the equilibrium v l is such that v l B v , which requires J T . In fact, under (2.7), J > T implies: Therefore, the contraposition leads to: We fix all the parameters but J. Then, we wish to rewrite this condition by giving it in terms of J. In this way, we have : and there is an equilibrium only when J ∈]J b , J h [, where : Knowing that J is nonnegative, there are only two possible cases : 1. 0 < J < J h with one steady-state which is a node, 2. J > J h with no steady-state.
A therapeutic perspective is to have the steady state outside the viability domain [8]. Therefore playing on cell lactate intake could be worth exploring: a large J involves an unbounded cell lactate concentration which leads to an exit of the cell viability domain and, finally, the glioma cell death. 4. The case ε = 0. We now study the limit system for ε = 0, together with the initial condition : We first give some preliminary results and then establish bounds on the solutions and study the well-posedness of the system. We finally compare the original system (with ε > 0) with this limit system (with ε = 0).
Preliminaries. The function v 0 given by (4.2) is defined as long as v 0 (t) belongs We define the function ϕ c , for any constant c > 0, by : It is easy to see that ϕ c is a monotone increasing function with ϕ c (0) = 0. We also define an inverse function of ϕ c : where ψ c is a bijection from ] − c, +∞[ onto R. It can also be a bijection from R + onto itself. Its derivative reads : Employing (4.2), we have : We rewrite (4.1)-(4.2) as : ) := G(t, u 0 (t)), (4.4) and set, for y ∈ [0, +∞[ : A priori bounds on the solutions. Thanks to (4.4), we have : and the system is quasipositive : for an initial conditionū 0 0, there holds u 0 (t) 0, ∀t ∈ R + . A fortiori, we have v 0 (t) = Ψ(u 0 (t)) 0. Using (4.2), we find an upper bound on v 0 , ∀t ∈ R + : We can also obtain an upper bound on u 0 using Theorem 2.1. When it exists, we call it B u,0 .

Lemma 1.
The function Ψ is Lipschitz continuous in u 0 , i.e. there exists K L > 0 such that for all u 1 , u 2 ∈ [0, +∞], Proof. Let u ∈ R + . We know that : and : It follows from the above that Ψ(u) = v, with : Therefore, Ψ(u) = v ∈ [0, B v,0 ] and : Consequently, G(t, u 0 ) is Lipschitz continuous in u 0 . Therefore, thanks to the Cauchy-Lipschitz theorem, we have the existence and uniqueness of the solution to (4.4), ∀t ∈ R + . Finally we have existence and uniqueness for v 0 , recalling that Ψ is a bijection.
Stability of the equilibrium. As proved in [3], (4.1)-(4.2) can have at most one equilibrium given under the above parameters condition. The Jacobian of the system at this point gives the eigenvalue : Therefore, u l is locally stable. Moreover, ∀t ∈ R + , v(t) = Ψ(u(t)), where Ψ is a bijective function. Setting z l = J T + v l k +v l , we have : Thus : and the stationnary point v l is locally stable.

Remark 5.
In the particular casev 0 = Ψ(ū 0 ), we have ∀t ∈ [0, t m ] : Remark 6. Setting v 0 (0) =ṽ 1 ∈ I 1 the negative root of (4.2) for u 0 (t) =ū 0 , we define a second solution to (4.1)-(4.2) for every t ∈ R + . We strongly think that we can study this degenerate system with similar mathematical tools, but this does not make sense from a biological point of view.

5.
Numerical simulations and comparison with experimental data. In this section, we first present several numerical simulations with relevant values of our parameters. We also compare the numerical simulations with different values of J and ε. We finally give and study experimental data. These simulations have been done with the Matlab software.
Numerical illustration with nonconstant J and F . We first consider the system given by (2.1)-(2.2). We expect that a cell manages its lactate concentration by means of its amount but not of the experiment's duration. Therefore we assume that it is biologically relevant to take a function J that does not depend of t. Besides, we expect that a cell imports more lactate when its lactate concentration is low.
In other words J should be monotone decreasing in x. Under this hypothesis we choose : , containing a creation term, a consumption term and an import term. This function J is a bounded and Lipschitz continuous function which enjoys the mild (C 2.1 ) condition. We also define, as given in [1] : The parameters for these two functions are given in Table 1. We also consider the parameters given in [2] and [8]. In that case, (ū 0 ;v 0 ) = (1.15; 1) and the parameters values are given in Table 2. Table 2. Parameters values.
The solutions u ε and v ε remain nonnegative, as shown in Figure 3. Here, we have an upper bound on v ε and J enjoys (C 2.1 ), so that u ε has an upper bound too. Thus, the intracellular and capillary lactate concentrations match with the  Numerical simulations with different values of J and ε. We now assume that J and F are constant in order to compare (2.1)-(2.2) with (4.1)-(4.2). The parameters values are given in Table 3, following [2] and [8], and the simulations are displayed in Figure 4. . Intracellular and capillary lactate dynamics with constant functions J and F . The intracellular lactate trajectory (on the left) and the capillary lactate trajectory (on the right) are both upper bounded and reach the corresponding steady state. The trajectories for the original system (with ε > 0) are also lower bounded. On the right top corner the capillary lactate of the original system has an initial dip, while it does not exist on the capillary lactate curve of the limit system (right bottom corner).

Parameter Value
Unit Table 3. Parameters values.
Note that, in that case, we have shown the existence of an upper bound on u and the presence of a locally stable steady state. The limit system has almost the same dynamics as the initial one. The notable difference is the presence of a hard initial dip for the dynamics of the capillary lactate with ε > 0. This dip does not exist in the dynamics given by the limit system.
Using the parameters values given in Table 3, we perform numerical simulations with various values of ε. The results are given in Figure 5.
While the intracellular lactate trajectories seem not to differ a lot, the capillary lactate trajectories show different initial dynamics. The smaller ε is, the steeper the dip is up to ε = 0, where there is no dip. Remark 7. An approach based on singular perturbation theory has been made on this model by Lahutte-Auboin et al. [7,8]. There, the authors give a geometrical explanation for the initial lactate dip and prove the existence of a periodic solution of the fast-slow system under a repetitive sequence of identical stimuli. In addition, our approach gives an estimate on the rate of convergence with respect to the parameters, which is usually not the case with singular perturbation theory.
Using the parameters given in Table 3, we now test different values for J. The results are given in Figure 6.
There is a limit value of J, denoted by J lim , such that there is an equilibrium only for J < J lim . With these parameters values, we have : There are two types of dynamics : those with J < J lim for which there is a steady state and those with J > J lim for which the intracellular lactate trajectory explodes.
Experimental data. In this section we compare typical results obtained by using the model (with constant J and F ) with in vivo data. This model is known to give good results when fitted to experimental data in small time variations (seconds) [2]. We want here to test its robustness on larger time variations (days).
For each five patients we have four lactate concentration measures (only three for patient 1) separated from each other by more than 80 days. Biological data were collected from patients exhibiting low grade gliomas (WHO grade 2) histologically proven, using monovoxel proton MR spectroscopy sequences performed on a same whole body 3 Tesla magnet (Verio, Siemens Ag) using specific 32 channels head coil. Raw data have been performed under the JMRUI software for appropriate quantification of metabolites, especially lactate concentration. Figure 6. Dynamics for different values for J. On the right, the intracellular lactate trajectories are divided into two groups : for J ∈ {1, 0.1, 0.01}, the concentration seems to explode, while for J ∈ {0.001, 0.0001}, it seems more stable. On the right the capillary lactate trajectories are devided into these two groups. For the first one, we can see a dip, while, for the second one, the steady state is not quickly reached.
We are unable to distinguish between capillary lactate and intracellular lactate using imaging data. Therefore lactate concentration measures are the sum of the capillary lactate and intracellular lactate concentrations.
Because lactate concentrations variations are intrinsic and depend on lactate exchanges, we assume that it is relevant to adjust the initial values (ū 0 andv 0 ) and the exchanges with neighboring cells (J) only. We give the other parameters values in Table 4. Table 4. Parameters values.
The results are given in Figure 7. Fitted values ofū 0 ,v 0 and J for each patient are given in Table 5.  Table 5. Fitted values ofū 0 ,v 0 and J.
Our model (with constant J and F ) is consistent with biological data. it is able to predict what happens after a lactate spike, when the brain sets up a huge regulation to turn back to an acceptable lactate concentration. 6. Discussion. In this study we analyze a model for lactate kinetics first given in [2]. This model is a first step in view of a better understanding of lactate dynamics in the brain. Lactate has a key role in neuroenergetics. Therefore studying its dynamics in the brain is necessary to better understand the energetic breakdown which is observed, for example in tumors. This model is known to give good results when fitted with experimental data [2], [8]. However, to the best of our knowledge, no mathematical analysis has been made to show the existence and uniqueness of the solution and conditions for bounds on the solution, but also comparisons between the original and limit systems.
In this paper, we study the original and limit systems and obtain existence, uniqueness and bounds on the solutions for the two systems. We also give a condition on J which ensures the existence of a locally stable steady state. Moreover, we give an upper estimate on the difference between the solutions of the two systems. We also give several numerical simulations, for different values of ε and J. Finally we confront the model with in vivo data.
When confronted with imaging data from NMR spectroscopy and perfusion, the model provides good results. Because there are large time variations, we cannot ensure that all the parameters remain constant in time. Therefore constant parameters are not good for explaining hard changes on the lactate concentration dynamics such as shifting from a decreasing concentration to an increasing one. Despite this, they can explain what happens after a lactate spike.
Differences in lactate dynamics suggest that there are several glioma profiles with different typical kinetics. This could indicate that non-agressive low grade gliomas show an increasing lactate concentration and then move to more agressive form (WHO II+ to WHO III+). At this stage the glioma will exhibit angiogenesis, modified proteins and altered transporters, which leads to fluctuating lactate kinetics [6]. This point should be considered as a critical one because of its therapeutical management consequences. Yet the patient should be referred to more agressive therapeutics such as radiotherapy or chemotherapy. Also the imaging control frequency should be restrained as well.
As already mentioned in [10], an initial dip exists in the brain lactate dynamics. This dip could be mathematically explained by the different initial values of the capillary lactate concentration between the original and limit systems. Therefore the dip can be biologically explained by compartment volume modifications.
It cannot be excluded that J stands only because the model is build to explain a local brain part dynamics. Then lactate exchanges between cells and the extracellular space have to be better studied to find a biologically relevant function J or other models with a different way to take these interactions into account. Despite this, the function J is described here with minimal assumptions. Therefore, this allows us to use this analysis with further choices for J or modified systems.
Indeed one perspective is to build more complex and suitable models for brain metabolism. Adding oxygen and glucose dynamics to this model can be the next step in view of a more accurate description of energy dynamics in the brain. It could also be interesting to build a model with different cell types (such as astrocyte and neuron), for a better understanding of the brain fuel substrate fluxes.