A simple reaction–diffusion system as a possible model for the origin of chemotaxis

Abstract

Chemotaxis is a directed cell movement in response to external chemical stimuli. In this paper, we propose a simple model for the origin of chemotaxis – namely how a directed movement in response to an external chemical signal may occur based on purely reaction–diffusion equations reflecting inner working of the cells. The model is inspired by the well-studied role of the rho-GTPase Cdc42 regulator of cell polarity, in particular in yeast cells. We analyse several versions of the model to better understand its analytic properties and prove global regularity in one and two dimensions. Using computer simulations, we demonstrate that in the framework of this model, at least in certain parameter regimes, the speed of the directed movement appears to be proportional to the size of the gradient of signalling chemical. This coincides with the form of the chemical drift in the most studied mean field model of chemotaxis, the Keller–Segel equation.

Keywords:

• Chemotaxis
• Reaction–diffusion system
• Cdc42
• Gradient tracking
• Regularity

1. Introduction

Chemotaxis is a directed cell movement in response to external chemical stimuli. Chemotaxis is ubiquitous in biology; for example, it plays a role in organism morphology [19, 20], reproduction processes [8, 17, 18, 23] and workings of immune system [4, 21]. There are many mathematical models of chemotaxis, the most studied is the Keller–Segel equation and its variants. Virtually all of these models incorporate a transport term driven by the concentration of the external chemical, which may be produced by the cells themselves (see, e.g. [11, 15], where further references can be found). Yet we are not aware of any mathematical models that would aim to explain the origin of the transport based on reaction–diffusion processes taking place inside cells. The way chemotaxis happens, at least for eukaryotic cells, is that cells translate chemical environmental cues into amplified intracellular signalling that results in elongated cell shape, actin polymerization toward the leading edge and movement along the gradient. In this paper, instead of presenting chemotaxis as an explicit transport term, we explore model that aims to explain the origin of the chemotactic ability of cells. Inspired by [3], we look at sexual yeast reproduction and simplify the polarization process into understanding active rho-GTPase Cdc42 concentration in one yeast under chemical gradient produced by another yeast. This is certainly just an element of a more complex picture involved into producing chemotactic response in cells, but we limit consideration to this one stage. Our first goal is to explore the well-posedness properties of the model and its variants and to understand analytic features involved. Our second goal is to get more information on the nature of transport generated by the model reacting to external chemical stimuli. In particular, a natural question is how the speed of transport, which we measure via the coordinate moment of a density, depends on the gradient of the attractive chemical. This question we approach through numerical simulations and find that for certain reasonable ranges of parameters, this dependence is linear.

The two-species mass-conserved activator-substrate (MCAS) system (Figure 1) that is the basis of our model consists of two partial differential equations (PDEs) governing the kinetics of the slowly diffusing activator u (GTP-bound GTPase on the membrane) and the rapidly diffusing substrate v (GDP-bound GTPase in the cytoplasm). In general, this system has the following form in 1D (see [3]): (1) $\begin{array}{rl}\frac{\mathrm{\partial }u}{\mathrm{\partial }t}& =F(u,v)+k\mathrm{\Delta }u,\\ \frac{\mathrm{\partial }v}{\mathrm{\partial }t}& =-F(u,v)+k_v\mathrm{\Delta }v.\end{array}$(1) Here, k refers to the diffusion of u, $k_v$ refers to the diffusion of v. These two diffusion constants usually differ by two orders of magnitude. $F(u,v)$ describes the biochemical interconversions between u and v, given in the form: $F(u,v)=h(u)v-g(u)u.$

Figure 1. Local activation via positive feedback and depletion of the substrate in the cytosol generates an activator-enriched domain on the cortex [3].

The positive feedback (i.e. conversion from inactive GDPase to active GTPase with energy) is denoted by $f(u)v$, while the negative feedback (i.e. conversion from active GTPase to inactive GDPase without energy) is denoted by $g(u)u$. Examples of $F(u,v)$ include: the simplest $h(u)=\tilde{a}{u}^2,g(u)=b;$Goryachev's (see [6]) (2) $h(u)=\tilde{a}{u}^2+cu,g(u)=b;$(2) and (see [14]) $h(u)=1,g(u)=\frac{b}{1+u^2}.$In [3], Turing-type instability for these types of reaction have been analysed; it was also shown that steady states with more than one peak are unstable for many kinds of $F(u,v)$. This analysis is in agreement with experiments as one usually only observes one bud in yeast asexual production [3].

Several complicated computational models have been developed to mimic gradient-induced polarization toward the pheromone source [9, 22] and have shown the rate of movement is dependent on pheromone concentration [1]. Here, we propose a simpler system that is capable of capturing such gradient tracking ability, specifically, in the context of chemotactic reaction of a single yeast cell to an external pheromone signal. We apply a modification to the Turing-type model described above and add a pheromone density profile $\tilde{\alpha }f(x)$ that depends on location in the form similar to [22] -- and we obtain the following system: (3) $\begin{array}{rl}\frac{\mathrm{\partial }u}{\mathrm{\partial }t}& =(\tilde{a}{u}^2+\tilde{\alpha }f(x))v-bu+k\mathrm{\Delta }u,\\ \frac{\mathrm{\partial }v}{\mathrm{\partial }t}& =-(\tilde{a}{u}^2+\tilde{\alpha }f(x))v+bu+k_v\mathrm{\Delta }v.\end{array}$(3) In (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }u}{\mathrm{\partial }t}& =(\tilde{a}{u}^2+\tilde{\alpha }f(x))v-bu+k\mathrm{\Delta }u,\\ \frac{\mathrm{\partial }v}{\mathrm{\partial }t}& =-(\tilde{a}{u}^2+\tilde{\alpha }f(x))v+bu+k_v\mathrm{\Delta }v.\end{array}$(3) ), $\tilde{a},b,\tilde{\alpha },k,k_v$ are constants. $\tilde{a}$ is the reaction activation constant, b is the reaction depletion constant, $\tilde{\alpha }$ is the overall pheromone strength, k is the diffusion coefficient for u, and $k_v$ is the diffusion coefficient for v. $f(x)$ is a bounded smooth function that describes the pheromone level at different locations.

