Dmitry V Sorokin

Igor Peterlík

Vladimír Ulman

David Svoboda

Martin Maška

MODEL-BASED GENERATION OF SYNTHETIC 3D TIME-LAPSE SEQUENCES OF MOTILE CELLS WITH GROWING FILOPODIA

Keywords: Simulation, 3D time-lapse sequence, synthetic cell, cell deformation, filopodium evolution

The existence of benchmark datasets is essential to objectively evaluate various image analysis methods. Nevertheless, manual annotations of fluorescence microscopy image data are very laborious and not often practicable, especially in the case of 3D+t experiments. In this work, we propose a simulation system capable of generating 3D time-lapse sequences of single motile cells with filopodial protrusions, accompanied by inherently generated ground truth. The system consists of three globally synchronized modules, each responsible for a separate task: the evolution of filopodia on a molecular level, linear elastic deformation of the entire cell with filopodia, and generation of realistic, time-coherent cell texture. The capability of our system is demonstrated by generating a synthetic 3D time-lapse sequence of a single lung cancer cell with two growing filopodia, visually resembling its real counterpart acquired using a confocal fluorescence microscope.

INTRODUCTION

Filopodia are dynamic bundles of parallel actin filaments protruding forward from the periphery of motile cells [START_REF] Carlier | Actin-based Motility[END_REF]. Their typical function is to probe the surrounding environment for cues and promote the locomotion of crawling cells, thus contributing to important physiological processes, such as chemotaxis, angiogenesis, or embryonic development [START_REF] Mattila | Filopodia: Molecular architecture and cellular functions[END_REF].

Modern time-lapse confocal microscopy and fluorescent labeling is a particularly powerful approach to study such processes in vivo [START_REF] Ortiz-De-Solórzano | Toward a morphodynamic model of the cell[END_REF], calling for robust, quantitative, and ideally fully automatic image analysis [START_REF] Cardona | Current challenges in open-source bioimage informatics[END_REF][START_REF] Carpenter | A call for bioimaging software usability[END_REF]. However, an objective assessment of image analysis methods remains an open question. A thorough evaluation requires the availability of reference annotations. Being often prepared manually by human This work was supported by the Czech Science Foundation under the grant number GA16-03909Y. experts, this leads to a subjective, highly error-prone, and extremely laborious task, especially when it comes to manually annotating 3D time-lapse image data [START_REF] Coutu | Probing cellular processes by long-term live imaging-historic problems and current solutions[END_REF].

To address the question of objective evaluation, two recent benchmarks of particle and cell tracking methods [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF][START_REF] Maška | A benchmark for comparison of cell tracking algorithms[END_REF] exploited computer-simulated image datasets, accompanied by inherently generated ground truth. The simulation platforms used in these benchmarks are capable of generating realistic time-lapse sequences of intra-cellular and inter-cellular particles [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF] and of motile, nuclearly stained cells [START_REF] Svoboda | Mitogen: A framework for generating 3D synthetic time-lapse sequences of cell populations in fluorescence microscopy[END_REF]. To the best of our knowledge, there are no tools currently available capable of generating 3D time-lapse sequences of motile, cytoplasmically stained cells with filopodial protrusions, although they might help to stimulate the development of filopodium segmentation and tracking routines fully in 3D [START_REF] Tsygankov | CellGeo: A computational platform for the analysis of shape changes in cells with complex geometries[END_REF][START_REF] Barry | Open source software for quantification of cell migration, protrusions, and fluorescence intensities[END_REF].

In this paper, we propose a first simulation system for generating 3D time-lapse sequences of single motile cells with filopodial protrusions. The simulation of filopodium growth is based on a stochastic model implemented using the Gillespie algorithm. The deformation of the entire cell with filopodia is performed using an elastic model. The realistic cell texture is generated by rendering simulated fluorescent particles into the images, submitted afterwards to a virtual microscope. The output generated by the simulation system has the form of a 3D synthetic time-lapse sequence, accompanied by labeled masks of the cell body and individual filopodia. Such a composite output allows one to benchmark filopodium segmentation and tracking algorithms and to study how segmentation and tracking errors can affect quantitative parameters of motile cells with filopodial protrusions, such as the number, length, and lifetime of filopodia or their protrusion and retraction rates.

PROPOSED SIMULATION SYSTEM

In this section, we describe our simulation system, including the system overview, filopodia growth simulation, filopodia and cell deformation, and realistic texture generation.

Algorithm 1

The pseudocode of the simulation process. Inputs: global time step (∆T ), time limit (T lim ), initial parameters of the model (initParams). Outputs: annotated 3D timelapse sequence (frame i for i = 1, ..., T lim /∆T ).

1: [fList 0 , cell 0 , texture 0 ] ← initializeModel(initParams) 2: frame 0 ← renderImg(texture 0 , fList 0 , cell 0 ) 3: for (i = 0; i • ∆T <= T lim ; i + +) do [cell i+1 , fList def ] ← deformFilo(fList i , cell i , ∆T )

6: fList i+1 ← updGeometry(fList new , fList def ) 7: texture i+1 ← updTexture(texture i , fList i+1 , fList i , cell i+1 , cell i , ∆T ) 8:
frame i+1 ← renderImg(texture i+1 , fList i+1 , cell i+1 ) 9: end for

System overview

The simulation of a 3D time-lapse sequence of a single motile cell with filopodial protrusions consists of three independent parts: modeling of filopodium growth, simultaneous deformation of the cell body and individual filopodia, and generation of the realistic cell texture. Each of these three tasks is done by a separate system module: geometry module (GM), deformation module (DM), and texture module (TM). The modules are independent software components and they are executed from the main component, a supervisor, which is responsible for the proper exchange of the data between the modules, setting the parameters, and synchronizing the time.

The whole simulation process (see Alg. 1) starts by initializing the meshes of the cell body (cell 0 ) and individual filopodia (fList 0 ), setting their growth parameters and the positions where they are attached to the cell body, generating the initial cell texture (texture 0 ), and finally rendering the first frame of the time-lapse sequence (frame 0 ). In our data representation, cell i includes only its geometry stored as a tetrahedral mesh. Each filopodium in fList i is represented by its geometry and growth-related data. The filopodium geometry is stored as a 3D polyline, equipped with a scalar data defining the radius in each node. The growth-related data consists of the parameters and molecule concentrations described in Section 2.2. After the initialization, the simulation process continues iteratively. In each iteration, the morphology of each individual filopodium (i.e., its length and width) is evolved (fList new ) first. Next, the cell body (cell i ) and the filopodia (fList i ) from the previous iteration are deformed using the pressures applied to the cell body and forces applied to the filopodia. As a result, we obtain the deformed cell body (cell i+1 ) and the deformed filopodia (fList def ) that are subsequently integrated with their actual morphology (fList new ), resulting in the fully updated filopodia (fList i+1 ). In case a filopodium has grown, the corresponding segment is added to the tip of its mesh. The position of the new segment is obtained by extrapolating the positions of the filopodium nodes in fList def in the nodes corresponding to the new segment, thus resulting in a smooth continuation of the deformed filopodium. In case the filopodium has shortened, we remove the corresponding segments from its mesh. After the geometry of the whole cell is updated, we update the texture accordingly, which boils down to updating the positions of the simulated fluorescent particles. Finally, the current frame of the sequence (frame i+1 ) is created by rendering the model.

Since the modules are independent, they can work with their own internal time steps, regularly exporting the particular outputs. A user-defined global time interval ∆T , during which each module is operating, corresponds to the time step between the frames of the resulting time-lapse sequence. The simulation process is terminated after the total accumulated time reaches a pre-defined limit T lim .

