On Fokker-Planck Equations with In- and Outflow of Mass

TL;DR

This paper studies nonlinear Fokker-Planck equations with non-conserved mass, modeling neuron growth, proving exponential decay to equilibrium despite the lack of mass conservation, using entropy methods and numerical validation.

Contribution

It provides new analytical results on decay to equilibrium for Fokker-Planck equations with mass in- and outflow, overcoming challenges due to non-conservation of mass.

Findings

Proved exponential decay towards equilibrium in various scenarios.

Developed entropy-based methods applicable without classical Sobolev inequalities.

Validated results through extensive numerical simulations.

Abstract

Motivated by modeling transport processes in the growth of neurons, we present results on (nonlinear) Fokker-Planck equations where the total mass is not conserved. This is either due to in- and outflow boundary conditions or to spatially distributed reaction terms. We are able to prove exponential decay towards equilibrium using entropy methods in several situations. As there is no conservation of mass it is difficult to exploit the gradient flow structure of the differential operator which renders the analysis more challenging. In particular, classical logarithmic Sobolev inequalities are not applicable any more. Our analytic results are illustrated by extensive numerical studies.

Equations

$${\partial }_t\rho +\nabla \cdot (-\nabla \rho +f(\rho )\nabla V)=0,$$

$${\partial }_t\rho =\nabla \cdot (f(\rho )\nabla E^{\prime }(\rho ))=0,E(\rho )={\int }_{\Omega }F(\rho )-\rho Vdx,$$

$$E(\rho |{\rho }_{\infty })={\int }_{\Omega }F(\rho )-F({\rho }_{\infty })-F^{\prime }({\rho }_{\infty })(\rho -{\rho }_{\infty })dx.$$

$${\partial }_t\rho =\nabla \cdot J,J=-\nabla \rho -f(\rho )u,$$

$$J\cdot n=-a\frac{f(\rho )}{\rho }+b\rho ,$$

$${\Omega }_3=\{(x_1,x_2,x_3):x_1\in (0,1),(x_2,x_3)\in {\Omega }_2(x_1)\},$$

$$\rho (x,t)\approx {\rho }_1(x_1,t)q(x,t),$$

$${\nabla }_{23}\cdot ({\nabla }_{23}q+\frac{f({\rho }_{1}q)}{{\rho }_1}{\nabla }_{23}V)=0,$$

$$-({\nabla }_{23}q+\frac{f({\rho }_{1}q)}{{\rho }_1}{\nabla }_{23}V)\cdot n=-a\frac{f({\rho }_{1}q)}{{\rho }_{1}q}+b{\rho }_{1}q.$$

$${\partial }_t(\rho \overline{q})={\partial }_{x_1}(\overline{q}({\partial }_{x_1}{\rho }_1+g({\rho }_1){\partial }_{x_1}V))+A({\rho }_1)-\beta {\rho }_1$$

$$\overline{q}(x_1,t)$$

$$g({\rho }_1(x_1,t))$$

$$A(x_1,\rho (x_1,t),t)$$

$$\beta (x_1,t)$$

$$g({\rho }_1)={\rho }_1,A(x_1,{\rho }_1,t)={\int }_{\partial {\Omega }_2(x_1)}a(x)ds=\alpha (x_1).$$

$${\partial }_t(\rho \overline{q})={\partial }_{x_1}(\overline{q}({\partial }_{x_1}{\rho }_1+{\rho }_1{\partial }_{x_1}V))+\alpha -\beta {\rho }_1.$$

$$A(x_1,\rho (x_1,t),t)={\int }_{\partial {\Omega }_2(x_1)}a(x)(1-{\rho }_{1}q)ds=(a_0-a_1{\rho }_1).$$

$${\partial }_t\rho +\nabla \cdot J=0\text{ with }J=-\nabla \rho +\rho \nabla V,\text{ on }\Omega \times (0,T)$$

$$-J\cdot n$$

$$J\cdot n$$

$$\displaystyle\int_{\Omega}u^{2}~{}dx\leq C_{F}\Big{(}\int_{\Omega}|\nabla u|^{2}~{}dx+\int_{\Gamma}u^{2}~{}d\sigma\Big{)}$$

$$E(\rho |{\rho }_{\infty })=\frac{1}{2}{\int }_{\Omega }\frac{(\rho -{\rho }_{\infty }{)}^2}{{\rho }_{\infty }}dx,$$

$$||\rho -{\rho }_{\infty }|{|}_{L^1(\Omega )}^2\le 2\max \{{\rho }_{\infty }\}E({\rho }_0|{\rho }_{\infty })e^{-C_F{K}_{1}t}\forall t\ge 0,$$

$${\int }_{\Omega }{\partial }_t\rho \phi dx-{\int }_{\Omega }J\cdot \nabla \phi dx+{\int }_{{\Gamma }_{\text{out}}}\beta \rho \phi d\sigma ={\int }_{{\Gamma }_{\text{in}}}\alpha \phi d\sigma$$

$${\int }_{\Omega }{e}^V{\partial }_{t}u\phi dx+{\int }_{\Omega }{e}^V\nabla u\cdot \nabla \phi dx+{\int }_{{\Gamma }_{\text{out}}}\beta u{e}^V\phi d\sigma ={\int }_{{\Gamma }_{\text{in}}}\alpha \phi d\sigma .$$

$$B[u,\phi ]={\int }_{\Omega }{e}^V\nabla u\cdot \nabla \phi dx+{\int }_{{\Gamma }_{\text{out}}}\beta u{e}^V\phi d\sigma .$$

$$\begin{array}{rl}B[u,u]& ={\int }_{\Omega }{e}^V|\nabla u{|}^{2}dx+{\int }_{{\Gamma }_{\text{out}}}\beta e^V|u{|}^{2}d\sigma \\ & \ge K_2{\int }_{\Omega }|\nabla u{|}^{2}dx+K_2{\beta }_0{\int }_{{\Gamma }_{\text{out}}}|u{|}^{2}d\sigma =\eta \Vert u{\Vert }_{H^1(\Omega )}^2-\eta \Vert u{\Vert }_{L^2(\Omega )}^2,\end{array}$$

$${\partial }_t{\int }_{\Omega }{e}^V\frac{|u^{-}{|}^2}{2}dx+{\int }_{\Omega }{e}^V|\nabla u^{-}{|}^{2}dx+{\int }_{{\Gamma }_{\text{out}}}\beta e^V|u^{-}{|}^{2}d\sigma ={\int }_{{\Gamma }_{\text{in}}}\alpha u^{-}d\sigma .$$

$${\int }_{\Omega }{e}^V\frac{|u^{-}{|}^2}{2}dx\le {\int }_{\Omega }{e}^V\frac{|u_0^{-}{|}^2}{2}dx.$$

$$\nabla \cdot \left(-\nabla {\rho }_{\infty }+{\rho }_{\infty }\nabla V\right)=0\text{ on }\Omega ,$$

