text:
\section{Introduction}

Chemotaxis partial differential equations (PDEs) were introduced by Keller and Segel (KS \cite{keller1970initiation}) to describe the aggregation of the slime mold amoeba Dictyostelium discoideum due to an attractive chemical substance. Related random walk model by Patlak
was known earlier \cite{Patlak_53}, see \cite{Oth_97} for
an analysis of basic taxis behaviors (aggregation, blowup and collapse) based on reinforced random walks.
Recall a common form of KS model \cite{fatkullin2012study}:
\begin{align}
    \rho_{t} = \nabla \cdot (\mu \, \nabla \rho - \chi\, \rho\,  \nabla c ), \;\; \;\;
    \epsilon \, c_t  = \Delta\, c - k^2 \, c + \rho, \label{KS1}
    \end{align}
    where $\chi, \mu $ ($\epsilon, k$) are positive (non-negative) constants. The model is called elliptic if $\epsilon =0$ (when $c$ evolves rapidly to a local equilibrium), and parabolic if $\epsilon >0$.
The $\rho$ is the density of active particles (bacteria), and $c$ is the concentration of
chemo-attractant.
The bacteria diffuse with mobility $\mu$ and drift in the direction of $\nabla c$ with velocity $\chi \nabla c$ where $\chi$ is called chemo-sensitivity.
\medskip

In the simplest regime
 $(\epsilon,k)=\mathbf{0}$, the concentration equation

 becomes the Poisson equation
$  - \Delta c = \rho$. Subject to a suitable boundary condition of $c$, the classical integral representation
$c=- \mathcal{K}\ast \rho$ holds,
where $\mathcal{K}$ is the Green's function of the Laplacian, and $\ast$ denotes convolution.
The KS system
then reduces to a scalar non-local nonlinear advection-diffusion PDE governing the evolution of the density function $\rho$:
\begin{align}
\rho_t = \mu \, \Delta \, \rho + \chi \, \nabla \cdot \big(\rho \, \nabla( \mathcal{K}\ast \rho)\big). \label{KSequation1}
\end{align}

For modeling chemotaxis in a fluid environment such as ocean,  \cite{chemomix_12,chemomix_16,chemomix_yao,chemomix_21,chemomix_22} studied equation
(\ref{KSequation1}) with the advective Lie derivative $\rho$ on the left hand:
\begin{align}
\rho_t +\nabla \cdot(\rho \, \textbf{v})= \mu \, \Delta \, \rho + \chi \, \nabla \cdot \big(\rho \, \nabla( \mathcal{K}\ast \rho)\big). \label{KSequation2}
\end{align}
The mixing mechanism of the flow field $\mathbf{v}$ is known to slow down or smooth out blowup or aggregation in (\ref{KSequation1}), see analysis in  \cite{chemomix_12,chemomix_16,chemomix_21,chemomix_22} and references therein, and related convection induced smoothing in fluids
\cite{hou2006dynamic,hou2009stabconv}.
Eq. \eqref{KSequation2} is the macroscopic limit (McKean-Vlasov equation) of the interacting particle system below as $J \uparrow \infty$:
\begin{align}
    d X^{j} = -\,  {\chi M \over J} \nabla_{X^j} \,
    \sum_{i=1:J,\, i\not = j}\, \mathcal{K} (|X^j - X^i |) \, dt + \textbf{v}(X^j) dt
+ \sqrt{2\,\mu} \, d W^j, \;\; j=1,\cdots, J,   \label{IPS1}
\end{align}
where $M$ is the conserved total mass (integral of $\rho$), $W^j$'s are independent Brownian motions in
$R^d$.
\medskip

In this work,
$ \textbf{v}$ is a prescribed divergence
free vector field. We shall approximate
$\rho$ of equation (\ref{KSequation2}) numerically based on the associated interacting particle system in two and three spatial dimensions ($d=2,3$), and carry out a systematic deep learning study (a.k.a. DeepParticle \cite{DP_22}) on the $(\mu, \textbf{v})$ dependency of solutions.
As we are interested in studying  near singular KS solutions, the main challenge for training data collection is to approximate the fields $\rho(\textbf{x})$ and $c(\textbf{x})$ reliably as they intensify. Due to singular behavior of Green's function as particles come close to each other, our approach is to regularize $\mathcal{K}$ in  (\ref{IPS1}) for approximating  $\rho(\textbf{x})$ as particles aggregate,
in similar spirit to
the vortex blob method for fluids (\cite{krasny1991vortex} and references therein) and \cite{blob_aggregation_ck2016}.
\medskip

Deep learning tools have been applied broadly for scientific computing in recent years, such as solving PDEs and their inverse problems, see \cite{Osher_21,DP_22} and references therein.
DeepParticle \cite{DP_22} is based on a particle
method for solving a time-dependent physically parameterized PDE, whose
solution is approximated by the particle empirical measure (distribution).  A deep neural
network (DNN) with physical parameter dependence
learns the mapping from the initial particle distribution to
the particle distribution at time $T$ with training data over sampled
physical parameters provided by the particle method.
The trained DNN then generates approximate solutions at time $T$
for new physical parameters unseen in training.
DeepParticle has been successfully designed and trained for learning
and generating invariant measures of stochastic interacting particle
systems (at $T=\infty$) arising in reaction-diffusion front speeds in three dimensional chaotic flows \cite{DP_22}. In this paper, we further develop DeepParticle to learn and generate KS solutions exhibiting aggregation behavior
at a finite time $T$ before blow-up for a range of diffusivity
$\mu$ and advection amplitude values.
\medskip

The rest of the paper is organized as follows. In Section 2, we briefly review the blow up behavior in the KS model and the regularized  particle methods to solve the KS model. In Section 3, we present our DeepParticle method to learn the transport map from an input distribution to a target distribution. Moreover, we will discuss details of implementation of our method and how to learn the distributions in particle simulation of the KS model.  In Section 4, we show numerical results to demonstrate the performance of our method. Concluding remarks are made in section 5.

\section{Regularized interacting particle method of KS equation}

\subsection{Blow up behavior in KS model}
We start from the simplest KS model without advection, namely \eqref{KSequation1} with $\mu=\chi=1$. This system has been extensively studied by many authors, see survey article \cite{perthame2004pde}. The conservation of mass holds:
\begin{align}
\frac{d}{ dt}\int_{\mathbb{R}^{2}} \rho(x,t) d x = 0,
\end{align}
which is also true for  Eq. \eqref{KSequation2} when the advection field $\textbf{v}$ is divergence free.  If we set the total mass $ M:=\int_{\mathbb{R}^{2}} \rho(x,0) d x$, the second moment has fixed time derivative, i.e.,
\begin{align}\label{second_order_moment}
    \frac{d}{ dt} \int_{\mathbb{R}^{2}}|x|^{2} \rho(x,t)  x=\frac{M}{2 \pi}(8 \pi-M),
\end{align}
where $8 \pi$ is called critical mass of the system. Accordingly, it is well-known that: 1) if $M>8\pi$, the system has no global smooth solutions; 2) if $M=8\pi$, the system has a global smooth solution, which blow up as $t\to\infty$; 3) if $M<8\pi$, the system has a global smooth solution.
\medskip