The labeled masks of the cell body and individual filopodia, serving as the ground truth for a generated 3D time-lapse sequence, are formed by projecting and rasterizing the meshes of the cell body and individual filopodia to the image domain.

Modeling of filopodium growth

A filopodium is a bundle of N actin filaments enclosed by the cell membrane [START_REF] Carlier | Actin-based Motility[END_REF]. Its length is defined as the length of the longest actin filament h = max(h n ) where n is the actin filament index. Each actin filament consists of actin monomers connected with a half-monomer-size overlap δ = 2.7 nm [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF] (see Fig. 1). The evolution of each actin filament is affected by four processes: 1) diffusion of free actin monomers along the filopodium; 2) polymerization and depolymerization of actin filaments; 3) fluctuation of the membrane under load, slowing down the polymerization rate of actin filaments; and 4) a constant velocity retrograde flow that continuously pulls actin filaments into the cell body. We modeled these processes on a molecular level using the Gillespie algorithm [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF].

The filopodium volume is discretized by dividing it into L compartments, with a compartment length of l d = 50 nm. Let a l , l = 1, . . . , L, denote the number of actin monomers in the l-th compartment, being indexed from the filopodium tip. The diffusion of each free actin monomer is defined as 1D random walk [START_REF] Erban | Reactive boundary conditions for stochastic simulations of reaction-diffusion processes[END_REF] with reflective boundary conditions at the filopodium tip (P 1→0 = 0). At the filopodium base, the number of actin monomers is kept unchanged at the constant bulk value (a L+1 = 5) as the number of free actin monomers in the cell body is significantly higher than in the filopodium. Since free actin monomers move between neighboring compartments and stay in the same compartment with equal probability, the reaction time of free actin monomer diffusion is given as ∆t =

l 2 d 3D
, where D is the diffusion constant [START_REF] Erban | Reactive boundary conditions for stochastic simulations of reaction-diffusion processes[END_REF]. The actin filament polymerization rate k n depends on the thermal energy of the system and the membrane load force, defined for each individual filament depending on membrane fluctuations [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF]. The actin filament depolymerization rate k d does not depend on the membrane load force and has the same value for all filaments.

Taking into account the above mentioned processes, the motion equation for each filament length h n is defined as:

∆h n = -ν retr ∆t + ξ(a n , k n , ∆t)δ -η(a n , k d , ∆t)δ
where ∆h n is the change of filament length over the reaction time ∆t, ν retr is the retrograde flow velocity, and ξ and η are the binary indicators of polymerization and depolymerization events, respectively.

Every iteration of the Gillespie simulation involves three steps [START_REF] Erban | A practical guide to stochastic simulations of reaction-diffusion processes[END_REF]. First, two random numbers, r 1 and r 2 , uniformly distributed in the [0, 1] interval are generated. Next, the type of event to occur in a particular iteration is established as depolymerization (ξ = 0

; η = 1) if 0 ≤ r 1 < k d k d +kn , polymer- ization (ξ = 1; η = 0) if k d k d +kn ≤ r 1 < kn k d +kn
, or free actin monomer diffusion (ξ = η = 0) if r 1 ≥ kn k d +kn . Finally, depending on the event type, the number of free actin monomers in each compartment is updated, and the reaction time ∆t is calculated as either ∆t = l 2 d /(3D) in case of diffusion or ∆t = (a tip n (k d + k n )) -1 ln(1/r 2 ) in case of polymerization or depolymerization, with a tip n being the number of free actin monomers in the filopodium compartment that corresponds to the tip of the n-th filament, not necessarily being the first one.

The Gillespie simulations are performed individually for each actin filament over the global time step ∆T . The resulting geometry changes are transmitted to the deformation and texture modules of the simulation system.

The 3D filopodium shape is modeled as a Gaussian with σ corresponding to the number of filaments. The radii at the points between compartments (nodes) are set accordingly.

Deformation of cell body and filopodia

To mimic realistic deformations, the cell body and the filopodia are modeled as elastic objects coupled together via mechanical constraints. For the both types of objects, the corotational finite element formulation of the linear elasticity is employed. This approach is computationally efficient and at the same time, it allows for correct simulation of large deformations [START_REF] Felippa | A systematic approach to the element independent corotational dynamics of finite elements[END_REF]. Therefore, the mechanical response of both cell body and filopodia to given applied forces and pressures is determined by Young's modulus and Poisson's ratio, which can be set independently for each object instance. While the same mathematical formulation is used for all simulated objects, different types of elements are employed for their discretization. The cell body is discretized by a volume mesh composed of linear tetrahedra: each element is given by four nodes having three positional degrees of freedom. As in [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF], random pressures can be applied to the surface of this mesh, imposing motion and deformation of the cell body. Since the geometry of filopodia is given by line segments equipped with radii, we opted for linear beam elements based on the Timoshenko formulation [START_REF] Przemieniecki | Matrix structural analysis of substructures[END_REF]. While being computationally efficient, the beam elements allow for modeling of elongation and bending, which are two basic types of deformations needed for realistic simulation of filopodial movement. Unlike tetrahedral elements, the beam is given by two nodes, each having six degrees of freedom: positions and orientations necessary to model the bending. It is therefore possible to apply linear and rotational forces (torques) to the nodes of filopodia.

Each filopodium is mechanically attached to the cell body: it is supposed that initially the base of each filopodium is located inside a tetrahedron that belongs to the cell body mesh. The initial position and orientation of the base with respect to the corresponding tetrahedron is then kept constant during the entire simulation step ∆T similarly as in [START_REF] Peterlík | Modeling and realtime simulation of a vascularized liver tissue[END_REF]. Therefore, the motion of the cell body due to the applied surface pressures results in the motion of the filopodium via the motion of its base. Vice versa, a motion of the filopodium is propagated to the cell body via the tetrahedron associated to the filopodium base. Numerically, the coupling between the objects is implemented using bilateral constraints based on Lagrange multipliers: in every internal step of the deformation module, first the dynamics of each object (i. e., every filopodium and the cell body) is computed independently. Then, this predictive motion is corrected by applying the constraints that maintain the initial positions and orientations of filopodium bases with respect to the cell body tetrahedra. The implicit Euler solver is used for time integration, which requires solution of the local system of linear equations corresponding to each object. For the cell body, the solution is computed by a solver based on the LDL-decomposition, whereas a solver based on the Thomas algorithm for block tridiagonal matrices is employed for each filopodium [START_REF] Cotin | New approaches to catheter navigation for interventional radiology simulation[END_REF].

Generation of cell texture

The cell texture is represented as a set of many simulated fluorescent particles placed into a phantom image of a simulated cell. In particular, a floating-point precision position is stored with every particle to represent its real world coordinate, regardless of a desired resolution of created images. When rendering successive frames of the generated sequence, the particle positions follow the cell geometry changes by tracking and interpolating displacements of the mesh vertices. The intensity assigned to a phantom image voxel is proportional to the number of particles that are projected into the voxel volume. The initial distribution of particles is sampled from an image with an isotropic mosaic of Perlin noises, and is restricted to the cell volume. Every voxel introduces multiple particles at the same respective coordinate, the number of them is proportional to the voxel intensity. For each filopodium, the particles are placed on its skeleton in amounts proportional to the local radius of the filopodium.

To obtain the final image that is visually similar to a real microscopy one, the phantom image is submitted to a virtual microscope [START_REF] Svoboda | Generation of digital phantoms of cell nuclei and simulation of image formation in 3D image cytometry[END_REF]. In particular, the phantom image is blurred by an experimental point spread function of a real optical system, affected by an uneven illumination phenomenon, resampled to the resolution as that of the real microscope setup, and degraded by shot noise, dark current, and read-out noise.

RESULTS

