A control problem for a cross-diffusion system in a nonhomogeneous medium†

We study a model, motivated by a bioremediation process, describing a cross-diffusion movement of a bacteria population b attracted by a chemoattractant signal c, in a nonhomogeneous stratified medium with n layers. We assume that this reaction–diffusion process is characterized by a low rate of degradation and a low diffusion coefficient of the chemoattractant, expressed in the model by a small parameter ∊. The model consists of n systems of nonlinear parabolic equations with transmission conditions between layers. We prove a global-in-time solution for the asymptotic model setup with respect to the small parameter of the problem, for arbitrarily large initial data. Next, we deal with the control problem focusing mainly on the reduction of the chemoattractant concentration, by acting upon the initial distribution of the bacteria population b0. To this end, we prove the existence of a solution to the control problem and determine the optimality conditions. AMS Subject Classifications: 35K57; 30E25; 35B20; 35K60


Introduction
Chemotaxis is the biological process of directed movement of cells in response to a chemical signal emitted by a substance or by another population, called chemoattractant, and play an important role in the interaction of cells with their environment. Chemotaxis is a complex process which involves many aspects such as various species of cell populations b can interact; the movement can be either towards the higher concentration of the signal (positive chemotaxis) or away from it (negative chemotaxis); the chemical signal can be secreted by a population b itself, and not necessarily by an external source; the chemotactic process can lead to an aggregation of the attracted individuals and to a production of the chemoattractant, or contrary, it can determine the degradation of the chemoattractant. Consequently, the theoretical understanding of these processes determined a growing interest in their mathematical modelling, which at its turn raised challenging problems.
The origin of the fundamental model is given in the work of Patlak [24]. Later, Keller and Segel [18][19][20] introduced a similar model based on another assumptions. Since then, a rich mathematical literature on various versions of the model has been emerged, mainly focusing on the wellposedness of it and we refer the reader to a very comprehensive survey in [15].
In this paper, we shall denote the density of the cell population by b and the density of the population spreading the signal by c, both of them being functions of time t, and space variable ξ .
A chemotactic system with only one population of cells consists of two equations for b and c with initial and boundary conditions: . They may depend on the space variable, too. The function K characterizes the chemotactic sensitivity and the cross-diffusion term in the first equation is indeed able to enforce the spontaneous emergence of structures provided that the process of chemotactic migration is accompanied by a production of the signal substance by the cells themselves [14]. Thus, the cross-diffusion term and the kinetic term ϕ(b, c) in the second equation can lead to the blow-up of the solution, even if the growth-death rate f (b, c) is zero and D, δ, K are constant. In the literature, the chemotactic system has been approached in simplified versions, able to avoid blow-up and to allow global solutions. We cite some more recent results.
In the 1-D case, it has been shown that blow-up does not occur [16] for D = 1, K constant and δ = ε small. When the space dimension d is greater or equal to two, the solutions generally exhibit blow-up, this being influenced by the model parameters and the characteristics of the initial data [13,14,25]. For example, in [17], a chemotaxis motion with constant diffusion coefficients is studied by using a nonlocal gradient sensing term to model the effective sampling radius of the species. In [9], Dyson et al. use a nonlocal term to model the species-induced production of the chemoattractant, ϕ 2 (b, c, ξ), in order to prevent blow-up in the d-D space, considering that the diffusion coefficients are constant. They prove the existence of solutions, which exist globally, and are L ∞ -bounded on finite time intervals. In [28], the system (1)-(4) is considered with and it is shown that for arbitrarily large initial data, this problem admits at least one global weak solution for which there exists T > 0, such that (b, c) is bounded and smooth in (T , ∞) × . The paper [8] deals with two types of reaction-diffusion systems, one arising in chemotaxis and the other in angiogenesis. The first refers to Equations (1)-(4) with D(b, c), δ(b, c) and K(b, c) constant, f (b, c) = 0, ϕ(b, c) = b − αc and the equation for c stationary, that is, − c = b − αc. This equation is obtained as a case limit, for , and the diffusion coefficient δ is of the order of 1/ε, with ε small. For this system (as well as for that of angiogenesis type), it is shown in [8] that when the L d/2 norm of initial data is small enough (for d ≥ 2), then there is a global (in time) weak solution that stays in all L p spaces with max{1; d/2 − 1} ≤ p < ∞. In [4], the same system but with α = 0 is studied in R 2 and a detailed proof of the existence of weak solutions below the critical mass, above which any solution blows up in finite time in the whole Euclidean space, is given [26]. The stability of the stationary solutions to a chemotaxis system was proved in [11] for D = 1, δ(b, c) = 0 and a general function ϕ(b, c). In [23], existence and uniqueness are studied for slow and singular fast diffusion of the cells in the case with a stationary equation for the chemoattractant.
Besides applications in biology, a chemotaxis model can apply in environment bioremediation in the cases when specific bacteria b are injected into a polluted medium (soil or water) with the purpose of cleaning it from an inside spread pollutant c [5,10,27].
Our study is motivated by an application to environment bioremediation and focuses on the case in which the kinetic term and the diffusion coefficient of the chemoattractant (pollutant) c have a weak influence on the flow, meaning that the rate of degradation of the chemoattractant is slow and it diffuses very little (or not at all, as in the case of oil polluting an environment).
Roughly speaking we shall start from a model reading as where ε is a small parameter in front of the diffusive and the kinetic terms for the chemoattractant. Such a model is obtained by making dimensionless Equations (1)-(4). Moreover, we assume that at the initial time t = 0 the chemoattractant concentration c 0 is constant and that there exists a certain nonlinear dynamics growth-death of the cell population, indicated by f . The model will be implemented into a nonhomogeneous medium, fact evidenced later by a particular space dependence of the coefficients. The aim is to investigate the way in which the solution depends on the small parameter ε and not to study the limit model when ε → 0. Accordingly, we shall not pass to the limit, but use a perturbation technique [7], by which the solution is expanded in series with respect to the powers of the small parameter, and retain the systems of ε m -order of approximation, obtained by equating the coefficients of ε m in Equations (5).
The main goal of the paper is to study the possibility of controlling the environment cleaning by acting upon the initial distribution of the bacteria. More exactly, the intention is to design the initial distribution of the bacteria such that the chemoattractant mean concentration should decrease under a certain critical value c crt within a time period T .
A secondary task of the control problem is motivated by the fact that bioremediation can be problematic because sometimes it might become difficult to remove the microorganisms released into the environment in order to let them at acceptable levels in the soil [21]. Consequently, for enhancing a complete and efficient decontamination process and for avoiding higher costs of a further cleaning of the environment from the remained bacteria, we investigate by the same control problem whether a moderate bacteria proliferation, whose growth is limited by a prescribed density value b c , would be sufficient for achieving the main objective of reducing the pollutant concentration up to a value under the dangerous critical threshold c crt . To this end, the state system is studied in the framework of an appropriate functional setting and a global-in-time solution for the asymptotic system derived from Equations (5) is obtained. Then, the existence of at least a solution to the control problem is proved and the optimality conditions proving information about the initial distribution of the bacteria b 0 are computed.

