Dynamics in a Chemotaxis Model with Periodic Source

: We consider a system of two differential equations modeling chemotaxis. The system consists of a parabolic equation describing the behavior of a biological species “ u ” coupled to an ODE patterning the concentration of a chemical substance “ v ”. The growth of the biological species is limited by a logistic-like term where the carrying capacity presents a time-periodic asymptotic behavior. The production of the chemical species is described in terms of a regular function h , which increases as “ u ” increases. Under suitable assumptions we prove that the solution is globally bounded in time by using an Alikakos-Moser iteration, and it fulﬁlls a certain periodic asymptotic behavior. Besides, numerical simulations are performed to illustrate the behavior of the solutions of the system showing that the model considered here can provide very interesting and complex dynamics.


Introduction
Chemotaxis is the ability of some living organisms to direct their movement in response to the presence of a chemical gradient. This response can be either positive (chemoattractant) or negative (chemorepellent). Mathematical models for chemotaxis have been studied since 1970 when Keller and Segel proposed a system of two parabolic equations involving nonlinear second order terms in the form −∇· (χu∇v) in the u-equation. Since the publication of the model, an extensive mathematical literature has treated the topic, see also [1]. To present an exhaustive literature review is not the aim of this article, therefore we refer to the reader to the survey works of Horstmann [2,3], and Bellomo et al. [4] for more details (see also [5]).
The mathematical model that we study in this article describes the behavior of a biological species "u" in terms of a PDE of parabolic type. The problem is posed in a bounded domain Ω ⊂ R n , with a regular boundary ∂Ω as follows x ∈ Ω, t > 0, The equation includes the linear diffusion of "u" which also moves following the direction of the chemical gradient of a non-diffusive substance "v". The chemotactic coefficient χ is assumed to be constant and positive, i.e., the biological species moves to a higher concentration of "v". The logistic term includes a carrying capacity that limits the growth of "u" and it presents a spatial and time dependence, f (x, t). Then, the reaction part The global existence of solutions for the fully parabolic system is achieved by employing an iterative method based on the Alikakos-Moser iteration. By using an energy method through a Lyapunov functional, the convergence of the solution to a homogeneous in space and periodic in time function u * is given. The parabolic-elliptic case, i.e., for v satisfying the equation has been studied in Negreanu, Tello and Vargas [7], where the global existence and similar asymptotic behavior are done. In that case, the proof follows a sub-super solutions method already featured in Pao [8], Tello and Winkler [9], Galakhov, Salieva and Tello [10] and Negreanu and Tello [11,12] among others. In [7], the problem is addressed for a non constant function f and u * satisfying the ODE where f * is a periodic in time function such that sup x∈Ω | f (x, t) − f * (t)| → 0, as t → ∞.
In Issa and Shen [13] the logistic term is u a 1 (x, t) − a 2 (x, t)u − a 3 (x, t) Ω udx and the authors got the existence of periodic solutions when the coefficients a i (for i = 1, 2, 3) are periodic in time.
Parabolic-ODE systems with chemotactic terms have been considered from the last three decades, and after the pionering works of Levine and Sleeman [14] and Anderson and Chaplain [15] modeling tumor angiogenesis, a considerable number of authors have analyzed such models. In Othmer and Stevens [16] and Stevens [17], the authors address a Parabolic-ODE system of chemotaxis passing to the limit from a discrete to a continuous system of equations. Concerning angiogenesis, the model has been raised in Kubo and Suzuki [18], Suzuki [19] and Kubo, H. Hoshino and K. Kimura [20]. Mathematical analysis of these models with two equations can be found in Fontelos, Friedman and Hu [21], Friedman and Tello [22] and Negreanu and Tello [11] among others. Systems with three or more equations involving chemotaxis and diffusive or non-diffusive processes also appear in ecology and other biological applications (see [12,23]). In Ref. [24] the authors study a similar Parabolic-Parabolic-ODE system where the chemosensitivities χ 1 , χ 2 are non-constant. Global existence and convergence of the solution to a steady-state (1, 1,w) satisfying h(1, 1,w) = 0 are presented under suitable assumptions on the coefficients and the spatial dimension of the domain. The results in [24] have been improved in Mizukami and Yokota [25] for a larger range of parameters.
Also in the context of cancer dynamics, chemotactic systems with non-diffusive equations have been recently used to model cancer cell invasion in Stinner, Surulescu and Winkler [26] with a model consisting of six equations where the cancer cells behavior is described by a parabolic equation with chemotactic terms. Denoting by "u" the cancer cells density, by "v" the fibers of the extracellular matrix (ECM) and by "l", "y 1 " and "y 2 " the concentration of chemoattractant, integrins bound to ECM fibers and integrins bound to proteolitic residuals then, their model is the following The authors prove the existence of global weak solutions together with some boundedness properties. The proof is based on the properties of the functional Notice that, in this case, three of the variables leading the movement satisfy ordinary differential equations (see also Stinner, Surulescu and Uatay [27], Tao and Winkler [28], Zhigun, Surulescu and Hunt [29] and Zhigun, Surulescu and Uatay [30] for similar models). Throughout the article we use the notation Ω t = Ω × (0, t), for t ∈ (0, ∞], we assume, without loss of generality, that |Ω| = 1 and we denote by g the function We work under the following hypotheses 1. The positive initial data (u 0 , v 0 ) of (1) satisfy, (u 0 , v 0 ) ≡ (0, 0) and for some s > max{4, n} and 2.
There exists a periodic function f * verifying 3.
Moreover, there exists a positive constant c such that with some ε as in (7).