To demonstrate the capability of our simulation system, we focused on a motile lung cancer cell with two filopodial protrusions. As a reference, we considered real image data with the voxel size of 0.126 × 0.126 × 1.0 µm, being acquired using a spinning disk confocal microscope equipped with a water Plan-Apo 63×/1.20 objective lens [START_REF] Maška | Automatic quantification of filopodia-based cell migration[END_REF]. The initial geometry of the simulated cell body was taken from the real image data. The bases of both filopodia were randomly placed on the cell body surface. For both filopodia, their growth parameters were taken from [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF] and their growth dynamics were regulated by random retrograde flow velocities of 45 to 55 nm/s. The cell body was deformed by random pressures generated similar to [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF]. The filopodia were deformed by applying random forces of 0.5 to 1 mN to their tips. The Young's modulus for the cell body and filopodia was set to 1 kPa and 50 kPa, respectively, reflecting the fact that filopodia are stiffer than the cell body [START_REF] Carlier | Actin-based Motility[END_REF]. The Poisson's ratio for both the cell body and filopodia was fixed at 0.4. The global time step and time limit were fixed at ∆T = 2 s and T lim = 100 s, respectively, yielding a time-lapse sequence of 51 frames long. The synthetic intensity profiles were adapted to visually imitate the real ones by scaling and shifting constants. Furthermore, factoring the cell body volume by utilizing a mosaic of Perlin noises during the cell texture initialization allowed us to introduce the more pronounced cell structures found near the outer cell membrane, and the membrane itself. Approximately nine million particles for the cell body and a quarter of million for both filopodia contributed to each phantom image.

The resulting 3D time-lapse sequence generated with the described parameters is shown in Fig. 2, demonstrating the growth and motion of filopodia. To compare visually the synthetic cell texture with its real counterpart, an example of XY and YZ slices of both is shown in Fig. 3.

CONCLUSION

We have presented a novel approach to generating synthetic 3D time-lapse sequences of single motile cells with filopodial protrusions. The proposed approach utilizes a stochastic model for the simulation of filopodium growth on a molecular level. The motion and deformation of the cell is modeled using the linear elasticity. The cell texture is composed of a large number of simulated fluorescent particles, initially sampled from a mosaic of Perlin noises and updated according to the cell geometry changes. As a proof of the concept, we generated a 3D time-lapse sequence of a single lung cancer cell with two growing filopodia, showing coherent physic-laws constrained shape changes followed by coherently developing texture. In future work, we intend to adapt the simulation parameters across different cell phenotypes, including those with filopodium branching and interaction with one another.

4 :

 4 fList new ← evolveFilo(fList i , ∆T ) 5:

Fig. 1 .

 1 Fig. 1. Filopodium molecular growth model (l d is the compartment size, h n is the filament length).

Fig. 2 .

 2 Fig. 2. Two maximum intensity projections of the resulting 3D time-lapse sequence. The voxel size is 0.126 × 0.126 × 1.0 µm.

Fig. 3 .

 3 Fig. 3. An example of XY and YZ slices of real (R) and synthetic (S) images. The voxel size is 0.126 × 0.126 × 1.0 µm.


Alain Miranville

MATHEMATICAL ANALYSIS OF A PARABOLIC-ELLIPTIC MODEL FOR BRAIN LACTATE KINETICS

Keywords: 2010 Mathematics Subject Classification. 35M33, 35K57, 35K67, 35B45 Brain lactate kinetics, parabolic-elliptic system, reaction-diffusion system, limit system, nonnegative solutions, well-posedness, error estimates, linear stability. 1

Our aim in this paper is to study properties of a parabolic-elliptic system related with brain lactate kinetics. These equations are obtained from a reaction-diffusion system, when a small parameter vanishes. In particular, we prove the existence and uniqueness of nonnegative solutions and obtain error estimates on the difference of the solutions to the initial reaction-diffusion system and those to the limit one, on bounded time intervals. We also study the linear stability of the unique spatially homogeneous equilibrium.

Introduction

The following system of ODE's:

(1.1)

du dt + κ( u k + u - v k + v ) = J, κ, k, k , J > 0, (1.2) ε dv dt + F v + κ( v k + v - u k + u ) = F L, ε, F, L > 0,
where ε is a small parameter, was proposed and studied as a model for brain lactate kinetics (see [START_REF] Costalat | Mathematical modeling of metabolism and hemodynamics[END_REF], [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF], [START_REF] Lahutte-Auboin | Dip and buffering in a fast-slow system associated to brain lactate kinetics[END_REF] and [START_REF] Lahutte-Auboin | On a minimal model for hemodynamics and metabolism of lactate: application to low grade glioma and therapeutic strategies[END_REF]; see also [START_REF] Aubert | Interaction between astrocytes and neurons studied using a mathematical model of compartmentalized energy metabolism[END_REF]). In this context, u = u(t) and v = v(t) correspond to the lactate concentrations in an interstitial (i.e., extra-cellular) domain and in a capillary domain, respectively. Furthermore, the nonlinear term κ( u k+u -v k +v ) stands for a co-transport through the brain-blood boundary (see [START_REF] Keener | Mathematical physiology[END_REF]). Finally, J and F are forcing and input terms, respectively, assumed frozen (more generally, J depends on t and u and accounts for the interactions with a third intracellular compartment (which includes both neurons and astrocytes), while F = F (t) (an applied electrical stimulus; see [START_REF] Hu | A temporary local energy pool coupled to neuronal activity: fluctuations of extracellular lactate levels in rat brain monitored with rapid-response enzyme-based sensor[END_REF]) is piecewise linear and periodic). This model has essential applications to the therapeutic management of glioma (also called glial tumors); see [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF] for thorough discussions on this issue.

Let us assume that u(0) and v(0) are nonnegative (recall that u and v are concentrations and are thus expected to be nonnegative). Then, noting that, if u(0) = 0, then du dt (0) > 0 and, if v(0) = 0, then dv dt (0) > 0, it follows from Cauchy-Lipschitz theorem that, for t > 0 small, u and v exist and are nonnegative. This also yields that the solutions are defined and remain nonnegative on the whole interval R + ; indeed, it is not difficult to prove that they are bounded on finite time intervals. Furthermore, in [START_REF] Costalat | Mathematical modeling of metabolism and hemodynamics[END_REF], [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF], [START_REF] Lahutte-Auboin | Dip and buffering in a fast-slow system associated to brain lactate kinetics[END_REF] and [START_REF] Lahutte-Auboin | On a minimal model for hemodynamics and metabolism of lactate: application to low grade glioma and therapeutic strategies[END_REF], questions related to the stability of the unique equilibrium were addressed. This constitutes an essential point in the modeling, since, as discussed in [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF], a therapeutic perspective of such a result is to have the steady state outside the viability domain, where cell necrosis occurs. Finally, in [START_REF] Lahutte-Auboin | Dip and buffering in a fast-slow system associated to brain lactate kinetics[END_REF], justifications for the dip and buffering which are observed in experiments (see [START_REF] Hu | A temporary local energy pool coupled to neuronal activity: fluctuations of extracellular lactate levels in rat brain monitored with rapid-response enzyme-based sensor[END_REF]) were given, based on geometrical arguments and averaging theory on a slow manifold.

We can note that the above ODE's model does not account for spatial diffusion. Taking this into account would be relevant and desirable from a biological point of view. The simplest possible corresponding PDE's (reaction-diffusion) system, accounting for spatial diffusion, reads

(1.3) ∂u ∂t -α∆u + κ( u k + u - v k + v ) = J, α > 0, (1.4) ε ∂v ∂t -β∆v + F v + κ( v k + v - u k + u ) = F L, β > 0,
where u = u(x, t) and v = v(x, t), which we consider in a bounded and regular domain Ω of R N , N = 1, 2 or 3, together with Neumann boundary conditions,

