Mathematical modeling of brain metabolites variations in the circadian rhythm

1.   Introduction

Magnetic resonance human imaging has emerged 30 years ago and progressively imposed itself as the reference for in vivo noninvasive multiorgans investigation. Recent and huge developments of technical aspects and sequences have lead MRI to switch from "anatomic" to multiparametric powerful investigation system. Namely, multinuclear spectroscopy/imaging and perfusion sequences today provide numerous so called "metabolic" information which match with tremendous increase of human metabolism knowledge (see ^[2,5]).

Then, significant metabolic modifications have been demonstrated in pathologies usually presenting a normal "anatomic" appearance, such as psychiatric diseases (see ^[1,3]). The correlative (and underlying) question concerned the measurements' reliability. Yet circadian variations of metabolite concentrations, as they can be calculated from spectroscopic acquisitions, appeared to be of importance in this context.

We performed three successive acquisitions at "key" moments referring to the circadian rhythm, 7.30 am, 1.30 pm and 5.30 pm, on 30 healthy volunteers and on 12 different areas in the brain. Several fluctuations were found within acquisition time and location, thus suggesting modeling usefulness.

We propose hereafter to develop a mathematical analysis and modeling for describing and integrating these fluctuations. More precisely, we consider models of the form

with a symport nonlinear term (see ^[4]). Here, $u$ is a metabolite's concentration, e.g., lactate. We give a mathematical analysis of these models; the main difficulty is the (expected) nonnegativity of the concentration. We then present numerical simulations which we compare with the medical data mentioned above.

2.   Mathematical analysis

We consider the initial value problem

Noting that the function