Statement of the model and the perturbation technique
We consider a nonhomogeneous 3D right cylinder domain where the base is an open bounded subset of R 2 , with a boundary of class C 2 . The medium nonhomogeneity is modelled by a stratification of in n parallel layers along the Ox axis, the separation of the layers being determined by the different values that movement parameters may have in each layer. More exactly, a stratification can be put into evidence when certain parameters do not depend on x (here x is the stratification variable) or are constant in each layer (x i−1 , x i ) but have different values from a layer to the other.
Therefore, the domain consists of n subdomains i , having the boundaries The surfaces 0 and n are the external horizontal boundaries, while i with i = 1, . . . , n − 1 are the boundaries between layers. We denote In each layer i the chemotaxis process is modelled by two equations, one for b i and the other for the chemoattractant c i . The interaction between the layers is established by transmission conditions for b i , that is, the continuity of the solutions and fluxes. We assume that the system is closed for b i and c i , namely the fluxes across the exterior frontiers are zero. The chemoattractant c i may display jumps at the interface between layers, because the ε m -approximation systems will not be anymore of diffusion type, as we shall see. Therefore, there is no need to specify these conditions. The stratification is set by the values for D i , δ i , c i,0 assumed constant in each layer i, but different for two consecutive layers. Also, the expressions of the functions f i , K i , ϕ i , which are assumed not to depend explicitly on ξ , are different from one layer to another.
With these considerations, the dimensionless mathematical model for a chemotaxis movement in a nonhomogeneous stratified medium has n equations for the unknowns b i and c i (the density of the cell population and the chemoattractant concentration, respectively, in each layer i = 1, . . . , n) and reads [1] for all i = 1, . . . , n, where c i,0 (assumed constant) and b i,0 (ξ ) are initial conditions for c i and b i . At the interface between two layers, we have the conditions for i = 1, . . . , n − 1, and on the exterior horizontal and lateral boundaries, we set Here, ν is the unit outer normal to lat i and ∂/∂ν is the normal derivative. Generally, c i and its normal derivative δ i (∂c i /∂x) can have a jump when crossing the internal boundaries i . As we said, do not specify here these conditions because they will not intervene in the asymptotic model we shall study.
We stress that this model is presented in a dimensionless form and the constantsD,K,f ,δ,φ are dimensionless parameters. They appear after the dimensionless procedure by which all dimensional variables and function are scaled by certain characteristic values indicated by the subscript a (e.g. L a for length, T a for time, D a for diffusion coefficient, etc.). We do not present the dimensionless computation, but we assert that by this procedure the dimensionless parameters are given byD If by such a procedure the computed valuesδ andφ follow to have the same order of magnitude, much smaller than that of the other parameters, we set in Equations (6)-(17) and consider the other dimensional parameters of order O(1) as ε → 0. We agree that this model characterizes the chemotactic process with a chemoattractant (e.g. oil in a bioremediation process) that diffuses very slow, and it is degraded with difficulty by the bacteria. This process happens in a different way in each layer, due to the nonhomogeneity modelled by the various functions f i , K i , ϕ i , the structure of the initial data and different diffusion coefficients.

