MFC with FEM]Generalized optimal transport and mean field control problems for reaction-diffusion systems with high-order finite element computation
Fu]Guosheng Fu
Department of Applied and Computational Mathematics and
Statistics, University of Notre Dame, Notre Dame, IN, USA.
G. Fu's work is supported by NSF grant DMS-2012031.

Osher]Stanley Osher
Department of Mathematics, University of California, Los Angeles,
Los Angeles, CA, USA.
S. Osher's work is supported in part by AFOSR MURI FP 9550-18-1-502 and ONR grants:  N00014-20-1-2093 and N00014-20-1-2787.

Pazner]Will Pazner
Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR, USA.

Li]Wuchen Li
Department of Mathematics, University of South Carolina, Columbia, SC, USA.
W. Li's work is supported by AFOSR MURI FP 9550-18-1-502, AFOSR YIP award No. FA9550-23-1-0087, NSF DMS-2245097, and NSF RTG: 2038080.

35K57, 49N80, 49M41, 65N30.

We design and compute a class of optimal control problems for reaction-diffusion systems. They form mean field control problems related to multi-density reaction-diffusion systems. To solve proposed optimal control problems numerically, we first apply high-order finite element methods to discretize the space-time domain and then solve the optimal control problem using augmented Lagrangian methods (ALG2). Numerical examples, including generalized optimal transport and mean field control problems between Gaussian distributions and image densities, demonstrate the effectiveness of the proposed modeling and computational methods for mean field control problems involving reaction-diffusion equations/systems.

§ INTRODUCTION

Reaction-diffusion systems are essential classes of modeling dynamics [37], which have applications in tumor growth modeling [9], propagation of pandemic spread [24, 25], evolutionary games [19], etc. The reaction term reflects general nonlinear interacting behaviors of agents/particles in complex systems. Patterns of population behaviors often arise in the solution of reaction-diffusion equations, which model the collective behaviors of particles/agents.

In recent years, mean field control problems [20, 22] have been studied, which are optimal control problems for mean-field limits of infinitely many identical particles/agents. The problem models the complex behaviors of identical particles/agents interacting with each other. It is worth mentioning that a special example of a mean-field control problem forms the dynamical optimal transport problem [3]. It studies a particular optimal control problem in density space, which transfers an initial density toward the terminal density function. The optimal transport problem also introduces a functional distance, namely the Wasserstein distance. It helps study and compute a class of initial value evolutionary dynamics, namely, Wasserstein gradient flows [21]. They are also valuable for modeling interaction behaviors of particle dynamics.

Nowadays, optimal transport and mean field control problems have vast applications in modeling and computations. Classical studies of optimal transport and mean field control problems are often limited to a single-density function. The complex interaction between densities is a vital modeling factor that has not been systemically studied.

In this paper, we generalize a class of mean field control of reaction-diffusion equations/systems proposed in [27]. It forms an optimal control problem among multiple density functions interacting with each other with a nonlinear reaction vector function. After change of variables, we reformulate the optimal control problem into an optimization problem with a Fisher information type potential functional. We apply the high-order finite element method to discretize the spatial-time domain and use the augmented Lagrangian method, ALG2 from [15], to compute the discretized optimization problems. Numerical examples, including mean field control problems between Gaussian densities and a system with $12$ images, demonstrate the solution of the proposed generalized mean field control problems.

Generalized optimal transport and mean field control problems have been widely investigated in  [6, 7, 14, 35].
For example, multi-population mean field games were discussed in [4], and generalized optimal transport distances between vector densities were formulated in [11]. Meanwhile,  unnormalized and unbalanced optimal transport were proposed in [12, 23, 32], which allows to control densities with different total masses. In applications, one also applies mean field control problems to model the propagation of pandemics [24, 25].
They are all important examples of optimal control problems for reaction-diffusion systems. Modeling and computational mean-field control problems for general convection-reaction-diffusion systems are new research directions [28, 29]. They have potential applications in classical modeling dynamics based on reaction-diffusion equations and systems. This paper proposes a class of mean field control of reaction-diffusion systems. Meanwhile, the optimal control problem of gradient flows in Wasserstein space is widely studied, namely Schrödinger bridge problems [10, 13, 26, 31, 36]. In addition, the generalized gradient flows in multiple-density spaces have been studied [17, 18, 34, 35]. Our formulation extends the optimal control problem of gradient flows based on reaction-diffusion equations and systems. Thus, we study the mean field control using formalisms of controlling gradient flows. In simulations, we apply high-order finite element schemes to simulate the proposed generalized optimal transport and mean field control problems.

The paper is organized as follows. In section <ref>, we first review gradient flows and their induced metric distances in generalized optimal transport spaces. We then formulate and derive optimality conditions for generalized optimal transport and mean field control problems of reaction-diffusion systems. In section <ref>, we approximate the proposed mean field control problems using high-order finite element schemes and then use the ALG2 method to compute the discretized optimization problems. Numerical results, including scalar reaction-diffusion equations and two-species, twelve-species reaction-diffusion systems, are presented in section <ref>. Finally, we conclude with a discussion in section <ref>.

§ GENERALIZED OPTIMAL TRANSPORT AND MEAN FIELD CONTROL OF REACTION-DIFFUSION SYSTEMS

This section presents the main formulation of mean field control (MFC) problems for scalar reaction-diffusion equations and systems.

§.§ Background: Reaction-diffusion induced  metric distances

Before delving into the mean field control problems that will be discussed in this manuscript, we first review the definitions of two metric distances: one for scalars and one for systems, which are obtained from generalized optimal transport type gradient flow problems. The material in this subsection follows closely our previous work on variational time implicit schemes for reaction-diffusion systems in [17].
We only give the definition of these metric distances without further elaboration, but refer to the references [1, 26, 27, 34] for more details on
optimal transport type gradient flows, distances, mean-field control and related problems.

§.§.§ The metric distance: scalar case

The scalar metric distance is derived from the following reaction-diffusion equation [27, 17]:

\begin{equation}\label{dissipative}
\begin{split}
&\partial_t \rho=\nabla\cdot(V_1(\rho)\nabla\frac{\delta}{\delta \rho}\mathcal{E}(\rho))-V_2(\rho)\frac{\delta}{\delta\rho}\mathcal{E}(\rho), \text{ on }[0,T]\times\Omega,
\end{split}
\end{equation}

with homogeneous Neumann boundary conditions
$V_1(\rho)\nabla\frac{\delta}{\delta \rho}\mathcal{E}(\rho)\cdot\nnu|_{\partial \Omega} = 0$, where $\nnu$ is the unit outward normal direction on the boundary $\partial\Omega$. Here $\Omega\subset\mathbb{R}^d$ is the spatial domain, $\rho\colon [0,T]\times\Omega \rightarrow \mathbb{R}$ is a scalar non-negative density function satisfying

\begin{equation}
\label{d-space0}
\rho(t,\cdot) \in  \mathcal{M}=\{\rho\in H^1(\Omega): \rho\geq 0\},\quad \forall t\geq 0,
\end{equation}

$\mathcal{E}\colon \mathcal{M}\rightarrow\mathbb{R}$ is an energy functional, $V_1, V_2\colon \mathbb{R}_+\rightarrow\mathbb{R}_+$ are two positive mobility functions, and $\frac{\delta}{\delta\rho}$ is the  first variation operator in $L^2$ space.

A crucial property of the equation (<ref>) is that it
satisfies an energy-dissipation law:

\begin{equation}
\label{ener-law}
\frac{d}{dt}\mathcal{E}(\rho)=-\int_\Omega\Big[\|\nabla\tfrac{\delta}{\delta \rho}\mathcal{E}(\rho)\|^2V_1(\rho) + |\tfrac{\delta}{\delta \rho}\mathcal{E}(\rho)|^2V_2(\rho)\Big]dx\leq  0,
\end{equation}

where we use the fact that functions $V_1$, $V_2$ are non-negative.
In the literature, the above right hand side (dissipation rate)

\begin{align}
\label{fisher}
\mathcal{I}(\rho)=
\int_\Omega\Big[\|\nabla\tfrac{\delta}{\delta \rho}\mathcal{E}(\rho)\|^2V_1(\rho) + |\tfrac{\delta}{\delta \rho}\mathcal{E}(\rho)|^2V_2(\rho)\Big]dx,
\end{align}

is called the generalized Fisher information functional.
This Fisher information functional induces a metric in the space $\mathcal{M}$, which further defines a distance between two densities $\rho^0, \rho^1\in \mathcal{M}$.

Define a distance functional
\[\mathrm{Dist}_{V_1,V_2}\colon \mathcal{M}\times \mathcal{M}\rightarrow\mathbb{R}_+\]
as below. Consider the following optimal control problem:

\begin{equation}\label{Dis1}
\mathrm{Dist}_{V_1, V_2}(\rho^0, \rho^1)^2:=\inf_{\rho, \bm v_1, v_2}\quad\int_0^T\int_\Omega \Big[\|\bm v_1\|^2 V_1(\rho)+ |v_2|^2V_2(\rho)\Big]dxdt,
\end{equation}

where the infimum is taken among $\rho(t,x)\colon [0,T]\times\Omega\rightarrow\mathbb{R}_+$,  $\bm v_1(t,x)\colon  [0,T]\times \Omega \rightarrow\mathbb{R}^d$,
$v_2(t,x)\colon  [0,T]\times \Omega \rightarrow\mathbb{R}$,
such that $\rho$ satisfies a reaction-diffusion type equation with drift vector field $\bm v_1$, drift mobility $V_1$, reaction rate $v_2$, reaction mobility $V_2$, connecting initial and terminal densities $\rho^0$, $\rho^1\in\mathcal{M}$:

\begin{equation}\label{Dis2}
\left\{\begin{aligned}
&\partial_t \rho + \nabla\cdot( V_1(\rho) \bm v_1)=V_2(\rho)v_2,\quad \text{on } [0,T]\times \Omega,\\
&\rho(0, x)=\rho^0(x),\quad \rho(T,x)=\rho^1(x),
\end{aligned}\right.
\end{equation}

with no-flux boundary condition $V_1(\rho)\bm v_1\cdot\nnu|_{\partial\Omega}=0$

Introducing the flux function $\bmm(t,x)\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$ and source function $s(t, x)\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}    \bmm(t,x)=V_1(\rho(t,x))\bm v_1(t,x), \quad s(t,x)=V_2(\rho(t,x))v_2(t,x),
\end{equation*}

the distance in Definition <ref> is rewritten as
the following  optimization problem with a linear constraint:

\begin{equation}
\label{DisM1}
\mathrm{Dist}_{V_1, V_2}(\rho^0, \rho^1)^2
:=\inf_{\rho, \bmm, s}\quad\int_0^T\int_\Omega \Big[\frac{\|\bmm\|^2} {V_1(\rho)}+ \frac{|s|^2}{V_2(\rho)}\Big]dxdt,
\end{equation}

such that

\begin{equation}
\label{DisM2}
\begin{split}
&\partial_t \rho(t,x)+ \nabla\cdot \bmm(t,x) =s(t,x), \quad \text{in }[0,T]\times\Omega,\\
\bmm\cdot\nnu = 0, \quad\text{ on } [0,T]\times\partial\Omega,\\
\rho(0,x)=\rho^0(x),\; \rho(T,x)=\rho^1(x), \quad
\text{in }\Omega.
\end{split}
\end{equation}

The optimization functional in (<ref>)
is convex under the condition that both
$V_1$ and $V_2$ are positive concave functions.

§.§.§ The metric distance: system case

The system  metric distance is derived from the following reaction-diffusion system
$M$ species and $R$ reactions [17, 34]:

\begin{equation}\label{dissipativeS}
\begin{split}
\partial_t\rho_i =\;\nabla\cdot\left(
V_{1,i}(\rho_i)\nabla\frac{\delta}{\delta \rho}\mathcal{E}_i(\rho_i)\right) \; -\sum_{p=1}^RV_{2,p}(\bmr)\gamma_{i,p}\sum_{j=1}^M\gamma_{j,p}
\frac{\delta}{\delta \rho}\mathcal{E}_j(\rho_j),
\end{split}
\end{equation}

with homogeneous Neumann boundary conditions
$V_{1,i}(\rho_i)\nabla\frac{\delta}{\delta \rho}\mathcal{E}_i(\rho_i)\cdot\nnu|_{\partial\Omega}=0$
for $1\le i\le M$.
Here $\rho_i$ is the density such that
$\mathcal{E}_i\colon \mathcal{M}\rightarrow\mathbb{R}$ is the energy functional
and $V_{1,i}\colon \mathbb{R}_+\rightarrow\mathbb{R}_+$ is the positive mobility function
for $i$-th species, and $V_{2,p}\colon \mathbb{R}^M_+\rightarrow\mathbb{R}_+$ is the positive mobility function for $p$-th reaction.
The bolded density $\bm\rho=(\rho_1,\cdots, \rho_M)$ is simply the collection of $M$ densities.
Moreover, the coefficient matrix  $\Gamma=(\gamma_{i,p})\in\mathbb{R}^{M\times R}$ satisfies
$\sum_{i=1}^M\gamma_{i,p}=0$ for all $1\le p\le R$, which is related to the total mass conservation:
$\frac{d}{dt}\int_{\Omega}\sum_{i=1}^M\rho_i\,dx =0.

The equation (<ref>)
satisfies the energy-dissipation law:

\begin{equation}
\label{ener-lawS}
\frac{d}{dt}\sum_{i=1}^M\mathcal{E}_i(\rho_i)=-
\bm{\mathcal{I}}(\bm\rho)
\le 0,
\end{equation}

where $\bm{\mathcal{I}}(\bm\rho)$ denotes the generalized Fisher information functional

\begin{equation}
\label{fisherS}
\bm{\mathcal{I}}(\bm\rho):=
\int_\Omega\Big[\sum_{i=1}^M\|\nabla\frac{\delta\mathcal{E}_i}{\delta \rho}\|^2V_{1,i}(\rho_i) +
\sum_{p=1}^R\left|\sum_{j=1}^M\gamma_{j,p}\frac{\delta\mathcal{E}_j}{\delta \rho}\right|^2V_{2,p}(\bm\rho)\Big]dx.
\end{equation}

In the above, we use the fact that functions $V_{1,i}$, $V_{2,p}$, $1\leq i\leq M$, $1\leq p\leq R$, are non-negative. The functional $\bm{\mathcal{I}}$, again, induces a metric function in space $\mathcal{M}^M$, which defines distances between two system densities $\bm\rho^0, \bm\rho^1\in \mathcal{M}^M$.

Define a distance functional

\begin{align*}
\bm{\mathrm{Dist}}_{\bm V_1,\bm V_2}\colon \mathcal{M}^M\times \mathcal{M}^M\rightarrow\mathbb{R}_+
\end{align*}

as below.

\begin{equation}
\label{DisS}
\begin{aligned}
\mathrm{\mathbf{Dist}}_{V_1, V_2}(\bmr^0, \bmr^1)^2:=& \inf_{\bmr, \cmm, \bms} \Big\{\int_0^T\int_\Omega
\left(\sum_{i=1}^M\frac{|\bmm_i|^2}{V_{1,i}(\rho_i)}+
\sum_{p=1}^R\frac{|s_p|^2}{V_{2,p}(\bmr)}\right)
$\partial_t\rho_i+\nabla\cdot \bmm_i=
\sum_{p=1}^R\gamma_{i,p}s_p,\forall 1\le i\le M$,\\
$\bmr(0,\cdot)=\bmr^0, \;\;
\bmr(T,\cdot)=\bmr^1.$
\end{tabular}\Big\},
\end{aligned}
\end{equation}

where $\cmm=(\bmm_1,\cdots, \bmm_M)$ is the collection of fluxes, and
$\bms=(s_1,\cdots, s_R)$ is the collection of sources.

Again, the optimization functional in (<ref>)
is convex under the condition that
$V_{1,i}$ and $V_{2,p}$ are positive concave functions for all $1\le i\le M$ and $1\le p\le R$.

§.§ MFC for scalar reaction-diffusion

In this subsection, we introduce the following MFC problem. This can be viewed as a generalization of the scalar distance (<ref>).

Suppose that there exist two positive mobility functions $V_1,V_2:\mathbb{R}_+\rightarrow\mathbb{R}_+$,
a terminal time $T>0$, a non-negative regularization parameter $\beta\ge0$,
an energy functional $\mathcal{E}\colon \mathcal{M}\rightarrow\mathbb{R}$,
a potential functional $\mathcal{F}\colon \mathcal{M}\rightarrow\mathbb{R}$,
and a terminal functional
$\mathcal{G}\colon \mathcal{M}\rightarrow\mathbb{R}$.

\begin{equation}\label{smfcA1}
\begin{split}
&\inf_{\rho, \tilde{\bmm}, \tilde s}\quad\int_0^T\left[\int_\Omega \left(\frac{\|\tilde{\bmm}\|^2} {2V_1(\rho)}+ \frac{|\tilde s|^2}{2V_2(\rho)}\right)dx
\end{split}
\end{equation}

where the infimum is among all densities
$\rho$ with $\rho(t,\cdot)\in\mathcal{M}$ for $t\in[0, T]$, flux $\tilde{\bmm}\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and source $\tilde s\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, such that

\begin{equation}\label{smfcA2}
\partial_t \rho + \nabla\cdot \tilde{\bmm}-\tilde s=\beta
\left(\nabla\cdot(V_1(\rho)\nabla\frac{\delta}{\delta \rho}\mathcal{E}(\rho))-V_2(\rho)\frac{\delta}{\delta\rho}\mathcal{E}(\rho)\right),
\end{equation}

with  boundary condition

\begin{align}
\label{smfcA3}
\left.\left(\tilde{\bmm}-\beta V_1(\rho)\nabla\frac{\delta}{\delta \rho}\mathcal{E}(\rho)\right)\cdot\nnu\right|_{\partial\Omega} =0,
\end{align}

and fixed initial density $\rho(0,\cdot) = \rho^0$ in $\Omega$.

From the modeling perspective,  the variational problem in Definition <ref> models the movement of the density from the initial density under the dynamical constraint of reaction-diffusion equations. It aims to find the optimal choices of vector fields and reaction rate functions under general kinetic, potential, and terminal energies.

We will construct numerical scheme for problem (<ref>) based on the following change of variables.

Denote $\bmm\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
\bmm=\tilde\bmm-\beta V_1(\rho)\nabla\frac{\delta}{\delta \rho}\mathcal{E}(\rho),\quad s=\tilde s-\beta V_2(\rho)\frac{\delta}{\delta \rho}\mathcal{E}(\rho),
\end{equation*}

and denote

\begin{equation}
\label{v3-form}
V_3(\rho)=\frac{1}{V_1(\rho)|\frac{\delta^2}{\delta \rho^2}\mathcal{E}(\rho)|^2},
\end{equation}

where $\frac{\delta^2}{\delta \rho^2}\mathcal{E}(\rho)$ is the second variational derivative of the energy functional $\mathcal{E}(\rho)$.
Then, the scalar MFC problem (<ref>) is equivalent to the following optimization problem:

\begin{equation}\label{smfc2A1}
\begin{split}
&\inf_{\rho,\bmm, s} ~~ \int_0^T\int_\Omega \Big[\frac{\|\bmm\|^2}{2 V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\beta^2\frac{|\nabla \rho|^2}{2V_3(\rho)}+\frac{\beta^2}{2}{|\frac{\delta\mathcal{E}}{\delta \rho}|^2}{V_2(\rho)}\Big]dxdt\\
\mathcal{E}(\rho^0)],
\end{split}
\end{equation}

where the infimum is among all densities $\rho$ with $\rho(t,\cdot)\in\mathcal{M}$, for $t\in [0, T]$, flux $\bmm\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and source $s\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, such that

\begin{equation}\label{smfc2A2}
\partial_t \rho(t, x) + \nabla\cdot \bmm(t,x)-s(t,x)=0,
\end{equation}

with no-flux boundary condition

\bmm\cdot\nnu|_{\partial\Omega}=0,

and fixed initial density $\rho(0,\cdot)=\rho^0$.

The proof of the above proposition follows from a similar argument in <cit.>.  We present the detailed derivation in the Appendix. The main idea is to use integration by parts to show that

\begin{equation*}
\begin{split}
&\int_0^T\int_\Omega \Big[\frac{\|\tilde{\bmm}\|^2} {2V_1(\rho)}+\frac{\|\tilde{s}\|^2} {2V_2(\rho)}\Big]dxdt\\
\int_0^T\int_\Omega \Big[\frac{\|{\bmm}\|^2} {2V_1(\rho)}
+\frac{\|s\|^2} {2V_2(\rho)}
+\beta^2\frac{|\nabla \rho|^2}{2V_3(\rho)}
+\frac{\beta^2}{2}{|\frac{\delta\mathcal{E}}{\delta \rho}|^2}{V_2(\rho)}
\Big]dxdt\\
\mathcal{E}(\rho^0)].
\end{split}
\end{equation*}

For simplicity of discussions, we assume that there are two functions $F:[0,T]\times\Omega\times \mathbb{R}_+\rightarrow \mathbb{R}$
and $G:\Omega\times\mathbb{R}_+\rightarrow \mathbb{R}$, such that