In \cite{chemomix_yao}, an extra advection term is introduced to the KS equation, which reflects the organism movement in some prescribed fluid flows. Then the second moment identity in Eq.\eqref{second_order_moment} and the subsequent blowup vs. global evolution results are no long valid. By comparison principle, \cite{chemomix_yao} shows  that if the total mass is smaller than critical mass, the system has global smooth solution with smooth initial data. In the cases with supercritical mass, there are only numerical experiments suggesting that the advection, if sufficiently large,  prevents the solutions from blowing up. Later, \cite{chemomix_22} shows that when the flow exerts a  `stretching' effect, the advection field $\textbf{v}$ indeed suppresses the growth or the concentration of chemo-attractant and hence the solution has global regularity and exists for all time.  Examples of stretching flows include hyperbolic flows where $\mathbf{v}(\mathbf{x})=\left( x_{1}, -\frac{1}{d-1}\mathbf{x}_{-}\right)$ and laminar flows where $\mathbf{v}(\mathbf{x})=\left(v\left(\mathbf{x}_{-}\right), \mathbf{0}_{-}\right)$, with $\mathbf{x}_{-}=\left(x_{2}, \ldots, x_{d}\right)$ and $v$ is periodic. However the case of chaotic flows or when the amplitude of advection $\textbf{v}$ is not sufficiently large, remains open.

\subsection{Regularized interacting particle methods}

The singular behavior of the governing PDE and the associated Green's function that cause trouble for particle methods is a well-known problem. The vortex blob method \cite{chorin1973discretization} provides a regularization approach to extend vortex sheet motion past the singularity formation time into the physically important roll-up regime, in which the points representing the vortex sheet are replaced by vortices of prescribed and fixed shape. Numerical calculations show regular motion for the centers of the blobs even after the time when a curvature singularity on a vortex sheet is formed. Later, a special form of the vortex blob method \cite{krasny1986desingularization} is used to calculate the roll-up of a periodic vortex sheet resulting from the classical Kelvin-Helmholtz instability. The singular solutions are closely related to the ill-posedness of vortex sheet motion \cite{caflisch1989singular} . Nonetheless,
as the regularization parameter approaches zero, the regularized vortex sheet solution converges to a weak solution of the Euler equation \cite{liu1995convergence}.
 As in vortex blob methods, we formulate a regularized interacting particle method (denoted as IPM from here on) to solve KS chemotaxis systems.
\medskip

 Let us represent the density function (solution of Eq.\eqref{KSequation2}) with empirical distributions of particle positions. In SDE \eqref{IPS1} on particle positions,
the chemo-attractant term $\,  {\chi M \over J} \nabla_{X^j} \,
    \sum_{i=1:J,\, i\not = j}\, \mathcal{K} (|X^j - X^i |) \, dt$ causes numerical instabilities when particles tend to concentrate. To overcome this difficulty, we replace the singular kernel in \eqref{IPS1} by a smoothed approximation $K_{\delta}$,
such that $K_{\delta}(z)\rightarrow K(z)$ as $\delta \rightarrow 0$, where $\delta>0$ is a regularization parameter. For example, we define
\begin{align}
    K_{\delta}(z) = K(z)\frac{|z|^2}{|z|^2+\delta^2}.
\end{align}
\medskip

Equipped with kernel $K_{\delta}$, we
obtain a system of regularized SDEs for the particles as follows:
\begin{align}
    d X^{j} = -\,  {\chi M \over J} \nabla_{X^j} \,
    \sum_{i=1:J,\, i\not = j}\, K_{\delta}(|X^j - X^i |) \, dt+ \textbf{v}(X^j) dt
+ \sqrt{2\,\mu} \, \, d W^j, \;\; j=1,2,\cdots, J,
    \label{IPS2}
\end{align}
where $J$ is the number of particles, $X^{j}\in R^d$ is the position of the $j$-th particle, where $d W^j$ are mutually independent $d-$dimensional Brownian motions. The convergence of
a random particle blob method, similar to (\ref{IPS2}) yet for KS without advection field $\textbf{v}$, is analyzed in
\cite{rpbm_2017}. For a stochastic particle method using heuristic collision and splitting rules to bypass the singular behavior of Green's function, see \cite{havskovec2009stochastic}.
\medskip

Representing PDE solutions by particles belongs to the Lagrangian framework.  The Lagrangian methods have several advantages: (1) easy to implement; (2) spatially mesh-free and self-adaptive; and (3)  computational costs scale linearly with the dimension of spatial variables in the underlying stochastic dynamical systems. If we discretize KS system  \eqref{KSequation2} on mesh grids with e.g. finite element \cite{carrillo2019hybrid} and spectral \cite{shen2020unconditionally} methods, the number of mesh grids depends exponentially on the spatial dimension. The key benefit of the Lagrangian framework in computing KS models
is its natural capability to follow the KS solution when a singular behavior is emerging. As we shall see, the stochastic particle method based on (\ref{IPS2}) reproduces the well-known aggregation behavior and captures the KS dynamics during the potential blow-up stage of evolution. This is another step forward in our program of computing high gradient solutions in the Lagrangian framework, which has shown encouraging results for a range of multi-dimensional singularly-perturbed  advection-diffusion PDEs.
We refer interested readers to our recent progress in developing Lagrangian methods to compute effective diffusivites in chaotic or random flows \cite{WangXinZhang:18,IPM_2020,SharpMMS_21,TDChaoticFlow_22} and KPP front speeds in chaotic flows \cite{IPM_2021}.
There are also deterministic particle methods (\cite{blob_diff_2019,blob_aggregation_ck2016} and references therein) for a class of KS and degenerate diffusion equations that fall in our DeepParticle framework (section 3).
\medskip

Though the IPM in our study here is mesh-free and self-adaptive for solving multi-dimensional KS chemotaxis systems, the computational costs remain high if we want to study the systems under a variation of parameters (evolution time $T$ and amplitude of advection $A$). Also, the numbers of particles $J$ cannot be too large as the chemo-attractant term has $O(J^2)$ complexity in each step. The total complexity is then $O(J^2\frac{T}{\Delta t})$ where $\Delta t$ is the time step of discretization. On the other hand, our numerical simulation of Eq.\eqref{IPS2} shows that the distribution at finite time $T$, starting from the same initial distribution, may have continuous dependence on the physical parameters. Therefore in the next section, we will introduce a machine learning algorithm to learn the feature of continuous dependency on parameters and generate approximated samples with $O(J)$ complexity.

\section{Deep particle method}
In this section we introduce a DeepParticle algorithm to learn the features of the transport map from a trivial (input) distribution to a target (output) distribution. The mapping error is measured by the 2-Wasserstein distance.

\subsection{Discrete Wasserstein distance}
\noindent	Given distributions $\mu$ and $\nu$ defined on metric spaces $X$ and $Y$, let us construct a transport map $f^0:X \to Y$ such that
		$f^{0}_{*}(\mu)= \nu$, where star denotes the push forward of the map.
For any function $f:X \to Y$, the $2$-Wasserstein distance between  $f_*(\mu)$ and $\nu$ is defined by:

\begin{align}\label{def:Wdistance}
	W_{2}(f_*(\mu), \nu):=\left(\inf _{\gamma \in \Gamma(f_*(\mu), \nu)} \int_{Y \times Y} dist(y', y)^{2} \mathrm{~d} \gamma(y', y)\right)^{1 / 2},
\end{align}
where $\Gamma(f_*(\mu), \nu)$ denotes the collection of all measures on $Y \times Y$ with marginals $f_*(\mu)$ and $\nu$ on the first and second factors respectively and $dist$ denotes the metric (distance) on $Y$. A
straightforward derivation yields:

\begin{align}\label{def:Wdistance1}
    W_{2}(f_*(\mu), \nu)=\left(\inf _{\gamma \in \Gamma(\mu, \nu)} \int_{X \times Y} dist(f(x), y)^{2} \mathrm{~d} \gamma(x, y)\right)^{1 /2},
\end{align}
where $\Gamma(\mu, \nu)$ denotes the collection of all measures on $X \times Y$ with marginals $\mu$ and $\nu$ on the first and second factors respectively and still $dist$ denotes the metric (distance) on $Y$.
We approximate $\mu$ and $\nu$ by empirical distribution functions:

$	\mu=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_i}$,
$	\nu=\frac{1}{N}\sum_{j=1}^{N}\delta_{y_j}$, where $N$ is the number of samples of distribution. Under the setting of learning distribution from interacting particle methods, we take $N<J$ sub-samples from terminal time position of system \eqref{IPS1} to represent the distribution. We will re-sample them every $1000$ steps of training. The preceding technique is usually referred as \emph{mini-batch} in the machine learning literature.
\medskip

	Any joint distribution in $\Gamma(\mu,\nu)$ can be represented by an $N\times N$ doubly stochastic matrix \cite{sinkhorn1964relationship}, denoted as transition matrix,  $\gamma=(\gamma_{ij})_{i,j}$ satisfying:

	\begin{align}
	\label{def:bistochasticity}
		\gamma_{ij} \geq 0; \quad \quad
	 \forall j, ~\sum_{i=1}^{N}\gamma_{ij}=1;
	 \quad \quad
	 \forall i, ~\sum_{j=1}^{N}\gamma_{ij}=1.
	\end{align}

Then \eqref{def:Wdistance1} becomes
\begin{align}\label{def:Wdistance_approx}
\hat{W}(f):=\left(\inf_{\gamma \in \Gamma^N}\sum_{i,j=1}^N\, dist(f(x_i),y_j)^2 \gamma_{ij}\right)^{1/2}.
\end{align}
$\hat{W}(f)$ in \eqref{def:Wdistance_approx} has a simple intuitive interpretation: given a $\gamma \in \Gamma(\mu, \nu)$ and any pair of locations $(x,y)$, the value of $\gamma(x,y)$ tells us what proportion of $f_*(\mu)$ mass at $f(x)$ should be transferred to $y$, in order to reconfigure $f_*(\mu)$ into $\nu$. Computing the effort of moving a unit of mass from $y'$ to $y$ by $dist(f(x),y)^{2}$ yields the interpretation of $\hat{W}$ as the minimal effort (optimal transportation \cite{villani2021topics}) to reconfigure $f_*(\mu)$ mass distribution into that of $\nu$.

\subsection{Training data and network configuration}
  Note that given any fixed set of $\{x_i\}_{i=1}^N\subset R^d$ and $\{y_j\}_{j=1}^N\subset R^d$ (training data), we have derived Eq.\eqref{def:Wdistance_approx} to be minimized by gradient descent. In addition, we aim to find a network that can represent the change of target distribution over some physical parameters. In such scenario, more than one set of data ($\{x_i\}$ and $\{y_j\}$ consists one set of data) should be assimilated. More precisely, let the total number of data set be denoted as $N_{dict}$. Then we have $N_{dict}$ pairs of i.i.d. samples of input and output distribution, denoted as $\{x_{i,r}\}$ and $\{y_{j,r}\}$ for $r=1\cdots N_{dict}$. Associated with the $r$-th data set ($\{x_{i,r}\}$ and $\{y_{j,r}\}$), we assume that there is a physical parameter $\eta_r\in R^p$. To represent this in the network, we encode $\eta_r$ to each data in the set, i.e., the input of the network is $\{(x_{i,r},\eta_r)\}_i \subset R^{d+p}$. This procedure is also called \textit{padding} in the literature.
 The output of network is then denoted as $f_\theta(x;\eta)$ where $\theta$'s are all trainable parameters of the network.
 \medskip

 Between the input ($l_0$) and output ($l_6$), there are 5 latent layers where each layer is $30$ in width. The adjacent layers $l_i$ and $l_{i+1}$ are fully connected, i.e.,
 \begin{align}
     l_{i+1}=\tanh(W_i(l_{i})+b_i) \quad i=0,\cdots 4,
 \end{align}
 where $W_i$ is (weight) matrix width of $l_{i+1}$ and $l_{i}$; $b_i$ is the (bias) vector with same dimension as $l_{i}$; $\tanh$ is the activation function. For the output layer, the formula is similar except we do not apply the activation function.  In case of $d=2$ and $p=1$ (e.g. the first numerical example in computing blowup behaviour of KS without advection, see Sec.\ref{sec:eg1}), there are in total $4966$ parameters (weight and bias) to update by gradient descent. In 3D cases, we find that the network performs well even without altering the width of latent layers.

 \begin{remark}
 Comparing with par-net approach of \cite{DP_22}, we found that  \textit{padding} achieves similar performance when learning distribution in particle simulation of KS model.
 \end{remark}

  In our computation, by equipping $\mathbb{R}^d$ with Euclidean metric, the training loss function is
\begin{align}\label{def:penalty-1}
	\hat{W}^2(f_\theta):=\sum_{r=1}^{n_{\eta}}\left(\inf_{\gamma_r \in \Gamma^N}\sum_{i,j=1}^N|f_\theta(x_{i,r};\eta_r)-y_{j,r}|^2\gamma_{ij,r}\right).
\end{align}

The goal of deep particle algorithm is then to find $(\theta,\{\gamma_r\})$ to minimize
\begin{align}\label{def:goal}
	P(\theta,\{\gamma_r\}):=\sum_{r=1}^{n_{\eta}}\left(|f_\theta(x_{i,r};\eta_r)-y_{j,r}|^2\gamma_{ij,r}\right).
\end{align}

\subsection{Iterative method in finding transition matrix $\gamma$}
Notice that, with fixed $\theta$, finding the transition matrix $\gamma$ to calculate discrete Wasserstein distance in \eqref{def:penalty-1} is linear programming with degree of freedom $N^2$. When the number of particles ($N$) becomes large, it is expensive to go by a conventional algorithm, e.g. interior point algorithm or simplex algorithm.  We now present a {\it mini-batch linear programming algorithm} to find the best $\gamma$ for each
 inner sum of \eqref{def:penalty-1}, while suppressing $k_r$ dependence in $f_\theta$.

\medskip

The problem \eqref{def:penalty-1} is a linear program on the bounded convex set $\Gamma^N$ of vector space of real $N\times N$ matrices. By Choquet's theorem, this problem admits solutions that are extremal points of $\Gamma^N$. The set of all doubly stochastic matrix $\Gamma^N$ is viewed as Birkhoff polytope. The Birkhoff–von Neumann theorem \cite{schrijver2003combinatorial} states that such polytope is the convex hull of all permutation matrices, i.e., those matrices such that $\gamma_{ij}=\delta_{j,\pi(i)}$ for some permutation $\pi$ of $\{1,...,N\}$, where $\delta_{jk}$ is the Kronecker symbol.
 \medskip

 Our algorithm is defined iteratively. We start from a permutation matrix, e.g., $\gamma_{ij}=\delta_{ij}$. In each iteration, we randomly select columns and corresponding rows such that the submatrix is a permutation matrix. Then the entries of submatrix consist of  a linear programming problem under the constraint that maintains column-wise and row-wise sums equal to one. To be precise, we randomly choose  $\{i_k\}_{k=1}^M$, ($M\ll  N$)  from $\{1,2, \cdots, N\}$ without replacement. Then $j_k$ is selected such that $\gamma_{i_k,j_k}=1$.  The cost function of the sub-problem is
 \begin{align}\label{def:lp-sub-cost}
	C(\gamma^*):=\sum_{k,l=1}^M|f_\theta(x_{i_k})-y_{j_l}|^2\gamma^*_{i_kj_l}
\end{align}
subject to
\begin{align}\label{def:lp-sub-cons}
\left\{	\begin{array}{l}
		\sum_{k=1}^M \gamma^*_{i_k,j_l}=1\quad\forall l=1,\cdots,M\\
			\sum_{l=1}^M \gamma^*_{i_k,j_l}=1\quad\forall \, k=1,\cdots,M\\
	\gamma^*_{i_kj_l}\geq 0 \quad \forall \, k,l = 1,\cdots,M.
\end{array}\right.
\end{align}
Then $\gamma^*$ is again a permutation matrix.  The linear programming sub-problem of much smaller size is solved by the interior point method  \cite{wright1997primal}. In addition, as the goal is to find a permutation matrix, we set the tolerance of interior point to be relatively large and project the resulting approximation to a permutation matrix. In our approach, we terminate the iteration of sub-problem until $\min_k(\max_l \gamma_{i_k,j_l}^*)>0.5$. This is to ensure that there is a unique large-value entry in each column. As a projection, we update $\gamma$ by,
\begin{align}
    \gamma_{ij}=\left\{ \begin{array}{l}
         1 \quad \text{if}\ j=\arg\max_l \gamma^*_{il}, \\
         0 \quad \text{otherwise.}
    \end{array}\right.
\end{align}

We observe that the global minimum of $\gamma$ in  \eqref{def:penalty-1} is also the solution of sub-problems \eqref{def:lp-sub-cost} with arbitrarily selected rows and columns, subject to the row and column partial sum values of the global minimum. The selection of rows can be one's own choice. In our approach, in each step after gradient descent, we randomly select rows to solve the optimization problem iteratively.
\medskip

Comparing with Wasserstein GAN proposed in \cite{arjovsky2017wasserstein}, \eqref{def:goal} is a Min-Min optimization. Both Adam gradient descent of $\theta$ and mimi-batch optimization of $\gamma$ are iteratively defined. We then alternatively update $\theta$ and $\gamma$ to seek a global minimum of \eqref{def:goal}.
\medskip

The cost of finding optimal $\gamma$ increases as $N$ increases, however the network itself is independent of $\gamma$. After training, our network acts as a sampler from  some target distribution $\nu$
without assumption of closed-form distribution of $\nu$. At this stage, the input data is no longer limited by  training data, an arbitrarily large amount of samples approximately obeying $\nu$ can be generated through $\mu$ (uniform distribution).

  \subsection{Full Training Algorithm}
   The full training process is outlined in Alg. \ref{alg:1},
   and carried out on a quad-core CPU desktop with an RTX2070 8GB GPU at UC Irvine. The training data is collected from the first order explicit IPM which  solves the regularized SDE system (\ref{IPS2}) by Euler's method in time. IPM is also the reference solver for evaluating DeepParticle generations.

\begin{algorithm}[htbp]
	\SetAlgoLined
		Randomly initialize weight parameters  $\theta$ in network $f_\theta:\mathbb{R}^d\to \mathbb{R}^d$; \\
	\Repeat{given number of training mini-batches, $N_{dict}$}{
	\For{ physical parameter set $r\gets0$ \KwTo $n_{\eta}$ }{
	randomly select $\{x_{i,r}\}$, $\{y_{j,r}\}$, $i,j=1:N$ from i.i.d. samples of input and target distribution (generated by IPM) with respect to physical parameter $\eta_r$;\\
	$\gamma_{ij,r}=\delta_{i,j}$, i.e., initialize as permutation matrix;\\
	}
	\If{not the first training mini-batch}{
		\For{ physical parameter set $r\gets0$ \KwTo $n_{\eta}$ }{
		$\delta P_r=+\infty$\\
			\While{$|\delta P_r|_{fro}<tol$}{
		$P_r=\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}$;\\
		randomly choose $\{i_{k,r}\}_{k=1}^M$, from $\{1,2, \cdots, N\}$ without replacement; \\
		find $\{j_{k,r}\}_{k=1}^M$, such that $\gamma_{i_{k,r}j_{k,r},r}=1$\\
		solve the linear programming
		sub-problem \eqref{def:lp-sub-cost}-\eqref{def:lp-sub-cons}
		to get $\gamma_r^*$ another permutation matrix;\\
		update $\{\gamma_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$ with $\{\gamma^*_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$.\\
		$\delta P_r=\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}^*-P_r$
	}
	}
	}
	\Repeat{given steps for each training mini-batch}{
		$P=\sum_{r=1}^{N_r}\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}$;\\
		$\theta\leftarrow \theta - \delta_1\nabla_{\theta} P$, $\delta_1$ is the learning step size;\\
		\Repeat{given linear programming steps, $N_{LP}$}{
					\For{ physical parameter set $r\gets0$ \KwTo $n_{\eta}$ }{
		$P_r=\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}$;\\
		randomly choose $\{i_{k,r}\}_{k=1}^M$, from $\{1,2, \cdots, N\}$ without replacement; \\
		find $\{j_{k,r}\}_{k=1}^M$, such that $\gamma_{i_{k,r}j_{k,r},r}=1$\\
		solve the linear programming
		sub-problem \eqref{def:lp-sub-cost}-\eqref{def:lp-sub-cons}
		to get $\gamma_r^*$ another permutation matrix;\\
		update $\{\gamma_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$ with $\{\gamma^*_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$.\\
	}
	}
	}
	}
	$\mathbf{Return}$
	\caption{DeepParticle Learning}\label{alg:1}
\end{algorithm}
\section{Numerical Examples}
\subsection{2D KS Simulation and Generation in the Absence of Advection}\label{sec:eg1}
First we consider the KS model without advection, namely \eqref{KSequation1} with $\mu=\chi=1$.
A straightforward derivation shows if the initial mass $M>8\pi$ and has finite second moment, the system will blow up in finite time.
\medskip

As the first numerical example, we consider learning the change of distribution depending on evolution time $T$ starting from the initial distribution. The initial distribution is assumed to be a uniform distribution on a ball with radius $1$ centered at the origin. Assuming the total mass is $16 \pi$, we have then the initial second moment as $8\pi$, same as in \cite{chemomix_yao}. By Eq.\eqref{second_order_moment}, the system will blow-up when $t>0.125$. Before applying the deep learning algorithm, we first investigate the error of regularized method. In Fig.\ref{fig:eg1_training}(a), we show the second moment of IPM simulation  with various regularization factors. We see that except for the value  $\delta^2=10^{-2}$, the smaller  regularization factor $\delta$ does not affect the accuracy of the IPM simulation when $t\in[0,0.1]$.
\medskip

Then we turn to investigate the performance of the proposed DeepParticle algorithms. For training data, we applied IPM with regularization $\delta^2=1e-3$ for $T=0.1$ with $J=10000$ particles. We keep the snapshots of the empirical distribution every $0.05$ interval. During training, we consider a mini-batch of size $8\times 2000$, which means, we take $N_{dict}=8$ sets of training data at various time $t$ (which replaces the general notation for physical parameter $\eta$ in the network expression $f_\theta(x,\eta)$) and in each set (mini-batch) we have $N=2000$ subsamples from $J=10000$ samples from IPM. Adams gradient descent is applied to learn the parameters (weight and bias) in the network. We renew the mini-batch every $1000$ steps and renew $\gamma$ every $200$ steps.
In Fig.\ref{fig:eg1_training}(b), we show the training loss computed by Eq.\eqref{def:penalty-1}. Since we do not renew $\gamma$ at every steps of gradient descent of parameters, the training loss is not uniformly decreasing. However, the loss is on a decreasing trend overall, which shows that the network successively learns the feature of the distribution as training progresses.
\medskip

After training, we denote the trained parameter in the network as $\theta_1$ and evaluate the network $f_{\theta_1}(x,t)$ at various time $t$ with a larger batch $\{x_i\}_{i=1}^{J'}$ of size $J'=1M$.
\begin{remark}
Though $J'\gg J=10K$, the computational cost of the network evaluation, $O(J')$, is obviously smaller than that of IPM, $O(J^2 \frac{t}{\Delta t})$ as mentioned. In addition, the {\bf generated} output of the network is the image of a continuous map ($f$) acting on finer samples of input distribution. This way, it is not a direct replica of training data when $t$, $A$, or  $\mu$ overlaps with a value in the training set. Then the generated output of the network will not congregate on points in the training empirical distribution, see Fig.\ref{fig:Kflow}. Leveraging this feature, we can further use these generations as a warm-up to accelerate IPM simulation, similar idea appears in \cite{DP_22} when computing invariant measures.
\end{remark}
In Fig.\ref{fig:eg1_training}(c), we plot the second moment of the particles by IPM solver (denoted as reference) and by our DP (denoted as network) generation. We see that the slope of second moment by reference IPM solver deviates from the network output after $t=0.075$. This is due to regularization $\delta >0$ in IPM. For training data, we rely on IPM solver most of the time when it has good accuracy prior to the time when the $\delta$ effect kicks in.
\begin{figure}[htbp]
    \centering
     \begin{subfigure}{0.3\textwidth}
    \includegraphics[width=\linewidth]{figures/reg_KS.png}
    \caption{Second moment from IPM with various regularization factors}
    \end{subfigure}
     \begin{subfigure}{0.3\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-traning.png}
    \caption{Training loss}
    \end{subfigure}
     \begin{subfigure}{0.3\textwidth}
   \includegraphics[width=\linewidth]{figures/2DTDC02ndmomv2.png}
    \caption{Second moment comparison}
    \end{subfigure}
    \caption{Performance of IPM (reference) and DeepParticle (network) algorithms.}
    \label{fig:eg1_training}
\end{figure}
In Fig.\ref{fig:eg1_samples}, we plot the histogram of the network output (1M particles) and reference distribution (10K particles) at different time $t$. We see that the network learns the feature of particle concentration and even has more concentrated output than the reference solver near the blowup time.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T1e-2.png}
    \caption{$t=0.01$}
    \end{subfigure}
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T4e-2.png}
    \caption{$t=0.04$}
    \end{subfigure}\\
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T10e-2.png}
    \caption{$t=0.1$ (largest $t$ in training data)}
    \end{subfigure}
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T12e-2.png}
    \caption{$t=0.12$ (near blowup time)}
    \end{subfigure}
    \caption{Comparison of particle distributions from DeepParticle (network output) and IPM (reference) solver at different times and $A=0$.}
    \label{fig:eg1_samples}
\end{figure}
\subsection{KS Simulation and Generation in the Presence of 2D Laminar Flows}
Next we consider $\textbf{v}\not=0$, so in addition to chemotaxis gradient, the movement of organism is also driven by some given environmental fluid velocity field. In \cite{chemomix_yao}, blow-up behavior of KS model given various strength of environment velocity is investigated. Now we let
\begin{align}\label{2Dlaminer}
    \textbf{v}(x,y)= A \, \begin{pmatrix}\exp(-y^2)\\0\end{pmatrix},
\end{align}
which represents a laminar flow of amplitude $A$ traveling along $x$ direction with $y$-dependent speed. There are two parameters to learn in the model: the advection amplitude $A$ and the evolution time $T$.
\paragraph{Learning dependence on $A$} We start from learning amplitude dependence while fixing the evolution time. In this example, we used IPM to generate $J=10000$ samples of solution after $T=0.02$ with $A=10^{[0:0.2:2]}$.  The initial distribution of IPM for each $A$ is a uniform distribution on the unit ball. In Fig.\ref{fig:2DlaminerC} we see that the distribution turns to a $V$ shape as $A$ increases. This can be attributed to the stretching effect by the laminar flow \eqref{2Dlaminer}. Our network is shown to learn this feature. In addition, our network can also predict distribution when $A$ is slightly outside the range of training set, see the case when $A=150$ in Fig.\ref{fig:2DlaminerC}(d). More precisely, Fig.\ref{fig:2DlaminerC}(c) shows that most of outputs in the training set ($1\leq A \leq 100$) satisfy $x<3$ while at $A=150$, a reasonable proportion of particles is on the right side of $x=3$ and our network indeed predicts it.
\begin{figure}[htbp]
    \centering

    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig10.png}
    \caption{$A=10$}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig80.png}
    \caption{$A=80$ (interpolation)}
    \end{subfigure}\\
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig100.png}
    \caption{$A=100$ (largest $A$ in training set)}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig150.png}
    \caption{$A=150$ (prediction)}
    \end{subfigure}
    \caption{Learning particle aggregation at different $A$ values with $t=0.02$ fixed.}
    \label{fig:2DlaminerC}
\end{figure}
\paragraph{Learning dependence on $T$}
We turn to study the dependence on the evolution time at fixed $A=100$. To generate training data, we run IPM with $t\in[0,0.1]$ and $J=10000$ particles, and take snapshots of the empirical distribution at $t=0,0.01,0.02,\cdots,0.1$. In Fig.\ref{fig:2DlaminerT}, we compare our network generated output with IPM reference solution at various evolution times. We see that both IPM and DP methods reproduce the near blow-up behaviors consistent with  \cite{chemomix_yao}.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD2e-2.png}
    \caption{$t=0.02$}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD4e-2.png}
    \caption{$t=0.04$}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD10e-2.png}
    \caption{$t=0.1$ (largest $t$ in training set)}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD14e-2.png}
    \caption{$t=0.14$ (prediction)}
    \end{subfigure}
    \caption{Learning particle aggregation at different times with fixed flow amplitude $A=100$.}
    \label{fig:2DlaminerT}
\end{figure}
\subsection{KS Simulation and Generation in the Presence of 3D flows} In this subsection, we turns to KS model in three dimension. There are two kinds of flows under consideration, the first is 3D laminar flow which is a natural generalization of 2D laminar flow. The second is Kolmogorov flow which is a well-known example of chaotic flow \cite{CG95,KLX_21}.
\paragraph{3D Laminar Flow}
As the first 3D example, we consider advection field in the form:
\begin{align}
    \textbf{v}(x,y,z)=A\, \begin{pmatrix}\exp(-y^2-z^2)\\0\\0\end{pmatrix}. \label{Laminar}
\end{align}
It describes the organism travelling along the $x$-direction while its speed depends on the radial position of $y$ and $z$ variable.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.8\textwidth}
    \includegraphics[width=\textwidth]{figures/laminar3Dsig10.png}
    \caption{ $A=10$}
    \label{fig:my_label}
    \end{subfigure}
    \begin{subfigure}{0.8\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3Dsig100.png}
    \caption{$A=100$}
    \label{fig:my_label}
    \end{subfigure}
    \caption{Generated distributions in three cross sections at different $A$ values in a 3D laminar flow (\ref{Laminar}).}
    \label{fig:3Dlaminar}
\end{figure}
In Fig.\ref{fig:3Dlaminar} we show the histogram of the generated distribution of our deep learning algorithms with $A=10$ and $A=100$, which reproduces the distribution of corresponding IPM simulation. The configuration of learning $A$ dependence is the same as one in 2D cases, except the input and output is now in $3$-dimension. From the numerical experiments (not shown), there is no need to increase the width of our network. By comparing (a) and (b) in Fig.\ref{fig:3Dlaminar}, we see that the distribution becomes V shape in $xy$ projection and $xz$ projection as the amplitude $A$ increases. This is due to the laminar setting of the flow. In $yz$ projection, the distribution remains as radial.
\begin{figure}[htbp]
    \centering
     \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C10.png}
    \caption{$A=10$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C30.png}
    \caption{$A=30$}
    \label{fig:my_label}
    \end{subfigure}\\
      \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C100.png}
    \caption{$A=100$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C130.png}
    \caption{$A=130$}
    \label{fig:my_label}
    \end{subfigure}
    \caption{Network generated vs. reference distributions projected to $xy$ plane in a 3D laminar flow (\ref{Laminar}).  }
    \label{fig:3DlaminarC}
\end{figure}
In Fig.\ref{fig:3DlaminarC}, we show the $xy$ projection of the distribution with various $A$. It confirms that in addition to interpretation, our network is able to generalize the influence of amplitude on the distribution when $A$ is not far from the range in the training set.
 \paragraph{3D Kolmogorov Flow}
In the last example, we investigate the case when the organism travels and aggregates in chaotic streamlines given by Kolmogorov flow \cite{CG95,KLX_21}:

\begin{align}
    \textbf{v}(x,y,z)=A\, \begin{pmatrix}\sin(2\pi z)\\\sin(2\pi x)\\\sin(2\pi y)\end{pmatrix}. \label{Kflow}
\end{align}
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/Kflow_C100.png}
    \caption{Network output ($N=1M$) vs. training data ($N=10K$) in flow (\ref{Kflow}) at   $A=100$, $t=0.1$.}
        \label{fig:Kflow}
\end{figure}
\begin{figure}[htbp]
    \centering
     \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C10compare.png}
    \caption{$A=10$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C30compare.png}
    \caption{$A=30$}
    \label{fig:my_label}
    \end{subfigure}\\
      \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C100compare.png}
    \caption{$A=100$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C110compare.png}
    \caption{$A=110$}
    \label{fig:my_label}
    \end{subfigure}
    \caption{Network generated vs. reference distributions projected to $xy$ plane in a 3D chaotic flow (\ref{Kflow}).}
    \label{fig:KflowC}
\end{figure}
In Fig.\ref{fig:Kflow}, we compare the network generated and reference distributions when $A=100$ and $t=0.1$ and Fig.\ref{fig:KflowC} compares the distributions with various amplitude $A$ when projected to $xy$ plane. The distributions are in general a radial distribution and the radius of the distribution increases when $A$ increases (see also the second moment plot in Fig.\ref{fig:Kflow2nd}). The phenomenon differs from the one in laminar flow. This may result from the mixing mechanism of the chaotic flow that spreads and acts against chemotaxis aggregation.
\begin{figure}[h]
    \centering
    \includegraphics[width=0.3\linewidth]{figures/Kflow_2nd.png}
    \caption{The second moment of generated/reference distribution vs. $A$ of 3D chaotic flow (\ref{Kflow}).}
    \label{fig:Kflow2nd}
\end{figure}
\section{Conclusions and Future Works}

We proposed a regularized interacting particle method (IPM) to compute aggregation patterns and near singular solutions in multi-dimensional Keller Segel (KS) systems.  We then studied a DeepParticle (DP) method to learn and generate solutions for KS systems with dependence on physical parameters (e.g. the small diffusivity of chemo-attractant and the reciprocal of the flow amplitude in the advection-dominated regime) by minimizing the 2-Wasserstein distance between the source and target distributions. During the training stage, we seek a mapping in the form of a deep neural network from
source to target distributions and update network parameters based on a discretized 2-Wasserstein distance defined on finite distribution samples. Our method is general in the sense that we do not require target  distributions to be in closed form and the generation map to be invertible. Our method is fully data-driven and applicable to the fast generation of distributions for more general KS systems with physical parameter dependency. Our iterative divide-and-conquer algorithm reduces considerably the computational cost of finding the optimal transition matrix in the Wasserstein distance. We carried out numerical experiments to demonstrate the performance of our method for KS aggregation learning and generation
without and with laminar and chaotic advections.
In future work, we plan to study DP method to learn and generate pattern forming solutions of parabolic type KS systems ($\epsilon > 0$ in (\ref{KS1})) among other KS like (e.g. chemotaxis-haptotaxis) systems for modeling and predicting cancer cell evolution \cite{chaplain2006}.

\section*{CRediT authorship contribution statement}
\noindent
\textbf{Zhongjian Wang}: Conceptualization, Programming, Methodology, Writing-Original draft preparation and Editing. \textbf{Jack Xin}: Conceptualization, Methodology, Writing-Reviewing and Editing. \textbf{Zhiwen Zhang}: Conceptualization,  Writing-Reviewing and Editing.

\section*{Declaration of competing interest}
\noindent
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

\section*{Acknowledgements}
\noindent
The research of JX is partially supported by NSF grants DMS-1924548 and DMS-1952644. The research of ZZ is supported by Hong Kong RGC grant (Projects 17300318 and 17307921), National Natural Science Foundation of China  (Project 12171406), Seed Funding Programme for Basic Research (HKU) and Seed Funding for Strategic Interdisciplinary Research Scheme 2021/22 (HKU).
The authors would like to thank Prof. John Lowengrub at UC Irvine and Prof. Philip Maini
at Oxford University for helpful discussions of cell growth and cancer modeling.

\bibliographystyle{plain}

\section{Introduction}

Chemotaxis partial differential equations (PDEs) were introduced by Keller and Segel (KS \cite{keller1970initiation}) to describe the aggregation of the slime mold amoeba Dictyostelium discoideum due to an attractive chemical substance. Related random walk model by Patlak
was known earlier \cite{Patlak_53}, see \cite{Oth_97} for
an analysis of basic taxis behaviors (aggregation, blowup and collapse) based on reinforced random walks.
Recall a common form of KS model \cite{fatkullin2012study}:
\begin{align}
    \rho_{t} = \nabla \cdot (\mu \, \nabla \rho - \chi\, \rho\,  \nabla c ), \;\; \;\;
    \epsilon \, c_t  = \Delta\, c - k^2 \, c + \rho, \label{KS1}
    \end{align}
    where $\chi, \mu $ ($\epsilon, k$) are positive (non-negative) constants. The model is called elliptic if $\epsilon =0$ (when $c$ evolves rapidly to a local equilibrium), and parabolic if $\epsilon >0$.
The $\rho$ is the density of active particles (bacteria), and $c$ is the concentration of
chemo-attractant.
The bacteria diffuse with mobility $\mu$ and drift in the direction of $\nabla c$ with velocity $\chi \nabla c$ where $\chi$ is called chemo-sensitivity.
\medskip

In the simplest regime
 $(\epsilon,k)=\mathbf{0}$, the concentration equation

 becomes the Poisson equation
$  - \Delta c = \rho$. Subject to a suitable boundary condition of $c$, the classical integral representation
$c=- \mathcal{K}\ast \rho$ holds,
where $\mathcal{K}$ is the Green's function of the Laplacian, and $\ast$ denotes convolution.
The KS system
then reduces to a scalar non-local nonlinear advection-diffusion PDE governing the evolution of the density function $\rho$:
\begin{align}
\rho_t = \mu \, \Delta \, \rho + \chi \, \nabla \cdot \big(\rho \, \nabla( \mathcal{K}\ast \rho)\big). \label{KSequation1}
\end{align}

For modeling chemotaxis in a fluid environment such as ocean,  \cite{chemomix_12,chemomix_16,chemomix_yao,chemomix_21,chemomix_22} studied equation
(\ref{KSequation1}) with the advective Lie derivative $\rho$ on the left hand:
\begin{align}
\rho_t +\nabla \cdot(\rho \, \textbf{v})= \mu \, \Delta \, \rho + \chi \, \nabla \cdot \big(\rho \, \nabla( \mathcal{K}\ast \rho)\big). \label{KSequation2}
\end{align}
The mixing mechanism of the flow field $\mathbf{v}$ is known to slow down or smooth out blowup or aggregation in (\ref{KSequation1}), see analysis in  \cite{chemomix_12,chemomix_16,chemomix_21,chemomix_22} and references therein, and related convection induced smoothing in fluids
\cite{hou2006dynamic,hou2009stabconv}.
Eq. \eqref{KSequation2} is the macroscopic limit (McKean-Vlasov equation) of the interacting particle system below as $J \uparrow \infty$:
\begin{align}
    d X^{j} = -\,  {\chi M \over J} \nabla_{X^j} \,
    \sum_{i=1:J,\, i\not = j}\, \mathcal{K} (|X^j - X^i |) \, dt + \textbf{v}(X^j) dt
+ \sqrt{2\,\mu} \, d W^j, \;\; j=1,\cdots, J,   \label{IPS1}
\end{align}
where $M$ is the conserved total mass (integral of $\rho$), $W^j$'s are independent Brownian motions in
$R^d$.
\medskip

In this work,
$ \textbf{v}$ is a prescribed divergence
free vector field. We shall approximate
$\rho$ of equation (\ref{KSequation2}) numerically based on the associated interacting particle system in two and three spatial dimensions ($d=2,3$), and carry out a systematic deep learning study (a.k.a. DeepParticle \cite{DP_22}) on the $(\mu, \textbf{v})$ dependency of solutions.
As we are interested in studying  near singular KS solutions, the main challenge for training data collection is to approximate the fields $\rho(\textbf{x})$ and $c(\textbf{x})$ reliably as they intensify. Due to singular behavior of Green's function as particles come close to each other, our approach is to regularize $\mathcal{K}$ in  (\ref{IPS1}) for approximating  $\rho(\textbf{x})$ as particles aggregate,
in similar spirit to
the vortex blob method for fluids (\cite{krasny1991vortex} and references therein) and \cite{blob_aggregation_ck2016}.
\medskip

Deep learning tools have been applied broadly for scientific computing in recent years, such as solving PDEs and their inverse problems, see \cite{Osher_21,DP_22} and references therein.
DeepParticle \cite{DP_22} is based on a particle
method for solving a time-dependent physically parameterized PDE, whose
solution is approximated by the particle empirical measure (distribution).  A deep neural
network (DNN) with physical parameter dependence
learns the mapping from the initial particle distribution to
the particle distribution at time $T$ with training data over sampled
physical parameters provided by the particle method.
The trained DNN then generates approximate solutions at time $T$
for new physical parameters unseen in training.
DeepParticle has been successfully designed and trained for learning
and generating invariant measures of stochastic interacting particle
systems (at $T=\infty$) arising in reaction-diffusion front speeds in three dimensional chaotic flows \cite{DP_22}. In this paper, we further develop DeepParticle to learn and generate KS solutions exhibiting aggregation behavior
at a finite time $T$ before blow-up for a range of diffusivity
$\mu$ and advection amplitude values.
\medskip

The rest of the paper is organized as follows. In Section 2, we briefly review the blow up behavior in the KS model and the regularized  particle methods to solve the KS model. In Section 3, we present our DeepParticle method to learn the transport map from an input distribution to a target distribution. Moreover, we will discuss details of implementation of our method and how to learn the distributions in particle simulation of the KS model.  In Section 4, we show numerical results to demonstrate the performance of our method. Concluding remarks are made in section 5.

\medskip

\section{Regularized interacting particle method of KS equation}

\subsection{Blow up behavior in KS model}
We start from the simplest KS model without advection, namely \eqref{KSequation1} with $\mu=\chi=1$. This system has been extensively studied by many authors, see survey article \cite{perthame2004pde}. The conservation of mass holds:
\begin{align}
\frac{d}{ dt}\int_{\mathbb{R}^{2}} \rho(x,t) d x = 0,
\end{align}
which is also true for  Eq. \eqref{KSequation2} when the advection field $\textbf{v}$ is divergence free.  If we set the total mass $ M:=\int_{\mathbb{R}^{2}} \rho(x,0) d x$, the second moment has fixed time derivative, i.e.,
\begin{align}\label{second_order_moment}
    \frac{d}{ dt} \int_{\mathbb{R}^{2}}|x|^{2} \rho(x,t)  x=\frac{M}{2 \pi}(8 \pi-M),
\end{align}
where $8 \pi$ is called critical mass of the system. Accordingly, it is well-known that: 1) if $M>8\pi$, the system has no global smooth solutions; 2) if $M=8\pi$, the system has a global smooth solution, which blow up as $t\to\infty$; 3) if $M<8\pi$, the system has a global smooth solution.
\medskip

In \cite{chemomix_yao}, an extra advection term is introduced to the KS equation, which reflects the organism movement in some prescribed fluid flows. Then the second moment identity in Eq.\eqref{second_order_moment} and the subsequent blowup vs. global evolution results are no long valid. By comparison principle, \cite{chemomix_yao} shows  that if the total mass is smaller than critical mass, the system has global smooth solution with smooth initial data. In the cases with supercritical mass, there are only numerical experiments suggesting that the advection, if sufficiently large,  prevents the solutions from blowing up. Later, \cite{chemomix_22} shows that when the flow exerts a  `stretching' effect, the advection field $\textbf{v}$ indeed suppresses the growth or the concentration of chemo-attractant and hence the solution has global regularity and exists for all time.  Examples of stretching flows include hyperbolic flows where $\mathbf{v}(\mathbf{x})=\left( x_{1}, -\frac{1}{d-1}\mathbf{x}_{-}\right)$ and laminar flows where $\mathbf{v}(\mathbf{x})=\left(v\left(\mathbf{x}_{-}\right), \mathbf{0}_{-}\right)$, with $\mathbf{x}_{-}=\left(x_{2}, \ldots, x_{d}\right)$ and $v$ is periodic. However the case of chaotic flows or when the amplitude of advection $\textbf{v}$ is not sufficiently large, remains open.

