Global well-posedness and pattern formations of the immune system induced by chemotaxis

Abstract: This paper studies a reaction-diffusion-advection system describing a directed movement of immune cells toward chemokines during the immune process. We investigate the global solvability of the model based on the bootstrap argument for minimal chemotaxis models. We also examine the stability of nonconstant steady states and the existence of periodic orbits from theoretical aspects of bifurcation analysis. Through numerical simulations, we observe the occurrence of steady or timeperiodic pattern formations.


Introduction
Chemotaxis is a directed movement of biological individuals in response to a chemical signal. It is well-known that the chemotaxis can significantly affect the immune system (see [1,2]). To be more specific, during the immune process, a chemical which is called chemokine is secreted by the immune cell at sites of inflammation. Then, eukaryotic cells sense the gradient of chemokines by a polarized distribution of receptors, and they move toward relatively high concentrations of the chemical (chemoattraction), or in the opposite direction (chemorepulsion). Thus, chemotaxis seems to be a key contribution to the spatial properties of the solutions to immune systems, such as pattern formations or stabilization.
The immune system is regulated by various immune cells and dysregulation of the immune system is associated with some kinds of diseases. Immunodeficiency occurs when the immune response is less active, whereas hypersensitivity to the immune response is associated with many allergic diseases (see [3,4]). There have been several results to analyze the dynamics by using a mathematical model (see [5,6,7,8,9]). However, due to many variables and their complex coupling, it is not easy to understand the mechanism mathematically.
In this regard, Lee et al. [10] proposed a simplified PDE system involving only three main factors; antigens, chemokines and immune cells. The modelling is motivated by the pioneer work of Keller and Segel [11], which studies the aggregation of cellular slime mold toward a higher concentration of a chemical signal. Lee et al. obtained the stability and instability conditions for the uniform steady states and provided numerical observations. In this paper, we investigate the following reaction-diffusionadvection system introduced in [ x ∈ Ω, t > 0, C t = D C ∆C + s C MA − µ C C, x ∈ Ω, t > 0, x ∈ Ω, t > 0, A(x, 0) = A 0 (x), C(x, 0) = C 0 (x), M(x, 0) = M 0 (x), x ∈ Ω, where Ω is a smoothly bounded domain in R n , n = 1, 2. The unknowns A, C, and M represent the concentration of antigens, the concentration of chemokines, and the density of immune cells, respectively. Each equation of (1.1) contains a uniform diffusion with a constant coefficient and natural degradation with a constant decay rate. The constant antigen source s A > 0 represents the persistent infection. Since immune cells are produced in bone marrow and circulated throughout the body, we also consider the source of M as a positive constant s M > 0. During the immune process, immune cells and antigens are consumed by regulation and phagocytosis, which are denoted by λ M MA and λ A MA, respectively. On the other hand, s C MA indicates that chemokines are secreted by immune cells depending on the mass of antigens. The chemosensitivity χ > 0 in the advection term contributes to the movement of immune cells toward the gradient of chemokines. For the reader's convenience, we provide Figure 1 introduced in [10], which shows interactions between A, C and M. We supplement this system with the conditions of nonnegative initial values such that for p > n, A 0 (x), M 0 (x) ∈ W 1,p (Ω), C 0 (x) ∈ W 1,∞ (Ω), x ∈ Ω, and with the homogeneous Neumann boundary conditions where ν denotes the outward normal vector to the boundary ∂Ω.
The aim of this work is to verify the global solvability of (1.1)-(1.3) and show the occurrence of pattern formations. The rest of this paper is organized as follows. Section 2 is devoted to show the global existence and uniqueness of classical solutions. In the one-dimensional setting, the existence can be shown without any restriction on χ or initial conditions. However, in the case of two dimensions, due to the similarity to the minimal Keller-Segel model, the existence can be obtained under a restriction that depends on the initial conditions and coefficients. One of the main theorems (Theorem 2.1) will be followed by some Lemmas and a priori estimates. In Section 3, we present the stability of nonconstant steady states and existence of periodic orbits. Numerical results are given in Section 4. We close the paper in Section 5 with discussions and conclusions.
Throughout this paper, c denotes a positive generic constant which may differ from line to line.

Global existence
In this section, we prove the existence of classical solutions to (1.1). Looking at the first equation of (1.1), we can see that the nonlinear phagocytosis effect −λ A MA, due to its negative sign, does not play any role in the blow-up of A and thus A is uniformly bounded in time unless s A is singular. From the second equation of (1.1), the production of chemokines can not exceed a linear proportional to M, and we thereby find a similar structure to the Keller-Segel model. Similarly to −λ A MA, the nonlinear regulation −λ M MA in the third equation of (1.1) can not be a crucial obstacle to guarantee the boundedness of M. In general, the qualitative properties of the solutions to the minimal Keller-Segel model like boundedness or blow-up phenomena highly depend on the spatial dimension (see [12,13,14,15,16,17]). For this reason, we shall provide the proofs for each dimension.

Preliminaries
First, we introduce some useful estimates related to the heat semigroup (see [18,17]). Lemma 2.1. Let {e ∆t } t≥0 be the heat semigroup in a bounded domain with the homogeneous Neumann boundary condition and let λ 1 > 0 be the first positive eigenvalue of the Laplace operator −∆. Then, for any 1 ≤ q ≤ p ≤ ∞, there exists c = c(p, q, Ω) > 0 such that for all w ∈ L q (Ω); for all w ∈ L q (Ω).
Next, we assert the local existence of solutions to (1.1) and boundedness of L 1 -norm. Due to the presence of reaction terms, the total mass is not conservative. But we can still obtain the uniform boundedness of the total mass in time.

Lemma 2.2.
Let Ω be a smooth bounded domain of R n , n ≥ 1. For any (A 0 , C 0 , M 0 ) satisfying (1.2), we have: (1) There exists a maximal time of existence T max > 0 such that the system (1.1) with (1.3) has a unique nonnegative classical solution satisfying (3) There exists a positive constant C 1 such that In particular, we have Proof. We let U = (M, C, A) and rewrite (1 x ∈ Ω, The real parts of eigenvalues of A(U) are all positive and hence (2.8) is uniformly parabolic. Then, by Amman's well-established parabolic theory introduced in [19, Theorem 7.3] (see also see [20,21] for similar approaches), we obtain the local existence, uniqueness of solutions fulfilling (2.3) and blow-up criterion (2.4). The maximum principle ensures the nonnegativity of the solution. More precisely, we observe that the first and third component of F(U) are nonnegative, due to s M , s A > 0, when M = 0, A = 0. Then, the nonnegative initial values (A 0 , M 0 ) imply that (A, M) is also nonnegative. From the nonnegativity of (A, M) and the initial value C 0 , we obtain the nonnegativity of C. To show (2.5), we multiply the first equation, the third equation of (1.1) by s C and the second equation by λ A + λ M , respectively, and add them up. Then, integrating the result over Ω gives Letting s C Ω A(·, t) + (λ A + λ M ) Ω C(·, t) + s C Ω M(·, t) := y(t), we infer from (2.9) that there exist positive constants c 1 , c 2 such that which guarantees the uniform boundedness (2.5). As to (2.6), we integrate the third equation of (1.1) over Ω to obtain Solving the ordinary differential inequality (2.10) with M 0 , we obtain (2.6). For the last, we consider the following parabolic auxiliary problem for the equation of A (2.11) Now, we apply (2.1) to the Duhamel formulation of (2.11) as Then, we have By the standard comparison principle for parabolic equation (see [22]), we obtain 0 ≤ A(·, t) ≤Ã(·, t) for t ∈ (0, T max ).
This completes the proof.

One dimensional case
In this subsection, we consider an one-dimensional problem. In this case, gradient estimate of C is essential to show the global existence.
where the positive constant c is independent of T max .
Proof. We apply the estimates of Lemma 2.1 to the Duhamel formulation of second equation of (1.1) as the following equation (2.14) Using (2.1), (2.6) and (2.7), we deduce from (2.14) that for any t ∈ (0, T max ) where c is independent of t and λ 1 is the first positive Neumann eigenvalue of the Laplace operator −∆ on Ω. Since the two integral terms in (2.15) are uniformly bounded in t, this completes the proof.
The boundedness of M L p follows from (2.13).
where the positive constant c depends on p.
Proof. Testing the last equation of (1.1) by M p−1 for any p > 1, we have (2.17) Using (2.13), the Gagliardo-Nirenberg inequality, and Hölder's inequality, we compute where c 3 and c 4 are positive constants independent of p. Combining (2.17) and (2.18), we see that where c 5 = c 4 + 1. Note from the Gagliardo-Nirenberg inequality and (2.6) that holds for some positive constants c 6 , c 7 , c 8 and where the positive constant c depends on p.
By the application of bootstrap arguments, we prove the uniform boundedness of A, C and M as desired.
where the positive constant c is independent of T max .
Proof. Let 0 < T < T max . The boudnedness of A for t ∈ (0, T ) is a direct consequence of (2.7). As to the estimate of C, applying the estimates of Lemma 2.1 to (2.14) with (2.7) and (2.16), we have for where c is independent of T . Now, we apply the Moser-Alikakos iterative technique [23] to show the boundedness of M. Adding p 3 (p − 1) Ω M p dx to the both sides of (2.19), we have where c 10 = c 5 + 1. Interpolation inequality entails where c 11 is a positive constant independent of p. Therefore, (2.21) turns into where c 12 is a positive constant independent of p. We define Thus, we rewrite (2.23) as where c 13 is a positive constant independent of p.
where c 14 is independent of t. Since T < T max is arbitrary, we complete the proof.

Two dimensional case
In this subsection, we deal with the case of two spatial dimensions. First, we derive the uniform boundedness of M log M L 1 (Ω) and ∇C L 2 (Ω) under a smallness condition on M 1 . When the total mass of M is conserved, a Lyapunov functional is a helpful tool to prove the global existence (see [15]). However, due to reaction terms, the construction of a Lyapunov functional for this problem is very challenging. Instead, we use the Poincare inequality and the Cauchy-Schwarz inequality. Proof. Testing the third equation of (1.1) by log M, the integration by parts implies (2.25) where we used Young's inequality. Using (2.7), we find that there exists a positive constant c 15 > 0 such that On the other hand, we test the second equation of (1.1) by −∆C to obtain (2.28) Adding (2µ C − µ M ) Ω M log M to the both sides of (2.28), we obtain For the other case µ M > 2µ C , due to the positive uniform bound of −x log x for x > 0, we find a positive constant c 16 such that Now letting M ave (t) := 1 |Ω| Ω M(·, t), the Poincare inequality gives where c Ω is a positive constant satisfying By the Cauchy-Schwarz inequality, we deduce from (2.32) that Solving this inequality directly, we completes the proof.
Thanks to the work done by Bellomo et al. [24], the boundedness of Ω M log M and Ω |∇C| 2 are sufficient to obtain the global existence. Proof. Suppose that T max < ∞. According to Lemma 2.5 and Lemma 2.7, there exists c > 0 such that where c is independent of T max . This is contrary to Theorem 2.2 which implies T max = ∞. We complete the proof.

Stability and bifurcation
In this section, we analyze the bifurcation of the system (1.1). For simplicity, we consider the one-dimensional case : where the domain length L is a positive constant. Recalling the work in [10], we find that the system (1.1) has the unique positive equilibrium (A, C, M) such that where We the following functions spaces: From the equilibrium (3.2), we can rewrite the system (3.1) in the following abstract form: To obtain the principle of exchanging stability, we focus on the linearized system of (3.3) and then the linearized system (3.4) is equivalent to the following eigenvalue problem: where M k is a 3 × 3 matrix such that The principle of exchanging stability is determined by the sign of eigenvalues of each M k , and we have the following proposition (see [10, Lemma 1]). (3.8) Proof. Let Λ 1 k , Λ 2 k , and Λ 3 k be eigenvalues of M k . By the Routh-Hurwitz criterion, ReΛ i k < 0 for all k = 1, 2, 3, · · · and j = 1, 2, 3 if and only if Then, one can easily check that (3.9) completes the proof. The details can be found in [10, Lemma 1].

14)
where ξ 1 k 0 is given as in (3.13). And if K 3 ≤ 0, then there are no bifurcated solutions.
Proof. We use a similar method introduced in [25, Theorem 4.2] based on the center manifold theory (see e.g., [26]). Let Φ be the center manifold function of the equation (3.3) at χ = χ c . Since M k 0 | χ=χ c has a simple zero eigenvalue with the eigenvector ξ 1 k 0 in (3.14), we can express z as and the equation (3.3) is decomposed as the following: where P 1 : H → span{z 1 k 0 } ⊂ H and P 2 : H → span{z 1 k 0 } ⊥ are the natural projections. Then the bifurcation type of the equation (3.3) completely depends on the first equation of (3.15). Let z 1 * k 0 = ξ 1 * k 0 cos k 0 πx L , where ξ 1 * k 0 is the eigenvector of M * k 0 |χ = χ c corresponding to the zero eigenvalue. Taking the H-inner product to the both sides of the first equation of (3.15) by z 1 * k 0 , we can get the following equation: where (·, ·) denotes the H-inner product. Since we consider χ near χ S k 0 , G χ can be approximated by [26]), we can express g as g(y) = K 2 y 2 − K 3 y 3 + o(y 3 ).