Main Results
Our particular analysis will address the initial-boundary value problem (1) in a bounded open domain Ω ⊂ R n , where the initial data are as in (4). The issue of the global solvability is presented in the following theorem.
Afterwards, we study the asymptotic properties of the solutions. We introduce the function u * as set out by for u * 0 defined by and f * as in (6). Notice that u * satisfies Equation (2) and it is an homogeneous in space and periodic in time function. We denote by v * (t), the solution of the ordinary differential The following assertion is the main result on the asymptotic behavior of solutions of (1). Theorem 2. Assume (4)-(11) and let us denote by (u, v) the corresponding solution of (1) from Theorem 1. Then, (u, v) has the following asymptotic behavior for any p ∈ [1, ∞), where u * and v * are given by (12) and (13), respectively.
The paper is organized as follows: in Section 3 we proof the existence of a unique pair of classical solutions. A first key step consists in findining a maximal weak solution following [31], and we then get the boundedness of the solution. As a crucial ingredient in our derivation of a L ∞ (Ω) bound for u, we employ a Alikakos-Moser-type iterative procedure [32]. By means of these and some further higher regularity properties will assert the statements on global existence and boundedness of u and v from Theorem 1. Our collection of estimates of Section 3 will moreover turn out to be sufficient to derive the stabilization result from Theorem 2 in Section 4 through an analysis into two steps. First, we prove that the solutions converge to their respective averages, i.e., Ω udx, Ω vdx, using energy estimates to conclude that these averages converge to the functions u * and v * , respectively. Finally, in Section 5 we perform a brief numerical study of the system under consideration. Some of the results presented in this paper were announced in [33].

Global Existence of Solutions
The present section is devoted to the proof of Theorem 1. We study the local-intime existence of classical solutions to (1) and we prove some preliminary technical facts. In order to prove the global existence of the solutions, we first obtain the local existence using classical results on partial differential equations and then we conclude the proof by constructing uniform bounds.