\subsection{Regularized interacting particle methods}

The singular behavior of the governing PDE and the associated Green's function that cause trouble for particle methods is a well-known problem. The vortex blob method \cite{chorin1973discretization} provides a regularization approach to extend vortex sheet motion past the singularity formation time into the physically important roll-up regime, in which the points representing the vortex sheet are replaced by vortices of prescribed and fixed shape. Numerical calculations show regular motion for the centers of the blobs even after the time when a curvature singularity on a vortex sheet is formed. Later, a special form of the vortex blob method \cite{krasny1986desingularization} is used to calculate the roll-up of a periodic vortex sheet resulting from the classical Kelvin-Helmholtz instability. The singular solutions are closely related to the ill-posedness of vortex sheet motion \cite{caflisch1989singular} . Nonetheless,
as the regularization parameter approaches zero, the regularized vortex sheet solution converges to a weak solution of the Euler equation \cite{liu1995convergence}.
 As in vortex blob methods, we formulate a regularized interacting particle method (denoted as IPM from here on) to solve KS chemotaxis systems.
\medskip

 Let us represent the density function (solution of Eq.\eqref{KSequation2}) with empirical distributions of particle positions. In SDE \eqref{IPS1} on particle positions,
the chemo-attractant term $\,  {\chi M \over J} \nabla_{X^j} \,
    \sum_{i=1:J,\, i\not = j}\, \mathcal{K} (|X^j - X^i |) \, dt$ causes numerical instabilities when particles tend to concentrate. To overcome this difficulty, we replace the singular kernel in \eqref{IPS1} by a smoothed approximation $K_{\delta}$,
