Efficient Numerical Simulation of Biochemotaxis Phenomena in Fluid Environments

A novel dimension splitting method is proposed for the efficient numerical simulation of a biochemotaxis model, which is a coupled system of chemotaxis–fluid equations and incompressible Navier–Stokes equations. A second-order pressure correction method is employed to decouple the velocity and pressure for the Navier–Stokes equations. Then, the alternating direction implicit scheme is used to solve the velocity equation, and the operator with dimension splitting effect is used instead of the traditional elliptic operator to solve the pressure equation. For the chemotactic equation, the operator splitting method and extrapolation technique are used to solve oxygen and cell density to achieve second-order time accuracy. The proposed dimension splitting method splits the two-dimensional problem into a one-dimensional problem by splitting the spatial derivative, which reduces the computation and storage costs. Finally, through interesting experiments, we show the evolution of the cell plume shape during the descent process. The effect of changing specific parameters on the velocity and plume shape during the descent process is also studied.


Introduction
Biochemotaxis models are used to describe oxidation patterns of some cells in nature.These cells tend to inhabit viscous fluids.The cells and chemical attractants are moved about with the fluid.The gravity produced by the cell aggregation then has an impact on the fluid's flow.Therefore, the biochemotaxis equation is usually coupled to the incompressible Navier-Stokes (NS) equation that governs fluid motion.The aim of this study is to develop an effective second-order dimension splitting scheme for the chemotaxis-Navier-Stokes equations proposed by [1,2].This scheme is designed to achieve the accurate and fast numerical simulation of cell biochemotaxis.The chemotaxis-Navier-Stokes model includes a reaction-diffusion equation that involves changes in oxygen concentration, a convectiondiffusion equation that describes how cells travel along an oxygen concentration gradient, and an incompressible Navier-Stokes equation that governs the motility of the fluid.
We now introduce the dimensionless chemotaxis-Navier-Stokes equation and define the dimensionless parameters.The equation is defined in terms of non-negative parameters, namely α, β, γ, and ζ, as well as the Schmidt number S c .We briefly outline these parameters and their meaning in the equation: where L and q r are characteristic length and characteristic cell density, respectively.c air is the air's oxygen concentration.V b is the volume of the cell, g = 9.8 m/s 2 is the gravitational acceleration, and ρ b is the density of the cell [3].In addition, the diffusion coefficients α n , α c and the chemotactic sensitivity χ q are non-negative.In a smooth, open, bounded, and square two-dimensional domain Ω ⊂ R 2 and a time interval (0, T], the dimensionless chemotaxis-Navier-Stokes system is as follows: ∂c ∂t + u∇c − ζ∆c + βr(c)q = 0, in Ω × (0, T], where the variable q = q(x) denotes the cell density with x = (x 1 , x 2 ) and the oxygen concentration is indicated by the variable c = c(x).The variables u = (v 1 , v 2 ) T , and p = p(x) are defined as the velocity and the pressure of the incompressible NS equation with a density ρ and viscosity η, respectively.Moreover, f(q) = (0, S c γq) T and: where > 0 is a small constant.Let ∂Ω be the boundary of computational domain Ω.The boundary ∂Ω top represents the fluid-air interface.∂Ω top is the boundary at the bottom of computational domain Ω and ∂Ω lr is the boundary at the left and right of computational domain Ω.On the ∂Ω top , there is no cell flux and the air's oxygen concentration c air has saturated.In this paper, the numerical simulation is considered in a closed area with only the top contacting air, so the boundary conditions are set as: In addition, the initial conditions are: As shown above, the chemotaxis-Navier-Stokes system is a non-linear, coupled system with multiple variables.To efficiently solve this system, it is important to decouple the equations.Additionally, a fine-mesh discretization is necessary for achieving high-resolution results in the two-dimensional case.However, it may result in high storage requirements and computational complexity.Thus, achieving high resolution while maintaining stability and computational efficiency remains an important research field.
In response to numerical instabilities and maintaining algorithmic accuracy, researchers have developed many effective numerical methods.In [4], a high-resolution hybrid finitevolume finite-difference method is developed for the chemotaxis-fluid system.The method takes two different approaches to study two different parts of the system (chemotaxis and fluid).In [3], based on the linear finite element method, Huang et al. applied the fluxcorrected transport algorithm to the chemotactic equation to ensure high resolution and stability.In [5], the chemotactic stage uses a high-resolution method (from CLAWPACK) that captures steep gradients, whereas the diffusion step employs the stable TR-BDF2 methods.In [6], to achieve stability, a finite volume approach based on a second-order central upwind structure is used.Recently, Wang et al. established a new convergence-fluid diffusion domain (cf-DD) model, which is derived using a diffusion domain, and developed a second-order hybrid finite-volume finite-difference method to ensure correctness and high resolution [7].The above works mainly consider the stability and high precision of numerical solutions.For two-dimensional problems, standard spatial discretization is performed directly in two-dimensional space, which has high computational complexity and storage requirements.In this paper, the dimension splitting method is used to solve the computational efficiency problem.
In response to the huge storage and computational effort required for two-dimensional discretization, the dimension splitting method is developed.The key is to split the operator according to the direction of the spatial derivative.An early dimension splitting method was the alternate direction implicit (ADI) method, which relied heavily on the finite difference method [8,9].Operator splitting can be used in combination with most common numerical methods, so it is widely chosen as a dimension reduction method.For example, the operator splitting method is used for the convection Cahn-Hilliard equation [10], the Allen-Cahn-type equation [11], and the phase field crystal equation [12].
Our innovation is to establish a fast and low-cost numerical method for the chemotaxis-Navier-Stokes equation by using the dimension splitting method.First, a second-order pressure correction scheme with splitting effect is used to decouple the NS equation.
The key is to use the operator A = (I − ∂ xx ) I − ∂ yy instead of the elliptic operator to correct the pressure [13], so that the decoupled two-dimensional pressure equation can be split into one-dimensional problem for computation.The ADI scheme is then used to solve the decoupled velocity field equation.For the chemotaxis equation, by dividing the second-order spatial partial derivatives, the sequential splitting method is employed to the dimensionality reduction, and then the time accuracy is improved to the second-order by an extrapolation method.The proposed method decomposes the two-dimensional equation into one-dimensional PDE subproblems.The computation of one-dimensional subproblems in the same direction is independent, so that parallel computation can be performed, which greatly reduces the computation and storage effort.For space discretization, the MAC finite difference scheme is considered [14,15].
The rest of the paper is organized as follows.In Section 2, we show the dimension splitting method for the coupled system of chemotaxis-fluid equations and incompressible Navier-Stokes equations.In addition, the marker-and-cell (MAC) finite difference scheme is introduced.In Section 3, we show some numerical experiments to demonstrate the efficacy of the suggested strategy.The concluding Section 4 of the article provides a comprehensive summary of the findings and discussions from the preceding sections, encapsulating the key takeaways and observations.

