From Kinetic Theory of Multicellular Systems to Hyperbolic Tissue Equations: Asymptotic Limits and Computing

From Kinetic Theory of Multicellular Systems to Hyperbolic Tissue Equations: Asymptotic Limits and Computing

Nisrine Outada, Nicolas Vauchelet,
Thami Akrid and Mohamed Khaladi
Département de Mathématiques, Faculté des Sciences Semlalia,
Laboratoire LMDP, Université Cadi Ayyad,
B.P. 2390, 40000 Marrakesh, Morocco
Sorbonne Universités, UPMC Univ Paris 06, UMR 7598,
Laboratoire Jacques-Louis Lions, Paris, France
UMI 209 UMMISCO, 32 Avenue Henri Varagnat,
F-93 143 Bondy Cedex, France
INRIA-Paris-Rocquencourt, EPC MAMBA, Domaine de Voluceau,
BP105, 78153 Le Chesnay Cedex, France

This paper deals with the analysis of the asymptotic limit toward the derivation of macroscopic equations for a class of equations modeling complex multicellular systems by methods of the kinetic theory. After having chosen an appropriate scaling of time and space, a Chapman-Enskog expansion is combined with a closed, by minimization, technique to derive hyperbolic models at the macroscopic level. The resulting macroscopic equations show how the macroscopic tissue behavior can be described by hyperbolic systems which seem the most natural in this context. We propose also an asymptotic-preserving well-balanced scheme for the one-dimensional hyperbolic model, in the two dimensional case, we consider a time splitting method between the conservative part and the source term where the conservative equation is approximated by the Lax-Friedrichs scheme.

Dedicated to Abdelghani Bellouquid who prematurely passed away on August 2015.

Keywords. Kinetic theory; multicellular systems; hyperbolic limits; chemotaxis; asymptotic preserving scheme; Lax-Friedrichs flux.

1 Introduction

The aim of this paper is the derivation of macroscopic hyperbolic models of biological tissues from the underlying description at the microscopic scale delivered by kinetic theory methods. We consider the hyperbolic asymptotic limit for microscopic system that connect the biological parameters, at the level of cells, involved in this level of description.

The first step of the derivation of macroscopic models in biology from the underlying description at the microscopic scale is arguably due to Alt [1] and Othmer, Dunbar and Alt [27], who introduced a new modeling approach by perturbation of the transport equation by a velocity jump-process, which appears appropriate to model the velocity dynamics of cells modeled as living particles. This method has been subsequently developed by various authors, among others, we cite [3, 4, 5, 6, 7, 9, 12, 14, 15, 19, 20, 21, 23, 28, 30, 32, 35]. The survey on mathematical challenges on the qualitative and asymptotic analysis of Keller and Segel type models [6] reports an exhaustive bibliography concerning different mathematical approach on the aforementioned topics. The interested reader can find a further updating of the research activity on the study of Keller-Segel models and their developments in [2, 8, 25, 37, 38], as well as on applications to population dynamics with diffusion [34] and pattern formation in cancer [33].

Different time-space scalings lead to equations characterized by different parabolic or hyperbolic structures. Different combinations of parabolic and hyperbolic scales also are used, according to the dispersive or non-dispersive nature of the biological system under consideration. The parabolic (low-field) limit of kinetic equations leads to a drift-diffusion type system (or reaction-diffusion system) in which the diffusion process dominate the behavior of the solutions[16, 36]. On the other hand, in the hyperbolic (high-field) limit the influence of the diffusion terms is of lower (or equal) order of magnitude in comparison with other convective or interaction terms.

Possible applications refer to modeling cell invasion, as well as chemotaxis and haptotaxis phenomena and related pattern formation[1, 5, 14]. Models with finite propagation speed appear to be consistent with physical reality rather than parabolic models. This feature is also induced by the essential characteristics of living organisms who have the ability to sense signals in the environment and adapt their movements accordingly.