such that $K_{\delta}(z)\rightarrow K(z)$ as $\delta \rightarrow 0$, where $\delta>0$ is a regularization parameter. For example, we define
\begin{align}
    K_{\delta}(z) = K(z)\frac{|z|^2}{|z|^2+\delta^2}.
\end{align}
\medskip

Equipped with kernel $K_{\delta}$, we
obtain a system of regularized SDEs for the particles as follows:
\begin{align}
    d X^{j} = -\,  {\chi M \over J} \nabla_{X^j} \,
    \sum_{i=1:J,\, i\not = j}\, K_{\delta}(|X^j - X^i |) \, dt+ \textbf{v}(X^j) dt
+ \sqrt{2\,\mu} \, \, d W^j, \;\; j=1,2,\cdots, J,
    \label{IPS2}
\end{align}
where $J$ is the number of particles, $X^{j}\in R^d$ is the position of the $j$-th particle, where $d W^j$ are mutually independent $d-$dimensional Brownian motions. The convergence of
a random particle blob method, similar to (\ref{IPS2}) yet for KS without advection field $\textbf{v}$, is analyzed in
\cite{rpbm_2017}. For a stochastic particle method using heuristic collision and splitting rules to bypass the singular behavior of Green's function, see \cite{havskovec2009stochastic}.
\medskip