∂u ∂ν = ∂v ∂ν = 0 on Γ,
where Γ = ∂Ω and ν is the unit outer normal vector. Note that the terms -α∆u and -β∆v correspond to random motions. Note however that more precise models should account for the geometry, i.e., the different compartments (interstitial, capillary), so that (1.3)-(1.4) should be viewed as a very first step towards PDE's models for brain lactate kinetics. We will consider more realistic models elsewhere. We studied in [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF] the existence, uniqueness and regularity of nonnegative solutions to (1.3)-(1.4) (note that the mathematical analysis of (1.3)-(1.4) (and, in particular, the wellposedness) appears to be challenging, due to the coupling terms, especially for negative initial data (though biologically irrelevant, this makes sense from a mathematical point of view); this is also the case for the ODE's model (1.1)-(1.2)). We further established the linear (exponential) stability of the unique spatially homogeneous equilibrium. We also mention [START_REF] Miranville | A singular reaction-diffusion equation associated with brain lactate kinetics[END_REF] in which we proved the existence, uniqueness and regularity of the solutions to the following singular reaction-diffusion equation:

(1.5) ∂u ∂t -∆u + F u + κ u k + u = f (x, t), F ≥ 0,
corresponding to the case where either u or v is known in (1.3) and (1.4); we can also think of (1.5) as an equation in each compartment, assuming that the lactate concentration is known in the other one. Our aim in this paper is to study the limit system

(1.6) ∂u ∂t -α∆u + κ( u k + u - v k + v ) = J, (1.7) -β∆v + F v + κ( v k + v - u k + u ) = F L,
corresponding to ε = 0 in (1.4). We prove the existence and uniqueness of nonnegative solutions to (1.6)-(1.7). We then prove that the solutions to the initial reaction-diffusion system converge to those to the limit parabolic-elliptic one, on finite time intervals, and provide an error estimate in terms of ε. We finally study the linear stability of the unique spatially homogeneous equilibrium. We can note that a similar analysis would also be relevant in the context of the ODE's model (1.1)-(1.2). Though some of the results obtained here could apply (in a simpler way) to this system, this will be considered in more details elsewhere.

Notation: We denote by ((•, •)) the usual L 2 -scalar product, with associated norm • . More generally, • X denotes the norm on the Banach space X and, if X is a Hilbert space, ((•, •)) X denotes the associated scalar product. Throughout the paper, the same letters c, c and c denote positive constants which may vary from line to line. Similarly, the same letter Q denotes continuous and monotone increasing (with respect to each argument) functions which may vary from line to line.

The case ε > 0

We consider the following initial and boundary value problem:

(2.1) ∂u ∂t -α∆u + κ( u k + u - v k + v ) = J, (2.2)
ε ∂v ∂t -β∆v + F v + κ( v k + v - u k + u ) = F L, ε > 0, (2.3) ∂u ∂ν = ∂v ∂ν = 0 on Γ, (2.4) u| t=0 = u 0 , v| t=0 = v 0 .
Note that (2.1)-(2.2) are equivalent to

(2.5) ∂u ∂t -α∆u + κ( k k + v - k k + u ) = J, (2.6) ε ∂v ∂t -β∆v + F v + κ( k k + u - k k + v ) = F L.
We assume that

(2.7) (u 0 , v 0 ) ∈ H 2 N (Ω) 2 , u 0 ≥ 0, v 0 ≥ 0 a.e.
x, where

H 2 N (Ω) = {w ∈ H 2 (Ω), ∂w ∂ν = 0 on Γ}.
We proved in [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF] the Theorem 2.1. We assume that (2.7) holds. Then, (2.1)-(2.4) possesses a unique strong solution (u, v) such that

(2.8) u ≥ 0, v ≥ 0 a.e. (x, t)
and, ∀T > 0,

(u, v) ∈ L ∞ (0, T ; H 2 N (Ω) 2 ) ∩ L 2 (0, T ; H 3 (Ω) 2 ), ( ∂u ∂t , ∂v ∂t ) ∈ L ∞ (0, T ; L 2 (Ω) 2 ) ∩ L 2 (0, T ; H 1 (Ω) 2 ).
Furthermore,

u(t) L ∞ (Ω) ≤ u 0 L ∞ (Ω) + (J + κ)t, t ≥ 0, and
v(t) L ∞ (Ω) ≤ e -F ε t v 0 L ∞ (Ω) + F L + κ F , t ≥ 0.
Finally, if M ≥ F L+κ F and 0 ≤ v 0 ≤ M a.e. x, then 0 ≤ v ≤ M a.e. (x, t).

Remark 2.2. As far as the above regularity is concerned, the corresponding constants in [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF] depend on ε, i.e., they are not bounded uniformly with respect to ε as this quantity goes to 0. However, it is not difficult, reading the details, to see that most constants can be made independent of ε, yielding regularity estimates on u, ∂u ∂t and v which are uniform with respect to ε as ε → 0. Now, we have not been able to derive, at least in a straightforward way, such uniform estimates on ∂v ∂t which would allow us to pass to the limit in (2.2) (say, in a weak (variational) form) to deduce the existence of a solution to the limit problem corresponding to ε = 0 (see however Section 4). We will thus give a direct proof of existence for the limit problem which also has an interest on its own. Remark 2.3. (i) It follows from the above that the capillary lactate concentration is uniformly (with respect to time) bounded. However, we have not been able to derive a similar upper bound on the interstitial lactate concentration u. We can note that, in the biological model, outside a bounded viability domain, cell necrosis occurs (see [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF]), meaning that one expects viable trajectories to be uniformly bounded. (ii) Multiplying (2.1) by u + k, integrating over Ω and by parts, we obtain

dE dt + α ∇u 2 + κ u L 1 (Ω) ≤ ((J + κv k + v , u + k)),
where

E = 1 2 u 2 + ku.
Noting that v is uniformly bounded (we assume that, say, 0 ≤ v 0 ≤ F L+κ F ), we take, for κ, J, F and L given, J small enough and k large enough such that

J + κv k + v < κ.
We thus deduce that

dE dt + α ∇u 2 + c u L 1 (Ω) ≤ c , c > 0,
which yields, noting that

α ∇u 2 + c u L 1 (Ω) ≥ c ( ∇u + u L 1 (Ω) ) -c ≥ c u -c , the differential inequality (2.9) dE dt + c √ E ≤ c , c > 0. Set E * = ( c c ) 2
, where c and c are the same constants as in (2.9), so that

dE * dt + c √ E * = c .

It then follows from comparison arguments that

(2.10)

E(t) ≤ max(E(0), E * ), t ≥ 0,
and we finally deduce that the L 2 -norm of u is uniformly bounded.

The case ε = 0

We consider in this section the following initial and boundary value problem:

(3.1) ∂u ∂t -α∆u + κ( k k + v - k k + u ) = J, (3.2) -β∆v + F v + κ( k k + u - k k + v ) = F L, (3.3) ∂u ∂ν = ∂v ∂ν = 0 on Γ, A. MIRANVILLE (3.4) u| t=0 = u 0 .
We assume that

(3.5) u 0 ∈ H 2 N (Ω), u 0 ≥ 0 a.e. x. Remark 3.1. It follows from (3.2) that -β∆v(0) + F v(0) - k k + v(0) = F L - k k + u 0 .
We will see below that this allows to define in a unique way v(0) such that v(0) ≥ 0 a.e.

x.

3.1.

Existence and uniqueness of solutions to an auxiliary problem. We consider the following modified initial and boundary value problem:

