Chemotaxis Model for Drug Delivery Using Turing’s Instability and Non-Linear Diffusion

Featured Application: this work has a potential application in several domains especially in the ﬁeld of drug delivery control and optimization for nanomedicine. Abstract: This paper is devoted to the study of the chemotaxis model for drug delivery purposes. The pattern formation for a volume-ﬁlling with nonlinear diffusive terms is investigated. The proposed mathematical model is governed by a reaction–diffusion system modeling the interaction between the cell density and the concentration of the chemoattractant. We investigate the pattern formation for the model using Turing’s principle and linear stability analysis. An asymptotic expansion is used to linearize the nonlinear diffusive terms. Next, we introduce an implicit ﬁnite volume scheme; it is presented on a triangular mesh satisfying the orthogonality condition. Finally, numerical results showing the formation of the spatial pattern for the chemotaxis model are presented and analyzed. The results demonstrate promising progress in understanding the process of controlling and designing targeted drug delivery.


Introduction
In many societies, cancer is the leading cause of death and a serious health burden. Different anticancer and chemotherapy medicines have been created to combat these risky diseases [1]. The majority of chemotherapy drugs do not have a selective effect on cancer cells and can destroy healthy cells as well. As a result, supplying the medication directly to the tumor without affecting healthy tissues is a significant obstacle in cancer care [2,3].
Various attempts, such as "smart drug delivery" using nano/submicron droplets [4], have been made to resolve this problem. Researchers looking into smart drug distribution are having trouble figuring out how to get the drug to the right place at the right time.
Drug delivery is one of the most focused research areas for treating such complex diseases as cancerous tumors in various body organs such as the brain, breast, lungs, kidney, and heart. Nanoparticles such as iron oxide are used in nanomedicine to diagnose and treat diseases, and MRI imaging with iron oxide was recently used to track drug delivery in mice. Bacteria were used to transport drugs to tumor sites [5], with chemoattractant gradients controlling the transport path.
Chemotaxis is a basic cell and organism guidance system that attracts microbes to sustenance, immune cells to infection sites. The Patlak-Keller-Segel (PKS) method is a good choice for modelers and analysts alike, and it is a portion of the capstone of biological modeling. Chemotaxis is the tendency of cells and living organisms to move in response to Appl. Sci. 2021, 11, 4979 2 of 10 an external stimulus in their environment. Due to its critical function in cancer metastasis, embryogenesis, cell growth, inflammation, and wound healing [6], chemotaxis has been extensively studied in the last two decades [7,8].
Macrophages are known to migrate to areas of inflammation, which are typically where hypoxic tumors are found. In [9], a mathematical model of macrophage migration toward hypoxic tumor sites is suggested as a drug delivery vehicle. To demonstrate macrophage penetration into tumors, the authors propose a rough modeling process. They also look into the action and development of chemotaxis and chemokines, as well as the efficacy of drug-loaded macrophages for drug delivery. A growing avascular tumor spheroid serves as the basis for the model.
Chemotaxis is a subset of the reaction-diffusion systems method. Through selforganizing/pattern formation, reaction-diffusion systems create a wide range of spatiotemporal patterns, including propagating circular, spiral, and periodic waves. For the first time, such spatiotemporal patterns are observed in the Belousov-Zhabotinsky reaction system [10]. A standard model for explaining the pattern-forming mechanism is the reaction-diffusion system. Indeed, Alan Turing developed in his seminal paper [11] the idea that chemicals can react and diffuse in such a way that they create steady-state nonhomogeneous spatial patterns of chemical concentrations under certain conditions. Although a diffusion process generally brings uniform distribution of a substance, the scenario of Turing presents a non-uniform pattern as a stationary state in a reaction-diffusion system.
We can find a wide variety of reaction-diffusion systems in the fields of physics, chemistry, and biology [12]. One of the most popular reaction-diffusion systems is the chemotaxis model.
Chemotaxis is the feature movement of a cell along a chemical concentration gradient either toward the chemical stimulus, in which case the chemical is called chemoattractant, or away from the chemical stimulus, in which case the chemical is called chemo-repellent.
Theoretical and mathematical modeling of chemotaxis dates to the pioneering work of Patlak in the 1950s [13] and Keller and Segel in the 1970s [14,15].
The review article by Horstmann [16] provides a detailed introduction into the mathematics of the Keller-Segel (KS) model for chemotaxis. The general form of the Keller-Segel model in which we are interested is given by the system In the above model, u denotes the cell density on a given bounded domain ⊂ R 2 for a time t ∈ (0, T ), T > 0, and v describes the concentration of the chemoattractant. D(u) denotes the diffusivity of the cells (sometimes also called motility), while χ (u, v) stands for the chemo-sensitivity; however, the coefficient D v represents the diffusivity of the chemicals. The function f (u, v) describes cell growth and cell death, while g (u, v) describes production and degradation of the chemoattractant. A key property of the above equations is their ability to give rise to spatial pattern formation when the chemical signal v acts as an auto-attractant, that is, when cells both produce and migrate up gradients of the chemical signal. The Keller-Segel model (1) is studied by many authors, e.g., see the review articles of Horstmann [16,17] and the paper [18]. From a numerical point of view, we motion that authors in [19] analyzed a finite volume method for the Keller-Segel model (1).
In this paper, we are interested in a squeezing probability modeling the elasticity feature of the cell, we study the pattern formation for the corresponding volume-filling chemotaxis model with nonlinear diffusive terms, and finally, we investigate numerically the generation of spatial patterns in a 2-D bounded domain.