Representing PDE solutions by particles belongs to the Lagrangian framework.  The Lagrangian methods have several advantages: (1) easy to implement; (2) spatially mesh-free and self-adaptive; and (3)  computational costs scale linearly with the dimension of spatial variables in the underlying stochastic dynamical systems. If we discretize KS system  \eqref{KSequation2} on mesh grids with e.g. finite element \cite{carrillo2019hybrid} and spectral \cite{shen2020unconditionally} methods, the number of mesh grids depends exponentially on the spatial dimension. The key benefit of the Lagrangian framework in computing KS models
is its natural capability to follow the KS solution when a singular behavior is emerging. As we shall see, the stochastic particle method based on (\ref{IPS2}) reproduces the well-known aggregation behavior and captures the KS dynamics during the potential blow-up stage of evolution. This is another step forward in our program of computing high gradient solutions in the Lagrangian framework, which has shown encouraging results for a range of multi-dimensional singularly-perturbed  advection-diffusion PDEs.
We refer interested readers to our recent progress in developing Lagrangian methods to compute effective diffusivites in chaotic or random flows \cite{WangXinZhang:18,IPM_2020,SharpMMS_21,TDChaoticFlow_22} and KPP front speeds in chaotic flows \cite{IPM_2021}.
There are also deterministic particle methods (\cite{blob_diff_2019,blob_aggregation_ck2016} and references therein) for a class of KS and degenerate diffusion equations that fall in our DeepParticle framework (section 3).
\medskip