Lemma 1.
Let Ω ⊂ R n be an open bounded set with smooth boundary, assume that the initial data, (u 0 , v 0 ), are as in (4) and hypotheses of Theorem 1 hold. Then, there exists a maximal existence time T max > 0 such that system (1) has an unique non-negative classical solution Proof. We consider the system (6.2) of [31] where and We can rewrite then (1) as follows with the same initial data as (1). We apply Theorem 6.4 in [31] and consider maximal interval of existence. So, the local-in-time existence for (1) is proved. In order to see the non-negativity of u we introduce the following change of variables: Then we can rewrite the first equation in (1) as Now, deriving with respect to the spatial variable in the previous equation we get Then, the first equation of (1) becomes we multiply by e −χv to get Notice that the equation for v remains as an ordinary differential equation So, the original system (1) becomes (16) and (17) together with the initial datã and the Neumann boundary condition Finally, the Maximum Principle for parabolic equations and the regularity of h prove the non-negativity of u, taking into account that Hypotheses (9) and (10) on h and the Maximum Principle applied to (17) prove This completes the proof.
Let us now collect some basic properties thereof which in our subsequent analysis will play important roles not only by providing some useful fundamental regularity features, but also by establishing the first quantitative information (18) on large time behavior. We remember that Ω ∞ = Ω × (0, ∞) for the next matches, and also, · L ∞ (Ω ∞ ) = sup t>0 · L ∞ (Ω) .

Lemma 2.
Under Hypotheses (4)-(11), the total mass of the component u(x, t) of the solution to (1) is bounded: Proof. We integrate the first equation of (1) directly over Ω to get Applying the Cauchy-Schwarz inequality, since |Ω| = 1, we directly obtain Finally, (18) is a consequence of the Maximum Principle applied to (20), i.e.,

Lemma 3.
Under the same assumptions of Lemma 2, the solution to (1) satisfies By the previous lemma it follows Finally, since t 0 ≤ 1 we have thereby completes the proof.

Lemma 4.
Under hypotheses of Theorem 1, the following assertion is verified: there exists a positive constant c 3 defined by We are now prepared to perform an iterative argument of Alikakos-Moser type in order to derive L ∞ (Ω) bounds for u and v.
The proof starts with the following lemma.

Lemma 5.
Let g(v) be defined by (15), then, for p ≥ 2 the following estimate holds where c is the constant given in assumption (10).

Proof.
We proceed by induction in p, then, for p ≥ 2 we have For the first integral in (23) we infer that From the expression of the above identity, we deduce We look now at the last term of (23). By the Mean Value Theorem and assumption (10) we have Moreover, for the restant term of (23), we have We now replace (24)- (26) in (23) to get which yields (22) and the proof ends.

Lemma 6.
Let us consider p ≥ 1 and g(v) as in (15). Let be a positive constant defined by then, there exist c 4 > 0 and c 5 > 0 given by and Proof. For p = 1 the result is a consequence of Lemmas 2 and 3. For p ≥ 2 we proceed by induction and assume the result for p − 1, i.e., Taking p ≥ 2, thanks to (22), we have We first recall the Young's inequality: multiplying it by g 1−p we get 1 u p g 1−p ≤ u p+1 g 1−p + 1 (p + 1) p+1 g 1−p , which is equivalent to −u p+1 g 1−p ≤ − 1 u p g 1−p + 1 (p + 1) p+1 g 1−p , we integrate in space over Ω, and in view of g 1−p ≤ 1, we get Thanks to the definition (27) of , we have We replace (32) into (31) to get By solving the differential Equation (33) after integration in time, we obtain Dropping the nonpositive term and making use of a favorable cancellation, it yields and t 0 e s−t By definition of c 4 it follows and due to (35), the following inequality holds

Lemma 7.
Under the assumptions of Theorem 1, we have where c 5 has been defined in Lemma 6.
Proof. According to Lemma 6 we have that Since Ω u p g −p dx we take limits when p → ∞, to obtain (36).
Proof. By contradiction, we assume that for any v > v 0 L ∞ (Ω) there exists t 0 > 0 such that v(x, t 0 ) = v which is the first t 0 fulfilling this condition. Since by assumption (4) v 0 < v, v must be an increasing function in a neighborhood of t 0 . Then, by applying (11), we obtain Thanks to assumption (11) we have that for v large enough v t (t 0 ) < 0, which is a contradiction and the proof ends.
The above results entail the claimed qualitative properties of u: Proof. The result is a consequence of Lemmas 7 and 8.

Proof of Theorem 1.
The global existence of (u, v) over Ω × (0, ∞) is a direct consequence of the local existence (Lemma 3.1, Theorem 6.4 in [31]) and the uniform boundedness of (u, v) in L ∞ (Ω) established in the previous Lemmas.