$${\int }_{\Omega }{e}^V\nabla u_{\infty }\cdot \nabla \phi dx+{\int }_{{\Gamma }_{\text{out}}}\beta u_{\infty }{e}^V\phi d\sigma ={\int }_{{\Gamma }_{\text{in}}}\alpha \phi d\sigma$$

$${\int }_{{\Gamma }_{\text{out}}}\beta u_{\infty }{e}^{V}d\sigma ={\int }_{{\Gamma }_{\text{in}}}\alpha d\sigma .$$

Full text

Keywords: Fokker-Planck Equations, Entropy Methods, Exponential Decay, Mass Evolution

1 Introduction

A lot of recent research effort was devoted to gain understanding of the develoment of neuron cells, which are highly relevant for brain function and disfunction, but still their growth process is far from being completely understood. Takano et al. stated that neurons develop ”structurally and functionally distinct processes called axons and dendrites” [TXF*+*15, p. 1] and while in the beginning all neurons look the same, over time they polarize and generate only one axon but multiple dendrites, see Figure 1. Currently, most mathematical models deal with the influence of proteins in this process, modeled by systems of reaction diffusion equations. Here, we are interested in a different aspect, namly in the transport of vesicles from the cell nucleus to the neurite tips. Our model is based on the Fokker-Planck equation

[TABLE]

where $\rho=\rho(x,t)$ denotes the density of vesicles and $V=V(x)$ is a given potential. The function $f(\rho)$ is chosen to be either $f(\rho)=\rho$ or $f(\rho)=\rho(1-\rho)$. While the first case is just linear transport, the second choice enforces the bound $\rho\leq 1$ on the vesicle density. To model the in- and outflux of vesicles, we will supplement (1) either with flux boundary conditions or reaction terms. For the first approach we divide the boundary of the domain $\Omega$ into three parts: An inflow region, an outflow region and an insulated part (cf. [BP16]). On the inflow part, we prescribe a fixed flux while on the outflow region, the flux is proportional to the density. For the second model, we add no-flux boundary conditions and reaction terms for the in- and outflow of vesicles. This can be seen as an averaging of in- and outflow boundary conditions in a thin higher-dimensional structure.

In all cases, the mass of $\rho$ changes during its evolution making it difficult to exploit the formal gradient flow structure of the equations (cf. [JKO98, Ott01]),

[TABLE]

where $F^{\prime}(\rho)=\frac{1}{f(\rho)}.$ We will partly exploit the underlying gradient flow structure by considering a relative entropy (or Bregman distance) to a stationary solution $\rho_{\infty}$ instead, namely

[TABLE]

After verifying existence and uniqueness of a stationary solutions, we can use the dissipation of the relative entropy to show that solutions of the PDE decay exponentially to the respective stationary or equilibrium solution (we call a stationary solution an equilibrium if the flux vanishes, i.e. $f(\rho_{\infty})\nabla E^{\prime}(\rho_{\infty})=0$).

Entropy methods are very convenient to analyse the long-time behaviour of linear and non-linear partial equations and are strongly related to functional inequalities like the logarithmic-Sobolev inequality, see [MV00]. In case of in- and outflow terms in the bulk we can directly exploit the dissipation generated by the reaction terms in order to show exponential convergence to equilibrium, similar to recent approaches for reaction-diffusion equations (cf. [Mie11, LM13, DF06, FK17, HHMM18]). In the case of in- and outflow boundary condition, exponential convergence needs to be shown using the bulk dissipation by the diffusion and transport. However, standard logarithmic Sobolev inequalities do not apply in our case as the total mass is not conserved. In particular, using a scaling argument, we can give a counter example in our setting. This is in contrast to many similar models in the literature, (see [MV00, ACD*+*04]). However, in the case when the differential operator is linear, we can resort to a variation of Friedrichs’ inequality taking the boundary values into account. This allows us to bound the entropy dissipation in terms of the relative entropy by which, together with a Gronwall-type argument, we recover exponential decay of the relative entropy. Combining this with a Csiszár-Kullback inequality, one also obtains decay in the $L^{1}$-norm.

This paper is organized as follows: In section 2 we explain the biological background of the models. In section 3, we present the linear model where in- and outflow terms are modeled by boundary in- and outflow. In section 4 we investigate a model with spatially distributed in- and outflow and in section 5 we combine this model with a density constraint.

2 Biological Background and Modelling

Neurons are the major part of the central nervous system receiving and transmitting information through the human body. A typical neuron consists of a cell nucleus and two types of ‘arms’ originating in the cell body (see Figure 1).

These arms are called dendrites and axons. Each neuron has several dendrites but only one axon. “An axon is typically a single long process that transmits signals to other neurons by the release of neurotransmitters. Dendrites are composed of multiple branches processes and dendritic spines, which contain neurotransmitter receptors to receive signals from other neurons” [TXF*+*15, p. 1]. The formation of these different processes, called polarization, is crucial for a proper functionality of the central nervous system.

The typical polarization of a neuron in vitro is divided into five stages. In the first stage the neuron extends filopodias around the cell body, which are “thin, actin-rich plasma-membrane protrusions that function as antennae for cells to probe their environment” [ML08, p. 1]. In stage two these filopodia develop into neurites which in the beginning are all equivalent and seem to grow and shrink randomly. The actual polarization starts in stage three where one minor neurite grows quicker than the others and develops into the future axon. In stage four the remaining neurites shrink into dendrites and in stage 5 the polarized neuron matures (see also the poster on neuronal polarization in [TXF*+*15]).

For years biologists have been trying to understand the molecular machinery hidden behind this procedure developing many explanation approaches, see [NMY], which are mainly based on the different concentration of certain proteins in the neurites see [NFN*+*15]. It is conjectured that the neurite growth is mainly driven through a certain vesicle flux. These vesicles are believed to merge with the cell membrane at the neurite tips making the neurite grow. On the other hand it is conjectured that a part of the membrane can also be separated forming a vesicle and making the neuron shrink. The flux of theses vesicles can be measured by special microscopes.

To better understand the growth of the neurite in stage two which seems randomly we suggest a mathematical model describing the vesicle flux in the axons. In the most general case we consider the three-dimensional neurite with several types of different boundary conditions. Thus, in the three-dimensional domain $\Omega_{N}$ modelling the neurite the vesicle density satisfies a Fokker-Planck equation

[TABLE]

with a three-dimensional velocity field $u$. This is complemented by a boundary condition on the flux $\partial\Omega_{N}$ (with $n$ denoting the outward unit normal):

[TABLE]

where the nonnegative functions $a$ and $b$ model rates of in- and outflow, respectively. Typically the supports of $a$ and $b$ do not intersect and we find up to three different regions on the boundary, namely the inflow part where $a$ is positive, the outflow part where $b$ is positive, and the isolated part where $a=b=0$. The influx term in our model therefore corresponds to vesicles entering the neurite, which happens mainly at the part of the boundary that is an interface to the cell nucleus, while outflux corresponds to vesicles merging with the cell membrane, which typically happens at the neurite tip. Note that in its present form, our model does not yet explicitly account for the growth of the neurite, this might however be encoded in the transport terms when rescaling the domain. Moreover, since frequently the directed transport along microtubuli dominates over intracellular fluid transport in neurites, it seems reasonable to assume a potential force $u=-\nabla V$.