(3.6) ∂u ∂t -α∆u + κ( u k + |u| - v k + |v| ) = J, (3.7)
-β∆v + F v + κ( v k + |v| - u k + |u| ) = F L, (3.8
) ∂u ∂ν = ∂v ∂ν = 0 on Γ, (3.9) u| t=0 = u 0 .
We associate with (3.6)-(3.9) the following weak/variational formulation, for T > 0 given:

Find (u, v) : [0, T ] → H 1 (Ω) 2 such that (3.10) d dt ((u, φ)) + α((∇u, ∇φ)) + ((ϕ k (u), φ)) -((ϕ k (v), φ)) = ((J, φ)), ∀φ ∈ H 1 (Ω), (3.11) β((∇v, ∇ψ)) + F ((v, ψ)) + ((ϕ k (v), ψ)) -((ϕ k (u), ψ)) = ((F L, ψ)), ∀ψ ∈ H 1 (Ω),
in the sense of distributions, and

(3.12) u(0) = u 0 in L 2 (Ω),
where we have set, for c > 0 given,

ϕ c (s) = κs c + |s| , s ∈ R.
We can note that ϕ c is bounded (with |ϕ c | ≤ κ) and of class C 1 , with ϕ c (s) = κc (c+|s|) 2 , so that ϕ c is also Lipschitz continuous, with Lipschitz constant κ c .

Let then 0 = λ 1 < λ 2 ≤ • • • be the eigenvalues of the minus Laplace operator associated with Neumann boundary conditions and w 1 , w 2 , • • • be associated eigenvectors such that the w j 's form an orthonormal in L 2 (Ω) and orthogonal in H 1 (Ω) basis. Setting

V m = Span(w 1 , • • •, w m ), m ∈ N,
we consider the following approximated problem, for T > 0 given:

Find (u m , v m ) : [0, T ] → V m × V m such that (3.13) d dt ((u m , φ)) + α((∇u m , ∇φ)) + ((ϕ k (u m ), φ)) -((ϕ k (v m ), φ)) = ((J, φ)), ∀φ ∈ V m , (3.14) β((∇v m , ∇ψ))+F ((v m , ψ))+((ϕ k (v m ), ψ))-((ϕ k (u m ), ψ)) = ((F L, ψ)), ∀ψ ∈ V m ,
in the sense of distributions, and

(3.15) u m (0) = u 0m ,
where u 0m = P m u 0 , P m being the orthogonal projector (for the

L 2 -norm) from L 2 (Ω) onto V m .
For w ∈ V m given, we consider the following elliptic problem:

Find z ∈ V m such that (3.16) a(z, φ) + ((ϕ k (z), φ)) = ((F L + ϕ k (w), φ)), ∀φ ∈ V m ,
where

a(•, •) = β((∇•, ∇•)) + F ((•, •)
) is bilinear, symmetric, continuous and coercive on V m (and also on H 1 (Ω)). Let then R = R m be the operator defined by

R : V m → V m , z → R(z), where ((R(z), φ)) H 1 (Ω) = a(z, φ) + ((ϕ k (z), φ)) -((F L + ϕ k (w), φ)), ∀φ ∈ V m .
It is clear that this operator is well defined and continuous (since ϕ k is Lipschitz continuous). Furthermore, there holds, for z ∈ V m ,

((R(z), z)) H 1 (Ω) = a(z, z) + ((ϕ k (z), z)) -((F L + ϕ k (w), z)) ≥ c z 2
H 1 (Ω) -c z , c > 0 (note indeed that w is given and recall that ϕ k and ϕ k are bounded). Therefore,

((R(z), z)) H 1 (Ω) ≥ c z 2 H 1 (Ω) -c , so that ((R(z), z)) H 1 (Ω) ≥ 0 whenever z H 1 (Ω) ≥ c c .
It thus follows from the Brouwer fixed point theorem that there exists [START_REF] Temam | Navier-Stokes equations: theory and numerical analysis[END_REF]), which is equivalent to (3.16). Note that all constants here (and also below) are independent of m. This thus defines a mapping F = F m ,

z ∈ V m , z H 1 (Ω) ≤ c c , such that R(z) = 0 in V m (see, e.g.,
F : V m → V m , w → z = F(w). Let then (w 1 , w 2 ) ∈ V m × V m and set z i = F(w i ), i = 1, 2. We have, setting z = z 1 -z 2 and w = w 1 -w 2 , (3.17) a(z, φ) + ((ϕ k (z 1 ) -ϕ k (z 2 ), φ)) = ((ϕ k (w 1 ) -ϕ k (w 2 ), φ)), ∀φ ∈ V m .
Taking φ = z and noting that ϕ k is monotone increasing and ϕ k is Lipschitz continuous, we obtain

z 2 H 1 (Ω) ≤ c w z , whence (3.18) F(w 1 ) -F(w 2 ) H 1 (Ω) ≤ c w 1 -w 2 ,
which yields that F is Lipschitz continuous on V m (both for the L 2 and H 1 -norms); this also yields that F is indeed a mapping, since

w 1 = w 2 implies z 1 = z 2 .
It follows from the above that (3.13)-(3.15) is equivalent to Find u m : [0, T ] → V m such that

(3.19) d dt ((u m , φ)) + α((∇u m , ∇φ)) + ((ϕ k (u m ), φ)) -((ϕ k • F(u m ), φ)) = ((J, φ)), ∀φ ∈ V m , in the sense of distributions, (3.20) u m (0) = u 0m
and then set v m = F(u m ). Since ϕ k and ϕ k are Lipschitz continuous on R and F is Lipschitz continuous on V m with respect to the L 2 -norm, it is easy to prove that (3.19)-(3.20) possesses a (unique) solution u m ∈ L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)) (see, e.g., [START_REF] Temam | Infinite-dimensional dynamical systems in mechanics and physics[END_REF]), whence, setting v m = F(u m ), the existence of a solution (u m , v m ) to (3.13)-(3.15) such that v m ∈ L ∞ (0, T ; H 1 (Ω)). We also note that ∂um ∂t ∈ L 2 (0, T ;

H -1 (Ω)), so that u m ∈ C([0, T ]; L 2 (Ω)). Writing u m (t) = m i=1 d i,m (t)w i , taking φ = λ i w i in (3.19
), multiplying the resulting equality by d i,m and summing over i, we obtain 1 2

d dt ∇u m 2 + α ∆u m 2 -((ϕ k (u m ), ∆u m )) + ((ϕ k • F(u m ), ∆u m )) = 0,
which yields, recalling that ϕ k and ϕ k are bounded,

d dt ∇u m 2 + α ∆u m 2 ≤ c,
whence estimates on u m in L ∞ (0, T ; H 1 (Ω)) and L 2 (0, T ; H 2 (Ω)). It thus follows that

∂um ∂t ∈ L ∞ (0, T ; H -1 (Ω)) ∩ L 2 (0, T ; L 2 (Ω)).
Since the above regularity estimates are uniform with respect to m, we deduce from classical Aubin-Lions compactness theorems that, at least for a subsequence which we do not relabel,

u m → u in L ∞ (0, T ; H 1 (Ω)) weak star, in L 2 (0, T ; H 2 (Ω)) weak, in C([0, T ]; L 2 (Ω)) and a.e. (x, t) ∈ Ω × (0, T ), v m → v in L ∞ (0, T ; H 1 (Ω)
) weak star, for some functions u and v. Actually, since v m = F(u m ) and F is Lipschitz continuous with respect to the L 2 -norm, we can see that