We assume that rho-GDPase diffuses infinitely fast, i.e. $k_v$ approaches ∞. Since the total mass of rho-GTPase and rho-GDPase is conserved, $M=\int (u(\cdot ,t)+v(\cdot ,t))\mathrm{d}x$ is a constant. Then we can obtain the following equation (4(4) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(\tilde{a}{u}^2+\tilde{\alpha }f(x))\frac{1}{|{\mathbf{T}}^d|}(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-bu+k\mathrm{\Delta }u.$(4) ) that describes the activator–substrate reaction between these two substances. The setting we have is $x\in {\mathbf{T}}^d$ when dimension d = 1, 2, with periodic boundary condition: (4) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(\tilde{a}{u}^2+\tilde{\alpha }f(x))\frac{1}{|{\mathbf{T}}^d|}(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-bu+k\mathrm{\Delta }u.$(4) In (4(4) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(\tilde{a}{u}^2+\tilde{\alpha }f(x))\frac{1}{|{\mathbf{T}}^d|}(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-bu+k\mathrm{\Delta }u.$(4) ), $|{\mathbf{T}}^d|$ is the measure of the domain, M is the total mass. We are interested in the non-negative solution u with $\int u\mathrm{d}x\le M$ for all time. By rescaling space and time, we can simplify Equation (4(4) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(\tilde{a}{u}^2+\tilde{\alpha }f(x))\frac{1}{|{\mathbf{T}}^d|}(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-bu+k\mathrm{\Delta }u.$(4) ) as follows: (5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) Here depletion rate and diffusion coefficient are normalized to 1, and $\frac{1}{|{\mathbf{T}}^d|}$ gets absorbed into $\tilde{a}$ and $\tilde{\alpha }$ (Figure 2).

Figure 2. The interconversion of Rho-GTPases between active and inactive forms can be modelled as a reaction–diffusion equation governing the dynamics of the slowly-diffusing activator u and the infinitely-diffusing substrate v.

Our main results are as follows. On the rigorous level, we are able to establish global regularity results for Equation (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ) in one and two dimensions for all non-negative initial data. To better understand the structure of the equation, we consider (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ) in the absence of regularizing diffusion and prove that for all finite $0\le t\le T\le \mathrm{\infty }$, the profile for active rho-GTPase u stays smooth even when there is no diffusion term. When we re-introduce diffusion in one dimension, uniform in time global bound on derivatives of u is shown. With diffusion in two dimension, global regularity with possible growth is proved. In our numerical experiments, we observe the profile of active rho-GTPase u move towards higher concentration of pheromone, and stops moving once it reaches the location with maximum pheromone concentration. In addition, we explore the speed of such movement through tracking the centre of mass of rho-GTPase profile. If pheromone concentration is linear, the centre of mass moves with a constant speed towards the pheromone peak. More importantly, we verify that the movement speed depends linearly on the pheromone gradient in a natural parameter range similar to that used in [3]. Note that such linear dependence of chemotactic drift on the gradient of the attractive chemical density $f(x)$ is a common feature of chemotaxis models, including the most studied Keller–Segel equation which in its simplest form reads (see, e.g. [15]) (6) $\mathrm{\partial }\rho -\mathrm{\Delta }\rho +\alpha \mathrm{\nabla }\cdot (\rho \mathrm{\nabla }f)=0.$(6) The emergence of transport mean field equations such as (6(6) $\mathrm{\partial }\rho -\mathrm{\Delta }\rho +\alpha \mathrm{\nabla }\cdot (\rho \mathrm{\nabla }f)=0.$(6) ) from kinetic equations has been extensively studied (see, e.g. [10, 12, 13, 16]). However the existence of chemotactic transport is already built in on the kinetic level. As far as we know, Equation (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ) is the first simple reaction–diffusion model that aims to rigorously analyse the emergence of chemotaxis from the inner cell workings, even if it is focused on just one stage of the process that can be quite complex.

This paper is organized as follows: in the following section, we introduce the general set up and key parameters of the model in more detail and present our numerical scheme. We then proceed to describe results of the numerical experiments. After this we state the rigorous results that we are able to prove and proceed with the proofs.

2. General set up and numerical scheme

We want to explore the origin of chemotactic ability of cells with simulation in $1D$ using the following equation (dropping the scripts in (4(4) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(\tilde{a}{u}^2+\tilde{\alpha }f(x))\frac{1}{|{\mathbf{T}}^d|}(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-bu+k\mathrm{\Delta }u.$(4) ) and denote the total mass of $u(x,t)$ as $U(t):={\int }_{{\mathbf{T}}^d}u(x,t)\mathrm{d}x$): (7) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))\frac{1}{|{\mathbf{T}}^d|}(M-U(t))-bu+k\mathrm{\Delta }u.$(7) The parameters we use are the same as in [3] and are shown in Table 1 with some basic conversions. We used the method of lines to turn spatially discretized PDE into a system of ODEs, then we use a robust ODE solver ODE15s in Matlab to solve. Note that since we assume rho-GDPase v is rapidly diffusing, we use Simpson's method to numerically integrate rho-GTPase to obtain $\int u(x,t)\mathrm{d}x$ and calculate rho-GDPase as follows: (8) $v(t)=\frac{1}{|{\mathbf{T}}^d|}(M-U(t)).$(8) For the computational part, we restrict ourselves to one-dimensional surface and assume the pheromone profile is generated by another yeast cell. If this cell is some distance away, one reasonable model of the two-dimensional pheromone distribution is a solution to the heat equation ${\mathrm{\partial }}_t\omega -\mathrm{\Delta }\omega =0$ with a δ function initial data, that is, as a fundamental solution of 2D heat equation. Then we can derive the pheromone profile on the cell boundary as shown in Figure 3, and in general, it has the form: (9) $f_h(x)=\frac{\beta }{4\pi \gamma t}\left(\exp (-\frac{1}{4\gamma t}(L^2+r^2-2Lr\cos \frac{x-x_{\mathrm{peak}}}{r}))\right),$(9) where $\varphi =\frac{x-x_{\mathrm{peak}}}{r}$ and $\varphi \in [-\pi ,\pi )$. We use γ to denote the diffusion coefficient for the source, β as the response strength to the source, and without loss of generality, we assume t = 1. With $\gamma =10$, t = 1, L = 10, $\beta =1500$, we can plot $f_h(x)$ in Figure. 4. We also plot a linear pheromone profile $f(x)$ that is similar to $f_h(x)$ with a peak at $x_{\mathrm{peak}}$ defined below in Figure 4 as well. The reason for this choice will be explained later. (10) $f(x)=\{\begin{array}{ll}\frac{2}{5}x,& \mathrm{if}0\le x\le 5,\\ 4-\frac{2}{5}x,& \mathrm{if}5\le x\le 10.\end{array}$(10)

Figure 3. Derivation of a pheromone profile generated from a heat equation.

Figure 4. Pheromone profile generated from a heat equation $f_h(x)$ and a similar linear pheromone profile $f(x)$.

Table 1. Parameters from [3].

Parameter	Description	Value/Range in simulation	Unit
k	Diffusion coefficient for u	0.01	$(\text{μ}\mathrm{m}{)}^2({\mathrm{s}}^{-1})$
α	Pheromone strength	0.5−3	${\mathrm{s}}^{-1}$
$|{\mathbf{T}}^d|$	Cell size	10	$\text{μ}\mathrm{m}$
M	Total mass	10	–
a	Strength of activation	1	$(\text{μ}\mathrm{m}{)}^2$
b	Strength of substrate	1	${\mathrm{s}}^{-1}$

To obtain the initial profile of u, we start with a uniformly distributed v and a bump function u, and we run the simulation without pheromone until it stabilizes according to (11) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=a{u}^2\frac{1}{|{\mathbf{T}}^d|}(M-U(t))-bu+k\mathrm{\Delta }u.$(11) Switching on the pheromone, we track the movement using centre of mass since the profile of u is relatively stable over time. The centre of mass as a function of time is defined using (12) ${\mathrm{CM}}_u(t)=\frac{\int xu(x,t)\mathrm{d}x}{U(t)}.$(12) We also measure the movement of the profile of u with the time derivative of the centre of mass: $\frac{\mathrm{d}C{M}_u(t)}{\mathrm{d}t}$.

3. Numerical results: pheromone induced movement

While there is no explicit transport term in (7(7) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))\frac{1}{|{\mathbf{T}}^d|}(M-U(t))-bu+k\mathrm{\Delta }u.$(7) ), we observe movement of the rho-GTPase u profile over time in response to ‘reallocation of resources’ to more favourable reaction regions induced by pheromone $\alpha f_h(x)$ and $\alpha f(x)$ as shown in Figure 5. As our numerical simulations show, at least in certain parameter ranges, this transport appears quite similar to the Keller–Segel-type transport with speed proportional to the concentration gradient.

Figure 5. Numeric solution to Equation (7(7) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))\frac{1}{|{\mathbf{T}}^d|}(M-U(t))-bu+k\mathrm{\Delta }u.$(7) ) with α = 2 and other parameters given in Table 1: (a) pheromone profile given by $f_h(x)$ and (b) pheromone profile given by $f(x)$.