If we consider the neurite as an almost axisymmetric structure with small diameter, i.e.,

[TABLE]

for some $\Omega_{2}(x_{1})\subset\mathbb{R}^{2}$ with diameter much smaller than one, we can make further approximations. In particular the equilibration orthogonal to the axis, which we assume to be the $x_{1}$ direction, will be fast, hence

[TABLE]

where $q$ is a stationary solution of

[TABLE]

where $\nabla_{23}$ denotes the gradient with respect to $(x_{2},x_{3})$. On the boundary, $q$ satisfies

[TABLE]

Now, taking an average orthogonal to the axis in the small cross-section $\Omega_{2}(x_{1})$, we find with Gauss’ theorem

[TABLE]

with

[TABLE]

Hence, the in- and outflow boundaries naturally lead to analogous reaction terms in the bulk, which motivates the study of such models as well. We mention that we can obtain a two-dimensional version of the equations in geometries approximating a thin sheet as well.

Let us mention that the models simplify in the special cases of functions $f$ we consider. In the linear case, with $f$ being the identity, we find

[TABLE]

Thus, the resulting equation is simply

[TABLE]

In the crowded case, we still have

[TABLE]

Hence, in both cases, the reaction terms have the same shape as the boundary terms and it is consequently natural to study the following cases:

• •

The linear Fokker-Planck equation $f(\rho)=\rho$ with in- and outflow boundary conditions, which will be the subject of section 3.

• •

The crowded Fokker-Planck equation $f(\rho)=\rho(1-\rho)$ with in- and outflow boundary conditions, which was done in [BP16].

• •

The linear Fokker-Planck equation $f(\rho)=\rho$ with reaction terms of the form $a-b\rho$, which is the subject of section 4.

• •

The crowded Fokker-Planck equation $f(\rho)=\rho(1-\rho)$ with reaction terms of the form $a(1-\rho)-b\rho$, which is the subject of section 5.

3 Linear Model with Boundary In- and Outflux

We start by considering the linear Fokker-Planck equation

[TABLE]

for given $T\geq 0$, $x\in\Omega\subset\mathbb{R}^{n}$, and where $V=V(x)$ is a given potential. The unknown function $\rho=\rho(x,t)$ describes the density of vesicles and we supplement the equation by the flux boundary conditions

[TABLE]

where $n$ denotes the outward normal. We make the following assumptions:

• (A1)

The connected and bounded domain $\Omega\subset\mathbb{R}^{n}$ has Lipschitz boundary $\partial\Omega$.

• (A2)

The potential satisfies $V(x)\in W^{1,\infty}(\Omega)$.

• (A3)

The initial concentration $\rho_{0}$ is non negative and fulfills $\rho_{0}\in L^{2}(\Omega)$.

• (A4)

The subsets $\Gamma_{\text{in}},\Gamma_{\text{out}}\subset\partial\Omega$ of the boundary are open, disjoint and $\Gamma_{\text{out}}$ is nonempty.

• (A5)

The functions $\alpha$ and $\beta$ satisfy $\alpha\in L^{\infty}(\Gamma_{\text{in}})$ with $\alpha\geq\alpha_{0}>0$ and $\beta\in L^{\infty}(\Gamma_{\text{out}})$ with $\beta\geq\beta_{0}>0$ for some $\alpha_{0},\beta_{0}>0$.

In this setting vesicles enter the domain at $\Gamma_{\text{in}}$ with the rate $\alpha$ and leave it at $\Gamma_{\text{out}}$ with rate $\beta\rho$. An additional insulated part of the domain $\partial\Omega\setminus\{\Gamma_{\text{in}}\cup\Gamma_{\text{out}}\}$ may exist, see also Figure 2 for a sketch of such a geometry in 2D. Note that in (10) the term $\beta\rho=0$ is zero if and only if $\rho=0$, which means that as soon as vesicles reach the exit, they can leave the domain with rate $\beta\rho$. In particular due to the boundary conditions (9) - (11) there is no mass conservation, i.e. the spatial integral over $\rho(x,t)$ changes in time.

Remark 3.1**.**

Choosing the influx parameter $\alpha=0$, the stationary solution is obviously $\rho_{\infty}=0$ as particles only leave but never enter the domain. As exponential convergence can be shown for the classical Fokker-Planck equation, as for example done in [MV00], it clearly hold in the ”no influx case“, too. But as the quadratic relative entropy, see below, is undefined in this case, we will only consider $\alpha>0$ from now on.

Exponential convergence of the relative entropy is based on the following version of Friedrichs’ inequality:

Theorem 3.2**.**

[Dud15*, p. 4]**

For a bounded domain $\Omega$ and functions $u\in H^{1}(\Omega)$ the Friedrichs’ inequality with boundary values holds, i.e. for any non empty subset $\Gamma\subseteq\partial\Omega$ of the boundary*

[TABLE]

holds for a constant $C_{F}>0$, that we call Friedrich’s constant.

More precisely, we consider the quadratic relative entropy

[TABLE]

where $\rho_{\infty}$ denotes the solution to the stationary version of (8)–(11) (this will be made precise below). The main result of this section is the following theorem:

Theorem 3.3**.**

Let (A1)-(A5) hold, then the solution $\rho$ to equation (8) with initial datum $\rho_{0}$ and boundary conditions (9) - (11) obeys the following exponential decay towards equilibrium

[TABLE]

where $K_{1}=\min\{\frac{\beta_{0}}{2},1\},\,C_{F}$ being the constant of Friedrich’s inequality and $\rho_{\infty}$ being the stationary solution.

3.1 Existence and Uniqueness for the Time Dependent Problem

As $\rho=\rho(x,t)$ describes the density of vesicles, one naturally expects $\rho\geq 0$, which is proven below. But first we introduce the notion of weak solution and prove an existence result.

Definition 3.4**.**

We say that a function $\rho\in L^{2}(0,T;H^{1}(\Omega))$ with $\partial_{t}\rho\in L^{2}(0,T;H^{-1}(\Omega))$ is a weak solution to equation (8) supplemented with the boundary conditions (9)-(11) if the identity

[TABLE]

holds for all $\varphi\in H^{1}(\Omega)$ and a. e. time $0\leq t\leq T$. We can rewrite the weak formulation in terms of the Slotboom variable $u\coloneqq\rho e^{-V}$. As the transformed flux is given as $J=-e^{V}\nabla u$, we obtain

[TABLE]

Lemma 3.5**.**

Let (A1) - (A5) hold. Then there exists a unique weak solution to equation (8) in the sense of definition 3.4.

Proof.

We use the Slotboom formulation of the problem and use the second and third term on the left side of (12) to define the continuous but non-coercive bilinear form