Hypotheses
To approach system (6)-(17) we assume the following hypotheses for all i = 1, . . . , n : are of class C 2 with respect to r 1 and |ϕ i (r 1 , r 2 )|, |(∂ϕ i /∂r 1 )(r 1 , r 2 )|, |(∂ 2 ϕ i /∂r 2 1 )(r 1 , r 2 )| are bounded. We remark that equations with nonlinear terms f i do not generally admit global solutions in time [12], but by the perturbation technique we shall deduce a global solution for this asymptotic model, even under the assumption of a polynomial form for f i . Thus, for all i = 1, . . . , n, we assume and d is the space dimension.

ε 0 -order and ε 1 -order approximations
We write the series expansions of all functions, with respect to the small parameter ε =φ =δ.
where φ m i corresponds to the ε m -approximation of φ i (m = 0, 1, . . .), and We replace these series in the system (6)- (17) and by equating the coefficients of the mth powers of ε, we deduce the systems corresponding to ε m -order approximation.
The ε 0 -order approximation for c i reads for each i = 1, . . . , n, implying that c 0 i is a positive constant in each layer, For b 0 i we get Further, identifying the coefficients of ε 1 and using Equation (27), we obtain the system for the ε 1 -order approximation, where we have denoted ξ := (y, z) ∈ , for i = 1, . . . , n, for i = 1, . . . , n − 1, After solving the system for the ε 0 -order approximation, c 1 i follows by Equations (35) and (36) and so the functions F i (t, ξ), G i (t, x i , ξ ) are known.
The next approximations (for m ≥ 2) lead to systems having similar forms as that for the ε 1order approximation, and so they pose the same mathematical problems. That is why we do no longer write them.
Definition 2.1 We call an asymptotic solution to Equations (6)- (17), up to the order of approximation ε 2 , a pair of functions where c 0 i and c 1 i are the solutions to the differential equations (25) and (26) and (35)  We shall rigorously explain the definition of the solutions b 0 i and b 1 i in the next sections, in which the well-posedness of these systems will be studied.