As we can see in Figure 5, one does not expect that the solution will be an exact traveling pulse since the background level of the pheromone affects the shape of the bump of u density. With time t and recorded $C{M}_u(t)$, we can compute the movement speed of the centre of mass, $\frac{\mathrm{d}C{M}_u(t)}{\mathrm{d}t}$ and the corresponding pheromone profile $f_h$ and f at the centre of mass. From Figures 6 and 7, we propose the hypothesis that the movement speed $\frac{\mathrm{d}C{M}_u(t)}{\mathrm{d}t}$ depends on the derivative of the pheromone.

Figure 6. Movement speed $\frac{{\mathrm{CM}}_u(t)}{\mathrm{d}t}$ as a function of pheromone profile slope $\frac{f{(}_h{\mathrm{CM}}_u(t))}{\mathrm{d}t}$. $\frac{\mathrm{d}{f}_h(x)}{\mathrm{d}x}>0$ for $x\in (0,2.5)$.

Figure 7. Movement speed $\frac{{\mathrm{CM}}_u(t)}{\mathrm{d}t}$ as a function of pheromone profile slope $\frac{f({\mathrm{CM}}_u(t))}{\mathrm{d}t}$. $\frac{\mathrm{d}f(x)}{\mathrm{d}x}$ is a constant for $x\in (0,2.5)$.

