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

Abstract

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 b₀. To this end, we prove the existence of a solution to the control problem and determine the optimality conditions.

Keywords:

• reaction–diffusion system
• chemotaxis
• perturbation method
• boundary value problems for nonlinear parabolic PDE

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–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 well-posedness 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:

Here, the time runs in (0, T), T is finite, which is an open bounded subset of ℝ^d and . In these equations, D(b, c) and δ (b, c) represent the diffusion coefficients of the attracted population b and chemoattractant c, respectively; g(b, c) and h(b, c) are functions describing the rates of growth and death of b; and and stand for the production and degradation of the chemoattractant. Further, we denote the right-hand sides by and 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 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, 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 f(b, c)=0 and 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 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, and the equation for c stationary, that is, This equation is obtained as a case limit, for of Equation (3) in which and the diffusion coefficient δ is of the order of 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 . In [4], the same system but with α=0 is studied in ℝ² 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, and a general function 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₀ 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 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 -order of approximation, obtained by equating the coefficients of 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₀ are computed.

2. 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 ℝ², with a boundary of class C². 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 i=1, …, n, where are the lateral boundary of Ω_i and , i=0, …, n. The surfaces Γ₀ and Γ_n are the external horizontal boundaries, while Γ_i with 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 -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 , and on the exterior horizontal and lateral boundaries, we set

Here, ν is the unit outer normal to and is the normal derivative.

Generally, c_i and its normal derivative 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 constants 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 by 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 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.

2.1 Hypotheses

To approach system (6)–(17) we assume the following hypotheses for all i=1, …, n:

• (i₁) c_{i, 0} is constant, c_{i, 0}≥0, and there exist an i such that c_{i, 0}>0;

• (i₂) b_{i, 0}≥0, and there exist an i such that b_{i, 0}>0;

• (i₃) , ;

• (i₄)  are of class C¹ with respect to r₁, are bounded, ;

• (i₅)  are of class C² with respect to r₁ and 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

with positive for i=1, …, n, where and d is the space dimension.

2.2 ϵ⁰-order and ϵ¹-order approximations

We write the series expansions of all functions, with respect to the small parameter Denoting generically and we have

where corresponds to the -approximation of φ_i (m=0, 1, …), and , represent the derivatives of Φ_i with respect to b_i and c_i.

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 -order approximation.

The -order approximation for c_i reads

for each i=1, …, n, implying that is a positive constant in each layer,

For we get

Further, identifying the coefficients of and using Equation (27), we obtain the system for the -order approximation,

where we have denoted

for i=1, …, n,

for i=1, …, n−1,

After solving the system for the -order approximation, follows by Equations (35) and (36) and so the functions F_i(t, ξ), are known.

The next approximations (for m≥2) lead to systems having similar forms as that for the -order 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 a pair of functions i=1, …, n,

where and are the solutions to the differential equations (25) and (26) and (35) and (36), respectively, and and are the solutions in the sense of distributions to Equations (28)–(34) and (37)–(43).

We shall rigorously explain the definition of the solutions and in the next sections, in which the well-posedness of these systems will be studied.

2.3 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, 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 concentration 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 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 -order approximation 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. . 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.