Asymptotic Behavior
The main propose of this section is to demonstrate Theorem 2, i.e., to obtain the convergence of the solution (u, v) to (u * , v * ). We proceed in two steps: first of all we get the convergence of the solution u to its average Ω u, to get later the convergence of the average to the periodic function u * given by (12). For it, we need to prove the boundedness of |∇v| in L 2 (Ω). The result is enclosed in the following lemma.

Lemma 10.
Suppose that the assumptions of Theorem 2 hold. Then, there exists c 6 > 0, independent of t, such that where v is the solution of (1). (22), for p = 2, and integrate over (0, t) to obtain, after routinary computations and thanks to Lemmas 4, 6 and 7 t 0 e (s−t)

Proof. We consider Equation
for any > 0. Recalling that v satisfies then taking gradients we get Now, we multiply (38) by ∇v and integrate over Ω to obtain, in view of assumptions (9) d dt and therefore, by the Young's inequality d dt After integration in time we get and due to (37), we conclude the lemma.
Similar results can be found in Tao and Winkler [28] (Theorem 1.1) and [34] for parabolic-elliptic and fully parabolic systems.

Proof.
We divide by u * in (2) and integrate over (0, t) to obtain the result.
We now define the positive function thus we achieve the following.

Lemma 13.
Under the assumptions of Theorem 2, there exists a positive constant c 7 independent of t such that the following estimate holds Proof. We integrate the first equation of (1) over Ω and in view of Since f and f * are uniformly bounded, we have for any δ > 0. We take δ = 1/4 and then We divide by Ω udx to get Now, we consider the following functions Functionals of quite a similar form have previously been used in several works on related chemotaxis problems, e.g., in [35]. Notice that F 1 ≥ 0 and F 2 ≥ c 0 . Let c 1 be defined in (18), then and also d dt We take gradients in the equation of v, multiply the obtained equation by λ∇v and integrate over Ω to get

Now we add both expressions to obtain
We apply the Cauchy-Swartz inequality to the term ∇v∇u[χ + λuh u ]/u then, operating we achieve which is reduced to Due to the discriminant of the polynomial is given by 16h v (h v + χuh u ), which is positive, we have two different roots λ + and λ − that are both positive. Since Then, we take λ = λ − to obtain Through the inequality (43) it results After integration over (0, T) and taking limits when T → ∞ we conclude the lemma.

Lemma 14.
Under the assumptions of Theorem 2 the following estimate holds with c 8 a positive constant.
Proof. We first notice that p(λ) defined in (45) achieves its minimum at ] is a compact set of R 2 . Due to (44) we get We now proceed as in Lemma 13 and we obtain After integration over (0, ∞), in view of we end the proof.
We have the following boundedness Proof. After integration in the time variable the expression with c 10 > 0. In view of the boundedness of u we have and the proof ends. In Negreanu, Tello and Vargas [6], a similar problem is studied for the fully parabolic system.
Proof. The following relations hold: By applying the Young's inequality we have The boundedness of u and Lemma 14 imply the result.
The following lemma is used to prove the behavior of the solution. The proof follows Lemma 5.1 in Friedman-Tello [22], where k is uniformly bounded, i.e., |k | < c. Here, the boundedness of k is replaced by a weaker assumption given in (iii). Lemma 17. Let k : R + → R a function satisfying (i) k(t) ≥ 0 for any t ≥ 0, (iii) k ≤ c for any t ≥ 0, then, k(t) → 0 as t → ∞.
Proof. By contradiction, we assume that there exists a sequence t n such that t n → ∞ and lim sup n→∞ k(t n ) ≥ > 0. Then, there exist a subsequence t n such that t n ≥ t n−1 + 1 and Then, k ≥ 4 in the interval [t n − a, t n ] for a := min{1, c /4}. So Proof. We consider k 1 defined in (39), then, thanks to Lemmas 13 and 16 we have Now, we define k 2 as follows By recalling the definition of F 2 (t) as in (42) due to (39), we get We multiply by F 2 and due to the Mean Value Theorem we claim Notice that Lemma 2 implies k 2 ≤ c 13 F 2 2 for some positive constant c 13 . Therefore, there exists c 14 > 0 such that In view of Lemma 2, assumption (6) and Lemma 9 it is easy to see that Now, by Lemma 17,(48) and (49) we obtain Since by taking into account (47) and (50), we get and the proof ends.
In order to obtain v − v * L 2 (Ω) → 0 as t → ∞, we proceed as before in the following lemma.