We have tried some other pheromone profiles with similar results. While the graphs on Figure 7 appear close to lines, they are not quite lines – but perhaps because the density bump of u has spatial scale, and so is exposed to a range of concentration slopes (we take the slope at the centre of mass as the basis of the functional relationship pictured in this figure). As one of our goals is to study the dependence of the movement speed of the centre of mass $C{M}_u(t)$ and on the gradient of the pheromone concentration $\alpha f^{\prime }(x)$, we will also consider the piecewise linear pheromone profile $f(x)$: it has extended regions with the constant slope that makes it easier to capture the effect more precisely. We illustrate this transport picture in Figure 8 by calculating the centre of mass using (12(12) ${\mathrm{CM}}_u(t)=\frac{\int xu(x,t)\mathrm{d}x}{U(t)}.$(12) ) of the initial profile of u and the profile at $t=\mathrm{10,000}\mathrm{s}$ pheromone profile $f(x)$ with pheromone strength $\alpha =2$.

Figure 8. Initial profile for rho-GTPase u, rho-GDPase v, and pheromone profile $f(x)$ and pheromone strength $\alpha =2$: (a) initial profile t = 0s and (b) profile when $t=10,000s$.

As one can expect, the profile of u slows down once its centre of mass starts to approach $x=x_{\mathrm{peak}}=5\text{μ}\mathrm{m}$. We can plot the centre of mass a function as a time of time for different pheromone strength α. As presented in Figure 9, the centre of mass of u stays at $x=x_{\mathrm{peak}}=5\text{μ}\mathrm{m}$ after $t=7000\mathrm{s}$ for the pheromone strength $\alpha =3$. In fact, if we run the simulation long enough, the centre of mass of u appears to get arbitrarily close to $x=x_{\mathrm{peak}}=5\text{μ}\mathrm{m}$ for all $\alpha >0$.

Figure 9. Centre of mass position ${\mathrm{CM}}_u(t)$ as a function of time with pheromone profile $f(x)$. In regions away from $x_{\mathrm{peak}}$, ${\mathrm{CM}}_u(t)$ changes linearly with time.

We then continue to explore the movement speed of the centre of mass as a function of pheromone strength α, which controls the slope of the pheromone concentration. From Figure 9, we can see a constant movement speed of the centre of mass when the profile of u is far away from $x_{\mathrm{peak}}$. Moreover, if we plot the movement speed (before the profile of u is too close to $x_{\mathrm{peak}}$) as a function of α as shown in Figure 10, we observe linear dependence of transport speed on pheromone strength. Such linear dependence corresponds to the mean field chemotaxis models in [11, 15].

Figure 10. Centre of mass movement speed $\frac{d{\mathrm{CM}}_u(t)}{\mathrm{d}t}$ as a function of pheromone strength α with pheromone profile $f(x)$. In regions away from $x_{\mathrm{peak}}$, $\frac{d{\mathrm{CM}}_u(t)}{\mathrm{d}t}$ changes linearly with α.

4. Mathematical analysis: global regularity

Our goal in this section is to initiate rigorous mathematical analysis of Equation (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ). There is much literature on regularity of solutions for reaction–diffusion type equations, but we could not locate references dealing directly with (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ) due to nonlocality produced by our assumption of infinite diffusion for v. In [2], the author presented global existence result to a similar equation for sufficiently small and non-negative initial data. In this section, we will present regularity result for (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ) in ${\mathbf{T}}^d$, with d = 1, 2 for any non-negative smooth initial data with diffusion and without diffusion. When there is no diffusion with d = 1, 2, we can still show that u is smooth for all times $0\le t<\mathrm{\infty }$. With diffusion in 1D, uniform in time global bound is proven while with diffusion in 2D, global regularity with possible growth is shown. At this time, we are unable to prove rigorous results that would provide more qualitative features of the evolution, and in particular establish the conjecture on linear dependence of the transport speed on pheromone gradient at least in some regimes. The part of the difficulty is that, as we mentioned before, the solution cannot be expected to take form of a traveling pulse of the fixed shape on an unbounded domain, and lack of such framework complicate analysis. This paper can be viewed as creating a foundation for such further investigation that remains an intriguing open problem.

4.1. Without diffusion

We start with the case without diffusion. Consider the following equation: (13) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-U(t))-u,U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(13)

Theorem 4.1

Suppose $u(t,x)$ is a non-negative solution to (13(13) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-U(t))-u,U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(13) ) with dimension d = 1, 2 and periodic boundary condition; $a,\alpha ,M$ are constant parameters, and $f(x)$ is a smooth non-negative function. If $u_0(x)$ is a smooth initial profile, then $u(t,x)$ stays smooth for all times $0\le t<\mathrm{\infty }$.

Proof.

