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 well-posedness 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 $\varepsilon$ 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_{\varepsilon}$ the intracellular lactate concentration and $v_{\varepsilon}$ the capillary lactate concentration where $\varepsilon$ stands for the volume separating the compartments. They are given in mM. One of the main parameters in the model is $\varepsilon$. 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 $\varepsilon$ 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_{\varepsilon}$ 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 \in \mathbb{R}^+$:

$\begin{align} u_{\varepsilon}'(t) & = J(t, u_{\varepsilon}(t))-T(\frac{u_{\varepsilon}(t)}{k+u_{\varepsilon}(t)}-\frac{v_{\varepsilon}(t)}{k'+v_{\varepsilon}(t)}), \label{u}\end{align}$

(2.1)

$\begin{align} \varepsilon v_{\varepsilon}'(t) & = F(t)(L-v_{\varepsilon}(t))+T(\frac{u_{\varepsilon}(t)}{k+u_{\varepsilon}(t)}-\frac{v_{\varepsilon}(t)}{k'+v_{\varepsilon}(t)}).\label{v} \end{align}$

(2.2)

The initial condition is given by :

$(u_{\varepsilon}(0);v_{\varepsilon}(0)) = (\bar{u}_0;\bar{v}_0) \in \mathbb{R}^+\times \mathbb{R}^+.$

The model is biologically described in [2], to which we refer the interested readers for a better understanding of this process.

---

3. The case $\varepsilon>0$

Well-posedness. Recall that an ODE system $x'(t) = f(t, x(t))$ on $\mathbb{R}^n$, $x = [x_1, \cdots, x_n]$, $f = [f_1, \cdots, f_n]$, is called quasipositive if the condition :

$x \geqslant 0, x_k = 0 \Rightarrow f_k(t, x)\geqslant 0$

is verified for all $k = 1, \cdots, n$. System (2.1)-(2.2) obviously is quasipositive. Hence solutions with nonnegative initial data $(\bar{u}_0;\bar{v}_0)$ remain in $(\mathbb{R}^{+})^{2}$ for all positive times.

Since we have for nonnegative $u_1$ and $u_2$ :

$|{\frac{u_1}{k+u_1}-\frac{u_2}{k+u_2}}| = \frac{k|{u_2-u_1}|}{(k+u_1)(k+u_2)} \leqslant \frac{|{u_2-u_1}|}{k},$

then we can rewrite (2.1)-(2.2), setting

$\begin{align*} &X(t): = (u_{\varepsilon}(t);v_{\varepsilon}(t)), \end{align*}$

to have $\forall t \in \mathbb{R}^+$ :

$X'(t) = H(t, X(t)), ~~~~~~~~~ X(0) = X_0,$

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 $\forall t \in \mathbb{R}^+$.

Bounds on the solution. By means of (2.2), we have $\forall t \in \mathbb{R}^+$,

$v_{\varepsilon}'(t) \leqslant -\frac{F_1 v_{\varepsilon}(t)}{\varepsilon}+\frac{F_2L}{\varepsilon}+\frac{T}{\varepsilon},$

which implies, using Gronwall's lemma, that

$v_{\varepsilon}(t) \leqslant \exp (\frac{-F_1 t}{\varepsilon})\bar{v}_0+\int_0^t \exp (\frac{-F_1 (t-s)}{\varepsilon})\frac{T+F_2L}{\varepsilon}ds,$

or equivalently,

$v_{\varepsilon}(t) \leqslant \exp (\frac{-F_1 t}{\varepsilon})\bar{v}_0+\frac{T+F_2L}{F_1}(1- \exp (\frac{-F_1 t}{\varepsilon})).$

One can see, using the above formula, that we have $\forall t \in \mathbb{R}^+$,

$v_{\varepsilon}(t) \leqslant \max(\bar{v}_0, \frac{T+F_2L}{F_1}): = B_v. \label{Bv}$

Theorem 3.1. We can exhibit a sufficient, but not necessary, condition to ensure a bound on $u_{\varepsilon}$. When it exists, let $B_J \in \mathbb{R}$ be such that :

$\forall (t, x) \in \mathbb{R}^2, \ J(t, x) \leqslant B_J$

and :