[TABLE]

This form fulfills a Gårding inequality, see e.g. [Eva10, 6.2.2 Theorem 2], so that for all $u\in H^{1}(\Omega)$ we have

[TABLE]

where $K_{2}=\inf\{e^{V}\},~{}\eta=K_{2}C_{F}^{-1}\min\{1,\beta_{0}\}$ and $C_{F}$ is the constant coming from the Friedrich’s inequality with boundary values, see theorem 3.2. Existence of a unique solution $u$ then follows directly form the ideas stated in [Eva10, 7.1.2] together with the trace theorem [Eva10, 5.5 Theorem 1], applied to the right side of (12). ∎

Lemma 3.6**.**

The weak solution $\rho$ to equation (8) is non-negative, i.e. $\rho\geq 0$ in a.e. $x\in\Omega$ and every $t\in T$.

Proof.

If we use the Slotboom-formulation of the problem, $u$ is non-negative if and only if $\rho$ is non-negative. Choosing the test function $u^{-}=\min\{u,0\}$, the weak formulation yields

[TABLE]

Omitting the non-positive right hand side as well as the non-negative second and the third term on the left hand side, and integration with respect to time gives

[TABLE]

As $u_{0}$ is assumed to be non-negative this yields the assertion. ∎

3.2 Stationary Solutions

We denote by $\rho_{\infty}$ the (weak) solution to

[TABLE]

supplemented with boundary conditions (9)–(11). Our first result is the following.

Lemma 3.7**.**

Any stationary solution $\rho_{\infty}\in H^{1}(\Omega)$ is bounded and strictly positive, i.e. there exists constants $C_{0}>0$ and $C>0$, such that $C_{0}<\rho_{\infty}\leq C$ for a.e. $x\in\Omega$.

Proof.

As $V$ is bounded, $\rho_{\infty}$ is bounded if and only if the stationary Slotboom-variable $u_{\infty}$, which satisfies

[TABLE]

for every $\varphi\in H^{1}(\Omega)$, is bounded. First we choose $\varphi\equiv 1$ in (15) obtaining

[TABLE]

Next we choose $\varphi=(u_{\infty}-C)_{+}$ in (15) for some arbitrary constant $C>0$ and obtain after adding an appropriate trivial term

[TABLE]

Applying Friedrich’s inequality from theorem 3.2 then results in

[TABLE]

with $K=\min\{1,\beta\}$. As $\alpha$ is assumed to be bounded, we can reformulate the right hand side of this inequality by using (weighted) Cauchy’s inequality, see [Eva10, p. 622], with constant $1/2\gamma$ and using the trace inequality for $H^{1}$-functions, see [Eva10, 5.5 Thm.1], to estimate

[TABLE]

Combining (17) and (18) and using (16) we gain

[TABLE]

If we choose $\gamma$ large enough so that $A$ is positive and $C\geq\frac{\alpha}{2\gamma}$ such that $B$ is positive, we can omit $B$ obtain that a.e. $u_{\infty}\leq C$.

To show that the stationary solution is strictly positive, we first note that the argument of Lemma 3.6 also holds for the stationary solution and thus $u_{\infty}\in L^{\infty}(\Omega)$. Together with $u\in H^{1}(\Omega)$, this allows us to apply [LS02, Theorem 4] to (15) to conclude strict positivity. ∎

Next we show that such a stationary solution actually exists and that it is unique, closely following the proof of Proposition 4.1 in [BP16].

Lemma 3.8** (Existence of the stationary solution).**

Let (A1)-(A5) hold, then there exists exactly one stationary solution $\rho_{\infty}\in H^{1}(\Omega)$ to (14) with boundary conditions (9)–(11).

Proof.

Using the Gårding inequality (13), existence follows from the standard theory for elliptic equations, see [Eva10, Section 6.2.]. Now let $\rho_{1}$ and $\rho_{2}$ be two stationary solutions, then $\omega=\rho_{1}-\rho_{2}$ satisfies

[TABLE]

with boundary conditions

[TABLE]

Now let $\nu=\omega e^{-V}$, then $\nu$ is the weak solution of $\nabla\cdot(e^{V}\nabla\nu)=0$ in $\Omega$ with boundary conditions

[TABLE]

Using the weak formulation of this boundary value problem with test function $\nu$ implies

[TABLE]

which yields $\nu=0=\omega$ (using Friedrichs’ inequality with boundary values) and thus uniqueness of the solution. ∎

3.3 Entropy Dissipation with Logarithmic Entropy

The standard approach to prove exponential convergence of $\rho(x,t)$ as $t\to\infty$ would be choosing a logarithmic entropy functional and using a logarithmic-Sobolev inequality to bound the dissipation by the relative entropy. But this approach fails in our case, due to the lack of mass conservation. More precisely, a scaling argument shows that the desired logarithmic-Sobolev inequality fails to hold. Indeed, consider

[TABLE]

Calculating the dissipation of this functional yields

[TABLE]

where we used $u\nabla\log(u)=\nabla u$ in the penultimate line and Gauss’ theorem in the last equation. To show the entropy entropy dissipation inequality we need the following lemma:

Lemma 3.9**.**

The following sum of integrals is positive:

[TABLE]

Proof.

The function $F(t)=t\log(t)$ is convex for $t>0$, so $F(t)-F(s)-F^{\prime}(s)(t-s)\geq 0$. Thus we have $t\log\frac{t}{s}-t+s\geq 0$, which we denote by ($*$) for the future. We want to add the following sum onto $B$ which is zero because of Gauss’s theorem:

[TABLE]

So we gain by addition of zero, ($*$) and $\log(x)-x+1\leq 0$ for $x\geq 0$ the following:

[TABLE]

∎

Applying Lemma 3.9 gives us

[TABLE]

If we define $\varphi=\sqrt{\frac{\rho}{\rho_{\infty}}}$, we would need the following inequality, weighted by the strictly positive function $\rho_{\infty}$,

[TABLE]

to gain the desired result. Surprisingly this inequality is badly scaled. To see this, assume that there exists a function $\varphi$ having trace zero on $\Gamma_{\text{out}}$ that satisfies the inequality. Replacing $\varphi$ by $K\varphi$ with constant $K\geq 0$ gives

[TABLE]

and if we chose $K$ small enough, the inequality becomes wrong.

3.4 Entropy dissipation with Quadratic Entropy

Although the logarithmic entropy is physically more natural, the preceding discussion showed that it is not suitable in our setting. However, as the problem is linear, the quadratic relative entropy defined as

[TABLE]

yields the desired exponential decay towards equilibrium. Its dissipation is given as

[TABLE]

To proceed we test the equation $-\nabla\cdot J_{\infty}=0$ with $\varphi=\frac{\rho-\rho_{\infty}}{\rho_{\infty}}$ which yields

[TABLE]

Adding the last two equations then results in

[TABLE]