We first show some a-priori bounds on u, namely that all Sobolev norms of u are controlled by $L^{\mathrm{\infty }}$ norm. It is clear that on a bounded domain, $L^{\mathrm{\infty }}$ norm controls all other $L^p$ norms. In the estimates that follow, D stands for any partial derivative (just ${\mathrm{\partial }}_x$ in one dimension), and $W^{k,p}$ is the usual Sobolev space. Multiplying (13(13) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-U(t))-u,U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(13) ) by $(-\mathrm{\Delta }{)}^{s}u$ and integrating, we get (14) ${\mathrm{\partial }}_t||u|{|}_{H^s}^2\le C\underset{I}{\underbrace{|\int u^2\cdot (-\mathrm{\Delta }{)}^{s}udx|}}-\underset{||u|{|}_{H^s}^2}{\underbrace{\int u(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x}}+\alpha M|\int f(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x|.$(14) The last term can by controlled by $\alpha M||f|{|}_{W^{2s,1}}||u|{|}_{L^{\mathrm{\infty }}}$ (moving all derivatives to f). Term I can be represented by a sum of integrals of the type $\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x$, where $l=0,\dots ,s$. Then with Hölder's inequality and Gagliardo–Nirenberg interpolation inequality in dimension d = 1, 2, we can bound them by $\begin{array}{rl}\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x& \le C||D^{l}u|{|}_p||D^{s-1}u|{|}_q||D^{s}u|{|}_2\mathrm{with}\frac{1}{p}+\frac{1}{q}+\frac{1}{2}=1;\\ ||D^{l}u|{|}_p& \le C||D^{s}u|{|}_2^{\alpha }||u|{|}_{\mathrm{\infty }}^{1-\alpha };\\ ||D^{s-l}u|{|}_q& \le C||D^{s}u|{|}_2^{\beta }||u|{|}_{\mathrm{\infty }}^{1-\beta }.\end{array}$In $1D$, $\alpha =\frac{2(-1+lp)}{p(-1+2s)}$ and $\beta =\frac{2(-1+(s-l))q}{q(-1+2s)}$, and $\alpha +\beta =2(\frac{1}{p}+\frac{1}{q})=1$.

In $2D$, $\alpha =\frac{-2+lp}{p(-1+s)}$ and $\beta =\frac{-2+(s-l)q}{q(-1+s)}$, and $\alpha +\beta =2(\frac{1}{p}+\frac{1}{q})=1$. Then $\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x\le C||D^{l}u|{|}_p||D^{s-1}u|{|}_q||D^{s}u|{|}_2\le C||u|{|}_{H^s}^2||u|{|}_{\mathrm{\infty }}.$Substituting back into (14(14) ${\mathrm{\partial }}_t||u|{|}_{H^s}^2\le C\underset{I}{\underbrace{|\int u^2\cdot (-\mathrm{\Delta }{)}^{s}udx|}}-\underset{||u|{|}_{H^s}^2}{\underbrace{\int u(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x}}+\alpha M|\int f(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x|.$(14) ) gives ${\mathrm{\partial }}_t||u|{|}_{H^s}^2\le (C||u|{|}_{\mathrm{\infty }}-1)||u|{|}_{H^s}^2+\alpha M||f|{|}_{W^{2s,1}}||u|{|}_{L^{\mathrm{\infty }}}.$Applying Gronwall inequality [7], we see that to show global regularity, it suffices to prove that ${\int }_0^T||u(\cdot ,t)|{|}_{\mathrm{\infty }}\mathrm{d}t$ remains bounded. We will show a stronger constraint that $\Vert u(\cdot ,t){\Vert }_{\mathrm{\infty }}$ remains finite for all times. Via contradiction, denote T the first time of blow up of $||u|{|}_{\mathrm{\infty }}$.

Consider $\rho =\frac{e^{-t}}{1+u}$, then ρ satisfies the following equation: (15) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-(a{e}^{-t}-2a\rho +a{\rho }^2{e}^t+\alpha f(x){\rho }^2{e}^t)(M-U(t))-{\rho }^2{e}^t.$(15) Then at time T, ρ will reach 0 at some point. For simplicity, let us focus on d = 1; the argument for d = 2 is similar. From (15(15) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-(a{e}^{-t}-2a\rho +a{\rho }^2{e}^t+\alpha f(x){\rho }^2{e}^t)(M-U(t))-{\rho }^2{e}^t.$(15) ), we see that ρ is decreasing for all times, therefore ρ is bounded from above by $||{\rho }_0|{|}_{\mathrm{\infty }}$. Now we take a derivative with respect to x to obtain ${\mathrm{\partial }}_x{\mathrm{\partial }}_t\rho =-(M-U(t))(-2a{\mathrm{\partial }}_x\rho +2a{e}^t\rho {\mathrm{\partial }}_x\rho +\alpha {\rho }^2{e}^t{\mathrm{\partial }}_{x}f+2\alpha \rho e^{t}f(x){\mathrm{\partial }}_x\rho )-e^t\rho {\mathrm{\partial }}_x\rho .$Then it is clear that for all x, ${\mathrm{\partial }}_t|{\mathrm{\partial }}_x\rho |\le C_0(T)+C_1(T)|{\mathrm{\partial }}_x\rho |.$By Gronwall inequality, we can see that $|{\mathrm{\partial }}_x\rho |$ is bounded for any finite time $t\le T$. We can take another derivative in x and obtain $\begin{array}{rl}{\mathrm{\partial }}_{xx}{\mathrm{\partial }}_t\rho =& -(M-U(t))(-2a{\mathrm{\partial }}_{xx}\rho +2a{e}^t\rho {\mathrm{\partial }}_{xx}\rho +2a{e}^t({\mathrm{\partial }}_x\rho {)}^2+\alpha {\rho }^2{e}^t{\mathrm{\partial }}_{xx}f+\alpha \rho e^t{\mathrm{\partial }}_{x}f{\mathrm{\partial }}_x\rho \\ & +2\alpha \rho e^{t}f(x){\mathrm{\partial }}_{xx}\rho +2\alpha e^{t}f(x)({\mathrm{\partial }}_x\rho {)}^2+2\alpha \rho e^t{\mathrm{\partial }}_{x}f{\mathrm{\partial }}_x\rho )-e^t\rho {\mathrm{\partial }}_{xx}\rho -e^t({\mathrm{\partial }}_x\rho {)}^2.\end{array}$The boundedness of ρ along with control of $|{\mathrm{\partial }}_x\rho |$ gives ${\mathrm{\partial }}_t|{\mathrm{\partial }}_{xx}^2\rho |\le C_2(T)+C_3(T)|{\mathrm{\partial }}_{xx}\rho |.$By Gronwall inequality again, one gets that $|{\mathrm{\partial }}_{xx}^2\rho |$ is bounded for any finite time $t\le T$. One can effectively continue this calculation and get that all derivatives in space are bounded for $t\le T$. Since blow up happens for the first time at time T, then $\rho (x_B,T)=0$ at some point $x_B$. There can be many such points, but let us focus on one of them. Due to control of ${\mathrm{\partial }}_x|\rho |$ and ${\mathrm{\partial }}_{xx}^2|\rho |$, we have $\rho (x,T)\le C(T)|x-x_B{|}^2$ by Taylor expansion. Observe that $\rho (x,t)\to \rho (x,T)$ monotonically for every x. Therefore, as $u=\frac{1}{e^t\rho }-1$, we have $u(x,t)\to u(x,T)$ (including when $u(x,t)=\mathrm{\infty }$). Then we have $u(x,T)=\frac{1}{e^T\rho }-1\ge \frac{1}{e^T|x-x_B{|}^2}-1$. Then by Fatou's lemma, we have (16) $\begin{array}{rl}M& \ge \underset{t->T}{lim inf}\int u(x,t)\mathrm{d}x\ge \int \underset{t->T}{lim inf}u(x,t)\mathrm{d}x=\int u(x,T)\mathrm{d}x\\ & \ge \int C|x-x_B{|}^{-2}\mathrm{d}x=\mathrm{\infty },\end{array}$(16) which is a contradiction. Therefore, we cannot have such finite time blow up. Note that the argument above works both in 1D and 2D, only the computation yielding control of the derivatives of ρ needs a minor adjustment. Note that the size of the Sobolev norms of the solution may depend on time.

4.2. With diffusion in 1D

Now we turn out attention to the system with diffusion in $1D$: (17) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-U(t))-u+\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2},U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(17) We will prove global regularity for (17(17) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-U(t))-u+\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2},U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(17) ) as well.