Volume-Filling Chemotaxis Model
In this section, we consider the volume-filling chemotaxis model (1) for which we introduce the squeezing probability q of a cell finding space at its neighboring location. q is a non-increasing function concerning the cell density u, and also, it is nonlinear greater than the linear one that corresponds to the interpretation that the cells are solid blocks [8,20,21] and equals to zero when the maximum number of cells exceeds a certain number u.
The number u represents the total number of cells that can be accommodated at any site, also called the crowding capacity. We take the squeezing probability q reflecting the elasticity characteristic of the cells according to positive parameter γ, introduced in [22]: In the following, the diffusion term and the cross-diffusion are given by where q is the squeezing probability for the cell finding space defined in (2), and d 1 is the diffusion coefficient. In the sequel, for the simplicity, we take d 2 = D v . Substituting these coefficients into system (1), we obtain with zero flux boundary conditions on Σ T = ∂ × (0, T) where n is the unit normal vector at the boundary ∂ and outward to . We assume that the initial conditions are non-negative, for all x ∈ Under realistic assumptions on the functions f , g, and q, the existence of global bounded classical solutions of system (3) object to the boundary conditions (4) that has been established in [22] by the maximum principle and Amman's theory [23].

Pattern Formation
Pattern formation in mathematics is the process that, by changing a bifurcation parameter, the spatially homogeneous steady states lose stability due to spatially inhomogeneous perturbations, and stable inhomogeneous solutions arise.
In this section, we are interested in the pattern formation for the system (3) and (4). We derive here the necessary conditions for diffusion-driven instability (also called Turing instability) of the steady-state and the initiation of spatial pattern for the above-mentioned system. We look for the spatially homogeneous steady states by setting the kinetics on the right-hand side of Equation (3) to be zero: We use the same technique as Murray [12,24]. We suppose that (u s , v s ) is a positive solution of Equation (5); where (u s , v s ) is a steady state of the system where A is the stability matrix (i.e., the Jacobian matrix of system (6) at the steady state (u s , v s ).
From now on, we take the derivatives of f and g to be evaluated at the steady state (u s , v s ) unless stated otherwise. Therefore, the conditions for which the steady state (u s , v s ) is linearly stable; that is, Re λ < 0, where λ is the eigenvalue of the matrix A, is guaranteed if Now, we consider the full chemotaxis model (3). We examine a small perturbation from the homogeneous steady state (u s , v s ) of the form where ε 1 is a small parameter strictly positive. f is a smooth function, then using an asymptotic expansion we obtain The system (3) transforms into the linearized system in (9) after performing the following steps: Asymptotic expansion for the functions g, q, and χ around (u s , v s ); • Taking into account Equation (5); • Equating first-order terms with respect to ε; • Neglecting higher-order terms, and dropping the tilde for the convenience.
For simplicity, we denote by ζ s = q(u s ) -q (u s )u s , ψ s = u s q(u s ) χ(v s ) and by χ s = χ(v s ) with the given functions in (14), one can deduce that ζ s > 0 and ψ s < 0. In what follows, we perform the standard argument (e.g., [24], Chapter II). We obtain the dispersion relation associated with system (9) where a = a k 2 and b = b k 2 are two functions given by For the steady state to be unstable to spatial perturbations, we require Re(λ) > 0 for some k = 0. We remark that the sum of the eigenvalues is strictly negative since the function a k 2 is strictly positive due to the conditions (7). Thus, the instability can only happen if b k 2 becomes negative for some k = 0, i.e., the dispersion relation (10) has two roots with opposite signs. Carrying the function b k 2 , then the last condition requires that equation (11) be negative.
Since we required that the determinant |A| = f u g v − f v g u > 0 from the conditions (7), the only possibility for inequality (11) to hold, is that the sum of the roots is positive. That leads to Equation (12) Referring to the condition (11), and since ζ s d 1 d 2 is positive, it is necessary that the discriminant of b k 2 = 0 be positive. The discriminant of b k 2 is represented by the following equation: Applying inequality (12), we obtain from this inequality that In this paper we take specific functions in system (3), where f represents the cell density kinetics, the logistic growth, while g represents the chemo-attractant kinetics. A linear function depending on each of the variables u and v represents the chemoattractant growth with rate ν and the decay with rate δ due to dilution. They are given by the form Applying the same pattern formation analysis as before and for specific functions (14), we obtain the only necessary condition in Equation (15) for the generation of spatial patterns: where ϕ(u) = uq(u). In the sequel of this section, we investigate the range of unstable wave numbers for which, in addition to condition (15), we can expect pattern formation for system (3). To do this, we study the bifurcation with the chemo-sensitivity χ, the growth rate ν, and the decay rate δ of the chemoattractant. For the bifurcation parameter χ, we compute the bifurcation value χ c , and we obtain the following: Furthermore, the range of unstable wave numbers k 2 1 < k 2 < k 2 2 is obtained from the zeros of the equation b k 2 = 0, and k 2 1 and k 2 2 are defined as follows: where S = −µd 2 − δD(u c ) + χvϕ(u c ). Figure 1 shows how the eigenvalue l varies as a function of k 2 for various χ. By analyzing this figure, we can deduce that the pattern formation can be expected when the values of χ exceed the critical values χ c , since there exists a range of wave numbers k 2 1 < k 2 < k 2 2 when the dispersion relation (10) (10). When > , there is a range of unstable wave numbers < < for which the homogeneous steady state is linearly unstable.

Numerical Scheme and Simulations
In this section, we use a finite volume scheme to discretize system (3). Following [19], assuming that Ω is an open bounded polygonal domain for which we consider an admissible mesh consisting of a triplet ( , ℰ , ), where is a family of disjoint open polygonal convex subsets of Ω called control volumes, ℰ is a family of subsets of Ω containing the edges of the control volumes, and is a family of points of Ω satisfying the orthogonality condition; that is, the line segment joining the centers of two neighboring control volumes is orthogonal to the interface between these control volumes (see, e.g., [25]). We give some notations on the mesh: The set of neighbors of is denoted by ( ), which is ( ) = { ∈ | ∃ ∈ ℰ , σ̅ = ∩ }, where ℰ is the subset of ℰ containing the edges of . We denote by ɳ , the unit normal vector to , = ∩ outward to . Integrating the equations of system (3) over each control volume , and using the Green-Gauss formula, one can deduce the following discretization for system (3) for all ∈ : With the discrete unknowns = ( ) ∈ and = ( ) ∈ , ∈ ⟦0, … , ⟧.   (10). When χ > χ c , there is a range of unstable wave numbers k 2 1 < k 2 < k 2 2 for which the homogeneous steady state is linearly unstable.
Performing the same analysis as for the bifurcation parameter χ, we obtain evidence that that the pattern formation can also be supported by increasing the value of v above the critical value ν c or by decreasing the value of δ below the critical value δ c .

Numerical Scheme and Simulations
In this section, we use a finite volume scheme to discretize system (3). Following [19], assuming that is an open bounded polygonal domain for which we consider an admissible mesh consisting of a triplet (T , E , P ), where T is a family of disjoint open polygonal convex subsets of called control volumes, E is a family of subsets of Ω containing the edges of the control volumes, and P is a family of points of satisfying the orthogonality condition; that is, the line segment joining the centers of two neighboring control volumes is orthogonal to the interface between these control volumes (see, e.g., [25]). We give some notations on the mesh: The set of neighbors of K is denoted by N (K), which is N (K) = {L ∈ T ∃σ ∈ E K , σ = K ∩ L , where E K is the subset of E containing the edges of K. We denote by η K,L the unit normal vector to σ K,L = ∂K ∩ ∂L outward to K.
Integrating the equations of system (3) over each control volume K, and using the Green-Gauss formula, one can deduce the following discretization for system (3) for all K ∈ T : for all K ∈ T and n ∈ [0, . . . , N]: where c + = max(c, 0), c − = max(−c, 0), and Φ = χφ(u). Φ ↓ represents the nonincreasing part of Φ, and Φ ↑ represents the nondecreasing part of Φ. We refer to [19] for details about the convergence analysis of the discretized problem (17) and (18).
In the next section, we show the numerical simulations of the model (3) in twodimensional space to illustrate our experimental results and to demonstrate the patterns generated by the system. The Newton algorithm is used to approach the solution of the nonlinear system defined by the first Equation of (18). This algorithm is coupled with a bigradient method to solve linear systems arising from the Newton algorithm process. The simulation of the model is applied in two-dimensional domain = (0, 10) × (0, 10) for which a non-uniform admissible grid is considered. The following data are considered for the numerical test: Further, our experimental simulation considers a small time step ∆t = 0.05, and the numerical flux G is defined by Equation (19); we consider Φ ↑ = Φ(min{z, u m }) and . The aforementioned parameters yield to the homogeneous steady state, (u s , v s ) = (0.25, 0.25), and to the range of unstable wave numbers 0.11 < k 2 < 426.36.
Initially the cell density is located with as a small random perturbation around the homogeneous steady state; in our test, we take u (0, x, y) = u s + rand(0, 1). Meanwhile, for the chemoattractant, we take v(0, x, y) = v s . Figure 2 shows the plot of the initial condition for the cell density distribution; the red regions correspond to the highest cell density so the cells are aggregated into a "random" distribution. A zero-flux boundary condition is imposed for the boundary conditions. where = max(c, 0), = max(− , 0) , and Φ = ( ). Φ ↓ represents the nonincreasing part of Φ, and Φ ↑ represents the nondecreasing part of Φ. We refer to [19] for details about the convergence analysis of the discretized problem (17) and (18).
In the next section, we show the numerical simulations of the model (3) in two-dimensional space to illustrate our experimental results and to demonstrate the patterns generated by the system. The Newton algorithm is used to approach the solution of the nonlinear system defined by the first Equation of (18). This algorithm is coupled with a bigradient method to solve linear systems arising from the Newton algorithm process. The simulation of the model is applied in two-dimensional domain Ω = (0, 10) × (0, 10) for which a non-uniform admissible grid is considered. The following data are considered for the numerical test:  Initially the cell density is located with as a small random perturbation around the homogeneous steady state; in our test, we take (0, , ) = + rand(0, 1). Meanwhile, for the chemoattractant, we take (0, , ) = . Figure 2 shows the plot of the initial condition for the cell density distribution; the red regions correspond to the highest cell density so the cells are aggregated into a "random" distribution. A zero-flux boundary condition is imposed for the boundary conditions. Figures 3 and 4 show the evolution of the cell density in space for different time moments. It shows that the random distribution of the cell density leads to a merging process in all directions of the space at = 0.05 s. Then after = 0.85 s, the spatial patterns begin to form, and the aggregations start to merge and emerge and then to form uniform spot patterns as at time = 14 s. Next, after time = 30 s, these spot patterns stabilize (with a stationary -norm of solutions less than 10 ) to finally form 33 spatially inhomogeneous stable patterns as shown in the snapshot to the bottom of Figure 4.

Discussion
The mathematical analysis of chemotaxis models shows a plenitude of spatial patterns such as the chemotaxis models applied to skin pigmentation patterns [7,26,27]. That leads to aggregations of one type of pigment cell into a banded spatial pattern and to other models applied to the aggregation patterns in an epidemic disease [28], tumor growth [29], angiogenesis in tumor progression [30], and mathematical modeling of macrophages as vehicles for drug delivery to the tumor, as was used in [9]. Recently, direct transport of bacteria-based drug delivery vehicles was developed in [31], and iron oxide (IO) nanoparticles are used to track the chemotaxis of macrophages for the drug delivery control process using MRI imaging and visualization [5]. Here we investigated numerically the Turing instability and pattern formation in two dimensions for the chemotaxis model to show that the cell density stabilizes at a certain moment, thus verifying the theoretical results about pattern formation using an appropriate numerical scheme. Indeed, Figure 5 shows the time evolution of the cell density at the point (3.275; 3.1) in the red line, (8; 7.5) in the blue line, and (4.975; 2.825) in the cyan line. We observe that the cell density at the three points varies in response to the gradient of the chemoattractant, which plays an essential role to stop the aggregation of the cells. Next, the cell density stabilizes at the three points when the time t is greater or equal to 23 s, and we obtain the same result for the

Discussion
The mathematical analysis of chemotaxis models shows a plenitude of spatial patterns such as the chemotaxis models applied to skin pigmentation patterns [7,26,27]. That leads to aggregations of one type of pigment cell into a banded spatial pattern and to other models applied to the aggregation patterns in an epidemic disease [28], tumor growth [29], angiogenesis in tumor progression [30], and mathematical modeling of macrophages as vehicles for drug delivery to the tumor, as was used in [9]. Recently, direct transport of bacteria-based drug delivery vehicles was developed in [31], and iron oxide (IO) nanoparticles are used to track the chemotaxis of macrophages for the drug delivery control process using MRI imaging and visualization [5]. Here we investigated numerically the Turing instability and pattern formation in two dimensions for the chemotaxis model to show that the cell density stabilizes at a certain moment, thus verifying the theoretical results about pattern formation using an appropriate numerical scheme. Indeed, Figure 5 shows the time evolution of the cell density at the point (3.275; 3.1) in the red line, (8; 7.5) in the blue line, and (4.975; 2.825) in the cyan line. We observe that the cell density at the three points varies in response to the gradient of the chemoattractant, which plays an essential role to stop the aggregation of the cells. Next, the cell density stabilizes at the three points when the time t is greater or equal to 23 s, and we obtain the same result for the

Discussion
The mathematical analysis of chemotaxis models shows a plenitude of spatial patterns such as the chemotaxis models applied to skin pigmentation patterns [7,26,27]. That leads to aggregations of one type of pigment cell into a banded spatial pattern and to other models applied to the aggregation patterns in an epidemic disease [28], tumor growth [29], angiogenesis in tumor progression [30], and mathematical modeling of macrophages as vehicles for drug delivery to the tumor, as was used in [9]. Recently, direct transport of bacteria-based drug delivery vehicles was developed in [31], and iron oxide (IO) nanoparticles are used to track the chemotaxis of macrophages for the drug delivery control process using MRI imaging and visualization [5]. Here we investigated numerically the Turing instability and pattern formation in two dimensions for the chemotaxis model to show that the cell density stabilizes at a certain moment, thus verifying the theoretical results about pattern formation using an appropriate numerical scheme. Indeed, Figure 5 shows the time evolution of the cell density at the point (3.275; 3.1) in the red line, (8; 7.5) in the blue line, and (4.975; 2.825) in the cyan line. We observe that the cell density at the three points varies in response to the gradient of the chemoattractant, which plays an essential role to stop the aggregation of the cells. Next, the cell density stabilizes at the three points when the time t is greater or equal to 23 s, and we obtain the same result for the other spot patterns. Consequently, this stabilizing proves that the volume-filling chemotaxis model (3) generates stationary spatial patterns. These results demonstrate that our proposed chemotaxis model has a great potential to improve drug delivery processes.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 10 other spot patterns. Consequently, this stabilizing proves that the volume-filling chemotaxis model (3) generates stationary spatial patterns. These results demonstrate that our proposed chemotaxis model has a great potential to improve drug delivery processes.

Conclusions
In this article, we propose a mathematical model and a finite volume numerical method to simulate the pattern formation for the volume-filling chemotaxis model. We carry out conditions of pattern formation for the general volume-filling chemotaxis model with nonlinear diffusive terms. In the diffusive term, we apply a particular squeezing probability to describe the elasticity feature of the cells and to take into account that the cells are not solid blocks. In addition, we consider logistical growth to describe the cell kinetics. A two-dimensional numerical simulation is presented for nonzero kinetics to prove the generation of spatial inhomogeneous stable patterns for the volume-filling chemotaxis model. We have identified the mechanism for destabilizing the stationary solutions. We have found the conditions linking the biological parameters of the system in order to locate a tumor in space; consequently, the targeted drug delivery can be applied rigorously. These results need to be compared and completed by some experiments, which will be the subject of a future work.

Conclusions
In this article, we propose a mathematical model and a finite volume numerical method to simulate the pattern formation for the volume-filling chemotaxis model. We carry out conditions of pattern formation for the general volume-filling chemotaxis model with nonlinear diffusive terms. In the diffusive term, we apply a particular squeezing probability to describe the elasticity feature of the cells and to take into account that the cells are not solid blocks. In addition, we consider logistical growth to describe the cell kinetics. A two-dimensional numerical simulation is presented for nonzero kinetics to prove the generation of spatial inhomogeneous stable patterns for the volume-filling chemotaxis model. We have identified the mechanism for destabilizing the stationary solutions. We have found the conditions linking the biological parameters of the system in order to locate a tumor in space; consequently, the targeted drug delivery can be applied rigorously. These results need to be compared and completed by some experiments, which will be the subject of a future work.