(v m ) is a Cauchy sequence in C([0, T ]; L 2 (Ω)) (note indeed that, if m ≥ m, then V m ⊂ V m and that the constant c in (3.18) is independent of m), so that v m → v in C([0, T ]; L 2 (Ω)).
Recalling finally that ϕ k and ϕ k are Lipschitz continuous, this is sufficient to pass to the limit in (3.13)-(3.15) and deduce the existence of a solution (u, v) to (3.10)-(3.12) (note that the initial condition u 0 = u(0) makes sense; actually, v(0) also makes sense). Indeed, we need to pass to the limit in relations of the form 2 and θ ∈ D(0, T ). More precisely, we have the Theorem 3.2. We assume that u 0 ∈ H 1 (Ω). Then, (3.10)-(3.12) possesses a unique solution (u, v) such that, ∀T > 0,

T 0 [-((u m , φ))θ (t) + α((∇u m , ∇φ))θ(t) + (ϕ k (u m ), φ))θ(t) -((ϕ k (v m ), φ))θ(t) -((J, φ))θ(t)] dt = 0 and T 0 [β((∇v m , ∇ψ)) + F ((v m , ψ)) + ((ϕ k (v m ), ψ)) -((ϕ k (u m ), ψ)) -((F L, ψ))]θ(t) dt = 0, for (φ, ψ) ∈ H 1 (Ω)
u ∈ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)) ∩ C([0, T ]; L 2 (Ω)), v ∈ L ∞ (0, T ; H 1 (Ω)) ∩ C([0, T ]; L 2 (Ω)) and ∂u ∂t ∈ L ∞ (0, T ; H -1 (Ω)) ∩ L 2 (0, T ; L 2 (Ω)).
Proof. There remains to prove the uniqueness. Let thus (u 1 , v 1 ) and (u 2 , v 2 ) be two such solutions, with initial data u 0,1 and u 0,2 , respectively. We have, setting (u, v)

= (u 1 -u 2 , v 1 -v 2 ) and u 0 = u 0,1 -u 0,2 , (3.21) d dt ((u, φ)) + α((∇u, ∇φ)) + ((ϕ k (u 1 ) -ϕ k (u 2 ), φ)) -((ϕ k (v 1 ) -ϕ k (v 2 ), φ)) = 0, ∀φ ∈ H 1 (Ω), (3.22) β((∇v, ∇ψ)) + F ((v, ψ)) + ((ϕ k (v 1 ) -ϕ k (v 2 ), ψ)) -((ϕ k (u 1 ) -ϕ k (u 2 ), ψ)) = 0, ∀ψ ∈ H 1 (Ω), (3.23) u(0) = u 0 .
Taking φ = u and ψ = v, we obtain, recalling that ϕ k and ϕ k are monotone increasing and Lipschitz continuous, We deduce from (3.26) and (3.28) the uniqueness, as well as the continuous dependence with respect to the initial data in the L 2 -norm.

3.2.

Existence and uniqueness of nonnegative solutions. We first prove additional regularity results on the solutions to (3.10)-(3.11), assuming that (3.5) holds. This can be fully justified within the Galerkin scheme considered above. Taking ψ = -∆v in (3.11), we have

β ∆v 2 + F ∇v 2 = ((ϕ k (v), ∆v)) -((ϕ k (u), ∆v
)), which yields, recalling that ϕ k and ϕ k are bounded,

β 2 ∆v 2 + F ∇v 2 ≤ c, whence estimates on v in L ∞ (0, T ; H 2 (Ω)), ∀T > 0.
Taking then φ = ∆ 2 u in (3.10) and ψ = ∆ 2 v in (3.11), we obtain 1 2

d dt ∆u 2 + α ∇∆u 2 = -((ϕ k (u), ∆ 2 u)) + ((ϕ k (v), ∆ 2 u))
and

β ∇∆v 2 + F ∆v 2 = -((ϕ k (v), ∆ 2 v)) + ((ϕ k (u), ∆ 2 v)), respectively. Noting that |((ϕ k (u), ∆ 2 u))| = |((ϕ k (u)∇u, ∇∆u))| ≤ c ∇u ∇∆u ,
we find, proceeding in a similar way for the other terms,

d dt ∆u 2 + α ∇∆u 2 ≤ c( ∇u 2 + ∇v 2 )
and

β 2 ∇∆v 2 + F ∆v 2 ≤ c( ∇u 2 + ∇v 2 ),
whence estimates on u and v in L ∞ (0, T ; H 2 (Ω)) ∩ L 2 (0, T ; H 3 (Ω)) and L ∞ (0, T ; H 3 (Ω)), respectively, ∀T > 0.

Remark 3.3. This yields that ∂u ∂t ∈ L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)), ∀T > 0. We further note that the solution to (3.10)-(3.12) is strong, i.e., (3.6)-(3.9) are satisfied almost everywhere.

We can now prove the Theorem 3.4. We assume that (3.5) holds. Then, (3.1)-(3.4) possesses a unique strong solution (u, v) such that u ≥ 0, v ≥ 0 a.e. (x, t) and, ∀T > 0,

u ∈ L ∞ (0, T ; H 2 N (Ω)) ∩ L 2 (0, T ; H 3 (Ω)), v ∈ L ∞ (0, T ; H 3 (Ω) ∩ H 2 N (Ω)) ∩ C([0, T ]; L 2 (Ω)) and ∂u ∂t ∈ L ∞ (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)).
Proof. Let (u, v) be the unique strong solution to (3.6)-(3.9). Multiplying (3.6) by -u - and (3.7) by -v -, where x -= max(0, -x), we have

(3.29) 1 2 d dt u -2 + α ∇u -2 + κ Ω |u -| 2 k + |u| dx + κ Ω vu - k + |v| dx ≤ 0 and (3.30) β ∇u -2 + F u -2 + κ Ω |v -| 2 k + |v| dx + κ Ω uv - k + |u| dx ≤ 0,
respectively. Writing v = v + -v -, where x + = max(0, x), we deduce from (3.29) that 1 2

d dt u -2 ≤ κ Ω u -v - k + |v| dx, whence (3.31) d dt u -2 ≤ c u -v -.
Proceeding in a similar way for (3.30), we find

F v -2 ≤ c u -v -, whence (3.32) v -≤ c u -.
Injecting this into (3.31), we deduce that

d dt u -2 ≤ c u -2 ,
which yields, owing to Gronwall's lemma,

(3.33) u -(t) ≤ e ct u -(0) , t ≥ 0,
whence, since u -(0) = 0, u ≥ 0 a.e. (x, t). This, together with (3.32), yields that v ≥ 0 a.e. (x, t). Consequently, (u, v) is a strong solution to (3.1)-(3.4), which finishes the proof.

Remark 3.5. Proceeding exactly as in [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF], we can prove that

u(t) L ∞ (Ω) ≤ u 0 L ∞ (Ω) + (J + κ)t, t ≥ 0. Furthermore, it follows from (3.2) that (3.34) -β∆v + F v ≤ F L + κ.
Multiplying (3.34) by v m+1 , m ∈ N, we have

β(m + 1) Ω v m |∇v| 2 dx + F v m+2 L m+2 (Ω) ≤ (F L + κ) Ω v m+1 dx, which yields F v m+2 L m+2 (Ω) ≤ (F L + κ)Vol(Ω) 1 m+2 v m+1 L m+2 (Ω) , whence (3.35) v L m+2 (Ω) ≤ F L + κ F Vol(Ω) 1 m+2 .
Passing to the limit m → +∞ in (3.35), we finally obtain (see, e.g., [START_REF] Alikakos | L p bounds of solutions to reaction-diffusion equations[END_REF])

v(t) L ∞ (Ω) ≤ F L + κ F , t ≥ 0,
meaning that the capillary lactate concentration is again uniformly (with respect to time) bounded. Also note that (2.10) still holds.

