Convergence Analysis of a Stochastic Interacting Particle-Field Algorithm for 3D Parabolic-Parabolic Keller-Segel SystemsThanks: To appear in Mathematics of Computation.
Abstract
Chemotaxis models describe the movement of organisms in response to chemical gradients. In this paper, we propose a stochastic interacting particle-field algorithm integrated with random batch approximation (SIPF-) for the three-dimensional (3D) parabolic-parabolic Keller-Segel (KS) system, also referred to as the fully parabolic KS system. The SIPF- method approximates the KS system by coupling particle-based density representations with a smooth field variable computed via spectral methods. By incorporating the random batch method (RBM), we bypass the mean-field limit and significantly reduce computational complexity. Under mild assumptions on the regularity of the original KS system and the boundedness of numerical approximations, we prove that the empirical measure associated with the SIPF- particle system converges in the -Wasserstein distance to the exact measure of the limiting McKean-Vlasov process with high probability. Finally, we present numerical experiments to validate the theoretical convergence results, and demonstrate the efficiency and robustness of the SIPF- method as a diagnostic tool for density concentration and potential finite-time singularity in 3D parabolic-parabolic KS systems. Specifically, the SIPF- method identifies the critical threshold of initial mass leading to blow-up under specific initial data configurations.
keywords
Fully parabolic Keller-Segel system, stochastic interacting particle-field (SIPF) algorithm, random batch method, 3D simulations, convergence analysis.MSC
35K51, 65C05, 65M12, 65M75, 65T50.1 Introduction
Chemotaxis is a biological phenomenon involving the movement of organisms (e.g., bacteria) in response to signals (e.g., chemoattractants), which can be produced by the organisms themselves. Theoretical and mathematical modeling of this phenomenon was initiated by Patlak [40] and by Keller and Segel [26]. In this work, we focus on the fully parabolic KS system as follows:
| (1) |
where are positive constants, and are non-negative constants. The model is called parabolic-elliptic if , and fully parabolic if . Here, denotes the density of active particles (bacteria), and represents the concentration of chemoattractants emitted by these bacteria. The model finds broad applications in biology, ecology, and medicine [41, 39, 2, 46].
The nonlinear, potentially singular behavior of KS equations, especially blow-up phenomena [35, 17, 1, 12], makes numerical methods essential. Mesh-based methods such as finite difference [5, 42, 11], finite element [43, 9, 44], and finite volume schemes [13, 6, 54] are widely employed. Recent advances include local discontinuous Galerkin methods with optimal convergence rates [30] and semi-discrete symmetrization-based schemes that avoid nonlinear solvers [31]. Despite their success, challenges remain in ensuring stability, convergence, and the effective handling of singularities, making numerical studies on the KS model an active and evolving field of research.
In addition to mesh-based methods, particle-based approaches have also been developed to address the challenges posed by the KS system, providing a complementary perspective. Stevens [48] developed an -particle system and established its convergence for the fully parabolic case. Haškovec and Schmeiser [16] proposed a convergent regularized particle system for the 2D parabolic-elliptic KS model. Moreover, Godinho and Quininao [14] showed well-posedness results and the propagation of chaos property in a subcritical KS equation. Craig and Bertozzi [7] proved the convergence of a blob method for the related aggregation equation. Liu and Yang [32] introduced a random particle blob method with a mollified kernel for the parabolic-elliptic case, proving its convergence when the macroscopic mean-field equation possesses a global weak solution [33]. Due to the potential finite-time blow-up of the KS system, deriving a priori error estimates inherently requires sufficient regularity of the exact solution. Consequently, theoretical convergence analyses in the literature are typically restricted to the pre-blow-up regime [48, 30, 32, 54]. For the singular regime, numerical studies primarily focus on structural preservation (e.g., positivity and energy dissipation) and empirical robustness [6, 10, 45].
In [51], we proposed a novel stochastic interacting particle-field (SIPF) algorithm for the fully parabolic KS system (1) in 3D. The SIPF method approximates the KS solution as the empirical measure of particles (see Eq. (2)) coupled with a smoother field variable computed using the spectral method (see Eq. (3)). The algorithm employs an implicit Euler discretization and a one-step recursion based on Green’s function, which efficiently captures concentration behavior and potentially finite-time blow-up with only dozens of Fourier modes in each dimension.
Despite the algorithm’s numerical efficiency and stability as observed in [51], a complete convergence analysis remains open. This paper fills that gap by establishing and validating convergence estimates for the SIPF- method, a random-batch variant designed to reduce computational cost. Our main result, presented in Theorem 4, establishes the convergence of the solution obtained from the SIPF- method to the exact solution in the pre-blow-up regime under mild assumptions. Specifically, the -Wasserstein distance between the SIPF- and exact density distributions, denoted as , depends on the time step , the number of Fourier modes , and the number of particles , and scales as . The proof, detailed in Section 3, builds upon key lemmas that quantify single-step update errors and analyze their propagation over time. Furthermore, we demonstrate that the spectral truncation inherently regularizes the singular kernel, ensuring the unconditional stability of the algorithm even as the exact solution approaches the finite-time singularity.
The error between the chemical concentration in the SIPF- method and the exact solution originates from Fourier truncation and implicit Euler time discretization. Under appropriate regularity assumptions on , the truncation error decays as the number of Fourier modes increases. The temporal error can be expressed through differences between Fourier coefficients of and , which further couples with the particle trajectory error accumulated from the preceding step. Here, and are the exact dynamics of the system associated with the density and its numerical approximations obtained via our methods, respectively.
At the numerical discretization level, the RBM [24, 25, 23, 4] is incorporated into the SIPF- method, which renders the particles fully independent and identically distributed (i.i.d.), thereby effectively circumventing the need to address the propagation of chaos [33]. At each time step, small random batches of particles are selected with replacement for particle interactions. In the error estimate between and , the i.i.d. property, together with the complex-valued mean value theorem [34] and Bernstein’s inequality, bound the deviations of the empirical mean of particle trajectories from its expectation. This deviation accounts for the uncertainty described in Theorem 4. The numerical experiments in Section 4 further demonstrate that, with the introduction of the RBM, the numerical examples maintain a high level of precision.
The error between and is influenced by the gradient of and , reflecting the dependence of the particle trajectories on the interaction potential. Using Parseval’s identity [27], we relate to Fourier-coefficient errors, establishing two coupled recursive inequalities between and . The first inequality relates to ; the second updates by incorporating a term (see Eqs. (70)-(71)). Substituting and decoupling these aforementioned recursive inequalities yields a bound for that depends only on earlier errors. Through the natural coupling , this error estimate translates to the -Wasserstein distance . Applying the discrete Gronwall inequality together with Fourier-coefficient error estimates leads to the global convergence result in Theorem 4.
In addition to the theoretical analysis, Section 4 provides comprehensive numerical experiments to validate the performance of the SIPF- method. These experiments verify the predicted convergence rates, confirm the validity of the regularity assumptions, and explore the complex dynamics of the KS system. A particular focus is placed on identifying the initial-data-dependent threshold for finite-time blow-up in three dimensions under specific initial data configurations. In Subsection 4.2, the SIPF- algorithm identifies blow-up regimes and captures severe concentration toward finite-time singularities (see Figures 4 and 6). Notably, the method remains robust in detecting blow-up even under practical discretization parameters (, , ) that are much coarser than those required for the convergence analysis (see Figure 5).
The rest of the paper is organized as follows. In Section 2, we review the SIPF method for solving the fully parabolic KS system and present the derivation of the SIPF- method, which incorporates the RBM to compute the particle interaction. In Section 3, under certain assumptions, we provide a detailed convergence analysis of the SIPF- method, which decomposes the proof into several lemmas. In Section 4, we present numerical results to validate the necessity of the assumptions, demonstrate the accuracy of the SIPF- method, and confirm the theoretical convergence rate derived in our analysis. Section 5 discusses extensions of the SIPF- framework to a wider class of mathematical-biological chemotaxis models, and Section 6 concludes the paper.
2 Derivation of the SIPF- Method
In this section, we present the SIPF- method for solving the fully parabolic KS model. Since the evolutionary phenomena of interest occur in the interior of the domain, we restrict the system (1) to a large domain and assume Dirichlet boundary conditions for particle density and Neumann boundary conditions for chemical concentration .
Throughout this section, we use the standard notation , , etc., to represent the exact solutions of the fully parabolic KS model. For the variables computed or approximated using the SIPF- algorithm, we instead use the notations , , etc.
As a numerical algorithm, we assume that the temporal domain is partitioned by with and . We approximate the density at by empirical particles , i.e.,
| (2) |
where denotes the conserved total mass (integral of over the domain ). For the chemical concentration , we adopt a Fourier basis approximation. Its unnormalized Fourier coefficients are defined by
The corresponding inverse Fourier representation is
| (3) |
where denotes the index set, i.e.,
| (4) |
and . Defining in the same way, the exact solution can also be approximated by a truncated spatial Fourier series expansion as follows:
| (5) |
Remark 1.
The choice of the Fourier basis over Hermite polynomials for approximating the chemical concentration is based on the fact that, since the blow-up phenomenon is localized in the domain interior, periodic boundary conditions effectively emulate an infinite spatial domain under this configuration. When the spatial location of the singularity remains distant from domain boundaries, its interaction with these artificial edges becomes negligible.
Then at , we generate empirical samples according to the initial condition of and set up using the Fourier series of . For ease of presenting our algorithm, with a slight abuse of notation, we use , and
| (6) |
to represent density and chemical concentration at time .
Considering the time-stepping for the system (1) from to , with and known, our algorithm, inspired by the operator splitting technique, consists of two sub-steps: updating chemical concentration and updating organism density .
Updating chemical concentration
Let be the time step. We discretize the equation of (1) in time by an implicit Euler scheme
| (7) |
From Eq. (7), we obtain the explicit formula for as follows:
| (8) |
It follows that
| (9) |
where is the Green’s function of the operator and represents an approximation of spatial convolution that differs from the continuous setup, as is computed using truncated Fourier basis functions and is given by a discrete particle representation. Unless otherwise stated, all subsequent norms will refer to the norms. In the case of , the Green’s function reads as follows:
| (10) |
Green’s function admits a closed-form Fourier transform,
| (11) |
For the term in Eq. (9), by Eq. (11) it is equivalent to modify Fourier coefficients to .
For the second term , we first approximate using a cosine series expansion. Then, we use the particle representation of given in Eq. (2) to derive
| (12) |
where the factor arises from the redefinition of the Fourier computational domain from to .
Finally, we summarize the one-step update of the Fourier coefficients of the chemical concentration in Algorithm 1, which follows the same procedure as in the original SIPF method [51].
Updating density of active particles
In the one-step update of density represented by particles , we apply the Euler-Maruyama scheme to solve the stochastic differential equation (SDE)
| (13) |
where the variables are i.i.d. standard normal random variables corresponding to the Brownian paths in the SDE formulation. For , substituting Eq. (9) in Eq. (13) gives:
| (14) |
from which is constructed via Eq. (2).
In this particle formulation, the computation of the spatial convolution differs slightly from that in the update of (i.e., Eq. (9)).
For the term , to avoid the singular points of , we evaluate the integral with quadrature points that are away from . Precisely, denote the standard quadrature point in as
| (15) |
where , , are integers ranging from to . When computing , we evaluate at , where a small spatial shift is defined as and at correspondingly. The latter is computed by inverse Fourier transform of the shifted coefficients, with modified to , where denotes the -th component of .
Motivated by mini-batch sampling [15, 47, 49, 50, 52, 53] and the random batch method (RBM) [24, 25, 4, 23], for each particle , we choose a small batch of size randomly with replacement. We restrict the interaction of to particles within this batch; specifically, we approximate the interaction term by the random batch sum .
We summarize the one-step update (for ) of the density in the SIPF- method as in Algorithm 2.
The spatial shift acts as a mild perturbation to avoid the singular attractive force induced by the Green’s function . While it may technically introduce bias, due to the even structure of and the regularity of , the influence of this shift is small when evaluating . Consequently, particle aggregation remains driven by the intrinsic nonlinear dynamics rather than the grid structure. Combined with the stochasticity of the random batch method, this ensures that no systematic directional bias is introduced, which allows for robust singularity detection independent of specific grid points.
Combining Eq. (9) and Eq. (14), we conclude that the recursion from
to
is thus fully defined. We summarize the SIPF- method in the following Algorithm 3.
Computational Complexity
We briefly analyze the complexity of the proposed SIPF- method. The memory usage is , where denotes the total number of Fourier modes. Regarding the computational cost, the original SIPF method typically requires operations for pairwise particle interactions. In contrast, by incorporating the RBM, the interaction cost in the SIPF- method is reduced to , where is the batch size.
Particle-wise Independence due to RBM
In the above derivation, are i.i.d. samples with distribution and independent of . The one-step trajectories follow the discrete-time rule:
| (16) |
where is computed via Eq. (9), and denotes the Brownian motion. It is worth noting that, for the updated position of the -th particle by Eq. (13), the interaction term, is computed by , where the selection of is independent of and hence can be viewed as i.i.d. samples of independent of and . Together with the independent Brownian motion term, we can deduce the independence of .
Correspondingly, we denote the exact dynamics of the system by , a -distributed random variable evolving continuously in time:
| (17) |
where is the exact concentration field, and the integral describes how the gradient field evolves in continuous time. Both processes share the same Brownian motion , indicating that both processes are driven by the same source of randomness.
3 Convergence of the SIPF- Method to Smooth Solutions
We now prove the convergence of the SIPF- method to classical solutions of the 3D parabolic-parabolic Keller-Segel equations. To ensure the validity of the following analysis, we introduce a set of assumptions that impose structure on the concentration fields and their gradients.
3.1 Assumptions and Theorem Statement
To establish the convergence analysis of the SIPF- method, we first introduce assumptions that ensure the boundedness of approximation errors for particles and gradients of the chemical concentration at any time within the temporal computational domain . Specifically, we make the following assumptions.
Assumption 1.
There exist constants such that for all and ,
| (18) | ||||
| (19) |
We note that this assumption only requires the errors to be bounded, but does not require them to converge to zero. The convergence of these errors to zero will be established in the subsequent analysis. Furthermore, the boundedness condition in Eq. (18) can be derived from Eqs. (16)-(17) and Assumption 2(c); Eq. (19) follows directly from the uniform bound in Assumption 2(c) detailed later.
Assumption 2.
Suppose both and satisfy Lipschitz continuity conditions in space and time, along with regularity and boundedness properties as follows:
(a) (Spatial Lipschitz Continuity) There exists a constant , depending on the regularity of and , as well as the parameters and in the system (1), such that for all and ,
| (20) |
This implies that the second derivatives (Hessian entries) exist almost everywhere and satisfy:
| (21) |
(b) (Temporal Lipschitz Continuity) There exists a constant , depending on the regularity of and the parameters and in the system (1), such that for any and ,
| (22) |
(c) (Uniform Boundedness) There exists a constant , depending on the regularity of and the parameters and in the system (1), such that for all and :
| (23) |
(d) (Regularity of Time Derivatives) The exact solution is assumed to be sufficiently smooth in time and space such that both and are bounded. There exists a constant such that for all :
| (24) |
Remark 3.
These regularity and boundedness conditions reflect inherent properties of the SIPF- method rather than external constraints on the KS system. As analyzed in Subsection 3.4, the discrete Fourier representation and Brownian diffusion naturally impose spectral regularity and suppress singular clustering, ensuring remains well-behaved throughout the computation. The boundedness conditions in Assumption 2 are validated through numerical simulations in Subsection 4.1.1.
Assumptions 1 and 2 characterize a regular regime where the exact solution remains smooth over the time interval . We now state our main theorem, which quantifies the convergence of the SIPF- method on .
Suppose that the exact solutions and the numerical solutions obtained by the SIPF- method in satisfy Assumptions 1 and 2 uniformly for all . Let , , , and denote the number of Fourier modes, the number of particles, the batch size, and the uniform time step in the SIPF- method, respectively. Then, for all discrete time levels , the following error estimates hold with high probability:
For all , the -Wasserstein distance (defined in Eq. (78)) between and satisfies , and the maximum error in the truncated Fourier coefficients of and satisfies that . More specifically, for all , the errors are bounded with high probability by:
| (25) | ||||
where , , denote positive constants specified in Eqs. (79)-(80).
A direct consequence of Theorem 4 is that the SIPF- solution converges to the exact solution as and . Specifically, the density converges in the 1-Wasserstein metric even without imposing specific scaling relations between the discretization parameters, as all error terms in the first inequality of Eq. (25) naturally vanish. The convergence of the chemical concentration requires the scaling conditions and to control the statistical and discretization errors, respectively.
Remark 5 (Practical convergence and scaling).
Theorem 4 establishes error estimates up to time , the maximal time at which Assumptions 1-2 hold. Beyond , although the theoretical bounds no longer apply, the Fourier spectral cutoff and the parabolic structure of the equation for act as a natural low-pass filter that prevents the numerical solution from blowing up. Note that this numerical upper bound grows with . By repeating the experiment with increasing values of , we obtain a reliable detection of the time and location of the singularity.
Regarding parameter scaling, the stochastic error term is negligible relative to the leading-order discretization errors. Hence, a moderate batch size (e.g., ) is sufficient in practice, as the total error is dominated by and . Furthermore, while the theoretical gradient error bound scales with due to differentiation, the particle density error benefits from the smoothing effect of time integration. Numerical experiments in Section 4 verify that moderate parameters (, , and ) are sufficient to guarantee convergence and blow-up detection for our numerical examples, which demonstrates the efficiency of the proposed method.
3.2 Several lemmas
The result of Theorem 4 relies on the following lemmas concerning the change in single-step update error of the SIPF- method. We first prove several lemmas in this subsection and leave the complete proof of Theorem 4 to the next subsection.
From Eqs. (3)-(5), the error between and can be decomposed into two components: the error in their Fourier coefficients and the truncation error of . As the number of Fourier modes tends to infinity, and given the smoothness of , the truncation error becomes negligible and can be omitted from the analysis. We now focus on the error analysis between the Fourier coefficients and of and , as presented in the following lemma.
Lemma 6.
Proof.
We first write the frequency , and define the notation , which satisfies as . Recall the update formula for the numerical approximation of the Fourier coefficients (from Eq. (9)):
| (26) |
where represents the Fourier coefficient of at the frequency .
The exact solution satisfies the continuous equation:
| (27) |
Integrating this equation from to and applying the Taylor expansion with integral remainder to the time derivative term yields . This allows us to express the exact solution in a form compatible with the implicit Euler scheme:
| (28) |
By Assumption 2(d), is bounded by . The temporal truncation error satisfies
| (29) |
where is a constant. Taking the Fourier transform of the discrete relation for the exact solution and rearranging terms, we get
| (30) |
To bound the term , we first state a generalization of the mean value theorem to complex-valued functions.
Let be an open subset of , and let be a continuously differentiable function on . Fix points such that the line segment connecting and lies entirely within . There exists such that
| (32) |
The proof of (32) is direct. First, we define the function
| (33) |
Since is also a continuously differentiable function, the mean value theorem implies that there exist points such that
| (34) |
which implies Eq. (32). Applying this result to , we obtain:
| (35) |
Based on Eq. (3.2), we have
| (36) |
Let , where are i.i.d. random variables. This follows from the fact that the particles and are separately i.i.d. Specifically, the i.i.d. property of is ensured by the RBM described in Alg. 2. Based on Assumption 1, is bounded. The empirical mean is defined as and the expectation of is .
According to Bernstein’s inequality, for i.i.d. random variables with (from Assumption 1) almost surely, the probability that the empirical mean deviates from the expectation is bounded as
| (37) |
where . With high probability (e.g., for very small ), the following estimate holds
| (38) |
This implies that, with probability ,
| (39) |
The error estimate between and is more complex than that between and . To analyze this, we introduce an intermediate quantity . Using the frequency notation from Lemma 6, where , we define
| (40) |
where . From Alg. 2, it follows that
| (41) |
where . The error between and can be estimated by
| (42) |
To estimate the error between and , we divide the analysis into two parts:
| (43) |
The first part, involving and , focuses on the spatial discretization error of the continuous convolution on the periodic domain . Specifically, represents the continuous integral, while is its discrete approximation. Since is strictly represented by a truncated Fourier series with modes in , can be viewed as a discrete Fourier projection of the continuous convolution. Thanks to the spatial shift that regularizes the singular gradient kernel on the grid, the quadrature error introduced by this discrete evaluation is of the same order as the Fourier truncation error, and is thus controlled by the spectral projection .
To analyze the error introduced by this approximation, we rely on the following lemma.
Lemma 7.
Proof.
Define the continuous convolution operator . Its Fourier multiplier is , which satisfies for all . Thus, .
The discrete sum evaluates the convolution on a uniform grid with spacing . The spatial shift offsets the quadrature points by half a grid spacing, i.e., . Since the kernel is radially antisymmetric, this half-grid shift regularizes the origin singularity and induces symmetric cancellation. Consequently, the spatial discretization error is dominated by the Fourier truncation error of the spectral projection .
For the truncated modes outside , the frequency magnitude satisfies . Applying Parseval’s identity yields the explicit projection error bound:
| (45) |
Since , the term is controlled by the higher-order spatial regularity of stipulated in Assumption 2(a)(c). Combining these estimates, we obtain
| (46) |
where is a constant.
We now estimate in and . Using the RBM in Alg. 2, we replace with . We write
| (47) |
Lemma 8.
For all , , we have the estimate as follows:
| (48) |
where , is the conserved total mass, is the total number of particles, and is the batch size.
Proof.
Similar to Lemma 3.1 in [24], we rewrite
| (49) |
where denotes the number of times particle is selected in the batch . Since the sampling is with replacement, follows a Binomial distribution with , which indicates that .
Hence,
| (50) |
According to Jensen’s inequality, we obtain:
| (51) |
where . Since all particles are located at distinct positions in the SIPF- algorithm ( for ), there exists a minimum separation distance between any two distinct particles. Consequently, is bounded for all pairs of particles. This ensures that , which is the maximum of these kernel gradient norms, is finite.
3.3 Proof of the Main Theorem
Building on the lemmas established in the previous subsection, we now prove the main theorem in this section. For simplicity of notation in the proof below, we define:
| (52) | ||||
| (53) |
The following provides a bound on the error between and .
Lemma 9.
Proof.
According to Eqs. (16)-(17), it follows that
| (55) |
by the triangle inequality and Tonelli’s theorem. To bound the integrand, we decompose the gradient difference using the triangle inequality:
| (56) |
We now estimate the expectation of each term on the right-hand side of Eq. (3.3). Regarding the first term, recalling the error decomposition in Eqs. (42)-(43), along with the bounds from Lemmas 7 and 8, we have
| (57) |
For the second term, the spatial Lipschitz continuity in Assumption 2(a) yields
| (58) |
where is the Lipschitz constant defined in Assumption 2(a). By the Lipschitz conditions in Assumption 2(a)(b), the exact dynamics in Eq. (17), and the property of the 3D Brownian motion, the temporal variation of the exact gradient satisfies
| (59) |
Now we analyze the error between and as follows.
We next estimate the error between and .
Lemma 10.
For all , under Assumption 2, with high probability,
| (61) |
where , , and are positive constants independent of , , and .
Proof.
Subtracting Eq. (26) from Eq. (30), multiplying by , and using the notation defined in the proof of Lemma 6 give
| (62) |
Since the initial chemical fields coincide, iteration of Eq. (62) yields
| (63) |
Substituting the density-source term in Eq. (63) into the inverse Fourier representation in Eqs. (3)–(5) and evaluating the resulting gradient at gives its contribution to . Applying Eq. (65) together with the decay of the resolvent factor gives, for ,
| (66) |
The Fourier sum in Eq. (3.3) can be estimated by grouping the indices according to . For , the shell contains indices and . Hence
| (67) |
Thus, the series converges. Summing Eq. (3.3) over and using gives
| (68) |
where and .
Finally, we are ready to prove Theorem 4.
Proof of Theorem 4. From Lemma 9 and Lemma 10, we obtain the system of inequalities that couples and defined in Eqs. (52)-(53) as follows:
| (70) | ||||
| (71) |
From this coupled system, we can derive a general bound for . To be specific, substituting Eq. (71) into Eq. (70), we iteratively propagate and simplify the inequality to derive:
| (72) |
where we define .
To bound the deterministic accumulation of the error, we introduce a majorizing sequence that isolates the deterministic components, defined by and
| (73) |
Let be the running maximum of this deterministic sequence. The recursion implies:
| (74) |
By the discrete Gronwall inequality, we obtain the explicit bound for the deterministic accumulation:
| (75) |
where is a constant.
As established in Lemma 8, the RBM gradient error satisfies . This zero-mean property implies that the single-step position noise forms a martingale difference sequence with respect to the natural filtration of the particle system. Consequently, the noise terms are mutually uncorrelated, and thus orthogonal in the space. Therefore, the global accumulation of the RBM noise is bounded by the square root of the sum of its variances:
| (76) |
where is used. Combining this stochastic bound with the deterministic bound yields the refined optimal error bound that holds with high probability:
| (77) |
for all , where higher-order terms are omitted, and , and . Moreover, we point out that is a constant defined in Lemma 7, are constants defined in Lemma 10, are constants defined in Lemma 9.
According to the discrete and continuous dynamics defined in Eqs. (16)-(17), the -Wasserstein distance between the approximate and exact distributions at time is given by:
| (78) |
where the infimum is taken over all possible couplings of the two distributions.
Under the natural coupling induced by shared initial conditions and Brownian motion paths (i.e., and evolve via the same Wiener process ), we explicitly construct a joint distribution . This coupling allows us to bound the Wasserstein distance as:
| (79) |
where the constants are given explicitly by , , , and , with , defined in Lemmas 8-9. Specifically, tracing back to these lemmas shows that , while are independent of .
The inequality follows from the fact that the Wasserstein distance is defined as the infimum over all possible couplings, and our construction provides one such coupling. This step follows from the elementary norm inequality for vectors in , which is a direct consequence of the Cauchy-Schwarz inequality.
To derive the bound for the Fourier coefficients in Eq. (25), we explicitly incorporate the Fourier mode into the error analysis. Unlike the strong trajectory error , the Fourier coefficients represent macroscopic weak functionals (expectations of the test functions ). By the convergence theory of SDEs, the stochastic strong fluctuations of order cancel out in expectation. Instead, the Brownian increments and the zero-mean RBM noise contribute to the macroscopic bias only through second-order Taylor expansion terms, yielding a single-step bias of from temporal discretization and from the RBM approximation.
Accordingly, we refine the accumulation of the source term error over steps (i.e., up to the final time ) by replacing the stochastic part of with its weak bias, while retaining the deterministic spatial part . Applying the standard accumulation estimates along with the frequency bound , we arrive at the error bound:
| (80) |
where we define , , , . Here, , and are constants defined in Lemma 6. We also omit higher-order terms in the leading-order bound. This completes the proof of Theorem 4.
3.4 Numerical Unconditional Stability of the SIPF- Method
The SIPF- method exhibits inherent numerical stability regarding the chemical concentration field. This property stems directly from the spectral discretization and the implicit time-stepping scheme, a property that is independent of the particle distribution. We establish that both the chemical concentration and its gradient remain uniformly bounded for any fixed Fourier mode , thereby justifying the boundedness assumptions employed in the convergence analysis.
Lemma 11 (Unconditional stability of and ).
Let be the finite number of Fourier modes and be the spatial domain. For any time step and the fixed Fourier mode , the reconstructed concentration field and its gradient satisfy the uniform bounds:
| (81) |
and
| (82) |
where are constants depending only on the initial data, , , and the domain size , and is a constant that depends on , , , and , but is independent of the time step size and the particle positions .
Proof.
The update formula for the Fourier coefficients in the SIPF- algorithm is derived from the implicit Euler discretization. By rearranging the terms in Eq. (8) and applying the Fourier transform, the update relates to the previous state and the current empirical density via the amplification factor
| (83) |
The particle source term satisfies
| (84) |
Applying the recursive relation and the triangle inequality yields a bound via geometric series for the magnitude of the coefficients:
| (85) |
We next bound by approximating the partial Fourier sum. Using and treating the sum over as a Riemann sum, which is approximated by an integral in spherical coordinates:
| (86) |
where and are constants.
Regarding the gradient , the SIPF- method employs a specific discretization to handle the singularity of the Green’s function . According to Alg. 2, the gradient at a particle position is computed as:
| (87) |
where the spatial shift ensures that the evaluation point is bounded away from the singularity of Green’s function at grid points . Specifically, let and . The gradient of the Green’s function satisfies
| (88) |
Letting , and using the inequality , we can bound the term to be independent of (and thus independent of ):
| (89) |
Since is bounded (as shown above) and the sum over is finite for a fixed , the first term in Eq. (3.4) is uniformly bounded by a constant depending on but independent of .
For the second term in Eq. (3.4), the RBM excludes self-interaction (). For any fixed grid resolution , the corresponding effective kernel corresponds to a spectrally truncated approximation that is smooth. Thus, for any finite number of particles, this sum is finite.
Combining these results, is guaranteed by the algorithm’s design, where the -norm for the vector-valued gradient is defined by , and is a constant depending on , , , and .
This lemma confirms that the boundedness condition in Assumption 2 is not merely an external hypothesis but a property guaranteed by the SIPF- method itself. This demonstrates that the SIPF- method effectively regularizes the singular Keller-Segel kernel. While the exact solution may exhibit finite-time blow-up (where ), the numerical field remains finite for any fixed Fourier mode . This unconditional numerical stability ensures that the algorithm is robust: it enables the simulation of blow-up phenomena by capturing the solution’s growth trend as increases, while avoiding numerical breakdown at fixed resolutions.
4 Numerical Experiments
The numerical experiments are organized into three main parts to evaluate the SIPF- method. In Subsection 4.1, we validate the SIPF- method by quantifying its accuracy against a high-resolution radial finite difference benchmark and verifying its convergence rates with respect to the time step and batch size. Subsection 4.2 investigates the method’s capability to detect finite-time blow-up phenomena and critical mass thresholds under various conditions. Finally, Subsection 4.3 provides empirical verification of the key theoretical assumptions used in our analysis, specifically the spatial Lipschitz continuity of the concentration gradient.
4.1 Validation of the SIPF- Method
4.1.1 Comparison with FDM
We first demonstrate the accuracy of the SIPF- method. In the radially symmetric case, the fully parabolic KS system (1) in 3D can be expressed as and , where . The system is then rewritten as follows in 1D:
| (90) |
The radial representation in 1D allows us to compute a reference solution using a very fine mesh by the finite difference method (FDM). It will serve as a benchmark to quantify the accuracy of the SIPF- method in 3D. We define the relative error between the cumulative distribution functions (CDFs) obtained from the FDM and the SIPF- method as
| (91) |
where and represent the CDFs of computed via the SIPF- and FDM methods, respectively, and denotes the -th radial mesh point in the FDM, which are the discrete points along the radial direction starting from the origin. To ensure the relative error is well-defined, we set it to zero wherever .
Here, the initial distribution is assumed to be a uniform distribution over a ball centered at with radius 1. The model parameters are chosen as follows:
| (92) |
For the numerical computation, we use Fourier basis functions in each spatial dimension to discretize the chemical concentration and use particles to represent the approximated distribution , where the batch size in Algorithm 2 is . The computational domain is , where , and the total mass is chosen to be . The evolution of and is computed using Algorithm 3 with a time step size , up to the final simulation time .
In Fig. 1, we present the evolution of particles over time, showing the dynamic behavior of . Additionally, in Fig. 2, we compare the cumulative probability curves of obtained from the radial FDM and the SIPF- method at , with a mean relative error of 0.05512 as defined in Eq. (91). This comparison demonstrates that the SIPF- algorithm achieves high accuracy in approximating the true solution. These results validate the effectiveness of the SIPF- algorithm in capturing the behavior of the particle distribution.
4.1.2 Convergence of the SIPF- Method
In this subsection, we validate the convergence of the SIPF- method numerically. Based on Eq. (80), the error between and can be quantified by the error between their Fourier coefficients and . We adopt the same initial conditions as in Subsection 4.1.1. To eliminate the uncertainty introduced by the RBM, the reference solution is computed using the original SIPF method [51] with parameters , , and . Additionally, we set to ensure that the system remains free of singularities, as verified in Fig. 3 of [51].
To investigate the convergence with respect to the time step , we vary from to . Since Theorem 4 holds with high probability, we perform 100 independent experiments for each to empirically validate the algorithm’s accuracy. The mean error of the Fourier coefficients is computed over these 100 trials. As shown in Fig. 3(a), the slope of the mean error versus on a logarithmic scale indicates an approximate first-order convergence rate, with . This result aligns with the theoretical bound given in Eq. (25) of Theorem 4.
Furthermore, we examine the mean error of for varying batch sizes , while keeping . From Eq. (25), with other parameters fixed, the theoretical error of with respect to the batch size should scale as . This is empirically verified in Fig. 3(b), where the fitted convergence rate is , closely matching the theoretical prediction.
4.2 Numerical Detection of Finite-time Blow-up
While the error analysis in Section 3 assumes regular solutions, our algorithm can effectively detect finite-time blow-up phenomena under practical conditions. The parameters are the same as in Eq. (92), with the radial system (90) serving as a reference benchmark. We consider a more concentrated initial condition than in the preceding subsections: a uniform distribution within a ball of radius .
Fig. 4 illustrates the effect of varying the initial mass on solution behavior at . We note that unlike the 2D case, the 3D KS system has no universal critical mass, and the blow-up formation depends on both the total mass and the initial spatial distribution. For our prescribed initial configuration, defined as a uniform density supported over a spherical ball of fixed radius, we numerically observed a sharp blow-up transition. Figs. 4(a) and 4(b) show particle distributions for and , respectively, with the latter exhibiting clear concentration and potential blow-up. In Fig. 4(c), we present the maximum value of over space at for different initial masses, computed via the radial FDM as a benchmark. We denote
| (93) |
The FDM exhibits numerical instabilities (marked in red) for initial masses between and , which identifies the threshold for blow-up specific to this concentrated initial configuration. This agrees with the blow-up transition observed in our proposed method in Fig. 4(b). Remarkably, the concentration and potential blow-up are captured with only particles and a moderate value of , demonstrating the algorithm’s ability to detect singularities without the strict constraints imposed in the convergence analysis.
To further investigate blow-up detection with our method, we examine the maximum chemical concentration as a function of time for different discretization levels and masses . As shown in Fig. 5, when (Fig. 5(a)), the maximum remains stable and shows minimal variation across values, indicating a non-blow-up regime. In contrast, for (Fig. 5(b)), the maximum exhibits strong dependence on , with curves diverging significantly. Notably, this blow-up signature is visible even with coarse discretization ( vs ), demonstrating that singularity detection does not require high-resolution computations. The divergence of solutions at different levels provides a practical diagnostic for intense concentration and blow-up phenomena without demanding highly refined discretizations.
We then demonstrate that our method can also capture concentration and potential blow-up for non-radial initial data. We consider two non-overlapping spheres of radius with centers at and , each containing equal mass (total mass ). Fig. 6 shows the time evolution for two different mass regimes: the top row () remains stable at all times, while the bottom row () exhibits clear concentration towards blow-up formation at . This demonstrates that our algorithm effectively detects blow-up even for non-symmetric initial configurations, highlighting its robustness beyond the radially symmetric case.
4.3 Validation of Theoretical Assumptions
To verify the spatial Lipschitz continuity in Assumption 2(a), we adjust the spatial discretization by varying from 6 to 24. At the final time , we randomly select 1000 pairs of particle points from a total of 10,000 particles in each calculation. The spatial Lipschitz constant for is defined as the maximum ratio of the gradient difference to the spatial distance over all pairs of particle points :
| (94) |
The results, shown in Table 1, list the computed Lipschitz constant for each value of . The variation in these values is relatively small, confirming that the spatial Lipschitz continuity holds for computed by the SIPF- method.
| Fourier modes () | Spatial Lipschitz constant () |
|---|---|
| 6 | 0.002085 |
| 12 | 0.002106 |
| 18 | 0.002036 |
| 24 | 0.001957 |
In previous sections, we developed and analyzed the SIPF- method for the classical 3D fully parabolic KS system (1). The underlying SIPF- framework, however, is flexible and not limited to this classical model. Since the evolution of the chemical field and particle dynamics are updated separately, we can modify the SIPF- method to incorporate additional biological mechanisms and simulate more complicated mathematical biology models [21, 22]. In what follows, we illustrate this flexibility through three representative extensions: multi-species chemotaxis systems with volume-exclusion effects, nonlinear diffusion, and anisotropic cell motility, and we leave their detailed implementation, numerical validation, and convergence analysis for future work.
Multi-species chemotaxis system with volume-exclusion effects. The SIPF- framework is readily extendable to multi-species chemotaxis with volume-exclusion (often referred to as volume-filling) effects, which are essential for modeling realistic biological scenarios where cells of different species compete for finite physical space [19, 37, 3]. A representative two-species chemotaxis model with cross-volume exclusion is
| (95) | ||||
where denotes the total cell density, is the chemical concentration of species , , , and are positive constants, are the production rates of chemical by species , is the density-dependent diffusivity, and is the crowding factor (derived from the volume-filling probability in [37]) that modulates chemotactic mobility. This formulation naturally extends to species by redefining . Crucially, the system (95) reduces to the classical multi-species KS model when and .
Since the equations governing the chemical field remain linear, we apply Algorithm 1 to update the Fourier coefficients of the chemical field as follows:
| (96) |
where denotes the Fourier coefficient of the chemical field , is the time step, is the frequency vector associated with the multi-index (with denoting the domain size), (the index set defined in (4)), and is the Fourier coefficient of the empirical density of species with .
For the density update, species with total mass is represented at time using particles , with the corresponding empirical measure
| (97) |
The Fourier coefficients of the density source in Eq. (96) can be computed directly from this empirical measure as
| (98) |
The density of species is then reconstructed at the Fourier grid points via the inverse Fourier transform, and the total density is obtained by summing these reconstructed densities over all species:
| (99) |
Using the reconstructed total density, the functions and are first evaluated at the Fourier grid points , . The unnormalized Fourier coefficients of these functions are approximated as
| (100) | ||||
This coefficient computation extends Step 5 of Algorithm 1 to handle the density-dependent fields and . The values required for the particle update are obtained via Fourier interpolation:
| (101) | ||||
These quantities, determined by the density at time , are kept fixed during the update from to . Assuming that the frozen diffusivity is sufficiently smooth and nonnegative (after regularization if needed), the Fokker–Planck correspondence leads to the following update for the -th particle of species :
| (102) | ||||
where denote independent standard Gaussian vectors, and is the Green’s function for the operator . For each , the batch consists of indices sampled with replacement from ; when , any occurrence with is omitted from the sum. The bracketed expression approximates , with its density contribution evaluated by random batching. Using Eq. (97), the updated positions yield the empirical densities and hence the source coefficients for the subsequent chemical field update. If is constant and , Eq. (102) reduces to the standard SIPF- update.
Chemotaxis systems with nonlinear diffusion. Density-dependent diffusion is widely used in chemotaxis models to describe changes in cellular dispersal induced by the local population density, particularly when crowding leads to nonlinear motility effects [28, 8, 18, 29, 20]. For the porous-medium diffusion term with , the diffusivity at time can be approximated by . The values of and evaluated at particle positions are then computed from the density reconstruction in Eq. (99) and the Fourier interpolation in Eq. (101). After freezing the regularized over one substep, we formally represent the diffusion equation using the associated Itô SDE, discretize it via the Euler–Maruyama scheme, and then apply the original chemotactic update.
Chemotaxis systems with anisotropic cell motility. Anisotropic motility arises in structured environments such as extracellular-matrix fibers and neural tracts [36, 38]. It can be modeled as
| (103) |
where denotes a sufficiently smooth, symmetric positive semidefinite motility tensor. If chemical diffusion remains isotropic, the spectral update of remains unchanged. The Fokker–Planck correspondence yields the formal particle update
| (104) |
where , denotes the matrix square root, and denote independent standard Gaussian vectors. Thus, anisotropy modifies the particle solver, whereas the field update and random-batch evaluation of preserve the structure of Algorithm 2. We leave the interpolation of and its divergence, together with the convergence analysis of the resulting scheme, for future work in a separate implementation study.
6 Conclusions
In this paper, we introduced a random batch variant [24] of the original SIPF method [51] to simulate the 3D fully parabolic KS system. This modification leverages the randomness in batch sampling to bypass the mean-field limit, which reduces computational complexity without sacrificing accuracy. We established a comprehensive convergence analysis for the fully coupled particle-field system. Specifically, we proved the convergence with high probability for both the density and the concentration field to their respective exact solutions and . The error bounds revealed a dependence on , , and , with the density and concentration fields exhibiting distinct but interrelated convergence behaviors.
Finally, we performed numerical experiments to verify these theoretical rates and confirm the robustness of the SIPF- method. Notably, the algorithm effectively captures density concentration and potential finite-time blow-up phenomena subject to the critical threshold of initial mass in three space dimensions under specific initial data configurations, even with discretization parameters milder than those required by the worst-case theoretical bounds.
Future work will focus on refining the numerical analysis of the classical model and extending the SIPF- method to broader biological applications. From a theoretical perspective, we will refine the error estimates, particularly to weaken the strong dependence on the Fourier mode , and improve the efficiency of the algorithm in high-dimensional settings. Beyond refinements to the classical model, we will build on the formulations presented in Section 5 to develop SIPF- schemes for multi-species chemotaxis systems with volume-exclusion effects, nonlinear diffusion, and anisotropic cell motility, and to investigate their numerical performance and convergence properties.
Acknowledgements
ZZ was partially supported by the National Natural Science Foundation of China (Projects 92470103 and 12171406), the Hong Kong RGC grant (projects 17304324 and 17300325), the Seed Funding Programme for Basic Research (HKU), and the Hong Kong RGC Research Fellow Scheme 2025. ZW was partially supported by NTU SUG-023162-00001 and MOE AcRF Tier 1 Grant RG17/24. JX was partially supported by NSF grant DMS-2309520, the Swedish Research Council grant no. 2021-06594 at the Institut Mittag-Leffler in Djursholm, Sweden, and the E. Schrödinger Institute in Vienna, Austria, during his stay in the Fall of 2025. The simulations were performed on the research computing facilities of the Information Technology Services at the University of Hong Kong, and the Greenplanet Cluster at UC Irvine.
References
- [1] (2015) Toward a mathematical theory of Keller-Segel models of pattern formation in biological tissues. Mathematical Models and Methods in Applied Sciences 25 (09), pp. 1663–1763. Cited by: §1.
- [2] (2019) A chemotaxis-based explanation of spheroid formation in 3D cultures of breast cancer cells. Journal of Theoretical Biology 479, pp. 73–80. Cited by: §1.
- [3] (2006) The Keller–Segel model for chemotaxis with prevention of overcrowding: linear vs. nonlinear diffusion. SIAM Journal on Mathematical Analysis 38 (4), pp. 1449–1475. Cited by: §5.
- [4] (2024) Convergence of random batch method with replacement for interacting particle systems. arXiv preprint arXiv:2407.19315. Cited by: §1, §2.
- [5] (2018) High-order positivity-preserving hybrid finite-volume-finite-difference methods for chemotaxis systems. Advances in Computational Mathematics 44, pp. 327–350. Cited by: §1.
- [6] (2008) A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models. Numerische Mathematik 111, pp. 169–205. Cited by: §1, §1.
- [7] (2016) A blob method for the aggregation equation. Mathematics of Computation 85 (300), pp. 1681–1717. Cited by: §1.
- [8] (2001) A new deterministic spatio-temporal continuum model for biofilm development. Computational and Mathematical Methods in Medicine 3 (3), pp. 161–175. Cited by: §5.
- [9] (2009) Fully discrete analysis of a discontinuous finite element method for the Keller-Segel chemotaxis model. Journal of Scientific Computing 40, pp. 211–256. Cited by: §1.
- [10] (2009) Discontinuous Galerkin methods for the chemotaxis and haptotaxis models. Journal of Computational and Applied Mathematics 224 (1), pp. 168–181. Cited by: §1.
- [11] (2012) Upwind-difference potentials method for Patlak-Keller-Segel chemotaxis model. Journal of Scientific Computing 53, pp. 689–713. Cited by: §1.
- [12] (2013) A study of blow-ups in the Keller–Segel model of chemotaxis. Nonlinearity 26 (1), pp. 81–94. Cited by: §1.
- [13] (2006) A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. Numerische Mathematik 104, pp. 457–488. Cited by: §1.
- [14] (2015) Propagation of chaos for a subcritical Keller-Segel model. In Annales de l’IHP Probabilités et Statistiques, Vol. 51, pp. 965–992. Cited by: §1.
- [15] (2016) Deep learning. MIT press. Cited by: §2.
- [16] (2009) Stochastic particle approximation for measure valued solutions of the 2D Keller-Segel system. Journal of Statistical Physics 135, pp. 133–151. Cited by: §1.
- [17] (1998) Self-similar blow-up for a reaction-diffusion system. Journal of Computational and Applied Mathematics 97 (1-2), pp. 99–119. Cited by: §1.
- [18] (2009) A user’s guide to PDE models for chemotaxis. Journal of Mathematical Biology 58 (1), pp. 183–217. Cited by: §5.
- [19] (2001) Global existence for a parabolic chemotaxis model with prevention of overcrowding. Advances in Applied Mathematics 26 (4), pp. 280–301. Cited by: §5.
- [20] (1995) Dictyostelium discoideum: cellular self-organization in an excitable biological medium. Proceedings of the Royal Society of London. Series B: Biological Sciences 259 (1356), pp. 249–257. Cited by: §5.
- [21] (2024) A stochastic interacting particle-field algorithm for a haptotaxis advection-diffusion system modeling cancer cell invasion. arXiv:2407.05626. Cited by: §5.
- [22] (2026) A novel stochastic particle-field algorithm for a reaction-diffusion-advection cancer invasion model. arXiv:2605.20140. Cited by: §5.
- [23] (2025) Mean field error estimate of the random batch method for large interacting particle system. ESAIM: Mathematical Modelling and Numerical Analysis 59 (1), pp. 265–289. Cited by: §1, §2.
- [24] (2020) Random batch methods (RBM) for interacting particle systems. Journal of Computational Physics 400, pp. 108877. Cited by: §1, §2, §3.2, §6.
- [25] (2021) Random batch methods for classical and quantum interacting particle systems and statistical samplings. In Active Particles, Volume 3: Advances in Theory, Models, and Applications, pp. 153–200. Cited by: §1, §2.
- [26] (1970) Initiation of slime mold aggregation viewed as an instability. Journal of Theoretical Biology 26 (3), pp. 399–415. Cited by: §1.
- [27] (2022) Fourier analysis. Cambridge University Press. Cited by: §1.
- [28] (2005) Preventing blow-up in a chemotaxis model. Journal of Mathematical Analysis and Applications 305 (2), pp. 566–588. Cited by: §5.
- [29] (2007) Global attractors for cross diffusion systems on domains of arbitrary dimension. The Rocky Mountain Journal of Mathematics, pp. 1645–1668. Cited by: §5.
- [30] (2017) Local discontinuous Galerkin method for the Keller-Segel chemotaxis model. Journal of Scientific Computing 73 (2), pp. 943–967. Cited by: §1, §1.
- [31] (2018) Positivity-preserving and asymptotic preserving method for 2D Keller-Segel equations. Mathematics of Computation 87 (311), pp. 1165–1189. Cited by: §1.
- [32] (2017) A random particle blob method for the Keller-Segel equation and convergence analysis. Mathematics of Computation 86 (304), pp. 725–745. Cited by: §1, §1.
- [33] (2019) Propagation of chaos for the Keller-Segel equation with a logarithmic cut-off. Methods and Applications of Analysis 26 (4), pp. 319–348. Cited by: §1, §1.
- [34] (2013) Theory of functions of a complex variable. American Mathematical Society. Cited by: §1.
- [35] (1995) Blow-up of radially symmetric solutions to a chemotaxis system. Advances in Mathematical Sciences and Applications 5, pp. 581. Cited by: §1.
- [36] (2002) The diffusion limit of transport equations II: chemotaxis equations. SIAM Journal on Applied Mathematics 62 (4), pp. 1222–1250. Cited by: §5.
- [37] (2002) Volume-filling and quorum-sensing in models for chemosensitive movement. Canadian Applied Mathematics Quarterly 10 (4), pp. 501–543. Cited by: §5, §5.
- [38] (2009) Continuous models for cell migration in tissues and applications to cell sorting via differential chemotaxis. Bulletin of Mathematical Biology 71 (5), pp. 1117–1147. Cited by: §5.
- [39] (2019) Mathematical models for chemotaxis and their applications in self-organisation phenomena. Journal of Theoretical Biology 481, pp. 162–182. Cited by: §1.
- [40] (1953) Random walk with persistence and external bias. The Bulletin of Mathematical Biophysics 15, pp. 311–338. Cited by: §1.
- [41] (2006) Transport equations in biology. Springer Science & Business Media. Cited by: §1.
- [42] (2005) Notes on finite difference schemes to a parabolic-elliptic system modelling chemotaxis. Applied Mathematics and Computation 171 (1), pp. 72–90. Cited by: §1.
- [43] (2007) Conservative upwind finite-element method for a simplified Keller-Segel system modelling chemotaxis. IMA Journal of Numerical Analysis 27 (2), pp. 332–365. Cited by: §1.
- [44] (2011) Error analysis of a conservative finite-element approximation for the Keller-Segel system of chemotaxis. Communications on Pure and Applied Analysis 11 (1), pp. 339–364. Cited by: §1.
- [45] (2020) Unconditionally bound preserving and energy dissipative schemes for a class of Keller-Segel equations. SIAM Journal on Numerical Analysis 58 (3), pp. 1674–1695. Cited by: §1.
- [46] (1994) Chemotaxis and chemokinesis in eukaryotic cells: the Keller-Segel equations as an approximation to a detailed model. Bulletin of Mathematical Biology 56 (1), pp. 129–146. External Links: ISSN 0092-8240 Cited by: §1.
- [47] (2019) Optimal transport: fast probabilistic approximation with exact solvers. Journal of Machine Learning Research 20 (105), pp. 1–23. Cited by: §2.
- [48] (2000) The derivation of chemotaxis equations as limit dynamics of moderately interacting stochastic many-particle systems. SIAM Journal on Applied Mathematics 61 (1), pp. 183–212. Cited by: §1, §1.
- [49] (2022) DeepParticle: learning invariant measure by a deep neural network minimizing Wasserstein distance on data generated by an interacting particle method. Journal of Computational Physics 464, pp. 111309. Cited by: §2.
- [50] (2024) A DeepParticle method for learning and generating aggregation patterns in multi-dimensional Keller-Segel chemotaxis systems. Physica D: Nonlinear Phenomena 460, pp. 134082. Cited by: §2.
- [51] (2025) A novel stochastic interacting particle-field algorithm for 3D parabolic-parabolic Keller-Segel chemotaxis system. Journal of Scientific Computing 102 (3), pp. 1–23. Cited by: §1, §1, §2, §4.1.2, §6.
- [52] (2024) Randomized methods for computing optimal transport without regularization and their convergence analysis. Journal of Scientific Computing 100 (2), pp. 37. Cited by: §2.
- [53] (2026) A Bidirectional DeepParticle method for efficiently solving low-dimensional transport map problems. Journal of Computational Physics 561, pp. 114983. Cited by: §2.
- [54] (2017) Finite volume methods for a Keller-Segel system: discrete energy, error estimates and numerical blow-up analysis. Numerische Mathematik 135 (1), pp. 265–311. Cited by: §1, §1.