Lemma 19. Under assumptions (4)-(8), the solution v fulfills
We call z = v − v * , by multiplying by z the above equation and after integrating over Ω, it yields where we have applied the Hölder inequality to the last term. Now, by assumption (9) it results 1 2 where θ(t) is uniformly bounded. We obtain the result by solving the differential inequality.
Proof. (End of the proof of Theorem 2.) The asymptotic behavior (14) of (u, v) is a direct consequence of Lemmas 18 and 19 and the uniform bounds of u and v established therein.

Numerical Tests
Now we show some numerical results for the purpose of further clarifying that all conditions in the statement of the previous theorems play a relevant role in the behavior of the solution of (1). The suppression of some of the above conditions, together with the election of the initial data, may end up in the existence of blow-up of the solutions. We illustrate numerically the uniform boundedness and the convergence for the solution (u, v) to obtain a numerical validation of Theorem 2. We use the Generalized Finite Differences Method for the space discretization and we performed several tests showing the explosion of solutions in the case that certain hypotheses of Theorem 2 are not verified as thus the asymptotic behavior of solutions.

Example 1: Uniform and Periodic Asymptotic Behavior
For our purpose, let us consider µ = 1 and χ = 0.3. The initial data used are and the function f is which fulfils assumptions (6) and (7). Direct calculations lead The function h(u, v) = ue −χv − v, as it is easy to check, is in accordance with the hypothesis of the theorem. We find the solution of the second ODE of (2) and (13) by numerical integration using the ode45 function of Matlab R2019a. In Figure 1, the u, v−solution is presented for 0, 0.5, 1 and 20 s. Table 1 shows the l ∞ norm of the discrete (u, v)−solution and the value of (u * (t), v * (t)) for different times. In Figure 2 we illustrate the asymptotic solution u * , v * (solid line) and the value of the discrete solution for times in [0, 20].   Our simulations are in keeping with the theoretical results about global existence and boundedness of solutions to (1).

Example 2: Blow-Up Solutions
Next in order we present the results for the explosion of solutions in the case that certain hypotheses of Theorem 2 are not verified. We choose µ = 1 and χ = π. The initial data used are u 0 (x, y) = 8, v 0 (x, y) = 0.7e −10[(x−1.2) 2 +(y−0.8) 2 ] and the same function f of the previous example. Now, we consider h(u, v) = u − v, clearly it does not fulfil the assumptions. As we see from Table 2 and Figure 3, the solutions become unbounded before 0.40 s.  The formation of various patterns due to the effect of chemotaxis rate, domain size, initial data, the nature of the functions f and h or the complexity arisen in the solutions for large values of the chemotactic term are our goal for immediate studies.

Conclusions
We have obtained, under suitable assumptions, that the solution of a chemotaxis system is globally bounded in time by using an Alikakos-Moser iteration, and it fulfills a periodic asymptotic behavior. A possible future work is the consideration of the nonconstant chemosensitivity, χ(u, v), as in [23][24][25][26][27][28][29][30][31][32][33][34][35][36], and to also consider biological systems with two species, two chemotactic terms and one chemical substance verifying a similar equation as in (1). Furthemore, a parabolic-parabolic-ordinary system with periodic terms serves as a model for some chemotaxis phenomena and appears naturally in the interaction of two biological species and a chemical. The presence of the periodic terms has a strong impact on the behavior of the solutions. We would find conditions on the system's data that guarantee the global existence of solutions, the convergence to some periodical solutions of an associated ODE's system. We got a similar result in [37] for a parabolic-parabolic-elliptic system.

Funding:
The authors acknowledge the support of the Escuela Técnica Superior de Ingenieros Industriales (UNED) of Spain, project 2021-IFC02. This work is also partially support by the Project MTM2017-83391-P DGICT, Spain.