Remark 3.6. (i) As mentioned in the introduction (for the case ε > 0, but the situation is the same here), the existence of solutions for negative initial data is a challenging issue. However, we can prove the following partial result (see also [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF] for the case ε > 0). Let δ 1 and δ 2 be two positive constants such that k -δ 1 > 0 and k -δ 2 > 0 and assume that u 0 ≥ -δ 1 a.e. x. We then consider the following modified initial and boundary value problem:

(3.36) ∂u ∂t -α∆u + κ( u k -δ 1 + |u + δ 1 | - v k -δ 2 + |v + δ 2 | ) = J, (3.37)
-β∆v + F v + κ( v k -δ 2 + |v + δ 2 | - u k -δ 1 + |u + δ 1 | ) = F L, (3.38) ∂u ∂ν = ∂v ∂ν = 0 on Γ, (3.39) u| t=0 = u 0 .
The existence and uniqueness of the solution to (3.36)-(3.39) can be proved by arguing as above. Next, we set ũ = u + δ 1 and ṽ = v + δ 2 . These functions are solutions to

(3.40) ∂ ũ ∂t -α∆ũ + κ( ũ k -δ 1 + |ũ| - ṽ k -δ 2 + |ṽ| ) = J, (3.41) -β∆ṽ + F ṽ + κ( ṽ k -δ 2 + |ṽ| - ũ k -δ 1 + |ũ| ) = F , (3.42) ∂ ũ ∂ν = ∂ṽ ∂ν = 0 on Γ, (3.43) ũ| t=0 = u 0 + δ 1 ,
where

J = J + κ( δ 1 k -δ 1 + |ũ| - δ 2 k -δ 2 + |ṽ| ) and F = F (L + δ 2 ) -κ( δ 1 k -δ 1 + |ũ| - δ 2 k -δ 2 + |ṽ| ).
Choosing δ 1 and δ 2 such that J ≥ 0 and F ≥ 0 (in particular, these hold when δ 1 and δ 2 are small enough) and noting that ũ(0) ≥ 0 a.e. x, we can prove, as in the proof of Theorem 3.4, that ũ(x, t) ≥ 0 and ṽ(x, t) ≥ 0 a.e. (x, t), so that (u, v) is solution to (

, t) ≥ -δ 1 and v(x, t) ≥ -δ 2 a.e. (x, t).

(ii) Similarly, we can prove that, if δ 1 and δ 2 are positive and small enough, with u 0 ≥ δ 1 a.e. x, then u(x, t) ≥ δ 1 and v(x, t) ≥ δ 2 a.e. (x, t).

It follows from the above that we can actually define the Lipschitz continuous (for the L 2 and H 1 -norms) mapping

F : H 1 (Ω) → H 1 (Ω), w → z = F(w),
where z is the unique solution to the following elliptic problem:

(3.44) a(z, φ) + ((ϕ k (z), φ)) = ((F L + ϕ k (w), φ)), ∀φ ∈ H 1 (Ω).
We then have the Proposition 3.7. The mapping F is differentiable with respect to the L 2 and H 1 -norms.

Proof. Let w 0 and w belong to H 1 (Ω) and set z 0 = F(w 0 ) and z = F(w). We then have

(3.45) a(z -z 0 , φ) + ((ϕ k (z) -ϕ k (z 0 ), φ)) = ((ϕ k (w) -ϕ k (w 0 ), φ)), ∀φ ∈ H 1 (Ω).
Taking φ = z -z 0 and reacalling tht ϕ k is monotone increasing, this yields

z -z 0 H 1 (Ω) ≤ c ϕ k (w) -ϕ k (w 0 ) , whence (3.46) z -z 0 H 1 (Ω) ≤ c w -w 0 .
Let then Z ∈ H 1 (Ω) be the solution to the linear elliptic problem (recall that ϕ k is nonnegative)

(3.47) a(Z, φ) + ((ϕ k (z 0 )Z, φ)) = ((ϕ k (w 0 )(w -w 0 ), φ)), ∀φ ∈ H 1 (Ω).
Setting h = w -w 0 , we can see that

(3.48) a(z -z 0 -Z, φ) + ((ϕ k (z) -ϕ k (z 0 ) -ϕ k (z 0 )Z, φ)) = ((ϕ k (w) -ϕ k (w 0 ) -ϕ k (w 0 )h, φ)), ∀φ ∈ H 1 (Ω).
Writing

ϕ k (z) -ϕ k (z 0 ) -ϕ k (z 0 )Z = ϕ k (z 0 )(z -z 0 -Z) + o( z -z 0 )
and

ϕ k (w) -ϕ k (w 0 ) -ϕ k (w 0 )h = o( h ),
we obtain, taking φ = z -z 0 -Z and employing (3.46) (also recall that ϕ k ≥ 0),

a(z -z 0 -Z, z -z 0 -Z) ≤ |((o( h ), z -z 0 -Z))|, whence z -z 0 -Z H 1 (Ω) = o( h ).
This yields that F is differentiable at w 0 , with F (w 0 ) • h = Z, F denoting the differential of F.

We deduce from Proposition 3.

(Ω)) ≤ Q(T, u 0 H 1 (Ω) ) ∂u ∂t L 2 (0,T ;H 1 (Ω)) .

Convergence to the limit problem

All constants c and c and functions Q in this section are independent of ε. Let (u ε , v ε ) and (u 0 , v 0 ) be the unique strong solutions to the initial and limit problems, respectively, as given in Theorems 2.1 and 3.4, where v 0 = v 0 (0), i.e., (4.1) We have the Theorem 4.1. The following error estimates hold, ∀T > 0:

∂u ε ∂t -α∆u ε + κ( k k + v ε - k k + u ε ) = J, (4.2) ε ∂v ε ∂t -β∆v ε + F v ε + κ( k k + u ε - k k + v ε ) = F L, (4.3) ∂u ε ∂ν = ∂v ε ∂ν = 0 on Γ, (4.4) u ε | t=0 = u 0 , v ε | t=0 = v 0 (0), and
(4.5) ∂u 0 ∂t -α∆u 0 + κ( k k + v 0 - k k + u 0 ) = J, (4.6) -β∆v 0 + F v 0 + κ( k k + u 0 - k k + v 0 ) = F L,
u ε (t) -u 0 (t) H 1 (Ω) ≤ Q(T, u 0 H 1 (Ω) )ε, v ε (t) -v 0 (t) H 1 (Ω) ≤ Q(T, u 0 H 1 (Ω) ) √ ε, t ∈ [0, T ],
and