The control problem
As we have already explained, the aim of the control problem is to provide information about the necessary initial density of bacteria, b 0 (ξ ) = (b i,0 (ξ )) i=1,...,n , such that two objectives would be achieved. The main one is related to the environment cleaning that is, to force the decrease of the chemoattractant mean concentration under a critical threshold c crt . This means to minimize in the cost functional the mean positive part of the difference between the concentratioñ c(t, ξ) = c 0 i (t, ξ) + εc 1 i (t, ξ) given by Equation (48) and c crt . The second objective aims at realizing the necessary pollutant concentration decrease by a process limiting a too large proliferation of the bacteria. Consequently, in the cost functional we add a term expressing the restriction of the bacteria growth, by minimizing the mean positive part of the difference between b 0 i (t, ξ) and a constant prescribed value b c . We introduce this control problem to give a quick response (rather than an extremely accurate one) and to indicate a first decision for a process evolving in real time. That is why we accept that the consideration of the ε 0 -order approximation b 0 i only, which is the dominant term in the asymptotic expansion of b i , is motivated. A more accurate computation, if necessary, may take into account the further approximations.
Mathematically, we have to minimize the cost functional where the superscript '+' means the positive part (i.e. ϕ + (ξ ) = max{ϕ(ξ ), 0}). A certain choice of the constant σ can induce a lower (or a higher) influence of the first term with respect to the second one in the cost functional. In particular, if the limiting condition of the proliferation of the bacteria is disregarded, the constant σ can be set zero.

G. Marinoschi
Therefore, the control problem is We recall that in our model, the initial distribution of the density of the bacteria may depend on the

Well-posedness of the state system
In order to approach the control problem we need some information about the ε 0 and ε 1 -state systems. We begin by a few preliminaries.

Preliminaries and functional setting
Since by our assumptions, c i,0 remains constant in each layer, one notes that in Equation (28) f i depends only on b 0 i and so one can denote The assumptions (19)- (23) imply that the functions μ i ∈ C 1 (R), with p ∈ [0, 2) if d = 3 and p ≥ 0 if d = 1 or d = 2. We rewrite the system (28)-(34) for the ε 0 -order approximation without indicating the superscript '0' We introduce now the global functions defined in the following way (indicated by a generic notation φ, φ 0 ): (for a function which is constant in each layer, as c 0 or D), function depending on b and c), where we recall that ξ = (x, y, z). We mention that we include the constantD in the definition of D(x).
We note that assumptions (50)-(52) and (i 1 )-(i 5 ) imply similar properties for the functions defined before. Namely, we have uniformly with respect to x, where C 1 =f max i=1,...,n {C i 1 }, and μ r is the derivative with respect to r, We recall that we consider p ∈ [0, 2) if d = 3 and p ≥ 0 if d = 1 or d = 2. System (53)-(59) will be treated in the functional framework of the Sobolev space V = H 1 ( ) endowed with the standard norm ψ V = ( ψ 2 + ∇ψ 2 ) 1/2 , and its dual V , with the pivot where ·, · V ,V represents the duality between V and V .