Our analysis is quite general as it can be applied to different species in response to multiple (chemo)tactic cues [10, 11, 22, 29]. Therefore, the derivation of hyperbolic models can contribute to further improvements in modeling biological reality. In fact, it seems that the approach introduced by Patlak [29] and Keller-Segel [24] is not always sufficiently precise to describe some structures as the evolution of bacteria movements, or the human endothelial cells movements on matrigel that lead to the formation of networks interpreted as the beginning of a vasculature [30]. These structures cannot be explained by parabolic models, which generally lead to pointwise blow-up [6], moreover the numerical experiments show the predictability of the hyperbolic models in this context.

We now briefly describe the contents of this paper. Section 2, presents the kinetic model and the scaling deemed to provide the general framework appropriate to derive, by asymptotic analysis, models at the macroscopic scale. Section 3, referring to [3], presents the general kinetic framework to be used toward the asymptotic analysis. Section 4, shows how specific models can be derived by the approach of our paper. Section 5, presents some computational simulations to show the predictive ability of the models derived in this paper and looks ahead to research perspectives.

2 Kinetic mathematical model

Let us consider a physical system constituted by a large number of cells interacting in a biological environment. The microscopic state is defined by the mechanical variable , where , . The statistical collective description of the system is encoded in the statistical distribution , which is called a distribution function. We also assume that the transport in position is linear with respect to the velocity. In this paper, we are interested in the system of different species in response to multiple chemotactic cues. The model, for reads:


where and denotes respectively the density of cells and the density (concentration) of multiple tactic cues and .

The operators and model the dynamics of biological organisms by velocity-jump process. The set of possible velocities is denoted by , assumed to be bounded and radially symmetric. The operators and describe proliferation/destruction interactions. The dimensionless time indicates that the spatial spread of and are on different time scales. The case corresponds to a steady state assumption for .

The problem of studying the relationships between the various scales of description, seems to be one of the most important problems of the mathematical modelling of complex systems . Different structures at the macroscopic scale can be obtained corresponding to different space–time scales. Subsequently, more detailed assumptions on the biological interactions lead to different models of pattern formation. However, a more recent tendency been the use hyperbolic equations to describe intermediate regimes at the macroscopic level rather than parabolic equations, for example [3, 4, 15, 21].

The next section deals with the derivation of macroscopic equations using a Champan-Enskog type perturbation approach for (2.1) and a closure by minimization method for (2.1). Our purpose is to derive hyperbolic-hyperbolic macroscopic model. The first approach consists in expanding the distribution function in terms of a small dimensionless parameter related to the intermolecular distances (the space scale dimensionless parameter). In [15], a hydrodynamic limit of such kinetic model was used to derive hyperbolic models for chemosensitive movements. While the closure method consists that the (m+1)-moments of the minimizer approximate the (m+1)-moments of the true solution.

3 Asymptotic analysis toward derivation of hyperbolic systems

3.1 The kinetic framework

Let us now consider the first equation in (2.1). We assume a hyperbolic scaling for this population it means that we scale time and space variables and , where is a small parameter which will be allowed to tend to zero, see [3] for more details. We deal also with the small interactions i.e . Then, we obtain the following transport equation for the distribution function


where the position and the velocity . In addition, the analysis developed is based on the assumption that admits the following decomposition:


with in the form


The operator represents the dominant part of the turning kernel modeling the tumble process in the absence of chemical substance and is the perturbation due to chemical cues. The parameter is a time scale which here refers to the turning frequency. The equation (3.2) becomes


The most commonly used assumption on the perturbation turning operators , and is that they are integral operators and read:




The turning kernels , and describe the reorientation of cells, i.e. the random velocity changes from the previous velocity to the new .

The following assumptions on the turning operators are needed to develop the hyperbolic asymptotic analysis:
Assumption H0: For all , the turning operators , and conserve the local mass:


Assumption H1: The turning operator conserve the population flux:


Assumption H2: For all and , there exists a unique function such that


It is clear, from (3.6)-(3.8), that and satisfy the assumption

The following lemma, whose proof can be found in [22], will be used a few times,

Lemma 3.1

Assume that , , which corresponds to the assumption that any individual of the population chooses any velocity with a fixed norm s(speed). Then,

where and denotes the Kronecker symbol, and the notation corresponds to the unit sphere in dimension d.

3.2 Hydrodynamic limit