Theorem 4.2

Suppose u is a non-negative solution to (17(17) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-U(t))-u+\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2},U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(17) ) in dimension d = 1 and periodic boundary condition. Let $a,\alpha ,M$ be constant parameters, and $f(x)$ a smooth function. If $u_0(x)$ is a smooth initial profile, then $u(x,t)$ stays smooth for all time; in particular, all Sobolev norms $||u|{|}_{H^s}$ with s>0 are bounded uniformly for all time.

Proof.

First we show that $||u|{|}_2$ is bounded: multiplying both sides by u and integrating, we have $\int u\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\mathrm{d}x=\int u(a{u}^2+f(x))(M-U(t))\mathrm{d}x-\int u^2\mathrm{d}x+\int u\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}\mathrm{d}x.$Therefore, (18) $\frac{1}{2}{\mathrm{\partial }}_t\int u^2\mathrm{d}x\le M\int a{u}^3\mathrm{d}x-\int u^2\mathrm{d}x-\int (u_x{)}^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}.$(18) Using Gagliardo–Nirenberg–Sobolev inequality (see, e.g. [5]) gives (19) $||u|{|}_{L^3}\le C||u|{|}_{L^1}^{5/9}||u_x|{|}_{L^2}^{4/9}.$(19) Substituting (19(19) $||u|{|}_{L^3}\le C||u|{|}_{L^1}^{5/9}||u_x|{|}_{L^2}^{4/9}.$(19) ) into (18(18) $\frac{1}{2}{\mathrm{\partial }}_t\int u^2\mathrm{d}x\le M\int a{u}^3\mathrm{d}x-\int u^2\mathrm{d}x-\int (u_x{)}^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}.$(18) ) yields (note that the constant C changes from line to line and may depend on M): (20) $\begin{array}{rl}\frac{1}{2}{\mathrm{\partial }}_t\int u^2\mathrm{d}x& \le C{(\int u\mathrm{d}x)}^{5/3}{(\int u_x^2\mathrm{d}x)}^{2/3}-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le \frac{C}{3}{(\int u\mathrm{d}x)}^5+\frac{2}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le -\frac{1}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}+C(M).\end{array}$(20) Note that in the second inequality above, we used Young's inequality $ab\le \frac{a^p}{p}+\frac{b^q}{q}$ with $p=3,q=\frac{3}{2}$. The calculation (20(20) $\begin{array}{rl}\frac{1}{2}{\mathrm{\partial }}_t\int u^2\mathrm{d}x& \le C{(\int u\mathrm{d}x)}^{5/3}{(\int u_x^2\mathrm{d}x)}^{2/3}-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le \frac{C}{3}{(\int u\mathrm{d}x)}^5+\frac{2}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le -\frac{1}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}+C(M).\end{array}$(20) ) implies that $||u|{|}_{L^2}$ is globally bounded by Gronwall inequality. Moreover, from (20(20) $\begin{array}{rl}\frac{1}{2}{\mathrm{\partial }}_t\int u^2\mathrm{d}x& \le C{(\int u\mathrm{d}x)}^{5/3}{(\int u_x^2\mathrm{d}x)}^{2/3}-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le \frac{C}{3}{(\int u\mathrm{d}x)}^5+\frac{2}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le -\frac{1}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}+C(M).\end{array}$(20) ), we can see that $||u|{|}_2$ is in fact uniformly bounded by $M^2||f|{|}_{\mathrm{\infty }}+C(M)$ since if $||u|{|}_2$ ever crosses this value for the first time at $t_0$, ${\mathrm{\partial }}_t||u|{|}^2$ becomes negative, which implies that before $t_0$, $||u|{|}_2$ lies above $M^2||f|{|}_{\mathrm{\infty }}+C{M}^5$ already. We arrive at a contradiction. Therefore, $||u|{|}_2$ is uniformly bounded for all time.