G. Marinoschi
We specify that for the writing simplicity we shall denote the scalar product and the norm in L 2 ( ) by (·, ·) and · .
We define the operator A 0 : V → V by and restrict it to L 2 ( ) by the operator A : Thus, the functional abstract setting of our problem is

Existence for the ε 0 -order approximation
The well-posedness for Equations (71) and (72) is concentrated in the following. The solution satisfies the estimate where the constant C V depends on the problem data. Moreover, Proof The proof is split into two steps and is given in [1]. Further, we indicate the arguments. We assert by the properties (62)-(64) and by handling appropriate Sobolev inequalities that the operator b −→ μ(b, x) turns out to be locally Lipschitz from V to L 2 ( ), uniformly in x. In the first step, according to the method presented in [3], we reduce the problem to a case with a globally Lipschitz operator, by approximating μ(·, x) by for each N natural. This truncated operator is Lipschitz from V to L 2 ( ) for each N fixed, and we consider the approximating problem where A N is the operator A (defined by Equation (70) Some further necessary estimates are deduced for this solution, that is, implying Equation (74), with C 0 depending on fixed data of the problem (p, T , D 0 , D ∞ ), and db N dt 2 L 2 (Q) The proof of the m-accretiveness of A N , as well as all estimates written before involve some technical computations given in detail in [1].
The positiveness is proved on the basis of Stampacchia's lemma and Gronwall's lemma and show that the solution falls within the accepted physical domain of positive densities.
The regularity H 2 (a.e. t) of the restriction of b to each layer b i comports a fine technique because it must adapt the known results [2,6] for a domain with a regular boundary ∂ to this transmission problem with n − 1 interfaces. These arguments are given in [1], too.
Finally, Equation (74) opens the ways to the second step by going back to Equations (71) and (72). In fact, if we choose R = C V and N large enough, N > R, we get that A N b N (t) = Ab N (t). Therefore, b N (t) with N large enough is the solution to problem (71) and (72), and further, it will be denoted by b.
For proving the uniqueness of the solution, we consider two solutions b andb corresponding to the same initial data b 0 . By the previous proof if is the solution to Equation (78) and (79).
It is obvious that in each layer the functions b i (meaning in fact the ε 0 -approximation b 0 i ) have the regularity induced by b restricted to Q i , for i = 1, . . . , n, that is,

Existence for the ε 1 -order approximation
The functions c 1 i are computed by Equations (35) and (36), and by the boundedness hypothesis (i 5 ) and Equations (83) and (75), we derive that Existence for the solution to Equations (37)-(43) for the ε 1 -order approximation b 1 i (t, ξ) is studied by the Lions' theorem for the time-dependent case [22] and the conclusions are [1] the following.
given byb In particular, the restrictions of the solution to each layer have the properties

The control problem
Using the notation (60) we can rewrite the control problem (P) as subject to Equations (53)-(59) (equivalently Equations (71) and (72)) and (84), where This means that b 0 is the solution to the ε 0 -approximation, Equations (53)-(59), with the initial datum b 0 and c is given by Equation (84).
..,n , and we note that U = n i=1 U i . We recall that c 0 (t, ξ) = c 0 and is constant in all layers. Because in (P) only the ε 0 -order approximation for b and only the ε 1 -order approximation for c are involved, for the writing simplicity we shall indicate them by b and c, without superscripts.

Existence of the optimal control
Theorem 4.1 Problem (P) has at least a solution.
Proof Since J(b 0 ) ≥ 0, its infimum exists and let us denote it by Let us take a minimizing sequence ..,n is the solution to Equations (53)-(59) with initial datum b k 0 and c k is given by Equation (84) corresponding to b k . We stress that the confusion between the sequence (b k ) k≥1 and the ε k -order approximation should be avoided.
Since b k 0 ∈ U we can select a subsequence (denoted still by k) such that Then, the Cauchy problem (71) 82) and (85). Therefore, selecting a subsequence we deduce that Then, b k → b * a.e. on Q and since ϕ is continuous it follows that a.e. on Q, implying that (c 0 + εc k − c crt ) + → (c 0 + εc * − c crt ) + a.e. in Q, as k → ∞. Moreover, ϕ is bounded by (i 5 ), so by the Lebesgue dominated convergence theorem, we get by Equations (84) that But, on the other hand, by Equation (89), the sequence (c 0 + εc k − c crt ) + k is bounded in L 2 (Q), so that it converges weakly on a subsequence to a limit in L 2 (Q), and in conclusion, by the limit uniqueness (c 0 + εc k − c crt ) + → (c 0 + εc * − c crt ) + weakly in L 2 (Q), as k → ∞.
Since (μ(b k )) k lies in a bounded subset of L 2 (Q) and μ(b k ) → μ(b * ) a.e. in Q, it follows that We also deduce by the Ascoli-Arzelà theorem that Moreover, b * is the solution to Equations (71) and (72) because, on the basis of the above convergencies, we can pass to the limit as k → ∞ in the weak form Finally, using the weakly lower semicontinuity property of the functions in Equation (89), we can pass to the limit and deduce which shows that J(b * 0 ) = d, hence (P) has at least a solution.