Numerical Schemes
In this section, we present the dimension splitting scheme of the chemotaxis-Navier-Stokes model and the second-order central difference scheme in the MAC grid configuration.

The Dimension Splitting for the Chemotaxis-Navier-Stokes System
We give a second-order pressure-stabilization perturbation scheme for the NS equation in (6): where A = (I − ∂ xx ) I − ∂ yy and 0 < χ ≤ 1.The operator A replaces the Poisson operator in the majority of the literature, so the two-dimensional pressure equation is decomposed into one-dimensional subproblems to compute.In addition, this method not only has the advantage of high efficiency, but also has the same stability and convergence of Poisson projection (reference [13]).
Let δt be the time step size, t n = nδt, n = 0, 1, . . .Let u n , c n , and q n be the approximations of u, c, and q at t n , respectively.p n− 1 2 is the approximation of p at t n− 1 2 .Setting p − 1 2 = ψ − 1 2 = 0, in each time interval t n−1 , t n .Therefore, the second-order and fast numerical solution of the chemotaxis-Navier-Stokes equation can be obtained by Steps 1-5 as follows.
Step 1: The ADI scheme for updating the velocity field: where , and q * = 3 2 q n−1 − 1 2 q n−2 .The variables ūn and u n fulfill the Dirichlet boundary conditions corresponding to the boundaries of the X-and Y-directions, respectively.
Step 2: The dimension splitting scheme for updating the pressure p n− 1 2 : where 0 < χ ≤ 1 is to ensure second-order accuracy.The variables ψn− 1 2 and ψ n− 1 2 satisfy the homogeneous Neumann boundary conditions corresponding to the boundaries of Xand Y-directions, respectively.
Step 3: The dimension splitting scheme for solving the oxygen concentration c based on the first-order operator splitting method: where the variables cn and c n satisfy the boundary conditions given by (4) corresponding to the X-and Y-directions, respectively.
Step 4: The dimension splitting scheme for solving the cell density q based on the firstorder operator splitting method: where the variables qn and q n satisfy the boundary conditions given by (4) corresponding to the X-and Y-directions, respectively.
Step 5: Since the numerical solutions q n and c n have only first-order accuracy in time, this paper uses an extrapolation method to improve the convergence order to secondorder.Using Step 3 and Step 4, the time step δt and 1 2 δt are, respectively, used in the time interval [t n−1 , t n ] to obtain the first-order solutions at t n as c n,1 , q n,1 and c n,2 , q n,2 .The numerical solution of u at time Then, we get second-order solutions by adopting the extrapolation c n = −c n,1 + 2c n,2 and q n = −q n,1 + 2q n,2 .