Direct calculations yield that
, which has the factor L 0 cos 3 k 0 πx L dx = 0. To calculate K 3 , we let and, with the aid of Theorem A.1.1 in [26], Therefore, ψ satisfies (3.12) and we have which is calculated as in (3.11). Hence, the system (3.1) has a pitchfork type bifurcation and we complete the proof as desired.

Hopf bifurcation
In this subsection, we consider the case χ c = min k≥1 χ H k < min k≥1 χ S k , which may give rise to the existence of periodic solutions. Let k 1 be the positive integer satisfying χ k 1 = χ c . Then, there exists at least one eigenmode k such that χ H k < χ S k . Recalling from Lemma 1 in [10], we calculate the characteristic polynomial of M k as Λ 3 + P k Λ 2 + (Q k,0 − χQ k,1 )Λ + (R k,0 − χR k,1 ) with P k , Q k,0 , Q k,1 , R k,0 and R k,1 given in (3.8). By the Routh-Hurwitz criterion, for a such mode k with χ H k < χ S k and the choice of χ = χ H k , M k has two pure imaginary and one negative eigenvalues: To obtain the existence of the Hopf bifurcation, we apply the bifurcation theory from [27] and we check the transversality condition. The similar proof can be found in [28, such that (z k (s), χ k (s)) is a nontrivial solution of (3.3) and z k (s) is periodic with period T k (s) ≈ 2π ζ k where ζ k = Q k,0 − χ H k Q k,1 and ξ 1,2 k are the eigenvectors of M k | χ=χ H k corresponding to Λ 1,2 k (χ H k ). Moreover, ρ k (s 1 ) ρ k (s 2 ) for all s 1 s 2 , s 1 , s 2 ∈ (−δ, δ) and all nontrivial periodic solutions of (3.3) around (0, χ H k ) must be on the orbit ρ k (s) in the sense that, if (3.3) has a periodic solutionz(x, t) with period T for some χ ∈ R around ρ k (s) such that for small > 0 |χ − χ H k | < , |T − 2π/ζ k | < and max t>0,x∈[0,L] |z(x, t)| < , then there exist numbers s ∈ (−δ, δ) and some θ > 0 such that (T, χ) = (T k (s), χ k (s)) andz(x, t) = z k (s, x, t + θ).

Numerical simulations
In this section, we provide numerical examples of stable and time-periodic solution to the system (3.1) for the one dimensional case. We use the upwind finite volume scheme introduced in [29,30]. The parameter sets are given in Table 1 introduced in [10]. Since the model is simplified, it is not easy to find all the parameters estimated from the experiments. The diffusion coefficients and decay/death rates are usually known (see [31,9] and references therein). For example, the diffusion coefficient D M is estimated in a range of 10 −11 -10 −7 cm 2 /h (see [32]). In this work, we choose D M = 5 × 10 −7 cm 2 /h in the range. We choose domain length L = (half of some integers) × 0.005m, which is relevant to contain a sufficient number of cells considering the size of the cells. Table 1. Parameter values used for numerical simulations. Values of some groups (diffusion coefficients and decay/death rates) are from typical biological ranges [31,9]

Nonconstant stable steady states
By direct computations with (3.2) and Proposition 3.1, we obtain a positive constant equilibrium (A, C, M) ≈ (1.86 × 10 −2 , 1.93 × 10 −5 , 2.59 × 10 −1 ) and Q k,0 R k,1 > Q k,1 R k,0 for all k ∈ N, which implies that χ S k < χ H k for all k ∈ N.  Table 1. We choose the domain size L = 0.005 and the initial condition First we observe the chemotactic effect in a fixed domain length L = 0.005 with varying χ. In this case, we have min k∈N χ S k = χ S 2 ≈ 1.097. Thus, if we choose χ > 1.097, then the system bifurcates from the equilibrium to another stable steady state solution. In Figure 2, we illustrate steady state solutions of A, C and M with χ taken to be 1.15, 1.6 and 2. We see that A has an interior spike, and C and M have boundary layers on each boundary. As χ increases, the solutions become more aggregated. In Figure 3, there are two steady state solutions with different structures. Theorem 3.1 indicates that the system (3.1) has two bifurcated stable steady solutions in case of χ c = min k∈N χ S k < min k∈N χ S k . Figure 3 shows that the structure of bifurcated solutions depends on the initial condition. Figure 4 illustrates the simulation in several cases of domain lengths, L = 0.0025, 0.005, 0.01, 0.015, and for each case, χ = 1.2, 1.15, 1.07, 1.065 are chosen. We see from Table 2 that the stable steady solutions have the profile cos k 0 πx L when χ is slightly bigger than χ S k 0 = min k∈N χ S k . These numerical results in Figure 3 and 4 agree with the theoretical results of Theorem 3.1. Table 2. k 0 and χ S k 0 corresponding to L where k 0 = argmin k∈N χ S k . L 0.025 0.005 0.01 0.015 k 0 1 2 3 5 χ S k 0 1.096 1.097 1.060 1.062 Figure 4. Space-time color maps (first row) and stable steady states (second row) of the antigen concentration A in the case of L = 0.005, 0.01, 0.015 respectively. In each case, we choose χ = 1.15, 1.07, 1.065, respectively, which are slightly bigger than χ S k 0 = min k∈N χ S k . The stable steady states have spatial profile cos k 0 πx L .

Time-periodic solutions
From Table 1, we find the positive constant equilibrium (A, C, M) = (2.00 × 10 −2 , 2.00 × 10 −5 , 2.50 × 10 −1 ) and min k∈N χ H k ≈ 1.04 < min k∈N χ S k ≈ 14. In this case, by Theorem (3.2), the system (3.1) posseses time-periodic solutions. Figure 5 illustrates time-periodic solutions to (3.1) with parameters of Table 1 and L = 0.005 when χ = 1.085 is chosen slightly bigger than min k∈N χ H k = χ H 2 ≈ 1.08. Letting ζ k (χ) = Q k,0 − χQ k,1 , we compute the period of the solution as T = 2π ζ k (χ H 2 ) ≈ 49.2. Figure 6 shows the spatial profile of time periodic solutions featuring cos k 1 πx L , where k 1 = argmin k∈N χ H k for varying domain size L = 0.0025, 0.005, 0.01, 0.0125. Then, k 1 and χ H k 1 = min k∈N χ H k are given in Table 3.   We observe that the solution of each case is time-periodic and has spatial profile cos k 1 πx L .

Conclusions and discussions
In this paper, we investigate the global solvability and bifurcation dynamics of the system (1.1) with the boundary condition (1.3) introduced in [10]. The system (1.1) has a minimal type of chemotactic sensitivity, i.e., constant χ = const. Due to a similar structure to the classical Keller-Segel model (see [11]), the global solvability of the problem depends on the dimension. In the one dimensional case, the global existence and uniform boundedness of solutions can be proved for any χ > 0. However, in the two dimensional case, the global solvability of this model is obtained under smallness assumptions (see [15] for Keller-Segel model case).
In view of the stability analysis, it is known that the system has a single positive constant equilibrium (A, C, M) described in (3.2) and there exists a threshold χ c = min k∈N {χ S k , χ H k } of χ such that the equilibrium (A, C, M) is asymptotically stable when χ < χ c and unstable when χ > χ c (see [10]). This study goes further from the work in [10] by means of local bifurcation theory. As χ increases passing through the threshold value χ c , the equilibrium (A, C, M) becomes unstable and two different type bifurcated solutions can occur; If χ S k 0 < χ H k for all k ∈ N, then there exist nonconstant steady states whose structure depends on initial values, and if χ H k 1 < χ S k for all k ∈ N, then we obtain time-periodic solutions. Moreover, we obtain the information of wave modes of bifurcated solutions in both cases. However, in Hopf bifurcation case, the information of wave modes is local which means that the solution can have different dynamics from the form of (3.18) when χ is much bigger than χ c = χ H k 1 . We see this dynamics, via the example illustrated in Figure 7. The global bifurcation will be investigated in the future.