The system in variations and the dual system
Let (b * 0 , b * , c * ) be an optimal pair in (P). The function b * is the solution to Equations (71) and (72) with the initial datum b * 0 and c * is given by Equation (84) corresponding to b * , Let λ > 0 and introduce b λ The corresponding solution to Equations (71) and (72) with the initial datum b λ 0 is denoted by b λ , and we define where we recall that w(ξ ) = v(ξ ) − b * 0 (ξ ). Combining Equations (106) and (107) we finally deduce which can be still written In this relation, In conclusion, we obtain relations (105) as claimed.
We recall that b m and b M represent a sequence of minimum and maximum values assigned for each layer i. Equivalently, the previous relation −p(0, ξ) ∈ N U (b * 0 ) means that −p i (0, ξ) ∈ N U i (b * i,0 ) and so the relations given by relations (105) follow.
We observe that the minimization problem admits at least a solution which can be realized not necessarily at b 0 = b m , but in a way expressed by relations (105). We conclude that the controller b * 0 (= b * i,0 ) may take either the value b i,m or b i,M on the subsets of i , where p i (0, ξ) has positive or negative values. In some layers, b * i,0 may reach also the value 0 (if e.g. p i (0, ξ) > 0 in i ), because b i,m was allowed to be nonnegative. This means that in some situations it might be possible not to place bacteria in some layers at the initial time. Finally, we assert that it is difficult to say if the system (98)-(104) can have a null solution p i (0, ξ) on a subset of i , but this analysis is beyond the scope of this paper. Therefore, we cannot exclude that in some circumstances it may happen that the minimum be reached for whatever b i,0 ∈ (b i,m , b i,M ).

Numerical results
Some graphics revealing the feature of the process with respect to the change of b 0 are presented. The simulations are made with Comsol Multiphysics (FLN License 1025226) for the system (6)- (17) in a 1D domain with the same values D, δ, χ , c 0 in all layers (in order to compare how results change for various b 0 ) with the following dimensionless data: In Figure 1 (right), we see that the values of c increase from x = 0 to x = 10. The medium cleaning is more efficient beginning with the first layers. In Figure 2 (right), the medium cleaning happens in a different way, being less efficient in the first layers (x ∈ [0, 2]) where the values of c remain high, but more efficient in the medium layers (x ∈ [4,6]).

Conclusion
We address a control problem related to a chemotactic motion of a bacteria towards a pollutant chemoattractant in a stratified medium. The control problem is focused on the reduction of the chemoattractant concentration, by acting upon the initial distribution of the bacteria population, b 0 . We prove that the control problem has at least a solution and provide the structure of the controller b 0 in each layer.