The MAC Finite Difference Scheme
When applying the standard finite difference scheme directly in a standard Cartesian grid for the incompressible Navier-Stokes equation, the resulting numerical pressure is unstable.To address this issue, we use the second-order central difference scheme in the MAC grid configuration to provide a fully discrete scheme for the dimension splitting method above.In the MAC scheme [15], the unknown variable (v 1 , v 2 , p) is placed in different positions, as illustrated in Figure 1.Set the computational domain as Ω = Ω x × Ω y = [x a , x b ] × [y a , y b ], and then a uniform Cartesian grid can be obtained: where N x and N y are given integers; x a , x b , and y a , y b are the start and end points of Ω in the x and y directions, respectively.For intermediate nodes, we can define as: To calculate the gradient of the pressure p in the decoupled velocity equation, we define: For the partial derivatives of variables p, q, and c at the cell center, we let U = (p, q, c), U i,j = U x i , y j , and U ȋ, j = U xi , yj .Define ∇ h U ȋ, j = δ x , δ y U ȋ, j, and ∆ h U ȋ, j = δ xx + δ yy U ȋ, j, where: Define the divergence of velocity on cell center as ∇h • u n ȋ, j = δx v 1, ȋ, j + δy v 2, ȋ, j, where: The first and second partial derivatives of v 1 at the position of v 1 are given by: The first and second partial derivatives of v 2 at the position of v 2 are given by: Then, we let yy .The full discretization for the chemotaxis-Navier-Stokes system can be summarized with the following five steps.

Numerical Simulations 4.1. Convergence Test
In this part, we test the method's convergence.The computational domain is set to be Ω = [−1, 1] 2 for 0 ≤ t ≤ T := 1.The exact solution for the convergence test is chosen as: The parameters α, β, γ, ζ, η, S c , and r(c) are set to be one.The value range of χ is 0 < χ ≤ 1.However, through our computational experience, most of the models coupled with incompressible NS equations have good numerical performance when χ = 1 2 .Therefore, this article fixes the χ as 1  2 based on computational experience.For simplicity, Dirichlet boundary conditions are used for velocity variables and homogeneous Newman boundary conditions for cell and oxygen variables.The spatial discretization size h x and h y is determined by taking {50, 75, 100, 125, 150, 175, 200} divisions in each coordinate direction.To make it easier to verify second-order convergence in space-time, we set the time step size as δt = 0.1 min{h x , h y }.The error is measured through discrete l 2 -norm, i.e., Figure 2a displays the l 2 errors for all unknown variables between the numerical solution and the precise solution at t = 1.It is clear that all unknown variables' error trends are parallel to the second-order reference line, proving that the suggested technique is second-order convergent.

Test with Random Initial Density of Cells
The initial density of cells in this paragraph is adjusted to q 0 (x, y) = 0.8 + 0.2 • rand, with rand being a random number uniformly distributed between 0∼1.Other initial values are consistent with (37).The simulation parameters are set as ζ = 5, α = 10, Sc = 500, η = 1, β = 50, and γ = 5000.In Figure 8, we can clearly see that our numerical approach accurately captures the descent dynamics of the plume.The convective instability occurs at t = 0.095.Apparently, the plumes are unevenly spaced and sink at different rates.At t = 0.141, two rapidly sinking plumes form mushroom shapes.At t = 0.144, all the plumes form mushroom shapes.

Conclusions
A fast dimension splitting scheme for chemotactic-Navier-Stokes systems is proposed.This method combines the pressure correction method with the dimension splitting effect, the ADI method, the operator splitting method, and the extrapolation method to achieve a fast second-order solution.The method decomposes the two-dimensional model into a series of one-dimensional subproblems, thus reducing computational complexity and storage space.The effectiveness and practicability of this method are demonstrated by numerical experiments simulating chemotaxis.Future research will focus on the following directions: (1) the dimension splitting schemes for three-dimensional chemotactic-Navier-Stokes systems: the method proposed in this article can be appropriately modified for solving the coupling system of the incompressible NS equation and the second-order parabolic equation.For example, the effective numerical results of the Allen-Cahn model of two-phase incompressible fluid, the surfactant model of binary fluid, etc., are what we will present in the future.(2) Due to the complexity and nonlinearity of the model, its stability analysis is a challenging problem, which will continue to be solved in the future.(3) The rapid numerical simulation of biochemotaxis in complex regions has more important practical application value.How to establish the dimensional splitting scheme of complex regions is the focus of our follow-up research.

Figure 1 .
Figure 1.A diagram of the location of variables (v 1 , v 2 , p, q, c).

Figure 2 .
Figure 2. The convergence and accuracy comparison results of the proposed method.(a) The proposed dimension splitting method.(b) Standard FD discretization without dimension splitting.