In this subsection, we use the last assumptions to derive an hyperbolic system on macroscopic scale for small perturbation parameter.

Let be solution of the equation (3.5) and consider the density of cells and the flux defined by:


To derive the equations for the moments in (3.12), we multiply (3.5) by and respectively, and integrate over to obtain the following system


Now let be a solution of the following (i)-equation


and set


To derive the equations for moments in (3.15), we multiply the equation (3.14) by and respectively and integrate over to obtain the following system


Finally, (3.13) and (3.16) yield the following system


In the following, we are interested to close the system (3.17). We start by two first equations of (3.17), we introduce such that

where the equilibrium distribution is defined by (3.11). Then, we deduce

Then, we assume the following asymptotic expansion in order 1 in ,


Replacing now by its expansion and using the first equality of (3.18), yields



where the pressure tensor is given by


Since conserves the local mass (3.9), the system (3.19) becomes

Remark 3.2

It is easy to see that the influence of the turning operator on the macroscopic equations (3.21) only comes into play through the stationary state in the computation of the right-hand side of the second equation in (3.21) and the pressure tensor . While the structure of the turning operator determines the effect of the chemical cues.

Taking into account the system (3.21) and using the second equality of (3.18) the system (3.17) reads now



It can be observed that system (3.22) is not yet closed. Indeed, it can be closed by looking for an approximate expression of . The approach consists in deriving a function which minimizers the -norm under the constraints that it has the same first moments, and , as . Once this function has been found, we replace by , and by in the others terms.

Toward this aim, we consider the set of velocities with and the unit sphere of . Let us introduce Lagrangian multipliers and respectively scalar and vector, and define the following operator:

The Euler-Lagrange equation (first variation) of reads . We use the constraints to define and . First, from the first equality in (3.15) one gets easily . Next, from Lemma 3.1 one obtains

then .



Consequently, using again lemma 3.1, the pressure tensor is

where denotes the identity matrix. Thus, the following nonlinear coupled hyperbolic model is derived:


with .

Remark 3.3

The second variation of is , then the extremum is a minimum.

4 Derivation of models

This section shows how the tools reviewed in the preceding section can be used to derive models. Let us consider the model defined by choosing the stationary state and the turning kernels. Consider as follows:


It is easy to check that satisfies the assumptions (3.10)-(3.11).
We take the turning kernel in (3.6) in the form

with a real constant, and consider that the turning kernel in (3.7) depends on the velocity , on the population , and on its gradient, defined by:

where is a mapping , are real constants and stands for the -mean of a function, i.e for .
Therefore, the turning operator is given by:


then, is a relaxation operator to .
While, the turning operator can be computed as follows:




Finally, take the turning kernel in (3.8) as follows:

with is a real constant.
Then, the turning operator is computed as follows:




Now we compute the pressure tensor . By using lemma 3.1, we have



Finally, the system (3.24) becomes, at first order with respect to ,


where and are defined in (3.23) and (4.25).

Theorem 4.1

If we consider for all , , and satisfy the assumption (3.18) then, we obtain the following system at first order with respect to ,


This theorem leads to some specific models which are presented in the next subsection.

4.1 A Cattaneo type model for chemosensitive movement

Taking in (4.31), one can derive the corresponding hyperbolic system for chemosensitive movement, at first order with respect to , as follows




In absence of interactions, the authors in [15] and [22] derived, respectively, the first two equations for by asymptotic analysis and moment closure. The system composed by the first two equations with is called the Cattaneo model for chemosensitive movement with density control [11, 21].

4.2 Derivation of Keller-Segel models

The approach proposed can be applied to derive a variety of models of Keller-Segel type. Indeed, by taking the system (4.31) and with different scalings this approach allows to derive various models.
From (3.23) and (4.25), we have and . To get our aim, we assume moreover in this subsection the following assumption,


and we set

Consequently, we have the following proposition,

Proposition 4.2

For , the system (4.31) becomes, with above assumptions (3.18) and (4.33), which are satisfied if and are bilinear


Let now and such that . Dividing the fourth equation of system (4.34) by and taking last limits, yields therefore the third equation of (4.34) writes