Though the IPM in our study here is mesh-free and self-adaptive for solving multi-dimensional KS chemotaxis systems, the computational costs remain high if we want to study the systems under a variation of parameters (evolution time $T$ and amplitude of advection $A$). Also, the numbers of particles $J$ cannot be too large as the chemo-attractant term has $O(J^2)$ complexity in each step. On the other hand, our numerical simulation of Eq.\eqref{IPS2} shows that the distribution at finite time $T$, starting from the same initial distribution, may have continuous dependence on the physical parameters. Therefore in the next section, we will introduce a machine learning algorithm to learn the feature of continuous dependency on parameters and generate approximate samples with $O(J)$ complexity.

\section{Deep particle method}
In this section we introduce a DeepParticle algorithm to learn the features of the transport map from a trivial (input) distribution to a target (output) distribution. The mapping error is measured by the 2-Wasserstein distance.

\subsection{Discrete Wasserstein distance}
\noindent	Given distributions $\mu$ and $\nu$ defined on metric spaces $X$ and $Y$, let us construct a transport map $f^0:X \to Y$ such that
		$f^{0}_{*}(\mu)= \nu$, where star denotes the push forward of the map.
For any function $f:X \to Y$, the $2$-Wasserstein distance between  $f_*(\mu)$ and $\nu$ is defined by:

\begin{align}\label{def:Wdistance}
	W_{2}(f_*(\mu), \nu):=\left(\inf _{\gamma \in \Gamma(f_*(\mu), \nu)} \int_{Y \times Y} dist(y', y)^{2} \mathrm{~d} \gamma(y', y)\right)^{1 / 2},
\end{align}
where $\Gamma(f_*(\mu), \nu)$ denotes the collection of all measures on $Y \times Y$ with marginals $f_*(\mu)$ and $\nu$ on the first and second factors respectively and $dist$ denotes the metric (distance) on $Y$. A
straightforward derivation yields:

\begin{align}\label{def:Wdistance1}
    W_{2}(f_*(\mu), \nu)=\left(\inf _{\gamma \in \Gamma(\mu, \nu)} \int_{X \times Y} dist(f(x), y)^{2} \mathrm{~d} \gamma(x, y)\right)^{1 /2},
\end{align}
where $\Gamma(\mu, \nu)$ denotes the collection of all measures on $X \times Y$ with marginals $\mu$ and $\nu$ on the first and second factors respectively and still $dist$ denotes the metric (distance) on $Y$.
We approximate $\mu$ and $\nu$ by empirical distribution functions:

$	\mu=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_i}$,
$	\nu=\frac{1}{N}\sum_{j=1}^{N}\delta_{y_j}$, where $N$ is the number of samples of distribution. Under the setting of learning distribution from interacting particle methods, we take $N<J$ sub-samples from terminal time position of system \eqref{IPS1} to represent the distribution. We will re-sample them every $1000$ steps of training. The preceding technique is usually referred as \emph{mini-batch} in the machine learning literature.
\medskip

	Any joint distribution in $\Gamma(\mu,\nu)$ can be represented by an $N\times N$ doubly stochastic matrix \cite{sinkhorn1964relationship}, denoted as transition matrix,  $\gamma=(\gamma_{ij})_{i,j}$ satisfying:

	\begin{align}
	\label{def:bistochasticity}
		\gamma_{ij} \geq 0; \quad \quad
	 \forall j, ~\sum_{i=1}^{N}\gamma_{ij}=1;
	 \quad \quad
	 \forall i, ~\sum_{j=1}^{N}\gamma_{ij}=1.
	\end{align}

Then \eqref{def:Wdistance1} becomes
\begin{align}\label{def:Wdistance_approx}
\hat{W}(f):=\left(\inf_{\gamma \in \Gamma^N}\sum_{i,j=1}^N\, dist(f(x_i),y_j)^2 \gamma_{ij}\right)^{1/2}.
\end{align}
$\hat{W}(f)$ in \eqref{def:Wdistance_approx} has a simple intuitive interpretation: given a $\gamma \in \Gamma(\mu, \nu)$ and any pair of locations $(x,y)$, the value of $\gamma(x,y)$ tells us what proportion of $f_*(\mu)$ mass at $f(x)$ should be transferred to $y$, in order to reconfigure $f_*(\mu)$ into $\nu$. Computing the effort of moving a unit of mass from $y'$ to $y$ by $dist(f(x),y)^{2}$ yields the interpretation of $\hat{W}$ as the minimal effort (optimal transportation \cite{villani2021topics}) to reconfigure $f_*(\mu)$ mass distribution into that of $\nu$.

\subsection{Training data and network configuration}
  Note that given any fixed set of $\{x_i\}_{i=1}^N\subset R^d$ and $\{y_j\}_{j=1}^N\subset R^d$ (training data), we have derived Eq.\eqref{def:Wdistance_approx} to be minimized by gradient descent. In addition, we aim to find a network that can represent the change of target distribution over some physical parameters. In such scenario, more than one set of data ($\{x_i\}$ and $\{y_j\}$ consists one set of data) should be assimilated. More precisely, let the total number of data set be denoted as $N_{dict}$. Then we have $N_{dict}$ pairs of i.i.d. samples of input and output distribution, denoted as $\{x_{i,r}\}$ and $\{y_{j,r}\}$ for $r=1\cdots N_{dict}$. Associated with the $r$-th data set ($\{x_{i,r}\}$ and $\{y_{j,r}\}$), we assume that there is a physical parameter $\eta_r\in R^p$. To represent this in the network, we encode $\eta_r$ to each data in the set, i.e., the input of the network is $\{(x_{i,r},\eta_r)\}_i \subset R^{d+p}$. This procedure is also called \textit{padding} in the literature.
 The output of network is then denoted as $f_\theta(x;\eta)$ where $\theta$'s are all trainable parameters of the network.
 \medskip

 Between the input ($l_0$) and output ($l_6$), there are 5 latent layers where each layer is $30$ in width. The adjacent layers $l_i$ and $l_{i+1}$ are fully connected, i.e.,
 \begin{align}
     l_{i+1}=\tanh(W_i(l_{i})+b_i) \quad i=0,\cdots 4,
 \end{align}
 where $W_i$ is (weight) matrix width of $l_{i+1}$ and $l_{i}$; $b_i$ is the (bias) vector with same dimension as $l_{i}$; $\tanh$ is the activation function. For the output layer, the formula is similar except we do not apply the activation function.  In case of $d=2$ and $p=1$ (e.g. the first numerical example in computing blowup behaviour of KS without advection, see Sec.\ref{sec:eg1}), there are in total $4966$ parameters (weight and bias) to update by gradient descent. In 3D cases, we find that the network performs well even without altering the width of latent layers.

 \begin{remark}
 Comparing with par-net approach of \cite{DP_22}, we found that  \textit{padding} achieves similar performance when learning distribution in particle simulation of KS model.
 \end{remark}

  In our computation, by equipping $\mathbb{R}^d$ with Euclidean metric, the training loss function is
\begin{align}\label{def:penalty-1}
	\hat{W}^2(f_\theta):=\sum_{r=1}^{n_{\eta}}\left(\inf_{\gamma_r \in \Gamma^N}\sum_{i,j=1}^N|f_\theta(x_{i,r};\eta_r)-y_{j,r}|^2\gamma_{ij,r}\right).
\end{align}

The goal of deep particle algorithm is then to find $(\theta,\{\gamma_r\})$ to minimize
\begin{align}\label{def:goal}
	P(\theta,\{\gamma_r\}):=\sum_{r=1}^{n_{\eta}}\left(|f_\theta(x_{i,r};\eta_r)-y_{j,r}|^2\gamma_{ij,r}\right).
\end{align}

\subsection{Iterative method in finding transition matrix $\gamma$}
Notice that, with fixed $\theta$, finding the transition matrix $\gamma$ to calculate discrete Wasserstein distance in \eqref{def:penalty-1} is linear programming with degree of freedom $N^2$. When the number of particles ($N$) becomes large, it is expensive to go by a conventional algorithm, e.g. interior point algorithm or simplex algorithm.  We now present a {\it mini-batch linear programming algorithm} to find the best $\gamma$ for each
 inner sum of \eqref{def:penalty-1}, while suppressing $k_r$ dependence in $f_\theta$.

\medskip

The problem \eqref{def:penalty-1} is a linear program on the bounded convex set $\Gamma^N$ of vector space of real $N\times N$ matrices. By Choquet's theorem, this problem admits solutions that are extremal points of $\Gamma^N$. The set of all doubly stochastic matrix $\Gamma^N$ is viewed as Birkhoff polytope. The Birkhoff–von Neumann theorem \cite{schrijver2003combinatorial} states that such polytope is the convex hull of all permutation matrices, i.e., those matrices such that $\gamma_{ij}=\delta_{j,\pi(i)}$ for some permutation $\pi$ of $\{1,...,N\}$, where $\delta_{jk}$ is the Kronecker symbol.
 \medskip

 Our algorithm is defined iteratively. We start from a permutation matrix, e.g., $\gamma_{ij}=\delta_{ij}$. In each iteration, we randomly select columns and corresponding rows such that the submatrix is a permutation matrix. Then the entries of submatrix consist of  a linear programming problem under the constraint that maintains column-wise and row-wise sums equal to one. To be precise, we randomly choose  $\{i_k\}_{k=1}^M$, ($M\ll  N$)  from $\{1,2, \cdots, N\}$ without replacement. Then $j_k$ is selected such that $\gamma_{i_k,j_k}=1$.  The cost function of the sub-problem is
 \begin{align}\label{def:lp-sub-cost}
	C(\gamma^*):=\sum_{k,l=1}^M|f_\theta(x_{i_k})-y_{j_l}|^2\gamma^*_{i_kj_l}
