BOUNDARY LAYERS AND STABILIZATION OF THE SINGULAR KELLER-SEGEL SYSTEM

. The original Keller-Segel system proposed in [23] remains poorly understood in many aspects due to the logarithmic singularity. As the chemical consumption rate is linear, the singular Keller-Segel model can be converted, via the Cole-Hopf transformation, into a system of viscous conservation laws without singularity. However the chemical diﬀusion rate parameter ε now plays a dual role in the transformed system by acting as the coeﬃcients of both diﬀusion and nonlinear convection. In this paper, we ﬁrst consider the dynamics of the transformed Keller-Segel system in a bounded interval with time-dependent Dirichlet boundary conditions. By imposing appropriate conditions on the boundary data, we show that boundary layer proﬁles are present as ε → 0 and large-time proﬁles of solutions are determined by the boundary data. We employ weighted energy estimates with the “eﬀective viscous ﬂux” technique to establish the uniform-in- ε estimates to show the emergence of boundary layer proﬁles. For asymptotic dynamics of solutions, we develop a new idea by exploring the convexity of an entropy expansion to get the basic L 1 estimate. We the obtain the corresponding results for the original Keller-Segel system by reversing the Cole-Hopf transformation. Numerical simulations are performed to interpret our analytical results and their implications.


(Communicated by Tao Luo)
Abstract. The original Keller-Segel system proposed in [23] remains poorly understood in many aspects due to the logarithmic singularity. As the chemical consumption rate is linear, the singular Keller-Segel model can be converted, via the Cole-Hopf transformation, into a system of viscous conservation laws without singularity. However the chemical diffusion rate parameter ε now plays a dual role in the transformed system by acting as the coefficients of both diffusion and nonlinear convection. In this paper, we first consider the dynamics of the transformed Keller-Segel system in a bounded interval with time-dependent Dirichlet boundary conditions. By imposing appropriate conditions on the boundary data, we show that boundary layer profiles are present as ε → 0 and large-time profiles of solutions are determined by the boundary data. We employ weighted energy estimates with the "effective viscous flux" technique to establish the uniform-in-ε estimates to show the emergence of boundary layer profiles. For asymptotic dynamics of solutions, we develop a new idea by exploring the convexity of an entropy expansion to get the basic L 1estimate. We the obtain the corresponding results for the original Keller-Segel system by reversing the Cole-Hopf transformation. Numerical simulations are performed to interpret our analytical results and their implications.
1. Introduction. The oriented movement of species up/down to the chemical concentration gradient is termed as chemotaxis which has been a significant mechanism to interpret abundant pattern formation and biological processes such as bacteria band formation and aggregation [38,48], slime mould formation [16], fish pigmentation patterning [41], angiogenesis in tumor progression [6,7,8], primitive streak formation [42], blood vessel formation [14], wound healing [44], and so on. Proposed by , the chemotaxis model has two prototypes according to the chemotactic sensitivity function. One was the linear sensitivity and the other was the logarithmic sensitivity. The former was derived in [25] to model the selfaggregation of Dictyostelium discoideum in response to cyclic adenosine monophosphate (cAMP), and the latter in [23] to model the wave propagation of bacterial chemotaxis. Compared to massive results on the Keller-Segel (KS) model with linear sensitivity, much less is known on the KS model with logarithmic sensitivity due to its singularity nature. However logarithmic sensitivity complies with the Webber-Fecher law and has many prominent applications in biology (cf. [2,3,10,22,23]) in addition to its indispensable role to reproducing the bacterial traveling bands (cf. [49]). This paper is concerned with the original KS model proposed in [23] where u(x, t) and w(x, t) denote the bacterial density and concentration of nutrient (chemical), respectively, at position x and time t. The parameter D > 0 is the diffusivity of bacterial, χ > 0 is referred to as the chemotactic coefficient measuring the intensity of chemotaxis, ε ≥ 0 is the chemical diffusion rate and m ≥ 0 is the consumption rate of nutrient. It has been shown (cf. [24,46,49]) that the KS model (1) will produce traveling bands (pulsating waves) if 0 ≤ m < 1, and fronts if m = 1 and no traveling waves if m > 1, where the logarithmic sensitivity is indispensable to generate traveling waves. In the case of 0 ≤ m < 1, the KS model (1) was employed by Keller and Segel to interpret the bacterial traveling band formation observed in the experiment by Adler [1]. When m = 1, (1) was first used by Nossal [40] to model the boundary movement of bacterial and later by Levine et al [27] to model the dynamics between vascular endothelial growth factor (VEGF) and vascular endothelial cells (VECs) in the initiation of tumor angiogenesis. Except the existence of traveling waves, the understanding of (1) with m = 1 is very poor due to the singularity of logarithm ln w (at w = 0), where in particular the stability of traveling waves remains an outstanding open question to date except some instability results [11,39]. However for the linear consumption case m = 1, the model can be understood to some extend since the logarithmic singularity can be resolved by a Cole-Hopf type transformation ( [26,35] which converts the KS model(1) into a non-singular system of conservation laws as follows Though the singularity no longer exists in (3), a quadratic nonlinear convection is generated. In multi-dimensions, v is a gradient vector and the curl of v is intrinsic required to be zero, namely curlv = ∇ × v = 0. A characteristic feature of the transformed system (3) distinct from other system of conservation laws (e.g. see [4,9,47]) is that the parameter ε plays a dual role: coefficient of viscosity (diffusion) and nonlinear convection. Hence it is hard to justify the parameter ε > 0 is "good" or "bad" for analysis, and how to find a balance between the nonlinear convection and viscosity with the curl-free condition becomes an art of analysis. Indeed the transformed chemotaxis model (3) has been well understood in one-dimension for both ε = 0 and ε > 0 from various aspects such as the traveling wave solutions (cf. [5,21,30,32,33,34,35]), global dynamics of large-data solutions in R (cf. [28,37]) or in the bounded interval subject to various boundary conditions (cf. [29,50,52]). However it still remains poorly understood in multi-dimensions except few results on the small-data solutions (cf. [12,15,43,51]) or radial solutions (cf. [53]). In addition to these works, there was another class of results by considering singular limits of solutions to (3) as ε → 0. Such a topic is of particular interest since the vanishing as ε → 0 occurs concurrently to both viscosity and quadratic nonlinear convection in the transformed system (3). It is also of relevance since the chemical diffusion rate ε > 0 was assumed to be zero in the analysis of many early works (cf. [23,24,27]) on the grounds of simplicity and hence it is desirable to reveal the role of ε. Next we shall first recall existing results connecting the limit problem of ε → 0 and then propose our new questions. If the spatial domain is unbounded (i.e. x ∈ R N , N ≥ 1), it has been shown in [43,49,51] that both traveling wave solutions (see [49]) and global solutions of the Cauchy problem (see [43,51]) are uniformly convergent in ε, namely the solutions with ε > 0 converges to those with ε = 0 as ε → 0 in L ∞ -norm. If the domain is an interval say (0, 1), and zero mixed Neumann-Dirichlet (ND) boundary conditions are prescribed: u x | x=0,1 = 0, v| x=0,1 = 0, ε ≥ 0, it was shown in [52] that the solution is still uniformly convergent in ε. However if the Dirichlet boundary conditions are imposed, one cannot impose the boundary conditions for v with ε = 0 since otherwise the problem may be over-determined. In this circumstance, boundary layers may arise due to the possible mismatch of boundary conditions. This was first observed and numerically verified in a recent work by Li and Zhao in [29], and later was justified in [18]. Considering that the boundary conditions are dynamic in vivo environment for tumor angiogenesis, in this paper we consider the system (3) with time-dependent Dirichlet boundary values, and for simplicity hereafter we assume χ = D = 1 since their specific values are not important for our analysis. Hence precisely we shall consider the initial-boundary value problem (3) for (x, t) ∈ [0, 1] × [0, ∞) as follows: where α(t) and β(t) are boundary data depending on t. In (4) we always assume ε > 0. The non-diffusive initial-boundary value problem associated with (4) is Since now the boundary conditions are time-dependent, the global existence and asymptotic behavior of solutions may become elusive due to time-variable boundary data. Whether the boundary layer profiles for constant Dirichlet boundary data can be destructed by time-varying boundary data is also concerned. Hence we set two goals to this paper. First we show that the global strong solutions of the initial-boundary value problem (4) and (5) exists and boundary layer profile will arise as ε → 0 under mild conditions on boundary data α(t) and β(t), where the solution component u converges in L ∞ , v converges in L 2 while diverges in L ∞ . Second, we prove under certain constraints, the time-dependent boundary data α(t) and β(t) will act as the asymptotic profiles of solutions to (4) approaching some constant states. We remark that the approaches and estimates developed in previous works [18,29] for constant boundary conditions are not adequate for our current problem with time-dependent boundary data and various delicate boundary estimates and uniform-in-ε estimates are desired. In this paper we shall introduce the so called "effective viscous flux" technique employed in the study of the Navier-Stokes equations (see [17,36]) to gain the desired estimates to achieve our first goal. For the second goal, we develop a new entropy-like energy framework and fully explore the convexity of the entropy expansion to establish a basic L 1 energy estimate, on which the results of the asymptotic behavior of solutions are built up. We shall state our main results in the next section.
2. Statement of main results. To proceed, we first specify some notations for clarity. In the sequel, H k [0, 1] denotes the usual k-th order Sobolev space on [0, 1] , where we simply denote · := · L 2 [0,1] . We also use · L ∞ to denote · L ∞ [0,1] . Unless otherwise specified, we use C to denote a generic positive constant and C(t) denotes a generic positive constant which depends on t. The values of the constants may vary line by line according to the context. The first result of this paper on the existence and uniform-in-ε boundedness of global solutions to (4) is stated as follows.
Theorem 2.1. Assume that the initial and boundary data satisfy where c 0 is a positive constant. Then for any ε ≥ 0, the initial boundary value problem (4) has a unique global solution (u, v), such that for any T > 0, there hold that (u, v) ∈ L ∞ (0, T ; H 2 (0, 1)) ∩ L 2 (0, T ; H 2 (0, 1)), u ≥ 0 and where C(T ) is a positive constant dependent on T but independent of ε.
Next we shall state the result on the asymptotic behavior of solutions to (4).

Remark 2.
The conditions on the time-dependent boundary data admit a family of functions approaching constant states with certain decaying/growth rates, such as algebraic or exponential, as time goes to infinity. Since the temporal integrability of the boundary data is not required, boundary functions which approach constant states with slow decaying/growth rates, such as α(t) = 2 ± 1 (1+t) for 0 < 1, are permitted. The long-time behavior result indicates that the solution decays asymptotically and its decay profile is determined by the boundary data. In addition, we require α(t) to be bounded from below and above away from zero, which is consistent with and generalizes the previous result [29] wherein the boundary data are constants.
Finally we reverse the results of the transformed system to the pre-transformed chemotaxis model (1) with m = 1. The counterpart of the initial-boundary value problem of (1) with m = 1 corresponding to (3) reads as Then we have following results for (10).
Remark 4. The result in Theorem 2.5 (ii) means that when α > εβ 2 (this condition is satisfied naturally in the case of β(t) = 0), the L ∞ -norm of w will exponentially decay to zero as time goes to infinity. However the result for the case α < εβ 2 is unclear. But our result implies that if the solutions diverge in this case, the divergence rate is not faster than an exponential rate.
The rest of this paper is organized as follows. In section 3, we shall establish the global existence of solutions of (4) and prove Theorem 2.1. In section 4, we explore the vanishing limits as ε → 0 of solutions (boundary layer solutions) and prove Theorem 2.3. The results on the asymptotic behavior of solutions (Theorem 2.4) will be shown in section 5, and the proof of Theorem 2.5 will be given in section 6. Finally we show the numerical simulations to illustrate boundary layer profiles and interpret our analytical results in section 7.
3. Proof of Theorem 2.1. In this section, we will prove Theorem 2.1. First, using the standard arguments (e.g. see [50]), one can show the local existence of solutions to (4).
Next we derive some a priori uniform-in-ε estimates of solutions, which not only extend the local solutions to global ones, but also play important parts in investigating the vanishing diffusion limit. We depart from the following boundary estimates on (u x , v x ). and Proof. By integrating (3) over (0, x) and using the boundary condition in (4), we have Then, integrating (17) with respect to x over (0, 1), yields (15). Similarly, (16) is obtained.
Let the assumptions in Theorem 2.1 hold. Then for any t > 0, there exists a positive constant C(t) which is dependent on t but independent of ε, such that Proof. To resolve the logarithmic singularity in the following estimates, inspired by [29], we make a technical treatment by introducing a change of variableũ = u + 1. Thus, problem (4) turns into Multiplying the first equation of (19) by lnũ and integrating the result by parts over [0, 1], we have d dt where η =ũ lnũ −ũ + 1 + R and R is a positive constant to be determined later.
Multiplying the second equation of (19) by v and integrating the result by parts over [0, 1], we have Adding (20) to (21) and integrating the result over (0, t) yield that Using the factũ ≥ 1 and Cauchy-Schwarz inequality, we get On the other hand, from the boundary conditions in (6), for any t > 0 there is a constant c 1 (t) which may depend on t such that Thus, using Lemma 3.2, (23), integration by parts and Cauchy-Schwarz inequality, we can estimate the third term on the right-hand side of (22) as follows: Noting thatũ x = u x , for the fourth term on the right-hand side of (22), we use Lemma 3.2, (6) and the integration by parts to get where we have used the fact that α(t) ≥ 0, (23) and ln(α(t) Inserting the above estimates into (22) yields which results in (18) by the Gronwall's inequality.
Lemma 3.4. Let the assumptions in Theorem 2.1 hold. Then for any t > 0, there exists a constant C(t) > 0 which is dependent on t but independent of ε, such that Proof. Multiplying the first equation of (4) by u, integrating the result by parts over [0, 1], and adding the resultant equality to (21), we have Integrating (26) with respect to t and using the boundary conditions in (4), we have For the second and third terms on the right-hand side of (27), by the Gagliardo-Nirenberg and Cauchy-Schwarz inequalities, we have where in the second inequality we have used (18). The last term on the right-hand side of (27) has been well estimated in (24). Therefore, we get from (18) and (24) that For the fourth term on the right-hand side of (27), by Lemma 3.2 and integration by parts, we have where we have used (18) and (23). Substituting these estimates into (27), we obtain which, together with Gronwall's inequality, yields (25).
Lemma 3.5. Let the assumptions in Theorem 2.1 hold. Then for any t > 0, it holds that where the constant C(t) is independent of ε but dependent on t.
Proof. We first multiply the first equation of (4) by u t and integrate the resulting equation over [0, 1] × [0, t] to get where in the last equality we have used the boundary conditions in (4). For the second term on the right-hand side of (29), by the Gagliardo-Nirenberg and Cauchy-Schwarz inequalities, (18) and (25), we have From Lemma 3.2 and integration by parts, the last term on the right-hand side of (29) can be estimated as follows: where (18) and (23) have been used. Substituting (30) and (31) into (29), we have Next, in order to obtain the estimate of Multiplying the first equation of (33) by u t and the second by v t , adding the results and integrating it over [0, 1] × [0, t], we have For I 1 , integrating by parts and using the boundary conditions in (4), we obtain Using (25), Gagliardo-Nirenberg and Cauchy-Schwarz inequalities, we can estimate the second term on the right-hand side of (35) as which updates (35) as dτ.
By Cauchy-Schwarz inequality, we have Integration by parts implies In order to estimate the boundary terms in (35) and (36), we follow the same procedure as in Lemma 3.2 and get and Using (23) and (37), integration by parts and Cauchy-Schwarz inequality, we can estimate the last term on the right-hand side of (35) as follows:

HONGYUN PENG, ZHI-AN WANG, KUN ZHAO AND CHANGJIANG ZHU
Similar to (24), we can estimate the last term on the right-hand side of (36) as where we have used (38). Next, we need to estimate I 4 . Integrating by parts, using the boundary conditions in (4), Gagliardo-Nirenberg and Cauchy-Schwarz inequalities, we obtain Substituting the estimates of I i (i = 1, 2, 3, 4) into (34), we get Using Gronwall's inequality and (25), we obtain which together with (32) leads to u x 2 + t 0 u t 2 dτ ≤ C(t). This, along with (39), leads immediately to (28).
The next lemma gives the estimate of L ∞ -norm of (u x , v). It turns out it is not easy to gain them by the routine procedure like the iteration method. Motivated by the studies for the Navier-Stokes equations (cf. [17,19,36]), we here introduce the following so-called "effective viscous flux G(x, t)": From the first equation of (3), it is easy to see that The quantity effective viscous flux G will play an important role deriving the L ∞norm of (u x , v).
Lemma 3.6. Let the assumptions in Theorem 2.1 hold. Then for any t > 0, there exists a constant C(t) > 0 which is independent of ε, such that Proof. Multiplying the second equation of (4) by 2nv 2n−1 (n ≥ 1 is an integer), integrating the result by parts over (0, 1), we obtain d dt where we have used the boundary conditions in (4) and the non-negativity of u and v 2n . Now, we need to control G L ∞ . Using Gagliardo-Nirenberg inequality, (40)-(41), (25) and (28), we get Using Lemma 3.2 and integration by parts, we have where we have used (23) and (25). Then it follows from (43)- (45) and Gronwall's inequality that Then, raising the power 1 2n to both sides of (46) and letting n → ∞, we obtain that v L ∞ ≤ C(t).
The following refined estimates of (u, v) will play an important role in the study of vanishing diffusion limit.
Lemma 3.7. Let the assumptions in Theorem 2.1 hold. Then for any t > 0, it holds that where the constant C(t) is independent of ε but depends on t.
Proof. Multiplying the first equation of (4) by −2εu xx in L 2 , using Cauchy-Schwarz inequality, (23) and Lemmas 3.4-3.6, we have Next, we differentiate the second equation of (3) with respect to x, and subtract the resulting equation from the first equation of (3), to get Multiplying (50) by 2εv x and integrating the result over (0, 1) yield I 5 can be estimated by Cauchy-Schwarz inequality and Lemmas 3.4-3.6 as For I 6 , we use integration by parts, Cauchy-Schwarz inequality, Sobolev embedding theorem and Lemmas 3.4-3.6, to get Noting that εv xx = v t − u x + ε(v 2 ) x , using (23), (42), the Gagliardo-Nirenberg and Cauchy-Schwarz inequalities, we have where we have used the following inequality derived from Gagliardo-Nirenberg inequality and Young inequality Substituting the estimates of I i (i = 5, 6, 7) into (51) and adding the resulting inequality to (49) yield Then the Gronwall's inequality leads to which immediately gives (48) and completes the proof of Lemma 3.7.

Vanishing diffusion limit and boundary layer solutions.
This section is concerned with the vanishing diffusion limit and boundary layer solutions. We first give the global existence of solutions to the non-diffusion problem (5).
Lemma 4.1. Assume that the initial and boundary data satisfy Then for any 0 < T < ∞, there exists a unique strong solution (u, v) to (5) in Proof. Noting that the energy estimates established in Theorem 2.1 still hold true for ε = 0, i.e., for any t > 0, there is a constant C(t) > 0, such that Next, we will give the estimate of v x . Differentiating the second equation of (5) with respect to x, then subtracting the resulting equation from the first equation Multiplying (53) by 2v x , integrating by parts over (0, 1) and using Cauchy-Schwarz inequality and (52), we deduce Applying Gronwall's inequality, we have This together with the first equation of (5) and (52) means Next differentiating (53) with respect to x, we have Multiplying (56) by 2v xx , integrating by parts over (0, 1), using Cauchy-Schwarz inequality, (52) and (54), we deduce Applying Gronwall's inequality and (52), we have v xx 2 ≤ C(t).

Proof of Theorem 2.3 (i).
Let (u ε , v ε ) and (u 0 , v 0 ) be the solutions to the initial boundary value problems (4) and (5), respectively. Let us set Then, by a straightforward calculation, we find that (ϕ ε , θ ε ) satisfies the following the initial boundary value problem: with initial data (ϕ ε , θ ε ) (x, 0) = (0, 0), and boundary condition: Assume that the assumptions listed in Theorem 2.1 and Lemma 4.1 are satisfied. Then for any t > 0, there exists a positive constant C(t) which is independent of ε, such that (61) Proof. Multiplying the first and second equations of (57) by 2ϕ ε and 2θ ε respectively, integrating the result by parts on [0, 1], using the boundary condition (59), we have d dt By Cauchy-Schwarz inequality and Theorem 2.1, we have
Finally, Theorem 2.3 is a consequence of Lemma 4.2.

Proof of Theorem 2.3 (ii).
Inspired by a recent work [19], we first establish the following lemma by the weighted L 2 -method dedicating to the boundary layer solutions.

5.
Large-time behavior. In this section, we prove Theorem 2.4. For the reader's convenience, we restate the initial-boundary value problem, which reads as The proof of Theorem 2.4 is divided into four steps contained in a series of subsections. First of all, we note that, due to the conditions of Theorem 2.4 and maximum principle, it holds that u(x, t) ≥ 0, provided that the solution exists. We depart with a basic estimate involving the logarithmic expansion of u.

Entropy estimates.
Lemma 5.1. Let the assumptions in Theorem 2.4 hold. Then there exists a constant C > 0 which is independent on t and ε, such that denotes the entropy expansion.
Proof. We divide the proof into three steps.
Step 1. By a direct calculation, we can show that By using the first equation of (70) and noting α depends only on t, we deduce that Then plugging (72) into (71), we find After integrating (73) over [0, 1], and using the boundary conditions we have Since β is independent of x, we derive from the second equation of (70) that Taking the L 2 inner product of (75) with v − β, we have Note that So we update (76) as By adding (77) to (74), we get that where We remark that in [29] the two terms on the right hand side of (78) vanish, due to the constant boundary conditions. The treatment of these non-constant terms is one of the major differences between this paper and [29].
Step 2. In this step, we derive an energy bound for the L 1 norm of u in terms of the entropy expansion defined by (79). We remark that under the Dirichlet type boundary conditions, the L 1 norm of u is not a conserved quantity. Hence, the energy method established in [31] for the mixed Neumann-Dirichlet boundary value problem can not be utilized for the Dirichlet boundary conditions. Luckily, such an issue was previously resolved in [29] for constant Dirichlet boundary data by developing a new approach through higher order nonlinear cancellation. Though such a technique also works for the time-dependent Dirichlet boundary conditions and can produce a uniform-in-time energy estimate for the low frequency part of the solution, the proof is lengthy and one needs more constraints on the boundary data to close the energy estimate. In this paper, we develop a very new approach (which has never appeared in any related work) to settle down the energy estimate for the low frequency part of the solution. The idea is to fully explore the convexity of the entropy expansion E(u, α) and compare it with a linear function. For this purpose, we set Then it can be readily checked that which imply that F α (u) ≥ 0 for any u ≥ 0. This leads to Step 3. By plugging (80) into (78), we see that where we used the first assumption of Theorem 2.4 and the Cauchy-Schwarz inequality. By applying the Gronwall's inequality to (81), we have By using the second assumption of Theorem 2.4, we deduce from (82) that where the constant C is independent of time and ε. By plugging (83) into (81), then integrating the resulting inequality with respect to time, we have in particular, where the constant C is independent of time and ε. This together with (83) completes the entropy estimate and hence the proof of Lemma 5.1.

L 2 estimates.
To perform further energy estimates, we let where (u, v) satisfies (4). Then (ũ,ṽ) satisfies Proof. Taking the L 2 inner product of the first equation of (85) withũ, we have Taking the L 2 inner product of the second equation of (85) withṽ yields Multiplying (87) by α, we have where we have applied integration by parts to the first term on the right hand side of (88). Adding (88) to (86), we have Now, we estimate the first term on the right hand side of (89) by using the L 1 estimate obtained from the previous subsection. To this end, we observe that where ũ 2 L ∞ can be estimated through the following procedure: Step 1. Note that for any x ∈ [0, 1] and t > 0, u(x, t) = x 0ũ y dy, which implies

HONGYUN PENG, ZHI-AN WANG, KUN ZHAO AND CHANGJIANG ZHU
Step 2. Sinceũ = u − α and α is independent of x, it holds thatũ x = u x . Then by Hölder's inequality and the positivity of u, we have Step 3. By applying (83) to (80) and using the first assumption of Theorem 2.4, we obtain Step 4. By applying (92) to the first term on the right hand side of (91), we obtain By plugging the preceding estimate into (90), we find which updates (89) as Note that So we update (94) as where we have used the first assumption of Theorem 2.4. Applying the Gronwall's inequality to (95) and using (84) and the second assumption of Theorem 2.4, we find that for some constant C which is independent of t and ε. Plugging (96) back into (95), then integrating the resulting inequality with respect to time, we conclude that where the constant C is independent of t and ε. This completes the energy estimate for the low frequency part of the solution.
Next, we shall move on to the estimation of the first order derivatives of the solution.
where the constant C is independent of t, but is inversely proportional to ε.
Proof. Taking the L 2 inner products of the first equation of (85) with −ũ xx , and the second with −ṽ xx , respectively, then adding the results, we have For the right hand side of (98), we first apply the basic Cauchy-Schwarz inequality to deduce For the L ∞ norms appearing in the above estimates, we note that since both the functionsũ andṽ equal zero on the boundary, it holds that and the same is true forṽ. Hence, we can update I 1 and I 3 as I 1 ≤ 1 4 ũ xx 2 + 8 ṽ x 2 ũ x 2 + 4α 2 ṽ x 2 + 4β 2 ũ x 2 ; 7. Simulations and implications. In this section, we numerically solve system (4) to illustrate the boundary layer profile (u ε , v ε ), verify our analytical results and discuss boundary effects. The model is solved in the interval [0, 1] with MATLAB based on the finite difference scheme with mesh size ∆x = 0.001, ∆t = 0.01. We first choose the initial and boundary data satisfying the requirements in (6) and implement numerical computations to the system (4). The solution profile (u, v)(x, t) at time t = 0.2 as ε → 0 is plotted in Fig.1. For the sake of comparison, we also numerically solve the non-diffusive problem (5) in the absence of boundary conditions for v and plot the numerical solution at t = 0.2 in Fig.1. We find from the simulations that the solution profile u(x, t) is convergent with respect to ε in [0, 1], whereas the solution profile v(x, t) becomes increasingly sharp near the boundary as ε → 0 and boundary layers arise. Outside the boundary layer (i.e. in the interior of [0, 1]) the solution profiles v(x, t) for small ε > 0 and ε = 0 match well.
In Fig.2, we proceed to plot the time evolution of the same solution solved in Fig.1 to observe the asymptotic profiles, where we find that large-time profiles of the solution is elusive. This is because the boundary data chosen in Fig.1 vary (oscillate) in time. But the simulations show that the boundary layer profiles (sharp transition near boundaries) persist in time given small ε > 0. However if we impose some decay properties to the boundary data, the results of Theorem 2.4 show that the asymptotic behavior of solutions may become tractable and converge to some constant states, where the decay profiles of solutions are determined by the boundary data. Here we numerically explore this analytical finding. For this, we choose the initial and boundary data (see the caption of Fig. 3) such that the decay of boundary data for u is exponential and for v is algebraic, as well as the initial data satisfying the compatibility conditions at the end points x = 0, 1, as required by Theorem 2.4. We plot the numerical solution profiles in Fig.3 at different times showing that the solution (u, v) will approach constant states as time evolves. In particular, we find that the convergence of u is much faster than that of v. This complies with our analytical results in Theorem 2.4 that the decay rates of u and v are same as the boundary data α(t) and β(t), respectively, where the former (exponential decay) is much fast than the latter (algebraic decay).  Finally we shall discuss some biological insights gained from our analytical and numerical results. In view of model (1) with m = 1 and the transformation (2), we see that the quantity v represents the velocity of chemotactic flux crossing the boundary. Therefore the results in Theorem 2.5 imply that if the chemical diffusion is small, although both cell density and chemical concentration have no boundary layers, the chemotactic flux (i.e. the term u(ln w) x = uv) may change drastically near the boundary since v has boundary layers. If the boundary data have oscillating properties, this phenomenon will persist in time. However if the boundary data have some decay properties, the boundary layer may vanish as time evolves. Therefore the nature of boundary date play an essential role in determining the solution behaviour near the boundary and large-time dynamics.