\begin{align*}
\int_{\Omega}F(t, x, \rho)dx = &\;
\mathcal{F}(\rho(t,\cdot))-\int_{\Omega}\frac{\beta^2}{2}{|\frac{\delta\mathcal{E}}{\delta \rho}(\rho(t,x))|^2}{V_2(\rho(t,x))}dx,\quad\forall \,t\geq 0,\\
\int_{\Omega}G(x,\rho)dx = &\;
\mathcal{G}(\rho)+\beta\mathcal{E}(\rho).
\end{align*}

Then the functional in (<ref>)
simplifies to the following formula:

\begin{align*}
\int_0^T\int_\Omega \Big[\frac{\|\bmm\|^2}{2 V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2|\nabla \rho|^2}{2V_3(\rho)}-F(t,x,\rho)\Big]dxdt+
\int_{\Omega}
G(x, \rho(T,x))dx.
\end{align*}

This is the form we use in our numerical scheme.

For completeness, we formulate the Karush–Kuhn–Tucker (KKT) condition, i.e., the critical point system, for the optimization problem (<ref>). The proof is presented in the Appendix; see also [27, 30].

Let $(\rho, \bmm, s)$ be the critical point of the optimization problem (<ref>). Then there exists a function $\phi\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that

\begin{equation*}
\frac{\bmm(t,x)}{V_1(\rho(t,x))}=\nabla \phi(t,x),\quad \frac{s(t,x)}{V_2(\rho(t,x))}=\phi(t,x),
\end{equation*}

\begin{equation}\label{SBPI}
\left\{\begin{aligned}
&\partial_t\rho(t,x) +\nabla\cdot(V_1(\rho(t,x))\nabla \phi(t,x))-V_2(\rho(t,x))\phi(t,x)=0,\\
&\partial_t\phi(t,x)+\frac{1}{2}\|\nabla \phi(t,x)\|^2 V_1'(\rho(t,x))+\frac{1}{2}|\phi(t,x)|^2V'_2(\rho(t,x))\\
&\hspace{6.3cm}+\frac{\delta}{\delta \rho}\Big[\mathcal{F}(\rho)-\frac{\beta^2}{2}\mathcal{I}(\rho)\Big]=0,
\end{aligned}\right.
\end{equation}

is the generalized Fisher information functional given in (<ref>), with initial and terminal time boundary conditions

\begin{equation}
\rho(0,x)=\rho^0(x),\quad \phi(T,x)=-\frac{\delta}{\delta \rho}\Big(\mathcal{G}(\rho(T,\cdot))+\beta \mathcal{E}(\rho(T,\cdot))\Big).
\end{equation}

§.§.§ Examples

Different choices of the functions
$V_1$, $V_2$, $V_3$, $F$, $G$,
along with the regularization parameter $\beta\ge0$ and initial density $\rho^0$,
lead to various scalar MFC problems (<ref>).
We present some examples of these functions in Table <ref> below.

$V_1(\rho)$:     $\alpha_1\rho^{\gamma_1}$ with $ 0\le \gamma_1\le 1$. $\alpha_1>0$.

$V_2(\rho)$:     $\alpha_2\rho^{\gamma_2}$ with $0\le \gamma_2\le 1$, or
$\alpha_2\frac{\rho-1}{\log(\rho)}$. $\alpha_2\ge 0$

$V_3(\rho)$:     $\alpha_3\rho^{\gamma_3}$ with $ 0\le \gamma_3\le 1$. $\alpha_3>0$.

$F(t,x,\rho)$:     $-\tau\rho(\log(\rho)+V(t,x))$, or
$-\tau(\rho^m+\rho V(t,x))$.
$\tau\ge 0$, $m> 1$.

$G(x,\rho)$:     $\gamma(\rho-\rho^1)^2$, or
$\gamma> 0$.

Example of various functions for scalar MFC (<ref>).

$V(t,x)$ for $F(t,x,\rho)$ is a given drift coefficient, and
$\rho^1$ for $G(x, \rho)$ is a given terminal density.
Taking the energy

\begin{align*}
\mathcal{E}(\rho)=\begin{cases}
\int_{\Omega}\frac{1}{\sqrt{\alpha_1\alpha_3}}\rho\log(\rho)\,dx, \quad \hspace{1.4cm}\text{ if }r_1=r_3=1,\\[1ex]
\int_{\Omega}\frac{4\rho^{2-(r_1+r_3)/2}}{(2-r_1-r_3)(4-r_1-r_3)\sqrt{\alpha_1\alpha_3}}\,dx, \quad \text{ else},\\
\end{cases}
\end{align*}

ensures that the relation (<ref>) holds.
Since $V_1, V_2, V_3$ in Table <ref>
are positive and concave, $F$ is concave, and $G$ is convex,  the objection functional in the MFC problem (<ref>) is convex under a linear constraint. Hence there exists a minimizer.

When the terminal density $\rho(T,\cdot)=\rho^1$ is prescribed (or when $\gamma\rightarrow +\infty$
for the terminal function in Table <ref>), the problem (<ref>) is called the planning problem for MFC. This problem with $\beta=0$ (no $V_3$) was consider in the work [27].

Taking the regularization parameter $\beta>0$ in MFC (<ref>) forms a Schrödinger bridge type problem  [10, 26]. A related planning problem with $\beta>0$ and no potential  $F=0$ was considered recently in [9], known as unbalanced regularized optimal mass transport (urOMT). It has applications in cancer imaging.
Therein, $V_1(\rho)=V_3(\rho)=\rho$, $V_2(\rho)=\alpha(t,x)\rho$ for $\alpha\ge 0$,
and the objective function has an additional term
$\int_0^T\int_{\Omega}\beta s\log(\rho)dxdt$ compared with our formulation (<ref>) under the same linear constraint (<ref>).

§.§ MFC for reaction-diffusion systems

In this subsection, we introduce a system of MFC problem.
It can be viewed as a generalization of the system distance (<ref>).

Suppose there exist $M\in\mathbb{N}_+$ species and $R\in\mathbb{N}_+$ reactions with reaction coefficient matrix $\Gamma=(\gamma_{i,p})\in\mathbb{R}^{M\times R}$. Let $V_{1,i}:\mathbb{R}_+\rightarrow\mathbb{R}_+$
for $1\le i\le M$,
for $1\le p\le R$
be positive mobility functions.
Take terminal time $T>0$, and non-negative regularization parameter $\beta\ge0$. Choose energy functionals $\mathcal{E}_i\colon \mathcal{M}\rightarrow\mathbb{R}$ for $1\le i\le M$,
an potential functional $\bm{\mathcal{F}}\colon \mathcal{M}^M\rightarrow\mathbb{R}$,
and a terminal functional
$\bm{\mathcal{G}}\colon \mathcal{M}^M\rightarrow\mathbb{R}$.

\begin{equation}\label{vmfcA1}
\begin{split}
&\inf_{\bmr, \tilde{\cmm}, \tilde\bso}\quad\int_0^T\left[\int_\Omega
\left(\sum_{i=1}^M\frac{|\tilde\bmm_i|^2}{2V_{1,i}(\rho_i)}+
\sum_{p=1}^R\frac{|\tilde s_p|^2}{2V_{2,p}(\bmr)}\right)dx
\end{split}
\end{equation}

where the infimum is among all densities
$\bmr=(\rho_1,\cdots,\rho_M)$ with $\bmr(t,\cdot)\in\mathcal{M}^M$, fluxes
$\tilde\cmm=(\tilde\bmm_1,\cdots,\tilde\bmm_M)$ with $\tilde{\bmm}_i\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and sources $\tilde \bso=(\tilde s_1, \cdots, \tilde s_M)\colon [0, T]\times \Omega \rightarrow\mathbb{R}^M$, such that

\begin{equation}\label{smfcA2}
\begin{split}
&\partial_t \rho_i + \nabla\cdot \tilde{\bmm}_i-
\sum_{p=1}^R\gamma_{i,p}\tilde s_p\\
\left(\nabla\cdot(V_{1,i}(\rho_i)\nabla\frac{\delta}{\delta \rho}\mathcal{E}_i(\rho_i))-
\sum_{p=1}^R\gamma_{i,p}V_{2,p}(\bmr)\sum_{j=1}^M\gamma_{j,p}\frac{\delta}{\delta\rho}\mathcal{E}_j(\rho_j)\right),
\end{split}
\end{equation}

with boundary condition

\begin{align}
\label{smfcA3}
\left.\left(\tilde{\bmm}_i-\beta V_{1,i}(\rho_i)\nabla\frac{\delta}{\delta \rho}\mathcal{E}_i(\rho_i)\right)\cdot\nnu\right|_{\partial\Omega} =0,
\end{align}

for all $1\le i\le M$,
and fixed initial density $\bmr(0,\cdot) = \bmr^0\in\mathcal{M}^M$.

Again, our numerical scheme is based on the following equivalent formulation after the change of variables, whose proof is given in the Appendix.

For $1\le i\le M$ and $1\le p\le R$, denote $\bmm_i\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$
and $s_p\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
\bmm_i=\tilde\bmm_i-\beta V_{1_i}(\rho_i)\nabla\frac{\delta}{\delta \rho}\mathcal{E}_i(\rho_i),\quad
s_p=\tilde s_p-
\beta V_{2,p}(\bmr)
\sum_{j=1}^M\gamma_{j,p}\frac{\delta}{\delta\rho}\mathcal{E}_j(\rho_j),
\end{equation*}

and denote

\begin{equation}
\label{v3-forms}
V_{3,i}(\rho_i)=\frac{1}{V_{1,i}(\rho_i)|\frac{\delta^2}{\delta \rho^2}\mathcal{E}_i(\rho_i)|^2}.
\end{equation}

Then, the system MFC problem (<ref>) is equivalent to the following optimization problem:

\begin{equation}\label{vmfc2A1}
\begin{split}
&\inf_{\rho,\bmm, s} ~~ \int_0^T\int_\Omega
\sum_{i=1}^M\left(\frac{|\bmm_i|^2}{2V_{1,i}(\rho_i)}+\beta^2\frac{|\nabla \rho_i|^2}{2V_{3,i}(\rho_i)}\right)dxdt\\
\sum_{p=1}^R\left(\frac{|s_p|^2}{2V_{2,p}(\bmr)}
+\frac{\beta^2}{2}{|\sum_{j=1}^M\gamma_{j,p}\frac{\delta\mathcal{E}_j}{\delta \rho}(\rho_j)|^2}{V_{2,p}(\bmr)}
\right)dxdt\\
\mathcal{E}_i(\rho_{i}^0)],
\end{split}
\end{equation}

subject to the constraints

\begin{equation}\label{vmfc2A2}
\partial_t \rho_i(t, x) + \nabla\cdot \bmm_i(t,x)-\sum_{p=1}^R\gamma_{i,p}s_p(t,x)=0,
\end{equation}

with no-flux boundary condition
\bmm_i\cdot\nnu|_{\partial\Omega}=0,
for all $1\le i\le M$,

and fixed initial density $\bmr(0,\cdot)=\bmr_0=(\rho_{1}^0,\cdots, \rho_{M}^0)\in\mathcal{M}^M$.

Note that when $M=1$ and $R=1$, the MFC problem (<ref>) corresponds to the scalar case considered in Definition (<ref>).

Similar to the scalar case, we assume that there are two functions $\bm F:[0,T]\times\Omega\times \mathbb{R}_+^M\rightarrow \mathbb{R}$
and $\bm G:\Omega\times\mathbb{R}_+^M\rightarrow \mathbb{R}$, such that

\begin{align*}
\int_{\Omega}\bm F(t, x, \bmr)dx = &\;
\bm{\mathcal{F}}(\bmr(t,\cdot))-\int_{\Omega}\sum_{p=1}^{R}\frac{\beta^2}{2}{|\sum_{j=1}^{M}\frac{\delta\mathcal{E}_j}{\delta \rho}(\rho_j)|^2}{V_{2,p}(\bmr)}dx,\\
\int_{\Omega}\bm G(x,\bmr)dx = &\;
\bm{\mathcal{G}}(\bmr)+\beta\sum_{i=1}^M\mathcal{E}_i(\rho_i).
\end{align*}

This simplifies the functional in (<ref>) as follows:

\begin{equation*}
\begin{split}
&\int_0^T\int_\Omega \left[
\sum_{i=1}^M\left(\frac{|\bmm_i|^2}{2V_{1,i}(\rho_i)}+\beta^2\frac{|\nabla \rho_i|^2}{2V_{3,i}(\rho_i)}\right)+ \sum_{p=1}^R\frac{|s_p|^2}{2V_{2,p}(\bmr)}-\bm{F}(t,x,\bmr)\right]dxdt\\
\int_{\Omega}
\bm{G}(x, \bmr(T,x))dx,
\end{split}
\end{equation*}

The KKT system for the optimization problem (<ref>) is given below. Its proof is again presented in the Appendix.

Let $(\bmr, \cmm, \bso)$ be the critical point of the optimization problem (<ref>). Then there exists a function $\bm{\phi}=(\phi_1,\cdots,\phi_M)\colon [0, T]\times\Omega\rightarrow \mathbb{R}^M$, such that

\begin{equation*}
\frac{\bmm_i(t,x)}{V_{1,i}(\rho_i(t,x))}=\nabla \phi_i(t,x),\quad \frac{s_p(t,x)}{V_{2,p}(\bmr(t,x))}=
\sum_{j=1}^M\gamma_{j,p}\phi_j(t,x),
\end{equation*}

\begin{equation}\label{SBPI}
\left\{\begin{aligned}
&\partial_t\rho_i +\nabla\cdot(V_{1,i}(\rho_i)\nabla \phi_i)-\sum_{p=1}^RV_{2,p}(\bmr)\gamma_{i,p}\sum_{j=1}^M\gamma_{j,p}\phi_j=0,\\
&\partial_t\phi_i+\frac{1}{2}\|\nabla \phi_i\|^2 V_{1,i}'(\rho_i)+\frac{1}{2}\sum_{p=1}^R|\sum_{j=1}^M\gamma_{j,p}\phi_j|^2
\frac{\partial}{\partial \rho_i}V_{2,p}(\bmr)\\
&\hspace{6.3cm}+\frac{\delta}{\delta \rho_i}\Big[\bm{\mathcal{F}}(\bmr)-\frac{\beta^2}{2}\bm{\mathcal{I}}(\bmr)\Big]=0,
\end{aligned}\right.
\end{equation}

where $\bm{\mathcal{I}}(\bmr)$
is the generalized Fisher information functional given in (<ref>),

with initial and terminal time boundary conditions

\begin{equation}
\rho_i(0,x)=\rho_i^0(x),\quad \phi_i(T,x)=-\frac{\delta}{\delta \rho_i}\Big(\bm{\mathcal{G}}(\bmr(T,\cdot))+\beta \mathcal{E}_i(\rho_i(T,\cdot))\Big).
\end{equation}

§.§.§ Examples

Let us consider examples with only pairwise interactions. This setting is similar to the reversible Markov chains on discrete states (symmetric weighted graph) [19, 34, 33]. We assume the mobility function $V_{2,p}$ only depends on two densities $\rho_{p_0}$ and $\rho_{p_1}$
for $p_0\not= p_1\in\{1,\cdots, M\}$ for $1\le p\le R$.
We take the reaction coefficient $\Gamma=(\gamma_{i,p})\in \mathbb{R}^{M\times R}$, such that

\begin{align}
\label{rr}
\gamma_{i,p}= \begin{cases}
1 & \text{if } i=p_0,\\
-1 & \text{if } i=p_1,\\
0 & \text{else}.
\end{cases}
\end{align}

These reactions naturally lead to a graph with $M$
vertices $\{X_1,\cdots, X_M\}$ and
$R$ edges $\{E_1,\cdots, E_R\}$ with
$E_p=(X_{p_0}, X_{p_1})$ connecting vertices $X_{p_0}$ and $X_{p_1}$, where each vertex represents a density species, and each edge represents a reaction between two densities.

Different choices of the functions
$V_{1,i}$, $V_{2,p}$, $V_{3,i}$, $\bm F$, $\bm G$,
along with the reaction coefficient matrix $\Lambda\in\mathbb{R}^{M\times R}$ in (<ref>),
regularization parameter $\beta\ge0$ and initial density $\bmr^0$,
lead to various system MFC problems (<ref>). We present some examples of these functions in Table <ref> below.

$V_{1,i}(\rho_i)$:     $\alpha_{1}\rho_i^{\gamma_1}$ with $ 0\le \gamma_1\le 1$. $\alpha_1>0$.

$V_{2,p}(\bmr)$:     $\alpha_2(\rho_{p_0}+\rho_{p_1})/2$, or

$\alpha_2\frac{\rho_{p_0}\rho_{p_1}}{\rho_{p_0}+\rho_{p_1}}$, or
$\alpha_2\frac{\rho_{p_0}-\rho_{p_1}}{\log(\rho_{p_0})-\log(\rho_{p_1})}$. $\alpha_2\ge 0$

$V_{3,i}(\rho_i)$:     $\alpha_3\rho_i^{\gamma_3}$ with $ 0\le \gamma_3\le 1$. $\alpha_3>0$.

$\bm F(t,x,\bmr)$:     $-
\sum_{i=1}^M\tau_i\rho_i(\log(\rho_i)+V_i(t,x))$,

$-\sum_{i=1}^M\tau_i(\rho_i^m+\rho_i V_i(t,x))$.
$\tau_i\ge 0$, $m> 1$.

$G(x,\rho)$:     $\gamma\sum_{i=1}^M(\rho_i-\rho_i^1)^2$, or
$\gamma> 0$.

Example of functions for system MFC (<ref>).

Here we can again take the energies

\begin{align*}
\mathcal{E}_i(\rho_i)=\begin{cases}
\int_{\Omega}\frac{1}{\sqrt{\alpha_1\alpha_3}}\rho_i\log(\rho_i)\,dx, \quad \hspace{1.4cm}\text{ if }r_1=r_3=1,\\[1ex]
\int_{\Omega}\frac{4\rho_i^{2-(r_1+r_3)/2}}{(2-r_1-r_3)(4-r_1-r_3)\sqrt{\alpha_1\alpha_3}}\,dx, \quad \text{ else },\\
\end{cases}
\end{align*}

so that the relations (<ref>) hold.
Here, the four choices for $V_{2,p}$ correspond to arithmetic, geometric, harmonic, and logarithmic averages of two densities; see [19, 34, 33].
Similar to the scalar case, $V_{1,i}, V_{2,p}, V_{3,i}$ in Table <ref>
are positive and concave, $\bm F$ is concave, and $\bm G$ is convex. Hence the objective functional in MFC problem (<ref>) is convex under a linear constraint. A minimizer is guaranteed.

This formulation naturally allows general nonlinear reactions among more than two densities. We leave the detailed modeling study of MFC for more general reaction-diffusion systems in future work.

§ HIGH ORDER DISCRETIZATIONS AND OPTIMIZATION ALGORITHMS

In this section, we present the high-order spatial-time finite element formulations for the proposed MFC problems (<ref>) and (<ref>).

§.§ Scalar case

In this subsection, we discretize the scalar model (<ref>).
By introducing a new approximation variable $\nn=\beta\nabla\rrho$, the optimization problem (<ref>) is rewritten as follows.

\begin{equation}\label{MS1}
\inf_{\rrho,\mm, \so, \nn} ~ \int_0^T\int_\Omega \Big[\frac{\|\mm\|^2}{2 V_1(\rrho)}+ \frac{|\so|^2}{2V_2(\rrho)}
-\FF(\rrho)\Big]dxdt +
\int_{\Omega}\GG(\rrho(T,x))dx,
\end{equation}

where the infimum is taken over density $\rrho\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, flux $\mm\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, source $\so\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, and
vector field
$\nn\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$,
such that

\begin{align}
\label{MS2}
\partial_t \rrho + \nabla\cdot \mm-\so&\;=0, \\
\label{MS3}
\nn - \beta \nabla \rrho &\;= 0,
\end{align}

on $[0,T]\times \Omega$,
with fixed initial density $\rrho(0,x)=\rrho^0(x)$ in $\Omega$
and no-flux boundary condition $\mm\cdot\nnu=0$ on $[0,T]\times \partial\Omega$.
Here we suppress the coordinate dependence of functions $F$ and $G$ for simplicity of presentation.

Following the high-order finite element scheme for MFC problems in [16], we discretize the constraint optimization problem (<ref>) using high-order finite element methods
and apply the ALG2 optimization algorithm to solve the discrete saddle-point problem.

§.§.§ The continuous saddle-point formulation

We introduce a scalar Lagrange multiplier $\pphi\in H^1(\Omega_T)$, where $\Omega_T:=[0,T]\times\Omega$ is the space-time domain,
for the constraint (<ref>),
and a vectorial Lagrange multiplier $\ssigma\in L^2([0,T])\otimes H(\mathrm{div}_0;\Omega)$ for the constraint (<ref>),  where

\[
\ttau\in H(\mathrm{div};\Omega):
\ttau\cdot\nnu|_{\partial\Omega} = 0\},
\]

in which the H(div)-conforming space

\[
\ttau\in [L^2(\Omega)]^d:
\nabla\cdot\ttau\in L^2(\Omega)
\}.
\]

We reformulate the constrained optimization problem (<ref>)
into the following saddle-point system:
Find the critical point of the system