Next, we estimate the higher order Sobolev norms. In fact, in 1D, this can be done using just control of the $L_1$ norm, but we are going to use $L_2$ norm for convenience as we have shown it remains bounded. Multiplying (5(5) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}=(a{u}^2+\alpha f(x))(M-{\int }_{{\mathbf{T}}^d}u\mathrm{d}x)-u+\mathrm{\Delta }u.$(5) ) by $(-\mathrm{\Delta }{)}^{s}u$ and integrating, we get (21) ${\mathrm{\partial }}_t||u|{|}_{H^s}^2\le C\underset{I}{\underbrace{|\int u^2\cdot (-\mathrm{\Delta }{)}^{s}udx|}}-\underset{||u|{|}_{H^s}^2}{\underbrace{\int u(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x}}-||u|{|}_{H^{s+1}}^2+\alpha M|\int f(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x|.$(21) The last term can by controlled by $\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2}$ (moving all derivatives to f). Term I can be represented by a sum of integrals of the type $\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x$, where $l=0,\dots ,s$. In 1D, $D={\mathrm{\partial }}_x$. Then with Hölder's inequality and Gagliardo–Nirenberg interpolation inequality with dimension d = 1, we can bound them by (22) $\begin{array}{rl}\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x& \le C||D^{l}u|{|}_p||D^{s-1}u|{|}_q||D^{s}u|{|}_2\mathrm{with}\frac{1}{p}+\frac{1}{q}+\frac{1}{2}=1\\ & \le C||u|{|}_{H^{s+1}}^{2-\frac{3}{2(s+1)}}||u|{|}_{L^2}^{1+\frac{3}{2(s+1)}}.\end{array}$(22) Here we deploy the Gagliardo Nirenberg inequalities $||D^{s}u|{|}_{L^2}=||u|{|}_{H^s}\le ||u|{|}_{L^2}^{\frac{1}{s+1}}||u|{|}_{H^{s+1}}^{\frac{s}{s+1}}$, $||D^{l}u|{|}_p\le C||u|{|}_{L^2}^{1-\alpha }||D^{s+1}u|{|}_{L^2}^{\alpha }$, with $\alpha =\frac{l-\frac{1}{p}+\frac{1}{2}}{s+1}$, and $||D^{s-l}u|{|}_q\le C||u|{|}_{L^2}^{1-\beta }||D^{s+1}u|{|}_{L^2}^{\beta }$, with $\beta =\frac{s-l-\frac{1}{q}+\frac{1}{2}}{s+1}$. Substituting gives (note that the constant C changes from line to line): (23) $\begin{array}{rl}{\mathrm{\partial }}_t||u|{|}_{H^s}^2& \le C||u|{|}_{H^{s+1}}^{2-\frac{3}{2(s+1)}}||u|{|}_{L^2}^{1+\frac{3}{2(s+1)}}-||u|{|}_{H^s}^2-||u|{|}_{H^{s+1}}^2+\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2}\\ & \le \frac{4s+1}{4s+4}||u|{|}_{H^{s+1}}^2-||u|{|}_{H^{s+1}}^2+\frac{3C}{4s+4}||u|{|}_{L^2}^{\frac{4s+10}{3}}-||u|{|}_{H^s}^2+\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2}\\ & \le -||u|{|}_{H^s}^2+\frac{3C}{4s+4}||u|{|}_{L^2}^{\frac{4s+10}{3}}+\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2},\end{array}$(23) where we use Young's inequality $ab\le \frac{a^p}{p}+\frac{b^q}{q}$ with $p=\frac{4s+4}{4s+1}$, and $q=\frac{4s+4}{3}$. Given that we proved $\Vert u{\Vert }_2$ is bounded uniformly in time, (23(23) $\begin{array}{rl}{\mathrm{\partial }}_t||u|{|}_{H^s}^2& \le C||u|{|}_{H^{s+1}}^{2-\frac{3}{2(s+1)}}||u|{|}_{L^2}^{1+\frac{3}{2(s+1)}}-||u|{|}_{H^s}^2-||u|{|}_{H^{s+1}}^2+\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2}\\ & \le \frac{4s+1}{4s+4}||u|{|}_{H^{s+1}}^2-||u|{|}_{H^{s+1}}^2+\frac{3C}{4s+4}||u|{|}_{L^2}^{\frac{4s+10}{3}}-||u|{|}_{H^s}^2+\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2}\\ & \le -||u|{|}_{H^s}^2+\frac{3C}{4s+4}||u|{|}_{L^2}^{\frac{4s+10}{3}}+\alpha M||f|{|}_{H^{2s}}||u|{|}_{L^2},\end{array}$(23) ) implies that $||u|{|}_{H^s}$ is also bounded uniformly for all time.

4.3. With diffusion in 2D

The equation that we are interested is given by (24) ${\mathrm{\partial }}_{t}u=a{u}^2(M-U(t))+\alpha f(x)(M-U(t))-u+\mathrm{\Delta }u,U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(24)

Theorem 4.3

Suppose u is a non-negative solution to (24(24) ${\mathrm{\partial }}_{t}u=a{u}^2(M-U(t))+\alpha f(x)(M-U(t))-u+\mathrm{\Delta }u,U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(24) ) with dimension d = 2 and periodic boundary condition. Let $a,\alpha ,M$ be constant parameters, and $f(x)$ a smooth function. If $u_0(x)$ is a smooth initial profile, then $u(x,t)$ stays smooth for any finite time, that is, Sobolev norms $||u|{|}_{H^s}$ with s>0 are bounded for any time $0\le t<\mathrm{\infty }$. The bound on the Sobolev norms may now depend on time.

Proof.