is $\mathcal C^1$ on $\mathbb{R} \times (-k', +\infty)$, it follows from the Cauchy–Lipschitz theorem that (2.1)-(2.2) possesses the maximal solution

We have the

Theorem 2.1. The maximal solution $y$ is defined on $\mathbb{R} ^+$ (i.e., $T_\star = +\infty)$.

Proof. Let $y_+$ be the maximal solution to the initial value problem

Noting that the function

is of class $\mathcal C^1$ and globally Lipschitz continuous on $\mathbb{R}$, it follows that $y_+$ is uniquely defined on $\mathbb{R} ^+$. Multiplying (2.3) by $-y^-_+$, where $y_+ = y^+_+-y^-_+$, $y_+^+ = \max (y_+, 0)$, $y_+^- = \max (-y_+, 0)$, we obtain

so that

Therefore

and

This yields that $y_+$ is the unique solution to

defined on $\mathbb{R} ^+$.

Let then $y_-$ be the maximal solution to

defined on $[0, T_\star ^-)$ (the existence and uniqueness of $y_-$ follows from the Cauchy–Lipschitz theorem). We can note that (2.8) possesses a unique equilibrium defined by

i.e.,

Setting $z = y_–y_-^e$, we have

where

with $h(x) = {{kx}\over {k'+x}}$. Therefore,

and

This yields that $y_-$ is bounded, so that $T_\star ^- = +\infty$ and

Now, it follows from the comparison principle that, $\forall t\in [0, T_\star)$,

In particular, if $T_\star < +\infty$, then $y$ remains bounded as $t\to T_\star$, which yields a contradiction. Thus, $T_\star = +\infty$, which finishes the proof.

Remark 2.2. We can actually prove, proceeding as above, that $y(t)\ge y_-^e$, $t\ge 0$. Indeed, we have, setting $z = y-y_-^e$,

Multiplying the above equality by $-z^-$, we obtain

which yields, noting that $g'\ge 0$ and $-zz^- = \vert z^-\vert ^2$,

and the result follows. Note however that (2.10) gives a better lower bound on $y$, in the sense that it yields that $y$ is nonnegative at least for some time. Note indeed that $y$ corresponds to a concentration and is expected to be nonnegative (see below).

We then have the

Theorem 2.3. We assume that $k > a$. Then the solution $y$ to (2.1)-(2.2) (defined on $\mathbb{R} ^+$) is bounded.

Proof. We again consider the solution $y_+$ to (2.6)-(2.7), defined on $\mathbb{R} ^+$. Note that, when $k > a$, (2.6) possesses a unique equilibrium $y_+^e$ defined by

i.e.,

Proceeding as in the proof of Theorem 2.1, we find

where

Therefore, $y_+$ is bounded and, recalling (2.12), it follows that $y$ is bounded, which finishes the proof of the theorem.

Actually, assuming that $y_0$ is large enough, we can also prove that $y$ is bounded on $\mathbb{R} ^+$. More precisely, we have the

Theorem 2.4. We assume that

Then, the solution $y$ to (2.1)-(2.2) is bounded.

Proof.

Let now $y_+$ be the solution (defined on $\mathbb{R} ^+$) to

Clearly,

so that

Next, we set $z = y-y_+$ and have

which we can rewrite as

Multiplying (2.20) by $z^+$, we obtain

which yields, noting/recalling that $k'+y > 0$ and $y_+\ge 0$,

so that

Employing once more the comparison principle, we deduce that

where $y_-$ is as defined in the proof of Theorem 2.1. The proof follows from (2.21).

Remark 2.5. One crucial question is whether the solution $y$ to (2.1)-(2.2) remains nonnegative, assuming that $y_0 > 0$ (note indeed that, if, for some $t\ge 0$, $y(t) = 0$, then

and can be negative; in particular, if $y_0 = 0$ and $\sin c < 0$, then $y'(0) < 0$); as already mentioned, $y$ corresponds to a concentration and is thus expected to be nonnegative. We thus see that $y$ can be negative when $y_0 = 0$ (and even when $y_0 > 0$ is small). We believe however that, when $y_0 > 0$ is large enough, then $y$ remains positive; this is confirmed by the numerical simulations below. Unfortunately, we have not been able to prove this. One difficulty is to find a nonnegative subsolution.

Remark 2.6. A model, ensuring nonnegativity, reads

As above, we have the existence and uniqueness of the solution, defined on $\mathbb{R} ^+$, such that $y(t) > -k'$, $t\ge 0$. Then, multiplying (2.22) by $-y^-$, we obtain

so that

and

The numerical simulations below show that this model gives satisfactory results for the lactate concentrations (the solution to the first model becomes indeed negative in that case). We can also note that, as far as the other metabolites concentrations are concerned, the first model matches better the experimental data. However, it is less precise for predicting the longtime behavior of the concentrations.

3.   Numerical simulations

3.1.   Using model (2.1)

We test model (2.1) with medical data, supplied from MRI scans by the Poitiers University Hospital (see the introduction). These data consist, for several patients, in the concentrations of brain metabolites (e.g., lactate, choline) at different brain areas and at different times of the day: 7.30 am, 1.30 pm and 5.30 pm. Basically, these data can be divided into four classes: the hat shape (the second dosage is maximal), the reverse hat shape (the second dosage is minimal), the increasing shape (the three dosages increase with respect to time) and the decreasing shape.

In Figure 1, we present, for each class of data, a comparison between the medical data and the simulations, during the ten hours (left) and during the three days (right) next to the first MRI scan. In the right columns of Figure 1, the blue circles correspond to actual measured data, and the red crosses correspond to predicted values, taking into account the circadian rhythm. The first row of the figure corresponds to a choline dosage, for a patient of the hat shape class. The corresponding model is $y'(t) = - \frac{0.007 y}{10.5 +y} + 0.5 \sin (0.265 t + 1.4)$. The second row corresponds to a choline dosage, for a patient of the reverse hat shape class, with the model $y'(t) = -\frac{0.01 y}{5 + y} + 0.5 \sin (0.265 t - 0.85)$. The third row corresponds to a choline dosage, for a patient of the increasing shape class, with the model $y'(t) = -\frac{0.004y}{0.8 +y } + 0.29 \sin (0.265 t + 1)$. Finally, the fourth row corresponds to a Cr dosage for a patient of the decreasing shape class, and the model is $y'(t) = - \frac{0.004y}{7+y} + 0.3 \sin (0.265 t + 2.8)$. Figure 1 shows that the model matches well the medical data.

3.2.   Using model (2.22)

Although model (2.1) seems convenient in many cases, it is not suitable for modeling the variations of lactate. Indeed, the concentrations being small for this metabolite, the simulations may become negative during the simulations. In Figure 2, we compare, in the case of a lactate concentration, two solutions given by the models (2.1) ($y'(t) = - \frac{0.11 y}{8 +y} + 0.5 \sin (0.265 t + 1.9)$) and (2.22) ($y'(t) = - \frac{2 y}{5.5 +y} + 0.9 \sin^2 (0.265 t + 0.9)$). With the first model, the solutions become negative after ten hours, while, with the second one, the solutions remain positive (see Remark (2.6)).

In Figure 3, we compare, for a choline concentration, two solutions : one given by model (2.1) ($y'(t) = - \frac{0.007 y}{10.5 +y} + 0.5 \sin (0.265 t + 1.4)$) and the second one given by model (2.22) ($y'(t) = - \frac{0.7 y}{10.7 +y} + 0.8 \sin^2 (0.1325 t + 1.4))$. As far as the first ten hours are concerned, the solutions are very similar (Figure 2, left). Furthermore, although both solutions match well the medical data, they are very different after the three days period (right). In particular, the range of variations of the solution of model (2.1) is significantly larger than the one for model (2.22) and may be less medically relevant.

Conflict of interest

The authors declare no conflicts of interest in this paper.