where the last inequality holds as $\nabla\cdot J_{\infty}=0,~{}\beta\geq\beta_{0}$ and the influx term is positive. With $K_{1}=\min\{\frac{\beta_{0}}{2},1\}$ and $C_{F}$ being the Friedrich’s constant, we gain

[TABLE]

as $\rho_{\infty}$ is strictly positive as proven in lemma 3.8. Combining this inequality with Gronwall’s lemma, we obtain

[TABLE]

and due to the Cauchy-Schwarz inequality the desired exponential decay in $L^{1}(\Omega)$ via

[TABLE]

3.5 Numerical Solutions

Finally we want to solve the one-dimensional problem numerically. After an appropriate scaling we assume $\Omega=[0,1]$ with $\partial\Omega=\{0,1\}$, so the one-dimensional version of (8) becomes

[TABLE]

with boundary condition

[TABLE]

Note that in this setting we can give an explicit characterization of the stationary solution. As the flux rates $\alpha$ and $\beta$ are constants we have $J=-\partial_{x}\rho_{\infty}+\rho_{\infty}\partial_{x}V=\alpha$ because of the boundary condition at $x=0$. Solving this ordinary differential equation, we obtain the strictly positive stationary solution

[TABLE]

For simplicity, we chose $V(x)=x$ in the following, i.e. the simplest case in which mass is transported from the left entrance to the right exit. We introduce an uniformly spaced grid with $n$ grid points $x_{0},..,x_{n-1}$, with $x_{0}=0$ and $x_{n-1}=1$ and discretize (24) using a fully explicit finite difference scheme. The boundary conditions are implemented by introducing two fictitious nodes $x_{-1}$ and $x_{n}$ outside of the physical domain. The algorithm is implemented in MATLAB, choosing the number of grid points $n=200$ and $dt=5\cdot 10^{-6}$.

3.5.1 Time Evolution of the Density

In Figure 3 the time evolution of the particle concentration $\rho$ solving (8) is compared to the stationary solution $\rho_{\infty}$ for $\alpha=1,\beta=0.9$ and initial concentration $\rho_{0}(x)=-0.1x+1.2$. For our choice $V(x)=x$ the stationary solution is of the form

[TABLE]

In Figure 3 (a) the direct comparison between the initial function $\rho_{0}(x)=-0.1x+1.2$ and the stationary solution is shown. We chose this initial concentration as its shape is very different to the stationary solution, i.e. the gradient has the opposite sign. In Figure 3 (b) the influence of the boundary conditions is strongly visible. At $x=0$ particles have entered the domain and at $x=1$ particles have left the domain. This effect is combined with the drift which leads to a maximum in the left part of the domain and to a minimum at the right. In Figure 3 (c) one sees that after some time the shape of the solution is similar to the one of the stationary solution and only the low frequency parts need longer to converge. Lastly Figure 3 (d) shows this long time behavior and not surprisingly there is no difference visible between $\rho$ and $\rho_{\infty}$ anymore.

3.5.2 Convergence Rates for the Relative Entropy

Next we compare the numerical rate of convergence to the analytical results of section 3.4. To this end we chose $\alpha=\beta=1$ which yields $\rho_{\infty}=1$. Starting again with initial concentration $\rho_{0}(x)=-0.1x+1.2$, we observe exponential convergence with rate $m_{n}=-2.33$ (corresponding to the case $\gamma=1$ in Figure 4).

Now we examine the calculations of section 3.4. Except for the application of Friedrich’s inequality, all other manipulations are equalities. Thus, comparing (22) and (23), we see that the analytic rate of convergence $m_{a}$ is determined by the constant in Friedrich’s inequality, i.e.

[TABLE]

To solve this minimization problem we define the functional

[TABLE]

is the functional Gâteaux-derivative in direction $\psi$ after integration by parts. It is zero if the three conditions

[TABLE]

hold. The function $\varphi(x)=a\sin(kx)+b\cos(kx)$ satisfies the first condition if $\lambda=k^{2}$. The third condition gives $b=2ak$ and for simplicity we chose $a=1$. Thus the second condition translates to

[TABLE]

where a numerical calculation shows that the smallest and positive $k$, which solves this equation is $k\approx 0.9602$. Inserting this in the first condition gives $m_{a}=2\lambda=\inf 2k^{2}\approx 1.8439$. So analytically the relative entropy can be bounded above by $\tilde{f}(x)=E(\rho_{0}|\rho_{\infty})e^{-1.8439t}$.

Surprisingly this value differs significantly from the numerically calculated slope $m_{n}=-2.33$. This can be explained by the fact that the operator

[TABLE]

is symmetric except for the drift term which is skew-symmetric. For any symmetric operator $A$ the spectral gap $\lambda$ determines the slowest possible convergence rate as

[TABLE]

by Gronwall’s lemma. The skew-symmetric part can however mix the eigenvalues which can result in a faster rate of convercence. To examine this phenomena in more detail, we consider only the the symmetric part of the operator (28). Its eigenvalue is determined by

[TABLE]

Elementary calculations yield that the function $u(x)=\cos(kx)$ is an eigenfunction with smallest eigenvalue $\lambda=k^{2}=0.8603^{2}$ where $k$ is the solution of $-k^{2}\cos(kx)+\lambda\cos(kx)=0$. Thus the convergence rate of the symmetric problem is $m_{s}=2\lambda=1.4802$, as calculating the infimum in (27) is equivalent to calculating the eigenvalue of the corresponding operator. Indeed, this is confirmed by our numerical calculations in the case $V=0$.

The entropy dissipation calculation in section 3.4 on the other hand is insensitive to the skew-symmetric part of the operator because we used integration by parts to gain (22) and to get rid of the potential-term. This explains the deviation between $m_{n}$ and $m_{a}$ as soon as $V\neq 0$.

Finally, we analyze numerically how the strength of the potential term influences the convergence rate. We consider

[TABLE]

for different values of $\gamma\leq 0$ and the same boundary conditions as in the initial situation, i.e. (9) - (11). This modification of the PDE effects the shape of the relative entropy and the (non exponential) relation between the scaling factor and the rate of convergence is shown in Figure 4 (c).

Remark 3.10**.**

Clearly $\alpha$ has no influence on the convergence rate as it is not part of the operator (28). It just influences the constant $E(\rho_{0}|\rho_{\infty})$ of the relative entropy as $\alpha$ influences $\rho_{\infty}$. The outflux term $\beta$ influences both the convergence rate and the constant.

Remark 3.11**.**

Note that depending on the choice of initial datum, the mass evolution due to the boundary conditions may dominate over the exponential convergence for short times. Clearly, asymptotically, one observes the exponential rate shown in theorem 3.3. This also motivates the proceeding subsection in which we study the mass evolution in more detail.

3.5.3 Evolution of the Total Mass

Most other works consider the case of an unbounded domain with confining potential or no-flux boundary conditions which yields a preservation of the total mass. This is not true in our case where we have

[TABLE]