Therefore, the control problem is subject to Equations (28)–(34) and (35)–(36), where We recall that in our model, the initial distribution of the density of the bacteria may depend on the space variable ξ (i.e. and only the chemoattractant initial concentration is considered constant (. A natural restriction for the initial densities of the bacteria b_{i, 0} ( representing the controllers in the control problem, is to impose their boundedness. The minimum and maximum bounds, b_{i, m} and b_{i, M}, are fixed numbers, . However, as we see later, this restriction does not necessarily imply that the minimum of J must be reached at b_{i, m}.

Finally, to be consistent with the requirement expressed by the first term in J, the maximum of the initial datum is chosen less than or equal to b_s, the appropriate upper threshold envisaged for the solution In the case with σ=0, the choice of b_{i, M} is free.

3. Well-posedness of the state system

In order to approach the control problem we need some information about the and -state systems. We begin by a few preliminaries.

3.1 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 and so one can denote

The assumptions (19)–(23) imply that the functions

with p∈[0, 2) if d=3 and p≥0 if d=1 or d=2.

We rewrite the system (28)–(34) for the -order approximation without indicating the superscript ‘0’

We introduce now the global functions defined in the following way (indicated by a generic notation φ, φ₀):

(for a function which is constant in each layer, as c₀ or D),

(for a function depending on b and c), where we recall that We mention that we include the constant in the definition of D(x).

We note that assumptions (50)–(52) and (i₁)–(i₅) imply similar properties for the functions defined before. Namely, we have

uniformly with respect to x, where and 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 endowed with the standard norm and its dual , with the pivot , such that . The value of at ψ∈V is

where represents the duality between V^′ and V.

We specify that for the writing simplicity we shall denote the scalar product and the norm in L²(Ω) by and .

We define the operator by

and restrict it to L²(Ω) by the operator

with .

Thus, the functional abstract setting of our problem is

3.2 Existence for the ϵ⁰-order approximation

The well-posedness for Equations (71) and (72) is concentrated in the following.

Theorem 3.1

Assume conditions (62)–(64), and let b₀∈D(A). Then, problem {\rm (71)} and {\rm (72)} has a unique solution

The solution satisfies the estimate

where the constant C_V depends on the problem data. Moreover,

Finally, let b₀≥0 a.e. in Ω. Then, the solution b to Equations {\rm (71)} and {\rm (72)} satisfies

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 turns out to be locally Lipschitz from V to L²(Ω), 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 by

for each N natural. This truncated operator is Lipschitz from V to L²(Ω) for each N fixed, and we consider the approximating problem

where A_N is the operator A (defined by Equation (70)) in which μ (b, x) is replaced by μ_N(b, x) and D(A_N)=D(A). For each N, it can be proved that A_N is quasi-m-accretive on L²(Ω) and then, for , problem (78) and (79) has a unique strong solution Some further necessary estimates are deduced for this solution, that is,

implying Equation (74), with C₀ depending on fixed data of the problem and

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² (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 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 and corresponding to the same initial data b₀. By the previous proof if then we get b(t)=b_N(t) and , where b_N(t) is the solution to Equation (78) and (79).

It is obvious that in each layer the functions b_i (meaning in fact the -approximation have the regularity induced by b restricted to Q_i, for i=1, …, n, that is,

3.3 Existence for the ϵ¹-order approximation

The functions are computed by Equations (35) and (36),

and by the boundedness hypothesis (i₅) and Equations (83) and (75), we derive that

Existence for the solution to Equations (37)–(43) for the -order approximation is studied by the Lions’ theorem for the time-dependent case [22] and the conclusions are [1] the following.

Theorem 3.2

System (37)–(43) has a unique solution

Corollary 3.3

Problem (6)–(17) admits a unique asymptotic solution up to the order of approximation

given by

In particular, the restrictions of the solution to each layer have the properties

4. 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⁰ is the solution to the -approximation, Equations (53)–(59), with the initial datum b₀ and c is given by Equation (84). Here, and , and we note that .

We recall that and is constant in all layers. Because in (P) only the -order approximation for b and only the -order approximation for c are involved, for the writing simplicity we shall indicate them by b and c, without superscripts.

4.1 Existence of the optimal control

Theorem 4.1

Problem (P) has at least a solution.

Proof Since J(b₀)≥0, its infimum exists and let us denote it by Let us take a minimizing sequence , such that

where is the solution to Equations (53)–(59) with initial datum and c^k is given by Equation (84) corresponding to b^k. We stress that the confusion between the sequence and the -order approximation should be avoided. Since we can select a subsequence (denoted still by k) such that Then, the Cauchy problem (71) and (72) has, by Theorem 3.1, a unique solution b^k, satisfying Equations (80)–(82) and (85). Therefore, selecting a subsequence we deduce that which implies by Lions’ theorem (since V is compact in that Then, a.e. on Q and since ϕ is continuous it follows that a.e. on Q, implying that a.e. in Q, as Moreover, ϕ is bounded by (\textbf i₅), so by the Lebesgue dominated convergence theorem, we get by Equations (84) that But, on the other hand, by Equation (89), the sequence is bounded in L²(Q), so that it converges weakly on a subsequence to a limit in L²(Q), and in conclusion, by the limit uniqueness Since lies in a bounded subset of L²(Q) and a.e. in Q, it follows that We also deduce by the Ascoli–Arzelà theorem that which implies that as 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 in the weak form getting Finally, using the weakly lower semicontinuity property of the functions in Equation (89), we can pass to the limit and deduce which shows that hence (P) has at least a solution.

4.2 The system in variations and the dual system

Let be an optimal pair in (P). The function b* is the solution to Equations (71) and (72) with the initial datum and c* is given by Equation (84) corresponding to b*,

Let λ>0 and introduce where The corresponding solution to Equations (71) and (72) with the initial datum is denoted by b^λ, and we define By a few computation, it can be shown that Y solves the system in variations (written for Y_i in each layer)

where with v_i∈U_i and We see that the system in variations can be treated in the same manner as the system for the -order approximation (being much simpler than it), and so it follows that it has a unique solution

We introduce now the system for the dual variable p, written for each layer,

It is obvious that it has a unique solution

4.3 Optimality conditions

Proposition 4.2

Let be optimal in (P). Then, the optimality conditions read

Proof We write that is optimal, that is, We take and have whence, passing to the right-hand side, performing some computations, dividing by λ and letting λ to go to 0, we get where generically indicates in each layer. This can still be written, after changing the order of integrations, as

Then, we test Equation (90) by p_i and integrate over (0, T), using the initial data and boundary conditions. We obtain writing globally

where we recall that Combining Equations (106) and (107) we finally deduce which can be still written

In this relation, is the subdifferential at of the indicator set of U and is the normal cone to U at 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 means that 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₀=b_m, but in a way expressed by relations (105). We conclude that the controller 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, 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 .

5. Numerical results

Some graphics revealing the feature of the process with respect to the change of b₀ 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 in all layers (in order to compare how results change for various b₀) with the following dimensionless data:

The figures represent the values of b and c along time, drawn at x fixed (0, 2, 4, 6, 8, 10) corresponding to two values of the initial datum b₀. In Figure 1, these graphics are plotted for b₀, the step function

while in Figure 2 the graphics are displayed for

Figure 1. Solutions b(t, ξ) (left) and c(t, ξ) (right) for b₀(x) given by (109).

Figure 2. Solutions b(t, ξ) (left) and c(t, ξ) (right) for b_0inv(x) given by (110).

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]).

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₀. We prove that the control problem has at least a solution and provide the structure of the controller b₀ in each layer.

Acknowledgement

The author thanks the reviewers for their comments and suggestions related to the control problem.

Funding

This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS –UEFISCDI [project number PN-II-ID-PCE-2011-3-0045].

Notes

†Dedicated to Mimmo Iannelli on the occasion of his 65th birthday.