\end{align}
subject to
\begin{align}\label{def:lp-sub-cons}
\left\{	\begin{array}{l}
		\sum_{k=1}^M \gamma^*_{i_k,j_l}=1\quad\forall l=1,\cdots,M\\
			\sum_{l=1}^M \gamma^*_{i_k,j_l}=1\quad\forall \, k=1,\cdots,M\\
	\gamma^*_{i_kj_l}\geq 0 \quad \forall \, k,l = 1,\cdots,M.
\end{array}\right.
\end{align}
Then $\gamma^*$ is again a permutation matrix.  The linear programming sub-problem of much smaller size is solved by the interior point method  \cite{wright1997primal}. In addition, as the goal is to find a permutation matrix, we set the tolerance of interior point to be relatively large and project the resulting approximation to a permutation matrix. In our approach, we terminate the iteration of sub-problem until $\min_k(\max_l \gamma_{i_k,j_l}^*)>0.5$. This is to ensure that there is a unique large-value entry in each column. As a projection, we update $\gamma$ by,
\begin{align}
    \gamma_{ij}=\left\{ \begin{array}{l}
         1 \quad \text{if}\ j=\arg\max_l \gamma^*_{il}, \\
         0 \quad \text{otherwise.}
    \end{array}\right.
\end{align}

We observe that the global minimum of $\gamma$ in  \eqref{def:penalty-1} is also the solution of sub-problems \eqref{def:lp-sub-cost} with arbitrarily selected rows and columns, subject to the row and column partial sum values of the global minimum. The selection of rows can be one's own choice. In our approach, in each step after gradient descent, we randomly select rows to solve the optimization problem iteratively.
\medskip

Comparing with Wasserstein GAN proposed in \cite{arjovsky2017wasserstein}, \eqref{def:goal} is a Min-Min optimization. Both Adam gradient descent of $\theta$ and mimi-batch optimization of $\gamma$ are iteratively defined. We then alternatively update $\theta$ and $\gamma$ to seek a global minimum of \eqref{def:goal}.
\medskip

The cost of finding optimal $\gamma$ increases as $N$ increases, however the network itself is independent of $\gamma$. After training, our network acts as a sampler from  some target distribution $\nu$
without assumption of closed-form distribution of $\nu$. At this stage, the input data is no longer limited by  training data, an arbitrarily large amount of samples approximately obeying $\nu$ can be generated through $\mu$ (uniform distribution).

  \subsection{Full Training Algorithm}
   The full training process is outlined in Alg. \ref{alg:1},
   and carried out on a quad-core CPU desktop with an RTX2070 8GB GPU at UC Irvine. The training data is collected from the first order explicit IPM which  solves the regularized SDE system (\ref{IPS2}) by Euler's method in time. IPM is also the reference solver for evaluating DeepParticle output.

\begin{algorithm}[htbp]
	\SetAlgoLined
		Randomly initialize weight parameters  $\theta$ in network $f_\theta:\mathbb{R}^d\to \mathbb{R}^d$; \\
	\Repeat{given number of training mini-batches, $N_{dict}$}{
	\For{ physical parameter set $r\gets0$ \KwTo $n_{\eta}$ }{
	randomly select $\{x_{i,r}\}$, $\{y_{j,r}\}$, $i,j=1:N$ from i.i.d. samples of input and target distribution (generated by IPM) with respect to physical parameter $\eta_r$;\\
	$\gamma_{ij,r}=\delta_{i,j}$, i.e., initialize as permutation matrix;\\
	}
	\If{not the first training mini-batch}{
		\For{ physical parameter set $r\gets0$ \KwTo $n_{\eta}$ }{
		$\delta P_r=+\infty$\\
			\While{$|\delta P_r|_{fro}<tol$}{
		$P_r=\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}$;\\
		randomly choose $\{i_{k,r}\}_{k=1}^M$, from $\{1,2, \cdots, N\}$ without replacement; \\
		find $\{j_{k,r}\}_{k=1}^M$, such that $\gamma_{i_{k,r}j_{k,r},r}=1$\\
		solve the linear programming
		sub-problem \eqref{def:lp-sub-cost}-\eqref{def:lp-sub-cons}
		to get $\gamma_r^*$ another permutation matrix;\\
		update $\{\gamma_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$ with $\{\gamma^*_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$.\\
		$\delta P_r=\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}^*-P_r$
	}
	}
	}
	\Repeat{given steps for each training mini-batch}{
		$P=\sum_{r=1}^{N_r}\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}$;\\
		$\theta\leftarrow \theta - \delta_1\nabla_{\theta} P$, $\delta_1$ is the learning step size;\\
		\Repeat{given linear programming steps, $N_{LP}$}{
					\For{ physical parameter set $r\gets0$ \KwTo $n_{\eta}$ }{
		$P_r=\sum_{i,j=1}^N|f_\theta(x_{i,r},\eta_r)-y_{j,r}|^2\gamma_{ij,r}$;\\
		randomly choose $\{i_{k,r}\}_{k=1}^M$, from $\{1,2, \cdots, N\}$ without replacement; \\
		find $\{j_{k,r}\}_{k=1}^M$, such that $\gamma_{i_{k,r}j_{k,r},r}=1$\\
		solve the linear programming
		sub-problem \eqref{def:lp-sub-cost}-\eqref{def:lp-sub-cons}
		to get $\gamma_r^*$ another permutation matrix;\\
		update $\{\gamma_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$ with $\{\gamma^*_{i_{k,r}j_{l,r}}\}_{k,l=1}^M$.\\
	}
	}
	}
	}
	$\mathbf{Return}$
	\caption{DeepParticle Learning}\label{alg:1}
\end{algorithm}
\section{Numerical Examples}
\subsection{2D KS Simulation and Generation in the Absence of Advection}\label{sec:eg1}
First we consider the KS model without advection, namely \eqref{KSequation1} with $\mu=\chi=1$.
A straightforward derivation shows if the initial mass $M>8\pi$ and has finite second moment, the system will blow up in finite time.
\medskip

As the first numerical example, we consider learning the change of distribution depending on evolution time $T$ starting from the initial distribution. The initial distribution is assumed to be a uniform distribution on a ball with radius $1$ centered at the origin. Assuming the total mass is $16 \pi$, we have then the initial second moment as $8\pi$, same as in \cite{chemomix_yao}. By Eq.\eqref{second_order_moment}, the system will blow-up when $t>0.125$. For training data, we applied IPM with regularization $\delta^2=1e-3$ for $T=0.1$ with $J=10000$ particles. We keep the snapshots of the empirical distribution every $0.05$ interval. During training, we consider a mini-batch of size $8\times 2000$, which means, we take $N_{dict}=8$ sets of training data at various time $t$ and in each set (mini-batch) we have $N=2000$ subsamples from $J=10000$ samples from IPM. Adams gradient descent is applied to learn the parameters (weight and bias) in the network. We renew the mini-batch every $1000$ steps and renew $\gamma$ every $200$ steps.
\medskip

Before applying the deep learning algorithm, we first investigate the error of regularized method. In Fig.\ref{fig:eg1_training}(a), we show the second moment of IPM simulation  with various regularization factors. We see that except for the value  $\delta^2=10^{-2}$, the smaller  regularization factor $\delta$ does not affect the accuracy of the IPM simulation when $t\in[0,0.1]$. In Fig.\ref{fig:eg1_training}(b), we show the training loss computed by Eq.\eqref{def:penalty-1}. Since we do not renew $\gamma$ at every steps of gradient descent of parameters, the training loss is not uniformly decreasing. However, the loss is on a decreasing trend overall, which shows that the network successively learns the feature of the distribution as training progresses. In Fig.\ref{fig:eg1_training}(c), we plot the second moment of the particles by IPM solver (denoted as reference) and by our DP (denoted as network) prediction. We see that the slope of second moment by reference IPM solver deviates from the network generated line after $t=0.075$. This is due to regularization $\delta >0$ in IPM. For training data, we rely on IPM solver most of the time when it has good accuracy prior to the time when $\delta$ effect kicks in.
\begin{figure}[htbp]
    \centering
     \begin{subfigure}{0.3\textwidth}
    \includegraphics[width=\linewidth]{figures/reg_KS.png}
    \caption{Second moment from IPM with various regularization factors}
    \end{subfigure}
     \begin{subfigure}{0.3\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-traning.png}
    \caption{Training loss}
    \end{subfigure}
     \begin{subfigure}{0.3\textwidth}
   \includegraphics[width=\linewidth]{figures/2DTDC02ndmomv2.png}
    \caption{Second moment comparison}
    \end{subfigure}
    \caption{Performance of IPM (reference) and DeepParticle (network) algorithms.}
    \label{fig:eg1_training}
\end{figure}
In Fig.\ref{fig:eg1_samples}, we plot the histogram of the output (1M particles) and reference distribution (10K particles) at different time $t$. We see that the network learns the feature of particle concentration and even has more concentrated output than the reference solver near the blowup time.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T1e-2.png}
    \caption{$t=0.01$}
    \end{subfigure}
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T4e-2.png}
    \caption{$t=0.04$}
    \end{subfigure}\\
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T10e-2.png}
    \caption{$t=0.1$ (largest $t$ in training data)}
    \end{subfigure}
    \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTDC0-T12e-2.png}
    \caption{$t=0.12$ (near blowup time)}
    \end{subfigure}
    \caption{Comparison of particle distributions from DeepParticle (prediction) and IPM (reference) solver at different times and $A=0$.}
    \label{fig:eg1_samples}
\end{figure}
\subsection{KS Simulation and Generation in the Presence of 2D Laminar Flows}
Next we consider $\textbf{v}\not=0$, so in addition to chemotaxis gradient, the movement of organism is also driven by some given environmental fluid velocity field. In \cite{chemomix_yao}, blow-up behavior of KS model given various strength of environment velocity is investigated. Now we let
\begin{align}\label{2Dlaminer}
    \textbf{v}(x,y)= A \, \begin{pmatrix}\exp(-y^2)\\0\end{pmatrix},
\end{align}
which represents a laminar flow of amplitude $A$ traveling along $x$ direction with $y$-dependent speed. There are two parameters to learn in the model: the advection amplitude $A$ and the evolution time $T$.
\paragraph{Learning dependence on $A$} We start from learning amplitude dependence while fixing the evolution time. In this example, we used IPM to generate $J=10000$ samples of solution after $T=0.02$ with $A=10^{[0:0.2:2]}$.  The initial distribution of IPM for each $A$ is a uniform distribution on the unit ball. In Fig.\ref{fig:2DlaminerC} we see that the distribution turns to a $V$ shape as $A$ increases. This can be attributed to the stretching effect by the laminar flow \eqref{2Dlaminer}. Our network is shown to learn this feature. In addition, our network can also predict distribution when $A$ is slightly outside the range of training set, see the case when $A=150$ in Fig.\ref{fig:2DlaminerC}(d). More precisely, Fig.\ref{fig:2DlaminerC}(c) shows that most of outputs in the training set ($1\leq A \leq 100$) satisfy $x<3$ while at $A=150$, a reasonable proportion of particles is on the right side of $x=3$ and our network indeed predicts it.
\begin{figure}[htbp]
    \centering

    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig10.png}
    \caption{$A=10$}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig80.png}
    \caption{$A=80$ (interpolation)}
    \end{subfigure}\\
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig100.png}
    \caption{$A=100$ (largest $A$ in training set)}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/laminar2Dsig150.png}
    \caption{$A=150$ (prediction)}
    \end{subfigure}
    \caption{Learning particle aggregation at different $A$ values with $t=0.02$ fixed.}
    \label{fig:2DlaminerC}
\end{figure}
\paragraph{Learning dependence on $T$}
We turn to study the dependence on the evolution time at fixed $A=100$. To generate training data, we run IPM with $t\in[0,0.1]$ and $J=10000$ particles, and take snapshots of the empirical distribution at $t=0,0.01,0.02,\cdots,0.1$. In Fig.\ref{fig:2DlaminerT}, we compare our network prediction with IPM generated reference solution at various evolution times. We see that both IPM and DP methods reproduce the near blow-up behaviors consistent with  \cite{chemomix_yao}.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD2e-2.png}
    \caption{$t=0.02$}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD4e-2.png}
    \caption{$t=0.04$}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD10e-2.png}
    \caption{$t=0.1$ (largest $t$ in training set)}
    \end{subfigure}
    \begin{subfigure}{0.48\textwidth}
    \includegraphics[width=\linewidth]{figures/2DTD14e-2.png}
    \caption{$t=0.14$ (prediction)}
    \end{subfigure}
    \caption{Learning particle aggregation at different times with fixed flow amplitude $A=100$.}
    \label{fig:2DlaminerT}