u ε -u 0 L 2 (0,T ;H 2 (Ω)) ≤ Q(T, u 0 H 1 (Ω) )ε, v ε -v 0 L 2 (0,T ;H 2 (Ω)) ≤ Q(T, u 0 H 1 (Ω) )ε. Proof. We have, setting u = u ε -u 0 and v = v ε -v 0 , (4.9) ∂u ∂t -α∆u + ϕ k (u ε ) -ϕ k (u 0 ) -ϕ k (v ε ) + ϕ k (v 0 ) = 0, (4.10) ε ∂v ∂t -β∆v + F v + ϕ k (v ε ) -ϕ k (v 0 ) -ϕ k (u ε ) + ϕ k (u 0 ) = -ε ∂v 0 ∂t , ε > 0, ( 4
ε d dt v 2 + c v 2 H 1 (Ω) ≤ c ( u v + ε 2 ∂v 0 ∂t 2 ), c > 0.
Combining (4.13) and (4.14), we have, owing to (3.49), (4.15)

d dt ( u 2 + ε v 2 ) + c( u 2 H 1 (Ω) + v 2 H 1 (Ω) ) ≤ c ( u 2 + ε v 2 + ε 2 ∂u ε ∂t 2 ), c > 0,
from which it follows, owing to Gronwall's lemma, (4.16)

u(t) 2 + ε v(t) 2 + c T 0 ( u(s) 2 H 1 (Ω) + v(s) 2 H 1 (Ω) ) ds ≤ Q(T )ε 2 ∂u ε ∂t 2 L 2 (0,T ;L 2 (Ω)) , c > 0, t ∈ [0, T ]. Multiplying now (4.1) by ∂u ε ∂t , we obtain d dt ∇u ε 2 + c ∂u ε ∂t 2 ≤ c , c > 0, whence (4.17) ∂u ε ∂t 2 L 2 (0,T ;L 2 (Ω)) ≤ Q(T, u 0 H 1 (Ω) ).
We finally deduce from (4.16)-(4.17) that

u(t) 2 + ε v(t) 2 + c T 0 ( u(s) 2 H 1 (Ω) + v(s) 2 H 1 (Ω) ) ds (4.18) ≤ Q(T, u 0 H 1 (Ω) )ε 2 , c > 0, t ∈ [0, T ].
Multiplying next (4.9) by -∆u and (4.10) by -∆v, we find, recalling that ϕ k and ϕ k are Lipschitz continuous,

d dt ∇u 2 + c ∆u 2 ≤ c ( u 2 + v 2 ), c > 0, and
ε d dt ∇v 2 + c ∆v 2 ≤ c ( u 2 + v 2 + ε 2 ∂v 0 ∂t 2 ), c > 0,
respectively. Summing these two inequalities, integrating over [0, T ] and proceeding as above, we have, adding the resulting inequality to (4.18),

u(t) 2 H 1 (Ω) + ε v(t) 2 H 1 (Ω) + c T 0 ( u(s) 2 H 2 (Ω) + v(s) 2 H 2 (Ω) ) ds ≤ c ( T 0 ( u(s) 2 + v(s) 2 ) ds + Q(T, u 0 H 1 (Ω) )ε 2 ), c > 0.
This yields, employing (4.18) to estimate the right-hand side,

u(t) 2 H 1 (Ω) + ε v(t) 2 H 1 (Ω) + c T 0 ( u(s) 2 H 2 (Ω) + v(s) 2 H 2 (Ω) ) ds ≤ Q(T, u 0 H 1 (Ω) )ε 2 , c > 0, t ∈ [0, T ],
which finishes the proof.

Remark 4.2. We have similar error estimates if we assume that u ε (0)-u 0 (0

) H 1 (Ω) ≤ cε and v ε (0) -v 0 (0) H 1 (Ω) ≤ c √ ε.

A stability result

As in [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF], (3.1)-(3.2) possesses a unique spatially homogeneous equilibrium (u, v) given by

v = L + J F > 0 and u = k( J κ + v k +v ) 1 -( J κ + v k +v )
.

Note that u is not necessarily positive. We thus assume in what follows that u > 0. The linearized (around (u, v)) system reads (5.1)

∂U ∂t -α∆U + κ( k (k + u) 2 U - k (k + v) 2 V ) = 0, (5.2) -β∆V + F V + κ( k (k + v) 2 V - k (k + u) 2 U ) = 0, (5.3) ∂U ∂ν = ∂V ∂ν = 0 on Γ, (5.4) U | t=0 = U 0 .
It is not difficult here to prove the existence, uniqueness and regularity of the solution to (5.1)-(5.4), assuming that U 0 is regular enough. Furthermore, proceeding as above, we can prove that, if U 0 ≥ 0 a.e. x, then U (x, t) ≥ 0 and V (x, t) ≥ 0 a.e. (x, t).

Multiplying now (5.1) by k (k+u) 2 U and (5.2) by k (k +v) 2 V , we obtain, summing the two resulting equalities, (5.5) 1 2 Multiplying next (5.2) by V , we easily find V ≤ c U , so that (5.7)

k (k + u) 2 d dt U 2 + αk (k + u) 2 ∇U 2 + βk (k + v) 2 ∇V 2 + F k (k + v) 2 V 2 + Ω ( k (k + u) 2 U - k (k + v) 2 V )
V (t) ≤ c U 0 , t ≥ 0.

We deduce from (5.6)-(5.7) that (u, v) is linearly stable with respect to the L 2 -norm. We can also prove the linear stability with respect to the H 1 -norm, proceeding in a similar way. Now, an important question is whether we also have a linear exponential stability as in [START_REF] Guillevin | On a reaction-diffusion system associated with brain lactate kinetics[END_REF] for the case ε > 0 (see also [START_REF] Costalat | Mathematical modeling of metabolism and hemodynamics[END_REF], [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF], [START_REF] Lahutte-Auboin | Dip and buffering in a fast-slow system associated to brain lactate kinetics[END_REF] and [START_REF] Lahutte-Auboin | On a minimal model for hemodynamics and metabolism of lactate: application to low grade glioma and therapeutic strategies[END_REF] for the ODE's model (1.1)-(1.2)). Indeed, as mentioned in the introduction, a therapeutic perspective of such a result is to have the (spatially homogeneous) steady state outside the viability domain, where cell necrosis occurs (see [START_REF] Lahutte-Auboin | Modélisation biomathématique du métabolisme énergétique cérébral : réduction de modèle et approche multi-échelle, application à l'aide à la décision pour la pathologie des gliomes[END_REF]).

We have, in this direction, the Theorem 5.1. The stationary solution (u, v) is linearly exponentially stable, in the sense that all eigenvalues s ∈ C associated with the linear system (5.1)-(5.2) satisfy Re(s) ≤ -ξ, for a given ξ > 0, Re denoting the real part.

Proof. We first note that it follows from (5.2) that (5.8)

V = k 1 (-β∆ + (F + k 2 )I) -1 U,
where k 1 = κk (k+u) 2 and k 2 = κk (k +v) 2 . Injecting this into (5.1), we obtain (5.9) ∂U ∂t -α∆U + k 1 U -k 1 k 2 (-β∆ + (F + k 2 )I) -1 U = 0.

We then look for solutions of the form U (x, t) = Û (x)e st , for s ∈ C, s = ζ + iη. Injecting this into (5.9), we find

β ∇v 2 +d dt u 2 ≤ c u 2 ,

 222 F v 2 ≤ c u v , respectively. In particular, it follows from (3.25) that (3.26) v ≤ c u , which, injected into (3.24), yields (3.27) whence, owing to Gronwall's lemma, (3.28) u(t) ≤ e ct u 0 , t ≥ 0.

  u 0 | t=0 = u 0 .

  U (t) ≤ U 0 , t ≥ 0.

(5. 10 ) 1 Ûαβ∆ 2 Û- k 1 F F + k 2 )

 101212 -α∆ Û + (s + k 1 ) Û -k 1 k 2 (-β∆ + (F + k 2 )I) --(βs + αF + βk 1 + αk 2 )∆ Û + ((F + k 2 )s + k 1 F ) Û = 0,where, owing to (5.10) and (5.11), 12) by the conjugate of Û , integrating over Ω and by parts and taking the real part, we have(5.14) αβ |∆ Û | 2 + (βζ + αF + βk 1 + αk 2 ) |∇ Û | 2 + ((F + k 2 )ζ + k 1 F ) | Û | 2 = 0.Therefore, when ζ ≥ 0, then, necessarily, Û ≡ 0. Furthermore, (5.14) can have nontrivial solutions only whenβζ + αF + βk 1 + αk 2 ≤ 0 or (F + k 2 )ζ + k 1 F ≤ 0.Therefore, necessarily,(5.15) ζ ≤ max(-αF + βk 1 + αk 2 β , < 0, which finishes the proof.