\begin{equation}\label{SD-scalarZ}
\begin{split}
\inf_{\rrho,\mm, \so, \nn, \rrho_T} \sup_{\pphi, \ssigma}~& \int_0^T\int_\Omega \Big[\frac{\|\mm\|^2}{2 V_1(\rrho)}+ \frac{|\so|^2}{2V_2(\rrho)}
-\FF(\rrho)\Big]dxdt \\
& \hspace{-0.8cm}
-\int_0^T\int_\Omega \Big[(\rrho\, \partial_t \pphi
+\mm\cdot\nabla \pphi+\so\pphi)
\Big]
& \hspace{-0.8cm}+
\int_\Omega \left(\rrho_T\, \pphi(T,\cdot)
\end{split}
\end{equation}

the variables $\mm,\nn\in [L^2(\Omega_T)]^d$, $\rrho, \so\in L^2(\Omega_T)$,
$\rho_T\in L^2(\Omega)$
with $\rrho,\rrho_T\ge 0$ a.e., $\pphi\in H^1(\Omega_T)$,
and $\ssigma\in L^2([0,T])\otimes H(\mathrm{div}_0;\Omega)$.
Notice that integration by parts is used to incorporate the linear constraints (<ref>)–(<ref>) to the saddle point problem (<ref>). Moreover, since no derivative information is needed for the physical variables $\mm, \nn, \rrho$ and $\so$, it is natural to approximate them using $L^2$-conforming spaces.

To simplify the notation, we collect the variables

\begin{align}
\label{sc1}
\uu:=\;(\rrho, \mm, \so, \nn), \text{ and }
\Pphi:=\;(\pphi, \ssigma),
\end{align}

and define the differential operator

\begin{equation}
\label{sc2}
\mathcal{D}(\Pphi):=
(\partial_t\pphi+\beta\nabla\cdot\ssigma, \nabla\pphi, \pphi, \ssigma)\in [L^2( \Omega_T)]^{2d+2}.
\end{equation}

With this notation, the above saddle-point problem simplifies to the following:
Find the critical point of

\begin{equation}\label{SD-scalar}
\begin{split}
\inf_{\uu, \rrho_T} \sup_{\Pphi}~& \int_0^T\int_\Omega\Big[\HH(\uu)
\,dxdt \\
& \hspace{-0.2cm}+
\int_\Omega \Big[\GG(\rrho_T)+\rrho_T\, \pphi(T,\cdot)
\end{split}
\end{equation}

the variables $\uu\in [L^2(\Omega_T)]^{2d+2}$, $\rho_T\in L^2(\Omega)$
with $\rrho, \rrho_T\ge 0 $ a.e., and
$\Pphi\in H^1(\Omega_T) \times  \Big(L^2([0,T])\otimes H(\mathrm{div}_0;\Omega)\Big)$,
the nonlinear functional $\HH(\uu)$ is given by

\[
\HH(\uu):=
\frac{\|\mm\|^2}{2 V_1(\rrho)}+ \frac{|\so|^2}{2V_2(\rrho)}
\]

and the inner-product is defined by

\[
\uu\cdot\mathcal{D}(\Pphi):=
\rrho (\partial_t \pphi+\beta\nabla\cdot\ssigma)
+\mm\cdot\nabla \pphi+\so\,\pphi+\nn\cdot\ssigma.
\]

§.§.§ The high-order finite element scheme and numerical integration

Following the work [16], we approximate the physical variables $\mm, \nn, \so$, $\rrho$ using (discontinuous) high-order integration rule spaces,
the variable $\pphi$ using a high-order $H^1$-conforming finite element space, and
the variable $\ssigma$ using a vectorial $H^1$-conforming finite element space.

Let $\Omega_h=\{T_\ell\}_{\ell=1}^{N_S}$ be a conforming mesh of the spatial domain $\Omega$ with $N_S$ elements, and
$I_h=\{I_j\}_{j=1}^{N_T}$ be a uniform discretization of the temporal domain $[0, T]$
with $N_T$ line segments. We denote the space-time mesh $\Omega_{T,h}:=I_h\otimes\Omega_h$.
For simplicity, we restrict ourselves to the case where the domain $\Omega=[0,1]^d$ is a unit hypercube and its spatial mesh $\Omega_h$ is a Cartesian mesh with uniform hypercubic elements.
Given a polynomial degree $k\ge 1$,
we denote the following $H^1$
and $L^2$ finite element spaces:

\begin{align}
\label{fes-V}
V_h^k:=&\;\{v\in H^1(\Omega_T): \;\; v|_{I_j\times T_\ell}
\in Q^k(I_j)\otimes Q^k(T_\ell)\;\;
\forall j, \ell\},\\
\label{fes-W}
W_h^{k-1}:=&\;\{w\in L^2(\Omega_T): \;\; w|_{I_j\times T_\ell}
\in Q^{k-1}(I_j)\otimes Q^{k-1}(T_\ell)\;\;
\forall j, \ell\},\\
\label{fes-M}
M_h^{k-1}:=&\;\{\mu\in L^2(\Omega):\;\;\; \;\; \mu|_{T_\ell}
\in Q^{k-1}(T_\ell)\;\;\forall \ell\},
\end{align}

where $Q^k(S)$ is the tensor-product polynomial space of degree no greater than $k$ in each coordinate direction.

Using the above finite element spaces, we obtain the following discrete saddle point problem:
Find the critical point of the discrete system

\begin{equation}\label{SD-scalarh}
\begin{split}
\inf_{\uu_h, \,\rrho_{T,h}} \sup_{\Pphi_h}~& \int_0^T\int_\Omega\Big[\HH(\uu_h)
\,dxdt \\
& \hspace{-0.2cm}+
\int_\Omega \Big[\GG(\rrho_{T,h})+\rrho_{T,h}\, \pphi_h(T,\cdot)
\end{split}
\end{equation}

where the variables
$\uu_h:=(\rrho_h, \mm_h, \so_h, \nn_h)$
with $\mm_h,\nn_h\in [W_h^{k-1}]^d$, $\rrho, \so\in W_h^{k-1}$,
$\rrho\ge0$ a.e., $\rho_{T,h}\in M_h^{k-1}$
with $\rrho_{T,h}\ge 0$ a.e.,
and $\Pphi_h:=(\pphi_h,\ssigma_h)$
with $\pphi_h\in V_h^k$
and $\ssigma_h\in [V_h^{k}]^d$
with $\ssigma_h\cdot\nnu|_{\partial\Omega}=0$.

In the scheme (<ref>), we use $L^2$-conforming spaces to approximate the
physical variables $\uu_h$ and $\rrho_{T,h}$, which is  natural as no derivative information is needed. Moreover, we use $H^1(\Omega_T)$-conforming finite element spaces to approximate the Lagrange multipliers $\Pphi_h$. The use of $H^1$-conforming space to approximate the scalar variable $\pphi_h$ is natural as the formulation (<ref>) requires space-time derivative information for $\pphi_h$.
On the other hand, the use of $H^1(\Omega_T)$-conforming space for the vectorial Lagrange multiplier $\ssigma_h$ is only for the simplicity of implementation.
A more natural finite element space for this variable is the following:

\[
[W_h^k]^d\cap \left(L^2([0,T])\otimes H(\mathrm{div}_0;\Omega)\right),
\]

which is discontinuous in time and $H(\mathrm{div})$-conforming in space.
We will implement this choice in future work.

Here we use polynomial degree of one order lower to approximate the physical variables $\uu_h$ than that for the Lagrange multipliers $\Pphi_h$.
Our lowest order scheme is a staggered grid scheme where the physical variables stay in cell centers and the Lagrange multipliers stay on mesh vertices.

In practice, the integrals in (<ref>) are discretized using
Gauss-Legendre quadrature rules with $k$ integration points per direction.
The total number of integration points for the space-time integral on $[0,T]\times \Omega$ is $(kN_T)\times (k^{d}N_S)$, and that for the spatial integral on $\Omega$
is $k^{d}N_S$, where  $N_T$ is the number of temporal cells and $N_T$ is the number of spatial cells.
Denoting these spatial/temporal Gauss-Legendre quadrature points as

\begin{equation}
\label{quad-pt}
\{\bm{\xi}_i\}_{i=1}^{k^dN_S}, \text{ and } \{\eta_j\}_{j=1}^{kN_T}
\end{equation}

with their associated quadrature weights as

\begin{equation}
\label{quad-wt}
\{\omega_i\}_{i=1}^{k^dN_S}, \text{ and } \{\zeta_j\}_{j=1}^{kN_T},
\end{equation}

we write the discrete integrals as follows:

\begin{align}
\label{d-int1}
\stint{f(t,x)}:=&\;
\sum_{j=1}^{kN_T}
\sum_{i=1}^{k^dN_S}f(\eta_j, \bm{\xi}_i)\omega_i\zeta_j,\\
\sint{f(x)}:=&\;
\sum_{i=1}^{k^dN_S}f(\bm{\xi}_i)\omega_i.
% \langle\!\langle f(t,x)\rangle\!\rangle_{h}
\label{d-int2}
\end{align}

Furthermore, we use the Gauss-Legendre basis for the $L^2$-conforming spaces $W_h^{k-1}$
and $M_h^{k-1}$ in the implementation, i.e., a function $w_h\in W_h^{k-1}$
is expressed as

\[
w_h =
\sum_{j=1}^{kN_T}
\sum_{i=1}^{k^dN_S} \mathsf{w}_{ij}\,\psi_j(t)\varphi_i(x),
\]

with coefficient $\mathsf{w}_{ij}\in \mathbb{R}$.
Here the tensor-product Gauss-Legendre basis function
$\psi_j(t)\varphi_i(x)\in W_h^{k-1}$ satisfies the nodal interpolation property

\[
\psi_j(\zeta_{j'})=\delta_{jj'}, \text{ and }
\varphi_i(\bm\xi_{i'})=\delta_{ii'},
\]

where $\delta_{ii'}=1$ if $i=i'$ and
$\delta_{ii'}=0$ if $i\not=i'$.
A function $v_h\in M_h^{k-1}$ is expressed as
v_h =
\sum_{i=1}^{k^dN_S} \mathsf{v}_{i}\varphi_i(x)
with coefficient $\mathsf{v}_i\in \mathbb{R}$.
Here $\mathrm{w}_{ij}=w_h(\zeta_j,\bm\xi_i)$ is the value of $w_h$
on the space-time integration point $(\zeta_j, \bm\xi_i)$, and
$\mathrm{v}_{i}=v_h(\bm\xi_i)$ is the value of $v_h$
on the spatial integration point $\bm\xi_i$.
Hence there holds

∑_i=1^k^dN_Sf(𝗐_ij)ω_iζ_j,        ∀w_h=∑_j=1^kN_T
∑_i=1^k^dN_S 𝗐_ij ψ_j(t)φ_i(x)∈W_h^k-1,

∑_i=1^k^dN_Sf(𝗏_i)ω_i,        ∀v_h
=∑_i=1^k^dN_S 𝗏_iφ_i(x)

With these notation, the fully discrete scheme (<ref>)
with numerical integration is given as follows:
Find the critical point of the fully discrete system

\begin{equation}\label{SD-scalarh1}
\begin{split}
\inf_{\uu_h, \,\rrho_{T,h}} \sup_{\Pphi_h}~& \stint{\HH(\uu_h)
\\
& \hspace{-0.2cm}+
\sint{\GG(\rrho_{T,h})+\rrho_{T,h}\, \pphi_h(T,\cdot)
\end{split}
\end{equation}

Here the unknowns $\rho_{T,h}=\sum_{i=1}^{k^dN_S}\uprho_{T,i}\varphi_i(x)$
with non-negative coefficient $\uprho_{T,i}\in\mathbb{R}_+$,
$\uu_h:=(\rrho_h, \mm_h, \so_h, \nn_h)$ with

_h=    ∑_j=1^kN_T
∑_i=1^k^dN_S_ij ψ_j(t)φ_i(x),  _h=       ∑_j=1^kN_T

_h=    ∑_j=1^kN_T
∑_i=1^k^dN_Ss_ijψ_j(t)φ_i(x),  _h=       ∑_j=1^kN_T

where coefficients $\uprho_{ij}\in \mathbb{R}_+$, $\bm{\mathrm{m}}_{ij},\bm{\mathrm{n}}_{ij}\in\mathbb{R}^d$, and $\mathrm{s}_{ij}\in\mathbb{R}$,
and $\Pphi_h:=(\pphi_h,\ssigma_h)\in [V_h^k]^{d+1}$
with $\ssigma_h\cdot\nnu|_{\partial\Omega}=0$.

We further note that the use of Gauss-Legendre integral rule along with the (nodal) Gauss-Legendre bases for the discontinuous spaces $W_h^{k-1}$ and $M_h^{k-1}$ significantly simplifies the optimization problem (<ref>), as no
explicit degrees of freedom coupling is introduced among the physical variables $\uu_h$ and
In particular, given a fixed Lagrange multiplier $\Pphi_h$, the nonlinear optimization problem (<ref>) for $\uu_h$ and
$\rrho_{T,h}$ can be solved in a pointwise fashion per integration point.
Moreover, positivity of the (space-time) density $\rrho_h$ and the terminal density
$\rrho_{T,h}$ is guaranteed on all quadrature points by design as the admissible set requires
their coefficients to be non-negative: $\uprho_{ij}, \uprho_{T,i}\ge0$.
The scheme (<ref>) does not prescribe the choice of basis functions for the $H^1$-conforming space $V_h^k$. The results are independent of the specific choice of bases for $V_h^k$. We use a nodal Gauss-Lobatto base in our implementation.

The scheme (<ref>) is our fully discrete high-order finite element discretization to the scalar constraint optimization problem (<ref>).
Next, we present its optimization solver using the ALG2 algorithm [15], which is also referred to as the alternating direction method of multipliers (ADMM) method [5].

§.§.§ Duality and augmented Lagrangian

We introduce the dual variables

\[
\uu_h^*=(\rrho_h^*, \mm_h^*, \so_h^*, \nn_h^*)\in [W_h^{k-1}]^{2d+2}, \text{ and }
\rho_{T,h}^*
\in M_h^{k-1}.
\]

To further simplify the notation, we denote

\[
\uu_h=\;\sum_{j=1}^{kN_T}
\sum_{i=1}^{k^dN_S}{\bm{\mathrm{u}}}_{ij}\,\psi_j(t)\varphi_i(x),\text{ and }
\uu_h^*=\;\sum_{j=1}^{kN_T}
\sum_{i=1}^{k^dN_S}{\bm{\mathrm{u}}}_{ij}^*\,\psi_j(t)\varphi_i(x),
\]

\[
\bm{\mathrm{u}}_{ij}:=
(\uprho_{ij}, \bm{\mathrm{m}}_{ij},\mathrm{s}_{ij}, \bm{\mathrm{n}}_{ij})
\in \mathbb{R}^{2d+2}, \text{ with }\uprho_{ij}\ge 0,
\]

\[
\bm{\mathrm{u}}_{ij}^*:=
(\uprho_{ij}^*, \bm{\mathrm{m}}_{ij}^*,\mathrm{s}_{ij}^*, \bm{\mathrm{n}}_{ij}^*)
\in \mathbb{R}^{2d+2}.
\]

Introducing the following Legendre transforms:

\begin{align}
\HH^*(\uui_{ij}^*):=&\; \sup_{\uui_{ij}\in \mathbb{R}^{2d+2},
\uprho_{ij}\ge 0} \uui_{ij}\cdot\uui_{ij}^*-
\HH(\uui_{ij}),\\
\GG^*(\rti^*):=&\; \sup_{\rti\in \mathbb{R}_+} \rti\cdot\rti^*-
\GG(\rti),
\end{align}

by duality, there holds

\begin{align}
\label{dual-H}
\HH(\uui_{ij}):=&\; \sup_{\uui_{ij}^*\in \mathbb{R}^{2d+2}} \uui_{ij}\cdot\uui_{ij}^*-
\HH^*(\uui_{ij}^*),\\
\label{dual-G}
\GG(\rti):=&\; \sup_{\rti^*\in \mathbb{R}} \rti\cdot\rti^*-
\GG^*(\rti^*).
\end{align}

Plugging the relations (<ref>) and (<ref>) back to the scheme (<ref>), we have the following dual formulation of (<ref>):
Find the critical point of

\begin{equation}\label{SD-scalarh2}
\begin{split}
\sup_{\uu_h, \,\rrho_{T,h}} \inf_{\Pphi_h,\uu_h^*, \,\rrho_{T,h}^*}~& \stint{\HH^*(\uu_h)
\\
& \hspace{-0.2cm}+
\sint{\GG^*(\rrho_{T,h}^*)-\rrho_{T,h}\, (\pphi_h(T,\cdot)+\rrho_{T,h}^*)
\end{split}
\end{equation}

It is clear that in the above formulation (<ref>),
the variable $\uu_h$ is the Lagrange multiplier for the constraint $\mathcal{D}(\Pphi_h)-\uu_h^*=0$, while the variable $\rrho_{T,h}$ is the Lagrange multiplier for the terminal constraint $\pphi_h(T,x)+\rho_{T,h}^*(x) = 0$.

Hence, the critical point of the system (<ref>) is equivalent to the
following augmented Lagrangian system

\begin{equation}\label{SD-scalarh3}
\begin{split}
\sup_{\uu_h, \,\rrho_{T,h}} \inf_{\Pphi_h,\uu_h^*, \,\rrho_{T,h}^*}~& \stint{\HH^*(\uu_h^*)
\\
& \hspace{-0.2cm}+
\sint{\GG^*(\rrho_{T,h}^*)-\rrho_{T,h}\, (\pphi_h(T,\cdot)+\rrho_{T,h}^*)
\frac{r}{2}\stint{(\mathcal{D}(\Pphi_h)-\uu_h^*)\cdot(\mathcal{D}(\Pphi_h)-\uu_h^*)}\\
\frac{r}{2}\sint{(\pphi_h(T,\cdot)+\rrho_{T,h}^*)^2},
\end{split}
\end{equation}

where $r>0$ is the augmented Lagrangian parameter.

§.§.§ The ALG2 algorithm and its efficient (modified) implementation

The ALG2 solves the optimization problem (<ref>) in a splitting fashion. One update from iteration level $\ell-1$ to $\ell$ contains three steps:

* Step A: Given data $\uu_h^{\ell-1},\rrho_{T,h}^{\ell-1},\uu_h^{*,\ell-1},\rrho_{T,h}^{*,\ell-1}$,
compute $\Phi_h^{\ell}$ by optimizing the target functional in (<ref>)
with respect to $\Phi_h$.

* Step B: Given data
compute $\uu_h^{*,\ell}$ and $\rrho_{T,h}^{*,\ell}$ by optimizing the target functional in (<ref>) with respect to $\uu_h^*$ and $\rrho_{T,h}^*$.

* Step C: Update the Lagrange multipliers
$\uu_h^{\ell},\rrho_{T,h}^{\ell}$ by the following explicit expressions:

\begin{align}
\label{sd-lag1}
\uu_h^{\ell}=&\;\uu_h^{\ell-1}+r(\mathcal{D}(\Pphi_h^{\ell})-\uu_h^{*,\ell}),\\
\label{sd-lag2}
\rrho_{T,h}^{\ell}=&\;\rrho_{T,h}^{\ell-1}-r(\pphi_h^{\ell}(T,\cdot)+\rrho_{T,h}^{*,\ell}).
\end{align}

Here Step A contains a global constant-coefficient linear system solver for the variables $\Pphi_h=(\pphi_h, \ssigma_h)$,
Step B is a pointwise nonlinear equation solved per quadrature point that involves the dual functions $\HH^*$ and $\GG^*$, and Step C is a simple pointwise update.
In practical implementation, we apply two further simplifications to the above ALG2 algorithm: (1) apply splitting to solve $\pphi_h$ and $\ssigma_h$ sequentially; (2) use duality to solve pointwise nonlinear problems for $\uu_h$ and $\rrho_{T,h}$, and then apply a simple pointwise update  for the dual variables $\uu_h^*$ and $\rrho_{T,h}^*$.
We note that the second modification avoids the
explicit computation of the dual functions, which has already been used in our previous work
on ALG2 for variational time implicit schemes for reaction-diffusion systems [17].
One iteration of the modified ALG2 algorithm is documented in Algorithm
<ref> below.
We skip the detailed derivation of each step. Since the whole algorithm is very similar to  Algorithm 2 in [17] with simple modifications.

One iteration of (modified) ALG2  for (<ref>).

[1]
Step A1.
Find $\pphi_h^{\ell}\in V_h^k$
such that, for all $\psi_h\in V_h^k$,

\begin{align*}
&    \stint{\partial_t\pphi_h^\ell\cdot\partial_t\psi_h+
\nabla\pphi_h^\ell\cdot\nabla\psi_h+\pphi_h^\ell\cdot\psi_h}
\stint{(\rrho_h^{*,\ell-1}-\frac{\rrho_h^{\ell-1}}{r}-\beta\nabla\cdot\ssigma_h^{\ell-1})\partial_t\psi_h+
\end{align*}

Step A2.
Find $\ssigma_h^{\ell}\in [V_h^k]^d$ with $\ssigma_h^\ell\cdot\nnu|_{\partial\Omega}=0$
such that

\begin{align*}
% \label{sd-sigma}
% \begin{split}
&    \stint{\beta^2(\nabla\cdot\ssigma_h^\ell)(\nabla\cdot\ttau_h)
\\
\stint{\beta(\rrho_h^{*,\ell-1}-\frac{\rrho_h^{\ell-1}}{r}-\partial_t\pphi_h^{\ell})\nabla\cdot\ttau_h+
% \end{split}
\end{align*}

for all $\ttau_h\in [V_h^k]^d$ with $\ttau_h\cdot\nnu|_{\partial\Omega}=0$.
Step A3. Evaluate $\mathcal{D}(\Pphi_h^\ell)$ on each space-time quadrature point $(\zeta_j, \bm\xi_i)$, and $\pphi_h^\ell(T,\cdot)$ on each spatial quadrature point $\bm\xi_i$.
Denote these values as

\begin{align*}
\mathcal{D}\Pphi_{ij}^\ell:=(\mathrm{d}\uprho_{ij}^\ell,
\bm{\mathrm{dm}}_{ij}^\ell,
\mathrm{ds}_{ij}^\ell,
\bm{\mathrm{dn}}_{ij}^\ell),\quad
\mathrm{d}\uprho_{T,i}^\ell=\pphi_h^\ell(T,\bm\xi_i).
\end{align*}

\overline{\bm{\mathrm{m}}}_{ij}^\ell,
\overline{{\mathrm{s}}}_{ij}^\ell,
\overline{\bm{\mathrm{n}}}_{ij}^\ell):=
\mathcal{D}\Phi_{ij}^\ell+\uui_{ij}^{\ell-1}/r,
\rti^{\ell-1}/r-\mathrm{d}\uprho_{T,i}^\ell.

Step B1. For each quadrature point $(\zeta_j,\bm\xi_i)$, find the nonnegative
density coefficient $\uprho_{ij}^\ell\in \mathbb{R}_+$ such that it minimizes the function
\begin{align*}
\frac{r^2|\overline{\bm{\mathrm{m}}}_{ij}^\ell|^2}{r+V_1(\uprho)}
\end{align*}

and find the nonnegative terminal density coefficient
$\uprho_{T,i}^\ell\in \mathbb{R}_+$ such that it minimizes the function
Step B2. Update the following coefficients:

\begin{align*}
\bm{\mathrm{m}}_{ij}^\ell=\;\frac{rV_1(\uprho_{ij}^\ell)}{r+V_1(\uprho_{ij}^\ell)}
\overline{\bm{\mathrm{m}}}_{ij}^\ell,\;\;
\overline{{\mathrm{s}}}_{ij}^\ell,\;\;
\bm{\mathrm{n}}_{ij}^\ell=\;\frac{rV_3(\uprho_{ij}^\ell)}{r+V_3(\uprho_{ij}^\ell)}
\overline{\bm{\mathrm{n}}}_{ij}^\ell.
\end{align*}

Step C. update the dual variables according to the following:

\begin{align*}
\uu_{ij}^{*,\ell} = \overline{\uu}_{ij}^{\ell}-\uu_{ij}^{\ell}/r, \quad\uprho_{T,i}^{*,\ell} =
\overline{\rrho}_{T,i}^{\ell}-\rrho_{T,i}^{\ell}/r.
\end{align*}

Step A1 involves solving a symmetric positive definite linear system
for a constant-coefficient reaction-diffusion equation, for which we use preconditioned conjugate gradient (PCG) method with a geometric multigrid preconditioner. It achieves a linear computational complexity $\mathcal{O}(N)$ with $N=\dim V_h^k$ being the total number of degrees of freedom.
Step A2 is an $H(div)$-elliptic linear system problem.

This system is well-conditioned for small parameter $\beta\ll 1$, which is usually
the case in our applications. We use the PCG method with a Jacobi preconditioner, which achieves optimal linear computational complexity $\mathcal{O}(dN)$.
Step B and C involve pointwise updates, which again have a linear complexity of $\mathcal{O}(N)$.
Hence the overall computational complexity of applying one iteration of Algorithm <ref>
is linear.
It is also clear that Step B1 guarantees the positivity of density on all quadrature points.

Unique solvability of the one-dimensional minimization problem in Step B1 of Algorithm <ref> can be guaranteed when the functions $L_{ij}$ and $L_{T,i}$
are strongly convex.
Strong convexity is ensured if the mobility functions $V_1, V_2$, and $V_3$ are positive and concave functions,  the potential function $F(\uprho)$ is concave, and
the function $G(\uprho)$ is convex, which are satisfied by the choices in Table <ref>.

§.§ System case

In this subsection, we discretize the model (<ref>) with $M$ species and $R$ reactions.
We note that due to our formulation, there are no additional technical difficulties for the system derivation compared with the scalar case. We still document details of the finite element scheme and its splitting optimization solver for the completeness of this paper.

Again, we introduce the vectors $\nn_i:=\beta\nabla\rrho_i$ to avoid derivative evaluation of the densities for the scheme.
To simplify the notation, we
denote the collection of densities $\brrho=(\rrho_1,\cdots,\rrho_M)$,
fluxes $\cmm=(\mm_1,\cdots,\mm_M)$, vectors $\bnn=(\nn_1,\cdots,\nn_M)$, and sources $\bso=(\so_1,\cdots, \so_R)$.
Then optimization problem (<ref>) is rewritten as follows:

\begin{equation}\label{MY1}
\begin{split}
\inf_{\brrho,\cmm, \bso, \bnn} ~& \int_0^T\int_\Omega \left[\sum_{i=1}^M\left(\frac{\|\mm_i\|^2}{2 V_{1,i}(\rrho_i)}
+ \sum_{p=1}^R\frac{|\sop|^2}{2V_{2,p}(\brrho)}\right] \, dxdt \\
\int_0^T\int_\Omega \bFF(\brrho)\,dxdt
\int_{\Omega}\bGG(\brrho(T,x))dx,
\end{split}
\end{equation}

subject to constraints

\begin{align}
\label{MY2}
\partial_t \rrho_i + \nabla\cdot \mm_i-\sum_{p=1}^R\gamma_{i,p}\so_p &\;=0, \\
\label{MY3}
\nn_i - \beta \nabla \rrho_i &\;= 0,
\end{align}

for all $1\le i\le M$ on $[0,T]\times \Omega$,
with fixed initial densities $\rrho_i(0,x)=\rrho_{i}^0(x)$ in $\Omega$
and no-flux boundary conditions $\mm_i\cdot\nnu=0$ on $[0,T]\times \partial\Omega$.

Similar to the scalar case in (<ref>), we reformulate the constrained optimization problem (<ref>)
into the following saddle-point system:
Find the critical point of

\begin{equation}\label{SY-scalar}
\begin{split}
\inf_{\buu, \brrho_T} \sup_{\bPphi}~& \int_0^T\int_\Omega\Big[\bHH(\buu)
\,dxdt \\
& \hspace{-0.2cm}+
\int_\Omega \Big[\bGG(\brrho_T)+\brrho_T\cdot \bpphi(T,\cdot)
\end{split}
\end{equation}

the variables

\begin{align}
\label{sy1}
&\buu:=\;(\brrho, \cmm, \bso, \bnn)\in [L^2(\Omega_T)]^{(2d+1)M+R}, \\
&\brrho_T:=\;(\rrho_{1,T}, \cdots, \rrho_{M,T})\in [L^2(\Omega)]^{M},\\
&\Pphi:=\;(\bpphi, \bssigma)\in [H^1(\Omega_T)]^M\times
\left[
L^2([0,T])\otimes H(\mathrm{div}_0;\Omega)
\right]^M,
\end{align}

\[
\bpphi=(\pphi_1,\cdots,\pphi_M), \text{ and }
\bssigma=(\ssigma_1,\cdots,\ssigma_M),
\]

and initial data $\brrho^0:=(\rrho_{1}^0,\cdots, \rrho_{M}^0)$.
Here the nonlinear function

\[
\bHH(\buu):=
\sum_{i=1}^M\left(\frac{\|\mm_i\|^2}{2 V_{1,i}(\rrho_i)}+\frac{|\nn_i|^2}{2V_{3,i}(\rrho_i)}\right)
+ \sum_{p=1}^R\frac{|\so_p|^2}{2V_{2,p}(\brrho)}
\]

and the operator

\[
\underline{\mathcal{D}}(\bPphi):=
(\partial_t\bpphi+\beta\nabla\cdot\bssigma, \nabla\bpphi, \Gamma^T\bpphi, \bssigma),
\]

\nabla\cdot\bssigma:=(\nabla\cdot\ssigma_1,\cdots, \nabla\cdot\ssigma_M)
is the component-wise divergence and
\nabla\bpphi:=(\nabla\pphi_1,\cdots, \nabla\pphi_M)
is the component-wise gradient, and
\Gamma = (\gamma_{i,p})\in \mathbb{R}^{M\times R}
is the reaction coefficient matrix.
Note that the inner-product
has the following component-wise form

\[
\uu\cdot\mathcal{D}(\Pphi)=
\sum_{i=1}^M\left(\rrho_i (\partial_t \pphi_i+\beta\nabla\cdot\ssigma_i)
+\mm_i\cdot\nabla \pphi_i+\nn_i\cdot\ssigma_i\right)
\]

Replacing the function spaces with appropriate (high-order) finite element spaces and applying numerical integration,
we derive the following fully discrete scheme:
Find the critical point of the fully discrete system

\begin{equation}\label{SY-scalarh1}
\begin{split}
\inf_{\buu_h, \,\brrho_{T,h}} \sup_{\bPphi_h}~& \stint{\bHH(\buu_h)
\\
& \hspace{-0.2cm}+
\sint{\bGG(\brrho_{T,h})+\brrho_{T,h}\cdot \bpphi_h(T,\cdot)
\end{split}
\end{equation}

Here the unknowns $\brrho_{T,h}=\sum_{\ii=1}^{k^dN_S}{\bm\uprho}_{T,\ii}\varphi_\ii(x)$
with non-negative coefficient ${\bm\uprho}_{T,i}\in\mathbb{R}_+^M$,
$\buu_h:=(\brrho_h, \cmm_h, \bso_h, \bnn_h)$ with

_h=    ∑_=1^kN_T
∑_=1^k^dN_S_ij ψ_j(t)φ_i(x),  _h=       ∑_j=1^kN_T

_h=    ∑_j=1^kN_T
∑_i=1^k^dN_Ss_ijψ_j(t)φ_i(x),  _h=       ∑_j=1^kN_T

where coefficients $\bm\uprho_{ij}\in \mathbb{R}_+^M$, $\underline{\bm{\mathrm{m}}}_{ij},\underline{\bm{\mathrm{n}}}_{ij}\in\mathbb{R}^{dM}$, and $\bm{\mathrm{s}}_{ij}\in\mathbb{R}^R$,
and $\bPphi_h:=(\bpphi_h,\bssigma_h)\in [V_h^k]^{(d+1)M}$
with $\bssigma_h\cdot\nnu|_{\partial\Omega}=0$.

We introduce the dual variables

\[
\buu_h^*=(\brrho_h^*, \cmm_h^*, \bso_h^*, \bnn_h^*)\in [W_h^{k-1}]^{(2d+1)M+R},
\]

\in M_h^{k-1}$, and denote

\[
\buu_h=\;\sum_{j=1}^{kN_T}
\sum_{i=1}^{k^dN_S}\underline{\bm{\mathrm{u}}}_{ij}\,\psi_j(t)\varphi_i(x),\text{ and }
\buu_h^*=\;\sum_{j=1}^{kN_T}
\sum_{i=1}^{k^dN_S}\underline{\bm{\mathrm{u}}}_{ij}^*\,\psi_j(t)\varphi_i(x),
\]

with coefficients

\[
\underline{\bm{\mathrm{u}}}_{ij}:=
(\bm\uprho_{ij}, \underline{\bm{\mathrm{m}}}_{ij},\bm{\mathrm{s}}_{ij}, \underline{\bm{\mathrm{n}}}_{ij})
\in \mathbb{R}^{(2d+1)M+R}, \text{ with }\bm\uprho_{ij}\ge 0,
\]

\[
\underline{\bm{\mathrm{u}}}_{ij}^*:=
(\bm\uprho_{ij}^*, \underline{\bm{\mathrm{m}}}_{ij}^*,\bm{\mathrm{s}}_{ij}^*, \underline{\bm{\mathrm{n}}}_{ij}^*)
\in \mathbb{R}^{(2d+1)M+R}.
\]

Using these dual variables,
we obtain the following augmented Lagrangian formulation of (<ref>):

\begin{equation}\label{SY-scalarh3}
\begin{split}
\sup_{\buu_h, \,\brrho_{T,h}} \inf_{\bPphi_h,\buu_h^*, \,\brrho_{T,h}^*}~& \stint{\bHH^*(\buu_h^*)
\\
& \hspace{-0.2cm}+
\sint{\bGG^*(\brrho_{T,h}^*)-\brrho_{T,h}\cdot (\bpphi_h(T,\cdot)+\brrho_{T,h}^*)}\\
\frac{r}{2}\stint{(\underline{\mathcal{D}}(\bPphi_h)-\buu_h^*)\cdot(\underline{\mathcal{D}}(\bPphi_h)-\buu_h^*)}\\
\frac{r}{2}\sint{|\bpphi_h(T,\cdot)+\brrho_{T,h}^*|^2}
\end{split}
\end{equation}

Finally, we introduce a splitting algorithm (modified ALG2) for
the saddle-point system (<ref>), where we sequentially compute each component of
$\bPphi_h$ in the linear elliptic update (Step A)
and of the densities $\brrho_h$ and $\brrho_{T,h}$ in the nonlinear update (Step B1); see also [17].

The coupled global elliptic linear system for $\bPphi_h=(\bpphi_h, \bssigma_h)$ in ALG2 takes the following form:
Find $(\bpphi_h, \bssigma_h)\in [V_h^{k}]^{(d+1)M}$ with
such that, for all
$(\bppsi_h, \bttau_h)\in [V_h^{k}]^{(d+1)M}$ with
$\bttau_h\cdot\nnu|_{\partial\Omega}=0$, there exists

\begin{align}
\label{eqn-Phi}
\begin{split}
\\
\end{split}
\end{align}

This system is solved sequentially in practice to drive down the overall computational cost.
The $i$-th component scalar reaction-diffusion solver for $\pphi_{i,h}^\ell$ at $\ell$-th iteration reads as follows:

\begin{align}
\label{eqn-phi1}
\begin{split}
\\
\\
\partial_t\psi_{i,h}}\\
\\
\textstyle\sum_{j=1}^{i-1}\gamma_{j,p}\pphi_{j,h}^\ell
\end{split}
\end{align}

The $i$-th component $H(\mathrm{div})$-elliptic solver for $\ssigma_{i,h}^\ell$ reads as follows

\begin{align}
\label{eqn-sigma1}
\begin{split}
\sum_{i=1}^M\stint{(\rrho_{i,h}^{*,\ell-1}-\frac{\rrho_{i,h}^{\ell-1}}{r}-\partial_t\pphi_{i,h}^\ell)
\beta\nabla\cdot\ttau_{i,h}
\end{split}
\end{align}

Meanwhile, the nonlinear update on each quadrature point for
the densities $\bm{\uprho}_{ij}^\ell=({\uprho}_{1,ij}^\ell,\cdots, {\uprho}_{M,ij}^\ell)\in \mathbb{R}_+^M$

takes the following form:
find $\bm{\uprho}_{ij}^\ell\in\mathbb{R}_+^M$ that minimizes the function

\begin{align}
\label{eqn-rho1}
\begin{split}
\sum_{\mathsf{m}=1}^M\frac{r^2|\overline{\bm{\mathrm{m}}}_{\mathsf{m},ij}^\ell|^2}{r+V_{1,\mathsf{m}}(\uprho_{\mathsf{m}})}
\end{split}
\end{align}

And the nonlinear update for the terminal densities
$\bm{\uprho}_{T,i}^\ell=({\uprho}_{1,T,i}^\ell,\cdots, {\uprho}_{M,T,i}^\ell)\in \mathbb{R}_+^M$
takes the following form:
find $\bm{\uprho}_{T,i}^\ell\in\mathbb{R}_+^M$ that minimizes the function

\begin{align}
\label{eqn-rho2}
\end{align}

Here the bar-values take the same form as in Step A3 of Algorithm <ref>,

\begin{align}
\label{barX}
\overline{\underline{\uui}}_{ij}^\ell:=(\overline{\bm\uprho}_{ij}^\ell,
\underline{\overline{\bm{\mathrm{m}}}}_{ij}^\ell,
\overline{\bm{\mathrm{s}}}_{ij}^\ell,
\overline{\underline{\bm{\mathrm{n}}}}_{ij}^\ell):=
\underline{\mathcal{D}}\bPphi_{ij}^\ell+\underline{\uui}_{ij}^{\ell}/r,
\;
\overline{\bm\uprho}_{T,i}^\ell:=
\bm{\uprho}_{T,i}^{\ell}/r-
\bm{\pphi}_h^\ell(T,\bm\xi_i).
\end{align}

The pointwise optimization problems in (<ref>) and (<ref>)
are $M$-dimensional problems. These problems are usually loosely coupled among the density components. We further reduce the computational cost by separately solving for each density component sequentially, which results in $M$ one-dimensional minimization problems per integration point. The positivity of the densities is guaranteed in our algorithm at the quadrature points.
After the densities $\bm{\uprho}_{ij}^\ell$ and $\bm{\uprho}_{T,i}^\ell$ have been computed, we then update the other physical variables on the quadrature point as follows:

m_𝗆,ij^ℓ=    rV_1,𝗆(_𝗆,ij^ℓ)/r+V_1,𝗆(_𝗆,ij^ℓ)
n_𝗆,ij^ℓ= rV_3,𝗆(_𝗆,ij^ℓ)/r+V_3,𝗆(_𝗆,ij^ℓ)

m_𝗆,ij^*,ℓ=    r/r+V_1,𝗆(_𝗆,ij^ℓ)
n_𝗆,ij^*,ℓ= r/r+V_3,𝗆(_𝗆,ij^ℓ)

_ij^*,ℓ=    _ij^ℓ-_ij^ℓ/r,
_T,i^*,ℓ= _T,i^ℓ-_T,i^ℓ/r,

s_p,ij^ℓ=    rV_2,p(_ij^ℓ)/r+V_2,p(_ij^ℓ)
s_p,ij^*,ℓ= r/r+V_2,p(_ij^ℓ)

For completeness, we collect the above procedures into the following algorithm.

One iteration of (modified) ALG2  for (<ref>).

[1]
Step A. For $i=1,\cdots, M$, solve the linear system (<ref>) for the unknown $\pphi_{i,h}^\ell$, and solve the linear system (<ref>)
for the unknown $\ssigma_{i,h}^\ell$.
Interpolate the operator $\underline{\mathcal{D}}(\bPphi_h^\ell)$
on quadrature points, and compute the bar-values (on quadrature points) in (<ref>).

Step B. For each quadrature point, sequentially compute the nonnegative density minimizers to (<ref>) and (<ref>).

Step C. Update the other variables on quadrature points according to (<ref>).

§ NUMERICAL RESULTS

In this section, we present one- and two-dimensional numerical results for the scalar scheme (<ref>)
and the system scheme (<ref>) using the finite element software MFEM [2].
The spatial domain $\Omega=[0,1]^d$ is taken to be either a unit line segment ($d=1$) or a unit square ($d=2$). The terminal time is  $T=1$.
Furthermore, we solve the mean field planning problem for all the numerical simulations, where the initial and terminal densities are prescribed. In this case, the terminal density is no longer an  unknown variable. Hence the scalar optimization problem (<ref>) simplifies as below

\begin{equation}\label{SD-scalarh3X}
\begin{split}
\sup_{\uu_h} \inf_{\Pphi_h,\uu_h^*}~& \stint{\HH^*(\uu_h^*)
\\
& \hspace{-0.2cm}+
\sint{-\rrho^1\, \pphi_h(T,\cdot)
\frac{r}{2}\stint{(\mathcal{D}(\Pphi_h)-\uu_h^*)\cdot(\mathcal{D}(\Pphi_h)-\uu_h^*)},
\end{split}
\end{equation}

and the system optimization problem (<ref>) simplifies as below

\begin{equation}\label{SY-scalarh3X}
\begin{split}
\sup_{\buu_h} \inf_{\bPphi_h,\buu_h^*}~& \stint{\bHH^*(\buu_h^*)
\\
& \hspace{-0.2cm}+
\sint{-\brrho^1\cdot \bpphi_h(T,\cdot)+\brrho^0\cdot\bpphi_h(0,\cdot)}\\
\frac{r}{2}\stint{(\underline{\mathcal{D}}(\bPphi_h)-\buu_h^*)\cdot(\underline{\mathcal{D}}(\bPphi_h)-\buu_h^*)}.
\end{split}
\end{equation}

We take the augmented Lagrangian parameter $r=1$ for all the simulations.

§.§ Scalar MFC for reaction-diffusion: $\beta=0$

We first consider a planning problem for  scalar MFC (<ref>) with regularization parameter $\beta=0$ (no Fisher information functional).
With $\beta=0$, we do not compute the $\nn_h, \nn_h^*$ and $\ssigma_h$ variables in the optimization problem (<ref>), since they are always zero and are decoupled from the other variables.

We solve the problem (<ref>) (with $\beta=0$) on
$\Omega_T=[0, 1]\times \Omega$ with $\Omega=[0,1]^d$ for $d=1,2$.
The initial (resp. terminal) densities are chosen to be:

\begin{equation*}
\rho^0(x)=\exp(-50|x-x_A|^2),\quad \rho^1(x)=\exp(-50|x-x_B|^2), \forall x=(x_1,\cdots,x_d)\in\Omega.
\end{equation*}

Here $x_A=0.25$, $x_B=0.75$ for $d=1$, and  $x_A=(0.25, 0.25)$, $x_B=(0.75,0.75)$ for $d=2$.
We fix $V_1(\rho)=\rho$ and potential function

\begin{equation}
\label{ener1}
F(t,x,\rho) = -(0.01\rho\log(\rho)+0.4\rho\cos(4\pi t)\prod_{i=1}^d\cos(4\pi x_i)).
\end{equation}

Moreover, we take the following four choices of $V_2(\rho)$ to highlight different reaction effects of the model:

\begin{align}
\label{V2s}
\begin{cases}
\text{Case 1}: V_2(\rho)=0 \\[.2ex]
\text{Case 2}: V_2(\rho)=20 \\[.2ex]
\text{Case 3}: V_2(\rho)=20\rho \\[.2ex]
\text{Case 4}: V_2(\rho)= 20\frac{\rho-1}{\log(\rho)}.
\end{cases}
\end{align}

The 1D results are computed on a $64\times 64$ uniform rectangular space-time mesh, and
the 2D results are computed on a $16\times 64\times 64$ uniform cubic space-time mesh.
We use polynomial degree $k=4$, and apply 400 ALG iterations for all the simulations.
The density contours for the $d=1$
on the space-time 2D domain $\Omega_T$ are shown in Figure <ref>.
The snapshots of density contours on $\Omega=[0,1]^2$ at
different times for the $d=2$ are shown in
Figure <ref>. We find the dynamics with and without reaction terms are completely different. The reaction function in Figure <ref> (b) and (d) shows that the reaction mobility function $V_2$ dominates the path between $\rho^0$ and $\rho^1$. While the transport mobility function $V_1$ dominates the path between $\rho^0$ and $\rho^1$ in Figure <ref> (a) and (c). In Figure <ref>, we observe a different pattern formulation for various choices of general nonlinear reaction mobility function $V_2$. This behaves very differently from the classical optimal transport problem with $V_2=0$.

[$V_2(\rho) = 0$.]

[$V_2(\rho) = 20$.]

[$V_2(\rho) = 20\rho$.]

[$V_2(\rho) = 20\frac{\rho-1}{\log(\rho)}$.]

Example <ref>.
Snapshots of density contour on the space-time domain $\Omega_T=[0,1]^2$.
The vertical axis represents time.

[$V_2(\rho) = 0$.]

[$V_2(\rho) = 20$.]

[$V_2(\rho) = 20\rho$.]

[$V_2(\rho) = 20\frac{\rho-1}{\log(\rho)}$.]

Example <ref>.
Snapshots of density contour at
$t=$ 0.1,0.3,0.5, 0.7, 0.9 (left to right).

§.§ Scalar MFC for reaction-diffusion: effect of  $\beta$

Here we use a similar setup as in Example <ref>, but study the effect of the regularization parameter $\beta$.
We take take mobility functions $V_1(\rho) = \rho$, $V_2(\rho) = 20$, and potential function

\begin{align}
\label{enerF1}
F(t,x,\rho) = -(0.005\rho\log(\rho)+0.4\rho\cos(4\pi t)\prod_{i=1}^d\cos(4\pi x_i)).
\end{align}

We further take the mobility function
$V_3(\rho) = \rho$, and vary $\beta$, namely,
$\beta=0.005$, and
The same discretization as in Example <ref> is used.
The density contours for $d=1$
on the space-time domain $\Omega_T$ are shown in Figure <ref>.
The snapshots of density contours at
different times for $d=2$ are shown in
Figure <ref>. In numerical examples, we show that increasing regularization $\beta$ leads to a more diffusive density evolution.

Example <ref>.
Snapshots of density contour on  $\Omega_T=[0,1]^2$.

[$V_2(\rho) = 20, \beta=0$. (no regularization)]

[$\beta=0.005$. (weak regularization)]

[$\beta=0.01$. (strong regularization)]

Example <ref>.
Snapshots of density contour at
$t=$ 0.1,0.3,0.5,0.7,0.9 (left to right).

§.§ System MFC for reaction-diffusion ($M=2$, $R=1$)

In this example, we consider a system model (<ref>) with $M=2$ species and $R=1$ reaction.
The initial/terminal densities for each component are given as:

\begin{align*}
\rho_1(0, x) = \exp(-50|x-x_A|^2), \;\;
\rho_1(1, x) = \exp(-50|x-x_B|^2), \\
\rho_2(0,x) = \exp(-50|x-x_B|^2), \;\;
\rho_2(1,x) = \exp(-50|x-x_A|^2).
\end{align*}

We take the mobilities $V_{1,i}(\rho_i)=V_{3,i}(\rho_i)=\rho_i$ for $i=1,2$,
and the mobility

\[
V_{2,1}(\rho_1, \rho_2)=20\frac{\rho_1-\rho_2}{\log(\rho_1)-\log(\rho_2)}.
\]

The potential function is given as follows

\begin{align*}
\bm F(t,x,\bmr) = &\;- \Big(0.01\rho_1\log(\rho_1)
+0.4\rho_1\cos(4\pi t)\prod_{i=1}^d\cos(4\pi x_i)\\
&\;\hspace{1cm} +0.005\rho_2\log(\rho_2)\Big)
\end{align*}

The potential function for the first component has a drift and an entropy term in $\rho_1$, while that for the second component only has a smaller entropy term in $\rho_2$. We use the same discretization as in the previous two examples.
The density contours for the 1D results ($d=1$) are shown in Figure <ref> for  $\beta=0$,
$\beta=0.005$, and
The snapshots of density contours at
different times for the 2D results ($d=2$) are shown in
Figure <ref>
for $\beta=0$ and
$\beta=0.01$. Without reaction mobility $V_2$, the second component density path $\rho_2$ will be similar to a Gaussian translation.
It is clear that the reaction mobility completely changes the density evolution. Moreover, increasing the regularization parameter $\beta$ leads to a diffusive density path, as expected.

[$\rho_1$. Left: $\beta=0$, middle: $\beta=0.005$, right: $\beta=0.01$.]

[$\rho_2$. Left: $\beta=0$, middle: $\beta=0.005$, right: $\beta=0.01$.]

Example <ref>.
Snapshots of $\rho_1$ (top) and $\rho_2$ (bottom).

[$\rho_1$. $\beta=0$. ]

[$\rho_1$. $\beta=0.01$.]

[$\rho_2$. $\beta=0$.]

[$\rho_2$.  $\beta=0.01$. ]

Example <ref>.
Snapshots of $\rho_1$ and $\rho_2$ at
$t=$ 0.1,0.3,0.5,0.7,0.9 (left to right).

§.§ 2D Scalar MFC for reaction-diffusion: image transfer

In this example, we consider a 2D scalar mean field planning problem (<ref>) with complex initial and terminal densities $\rho^0(x)$ and $\rho^1(x)$, and
a complex spatial coordinate dependent mobility $V_2(x, \rho)=20\rho^2(x)\rho$ and potential  $F(x,\rho)= -0.001\rho^3(x)\rho$.
Here the four non-negative functions $\rho^0(x)$, $\rho^1(x)$, $\rho^2(x)$, and $\rho^3(x)$ are normalized mascot images as shown in Figure <ref>, which are logos from University of Notre Dame, UCLA, Portland State university, and University of South Carolina, respectively.

Example <ref>.
The four functions $\rho^0(x)$, $\rho^1(x)$, $\rho^2(x)$, $\rho^3(x)$.
Initial density is $\rho^0$, terminal density is $\rho^1$, mobility $V_2(x,\rho) = 20\rho^2(x)\rho$
and potential  $F(x, \rho)=-0.001\rho^3(x)\rho$.

We further take $V_1(\rho)=V_3(\rho) = \rho$, and vary the regularization parameter $\beta=0$, $\beta=5\times 10^{-4}$, and
We apply the scheme (<ref>) with polynomial degree $k=4$ on a $16\times 64\times 64$ uniform mesh, and perform 2000 ALG iterations.
Snapshots of density contours at
different times for different $\beta$
are shown in Figure <ref>.
It is interesting to observe that the mobility coefficient $\rho^2(x)$ in $V_2$ and the interaction coefficient $\rho^3(x)$ in $F$
are imprinted in the density evolution.
We also observe a strong diffusion effect when $\beta=10^{-3}$.

[$\beta=5\times 10^{-4}$]

Example <ref>.
Snapshots of density contours at
$t=$ 0.2,0.4,0.6,0.8 (left to right) for different $\beta$.

§.§ 2D System MFC for reaction-diffusion ($M=2$, $R=1$): image transfer

In this example, we consider a system model (<ref>) with $M=2$ species and $R=1$ reaction.
We take the initial and terminal densities as the images in Figure <ref>.
Specifically, the two initial densities $\rho_1^0 = \rho^0(x)$, $\rho_2^0=\rho^3(x)$,
and the two terminal densities $\rho_1^1 = \rho^1(x)$ and $\rho_2^1=\rho^2(x)$.
We take $V_{1,i}(\rho)=\rho$,

\[
V_{2,1}(\bmr) = 20\frac{\rho_1+\rho_2}{2},
\]

no potential $\bm F(\bmr) = 0$, and no regularization $\beta=0$.
We use the same discretization as the previous example, and apply 2000 ALG iterations.
Snapshots of the density contour at different times are shown in Figure <ref>.
Here the reaction with mobility $V_{2,1}$ makes the density evolution different from a classical optimal transport path for each component.

Example <ref>.
Snapshots of $\rho_1$ (top) and $\rho_2$ (bottom) at
$t=$ 0.2,0.4,0.6,0.8 (left to right).

§.§ 2D System MFC for reaction-diffusion ($M=12$, $R=12$)

In the last example, we consider a system model (<ref>) with $M=12$ species and $R=12$ reaction.
The initial densities are the ancient Chinese calligraphy in seal script
for the 12 Chinese zodiac animals, which are downloaded from Richard Sears' website <https://hanziyuan.net/>, while the terminal densities are their associated (gray-scale) images, which are generated by
Baidu's text-to-image AI tool WenXin YiGe <https://yige.baidu.com/>. We rescale the images so that the maximal value is 1 and minimial value is 0; see Figure <ref>–<ref> for the 12 initial ($t=0$) and terminal ($t=1$) density approximations in gray scale.

Here we take $V_{1,i}(\rho)=\rho$,

\[
V_{2,p}(\bmr) = 20 \frac{\rho_{p}-\rho_{p+1}}{\log(\rho_{p})-\log(\rho_{p+1})},\quad \forall 1\le p\le 12,
\]

where the convention $\rho_{13}=\rho_1$ is used. The reaction patterns for this system are cyclic. Again, we set potential  $\bm F(\bmr) = 0$, and regularization $\beta=0$.
We use the same discretization as the previous example and apply 2000 ALG iterations.
The results at time $t=0.0, 0.2, 0.5, 0.8$, and $1.0$ are shown in Figure <ref>–<ref>.
The color range is gray-scale from 0 (black) to 1 (white).
We observe interesting and complex densities' evolutions from these figures.

For comparison purposes, we also plot the snapshots of densities at time $t=0.5$ for optimal transport without reaction ($V_{2,p}=0$)
in Figure <ref>. It is observed that the results in the middle row of Figure <ref>–<ref> for the reaction-diffusion system model are very different from the scalar optimal transport results in Figure <ref>.
These differences come from the nonlinear reaction mobility functions.

Example <ref>.
Snapshots of first 6 densities $\rho_1$ to $\rho_6$ (left to right) at times $t=0.0, 0.2, 0.5, 0.8, 1.0$ (top to bottom).

Example <ref>.
Snapshots of last 6 densities $\rho_7$ to $\rho_{12}$ (left to right) at times $t=0.0, 0.2, 0.5, 0.8, 1.0$ (top to bottom).

Example <ref>. (No reaction $V_{2,p} = 0$.)
Snapshots of densities at $t=0.5$. From left to right, top to bottom: $\rho_1$ to $\rho_{12}$.

§ DISCUSSIONS

In this paper, we model and compute a class of generalized optimal transport and mean field control problems for reaction-diffusion equations and systems. The control problems are constructed by general choices of transport and reaction mobility functions, such as $V_1$, $V_2$, and $V_3$, derived from entropy dissipation properties. We apply a high-order spatial-time finite element method to discretize the spatial-time domain and use the ALG2 algorithm to compute mean field control problems. Numerical examples in Section <ref>, including transporting two Gaussian distributions and a system of $12$ images, demonstrate the effectiveness of the proposed mean field control models and computations.

In future work, we shall study general modeling, computation, and inverse mean field control problems for reaction-diffusion systems. The generalized nonlinear reaction functions often represent complex behaviors between different populations, exhibiting patterns in social dynamical systems. The mean field control problem over transportation and reactions provides new patterns in population behaviors observed from our numerical examples. The analytical study of these new patterns could be a future research direction. We also expect that the mean field control problem of reaction-diffusion systems has vast applications in pandemic control, computer vision, and image processing problems. The other important direction is the parallel and high-order computation of generalized optimal transport and mean field control problems on three-dimensional spatial domains.

[1]
L. Ambrosio, N. Gigli, and G. Savaré, Gradient Flows in Metric
Spaces and in the Space of Probability Measures, Lectures in Mathematics
ETH Zürich, Birkhäuser, Basel, 2. ed ed., 2008.

[2]
R. Anderson, J. Andrej, A. Barker, and et al., MFEM: A modular
finite element methods library, Comput. Math. Appl., 81 (2021), pp. 42–74.

[3]
J.-D. Benamou and Y. Brenier, A numerical method for the optimal
time-continuous mass transport problem and related problems, Contemporary
mathematics, 226 (1999), pp. 1–12.

[4]
A. Bensoussan, T. Huang, and M. Laurière, Mean Field Control
and Mean Field Game Models with Several Populations,  (2018).

[5]
S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al., Distributed optimization and statistical learning via the alternating
direction method of multipliers, Foundations and Trends® in
Machine learning, 3 (2011), pp. 1–122.

[6]
P. Cardaliaguet, G. Carlier, and B. Nazaret, Geodesics for a class
of distances in the space of probability measures, Calculus of Variations
and Partial Differential Equations, 48 (2013), pp. 395–420.

[7]
J. Carrillo, S. Lisini, G. Savaré, and D. Slepčev, Nonlinear mobility continuity equations and generalized displacement
convexity, Journal of Functional Analysis, 258 (2010), pp. 1273–1309.

[8]
J. A. Carrillo, K. Craig, L. Wang, and C. Wei, Primal dual methods
for Wasserstein gradient flows, Found. Comput. Math., 22 (2022),
pp. 389–443.

[9]
X. Chen, H. Benveniste, and A. R. Tannenbaum, Unbalanced
Regularized Optimal Mass Transport with Applications to Fluid
Flows in the Brain,  (2023).

[10]
Y. Chen, T. T. Georgiou, and M. Pavon, On the Relation Between
Optimal Transport and Schrödinger Bridges: A Stochastic Control
Viewpoint, Journal of Optimization Theory and Applications, 169 (2016),
pp. 671–691.

[11]
Y. Chen, T. T. Georgiou, and A. Tannenbaum, Vector-Valued Optimal
Mass Transport, SIAM Journal on Applied Mathematics, 78 (2018),
pp. 1682–1696.

[12]
L. Chizat, G. Peyré, B. Schmitzer, and F.-X. Vialard, Unbalanced
optimal transport: dynamic and Kantorovich formulations, J. Funct. Anal.,
274 (2018), pp. 3090–3123.

[13]
G. Conforti, A second order equation for Schrödinger bridges
with applications to the hot gas experiment and entropic transportation
cost, Probability Theory and Related Fields, 174 (2019), pp. 1–47.

[14]
J. Dolbeault, B. Nazaret, and G. Savaré, A new class of
transport distances between measures, Calculus of Variations and Partial
Differential Equations, 34 (2009), pp. 193–231.

[15]
M. Fortin and R. Glowinski, Augmented Lagrangian methods, vol. 15
of Studies in Mathematics and its Applications, North-Holland Publishing Co.,
Amsterdam, 1983.
Applications to the numerical solution of boundary value problems,
Translated from the French by B. Hunt and D. C. Spicer.

[16]
G. Fu, S. Liu, S. Osher, and W. Li, High order computation of
optimal transport, mean field planning, and mean field games,
arXiv:2302.02308 [math.NA],  (2023).

[17]
G. Fu, S. Osher, and W. Li, High order spatial discretization for
variational time implicit schemes: Wasserstein gradient flows and
reaction-diffusion systems, arXiv:2303.08950 [math.NA],  (2023).

[18]
T. O. Gallouët and L. Monsaingeon, A JKO Splitting Scheme
for Kantorovich–Fisher–Rao Gradient Flows, SIAM Journal on
Mathematical Analysis, 49 (2017), pp. 1100–1130.

[19]
Y. Gao, W. Li, and J.-G. Liu, Master equations for finite state mean
field games with nonlinear activations, arXiv:2212.05675 [math.OC],  (2022).

[20]
M. Huang, R. P. Malhamé, P. E. Caines, et al., Large population
stochastic dynamic games: closed-loop mckean-vlasov systems and the nash
certainty equivalence principle, Communications in Information & Systems, 6
(2006), pp. 221–252.

[21]
R. Jordan, D. Kinderlehrer, and F. Otto, The variational formulation
of the fokker–planck equation, SIAM Journal on Mathematical Analysis, 29
(1998), pp. 1–17.

[22]
J.-M. Lasry and P.-L. Lions, Mean field games, Jpn. J. Math., 2
(2007), pp. 229–260.

[23]
W. Lee, R. Lai, W. Li, and S. Osher, Generalized unnormalized
optimal transport and its fast algorithms, arXiv preprint arXiv:2001.11530,

[24]
W. Lee, S. Liu, W. Li, and S. Osher, Mean field control problems for
vaccine distribution, Research in the Mathematical Sciences, 9 (2022),
p. 51.

[25]
W. Lee, S. Liu, H. Tembine, W. Li, and S. Osher, Controlling
Propagation of Epidemics via Mean-Field Control, SIAM Journal on
Applied Mathematics, 81 (2021), pp. 190–207.

[26]
F. Léger and W. Li, Hopf–Cole transformation via
generalized Schrödinger bridge problem, Journal of Differential
Equations, 274 (2021), pp. 788–827.

[27]
W. Li, W. Lee, and S. Osher, Computational mean-field information
dynamics associated with reaction-diffusion equations, J. Comput. Phys., 466
(2022), pp. Paper No. 111409, 30.

[28]
W. Li, S. Liu, and S. Osher, Controlling conservation laws ii:
Compressible navier–stokes equations, Journal of Computational Physics, 463
(2022), p. 111264.

[29]
height 2pt depth -1.6pt width 23pt, Controlling
conservation laws i: entropy-entropy flux, Journal of Computational Physics,
480 (2023).

[30]
W. Li, J. Lu, and L. Wang, Fisher information regularization schemes
for Wasserstein gradient flows, J. Comput. Phys., 416 (2020), pp. 109449,

[31]
W. Li, P. Yin, and S. Osher, Computations of Optimal Transport
Distance with Fisher Information Regularization, Journal of Scientific
Computing, 75 (2018), pp. 1581–1595.

[32]
M. Liero, A. Mielke, and G. Savaré, Optimal
Entropy-Transport problems and a new
Hellinger–Kantorovich distance between positive measures,
Inventiones mathematicae, 211 (2018), pp. 969–1117.

[33]
J. Maas, Gradient flows of the entropy for finite Markov
chains, Journal of Functional Analysis, 261 (2011), pp. 2250–2292.

[34]
A. Mielke, A gradient structure for reaction-diffusion systems and
for energy-drift-diffusion systems, Nonlinearity, 24 (2011), pp. 1329–1346.

[35]
A. Mielke, D. R. M. Renger, and M. A. Peletier, A Generalization
of Onsager's Reciprocity Relations to Gradient Flows with
Nonlinear Mobility, Journal of Non-Equilibrium Thermodynamics, 41
(2016), pp. 141–149.

[36]
L. Monsaingeon, L. Tamanini, and D. Vorotnikov, The dynamical
Schrödinger problem in abstract metric spaces,  (2020).

[37]
L. Onsager and S. Machlup, Fluctuations and Irreversible
Processes, Physical Review, 91 (1953), pp. 1505–1512.

§ APPENDIX

This section gives detailed proofs of Propositions <ref>–<ref>.

Denote $\mathcal{E}(\rho)=\int_\Omega E(\rho(x))dx$. Thus $\frac{\delta}{\delta\rho}\mathcal{E}(\rho)=E'(\rho)$.
By the change of variable formula, we obtain

\begin{equation*}
\bmm=\tilde\bmm-\beta V_1(\rho)\nabla E'(\rho),\quad {s}=\tilde s-\beta V_2(\rho)E'(\rho).
\end{equation*}

The constraint set (<ref>) satisfies

\begin{equation*}
\begin{split}
0=&\partial_t \rho+\nabla\cdot(\tilde\bmm-\beta V_1(\rho)\nabla E'(\rho))-(\tilde s-\beta V_2(\rho)E'(\rho))\\
=&\partial_t \rho+\nabla\cdot \bmm-{s}.
\end{split}
\end{equation*}

Moreover, the terms in objective functional (<ref>) satisfy

\begin{equation*}
\begin{split}
&\int_0^T\int_\Omega \Big[\frac{\|\tilde \bmm(t,x)\|^2}{ 2V_1(\rho(t,x))}+ \frac{|\tilde s(t,x)|^2}{2V_2(\rho(t,x))}\Big] dxdt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm+\beta V_1(\rho)\nabla E'(\rho)\|}{ 2V_1(\rho)}+ \frac{|s+\beta V_2(\rho)E'(\rho)|^2}{2V_2(\rho)}\Big] dxdt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2+\beta^2 V_1(\rho)^2\|\nabla E'(\rho)\|^2}{ 2V_1(\rho)}+ \frac{|s|^2+\beta^2 V_2(\rho)^2|E'(\rho)|^2}{2V_2(\rho)}\Big] dxdt\\
&+\beta\int_0^T\int_\Omega \Big[\bmm\cdot \nabla E'(\rho)-s\cdot E'(\rho)\Big]  dx dt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2}{ 2V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2}{2}\|\nabla E'(\rho)\|^2V_1(\rho)+\frac{\beta^2}{2}|E'(\rho)|^2V_2(\rho)\Big] dxdt\\
&+\beta\int_0^T\int_\Omega E'(\rho)\Big[-\nabla\cdot \bmm+s\Big]  dx dt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2}{ 2V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2}{2}\|\nabla E'(\rho)\|^2V_1(\rho)+\frac{\beta^2}{2}|E'(\rho)|^2V_2(\rho)\Big] dxdt\\
&+\beta\int_0^T\int_\Omega \Big[E'(\rho)\partial_t \rho\Big]  dx dt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2}{ 2V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2}{2}\|\nabla E'(\rho)\|^2V_1(\rho)+\frac{\beta^2}{2}|E'(\rho)|^2V_2(\rho)\Big] dxdt\\
&+\beta\int_\Omega \Big[E(\rho(T,\cdot))-E(\rho^0)\Big]dx,
\end{split}
\end{equation*}

where the last equality follows the fact that $\int_0^TE'(\rho)\partial_t \rho dt=\int_0^T\partial_t E(\rho) dt=E(\rho(T,\cdot))-E(\rho^0)$.
Thus we derive the variational problem (<ref>) using the definition of $V_3$ in (<ref>).

Denote the Lagrange multiplier of problem (<ref>) as $\phi\colon [0, T]\times\Omega\rightarrow\mathbb{R}$. Consider the following saddle point problem:

\begin{equation*}
\begin{split}
\inf_{\bmm, s, \rho, \rho_T}\sup_\phi \quad \mathcal{L}(\bmm, s, \rho, \rho_T, \phi),
\end{split}
\end{equation*}

\begin{equation*}
\begin{split}
\mathcal{L}(\bmm, s,\rho, \phi)=& \int_0^T \int_\Omega \Big[ \frac{\|\bmm\|^2}{2V_1(\rho)}+\frac{|s|^2}{2V_2(\rho)}+\phi\Big(\partial_t \rho+ \nabla\cdot \bmm-s\Big)\Big] dxdt\\
&+\int_{0}^T\Big[\frac{\beta^2}{2}\mathcal{I}(\rho)-\mathcal{F}(\rho)\Big] dt
\end{split}
\end{equation*}

Assume $\rho>0$. By solving the saddle point problem of $\mathcal{L}$, we obtain

\begin{equation*}
\left\{\begin{split}
&\frac{\delta}{\delta \bmm}\mathcal{L}=0,\\
& \frac{\delta}{\delta s}\mathcal{L}=0,\\
&\frac{\delta}{\delta \rho}\mathcal{L}=0,\\
&\frac{\delta}{\delta \phi}\mathcal{L}=0,\\
&\frac{\delta}{\delta \rho_T}\mathcal{L}=0,\\
\end{split}\right.\quad\Rightarrow\quad\left\{\begin{split}
&\frac{\bmm}{V_1}=\nabla \phi,\\
& \frac{s}{V_2}=\phi,\\
&-\frac{1}{2}\frac{\|\bmm\|^2}{V_1^2}V_1'-\frac{1}{2}\frac{|s|^2}{V_2^2}V_2'+\frac{\delta}{\delta \rho}\Big[\frac{\beta^2}{2}\mathcal{I}(\rho)-\mathcal{F}(\rho)\Big]-\partial_t\phi=0,\\
&\partial_t\rho+\nabla\cdot \bmm-s=0,\\
& \phi_T+\frac{\delta}{\delta\rho_T}\Big(\mathcal{G}(\rho_T)+\beta\mathcal{E}(\rho_T)\Big)=0.
\end{split}\right.
\end{equation*}

This finishes the proof.

The proof is similar to the scalar case. Denote $\mathcal{E}_i(\rho_i)=\int_\Omega E_i(\rho_i(x))dx$. Thus $\frac{\delta}{\delta\rho_i}\mathcal{E}_i(\rho_i)=E_i'(\rho_i)$.
By the change of variable formula, we obtain

\begin{equation*}
\bmm_i=\tilde\bmm_i-\beta V_{1,i}(\rho_i)\nabla E'_i(\rho_i),\quad {s}_p=\tilde s_p-\beta V_{2,p}(\bmr)\sum_{j=1}^M\gamma_{j,p}E'_j(\rho_j).
\end{equation*}

The constraint set (<ref>) satisfies

\begin{equation*}
\begin{split}
0=&\partial_t \rho_i+\nabla\cdot(\tilde\bmm_i-\beta V_{1,i}(\rho_i)\nabla E'_i(\rho_i))-(\sum_{p=1}^R\gamma_{i,p}\tilde s_p-\beta\sum_{p=1}^R\gamma_{i,p}V_{2,p}(\bmr)\sum_{j=1}^M\gamma_{j,p}E'_j(\rho_j))\\
=&\partial_t \rho_i+\nabla\cdot \bmm_i-\sum_{p=1}^R\gamma_{i,p}{s}_p.
\end{split}
\end{equation*}

Moreover, the terms in objective functional (<ref>) satisfies

\begin{equation*}
\begin{split}
&\int_0^T\int_\Omega \Big[\sum_{i=1}^M\frac{\|\tilde \bmm_i(t,x)\|^2}{ 2V_{1,i}(\rho_i(t,x))}+ \sum_{p=1}^R\frac{|\tilde s_p(t,x)|^2}{2V_{2,p}(\bmr(t,x))}\Big] dxdt\\
=&\int_0^T\int_\Omega \Big[\sum_{i=1}^M\frac{\| \bmm_i+\beta V_{1,i}(\rho_i)\nabla E'_i(\rho_i)\|}{ 2V_{1,i}(\rho_i)}+ \sum_{p=1}^R\frac{|s_p+\beta \sum_{j=1}^M \gamma_{j,p}V_{2,p}(\bmr)E'_j(\rho_j)|^2}{2V_{2,p}(\bmr)}\Big] dxdt\\
=&\int_0^T\int_\Omega \Big[\sum_{i=1}^M\frac{\| \bmm_i\|^2+\beta^2 V_{1,i}(\rho_i)^2\|\nabla E'_i(\rho_i)\|^2}{ 2V_{1,i}(\rho_i)}+ \sum_{p=1}^R\frac{|s_p|^2+\beta^2 |\sum_{j=1}^M\gamma_{j,p}V_{2,p}(\rho)E'_j(\rho_j)|^2}{2V_{2,p}(\bmr)}\Big] dxdt\\
&+\beta\int_0^T\int_\Omega \Big[\sum_{i=1}^M\bmm_i\cdot \nabla E'_i(\rho_i)-\sum_{p=1}^R s_p\cdot\sum_{j=1}^M\gamma_{j,p}V_{2,p}(\bmr)E'_j(\rho_j)\Big]  dx dt.
\end{split}
\end{equation*}

We only need to show that

\begin{equation*}
\int_0^T\int_\Omega \Big[\sum_{i=1}^M\bmm_i\cdot \nabla E'_i(\rho_i)-\sum_{p=1}^R s_p\cdot\sum_{j=1}^M\gamma_{j,p}E'_j(\rho_j)\Big]  dx dt=\int_\Omega \sum_{i=1}^M\Big[E_i(\rho_i(T,\cdot))-E_i(\rho_i^0)\Big]dx.
\end{equation*}

This is true from the following fact:

\begin{equation*}
\begin{split}
&\int_0^T\int_\Omega \Big[\sum_{i=1}^M\bmm_i\cdot \nabla E'_i(\rho_i)-\sum_{p=1}^R s_p\cdot\sum_{j=1}^M\gamma_{j,p}E'_j(\rho_j)\Big]  dx dt\\
=&\int_0^T \int_\Omega \sum_{i=1}^M  E'_i(\rho_i)\cdot \Big[-\nabla\cdot \bmm_i+\sum_{p=1}^R\gamma_{i,p}s_p\Big] dx dt\\
=&\int_\Omega \int_0^T \sum_{i=1}^M E'_i(\rho_i)\cdot\partial_t\rho_i dt dx\\
\int_\Omega \sum_{i=1}^M\Big[E_i(\rho_i(T,\cdot))-E_i(\rho_i^0)\Big]dx,
% \int_0^T \int_\Omega \Big[G(\bmr_T)-G(\bmr_0)\Big]dxdt,
\end{split}
\end{equation*}

where the last equality uses the integration by parts in the time variable. This finishes the proof.

The proof is similar to the scalar case. We derive the minimizer system for variational problem (<ref>). Denote $\phi_i\colon [0, T]\times \Omega\rightarrow \mathbb{R}$ as the Lagrange multiplier of constraint (<ref>), for $i=1,2,\cdots, M$. Write $\phi=(\phi_i)_{i=1}^M$. Consider the following saddle point problem

\begin{equation*}
\inf_{\mm,s,\rho, \rho_T}\sup_{\phi}\quad \mathcal{L}_1(\mm, s, \bmr, \bmr_T, \phi),
\end{equation*}

\begin{equation*}
\begin{split}
&\mathcal{L}_{1}(\mm, s, \bmr, \bmr_T, \phi)\\
=&\int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^M
\frac{\|\bmm_i\|^2}{V_{1,i}(\rho_i)}
% (\m_i, V_{1,i}(\rho_i)^{-1}\m_i)
&+\mathcal{G}(\bmr_T)+\beta \Big(\mathcal{E}(\bmr_T)-\mathcal{E}(\bmr_0)\Big)\\
&+\sum_{i=1}^M\int_0^T\int_\Omega \phi_i\cdot\Big\{\partial_t \rho_i + \nabla\cdot \m_i-\sum_{p=1}^R\gamma_{i,p}\s_{p}\Big\} dxdt.
\end{split}
\end{equation*}

Assume $\rho_i>0$, $i=1,2,\cdots, M$. By solving the saddle point problem of $\mathcal{L}_1$, we obtain

\begin{equation*}
\left\{\begin{split}
&\frac{\delta}{\delta \bmm_{i}}\mathcal{L}_1=0,\\
& \frac{\delta}{\delta \s_{p}}\mathcal{L}_1=0,\\
&\frac{\delta}{\delta \rho_i}\mathcal{L}_1=0,\\
&\frac{\delta}{\delta \phi_i}\mathcal{L}_1=0,\\
&\frac{\delta}{\delta \bmr_T}\mathcal{L}_1=0,
\end{split}\right.\quad\Rightarrow\quad\left\{\begin{split}
&{\bmm_i}=V_{1,i}(\rho_i)\nabla \phi_i,\\
& \s_{p}=V_{2,p}(\bmr)\sum_{j=1}^M \gamma_{j,p}\phi_p,\\
% -\frac{1}{2}(V_{1,i}(\rho_i)^{-1}\bmm_i, \frac{\partial}{\partial \rho_i}V_{1,i}(\rho_i) V_{1,i}(\rho_i)^{-1}\bmm_i)\\
&-\frac{1}{2}\sum_{p=1}^R\frac{|s_{p}|^2}{V_{2,p}(\bmr)^2}\frac{\partial}{\partial \rho_i}V_{2,p}(\bmr)\\
&+\frac{\partial}{\partial \rho_i}[\frac{\beta^2}{2}\bm{\mathcal{I}}(\bmr)-\mathcal{F}(\bmr)]-\partial_t \phi_i=0,\\
& \partial_t \rho_i + \nabla\cdot \bmm_i-\sum_{p=1}^R\gamma_{i,p}s_{p}=0,\\
&\phi_{T}+\frac{\delta}{\delta \bmr_{T}}\mathcal{G}(\bmr_T)+\beta \frac{\delta}{\delta\bmr_{T}}\mathcal{E}(\bmr_T)=0.
\end{split}\right.
\end{equation*}

This finishes the proof.

§.§ Case 2

In this subsection, we consider the other example of reaction-diffusion systems. We assume that equation (<ref>) satisfies a gradient flow in generalized optimal transport spaces; see details in [Onsager, Miekle]. Denote $C^2$ strictly convex functions $G_i\colon \mathbb{R}\rightarrow\mathbb{R}$, $i=1,2,\cdots, n$. Write

\begin{equation*}
\end{equation*}

Suppose equation (<ref>) can be reformulated as follows:

\begin{equation*}
\begin{split}
\partial_t u_i=&\nabla\cdot(V_{1,i}(u)\nabla  E'_i(u_i))-\sum_{j=1, j\neq i}^nV_{2,ij}(u)(G'_i(u_i)-G'_j(u_j)).
\end{split}
\end{equation*}

In the above formula, we assume that the following scalar mobility function $V_{1, i}(u)\in \mathbb{R}^{1}$ holds:

\begin{equation*}
\end{equation*}

for any $i=1, 2, \cdots, n$. And there exists a mobility function $V_{2,ij}(u)\in\mathbb{R}^{n\times n}$, with

\begin{equation*}
V_{2,ij}(u)=V_{2,ji}(u)\geq 0,
\end{equation*}

\begin{equation*}
\sum_{j=1, j\neq i}^nV_{2,ij}(u)(G_i'(u_i)-G_j'(u_j))=-R_i(u),
\end{equation*}

for any $i,j=1, 2, \cdots, n$.

As in the scalar equation case, we observe the following entropy dissipation property holds:

\begin{equation*}
\begin{split}
&\frac{d}{dt}\int_\Omega G(u(t,x))dx=\sum_{i=1}^n\int_\Omega (G'_i(u_i), \partial_t u_i) dx \\
=&\sum_{i=1}^n\int_\Omega G'_i(u_i)\Big(\nabla\cdot(V_{1,i}(u)\nabla  G'_i(u_i)-\sum_{j=1, j\neq i}^nV_{2,ij}(u)(G_i'(u_i)-G_j'(u_j))\Big)dx\\
=&-\int_\Omega \sum_{i=1}^n(\nabla G'_i(u_i), V_{1,i}(u) \nabla G'_i(u_i))+\frac{1}{2}\sum_{i,j=1}^n|G'_i(u_i)- G'_j(u_j)|^2V_{2,ij}(u)dx\\
:=&-\mathcal{I}(u)\leq 0,
\end{split}
\end{equation*}

where we apply the integration by parts in the second last equality.

Similarly, we denote $\mathcal{I}(u)$ as a generalized Fisher information functional.
We can derive a vector density space metric from the entropy dissipation property. It introduces the following optimal control problems for reaction-diffusion system (<ref>).

Define a terminal time $T>0$, an energy functional $\mathcal{F}\colon\mathcal{M}\rightarrow\mathbb{R}$, and a terminal functional $\mathcal{G}\colon\mathcal{M}\rightarrow\mathbb{R}$. Consider

\begin{equation}
\begin{split}
&\inf_{v_1, v_2, u}~~\int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^n(v_{1,i}(t,x), V_{1,i}(u(t,x))v_{1,i}(t,x))\\
&\hspace{2.5cm}+\frac{1}{4}\sum_{i,j=1, i\neq j}^n|v_{2,ij}(t,x)|^2V_{2,ij}(u(t,x))-\mathcal{F}(u(t,\cdot))\Big]dxdt+\mathcal{G}(u_T),
\end{split}
\end{equation}

where the infimum is among all densities $u_i\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, vector fields $v_{1,i}\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and anti-symmetric reaction rates function $v_{2,ij}\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, $v_{2,ij}=-v_{2,ji}$, for $i,j=1,2,\cdots, n$, such that

\begin{equation}\label{SOCII}
\begin{split}
&\partial_t u_i(t,x) + \nabla\cdot( V_{1,i}(u(t,x)) v_{1,i}(t,x))-\sum_{j=1}^nV_{2,ij}(u(t,x))v_{2,ij}(t,x)\\
=&\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x))\Big),
\end{split}
\end{equation}

with fixed initial densities $u_i(0,x)=u_{i0}(x)$.

As in the scalar case, there are two reformulations of optimal control problems of reaction-diffusion systems (<ref>).

Denote $m_i\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s_{ij}\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
m_i(t,x)=V_{1,i}(u(t,x))v_{1,i}(t,x), \quad s_{ij}(t,x)=V_{2,ij}(u(t,x))v_{2,ij}(t,x).
\end{equation*}

\begin{equation}\label{SOCRI}
\begin{split}
&\inf_{m, s, u} \quad \int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^n\frac{\|m_i(t,x)\|^2}{V_{1,i}(u(t,x))}+\frac{1}{4}\sum_{i,j=1, i\neq j}^n\frac{|s_{ij}(t,x)|^2}{V_{2,ij}(u(t,x))}\\
\end{split}
\end{equation}

where the infimum is among all densities $u_i\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, flux functions $m_i\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and anti-symmetric reaction functions $s_{ij}\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, $s_{ij}=-s_{ji}$ for $i,j=1,2,\cdots, n$, such that

\begin{equation}\label{SOCII}
\begin{split}
&\partial_t u_i(t,x) + \nabla\cdot m_i(t,x)-\sum_{j=1}^ns_{ij}(t,x)=\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x))\Big),
\end{split}
\end{equation}

with fixed initial densities $u_i(0,x)=u_{i0}(x)$.

For $i=1,2,\cdots, n$, denote $\bmm_i\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s_{ij}\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
\bmm_i=m_i-\beta V_{1,i}(u)\nabla G'_i(u_i),\quad s_{ij}=s_{ij}-\beta V_{2,ij}(u)(G_j'(u_j)-G_i'(u_i)).
\end{equation*}

\begin{equation*}
\end{equation*}

\begin{equation}
\begin{split}
&\inf_{\bmm, s, u} ~~ \int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^n\frac{\|\bmm_i(t,x)\|^2}{ V_{1,i}(u(t,x))}+\frac{1}{4}\sum_{i,j=1, i\neq j}^n\frac{|s_{ij}(t,x)|^2}{V_{2,ij}(u(t,x))}+\frac{\beta^2}{2}\sum_{i=1}^n\frac{|\nabla u_i(t,x)|^2}{V_{3,i}(u(t,x))}\\
&\hspace{2.0cm}+\frac{\beta^2}{4}\sum_{i,j=1, i\neq j}^n{|G_i'(u_i(t,x))-G_j'(u_j(t,x))|^2}{V_{2,ij}(u(t,x))}-\mathcal{F}(u(t,\cdot))\Big]dxdt\\
&\hspace{1cm}+\mathcal{G}(u_T)+\beta\int_\Omega \sum_{i=1}^n\Big[G_i(u_i(T,x))-G_i(u_i(0,x))\Big]dx,
\end{split}
\end{equation}

where the infimum is among all densities $u_i\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, $\bmm_i\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and anti-symmetric reaction functions $s_{ij}\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, $s_{ij}=-s_{ji}$, for $i,j =1,2,\cdots,n$, such that

\begin{equation}\label{SOCII}
\partial_t u_i(t,x) + \nabla\cdot \bmm_i(t,x)-\sum_{j=1}^ns_{ij}(t,x)=0,
\end{equation}

with fixed initial densities $u_i(0,x)=u_{i0}(x)$.

We also derive a proposition that presents the critical point system of variational problems (<ref>), (<ref>), and (<ref>). We demonstrate the equivalence among these optimal control problems of reaction-diffusion systems and their critical point systems.

The critical point of variation problem (<ref>) has the following two formulations. Let $i=1,2,\cdots, n$.

(i) There exists a function $\Phi_i\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that

\begin{equation*}
\begin{split}
\end{split}
\end{equation*}

And $(u, \Phi)$ satisfies the following system. When $t\in (0, T)$

\begin{equation}\label{SSBP}
\left\{\begin{aligned}
&\partial_tu_i(t,x) +\nabla\cdot(V_{1,i}(u(t,x))\nabla\Phi_i(t,x))\\
&\hspace{1.5cm}-\sum_{j=1,j\neq i}^nV_{2,ij}(u(t,x))(\Phi_i(t,x)-\Phi_j(t,x))\\
=&\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x))\Big),\\
&\partial_t\Phi_i(t,x)+\frac{1}{2}\sum_{j=1}^n(\nabla\Phi_j(t,x),  \frac{\partial}{\partial u_i}V_{1,j}(u(t,x))\nabla\Phi_j(t,x))\\
&\hspace{1.7cm}+\frac{1}{4}\sum_{j,k=1}^n|\Phi_j(t,x)-\Phi_k(t,x)|^2\frac{\partial}{\partial u_i}V_{2,kj}(u(t,x))+\frac{\delta}{\delta u_i}\mathcal{F}(u)\\
=&-\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla \Phi_i)-\sum_{j=1}^n(\nabla \Phi_j, \frac{\delta}{\delta u_i}A_j(u(t,x))\nabla u_j)\\
&\hspace{1cm}+\sum_{j=1}^n\Phi_j(t,x)\frac{\partial}{\partial u_i}R_j(u(t,x))\Big).
\end{aligned}\right.
\end{equation}

And both initial and terminal time boundary conditions satisfy

\begin{equation}
u_i(0,x)=u_{i0}(x),\quad \Phi_i(T,x)=-\frac{\delta}{\delta u_{iT}}\mathcal{G}(u_{iT}).
\end{equation}

(ii) There exists a function $S_i\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that

\begin{equation*}
\begin{split}
&\frac{\bmm_i(t,x)}{V_{1,i}(u(t,x))}=\nabla S_i(t,x),\\
\end{split}
\end{equation*}

And $(u, S)$ satisfies the following system. When $t\in (0, T)$,

\begin{equation}\label{SSBPI}
\left\{\begin{aligned}
&\partial_tu_i(t,x) +\nabla\cdot(V_{1,i}(u(t,x))\nabla S_i(t,x))\\
&\hspace{1.5cm}-\sum_{j=1,j\neq i}^nV_{2,j}(u(t,x))(S_i(t,x)-S_j(t,x))=0,\\
&\partial_tS_i(t,x)+\frac{1}{2}\sum_{j=1}^n(\nabla S_j(t,x),\frac{\partial}{\partial u_i} V_1(u(t,x)) \nabla S_j(x))\\
&\hspace{1.6cm}+\frac{1}{4}\sum_{j,k=1}^n|S_j(t,x)-S_k(t,x)|^2\frac{\partial}{\partial u_i}V_{2,jk}(u(t,x))\\
&\hspace{1.6cm}+\frac{\delta}{\delta u_i}\Big[\mathcal{F}(u)-\frac{\beta^2}{2}\mathcal{I}(u)\Big]=0,
\end{aligned}\right.
\end{equation}

where $\mathcal{I}(u)=\int_\Omega \sum_{i=1}^n(\nabla G'_i(u_i),V_{1,i}(u)\nabla G'_i(u_i))+\frac{1}{2}\sum_{i,j=1}^n|G_i'(u_i)-G_j'(u_j)|^2V_{2,ij}(u) dx$ is the generalized Fisher information functional. And both initial and terminal time boundary conditions satisfy

\begin{equation}
u_i(0,x)=u_{i0}(x),\quad S_i(T,x)=-\frac{\delta}{\delta u_{iT}(x)}\mathcal{G}(u_T)-\beta G_i'(u_i(T,x)).
\end{equation}

\begin{equation*}
S_i(t,x)=\Phi_i(t,x)-\beta G'_i(u_i(t,x)).
\end{equation*}

Then equation system (<ref>) is equivalent to equation system (<ref>).

The proof follows from similar derivations in Proposition <ref>. We still present the proof for the completeness of this paper. Variational problem (<ref>) is equivalent to problem (<ref>). Denote $\Phi_i\colon [0, T]\times \Omega\rightarrow \mathbb{R}$ as the Lagrange multiplier of constraint (<ref>), for $i=1,2,\cdots, n$. Consider the following saddle point problem

\begin{equation*}
\inf_{m,s,u, u_T}\sup_{\Phi}\quad \mathcal{L}_{32}(m, s, u, u_T, \Phi),
\end{equation*}

\begin{equation*}
\begin{split}
&\mathcal{L}_{32}(m, s, u, u_T, \Phi)\\
=&\int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^n(m_i, V_{1,i}(u)^{-1}m_i)+\frac{1}{4}\sum_{i,j=1, i\neq j}^n\frac{|s_{ij}|^2}{V_{2,ij}(u)}-\mathcal{F}(u(t,\cdot))\Big]dxdt+\mathcal{G}(u_T)\\
&+\sum_{i=1}^n\int_0^T\int_\Omega \Phi_i\cdot\Big\{\partial_t u_i + \nabla\cdot m_i-\sum_{j=1}^ns_{ij}-\beta\Big(\nabla\cdot(A_i(u)\nabla u_i)+R_i(u)\Big)\Big\} dxdt.
\end{split}
\end{equation*}

Again, assume $u_i>0$, $i=1,2,\cdots, n$. By solving the saddle point problem of $\mathcal{L}_1$, we obtain

\begin{equation*}
\left\{\begin{split}
&\frac{\delta}{\delta m_{i}}\mathcal{L}_{32}=0,\\
& \frac{\delta}{\delta s_{ij}}\mathcal{L}_{32}=0,\\
&\frac{\delta}{\delta u_i}\mathcal{L}_{32}=0,\\
&\frac{\delta}{\delta \Phi_i}\mathcal{L}_{32}=0,\\
&\frac{\delta}{\delta u_T}\mathcal{L}_{32}=0,
\end{split}\right.\quad\Rightarrow\quad\left\{\begin{split}
& s_{ij}=V_{2,ij}(u)(\Phi_i-\Phi_j),\\
&-\frac{1}{2}(V_{1,i}(u)^{-1}m_i, \frac{\partial}{\partial u_i}V_{1,i}(u) V_{1,i}(u)^{-1}m_i)\\
&-\frac{1}{4}\sum_{j,k=1, j\neq k}^n\frac{|s_{jk}|^2}{V_{2,jk}(u)^2}\frac{\partial}{\partial u_i}V_{2,jk}(u)-\frac{\partial}{\partial u_i}\mathcal{F}(u(t,\cdot))-\partial_t\Phi_i\\
&-\beta\nabla\cdot(A_i(u)\nabla\Phi_i)+\beta \sum_{j=1}^n(\nabla \Phi_j, \frac{\delta}{\delta u_i}A_j(u(t,x))\nabla u_j)\\
&-\beta \sum_{j=1}^n\frac{\partial}{\partial u_j}R_i(u)=0,\\
& \partial_t u_i + \nabla\cdot m_i-\sum_{j=1}^ns_{ij}-\beta\Big(\nabla\cdot(A_i(u)\nabla u_i)+R_i(u)\Big)=0,\\
&\Phi_T+\frac{\delta}{\delta u_T}\mathcal{G}(u_T)=0.
\end{split}\right.
\end{equation*}

This finishes the proof of (i).

We next prove (ii). We derive the minimizer system for variational problem (<ref>). Denote $S_i\colon [0, T]\times \Omega\rightarrow \mathbb{R}$ as the Lagrange multiplier of constraint (<ref>), for $i=1,2,\cdots, n$. Consider the following saddle point problem

\begin{equation*}
\inf_{\m,\s,u, u_T}\sup_{S}\quad \mathcal{L}_{42}(\m, \s, u, u_T, S),
\end{equation*}

\begin{equation*}
\begin{split}
&\mathcal{L}_{42}(\m, \s, u, u_T, S)\\
=&\int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^n(\m_i, V_{1,i}(u)^{-1}\m_i)+\frac{1}{4}\sum_{i,j=1, i\neq j}^n\frac{|\s_{ij}|^2}{V_{2,ij}(u)}+\frac{\beta^2}{2}\mathcal{I}(u)-\mathcal{F}(u(t,\cdot))\Big]dxdt\\
&+\mathcal{G}(u_T)+\beta \int_\Omega (G(u_T)-G(u_0))dx\\
&+\sum_{i=1}^n\int_0^T\int_\Omega S_i\cdot\Big\{\partial_t u_i + \nabla\cdot \m_i-\sum_{j=1}^n\s_{ij}\Big\} dxdt.
\end{split}
\end{equation*}

Assume $u_i>0$, $i=1,2,\cdots, n$. By solving the saddle point problem of $\mathcal{L}_1$, we obtain

\begin{equation*}
\left\{\begin{split}
&\frac{\delta}{\delta \m_{i}}\mathcal{L}_{42}=0,\\
& \frac{\delta}{\delta \s_{ij}}\mathcal{L}_{42}=0,\\
&\frac{\delta}{\delta u_i}\mathcal{L}_{42}=0,\\
&\frac{\delta}{\delta S_i}\mathcal{L}_{42}=0,\\
&\frac{\delta}{\delta u_T}\mathcal{L}_{42}=0,
\end{split}\right.\quad\Rightarrow\quad\left\{\begin{split}
&{\m_i}=V_{1,i}(u)\nabla S_i,\\
& \s_{ij}=V_{2,ij}(u)(S_i-S_j),\\
&-\frac{1}{2}(V_{1,i}(u)^{-1}\m_i, \frac{\partial}{\partial u_i}V_{1,i}(u) V_{1,i}(u)^{-1}\m_i)\\
&-\frac{1}{4}\sum_{j,k=1, j\neq k}^n\frac{|\s_{jk}|^2}{V_{2,jk}(u)^2}\frac{\partial}{\partial u_i}V_{2,jk}(u)\\
&+\frac{\partial}{\partial u_i}[\frac{\beta^2}{2}\mathcal{I}(u)-\mathcal{F}(u)]-\partial_t S_i=0,\\
& \partial_t u_i + \nabla\cdot m_i-\sum_{j=1}^ns_{ij}=0,\\
&S_{iT}+\frac{\delta}{\delta u_{iT}}\mathcal{G}(u_T)+\beta G'_i(u_{iT})=0.
\end{split}\right.
\end{equation*}

This finishes the proof of (ii).

We can prove (iii) by similar derivations in Proposition <ref>.

This section presents the main formulation of mean field control problems for scalar reaction-diffusion equations and systems.

§.§ Scalar reaction-diffusion equations

We first consider a scalar nonlinear reaction-diffusion equation

\begin{equation}\label{gheat}
\partial_t u(t,x)=\Delta F(u(t,x))+R(u(t,x)),
\end{equation}

where $x\in \Omega\subset\mathbb{R}^d$, $\Omega$ is a compact convex set, $u\colon [0,\infty)\times\Omega\rightarrow\mathbb{R}_+$ is a positive density function, $\Delta$ is the Euclidean Laplacian operator, and $F$, $R\colon \mathbb{R}_+\rightarrow\mathbb{R}$ are smooth functions. We assume the periodic boundary condition on $\Omega$. Denote

\begin{equation*}
\mathcal{M}=\{u\in C^{\infty}(\Omega)\colon u(x)\geq 0\}.
\end{equation*}

We briefly review that equation (<ref>) satisfies a gradient flow in generalized optimal transport spaces; see details in [LLO, M]. Denote a $C^2$ strictly convex function $G\colon \mathbb{R}\rightarrow\mathbb{R}$.
Then equation (<ref>) can be reformulated as follows:

\begin{equation*}
\begin{split}
\partial_t u
=&\nabla\cdot(F'(u)\nabla u)+R(u)\\
=&\nabla\cdot(V_1(\rho)\nabla G'(\rho))-V_2(\rho)G'(\rho).
\end{split}
\end{equation*}

where we use the fact $\nabla G'(\rho)=G''(u)\nabla u$ in the second equality, and define mobility functions

\begin{equation}\label{mob}
V_1(\rho):=\frac{F'(u)}{G''(u)}, \quad  V_2(\rho):=-\frac{R(u)}{G'(\rho)}.
\end{equation}

We observe the following entropy dissipation property holds:

\begin{equation*}
\begin{split}
\frac{d}{dt}\int_\Omega G(u(t,x))dx =&\int_\Omega G'(\rho)\partial_t u dx \\
=&\int_\Omega G'(\rho) \Big[\nabla\cdot(V_1(\rho)\nabla G'(\rho)-V_2(\rho)G'(\rho)\Big]dx\\
=&-\int_\Omega \|\nabla G'(\rho)\|^2V_1(\rho)+|G'(\rho)|^2V_2(\rho)\Big]dx\\
\end{split}
\end{equation*}

We call the term $\mathcal{I}(u)$ generalized Fisher information functional. The above equality illustrates that the dissipation of entropy functional equals the generalized Fisher information functional. In literature, it is known that the entropy dissipation property induces an infinite dimensional metric in $\mathcal{M}$ and optimal control problems for reaction-diffusion equations. See [AMS, LO] for detailed studies among metric spaces, entropy dissipation properties, and optimal control problems.

In this subsection, we directly study optimal control problems for scalar reaction-diffusion equation (<ref>).

Define mobility functions

\begin{equation*}
V_1(\rho)=\frac{F'(u)}{G''(u)}, \quad  V_2(\rho)=-\frac{R(u)}{G'(\rho)}.
\end{equation*}

Define an energy functional $\mathcal{F}\colon\mathcal{M}\rightarrow\mathbb{R}$, a terminal time $T>0$, and a terminal functional $\mathcal{G}\colon\mathcal{M}\rightarrow\mathbb{R}$. Consider

\begin{equation}
\begin{split}
&\inf_{v_1, v_2, u}~~\int_0^T\int_\Omega \Big[\frac{1}{2}\|v_1(t,x)\|^2 V_1(u(t,x))+ \frac{1}{2}|v_2(t,x)|^2V_2(u(t,x))-\mathcal{F}(u(t,\cdot))\Big]dxdt\\
\end{split}
\end{equation}

where the infimum is among all densities $u\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, vector fields $v_1\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and reaction rates function $v_2\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, such that

\begin{equation}\label{OCII}
\partial_t u(t,x) + \nabla\cdot( V_1(u(t,x)) v_1(t,x))-V_2(u(t,x))v_2(t,x)=\beta\Delta F(u(t,x))+\beta R(u(t,x)),
\end{equation}

with fixed initial density $u(0,x)=u_0(x)$.

We note that variational problem (<ref>) is a generalized optimal control problem motivated by the Benamou-Brenier formula in optimal transport.
E.g., let $V_1(\rho)=u$, $V_2(\rho)=0$, $\beta=0$. Then variational problem (<ref>) forms the classical optimal transport and mean field control problem. In general modeling perspectives, variational problem (<ref>) models that the density moves from the initial density under the dynamical constraint of reaction-diffusion equations. It aims to find the optimal choices of vector fields and reaction rate functions under general kinetic, potential, and terminal energies.

There are two equivalent reformulations of variational problem (<ref>). Both of them are useful in designing numerical schemes in later sections.

Denote $m\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
m(t,x)=V_1(u(t,x))v_1(t,x), \quad s(t,x)=V_2(u(t,x))v_2(t,x).
\end{equation*}

\begin{equation}\label{OCRI}
\begin{split}
&\inf_{m, s, u} \quad \int_0^T\int_\Omega \Big[\frac{\|m(t,x)\|^2}{ 2V_1(u(t,x))}+ \frac{|s(t,x)|^2}{2V_2(u(t,x))}\Big]dx-\mathcal{F}(u(t,\cdot))dt+\mathcal{G}(u_T),
\end{split}
\end{equation}

where the infimum is among all densities $u\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, $m\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and $s\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, such that

\begin{equation}\label{OCII}
\partial_t u(t,x) + \nabla\cdot m(t,x)-s(t,x)=\beta\Delta F(u(t,x))+\beta R(u(t,x)),
\end{equation}

with fixed initial density $u_0$.

Denote $\bmm\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
\bmm=m-\beta V_1(\rho)\nabla G'(\rho),\quad s=s-\beta V_2(\rho)G'(\rho).
\end{equation*}

\begin{equation*}
\end{equation*}

\begin{equation}
\begin{split}
&\inf_{\bmm, s, u} ~~ \int_0^T\int_\Omega \Big[\frac{\|\bmm(t,x)\|^2}{2 V_1(u(t,x))}+ \frac{|s(t,x)|^2}{2V_2(u(t,x))}\\
&\hspace{2.0cm}+\beta^2\frac{|\nabla u(t,x)|^2}{2V_3(u(t,x))}+\frac{\beta^2}{2}{|G'(u(t,x))|^2}{V_2(u(t,x))}\Big]dx-\mathcal{F}(u(t,\cdot))dt\\
&\hspace{1cm}+\mathcal{G}(u_T)+\beta\int_\Omega \Big[G(u(T,x))-G(u(0,x))\Big]dx,
\end{split}
\end{equation}

where the infimum is among all densities $u\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, $\bmm\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and $s\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, such that

\begin{equation}\label{OCII}
\partial_t u(t,x) + \nabla\cdot \bmm(t,x)-s(t,x)=0,
\end{equation}

with fixed initial density $u(0,x)=u_0(x)$.

We next present a proposition, which presents the critical point system of variational problems (<ref>), (<ref>), and (<ref>). We demonstrate the equivalence among these optimal control problems and their critical point systems.

The critical point of variation problem (<ref>) has the following two formulations.

(i) There exists a function $\Phi\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that

\begin{equation*}
v_1(t,x)=\frac{m(t,x)}{V_1(u(t,x))}=\nabla\Phi(t,x),\quad v_2(t,x)=\frac{s(t,x)}{V_2(u(t,x))}=\Phi(t,x).
\end{equation*}

And $(u, \Phi)$ satisfies, when $t\in (0, T)$

\begin{equation}\label{SBP}
\left\{\begin{aligned}
&\partial_tu(t,x) +\nabla\cdot(V_1(u(t,x))\nabla\Phi(t,x))-V_2(u(t,x))\Phi(t,x)\\
=&\beta\Big(\Delta F(u(t,x))+R(u(t,x))\Big),\\
&\partial_t\Phi(t,x)+\frac{1}{2}\|\nabla\Phi(t,x)\|^2 V_1'(u(t,x))+\frac{1}{2}|\Phi(t,x)|^2V'_2(u(t,x))+\frac{\delta}{\delta u}\mathcal{F}(u)\\
\end{aligned}\right.
\end{equation}

And initial and terminal time boundary conditions satisfy

\begin{equation}
u(0,x)=u_0(x),\quad \Phi(T,x)=-\frac{\delta}{\delta u_T}\mathcal{G}(u_T).
\end{equation}

(ii) There exists a function $S\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that $(u, S)$ satisfies

\begin{equation*}
\frac{\bmm(t,x)}{V_1(u(t,x))}=\nabla S(t,x),\quad \frac{s(t,x)}{V_2(u(t,x))}=S(t,x),
\end{equation*}

\begin{equation}\label{SBPI}
\left\{\begin{aligned}
&\partial_tu(t,x) +\nabla\cdot(V_1(u(t,x))\nabla S(t,x))-V_2(u(t,x))S(t,x)=0,\\
&\partial_tS(t,x)+\frac{1}{2}\|\nabla S(t,x)\|^2 V_1'(u(t,x))+\frac{1}{2}|S(t,x)|^2V'_2(u(t,x))\\
&\hspace{6.3cm}+\frac{\delta}{\delta u}\Big[\mathcal{F}(u)-\frac{\beta^2}{2}\mathcal{I}(u)\Big]=0,
\end{aligned}\right.
\end{equation}

where $\mathcal{I}(u)=\int_\Omega \|\nabla G'(\rho)\|^2V_1(\rho)+|G'(\rho)|^2V_2(\rho) dx$ is the generalized Fisher information functional. And initial and terminal time boundary conditions satisfy

\begin{equation}
u(0,x)=u_0(x),\quad S(T,x)=-\frac{\delta}{\delta u_T(x)}\Big(\mathcal{G}(u_T)+\beta \int_\Omega G(u(T,x))dx\Big).
\end{equation}

\begin{equation*}
S(t,x)=\Phi(t,x)-\beta G'(u(t,x)).
\end{equation*}

Then equation system (<ref>) is equivalent to equation system (<ref>).

We first prove (i). Denote $m=V_1v_1$ and $s=V_2v_2$. Then the variational problem (<ref>) forms problem (<ref>). Denote the Lagrange multiplier of problem (<ref>) as $\Phi\colon [0, T]\times\Omega\rightarrow\mathbb{R}$. Consider the following saddle point problem:

\begin{equation*}
\begin{split}
\inf_{m, s, u}\sup_\Phi \quad \mathcal{L}_1(m, s, u,\Phi),
\end{split}
\end{equation*}

\begin{equation*}
\begin{split}
\mathcal{L}_1(m, s,u,\Phi)=& \int_0^T \int_\Omega \Big\{ \frac{\|m\|^2}{2V_1(\rho)}+\frac{|s|^2}{2V_2(\rho)}-\mathcal{F}(u)\\
&\hspace{1cm}+\Phi\Big(\partial_t u+ \nabla\cdot m-s-\beta(\Delta F(u)+R(u))\Big)\Big\} dx dt.
\end{split}
\end{equation*}

Assume $u>0$. By solving the saddle point problem of $\mathcal{L}_1$, we obtain

\begin{equation*}
\left\{\begin{split}
&\frac{\delta}{\delta m}\mathcal{L}_1=0,\\
& \frac{\delta}{\delta s}\mathcal{L}_1=0,\\
&\frac{\delta}{\delta u}\mathcal{L}_1=0,\\
&\frac{\delta}{\delta \Phi}\mathcal{L}_1=0,
\end{split}\right.\quad\Rightarrow\quad\left\{\begin{split}
& \frac{s}{V_2}=\Phi,\\
&-\frac{1}{2}\frac{\|m\|^2}{V_1^2}V_1'-\frac{1}{2}\frac{|s|^2}{V_2^2}V_2'-\frac{\delta}{\delta u}\mathcal{F}(u),\\
&-\partial_t\Phi-\beta\Delta\Phi F'(u)-\beta\Phi R'(u)=0,\\
&\partial_tu+\nabla\cdot m-s=\beta(\Delta F(u)+R(u)).
\end{split}\right.
\end{equation*}

This finishes the proof of (i).

We next prove (ii). From the change of variable formula, we obtain

\begin{equation*}
\bmm=m-\beta V_1(\rho)\nabla G'(\rho),\quad s=s-\beta V_2(\rho)G'(\rho).
\end{equation*}

The constraint set (<ref>) satisfies

\begin{equation*}
\begin{split}
0=&\partial_t u+\nabla\cdot(m-\beta V_1(\rho)\nabla G'(\rho))-(s-\beta V_2(\rho)G'(\rho))\\
=&\partial_t u+\nabla\cdot \bmm-\beta s.
\end{split}
\end{equation*}

Moreover, the terms in objective functional (<ref>) satisfies

\begin{equation*}
\begin{split}
&\int_0^T\int_\Omega \Big[\frac{\|m(t,x)\|^2}{ 2V_1(u(t,x))}+ \frac{|s(t,x)|^2}{2V_2(u(t,x))}\Big] dxdt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm+\beta V_1(\rho)\nabla G'(\rho)\|}{ 2V_1(\rho)}+ \frac{|s+\beta V_2(\rho)G'(\rho)|^2}{2V_2(\rho)}\Big] dxdt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2+\beta^2 V_1(\rho)^2\|\nabla G'(\rho)\|^2}{ 2V_1(\rho)}+ \frac{|s|^2+\beta^2 V_2(\rho)^2|G'(\rho)|^2}{2V_2(\rho)}\Big] dxdt\\
&+\beta\int_0^T\int_\Omega \Big[\bmm\cdot \nabla G'(\rho)-s\cdot G'(\rho)\Big]  dx dt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2}{ 2V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2}{2}\|\nabla G'(\rho)\|^2V_1(\rho)+\frac{\beta^2}{2}|G'(\rho)|^2V_2(\rho)\Big] dxdt\\
&+\beta\int_0^T\int_\Omega G'(\rho)\Big[-\nabla\cdot \bmm+s\Big]  dx dt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2}{ 2V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2}{2}\|\nabla G'(\rho)\|^2V_1(\rho)+\frac{\beta^2}{2}|G'(\rho)|^2V_2(\rho)\Big] dxdt\\
&+\beta\int_0^T\int_\Omega \Big[G'(\rho)\partial_t u\Big]  dx dt\\
=&\int_0^T\int_\Omega \Big[\frac{\| \bmm\|^2}{ 2V_1(\rho)}+ \frac{|s|^2}{2V_2(\rho)}+\frac{\beta^2}{2}\|\nabla G'(\rho)\|^2V_1(\rho)+\frac{\beta^2}{2}|G'(\rho)|^2V_2(\rho)\Big] dxdt\\
&+\beta\int_0^T \Big[G(u_T)-G(u_0)\Big]dx,
\end{split}
\end{equation*}

where the last equality follows the fact that $\int_0^TG'(\rho)\partial_t u dt=\int_0^T\partial_t G(u) dt=G(u_T)-G(u_0)$.
Thus we derive the variational problem (<ref>).

Again, denote the Lagrange multiplier of problem (<ref>) as $S\colon [0, T]\times\Omega\rightarrow\mathbb{R}$. Consider the following saddle point problem:

\begin{equation*}
\begin{split}
\inf_{\bmm, s, u}\sup_S \quad \mathcal{L}_2(\bmm, s, u, S),
\end{split}
\end{equation*}

\begin{equation*}
\mathcal{L}_2(\bmm, s,u, S)= \int_0^T \int_\Omega \Big\{ \frac{\|\bmm\|^2}{2V_1(\rho)}+\frac{|s|^2}{2V_2(\rho)}+S\Big(\partial_t u+ \nabla\cdot \bmm-s\Big)\Big\} dx+\frac{\beta^2}{2}\mathcal{I}(u)-\mathcal{F}(u) dt.
\end{equation*}

Assume $u>0$. By solving the saddle point problem of $\mathcal{L}_2$, we obtain

\begin{equation*}
\left\{\begin{split}
&\frac{\delta}{\delta \bmm}\mathcal{L}_2=0,\\
& \frac{\delta}{\delta s}\mathcal{L}_2=0,\\
&\frac{\delta}{\delta u}\mathcal{L}_2=0,\\
&\frac{\delta}{\delta S}\mathcal{L}_2=0,
\end{split}\right.\quad\Rightarrow\quad\left\{\begin{split}
&\frac{\bmm}{V_1}=\nabla S,\\
& \frac{s}{V_2}=S,\\
&-\frac{1}{2}\frac{\|\bmm\|^2}{V_1^2}V_1'-\frac{1}{2}\frac{|s|^2}{V_2^2}V_2'+\frac{\delta}{\delta u}\Big[\frac{\beta^2}{2}\mathcal{I}(u)-\mathcal{F}(u)\Big]-\partial_tS=0,\\
&\partial_tu+\nabla\cdot \bmm-s=0.
\end{split}\right.
\end{equation*}

This finishes the proof of (i).

We last prove (iii). We first rewrite the dynamics (<ref>) in the Hamiltonian formalism:

\begin{equation*}
\partial_tu=\frac{\delta}{\delta \Phi}\mathcal{\tilde H}(u, \Phi), \quad \partial_t\Phi=-\frac{\delta}{\delta u}\mathcal{\tilde H}(u,\Phi),
\end{equation*}

where $\frac{\delta}{\delta}$ is the $L^2$ first variation operator, and $\tilde H$ is the Hamiltonian functional defined as

\begin{equation*}
\mathcal{\tilde H}(u, \Phi)=\int_\Omega \Big(\frac{1}{2}\|\nabla\Phi\|^2V_1(\rho)+\frac{1}{2}|\Phi|^2V_2(\rho)+\beta\Phi(\Delta F(u)+R(u))\Big)dx+\mathcal{F}(u).
\end{equation*}

Denote $S=\Phi-\beta G'(\rho)$. Then

\begin{equation*}
\begin{split}
&\mathcal{\tilde H}(u, \Phi)\\=&\int_\Omega \Big(\frac{1}{2}\|\nabla\Phi\|^2V_1(\rho)+\frac{1}{2}|\Phi|^2V_2(\rho)-\beta(\nabla\Phi, \nabla G'(\rho))V_1(\rho)-\beta(\Phi, G'(\rho))V_2(\rho)\Big)dx+\mathcal{F}(u)\\
=&\frac{1}{2}\int_\Omega \Big(\|\nabla\Phi-\beta\nabla G'(\rho)\|^2V_1(\rho)+|\Phi-\beta G'(\rho)|^2V_2(\rho)\\
&\hspace{2cm}-\beta^2\|\nabla G'(\rho)\|^2V_1(\rho)-\beta^2|G'(\rho)|^2V_2(\rho)\Big)dx+\mathcal{F}(u)\\
=&\frac{1}{2}\int_\Omega \Big(\|\nabla S\|^2V_1(\rho)+|S|^2V_2(\rho)-\beta^2\mathcal{I}(u)\Big)dx+\mathcal{F}(u)\\
:=&\mathcal{\tilde H}_1(u, S).
\end{split}
\end{equation*}

We simply check that $(u, S)$ satisfies the equation

\begin{equation*}
\partial_tu=\frac{\delta}{\delta S}\mathcal{\tilde H}_1(u, S),\quad \partial_t S=-\frac{\delta}{\delta u}\mathcal{\tilde H}_1(u, S),
\end{equation*}

This finishes the proof.

§.§ Reaction-diffusion systems

We next consider a system of $n$ number nonlinear reaction-diffusion equations:

\begin{equation}\label{gsys}
\partial_t u_i(t,x)=\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x)),
\end{equation}

where $i=1,2,\cdots, n$, $A_i(u)\in \mathbb{R}_+$ is a nonlinear positive scalar, and $R_i\colon \mathbb{R}^n \rightarrow\mathbb{R}$ is a reaction function. There are several different choices of reaction functions. We demonstrate them below. We design mean field control problems for reaction-diffusion systems for each type of reaction function.

We use a chemical reaction-diffusion system as a motivation example. Consider $n=M$ different chemical species
$X_1,\dots, X_M$ reacting according to $R$ mass-action laws:

\begin{align}
\label{mal}
\alpha_1^pX_1+\cdots+\alpha_M^pX_M\xrightleftharpoons[k_-^p]{k_+^p}
\beta_1^pX_1+\cdots+\beta_M^pX_M,
\end{align}

where $p=1, \cdots, R$, $R$ represents the total number of reactions, $\bma^p=(\alpha_1^p,\cdots,\alpha_M^p), \bmb^p=(\beta_1^p,\cdots,\beta_M^p)\in\mathbb{N}_0^M$ are the vectors of the stoichiometric coefficients, and
$k_+^p, k_-^p$ are the positive forward and backward reaction rates. From now on, we assume that $k_+^p=k_-^p=k^p>0$.

Consider the following equation:

\begin{align}
\label{rds}
\partial_t u_i= \gamma_i \Delta u_i -\sum_{p=1}^Rk^p(\alpha_i^p-\beta^p_i)(\bu^{\bma^p}-\bu^{\bmb^p}),
\end{align}

for $1\le i\le M$, where $\gamma_i=A_i(u)\geq 0$,
$\bu=(\rho_1, \cdots, \rho_M)$ and
the multi-index notation
$\bu^{\bma^p}:=\prod_{i=1}^M\bu_i^{\alpha_i^p}$ is used.

We next recast the above system (<ref>) into a gradient flow formulation, following appropriate mobility functions.
As in [34], consider

\begin{align}
\label{mob}
\ell(x,y)=\begin{cases}
\frac{x-y}{\log(x)-\log(y)}&\text{ for } x\not=y,\\[2ex]
y&\text{ for } x=y,
\end{cases}
\end{align}

and denote the following mobility functions:

\begin{align}
\label{mobs-V1}
V_{1, i}(u_i)= &\; \gamma_i u_i, \hspace{2.3cm}\forall 1\le i\le M, \\
\label{mobs-V2}
V_{2,p}(\bmr) = &\;k^p\,\ell\left(\bu^{\bma^p}, \bu^{\bmb^p}\right),\quad
\forall 1\le p\le R.
\end{align}

Denote an entropy function

\begin{equation*}
\end{equation*}

with $G_i(u_i):=u_i\log u_i-u_i$.
Using these notations, it can be shown that (<ref>) is equivalent to

\begin{align}
\label{rds-V1V2}
\partial_tu_i =&\;\nabla\cdot\left(
V_{1,i}(u_i)\nabla G_i'(u_i)\right),
\nonumber \\
&\; -\sum_{p=1}^RV_{2,p}(\bu)(\alpha_i^p-\beta^p_i)\sum_{j=1}^M(\alpha_j^p-\beta_j^p)G_j'(u_j).
\end{align}

where $G_i'(u_i)=\log u_i$, and $V_{1,i}(u_i)\nabla G_i'(u_i)=\gamma_iu_i\nabla\log u_i=\gamma_i\nabla u_i$. The formulation (<ref>) is a gradient flow in a generalized optimal transport type metric space.
Thus the following entropy dissipation property holds:

\begin{equation}
\label{ener-law2}
\begin{split}
\frac{d}{dt}\sum_{i=1}^M\int_\Omega G_i(u_i)dx
=&\;-\sum_{i=1}^M\int_\Omega\|\nabla G_i'(u_i)\|^2V_{1,i}(u_i)\,dx \\
:=& I(u),
\end{split}
\end{equation}

We also denote $\mathcal{I}(u)$ as a generalized Fisher information functional. As in the scalar case, we can derive a vector density space metric from the entropy dissipation property. This further introduces the following optimal control problems for reaction-diffusion system (<ref>).

Define a terminal time $T>0$, an energy functional $\mathcal{F}\colon\mathcal{M}\rightarrow\mathbb{R}$, and a terminal functional $\mathcal{G}\colon\mathcal{M}\rightarrow\mathbb{R}$. Consider

\begin{equation}
\begin{split}
&\inf_{v_1, v_2, u}~~\int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^M(v_{1,i}(t,x), V_{1,i}(u(t,x))v_{1,i}(t,x))\\
\end{split}
\end{equation}

where the infimum is among all densities $u_i\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, vector fields $v_{1,i}\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and anti-symmetric reaction rates function $v_{2,p}\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, $i=1,\cdots, M$, $p=1,\cdots, R$, such that

\begin{equation}\label{SOCIIcaseI}
\begin{split}
&\partial_t u_i(t,x) + \nabla\cdot( V_{1,i}(u(t,x)) v_{1,i}(t,x))-\sum_{p=1}^R(\alpha_j^p-\beta_j^p)V_{2,p}(u(t,x))v_{2,p}(t,x)\\
=&\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x))\Big),
\end{split}
\end{equation}

with fixed initial densities $u_i(0,x)=u_{i0}(x)$.

There are two reformulations of optimal control problems of reaction-diffusion systems (<ref>).

Denote $m_i\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s_{p}\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
m_i(t,x)=V_{1,i}(u(t,x))v_{1,i}(t,x), \quad s_{p}(t,x)=V_{2,p}(u(t,x))v_{2,p}(t,x).
\end{equation*}

\begin{equation}\label{SOCRI}
\begin{split}
&\inf_{m, s, u} \quad \int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^M\frac{\|m_i(t,x)\|^2}{V_{1,i}(u(t,x))}+\frac{1}{2}\sum_{p=1}^R\frac{|s_{p}(t,x)|^2}{V_{2,p}(u(t,x))}\Big]dx\\
&\hspace{3cm}-\mathcal{F}(u(t,\cdot)) dt+\mathcal{G}(u_T),
\end{split}
\end{equation}

where the infimum is among all densities $u_i\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, flux functions $m_i\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and reaction functions $s_{p}\colon [0, T]\times \Omega \rightarrow\mathbb{R}$ for $i=1,2,\cdots, M$, $p=1,\cdots, R$, such that

\begin{equation}\label{SOCII}
\begin{split}
&\partial_t u_i(t,x) + \nabla\cdot m_i(t,x)-\sum_{p=1}^R (\alpha_i^p-\beta_i^p)s_{p}(t,x)=\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x))\Big),
\end{split}
\end{equation}

with fixed initial densities $u_i(0,x)=u_{i0}(x)$.

For $i=1,2,\cdots, M$ and $p=1,\cdots, R$, denote $\bmm_i\colon [0,T]\times\Omega\rightarrow\mathbb{R}^d$,
and $s_{p}\colon [0,T]\times\Omega\rightarrow\mathbb{R}$, such that

\begin{equation*}
\bmm_i=m_i-\beta V_{1,i}(u)\nabla G'_i(u_i),\quad s_{p}=s_{p}-\beta V_{2,p}(u)\sum_{j=1}^M (\alpha_j^p-\beta_j^p) G_j'(u_j).
\end{equation*}

\begin{equation*}
\end{equation*}

\begin{equation}
\begin{split}
&\inf_{\bmm, s, u} ~~ \int_0^T\int_\Omega \Big[\frac{1}{2}\sum_{i=1}^n\frac{\|\bmm_i(t,x)\|^2}{ V_{1,i}(u(t,x))}+\frac{1}{2}\sum_{p=1}^R\frac{|s_{p}(t,x)|^2}{V_{2,p}(u(t,x))}+\frac{\beta^2}{2}\sum_{i=1}^n\frac{|\nabla u_i(t,x)|^2}{V_{3,i}(u(t,x))}\\
&\hspace{2.0cm}+\frac{\beta^2}{2}\sum_{p=1}^R{|\sum_{j=1}^M(\alpha_j^p-\beta_j^p) G_j'(u_j)|^2}{V_{2,p}(u(t,x))}\Big]dx-\mathcal{F}(u(t,\cdot))dt\\
&\hspace{1cm}+\mathcal{G}(u_T)+\beta\int_\Omega \sum_{p=1}^R\Big[G_i(u_i(T,x))-G_i(u_i(0,x))\Big]dx,
\end{split}
\end{equation}

where the infimum is among all densities $u_i\colon [0, T]\times \Omega\rightarrow\mathbb{R}_+$, $\bmm_i\colon [0, T]\times \Omega \rightarrow\mathbb{R}^d$, and reaction functions $s_{p}\colon [0, T]\times \Omega \rightarrow\mathbb{R}$, such that

\begin{equation}\label{SOCII}
\partial_t u_i(t,x) + \nabla\cdot \bmm_i(t,x)-\sum_{p=1}^R(\alpha_i^p-\beta_i^p)s_{p}(t,x)=0,
\end{equation}

with fixed initial densities $u_i(0,x)=u_{i0}(x)$.

We also derive a proposition that presents the critical point system of variational problems (<ref>), (<ref>), and (<ref>). We demonstrate the equivalence among these optimal control problems of reaction-diffusion systems and their critical point systems.

The critical point of variation problem (<ref>) has the following two formulations. Let $i=1,2,\cdots, n$, $p=1,\cdots, R$.

(i) There exists a function $\Phi_i\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that

\begin{equation*}
\begin{split}
&v_{2,p}(t,x)=\frac{s_{p}(t,x)}{V_{2,p}(u(t,x))}=\sum_{j=1}^R (\alpha_j^p-\beta_j^p)\Phi_j(u(t,x)).
\end{split}
\end{equation*}

And $(u, \Phi)$ satisfies the following system. When $t\in (0, T)$,

\begin{equation}\label{SSBPcaseI}
\left\{\begin{aligned}
&\partial_tu_i(t,x) +\nabla\cdot(V_{1,i}(u(t,x))\nabla\Phi_i(t,x))\\
=&\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla u_i(t,x))+R_i(u(t,x))\Big),\\
&\partial_t\Phi_i(t,x)+\frac{1}{2}\sum_{j=1}^n(\nabla\Phi_j(t,x),  \frac{\partial}{\partial u_i}V_{1,j}(u(t,x))\nabla\Phi_j(t,x))\\
&\hspace{1.7cm}+\frac{1}{2}\sum_{p=1}^R|\sum_{j=1}^M(\alpha_j^p-\beta_j^p)\Phi_j(t,x)|^2\frac{\partial}{\partial u_i}V_{2,p}(u(t,x))+\frac{\delta}{\delta u_i}\mathcal{F}(u)\\
=&-\beta\Big(\nabla\cdot(A_i(u(t,x))\nabla \Phi_i)-\sum_{j=1}^n(\nabla \Phi_j, \frac{\delta}{\delta u_i}A_j(u(t,x))\nabla u_j)\\
&\hspace{1cm}+\sum_{j=1}^n\Phi_j(t,x)\frac{\partial}{\partial u_i}R_j(u(t,x))\Big).
\end{aligned}\right.
\end{equation}

And both initial and terminal time boundary conditions satisfy

\begin{equation}
u_i(0,x)=u_{i0}(x),\quad \Phi_i(T,x)=-\frac{\delta}{\delta u_{iT}}\mathcal{G}(u_{iT}).
\end{equation}

(ii) There exists a function $S_i\colon [0, T]\times\Omega\rightarrow \mathbb{R}$, such that

\begin{equation*}
\begin{split}
&\frac{\bmm_i(t,x)}{V_{1,i}(u(t,x))}=\nabla S_i(t,x),\\
\end{split}
\end{equation*}

And $(u, S)$ satisfies the following system. When $t\in (0, T)$,

\begin{equation}\label{SSBPIcaseI}
\left\{\begin{aligned}
&\partial_tu_i(t,x) +\nabla\cdot(V_{1,i}(u(t,x))\nabla S_i(t,x))\\
&\partial_tS_i(t,x)+\frac{1}{2}\sum_{j=1}^n(\nabla S_j(t,x),\frac{\partial}{\partial u_i} V_1(u(t,x)) \nabla S_j(x))\\
&\hspace{1.6cm}+\frac{1}{2}\sum_{p=1}^R|\sum_{j=1}^M (\alpha_j^p-\beta_j^p) S_j|^2\frac{\partial}{\partial u_i}V_{2,p}(u(t,x))\\
&\hspace{1.6cm}+\frac{\delta}{\delta u_i}\Big[\mathcal{F}(u)-\frac{\beta^2}{2}\mathcal{I}(u)\Big]=0,
\end{aligned}\right.
\end{equation}

$$\mathcal{I}(u)=\int_\Omega \sum_{i=1}^n(\nabla G'_i(u_i),V_{1,i}(u)\nabla G'_i(u_i))+\frac{1}{2}\sum_{p=1}^R|(\alpha_j^p-\beta_j^p)G'_j(u)|^2V_{2,p}(u) dx$$
is the generalized Fisher information functional defined in entropy dissipation property (<ref>). And both initial and terminal time boundary conditions satisfy

\begin{equation}
u_i(0,x)=u_{i0}(x),\quad S_i(T,x)=-\frac{\delta}{\delta u_{iT}(x)}\mathcal{G}(u_T)-\beta G_i'(u_i(T,x)).
\end{equation}

\begin{equation*}
S_i(t,x)=\Phi_i(t,x)-\beta G'_i(u_i(t,x)).
\end{equation*}

Then equation system (<ref>) is equivalent to equation system (<ref>).