$\begin{align*} B_J &< T(1 - \frac{B_v}{k'+B_v}) \Leftrightarrow B_J(k'+B_v) < Tk'. &\mathit{(}{C_{{\rm{2}}.{\rm{1}}}}\mathit{)} \end{align*}$

In that case, we have, setting $z = \frac{B_v}{k'+B_v}+\frac{B_J}{T}$ and $\forall t \in \mathbb{R}^+$ :

$\begin{align} u_{\varepsilon}(t) \leqslant \max(\frac{kz}{1-z}, \bar{u}_0): = B_u. \end{align}$

(3.1)

Remark 1. Condition $(C_{2.1})$ is related to the equation $f(x) = 0$ with $f(x) = B_J-\frac{Tx}{k+x}+\frac{TB_v}{k'+B_v}$ for which a positive solution exists if and only if $B_J <T(1-\frac{B_v}{k'+B_v})$, 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.

Proof. Equation (2.1) gives $\forall t \in \mathbb{R}^+$ :

$u_{\varepsilon}'(t) \leqslant B_J+ T\frac{B_v}{B_v+k'}-T\frac{u_{\varepsilon}(t) }{k+u_{\varepsilon}(t) }.$

We set $z = \frac{B_v}{k'+B_v}+\frac{B_J}{T}$ and have, thanks to (C$_{2.1}$) :

$1-z = \frac{(k'+B_v)T}{(k'+B_v)T}-\frac{B_v(T+B_J)+k'B_J}{(k'+B_v)T} = \frac{k'T-B_J(k'+B_v)}{(k'+B_v)T} > 0.$

Let $t \in \mathbb{R}^+$ be such that:

$u_{\varepsilon}(t) > \frac{kz}{1-z}.$

Then,

$u_{\varepsilon}(t)(1-\frac{B_v}{k'+B_v}-\frac{B_J}{T}) > k(\frac{B_v}{k'+B_v}+\frac{B_J}{T}),$

which yields

$B_J+ T\frac{B_v}{B_v+k'}-T\frac{u_{\varepsilon}(t)}{k+u_{\varepsilon}(t)} < 0,$

hence

$u_{\varepsilon}'(t) < 0.$

We finally deduce that

$\begin{align*} u_{\varepsilon}(t) \leqslant \max(\frac{kz}{1-z}, \bar{u}_0). \end{align*}$

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_{\varepsilon}$ is large, we can exhibit a better upper bound satisfying (C$_{2.1}$). It is thus sufficient to write that $u_{\varepsilon}$ is bounded from above with a mild (C$_{2.1}$) condition. For example, we can take $J = J_{test}$ such that, $\forall (x, s) \in \mathbb{R}^+\times\mathbb{R}^+$ :

$J_{test}(s, x) = \underbrace{G_J}_{\text{creation}}-\underbrace{L_J}_{\text{ consumption}}+\underbrace{\frac{C_J}{\varepsilon_J+x}}_{\text{import}},$

for positive constants $G_J$, $L_J$, $C_j$ and $\varepsilon_J$ such that $G_J>L_J$ and $G_J <L_J+\frac{Tk'}{k'+B_v}$. 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 \geqslant \frac{C_j}{\frac{Tk'}{k'+B_v}-G_J+L_J} = N_{test}$, $J_{test}$ is bounded by $\frac{Tk'}{k'+B_v}$ and satisfies (C$_{2.1}$). We conclude that, setting $z = \frac{B_v}{k'+B_v}+\frac{B_J}{T}$, then $u(t) \leqslant \max(N_{test}, \frac{kz}{1-z}, \bar{u}_0)$, $\forall t \in \mathbb{R}^+$. 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_{\varepsilon}$ and $v_{\varepsilon}$ are nonnegative functions. We can exhibit lower bounds by using the same method. Indeed, (2.2) gives, $\forall t \in \mathbb{R}^+$ :

$v_{\varepsilon}'(t) \geqslant -\frac{F_2 v_{\varepsilon}(t)}{\varepsilon}+\frac{F_1L}{\varepsilon}-\frac{T}{\varepsilon}\frac{B_v}{k'+B_v}.$

Then, if $\frac{F_1L-T\frac{B_v}{k'+B_v}}{F_2} \geqslant 0$ and for $t \in \mathbb{R}^+$ such that

$v_{\varepsilon}(t) \leqslant \frac{F_1L-T\frac{B_v}{k'+B_v}}{F_2},$

we have $v_{\varepsilon}(t)'\geqslant 0$, so that, $\forall t \in \mathbb{R}^+$ :

$\begin{align*} v_{\varepsilon}(t) &\geqslant \min(\bar{v}_0, \frac{F_1L-T\frac{B_v}{k'+B_v}}{F_2}). \end{align*}$

If $\frac{F_1L-T\frac{B_v}{k'+B_v}}{F_2} \leqslant 0$, we cannot find a positive lower bound on $v_{\varepsilon}$ and we keep $v_{\varepsilon} \geqslant 0$. Finally :

$\begin{align} v_{\varepsilon}(t) &\geqslant \min(\bar{v}_0, \max(\frac{F_1L-T\frac{B_v}{k'+B_v}}{F_2}, 0)): = M_v. \label{Mv} \end{align}$

(3.2)

Similarly, $\forall t \in \mathbb{R}^+$ :

$\begin{align*} u_{\varepsilon}'(t)&\geqslant T(\frac{M_v}{k'+M_v}-\frac{u_{\varepsilon}}{k+u_{\varepsilon}}). \end{align*}$

Then $u_{\varepsilon}(t)\leqslant M_v\frac{k}{k'}$ leads to $u'(t)\geqslant 0$. This shows that

$\begin{align} u_{\varepsilon}(t) &\geqslant \min (\bar{u}_0, M_v\frac{k}{k'}): = M_u. \end{align}$

(3.3)

Remark 3. The upper bound $B_v$ on $v_{\varepsilon}$ can be derived even when the initial data $\bar{v}_0$ and $\bar{u}_0$ depend on $\varepsilon$ and are bounded with respect to this parameter. To do so, we use the same method, adapting the final step so that

$v_{\varepsilon}(t) \leqslant \max( \underset{\varepsilon>0}{\sup{\bar{v}_0}}, \frac{T+F_2L}{F_1}) .$

The lower bounds $M_v$ and $M_u$ can be obtained in the same way for initial data depending on $\varepsilon$.

Stability of the equilibrium. For constant $J$ and $F$, an equilibrium for (2.1)-(2.2) has been found in [4] :

$\begin{align} u_l&: = \frac{k(\frac{J}{T}+\frac{v_l}{k'+v_l})}{1-(\frac{J}{T}+\frac{v_l}{k'+v_l})}, \end{align}$

(3.4)

$\begin{align} v_l&: = L+\frac{J}{F}. \end{align}$

(3.5)

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 :

$\begin{align} \frac{J}{T}+\frac{LF+J}{F(k'+L)+J}&< 1 \Leftrightarrow J^2+JF(L+k')-TFk' <0. \end{align}$

(3.6)

Remark 4. We have already shown that $v_{\varepsilon}(t) \leqslant L+\frac{T}{F} = B_v$, $\forall t \in \mathbb{R}^+$. We verify that the equilibrium $v_l$ is such that $v_l\leqslant B_v$, which requires $J\leqslant T$. In fact, under $(2.7)$, $J > T$ implies:

$J^2+JF(L+k')-TFk'> T^2+TF(L+k')-TFk' = T(T+FL) > 0.$

Therefore, the contraposition leads to:

$J^2+JF(L+k')-TFk'\leqslant 0 \text{ implies } J \leqslant T.$

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 :

$\Delta_J = F^2(L+k')^2+4TFk'>0$

and there is an equilibrium only when $J \in]J_b, J_h[$, where :

$\begin{align*} J_b: = \frac12 (-F(L+k')-\sqrt{\Delta_J)}, \\ J_h: = \frac12 (-F(L+k')+\sqrt{\Delta_J)}. \end{align*}$

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 $\varepsilon=0$

We now study the limit system for $\varepsilon = 0$, $F(t): = F$, $J(t, x): = J$ constant, given $\forall~t \in \mathbb{R}^+$ by :

$\begin{align} u_0'(t) & = J-T(\frac{u_0(t)}{k+u_{0}(t)}-\frac{v_0(t)}{k'+v_0(t)}), \label{ul}\end{align}$

(4.1)

$\begin{align} 0 & = F(L-v_0(t))+T(\frac{u_0(t)}{k+u_0(t)}-\frac{v_0(t)}{k'+v_0(t)}), \label{vl} \end{align}$

(4.2)

together with the initial condition :

$u_0(0) = \bar{u}_0 \in \mathbb{R}^+.$

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 $\varepsilon>0$) with this limit system (with $\varepsilon = 0$).

Preliminaries. The function $v_{0}$ given by (4.2) is defined as long as $v_{0}(t)$ belongs to $I =]-\infty, -k'[~\cup~]-k', +\infty[: = I_1~\cup~I_2$ and is continuous. Taking $v_{0}(0) = \tilde{v}_0 \in I_2$, then $v_{0}(t) \in I_2$ $\forall~t \in \mathbb{R}^+$. It is biologically relevant to take $v_0(0) = \tilde{v}_0$ as the positive root of (4.2) for $u_0(0) = \bar{u}_0$.

We define the function $\varphi_c$, for any constant $c>0$, by :

${\varphi _c}\left\{ \begin{array}{l} ] - c, + \infty [ \to ] - \infty ,T[\\ \;\;\;\;\;\;\;\;\;\;\;\;s \mapsto \frac{{Ts}}{{c + s}}. \end{array} \right.$

It is easy to see that $\varphi_c$ is a monotone increasing function with $\varphi_c(0) = 0$. We also define an inverse function of $\varphi_c$ :

$\varphi _c^{ - 1}\left\{ \begin{array}{l} [0,T[ \to [0, + \infty [\\ \;\;\;\;\;z \mapsto \frac{{cz}}{{T - z}}. \end{array} \right.$

Furthermore, we introduce the function $\psi_c$ defined by :

${{\psi }_{c}}\left\{ \begin{align} & ]-c,+\infty [\to \mathbb{R} \\ & \ \ \ \ \ \ \ \ \ \ \ s\mapsto Fs+{{\varphi }_{c}}(s), \\ \end{align} \right.$

where $\psi_c$ is a bijection from $]-c, +\infty[$ onto $\mathbb{R}$. It can also be a bijection from $\mathbb{R}^+$ onto itself. Its derivative reads :

$\psi_c'(s) = F+\frac{Tc}{(c+s)^2}.$

Employing (4.2), we have :

$\psi_{k'}(v_0(t)) = FL+\varphi_k(u_0(t)).$

We rewrite (4.1)-(4.2) as :

$\begin{align} v_{0}(t)& = \psi_{k'}^{-1}(FL+\varphi_k(u_0(t))): = \Psi(u_{0}(t)), \end{align}$

(4.3)

$\begin{align} u_{0}'(t) & = J-T(\frac{u_0(t)}{k+u_0(t)}-\frac{\Psi(u_{0}(t))}{k'+\Psi(u_{0}(t))}): = G(t, u_{0}(t)), \label{Fu} \end{align}$

(4.4)

and set, for $y\in[0, +\infty[$ :

$\begin{align*} \Psi^{-1}(y)& = \varphi_k^{-1}(\psi_{k'}(y)-FL). \end{align*}$

A priori bounds on the solutions. Thanks to (4.4), we have :

$G(t, 0) = J+\frac{T\psi_{k'}^{-1}(FL)}{k'+\psi_{k'}^{-1}(FL)}\geqslant J,$

and the system is quasipositive : for an initial condition $\bar{u}_0\geqslant 0$, there holds $u_{0}(t)\geqslant 0$, $\forall t \in \mathbb{R}^+$. A fortiori, we have $v_{0}(t) = \Psi(u_{0}(t)) \geqslant 0$.

Using (4.2), we find an upper bound on $v_0$, $\forall t \in \mathbb{R}^+$ :

$\begin{align} v_{0}(t) &\leqslant L+\frac{T}{F}: = B_{v, 0}. \end{align}$

(4.5)

We can also obtain an upper bound on $u_0$ using Theorem 2.1. When it exists, we call it $B_{u, 0}$.

We now rewrite (4.2) as :

$v_{0}(t)^2+v_{0}(t)(k'-L+\frac{T}{F}-z)-k'(L+z) = 0,$

where $z = \frac{Tu_{0}(t)}{(k+u_{0}(t))F}\leqslant \frac{T}{F}$.

Noting that $v_{0}$ is positive, $\forall t \in \mathbb{R}^+$, we have :

$\begin{align} v_{0}(t) = \frac12 \big(z+L-\frac{T}{F}-k' + \sqrt{(\frac{T}{F}+k'-L-z)^2+4k'(L+z)} \big). \end{align}$

(4.6)

Well-posedness. Equation (4.4) gives :

$\begin{align*} u_{0}'(t) & = J-T(\frac{u_{0}(t)}{k+u_{0}(t)}-\frac{\Psi(u_{0}(t))}{k'+\Psi(u_{0}(t))}): = G(t, u_{0}(t)). \end{align*}$

Lemma 1. The function $\Psi$ is Lipschitz continuous in $u_0$, i.e. there exists $K_L>0$ such that for all $u_1, u_2 \in [0, +\infty],$

$|{\Psi(u_1)-\Psi(u_2)}| \leqslant K_L |{u_1-u_2}|.$

Proof. Let $u \in \mathbb{R}^+$. We know that :

$\begin{align*} \Psi'(u)& = \frac{1}{(\Psi^{-1})'(\Psi(u))} \end{align*}$

and :

$\begin{align*} |{(\Psi^{-1})'(\Psi(u))}| & = |{(\varphi_k^{-1})'(\psi_{k'}(\Psi(u))-FL)\psi_{k'}(\Psi(u)))}|\\ & = |{\frac{Tk}{(T-\psi_{k'}(\Psi(u))+FL)^2}(F+\frac{Tk}{(k+\Psi(u))^2})}|\\ & = |{\frac{Tk(Tk+F(k+\Psi(u))^2)}{(T+\varphi_k(u))^2(k+\Psi(u))^2}}|. \end{align*}$

It follows from the above that $\Psi(u) = v$, with :

$\begin{align*} 0 & = F(L-v)+T(\frac{u}{k+u}-\frac{v}{k'+v}) \Rightarrow v \leqslant FL+T = B_{v, 0}, \\ v &\geqslant 0. \end{align*}$

Therefore, $\Psi(u) = v \in [0, B_{v, 0}]$ and :

$\begin{align*} |{\Psi'(u)}|& = \frac{1}{|{(\Psi^{-1})'(\Psi(u))}|}\\ & = |{\frac{(T+\varphi_k(u))^2(k+\Psi(u))^2}{Tk(Tk+F(k+\Psi(u))^2)}}|\\ &\leqslant (T+\frac{Tu}{k+u})^2(k+\Psi(u))^2\\ &\leqslant 4T^2(k+B_{v, 0})^2 : = K_L. \end{align*}$

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), $\forall t \in \mathbb{R}^+$. Finally we have existence and uniqueness for $v_0$, recalling that $\Psi$ 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 :

$\lambda: = -T \frac{k}{(k+u_l)^2} <0.$

Therefore, $u_l$ is locally stable. Moreover, $\forall t \in \mathbb{R}^+$, $v(t) = \Psi(u(t))$, where $\Psi$ is a bijective function. Setting $z_l = \frac{J}{T}+\frac{v_l}{k'+v_l}$, we have :

$\begin{align*} T\frac{u_l}{k+u_l}& = T\frac{kz_l}{1-z_l}\frac{1-z_l}{k}\\ & = Tz_l\\ & = \frac{TJ}{T}+T\frac{v_l}{k'+v_l}\\ & = F(\frac{J}{F}+L-L)+T\frac{v_l}{k'+v_l}\\ & = F(v_l-L)+T\frac{v_l}{k'+v_l}. \end{align*}$

Thus :

$Fv_l+T\frac{v_l}{k'+v_l} = FL+T\frac{u_l}{k+u_l} \Leftrightarrow v_l = \Psi(u_l),$

and the stationnary point $v_l$ is locally stable.

Comparison between the original and the limit systems. We wish to bound the difference between $(u_{\varepsilon};v_{\varepsilon})$ solution to (2.1)-(2.2) with $F(t): = F$ and $J(t, x): = J$ and $(u_{0};v_{0})$ solution to (4.1)-(4.2), $\forall t \in \mathbb{R}^+$. To do so, we choose the same initial condition for $u_{\varepsilon}$ and $u_0$ :

$u_{\varepsilon}(0) = u_0(0) = \bar{u}_0.$

We set $u = u_{\varepsilon}-u_{0}$ and $v = v_{\varepsilon}-v_0$. Using (2.1), (2.2), (4.1) and (4.2), we have, $\forall~t \in \mathbb{R}^+$ :

$\begin{align} u'(t) & = T(\frac{k'v(t)}{(v_{\varepsilon}(t)+k')(v_{0}(t)+k')}-\frac{k u(t)}{(u_{\varepsilon}(t)+k)(u_0(t)+k)})\label{ue}, \end{align}$

(4.7)

$\begin{align} \varepsilon v'(t) & = -Fv(t)+T(\frac{k u(t)}{(u_{\varepsilon}(t)+k)(u_0(t)+k)}-\frac{k'v(t)}{(v_{\varepsilon}(t)+k')(v_0(t)+k')})\label{ve} -\varepsilon v_{0}'(t). \end{align}$

(4.8)

It follows from the above that, $\forall~t \in \mathbb{R}^+$ :

$\begin{align*} {u}_{0}'(t) & = J-T(\frac{u_{0}(t)}{k+u_{0}(t)}-\frac{v_{0}(t)}{k'+v_{0}(t)}) \in [J-T, J+T], \\ v_{0}(t) &\leqslant B_{v, 0}, \end{align*}$

Therefore, differentiating (4.2), we find :

$\begin{align*} F{v}_{0}'(t) = T(\frac{k{u}_{0}'(t)}{(k+u_{0}(t))^2}-\frac{k'{v}_{0}'(t)}{(k'+v_{0}(t))^2}) & \Rightarrow {v}_{0}'(t) (F+\frac{k'T}{(k'+v_{0}(t))^2}) = \frac{Tk {u}_{0}'(t)}{(k+u_{0}(t))^2}, \end{align*}$

hence $\forall~t \in \mathbb{R}^+$ :

$|{v_{0}'(t)}| \leqslant \frac{kT(J+T)}{(F+\frac{k'T}{(k'+B_{v, 0})^2})}: = \gamma.$

Next, multiplying $(4.7)$ by $u(t)$ and $(4.8)$ by $v(t)$ gives, $\forall~t \in \mathbb{R}^+$ :

$\begin{align} \frac12 \frac{d}{dt}(u^2(t)) &\leqslant \frac{T}{k'}|{u(t)}||{v(t)}|\label{sm1}, \end{align}$

(4.9)

$\begin{align}\varepsilon \frac12 \frac{d}{dt}(v^2(t))+Fv^2(t) &\leqslant \frac{T}{k}|{u(t)}||{v(t)}|+\varepsilon |{v(t)}|\gamma.\label{sm2} \end{align}$

(4.10)

Noting that

$\begin{align*} \frac{T}{k}|{u(t)}||{v(t)}|+\varepsilon |{v(t)}|\gamma & = (\frac{T}{k}|{u(t)}|\frac{2}{\sqrt{F}})(\frac{\sqrt{F}}{2}|{v(t)}|)+(|{v(t)}|\frac{\sqrt{F}}{2})(\frac{2}{\sqrt{F}}\varepsilon\gamma)\\ &\leqslant \frac{F}{2}v^2(t)+\frac{4T^2}{Fk^2}u^2(t)+\frac{4\gamma^2}{F}\varepsilon^2 \end{align*}$

and

$\begin{align*} \frac{T}{k'}|{u(t)}||{v(t)}| & = (\frac{T}{k'}|{u(t)}|\frac{\sqrt{2}}{\sqrt{F}})(\frac{\sqrt{F}}{\sqrt{2}}|{v(t)}|)\\ &\leqslant \frac{F}{2}v^2(t)+\frac{2T^2}{Fk'^2}u^2(t), \end{align*}$

summing (4.9) and (4.10) thus yields, $\forall~t \in \mathbb{R}^+$ :

$\begin{align*} \frac{d}{dt}(u^2(t)+\varepsilon v^2(t)) &\leqslant (\frac{8T^2}{Fk^2}+\frac{4T^2}{Fk'^2})(u^2(t)+\varepsilon v^2(t))+\frac{8\gamma^2}{F}\varepsilon^2. \end{align*}$

Noting finally that :

$u^2(0) = 0 \text{ and }u^2(0)+\varepsilon v^2(0) = \varepsilon (\bar{v}_0-\Psi(\bar{u}_0))^2,$

Gronwall's lemma gives, $\forall~t \in \mathbb{R}^+$ :

$\begin{align*} &u^2(t)+\varepsilon v^2(t) \\&\leqslant \exp \big(\frac{T^2t}{F}(\frac{8}{k^2}+\frac{4}{k'^2}) \big)(\varepsilon (\bar{v}_0-\Psi(\bar{u}_0))^2+\frac{k^2(J+T)^2}{(F+\frac{k'T}{(k'+L+\frac{T}{F})^2})^2}\frac{2\varepsilon^2}{(\frac{2}{k^2}+\frac{1}{k'^2})})\\ &-\frac{k^2(J+T)^2}{(F+\frac{k'T}{(k'+L+\frac{T}{F})^2})^2}\frac{2\varepsilon^2}{(\frac{2}{k^2}+\frac{1}{k'^2})}. \end{align*}$

Remark 5. In the particular case $\bar{v}_0 = \Psi(\bar{u}_0)$, we have $\forall t \in [0, t_m]$ :

$\begin{align*} u^2(t)+\varepsilon v^2(t) &\leqslant (\exp \big(\frac{T^2t_m}{F}(\frac{8}{k^2}+\frac{4}{k'^2}) \big)-1)\frac{2\gamma^2\varepsilon^2}{T^2(\frac{2}{k^2}+\frac{1}{k'^2})} \end{align*}$

which yields that, on the finite time interval $[0, t_m]$,

$\begin{align} |{u(t)}|&\leqslant C_{t_m}\varepsilon, \;\;\;& |{v(t)}|&\leqslant C_{t_m}\sqrt{\varepsilon}. \end{align}$

(4.11)

Remark 6. Setting $v_{0}(0) = \tilde{v}_1 \in I_1$ the negative root of (4.2) for $u_0(t) = \bar{u}_0$, we define a second solution to (4.1)-(4.2) for every $t \in \mathbb{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 $\varepsilon$. 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 :

$J\left\{ \begin{align} & {{\mathbb{R}}^{+}}\to {{\mathbb{R}}^{+}} \\ & \ \ \ x\mapsto {{G}_{J}}-{{L}_{j}}+\frac{Cj}{x+{{\varepsilon }_{j}}}, \\ \end{align} \right.$

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] :

$F\left\{ \begin{array}{l} {\mathbb{R}^ + } \to {\mathbb{R}^ + }\\ t \mapsto \left\{ \begin{array}{l} {F_0}(1 + {\alpha _f})\;\;\;{\rm{if}}\;\exists N \in\mathbb{N} /(N - 1){t_f} + {t_i}{\rm{ < }}t{\rm{ < }}N{t_f},\\ \;\;\;\;\;\;\;\;\;\;\;\;{F_0}\;\;{\rm{if}}\;\;{\rm{not}} \end{array} \right. \end{array} \right.$

The parameters for these two functions are given in Table 1.

We also consider the parameters given in [2] and [8]. In that case, $(\bar{u}_0;\bar{v}_0) = (1.15;1)$ and the parameters values are given in Table 2.

The solutions $u_{\varepsilon}$ and $v_{\varepsilon}$ remain nonnegative, as shown in Figure 3. Here, we have an upper bound on $v_{\varepsilon}$ and J enjoys (C$_{2.1}$), so that $u_{\varepsilon}$ has an upper bound too. Thus, the intracellular and capillary lactate concentrations match with the mathematical analysis. At the beginning, there is a dip for the capillary lactate concentration. The fluctuations on $F$ are rapidly damped out as time grows.

Numerical simulations with different values of $J$ and $\varepsilon$. 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.

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 $\varepsilon>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 $\varepsilon$. 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 $\varepsilon$ is, the steeper the dip is up to $\varepsilon = 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 :

$J_{lim} = 0.00851 \text{ mM.s$^{-1}$}.$

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.

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 ($\bar{u}_0$ and $\bar{v}_0$) and the exchanges with neighboring cells ($J$) only. We give the other parameters values in Table 4.

The results are given in Figure 7. Fitted values of $\bar{u}_0$, $\bar{v}_0$ and $J$ for each patient are given in Table 5.

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 $\varepsilon$ 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 Ⅱ+ to WHO Ⅲ+). 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.

---

Acknowledgments

The authors wish to thank an anonymous referee for her/his careful reading of the paper and helpful comments.

---