as vesicles can enter or leave the domain at $\Gamma_{\text{in}}$ or $\Gamma_{\text{out}}$. We want to examine this evolution numerically using two different initial conditions. As $\rho$ converges exponentially to some equilibrium state, its mass also converges exponentially to the mass of the stationary solution. Yet in contrast to the relative entropy, the evolution of the mass is not monotone in time. To shed light on this phenomena we are now looking for initial functions which enforce non monotone mass evolution. In Figure 5 (a) the initial function is

[TABLE]

The mass of this initial function is about 1.1863 which is more than the mass of the equilibrium state 1.0703, so in the long run the mass evolution should be a decreasing function. But at the beginning the mass of $\rho$ is increasing in time. This can be explained by the fact that particles are pumped into the domain with rate $\alpha$ at $x=0$ whereas no particles can leave the domain as they are no particles at the exit of the domain at $x=1$. Thus the mass increases until the drift term has transported (a substantial number of) particles to $x=1$. In Figure 5 (b) the initial function is

[TABLE]

which yields the opposite behavior. At first the mass of the particles concentration is decreasing and then it increases. This is due to the fact that the particle concentration at $x=1$ is 30 which is very much compared to the rest of the domain. As the outflux is proportional to the concentration at $x=1$ there are more particles leaving the domain than entering it. The mass of the initial function is 1.0711 and the mass of the equilibrium state is 1.0696. Note that one can see that the maximum of the mass of the initial function and the mass of the stationary solution is not an upper bound for the mass of $\rho$.

4 Linear Model with Spatially Distributed In- and Outflux

In our second model we consider the Fokker-Planck equation where influx and outflux is modeled by reaction terms, i.e.

[TABLE]

for $t\geq 0$ and $x\in\mathbb{R}^{n}$, $\alpha\in L^{\infty}(\Omega)$ with $\alpha\geq 0,~{}\beta\in L^{\infty}(\Omega)$ with $\beta\geq\beta_{0}>0$ for some $\beta_{0}>0$ and $V=V(x)$ is a smooth and bounded potential as in the previous section. Furthermore we assume a no flux condition i.e.

[TABLE]

We make the following assumptions:

• (B1)

The domain $\Omega\subset\mathbb{R}^{n}$ is a bounded with $\partial\Omega\in C^{1,1}$.

• (B2)

The initial function is non-negative and fulfills $\rho_{0}\in L^{2}(\Omega)$.

• (B3)

The potential $V$ is smooth and bounded, i.e. there are constants $K_{2}=\inf\{e^{-V}\}$ and $K_{3}=\sup\{e^{-V}\}$. Furthermore $\nabla V$ and $\Delta V$ are also bounded and $-\Delta V-\beta e^{V}\leq 0$ on $\Omega$.

• (B4)

The functions $\alpha$ and $\beta$ satisfy $\alpha\in L^{\infty}(\Omega)$ with $\alpha\geq\alpha_{0}>0$ and $\beta\in L^{\infty}(\Omega)$ with $\beta\geq\beta_{0}>0$ for some $\alpha_{0},\beta_{0}>0$.

The reaction terms in (31) translate to vesicles entering and leaving the domain at every point with rates $\alpha$ and $\beta\rho e^{-V}$, respectively. In this case, the stationary solution is given by $\rho_{\infty}=\frac{\alpha}{\beta}e^{V}$. The following theorem, which we will prove in 4.3, is the main result of this section:

Theorem 4.1**.**

Let (B1)-(B3) hold, then equation (31) with no flux boundary conditions obeys the following exponential decay towards equilibrium

[TABLE]

where $K_{1}$ depends on $\mathrm{max}\{||\rho_{\infty}||_{\infty},||\rho_{0}||_{\infty}\}$ only and $E(\rho|\rho_{\infty})$ is the logarithmic relative entropy defined in (19). Furthermore $K_{2}=\inf\{e^{-V}\}$ and $K_{4}$ come from the Cziszár-Kullback-Pinsker inequality in Lemma 4.8.

4.1 Existence of the Time Dependent Problem

Definition 4.2**.**

A function $\rho\in L^{2}(0,T;H^{1}(\Omega))$ with $\partial_{t}\rho\in L^{2}(0,T;H^{-1}(\Omega))$ is a weak solution to equation (31) supplemented with the boundary condition (32) if the identity

[TABLE]

holds for all $\varphi\in H^{1}(\Omega)$ and a. e. time $0\leq t\leq T$. Again we can rewrite the weak formulation in terms of the Slotboom variable $u\coloneqq\rho e^{-V}$ and obtain

[TABLE]

Lemma 4.3**.**

If $\nabla V\in L^{\infty}(\Omega)$ and $\rho_{0}\in L^{2}(\Omega)$, there exists a unique weak solution to equation (31) in the sense of definition 33.

Proof.

We use the Slotboom formulation, as in (34), which yields the bilinear form

[TABLE]

that obviously fulfills [Eva10, 6.2.2 Theorem 2] as $V$ is bounded and $\beta\geq\beta_{0}>0$. Next we apply the ideas stated in [Eva10, pp. 353 - 358]. ∎

4.2 Stationary Solutions

Lemma 4.4**.**

For $V,\,\nabla V$ and $\Delta V$ bounded and $-\Delta V-\beta e^{V}\leq 0$ on $\Omega$, the unique stationary solution $\rho_{\infty}\in H^{2}(\Omega)$ to

[TABLE]

with boundary conditions (32) is given by

[TABLE]

Proof.

Rewriting the problem as

[TABLE]

we see that our assumption on $V$ guarantees a sign on the lower order term and the assertion follows, e.g. from [Gri85, Thm. 2.4.2.6]. ∎

Remark 4.5**.**

If we choose $\beta=\alpha=0$, we lose uniqueness of the solution of (31). But if we demand the solution to be positive, we can restore uniqueness see [DV09, Theorem 1.1].

Having proven the uniqueness of the stationary solution, enables us to give an upper bound for $\rho$, which we need for the equilibrium property. i.e:

Lemma 4.6**.**

If (31) and (32) hold, $\rho\geq 0$ and there is an upper bound $L$ for $\rho$ i.e. for every $t\in(0,T)$ there holds

[TABLE]

Proof.

We define $\tilde{L}\coloneqq\mathrm{max}\{\alpha\beta^{-1},||\rho_{0}e^{-V}||_{\infty}\}$, which shall be an upper bound for $u$ reminding ourselves that $\rho_{0}e^{-V}=u_{0}$. As $e^{-V}\partial_{t}u=\partial_{t}\rho$ and $\tilde{L}$ is a constant we can write

[TABLE]

Using the weak formulation of this equation with test function $(u-\tilde{L})_{+}\in H^{1}(\Omega)$, which has the derivative

[TABLE]

we obtain

[TABLE]

by integration by parts, using $J\cdot n=(e^{-V}\nabla u)\cdot n=0$ on $\partial\Omega$. As the second integral is positive, we can omit it to achieve

[TABLE]

Now we use $\frac{e^{-V}}{2}\geq C$ as $V$ is bounded, and Gronwall’s lemma yields