\end{figure}
\subsection{KS Simulation and Generation in the Presence of 3D flows} In this subsection, we turns to KS model in three dimension. There are two kinds of flows under consideration, the first is 3D laminar flow which is a natural generalization of 2D laminar flow. The second is Kolmogorov flow which is a well-known example of chaotic flow \cite{CG95,KLX_21}.
\paragraph{3D Laminar Flow}
As the first 3D example, we consider advection field in the form:
\begin{align}
    \textbf{v}(x,y,z)=A\, \begin{pmatrix}\exp(-y^2-z^2)\\0\\0\end{pmatrix}. \label{Laminar}
\end{align}
It describes the organism travelling along the $x$-direction while its speed depends on the radial position of $y$ and $z$ variable.
\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.8\textwidth}
    \includegraphics[width=\textwidth]{figures/laminar3Dsig10.png}
    \caption{ $A=10$}
    \label{fig:my_label}
    \end{subfigure}
    \begin{subfigure}{0.8\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3Dsig100.png}
    \caption{$A=100$}
    \label{fig:my_label}
    \end{subfigure}
    \caption{Generated outputs in three cross sections at different $A$ values in a 3D laminar flow (\ref{Laminar}).}
    \label{fig:3Dlaminar}
\end{figure}
In Fig.\ref{fig:3Dlaminar} we show the histogram of the prediction of our deep learning algorithms with $A=10$ and $A=100$, which reproduces the distribution of corresponding IPM simulation. The configuration of learning $A$ dependence is the same as one in 2D cases, except the input and output is now in $3$-dimension. From the numerical experiments (not shown), there is no need to increase the width of our network. By comparing (a) and (b) in Fig.\ref{fig:3Dlaminar}, we see that the distribution becomes V shape in $xy$ projection and $xz$ projection as the amplitude $A$ increases. This is due to the laminar setting of the flow. In $yz$ projection, the distribution remains as radial.
\begin{figure}[htbp]
    \centering
     \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C10.png}
    \caption{$A=10$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C30.png}
    \caption{$A=30$}
    \label{fig:my_label}
    \end{subfigure}\\
      \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C100.png}
    \caption{$A=100$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.45\textwidth}
        \includegraphics[width=\textwidth]{figures/laminar3D_C130.png}
    \caption{$A=130$}
    \label{fig:my_label}
    \end{subfigure}
    \caption{Generated vs. reference distributions projected to the $xy$ plane in a 3D laminar flow (\ref{Laminar}).  }
    \label{fig:3DlaminarC}
\end{figure}
In Fig.\ref{fig:3DlaminarC}, we show the $xy$ projection of the distribution with various $A$. It confirms that in addition to interpretation, our network is able to generalize the influence of amplitude on the distribution when $A$ is not far from the range in the training set.
 \paragraph{3D Kolmogorov Flow}
In the last example, we investigate the case when the organism travels and aggregates in chaotic streamlines given by Kolmogorov flow \cite{CG95,KLX_21}:

\begin{align}
    \textbf{v}(x,y,z)=A\, \begin{pmatrix}\sin(2\pi z)\\\sin(2\pi x)\\\sin(2\pi y)\end{pmatrix}. \label{Kflow}
\end{align}
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{figures/Kflow_C100.png}
    \caption{Network output ($N=1M$) vs. training data ($N=10K$) in flow (\ref{Kflow}) at   $A=100$, $t=0.1$.}
        \label{fig:Kflow}
\end{figure}
\begin{figure}[htbp]
    \centering
     \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C10compare.png}
    \caption{$A=10$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C30compare.png}
    \caption{$A=30$}
    \label{fig:my_label}
    \end{subfigure}\\
      \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C100compare.png}
    \caption{$A=100$}
    \label{fig:my_label}
    \end{subfigure}
      \begin{subfigure}{0.48\textwidth}
        \includegraphics[width=\textwidth]{figures/Kflow_C110compare.png}
    \caption{$A=110$}
    \label{fig:my_label}
    \end{subfigure}
    \caption{Generated vs. reference distributions projected to the $xy$ plane in a 3D chaotic flow (\ref{Kflow}).}
    \label{fig:KflowC}
\end{figure}
In Fig.\ref{fig:Kflow}, we compare the network generated and reference distributions when $A=100$ and $t=0.1$ and Fig.\ref{fig:KflowC} compares the distributions with various amplitude $A$ when projected to $xy$ plane. The distributions are in general a radial distribution and the radius of the distribution increases when $A$ increases (see also the second moment plot in Fig.\ref{fig:Kflow2nd}). The phenomenon differs from the one in laminar flow. This may result from the mixing mechanism of the chaotic flow that spreads and acts against chemotaxis aggregation.
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.35\linewidth]{figures/Kflow_2nd.png}
    \caption{The second moment of generated/reference distribution vs. $A$ of 3D chaotic flow (\ref{Kflow}).}
    \label{fig:Kflow2nd}
\end{figure}
\section{Conclusions and Future Works}

We proposed a regularized interacting particle method (IPM) to compute aggregation patterns and near singular solutions in multi-dimensional Keller Segel (KS) systems.  We then studied a DeepParticle (DP) method to learn and generate solutions for KS systems with dependence on physical parameters (e.g. the small diffusivity of chemo-attractant and the reciprocal of the flow amplitude in the advection-dominated regime) by minimizing the 2-Wasserstein distance between the source and target distributions. During the training stage, we seek a mapping in the form of a deep neural network from
source to target distributions and update network parameters based on a discretized 2-Wasserstein distance defined on finite distribution samples. Our method is general in the sense that we do not require target  distributions to be in closed form and the generation map to be invertible. Our method is fully data-driven and applicable to the fast generation of distributions for more general KS systems with physical parameter dependency. Our iterative divide-and-conquer algorithm reduces considerably the computational cost of finding the optimal transition matrix in the Wasserstein distance. We carried out numerical experiments to demonstrate the performance of our method for KS aggregation learning and generation
without and with laminar and chaotic advections.
In future work, we plan to study DP method to learn and generate pattern forming solutions of parabolic type KS systems ($\epsilon > 0$ in (\ref{KS1})) among other KS like (e.g. chemotaxis-haptotaxis) systems for modeling and predicting cancer cell evolution \cite{chaplain2006}.

\section*{CRediT authorship contribution statement}
\noindent
\textbf{Zhongjian Wang}: Conceptualization, Programming, Methodology, Writing-Original draft preparation and Editing. \textbf{Jack Xin}: Conceptualization, Methodology, Writing-Reviewing and Editing. \textbf{Zhiwen Zhang}: Conceptualization,  Writing-Reviewing and Editing.

\section*{Declaration of competing interest}
\noindent
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

\section*{Acknowledgements}
\noindent
The research of JX is partially supported by NSF grants DMS-1924548 and DMS-1952644. The research of ZZ is supported by Hong Kong RGC grant (Projects 17300318 and 17307921), National Natural Science Foundation of China  (Project 12171406), Seed Funding Programme for Basic Research (HKU) and Seed Funding for Strategic Interdisciplinary Research Scheme 2021/22 (HKU).
The authors would like to thank Prof. John Lowengrub at UC Irvine and Prof. Philip Maini
at Oxford University for helpful discussions of cell growth and cancer modeling.

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv


text:
Wrestling Selected to Host NCAA Division III Northeast Regional
WORCESTER -- For the second time in six years, WPI has been tabbed as a host of an NCAA Division III Wrestling Regional Championship. The Crimson and Gray will host the two-day Northeast Regional in the Sports & Recreation Center on Saturday, February 24th and Sunday February 25th.
WPI is one of six regional sites, three of which were announced recently. Others revealed to be hosting include Augsburg (MN) in the Upper Midwest Region and Loras (IA) in the Lower Midwest Region. Hosts for the remaining three regionals (Central, Southeast and Mideast) will be determined in the upcoming weeks.
With the announcement of the hosts, the NCAA also released new regional alignments. The Northeast Regional will featured 17 teams - all 16 members of the New England Wrestling Association (NEWA) along with newcomer NYU. The realignment sees Cortland State, Ithaca and Oneonta State, all teams that took part in the 2013 regional at WPI, move to the Mideast and join Oswego State, which changed regions in 2014.
All 14 of the NEWA schools that competed in 2013 are back along with new additions Castleton and New England College, both of which have added wrestling in the last five years.
"We are very pleased and excited to once again be one of the six NCAA regional hosts," stated WPI Head Wrestling Coach, Steve Hall. "WPI last hosted in 2013, which was the inaugural year of the regional format, and it is fitting that we are hosting again with the first year of the new regional alignment."
The top three wrestlers from each weight class in each of the six regional sites will advance to the national championships, allowing for a total of 180 qualifiers. The year's NCAA Division III Championship will be hosted by Baldwin Wallace in Cleveland, Ohio on Friday, March 9th and Saturday, March 10th.
WPI kicks off its season on Saturday November 11th at the Ursinus (PA) Fall Brawl.
Regional Alignment Info

redpajama_set_name: RedPajamaCommonCrawl


text:
Mayan Debut "War on Terror" Music Video
Mayan have just debuted a new music video for their song "War On Terror." You can check out the video below.
The band -- which features Epica guitarist/songwriter Mark Jansen, Jack Driessen (ex-After Forever) and Frank Schiphorst -- released their debut album, Quarterpast, last month on Nuclear Blast Records.
Quarterpast features guest vocals from several special guests, including Floor Jansen of Revamp, Simone Simons of Eipca), Henning Basse of Sons of Seasons, and Italian opera singer Laura Macrì.


text: Q: USB not recognized in Virtualbox - Ubuntu 12.04 Host - Windows 7 Guest I have Ubuntu 12.04 host and cannot detect usb devices on windows 7 guest on virtualbox (guest-addon is installed), though  I configured vboxusers group and enabled it for my user account and still cannot detect any device :(

redpajama_set_name: RedPajamaStackExchange


text:
Q: How can I slow down a video (which I use as wallpaper) so it will look more realistic? I have an .mpg file of the Earth spinning. I use it as a background/wallpaper for my desktop.  However, it spins too fast.
How can I slow it down?

A: You can use Windows Live Movie Maker.  Import the video into the program, then go to the "Video Tools" tab.  There will be a box where you can adjust the speed of the video.  Then you can just save the movie and you'll be done.


Morgan Company Lectern. Height 44", Width 20", Depth 15". Weight 35 lbs.
The solid oak lectern is perfect for modern or traditional spaces. The simple design allows the natural beauty of the wood to take center stage. The design has straight sides and a recessed front panel with 2 petite feet at the base of the front. The top is supported by the 3 sides while the interior is left open for storage. Two shelves are included to hold books and other materials. The lectern is available as shown or with a solid cast brass cross and brass lamp.
The simple oak lectern is crafted by skilled wood workers and is carefully sanded to a satin finish and is stained with Enduracote III stain, which is catalyzed to lock in color and resist water. The stain offers long-lasting durability and years of beauty. Please note that due to variation in wood grain and tone, actual finished pieces may vary from what is shown.