First we derive an a-priori estimate. In the 2D case, the analog of the estimate (20(20) $\begin{array}{rl}\frac{1}{2}{\mathrm{\partial }}_t\int u^2\mathrm{d}x& \le C{(\int u\mathrm{d}x)}^{5/3}{(\int u_x^2\mathrm{d}x)}^{2/3}-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le \frac{C}{3}{(\int u\mathrm{d}x)}^5+\frac{2}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x-\int u_x^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}\\ & \le -\frac{1}{3}(\int u_x^2\mathrm{d}x)-\int u^2\mathrm{d}x+M^2||f|{|}_{\mathrm{\infty }}+C(M).\end{array}$(20) ) is not available, as the exponents do not allow to control $\Vert u{\Vert }_{L^2}$ uniformly in time in this way. Therefore, we need a more nuanced argument. Note that (25) ${\int }_0^T\int {\mathrm{\partial }}_{t}u\mathrm{d}x\mathrm{d}t={\int }_0^T\int (a{u}^2(M-U(t))+\alpha f(x)(M-U(t))-u+\mathrm{\Delta }u)\mathrm{d}x\mathrm{d}t$(25) gives (in the following calculations, the constant C changes from line to line): (26) ${\int }_0^T\int (a{u}^2(M-U(t)))\mathrm{d}x\mathrm{d}t\le U(T)+MT.$(26) Then multiplying (24(24) ${\mathrm{\partial }}_{t}u=a{u}^2(M-U(t))+\alpha f(x)(M-U(t))-u+\mathrm{\Delta }u,U(t)={\int }_{{\mathbf{T}}^d}u\mathrm{d}x.$(24) ) by $(-\mathrm{\Delta }{)}^{s}u$ and integrating in x, we obtain (27) ${\mathrm{\partial }}_t||u|{|}_{H^s}^2\le -||u|{|}_{H^{s+1}}^2-||u|{|}_{H^s}^2+a(M-\int u\mathrm{d}x)\underset{I}{\underbrace{\int u^2(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x}}+\alpha M||f|{|}_{H^s}||u|{|}_{H^s}.$(27) The integral I is a sum of terms of the form: $\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x$, and $l=0,\dots ,s$. We estimate it as follows: (28) $\begin{array}{rl}\int D^{l}u{D}^{s-l}u{D}^{s}u\mathrm{d}x& \le C||u|{|}_{H^s}||D^{l}u|{|}_p||D^{s-l}|{|}_q\mathrm{with}\frac{1}{p}+\frac{1}{q}+\frac{1}{2}=1\\ & \le C||u|{|}_{H^s}||u|{|}_{H^{s+1}}^{\alpha }||u|{|}_{L^2}^{1-\alpha }||u|{|}_{H^{s+1}}^{\beta }||u|{|}_{L^2}^{1-\beta }\\ & \le C||u|{|}_{L^2}||u|{|}_{H^s}||u|{|}_{H^{s+1}}.\end{array}$(28) Here we use Gagliardo Nirenberg inequalities for n = 2: $||D^{l}u|{|}_p\le C|u|{|}_{H^{s+1}}^{\alpha }||u|{|}_{L^2}^{1-\alpha }$ with $\alpha =\frac{l+1-\frac{2}{p}}{s+1}$, and $||D^{s-l}|{|}_q\le C||u|{|}_{H^{s+1}}^{\beta }||u|{|}_{L^2}^{1-\beta }$ with $\beta =\frac{s-l+1-\frac{2}{q}}{s+1}$, $\alpha +\beta =1$. Substituting these estimates back into (27(27) ${\mathrm{\partial }}_t||u|{|}_{H^s}^2\le -||u|{|}_{H^{s+1}}^2-||u|{|}_{H^s}^2+a(M-\int u\mathrm{d}x)\underset{I}{\underbrace{\int u^2(-\mathrm{\Delta }{)}^{s}u\mathrm{d}x}}+\alpha M||f|{|}_{H^s}||u|{|}_{H^s}.$(27) ) gives (29) $\begin{array}{rl}{\mathrm{\partial }}_t||u|{|}_{H^s}^2& \le -||u|{|}_{H^{s+1}}^2-||u|{|}_{H^s}^2+Ca(M-U(t))||u|{|}_{L^2}||u|{|}_{H^s}||u|{|}_{H^{s+1}}\\ & +\alpha M||f|{|}_{H^s}||u|{|}_{H^s}\\ & \le CM{a}^2(M-U(t))||u|{|}_{L^2}^2||u|{|}_{H^s}^2+\frac{1}{2}||u|{|}_{H^{s+1}}^2-||u|{|}_{H^{s+1}}^2\\ & -||u|{|}_{H^s}^2+\alpha M||f|{|}_{H^s}||u|{|}_{H^s}\\ & \le CM{a}^2(M-U(t))||u|{|}_{L^2}^2||u|{|}_{H^s}^2-\frac{1}{2}||u|{|}_{H^{s+1}}^2+\alpha M||f|{|}_{H^s}||u|{|}_{H^s},\end{array}$(29) where in the second line, we used the inequality: $2ab\le ϵa^2+\frac{1}{ϵ}{b}^2$. Then by Gronwall inequality, we get $||u|{|}_{H^s}\le \exp (CMa{\int }_0^T\int a{u}^2(M-U(t))\mathrm{d}x\mathrm{d}t)(||u(0)|{|}_{H^s}+{\int }_0^T\alpha M||f|{|}_{H^s}\mathrm{d}t).$From (26(26) ${\int }_0^T\int (a{u}^2(M-U(t)))\mathrm{d}x\mathrm{d}t\le U(T)+MT.$(26) ), we see that $||u|{|}_{H^s}$ is bounded for any finite $t\le T<\mathrm{\infty }$.

Acknowledgments

AK has been partially supported by the NSF-DMS award 2006372. The authors thank Daniel Lew for patiently teaching them some of the relevant biology (all inadequacies are our fault) and helpful discussions. We are grateful to anonymous referees for useful comments that helped improve the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Science Foundation  Division of Mathematical Sciences [2006372].