[TABLE]

so $u\leq\tilde{L}$ a.e. and therefore $\rho\leq\tilde{L}\sup\{e^{-V}\}=L$ a.e. in $\Omega$ and for every $t\in(0,T)$. Repeating the same argument with test function $u^{-}=\min\{u,0\}$ gives $0\leq\rho$. ∎

4.3 Long Time Behavior

Our proof of the exponential decay to equilibrium is closely related to ideas presented in [DF06]. A central tool in this paper which is used to give an upper bound on the relative entropy is the following lemma which we state for the sake of completeness:

Lemma 4.7**.**

[DF06*, Lemma 2.1]**

The function*

[TABLE]

is continuous on $(0,+\infty)^{2}$. For all $x>0$ the function $\varphi(x,\cdot)$ is strictly decreasing and respectively for all $y>0$ the function $\varphi(\cdot,y)$ is strictly increasing on $(0,+\infty)$. Finally it satisfies

[TABLE]

With this at hand, we proceed to the proof of theorem 4.1:

Proof of theorem 4.1: The previous lemma enables us to rewrite the logarithmic entropy (19) as

[TABLE]

where $K_{1}\coloneqq\max\{1,\varphi(L,0)\}>0$. Note that we just take the maximum with $1$ to ensure that $K_{1}$ is positive which will make further computations more easy. Furthermore the entropy dissipation $D(\rho|\rho_{\infty})$ is given by

[TABLE]

and by integration by parts with no flux boundary conditions we get

[TABLE]

where the penultimate inequality holds because of the elementary inequality $(a-b)(\log(a)-\log(b))\geq 4(\sqrt{a}-\sqrt{b})^{2}$. As $V$ is bounded, there exists a constant $K_{2}>0$ with $\inf\{e^{-V}\}=K_{2}$. Combining the estimates for the relative entropy and the entropy dissipation, we achieve

[TABLE]

Combining this with Gronwall’s lemma, we obtain

[TABLE]

To conclude the proof we use the following lemma which is a generalization of the Cziszar-Kullback-Pinsker inequality for functions which are not probability densities. $\qquad\square$

Lemma 4.8**.**

[HHMM18*, Lemma A.1]**

Let $\Omega$ be a measurable domain in $\mathbb{R}^{d}$. Let $\rho,\rho_{\infty}\colon\Omega\rightarrow\mathbb{R}^{+}$ be measurable. Then,*

[TABLE]

Because of lemma 4.6 the constant $K_{4}$ in lemma 4.8 is finite. Combining inequality (37) with (38) we conclude the desired exponential convergence, i.e.,

[TABLE]

Remark 4.9**.**

Combining the results of this section with the dissipation of the logarithmic entropy in section 3.3 it seems natural that one can also conclude the exponential convergence to equilibrium in the case of the Fokker-Planck equation with both in- and outflow in the bulk and on the boundary, i.e.

[TABLE]

with boundary conditions (9)-(11). However, this is not possible with the current techniques for reaction-diffusion equations that we used above, since they require that the stationary solution is of the form $\rho=ce^{-V}$, which is not the case with the flow boundary conditions unless $\alpha=\beta=0$.

Remark 4.10**.**

Replacing $\nabla V$ by a vector field $\vec{u}$ that is not necessarily the gradient of some vector field $V$ would be an interesting generalization of this model.

4.4 Numerical Solution

The following results are again based on a finite difference scheme with explicit time discretization, see section 3.5 for details.

4.4.1 Time Evolution of the Density

In Figure 7 one can see the time evolution of the concentration $\rho$ solving (31) in comparison to the calculated stationary solution $\rho_{\infty}=\frac{\alpha}{\beta}e^{V(x)}$ for $\alpha=1,~{}\beta=0.9$ and initial concentration $\rho_{0}(x)=-0.1x+1.2$. In contrast to the previous section there is no flow direction determined by the in- and outflux terms as particles enter the domain uniform in space.

4.4.2 Convergence Rates for the Relative Entropy

The numerical solution of the logarithm of the relative entropy (35) for the one dimensional version of (31) for $\alpha=1,~{}\beta=0.9$ and initial concentration $\rho(x)=-0.1x+1.2$ can be seen in Figure (8). For the data given above the relative entropy has roughly the shape of $0,22e^{-1.04t}$ in the interval $t\in[0,15]$. After that machine precision is reached. Here we did not use the explicit stationary solution $\rho_{\infty}=\alpha/\beta e^{x}$ as the reference value in the relative entropy but we calculated the stationary solution numerically up to machine precision.

5 Nonlinear Model with Spatially Distributed In- and Outflux

As vesicles have a positive volume, there naturally exists a maximal number that can fit into an axon. This motivates the next generalization of the Fokker-Planck equation given as

[TABLE]

with the same properties for $\rho,\alpha$ and $\beta$ as in the previous sections. The additional nonlinear term $\rho(1-\rho)$ and the modification of the inflow rate by $(1-\rho)$ ensure that the density stays within $[0,1]$ for all times, where $1$ corresponds to the scaled maximal density of vesicles. Again we assume no flux boundary conditions $J\cdot n=0$ on $\partial\Omega\times(0,T)$ and the following three assumptions:

• (C1)

The domain $\Omega\subset\mathbb{R}^{n}$ is a connected and bounded with $\partial\Omega\in C^{1,1}$.

• (C2)

The initial condition satisfies $\rho_{0}\in W^{2-2/p,p}(\Omega)$ for some fixed $2<p<3$ and the box constraints $0<\rho_{0}<1$.

• (C3)

The potential $V$ is smooth and bounded with $\nabla V\in L^{\infty}(\Omega)$.

Again we want to show exponential decay to equilibrium which is surprisingly simple although in contrast to the previous settings (40) is not a linear problem anymore. We chose the logarithmic entropy

[TABLE]

for equation (40) where $h(\rho)$ denotes the entropy density. We obtain the corresponding relative entropy

[TABLE]

motivated by viewing $\rho$ and $1-\rho$ as two different types of species. Our aim is the following result, which we will prove in section 5.3:

Theorem 5.1**.**

Let (C1)-(C3) and (B4) hold. Then every weak solution to equation (40) with no flux boundary conditions obeys the following exponential decay towards equilibrium:

[TABLE]

with $\tilde{C}=\alpha\min\{1,\inf\frac{1-\rho_{\infty}}{\rho_{\infty}}\}$.

5.1 Existence of the Time Dependent Problem

Definition 5.2**.**

We say that a function $\rho\in L^{2}(0,T;H^{1}(\Omega))$ with $\partial_{t}\rho\in L^{2}(0,T;H^{-1}(\Omega))$ is a weak solution to equation (40), supplemented with the boundary condition $J\cdot n=0$ on $\partial\Omega$, if the identity

[TABLE]

holds for all $\varphi\in H^{1}(\Omega)$ and a.e. $0\leq t\leq T$.

Lemma 5.3**.**

Let (C1)-(C3) and (B4) hold, then for every $\rho_{0}\in L^{1}(\Omega)$ with $0\leq\rho_{0}\leq 1$, there exists a weak solution to equation (40) in the sense of (42) which satisfies $0\leq\rho\leq 1$ a.e. in $\Omega$. With

[TABLE]

being the entropy variable, this solution fulfills the dissipation inequality

[TABLE]

Proof.

This proof is based one implicit Euler discretization, following e.g. [BP16, GSW18], to which we refer for more details. We start by rewriting (40), using the entropy variable $u$ and exploiting its formal gradient flow structure, as

[TABLE]

Now fix $N\in\mathbb{N}$ and consider a discretization of $(0,T]$ into subintervals $(0,T]=\cup_{k=1}^{N}[(k-1)\tau,k\tau]$ with time steps $\tau=\frac{T}{N}$, we obtain the following sequence of elliptic problems

[TABLE]

The existence of a solution $\rho_{k+1}$ (given $\rho_{k}$) to the nonlinear equation (44) can be proven via a fixed point argument, see [BP16, Theorem 3.5]. In particular, using the transformation $\rho_{k+1}=h^{\prime-1}(u_{k+1})$ enforces the bounds $0\leq\rho_{k+1}\leq 1$ (sometimes called boundedness by entropy).

In order to be able to pass to the limit $\tau\to 0$, we use the discrete entropy dissipation to get a priori bounds. As the entropy density $h(\rho)=\rho\log\rho-(1-\rho)\log(1-\rho)-\rho V$ is strictly convex for $\rho\in S^{0}$, with $S^{0}$ being the interior of $S=\{\rho\in\mathbb{R}~{}|~{}0\leq\rho\leq 1\}$ we obtain

[TABLE]

Now taking $\varphi=u_{k}\in H^{1}(\Omega)$ as test function in (44), and using $\rho_{k}=h^{\prime-1}(u_{k})$, we obtain the discrete entropy dissipation given as

[TABLE]

Solving the recursion then yields

[TABLE]

To pass to the limit $\tau\rightarrow 0$ we denote by $\rho_{k}$ a sequence of solutions to (45). We define $\rho_{\tau}(x,t)=\rho_{k}(x)$ for $x\in\Omega$ and $t\in((k-1)\tau,k\tau)]$. Then for $\tau\leq t\leq T$, the function $\rho_{\tau}$ solves the following problem

[TABLE]

where $\sigma_{\tau}$ denotes the shift operator, that is $(\sigma_{\tau}\rho_{\tau})(x,t)=\rho_{\tau}(x,t-\tau)$ and for all test functions $\varphi\in L^{2}(0,T;H^{1}(\Omega))$. Next the entropy dissipation inequality (46) becomes

[TABLE]

Following [GSW18, Appendix, Lemma 1], there exists a constant $C$ such that

[TABLE]

which gives the a-priori estimate $\|\rho_{\tau}\|_{L^{2}(0,T;H^{1}(\Omega))}\leq K$ when combined with (48). Thus, upon extraction of a subsequence, $\rho_{\tau}$ converges strongly in $L^{2}(\Omega)$. Together with the weak convergence of $\nabla u_{\tau}$, this is enough to pass to the limit in (46) in all terms but the first one. There we have to apply a special version of the Aubin-Lions lemma for piece-wise constant interpolations [DJ12, Thm 1] which allows us to take $\tau\to 0$. Finally, taking the limit in (48) yields the desired entropy dissipation inequality. ∎

5.2 Stationary Solution

Lemma 5.4**.**

There exists exactly one stationary solution $\rho_{\infty}\in H^{1}(\Omega)$ of equation (40) with no flux boundary conditions given by

[TABLE]

Proof.

Uniqueness of the stationary solution is a direct consequence of Theorem 5.1. Indeed, assuming that there are two different stationary solutions $\rho_{\infty}$ and $\tilde{\rho}_{\infty}$, inserting them into (50) yields

[TABLE]

As the left hand side of the inequality is a constant whereas the right side is a decreasing function in $t$, we obtain a contradiction for $t$ large enough. ∎

5.3 Long Time Behaviour

First we rewrite the reactions terms in equation (40) as

[TABLE]

Next using the analogue of (43) for the relative entropy, we see that the entropy dissipation is given by

[TABLE]

Neglecting the first non-negative part and introducing the function

[TABLE]

we can further estimate the dissipation from below by

[TABLE]

where we used the nonnegativity of the relative entropy and with $\tilde{C}=\alpha\min\{1,\inf\frac{1-\rho_{\infty}}{\rho_{\infty}}\}$. With Gronwall’s lemma and the Cziszár-Kullback-Pinsker inequality in lemma 4.8, we finally achieve

[TABLE]

5.4 Numerical Solution

Even though we are now dealing with a nonlinear equation, we again use the fully explicit scheme of section 3.5 and obtain the following results.

5.4.1 Analysis of the Time Evolution

In Figure 9 we show the time evolution of $\rho(x,t)$ solving (40) compared to the solution (49) with $\alpha=1,\beta=0.9$, initial concentration $\rho_{0}(x)=-(x-0.5)^{2}+1$, potential $V(x)=x$ and $\Omega=[0,1]$. We chose this particular initial particle concentration (see Figure 9 (a)) as it has a completely different shape compared to the stationary solution and secondly has the value 1 at $x=0.5$, so that at one point of the domain the density constraint of 1 is reached. We did not chose the same initial function as in the two previous sections as this initial function does not fulfill the box constraint. In Figure 9 (b) and (c) the effect of the diffusion becomes visible as it has flatten the particle concentration and the effect of the drift effect as there are more particles at the right part of the domain than in the left part. Finally in Figure 9 (d) there is no difference between the stationary solution and the concentration visible. In comparison to the results of the previous model one can see that the density constraint of 1 is never overstepped.

5.4.2 Analysis of the Relative Entropy

In Figure 10 the relative entropy for the one dimensional version of (40) for $\alpha=1,\beta=0.9$, initial concentration $\rho_{0}(x)=-(x-0.5)^{2}+1$ and potential $V(x)=x$ can be seen. To better compare the results with the one of the previous section we chose the same $\alpha$ and the same $\beta$. Comparing the convergence velocity of this model with the previous one, one can see that the case without a density constraint obeys a quicker exponential decay. This can be explained by the following intuition: Figuratively the relative entropy measures the distance to an equilibrium state and the convergence rate how quick this status is reached. In this setting the influx and the drift term are multiplied by the factor $(1-\rho)$ whereas in the previous section it was not, so they have less influence than in the previous model.

Acknowledgments: The authors acknowledge support by EXC 1003 Cells in Motion Cluster of Excellence, Münster, funded by the German science foundation DFG. MB was further supported by ERC via Grant EU FP7 - ERC Consolidator Grant 615216 LifeInverse. The authors would like to thank Andreas Püschel and Danila di Meo (WWU Münster) for details on the biological background.