The global existence and boundedness of solutions to quasi-linear reaction-diffusion systems are investigated. The system arises from compartmental models describing the spread of infectious diseases proposed in [30, 31], where the diffusion rate is assumed to depend on the total population, leading to quasilinear diffusion with possible degeneracy. The mathematical analysis of this model has been addressed recently in [1] where it was essentially assumed that all sub-populations diffuse at the same rate, which yields a positive lower bound of the total population, thus removing the degeneracy. In this work, we remove this assumption completely and show the global existence and boundedness of solutions by exploiting a recently developed $L^p$-energy method. Our approach is applicable to a larger class of systems and is sufficiently robust to allow model variants and different boundary conditions.

Classification AMS 2010: 35A01, 35K57, 35K59, 35Q92

Keywords: Quasi-linear reaction-diffusion systems;  Global existence; $L^p$-energy methods; Mass control.

§ INTRODUCTION

§.§ Mathematical analysis of quasilinear SEIRD models

Modelling the spatial spread of infectious disease is a classical topic and has recently attracted much more interest from the mathematical community due to the Covid-19 epidemic. Classical models include the susceptible-infected-removed (SIR) system and its variants, such as susceptible-exposed-infected-removed (SEIR) or susceptible-exposed-infected-removed-deceased (SEIRD). To account for spatial heterogeneity, a compartmental model was proposed in [30, 31] where the diffusion coefficients are highly heterogeneous and depend on the total population density, leading to a quasi-linear reaction-diffusion system. The simulations therein showed a strong qualitative agreement from the forecast using the model and collected data in Lombardy.

Let $\Omega \subset \mathbb R^N$, $N\ge 1$, be a bounded domain[The domain $\Omega$ was considered in $\R^2$ and $\R^3$ in [30] and [1], respectively. For the current work, we can deal with the case of arbitrary dimension.] with smooth boundary $\partial\Omega$ such that $\Omega$ lies locally on one side of $\pa\Omega$. Let $s(x,t)$, $e(x,t)$, $i(x,t)$, $r(x,t)$, and $d(x,t)$ denote the susceptible, exposed, infected, recovered, and deceased population densities, respectively, at spatial position $x\in\Omega$ and at time $t>0$, and let $n(x,t):=s(x,t)+e(x,t)+i(x,t)+r(x,t)$ be the total living population density. The proposed model in [30] reads as

\begin{equation}\label{sys-m}
\begin{cases}
\pa_t s= \alpha n - (1 - A_0/n) \beta_i si - (1 - A_0/n) \beta_e se - \mu s + \nabla\cdot(n \nu_s \nabla s)\\
\pa_t e= (1 - A_0/n) \beta_i si + (1 - A_0/n) \beta_e se - \sigma e - \phi_e e - \mu e + \nabla\cdot(n \nu_e \nabla e)\\
\pa_t i = \sigma e - \phi_d i - \phi_r i - \mu i + \nabla\cdot(n \nu_i \nabla i)\\
\pa_t r = \phi_r i + \phi_e e - \mu r + \nabla\cdot(n \nu_r \nabla r)\\
\pa_t d = \phi_d i \,,
\end{cases}
\end{equation}

where the positive diffusion coefficients $\nu_s, \nu_e, \nu_i, \nu_r$, the birth rate $\alpha$, the inverse of the incubation period $\sigma$, the asymptomatic recovery rate $\phi_e$, the infected recovery rate $\phi_r$, the infected mortality rate $\phi_d$, the asymptomatic contact rate $\beta_e$, the symptomatic contact rate $\beta_i$, and the general mortality rate $\mu$ are positive and may depend on time and space.

The mathematical analysis of (<ref>) was first studied in [1]. One distinguishing characteristic of this system, in comparison with many other variants of SIR models in the literature, is that the diffusion of all species depends on the total population density $n$, which makes (<ref>) a quasilinear reaction-diffusion system. This also brings possible degeneracy to the diffusion operators and makes the problem more challenging. The global existence and boundedness of solutions to (<ref>) was shown in [1] under the following assumptions and modifications:

(i) all diffusion rates are the same[In fact, [1] generalised the diffusion rate to $\kappa(n)$ for a continuous function $\kappa(\cdot)$, but it does not alter the fact that the diffusion rates must be the same so that the technique therein can be applied.], i.e. $\nu_s = \nu_e = \nu_i = \nu_r = \nu$,

(ii) the terms $(1-A_0/n)$ and $\phi_di$ are replaced by a non-singular function $A(n)$, for instance $A(n) = (1-A_0/n)_+$, and $\phi_dni$, respectively.

The replacement of $(1-A_0/n)$ avoids the singularity when $n$ gets close to zero. A closer examination later reveals that this is not necessary since one can estimate
$(si)/n \le s$ and $(se)/n\le e$. The assumptions of the same diffusion rate and replacement of $\phi_di$ are, on the other hand, essential as they help to obtain a key property that the total population density $n$ is bounded pointwise from below by a positive constant for all time, provided the initial density $n_0$ is also bounded from below. Indeed, summing the equations of $s$, $e$, $i$, and $r$, keeping in mind the aforementioned modifications, yields an equation for $n$ of the form

\begin{equation*}
\partial_tn - \nabla\cdot(\nu n\nabla n) = -\phi_din + (\alpha-\mu)n \red{\geq (\alpha-\mu)n-\phi_d n^2}
\end{equation*}

from which the lower bound for $n$ on finite time intervals follows from the lower bound of initial data. With this lower bound of $n$, the diffusion operators become non-degenerate and the analysis can be carried out in a standard way. This strategy used in [1] seems to break down   when the diffusion rates are different, but it still applies when the term $-\phi_d i$ stays in place, since in this case

\begin{equation*}
\partial_tn - \nabla\cdot(\nu n\nabla n) = -\phi_di + (\alpha-\mu)n \geq (\alpha-\mu-\phi_d)n.
\end{equation*}

So $n$ is bounded from below on finite time intervals if $n(x,0)\geq\alpha>0$ for all $x\in\Omega$.

In this paper, we set out to remove the assumptions (i) and (ii) above and show the global existence and boundedness of solutions to the original system (<ref>). In fact, we show the results for a much larger class of quasilinear reaction-diffusion systems which contains (<ref>) as a special case. Our key idea is to exploit a recently developed $L^p$-energy approach in e.g. [22, 12] which does not require any lower bound on $n$ except for its natural nonnegativity. In the next subsection, we provide the general setting while the main results and ideas are presented in subsection <ref>.

§.§ Problem setting

Let $1\leq N\in\mathbb{N}$, and $\Omega\subset\mathbb{R}^N$ be a bounded domain with Lipschitz boundary $\pa \Omega$. Let $2\leq m\in\mathbb{N}$. In this paper, we study the global existence and boundedness of the following quasi-linear reaction-diffusion system of concentrations $u=(u_1, \cdots, u_m)$, for any $i\in\{1,\cdots, m\}$,

\begin{equation}\label{system}
\begin{cases}
\pa_t u_i-\nabla\cdot(D_i(x,t)\Phi(u)\nabla u_i)=f_i(x,t,u),  &x\in \Omega, t>0\\
(D_i(x,t)\Phi(u)\nabla u_i)\cdot \eta=0,   &x\in \pa \Omega, t>0\\
u_i(x,0)=u_{i,0}(x),   &x\in \Omega,
\end{cases}
\end{equation}

where $\eta$ is the unit outward normal vector on $\pa\Omega$, initial date $u_{i,0}$ are bounded and non-negative,  the diffusion matrix $D_i$: $\Omega\times [0,\infty)\rightarrow\mathbb{R}^{N\times N}$ satisfies

\begin{equation}\label{D:cond-1}
\lambda |\xi|^2\leq \xi^\top D_i(x,t)\, \xi, \quad \forall (x,t)\in \Omega\times[0,\infty), ~\forall \xi\in\mathbb{R}^N, ~\forall i=1,\cdots,m
\end{equation}

for some $\lambda>0$, and for each $T>0$,

\begin{equation}\label{D:cond-2}
D_{i}\in L^{\infty}(\Omega\times(0,T); \R^{N\times N}), \quad\forall i=1,\cdots,m,
\end{equation}

and $\Phi(u)$ satisfies:

*  $\Phi(u): \R^m_+\rightarrow \R_+$ is continuous;

*  There is some $b\geq0$ and $M>0$, such that

\begin{equation*}
\Phi(u)\geq M \,u_i^b, \quad \forall  u \in \mathbb R_+^m \text{ and } i=1,\cdots,m;
\end{equation*}

*  There exist $\pi>0$ and $\widetilde{M}>0$ such that

\begin{equation*}
\Phi(u) \leq \widetilde{M} \bra{1+\sumi u_i^{\pi}}, \quad \forall  u\in \mathbb R_+^m.
\end{equation*}

The nonlinearities $f_i(x, t, u): \Omega\times\mathbb{R}_+\times\mathbb{R}^m_+\rightarrow\mathbb{R}$ satisfy the following conditions:

*  For any $i=1,\cdots,m$ and any $(x,t)\in\Omega\times\mathbb{R}_+$,  $f_i(x,t,\cdot): \mathbb{R}^m\rightarrow \mathbb{R}$  is locally Lipschitz continuous uniformly in $(x,t)\in\Omega\times(0,T)$ for any $T>0$;

*  For any $i=1,\cdots,m$ and any  $(x,t)\in\Omega\times\mathbb{R}_+$,  $f_i(x,t,\cdot)$ is quasi-positive, i.e., $f_i(x,t,u)\geq0$ for all $u\in \R ^m_+$ with $u_i=0$ for all $i=1,\cdots,m$;

*  There exists $c_1,\cdots,c_m>0$ and $K_1,K_2\in\R$ such that

\begin{equation*}
\sumi c_if_i(x,t,u)\leq K_1\sumi u_i+K_2, \quad \forall (x,t,u)\in\Omega\times\R_+\times\R^m_+;
\end{equation*}

*  There exist $K_3>0, r>0$, and a lower triangular matrix $A=(a_{ij})$ with positive diagonal entries, and nonnegative entries otherwise, such that, for any $i=1,\cdots,m$

\begin{equation*}
\sum^i_{j=1}a_{ij}f_j(x,t,u)\leq K_3 \bra{1+\sumi u_i^r}, \quad \forall (x,t,u)\in\Omega\times\R_+\times\R^m_+
\end{equation*}

(we call this assumption an intermediate sum of order $r$);

*  The nonlinearities are bounded above by a polynomial, i.e., there exists $l>0$ and $K_4>0$ such that

\begin{equation*}
f_i(x,t,u)\leq K_4 \bra{1+\sumi u_i^l},  \quad \forall (x,t,u)\in\Omega\times\R_+\times\R^m_+, ~\forall i=1,\cdots, m.
\end{equation*}

The local Lipschitz continuity (<ref>) of the nonlinearities implies the existence of a local solution to (<ref>) on a maximal interval $[0, T_{\max})$.  The quasi-positivity assumption (<ref>) assures that the solution to (<ref>) is non-negative (as long as it exists) if the initial data is non-negative,  which is a natural assumption since we consider here $u_i$ as concentrations or densities.  The assumption (<ref>) gives an upper bound on the total mass of the system.
Reaction-diffusion systems satisfying (<ref>), (<ref>) and (<ref>) appear naturally in modeling many real life phenomena, ranging from chemistry, biology, ecology, or social sciences. Remarkably, these natural assumptions are not enough to ensure global existence of bounded solutions as it was pointed out by counterexamples in [25] and [26] even for semilinear systems, i.e. $\Phi(u)\equiv 1$ and $D_i(x,t)\equiv D_i$ in (<ref>) where $D_i\ne D_j$ for some $i\ne j$. The study of global existence of semilinear systems, i.e. (<ref>) with $\Phi(u)\equiv 1$, has made considerable progress in the last decade, see e.g. [4, 8, 15, 17, 10, 11, 12] and the survey [24].

One can readily check that all assumptions (<ref>)–(<ref>), (<ref>)–(<ref>) are fulfilled for the SEIR system (<ref>) (excluding the equation of $d$ as it is uncoupled) provided that the diffusion rates are bounded from below by positive constants, and all other rates are nonnegative and bounded functions of $x$ and $t$. Indeed, by writing $u_1 = s$, $u_2 = e$, $u_3 = i$, $u_4 = r$, and $u = (u_1,\ldots,u_4)$, we have $\Phi(u) = u_1+u_2+u_3 + u_4$ satisfying (<ref>)–(<ref>) with $b = \pi = 1$. Now (<ref>)–(<ref>) and (<ref>) are obviously fulfilled. To check (<ref>), we observe the first nonlinearity

\begin{align*}
&\alpha n - (1 - A_0/n) \beta_i si - (1 - A_0/n) \beta_e se - \mu s\\
&= \alpha n - \beta_i si + A_0\beta_i\frac{si}{n} - \beta_ese + A_0\beta_e \frac{se}{n} - \mu s\\
&\le \alpha n  + (A_0\beta_i + A_0\beta_e)s
\end{align*}

where we used $si/n \le s$ and $se/n\le e$ and the nonnegativity of the rates. Thanks to the boundedness of the rates, we see that the first nonlinearity is bounded above by a linear combination of all components. It's immediate that the third and forth nonlinearities, as well as the sum of the first and second nonlinearities are bounded by a linear combination of all components, since all nonlinear terms are cancelled out when summing. Hence, by choosing $r = 1$ and the matrix

\begin{equation*}
A = \begin{pmatrix}
1 & 0 & 0 & 0\\
1 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1
\end{pmatrix}
\end{equation*}

we see that (<ref>) is satisfied. It is noted that all these assumptions also hold if we replace $-\phi_di$ in (<ref>) by $-\phi_d in$ as done in [1]. Therefore, (<ref>), as well as the modified one in [1], is indeed a special case of the general system (<ref>).

§.§ Main results

Let us start with the first main result about the global existence and boundedness of system (<ref>).

Assume (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>) with

\begin{equation}\label{growth-condition}
1\leq r<1+b+\frac{2}{N}.
\end{equation}

Then for any nonnegative, bounded initial datum $u_0\in (L^{\infty}(\Omega))^m$, there exists a global weak solution to (<ref>)
with $u_i\in L^{\infty}_{loc}(0,\infty; L^{\infty}(\Omega))$ for all $i=1,\cdots,m.$  In particular, if $K_1<0$ or $K_1=K_2=0$, then the solution is bounded uniformly in time, i.e.

\begin{equation}\label{Linftybound}
\sup_{t\geq0}\|u_i(t)\|_{\LO{\infty}}<+\infty, \quad \forall i=1,\cdots,m.
\end{equation}

* When $b=0$, we have $\Phi(u)\geq M$. By letting $\Phi(u)=1$, our results cover the result in <cit.>.

* It is noted that under the assumption $K_1<0$ or $K_1 = K_2 = 0$, we can get the $L^1(\Omega)$-norm of the solution bounded uniformly in time (see Lemma <ref>). This fact can be used to show the uniform-in-time bound (<ref>). As a result, if we can use other ideas or use other structure of the system to get the $L^1(\Omega)$-norm of the solution bound uniformly in time, we can remove the assumption $K_1<0$ or $K_1 = K_2 = 0$.

* It is noted that the key idea to prove global existence is to construct the $L^p$-energy function. Actually, the intermediate sum condition (<ref>) is essential to construct the $L^p$-energy function.

The proof of Theorem <ref> is based on an $L^p$-energy approach. The traditional idea has been to construct an energy function that is decreasing or at least bounded in time for the (<ref>) of the form

\begin{equation*}
\mathscr{L}[u] = \sumi \intO h_i(u_i) \,dx.
\end{equation*}

It is noted that if $h_i(z) \sim z^p$, this yields an $L^p$-estimate of the solution, and for $p$ large enough, using a bootstrapping procedure, one obtains an $L^\infty$-estimate which implies global existence. Unfortunately, this approach is very likely to fail under the general assumptions (<ref>)-(<ref>), except for some very special cases.
The duality method (see e.g. [13, 18, 19, 24, 4, 21]) is very efficient when dealing with systems with constant or smooth diffusion coefficients. Using this method, one gets an $L^{2+\eps}(\Omega\times(0,T))$-estimate from the mass control condition (<ref>).
Combining this initial estimate, another duality argument (see [21]) and bootstrap argument, using the intermediate sum condition (<ref>) and the growth assumption (<ref>), one can obtain an  $L^\infty(\Omega\times(0, T))$-estimate that ensures global existence (more details can be found in [21]).
It is noted that this method does not extend to the case of merely bounded measurable or quasi-linear diffusion coefficients, unless some additional regularity assumptions are given (see [6, 2]). Recently, an $L^p$-energy approach was given in [12, 22],
whose preliminary ideas have been used previously in, e.g., [23, 14]. It is remarked that a similar method has also been successfully applied to cross-diffusion systems, see [16]. The core idea of this $L^p$-energy method is to construct a generalized $L^p$-energy function of the following form:

\begin{equation*}
\mathscr{L}_p[u](t)=\int_{\Omega} \mathscr{H}_p[u](t) \,d x,
\end{equation*}

where $p\in \mathbb{N}$ with $p\geq 2$ and

\begin{equation*}
\mathscr{H}_p[u](t)=\sum_{\beta \in \mathbb{Z}_{+}^m,|\beta|=p}\left(\begin{array}{l}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u(t)^\beta
\end{equation*}

\begin{equation*}
u(t)^\beta = \prod_{i=1}^{m}u_i^{\beta_i}, \theta^{\beta^2}= \prod_{i=1}^{m}\theta_i^{\beta^2_i} \text{ and }
\left(\begin{array}{l}
p \\
\beta
\end{array}\right)=\frac{p !}{\beta_{1} ! \cdots \beta_{m} !},
\end{equation*}

where $\theta=\left(\theta_1, \ldots, \theta_m\right)$ and $\theta_1, \ldots, \theta_m$ are positive real numbers which will be determined later. The function $\mathscr L_p[u]$ contains all (mixed) multi-variable polynomials of order $p$ with carefully chosen coefficients. Thanks to the non-negativity of the solution, $(\mathscr{L}_p[u])^{1/p}$ is an equivalent $L^p$-estimate of $u$. Certainly, the difficulty is to choose these coefficients so that they are compatible with both diffusion and reactions. In our case, the assumptions on the quasilinear diffusion (<ref>) and the intermediate sum condition (<ref>) allow us construct such an energy function.

Comparing to [12, 21], this work extends this $L^p$-energy method to the case of degenerate quasi-linear diffusion coefficients.  Moreover, it is also shown that this method is sufficiently robust to model variants and different boundary conditions.

Conditionally, if we can, by using a specific structure, get a better a prior estimate of solutions, the exponent $r$ in the intermediate sum condition (<ref>) can be enlarged. This is contained in the following theorem.

Assume (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>).
Suppose that there exists a constant $a\geq 1$ such that

\begin{equation}\label{condi-u-1}
\|u_i\|_{L^{\infty}(0,T;\LO{a})}\leq \mathscr{F}(T), \quad {\forall i=1,\cdots,m}
\end{equation}

\begin{equation}\label{condi-growth1}
1\leq r<1+b+\frac{2a}{N}.
\end{equation}

Then for any nonnegative, bounded initial datum $u_0\in (L^{\infty}(\Omega))^m$, there exists a global weak solution to (<ref>)
with $u_i\in L^{\infty}_{loc}(0,\infty; L^{\infty}(\Omega))$ for all $i=1,\cdots,m.$  In particular, if $\sup_{T\geq0}\mathscr{F}(T)<\infty$, then the solution is bounded uniformly in time, i.e.

\begin{equation}
\sup_{t\geq0}\|u_i(t)\|_{\LO{\infty}}<+\infty, \quad \forall i=1,\cdots,m.
\end{equation}

When $b=0$,  our results cover the result in <cit.>.

Assume (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>).
Suppose that there exists a constant $q> 1$ such that

\begin{equation}\label{condi-u-2}
\|u_i\|_{L^{q}(0,T;\LO{q})}\leq \mathscr{F}(T), \quad {\forall i=1,\cdots,m}
\end{equation}

\begin{equation}\label{condi-growth2}
1\leq r< 1 + \frac{N}{N+2}b + \frac{2q}{N+2}.
\end{equation}

Then for any nonnegative, bounded initial datum $u_0\in (L^{\infty}(\Omega))^m$, there exists a global weak solution to (<ref>)
with $u_i\in L^{\infty}_{loc}(0,\infty; L^{\infty}(\Omega))$ for all $i=1,\cdots,m.$  In particular, if $\sup_{T\geq0}\mathscr{F}(T)<\infty$, then the solution is bounded uniformly in time, i.e.

\begin{equation}
\sup_{t\geq0}\|u_i(t)\|_{\LO{\infty}}<+\infty, \quad \forall i=1,\cdots,m.
\end{equation}

When $b=0$,  our results cover the result in <cit.>.

It is worth noting that other conditions can also lead to a prior estimates. For example, the entropy condition has been considered in many papers, especially when it involves chemical reactions, see e.g. [5, 3, 28, 29]. This condition means that there exist the scalars $\mu_i\in\R$ such that

\begin{equation*}
\sumi f_i(x,t,u)(\log u_i+\mu_i)\leq 0, \quad \forall (x,t,u)\in \Omega\times\R_+\times \R^m_+.
\end{equation*}

This condition guarantees an $\LO1$-estimate for $u_i\log u_i$.
Actually, this condition guarantees an $\LO1$-estimate for $H(u(\cdot,t))$, where

\begin{equation*}
H(u)=\sumi u_i(\log u_i-1+\mu_i), \quad \forall u_i\geq 0.
\end{equation*}

Moreover, we can assume there exists a set $M=\prod_{k=1}^m(\alpha_i,\beta_i)$, where $\alpha_i,\beta_i$ are extended real numbers such that $\alpha_i<\beta_i$ for each $i=1,\ldots, m$, and solutions to (<ref>) remain in $M$ if initial data lies in $M$, a function $H:M\to \mathbb{R}_+$ that is $C^2$ and has the form

\[
H(u)=\sumi h_i(u_i),
\]

where $h_i:(\alpha_i,\beta_i)\to\mathbb{R}_+$ satisfies

\begin{equation}\label{H1}\tag{H1}
\begin{gathered}
h_i''(z) \ge 0, \quad \forall z\in (\alpha_i,\beta_i);\\
h_i(z) \text{ is bounded implies } z \text{ is bounded};\\
\nabla H(u)\cdot F(x,t,u)\le K_5\sumi h_i(u_i)+K_6, \quad \forall (x,t,u) \in\Omega\times\mathbb{R}_+\times\mathbb{R}_+^m
\end{gathered}
\end{equation}

for some $K_5, K_6>0$.
Then, we would obtain an $L^1(\Omega)$-estimate for $H(u(\cdot,t))$ (this has been introduced in [18]). In addition, intermediate sum conditions could also be written in the form

\begin{equation}\label{H2}\tag{H2}
A\begin{pmatrix}h_1'(u_1)F_1(x,t,u)\\\vdots\\h_m'(u_1)F_m(x,t,u)\end{pmatrix}\le K_7\vec{1}\left(\sumi h_i(u_i)+1\right)^r, \quad \forall(x,t,u)\in\Omega\times\mathbb{R}_+\times\mathbb{R}_+^m,
\end{equation}

that would lead to results in the same manner as we obtain from (<ref>) and (<ref>) above.

Assume (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>). Assume moreover that

\begin{equation}\label{growth-condition-1}
1\leq r<1+b+\frac{2}{N}.
\end{equation}

Then for any nonnegative, bounded initial datum $u_0\in (L^{\infty}(\Omega))^m$, there exists a global weak solution to (<ref>)
with $u_i\in L^{\infty}_{loc}(0,\infty; L^{\infty}(\Omega))$ for all $i=1,\cdots,m.$  In particular, if $K_5<0$ or $K_5=K_6=0$, then the solution is bounded uniformly in time, i.e.

\begin{equation}
\sup_{t\geq0}\|u_i(t)\|_{\LO{\infty}}<+\infty, \quad \forall i=1,\cdots,m.
\end{equation}

When $b=0$,  our results cover the result in <cit.>.

Our method is sufficiently robust to extend to other boundary conditions (for example, Robin-type boundary conditions). The precise results are given in the following theorem.

Consider the system

\begin{equation}\label{Sys_other_bc}
\begin{cases}
\pa_t u_i - \nabla\cdot(D_i(x,t)\Phi(u)\nabla u_i) = f_i(x,t,u), &x\in\Omega, t>0,\\
D_i(x,t)\Phi(u)\nabla u_i(x,t)\cdot \eta + \alpha_i u_i(x,t) = 0, &x\in\pa\Omega, t>0,\\
u_{i}(x,0) = u_{i,0}(x), &x\in\Omega,
\end{cases}
\end{equation}

where $\eta$ is the unit  outward normal vector on $\pa\Omega$, and $\alpha_i\geq 0$ for all $i=1,\ldots, m$.

Assume (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>). Moreover, assume either (<ref>) or (<ref>) or (<ref>) with

\begin{equation*}
0\leq r< \begin{cases}
1 + b + \frac{2}{N}, &\text{ in case of } \eqref{A3},\\
1 + b  + \frac{2a}{N}, &\text{ in case of } \eqref{condi-u-1},\\
1 + \frac{N}{N+2}b + \frac{2q}{N+2}, &\text{ in case of } \eqref{condi-u-2}.
\end{cases}
\end{equation*}

Then for any non-negative, bounded initial data $u_0\in \LO{\infty}^m$, there exists a  global weak solution to (<ref>) (see Definition <ref>) with $u_i\in L^{\infty}_{loc}(0,\infty;\LO{\infty})$ for all $i=1,\cdots, m$.
In particular, if $K_1< 0$ or $K_1= K_2 = 0$ in case of (<ref>), or $\sup_{T\geq0}\mathscr{F}(T)<\infty$ in case of (<ref>) or (<ref>), then the solution is bounded uniformly in time, i.e.

\begin{equation*}
\sup_{t\geq0}\|u_i(t)\|_{\LO{\infty}}<+\infty, \quad \forall i=1,\cdots,m.
\end{equation*}

We believe Theorem <ref> can also be extended to the semilinear boundary conditions of the form

\begin{equation*}
D_i(x,t)\nabla u_i(x,t)\cdot \eta + \alpha_i u_i(x,t) = G_i(u), \quad x\in\pa\Omega,
\end{equation*}

where the nonlinearities $G_i$ also satisfy a quasi-positivity condition and an intermediate sum condition. The details are left for the interested reader. We refer to [22, 27] for a related work dealing with constant diffusion coefficients.

The paper is organised as follows.
In Section <ref>, we prove the existence of Theorem <ref> by first considering an approximate system where we regularize the nonlinearities to obtain global approximate weak solutions. Then we derive uniform estimates, applying the key idea of the $L^p$-energy functions, and pass to the limit to obtain global existence of (<ref>). The uniform-in-time boundedness is proved in Section <ref>. The proofs of extended Theorems <ref>, <ref>, <ref>, and <ref> are given in Section <ref>. In Section <ref>, we show application of our results to a Susceptible-Exposed-Infected-Recovered (SEIRD) model and its variants. Finally, in Appendix <ref>, we give  the proof of Lemma <ref>.

Notation. In this paper we use the following notation, some of which will be recalled from time to time:

* For $T>0$ and $p\in [1,\infty]$, $Q_T:= \Omega\times(0,T)$ and

\begin{equation*}
L^p(Q_T):= L^p(0,T;L^p(\Omega))
\end{equation*}

equipped with the usual norm

\begin{equation*}
\|f\|_{\LQ{p}}:= \bra{\int_0^T\intO |f|^pdxdt}^{1/p}
\end{equation*}

for $1\leq p < \infty$ and

\begin{equation*}
\|f\|_{\LQ{\infty}}:= \underset{(x,t)\in Q_T}{\text{\normalfont ess sup}}|f(x,t)|.
\end{equation*}

* For  $p\in [1,\infty]$, $\tau \geq 0$ and $\delta>0$, we denote by

\begin{equation*}
Q_{\tau,\tau+\delta}:= \Omega\times(\tau,\tau+\delta)
\end{equation*}

\begin{equation*}
L^p(Q_{\tau,\tau+\delta}):= L^p(\tau,\tau+\delta; L^p(\Omega)).
\end{equation*}

§ GLOBAL EXISTENCE OF BOUNDED WEAK SOLUTIONS

A vector of nonnegative state variables $u=(u_1,\cdots,u_m)$ is called a weak solution to (<ref>) on $(0,T)$ if

$$\pa_tu_i\in L^2(0,T;(H^1(\Omega))'), \quad \Phi(u)\nabla u_i\in L^2(0,T; L^2(\Omega)) $$

$$\nabla (u_i^{\frac{b+2}{2}})\in L^2(0,T; L^2(\Omega)) $$

with $u_i(\cdot, 0)=u_{i,0}(\cdot)$ for all $i=1,\cdots,m$, and  for any test function $ \red{\psi} \in L^2(0,T; H^1(\Omega))$
we have

\begin{equation*}
\begin{split}
\int^T_0 \intO \pa_t u_i  \psi\,dxdt+\int^T_0 \intO D_i(x,t)\Phi(u)\nabla u_i \cdot \nabla \psi \,dxdt
=\int^T_0 \intO f_i(u)\psi \,dxdt.
\end{split}
\end{equation*}

§.§ Approximate system

For $i=1,\cdots,m$ and $\eps>0$, consider the following approximating system for $u^{\eps}=(u_1^{\eps}, \cdots, u_m^{\eps})$,

\begin{equation}\label{appro-system}
\begin{cases}
\pa_t u_i^{\eps}-\nabla\cdot(D_i(x,t)\Phi(u^{\eps})\nabla u_i^{\eps})=f^{\eps}_i(u^{\eps}),  &x\in \Omega, t>0,\\
(D_i(x,t) \Phi(u^{\eps}) \nabla u^{\eps}_i)\cdot \eta=0,   &x\in \pa \Omega, t>0,\\
u^{\eps}_i(x,0)=u_{i,0}(x)+\eps,   &x\in \Omega,
\end{cases}
\end{equation}

\begin{equation*}
f^{\eps}_i(u^{\eps}):=\frac{f_i(u^{\eps})}{1+\eps \sum^m_{j=1}|f_j(u^{\eps})|}
\end{equation*}

and $u_{i,0}^{\eps}\in \LO{\infty}$ is non-negative and $u_{i,0}^{\eps}\stackrel{\eps\rightarrow0}{\longrightarrow}u_{i,0}$ in $\LO{\infty}$. With this approximation, it is easy to check that the approximated non-linearities $f_i^{\eps}$ still satisfy the assumptions (<ref>)-(<ref>).

For any fixed $\eps>0$, there exists a global bounded, nonnegative weak solution to (<ref>) on any finite time interval $(0,T)$, $T>0$.

Since for fixed $\eps>0$, the nonlinearities $f_i^{\eps}(u^{\eps})$ are Lipschitz continuous and bounded, i.e.,

\begin{equation*}
f^{\eps}_i(u^{\eps}):=\frac{f_i(u^{\eps})}{1+\eps \sum^m_{j=1}|f_j(u^{\eps})|}\leq \frac{1}{\eps}, \quad \forall (x,t)\in Q_T.
\end{equation*}

The global existence of a weak solution of (<ref>) is standard (Using Galerkin approach).

Next we prove the nonnegativity of $u^{\eps}$. Denote $u_{i,+}^{\eps}:=\max\{u_i^{\eps},0\}$ and $u_{i,-}^{\eps}:=\min\{u_i^{\eps},0\}$. We consider the auxiliary system of (<ref>)

\begin{equation*}
\pa_t u_i^{\eps}-\nabla\cdot(D_i(x,t)\Phi(u^{\eps})\nabla u_i^{\eps})=f^{\eps}_i(u_+^{\eps}),
\end{equation*}

where $u_+^{\eps}=(u_{i,+}^{\eps})_{i=1,\cdots,m}$.

By multiplying this auxiliary system by $u_{i,-}^{\eps}$ and using the quasi-positivity assumption (<ref>) (recall that this property also holds for $f_i^{\eps}$), we obtain

\begin{equation*}
\frac{1}{2}\intO |u_{i,-}^{\eps}|^2 \,dx+\lambda\intO \Phi(u^{\eps}) |\nabla u_{i,-}^{\eps}|^2\,dx \leq 0,
\end{equation*}

where we use $u_{i,0,-}^{\eps}=0$. Thus, we get $u_{i,-}^{\eps}=0$ a.e. in $Q_T$, this shows the desired nonnegativity.

§.§ Uniform-in-$\eps$ estimate

In this subsection, we prove crucial uniform-in-$\eps$ estimates for the solution to (<ref>). Moreover, we want to emphasize that all the constants in this subsection are independent of $\eps$.

The following bound in $L^{\infty}(0,T; L^1(\Omega))$ is immediate.

Assume (<ref>), (<ref>) and (<ref>). Then  for any $T>0$, there exists a constant $M_T$ depending on $T, \Omega, \|u_{i,0}\|_{L^1(\Omega)}$ and $c_1,\cdots,c_m$, $K_1, K_2$ in (<ref>) such that

\begin{equation*}
\sup_{t\in(0,T)}\|u_i^{\eps}(t)\|_{\LO1}\leq M_T, \quad \forall i=1,\cdots,m.
\end{equation*}

By summing the equation in (<ref>), integrating on $\Omega$ and using (<ref>) we have

\begin{equation*}
\frac{d}{dt}\sumi\intO c_iu_i^{\eps} \,dx \leq \intO \left(K_1 \sumi u_i^{\eps}(x,t) + K_2\right) \,dx.
\end{equation*}

The classical Gronwall inequality gives the desired estimate.

The following  $L^p$-energy function has been developed in [22, 12].
We write $\mathbb{Z}_{+}^m$ as the set of all $m$ tuples of nonnegative integers. Addition and scalar multiplication by nonnegative integers of elements in $\mathbb{Z}_{+}^m$ is understood in the usual manner. If $\beta=\left(\beta_1, \ldots, \beta_m\right) \in \mathbb{Z}_{+}^m$ and $p \in \mathbb{Z}_+$, then we define $\beta^p=\left(\left(\beta_1\right)^p, \ldots,\left(\beta_m\right)^p\right)$. Also, if $\alpha=\left(\alpha_1, \ldots, \alpha_m\right) \in$ $\mathbb{Z}_{+}^m$, then we define $|\alpha|=\sum_{i=1}^m \alpha_i$. Finally, if $z=\left(z_1, \ldots, z_m\right) \in \mathbb{R}_{+}^m$ and $\alpha=$ $\left(\alpha_1, \ldots, \alpha_m\right) \in \mathbb{Z}_{+}^m$, then we define $z^\alpha=z_1^{\alpha_1} \cdots z_m^{\alpha_m}$, where we interpret $0^0$ to be 1. For any $2\leq p \in \mathbb N$, we build our $L^p$-energy function of the form

\begin{equation}\label{def-Lp-energy-fun}
\mathscr{L}_p[u](t)=\int_{\Omega} \mathscr{H}_p[u](t) \,d x,
\end{equation}

\begin{equation}\label{def-H-p}
\mathscr{H}_p[u](t)=\sum_{\beta \in \mathbb{Z}_{+}^m,|\beta|=p}\left(\begin{array}{l}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u(t)^\beta
\end{equation}

\begin{equation*}
\left(\begin{array}{l}
p \\
\beta
\end{array}\right)=\frac{p !}{\beta_{1} ! \cdots \beta_{m} !}
\end{equation*}

and $\theta=\left(\theta_1, \ldots, \theta_m\right)$, where $\theta_1, \ldots, \theta_m$ are positive real numbers which will be determined later. For $p=0,1,2$, one can write these functions explicitly as

\begin{equation*}
\mathscr{H}_0[u](t)=1 \text { and } \mathscr{H}_1[u](t)=\sum_{j=1}^m \theta_j u_j(t)
\end{equation*}

\begin{equation*}
\mathscr{H}_2[u](t)=\sum_{i=1}^m \theta_i^4 u_i(t)^2+2 \sum_{i=1}^{m-1} \sum_{j=i+1}^m \theta_i \theta_j u_i(t) u_j(t) .
\end{equation*}

Thanks to the nonnegativity of the solution, we have

\begin{equation*}
\mathscr{L}_p[u](t) \sim \sum_{i=1}^m\|u_i(t)\|_{L^p(\Omega)}^p ~\text{ for } p\geq 1.
\end{equation*}

This will allow us to use $\mathscr{L}_p[u](t)$ to obtain a priori estimates on $u$ for each $2 \leq p \in \mathbb{N}$. We will need two technical lemmas.

Suppose $m \in \mathbb{N}, \theta=\left(\theta_1, \ldots, \theta_m\right)$, where $\theta_1, \ldots, \theta_m$ are positive real numbers, $\beta \in \mathbb{Z}_{+}^m$, and $\mathscr{H}_p[u]$ is defined in (<ref>). Then

\frac{\partial}{\partial t} \mathscr{H}_0[u](t)=0, \quad \frac{\partial}{\partial t} \mathscr{H}_1[u](t)=\sum_{j=1}^m \theta_j \frac{\partial}{\partial t} u_j(t)

and for $p \in \mathbb{N}$ such that $p \geq 2$,

\frac{\partial}{\partial t} \mathscr{H}_p[u](t)=\sum_{|\beta|=p-1}\left(\begin{array}{l}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u(t)^\beta \sum_{j=1}^m \theta_j^{2 \beta_j+1} \frac{\partial}{\partial t} u_j(t).

Suppose $m \in \mathbb{N}, \theta=\left(\theta_1, \cdots, \theta_m\right)$, where $\theta_1, \cdots, \theta_m$ are positive real numbers, and let $\mathscr{H}_p[u]$ be defined in (<ref>). If $p \in \mathbb{N}$ such that $p \geq 2$, then

\begin{equation*}
\begin{aligned}
& \sum_{|\beta|=p-1}\left(\begin{array}{l}
p \\
\beta
\end{array}\right) \theta^{\beta^2} \sum_{k=1}^m \theta_k^{2 \beta_k+1}\left(A_k \nabla u_k\right) \cdot \nabla u^\beta \\
& =\sum_{|\beta|=p-2}\left(\begin{array}{c}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u^{\beta} \sum_{k=1}^m \sum_{l=1}^m C_{k, l}(\beta)\left(A_k \nabla u_k\right) \cdot \nabla u_l,
\end{aligned}
\end{equation*}

\begin{equation*}
C_{k, l}(\beta)=\left\{\begin{array}{cc}
\theta_k^{2 \beta_k+1} \theta_l^{2 \beta_l+1}, & k \neq \red{l}, \\
\theta_k^{4 \beta_k+4}, & k= \red{l}.
\end{array}\right.
\end{equation*}

Next, we need the following functional inequality.

Suppose $\Omega\subset \mathbb{R}^N$ such that the Gagliardo-Nirenberg inequality is satisfied and basic trace theorems apply (for instance $\Omega$ has a Lipschitz boundary). Let $a\geq 1$, $p\geq 2a$, $w: \overline{\Omega}\rightarrow\mathbb{R}_+$ such that $w^{\frac{b+p}{2}}\in W^{1,2}(\Omega)$ and there exists $K\geq0$ such that $\|w\|_{L^a(\Omega)}\leq K$. If $1\leq r<1+b+ \frac{2a}{N}$ and $b\geq 0$, then there exists $C_{\varepsilon}\geq 0$ (depending on $p,\varepsilon, r, a, b, K, \Omega,$ but independent of $w$) such that

\begin{equation*}
\int_{\Omega}w^{p-1+r} \,dx +\int_{\Omega}w^{b+p}\,dx \leq \varepsilon \left(\int_{\Omega} w^{p-2+b}|\nabla w|^{2}\,dx+\int_{\Omega} w^{b+p}\,dx\right)+C_{\varepsilon}.
\end{equation*}

By using Sobolev's embedding we have

\begin{equation}\label{eq-L1-1}
\begin{split}
\int_{\Omega} &w^{p-2+b}|\nabla w|^{2}\,dx+\int_{\Omega} w^{b+p}\,dx\\
&=\left(\frac{2}{b+p}\right)^{2}\int_{\Omega} |\nabla w^{\frac{b+p}{2}}|^{2}\,dx+\int_{\Omega} (w^{\frac{b+p}{2}})^{2}\,dx\\
&\geq \min\left\{\left(\frac{2}{b+p}\right)^{2}, 1\right\}\|w^{\frac{b+p}{2}}\|^{2}_{H^{1}(\Omega)}.
\end{split}
\end{equation}

Thus we have

\begin{equation}\label{eq-l1}
\int_{\Omega}w^{b+p}\,dx\leq \eps\left(\int_{\Omega} w^{p-2+b}|\nabla w|^{2}\,dx+ \int_{\Omega} w^{b+p}\,dx\right)+C(p,\varepsilon,b,K).
\end{equation}

\begin{equation*}
\int_{\Omega}w^{p-1+r}\,dx=\|w^{\frac{p-1+r}{2}}\|^2_{\LO2},
\end{equation*}

if $r<b+1$. Then we have

\begin{equation*}
\int_{\Omega}w^{p-1+r} \,dx \leq \varepsilon \int_{\Omega} w^{b+p}\,dx +C_{\varepsilon}.
\end{equation*}

Next we consider $r\geq b+1$.
Thanks to the Gagliardo-Nirenberg inequality, we have

\begin{equation}\label{eq-L1-2}
\begin{split}
\int_{\Omega}w^{p-1+r}\,dx&=\|w^{\frac{b+p}{2}}\|_{L^{\frac{2(r+p-1)}{b+p}} (\Omega)}^{\frac{2(r+p-1)}{b+p}}\\
&\leq C \|\nabla(w^{\frac{b+p}{2}})\|_{L^{2}(\Omega)}^{\frac{2(r+p-1)}{b+p}\alpha} \|w^{\frac{b+p}{2}}\|_{L^{1} (\Omega)}^{\frac{2(r+p-1)}{b+p}(1-\alpha)},
\end{split}
\end{equation}

\begin{equation*}
\begin{split}
\frac{b+p}{2(r+p-1)}=\alpha(\frac{1}{2}-\frac{1}{N})+1-\alpha.
\end{split}
\end{equation*}

It follows that

\begin{equation*}
\begin{split}
\alpha=\frac{N[2(r+p-1)-(b+p)]}{(r+p-1)(N+2)}
\end{split}
\end{equation*}

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

Since $1\leq r<1+b+\frac{2a}{N}$, we find

\begin{equation*}
\alpha\frac{2(r+p-1)}{b+p}<2.
\end{equation*}

We can  use Young's inequality to estimate

\begin{equation}\label{eq-L1-3}
\begin{split}
\|w^{\frac{b+p}{2}}\|_{L^{\frac{2(r+p-1)}{b+p}}(\Omega)}^{\frac{2(r+p-1)}{b+p}}
&\leq \varepsilon \|\nabla(w^{\frac{b+p}{2}})\|_{L^{2}(\Omega)}^{2}+C_{\varepsilon}\|w^{\frac{b+p}{2}}\|_{L^{1} (\Omega)}^{\frac{(1-\alpha)2(r+p-1)}{b+p-(r+p-1)\alpha}},
\end{split}
\end{equation}

\begin{equation}\label{eq-L1-4}
\begin{split}
C_{\varepsilon}\|w^{\frac{b+p}{2}}\|_{L^{1} (\Omega)}^{\frac{(1-\alpha)2(r+p-1)}{b+p-(r+p-1)\alpha}}=C_{\varepsilon}\|w \|_{L^{\frac{b+p}{2}} (\Omega)}^{\frac{(1-\alpha)(b+p)(r+p-1)}{b+p-(r+p-1)\alpha}}.
\end{split}
\end{equation}

If $b+p=2a$, then this term is bounded by a constant depending on $K$, since $\|w\|_{L^a(\Omega)}\leq K$.
If $b+p>2a$, we use an interpolation inequality to have

\begin{equation}\label{eq-L1-5}
\begin{split}
C_{\varepsilon}\|w\|_{L^{\frac{b+p}{2}} (\Omega)}\leq C_{\varepsilon}\|w\|^{\theta}_{L^{p+r-1} (\Omega)}\|w\|^{1-\theta}_{L^{a} (\Omega)}\leq C_{\varepsilon, K} \|w\|^{\theta}_{L^{p+r-1} (\Omega)},
\end{split}
\end{equation}

where $\theta\in(0,1)$ satisfies

Note that

\begin{equation}\label{eq-L1-6}
\theta\frac{(1-\alpha)(b+p)(r+p-1)}{b+p-(r+p-1)\alpha}<p+r-1
\end{equation}

due to $r<b+1+ \frac{2a}{N}$ and $b\leq r-1<2r$.

From  (<ref>)-(<ref>) and Young's inequality, it follows that

\begin{equation}\label{eq-7}
\begin{split}
C_{\varepsilon}\|w \|^{\frac{(1-\alpha)(b+p)(r+p-1)}{b+p-(r+p-1)\alpha}}_{L^{\frac{b+p}{2}} (\Omega)}&\leq C_{\varepsilon, K}
\|w \|^{\frac{(1-\alpha)(b+p)(r+p-1)}{b+p-(r+p-1)\alpha}\theta}_{L^{p+r-1} (\Omega)}\\
&\leq \frac{1}{2} \| w \|^{p+r-1}_{L^{p+r-1}(\Omega)}+C(p,\varepsilon,b,K).
\end{split}
\end{equation}

Inserting this into (<ref>), we get

\begin{equation*}
\begin{split}
\|w\|^{p+r-1}_{L^{p+r-1}(\Omega)}\leq 2\varepsilon \|w^{\frac{b+p}{2}}\|^{2}_{W^{1,2}(\Omega)} +C(p,\varepsilon,b,K).
\end{split}
\end{equation*}

Combine this with (<ref>) and (<ref>) leads to the desired estimate

\begin{equation*}
\int_{\Omega}w^{p-1+r} +\int_{\Omega}w^{b+p}\,dx \leq \varepsilon \left(\int_{\Omega} w^{p-2+b}|\nabla w|^{2}\,dx+\int_{\Omega} w^{b+p}\,dx\right)+C(p,\varepsilon,b,K).
\end{equation*}

The following lemma shows results for the intermediate sum condition (<ref>), which is crucial for constructing the $L^p$ energy function.

Assume (<ref>). Then there exist componentwise increasing functions $g_i: \mathbb{R}^{m-i} \rightarrow \mathbb{R}_{+}$ for $i=1, \ldots, m-1$ such that if  $\theta=\left(\theta_1, \ldots, \theta_m\right) \in(0, \infty)^m$ satisfies $\theta_m>0$ and $\theta_i \geq g_i\left(\theta_{i+1}, \ldots, \theta_m\right)$ for all $i=1, \ldots, m-1$, then

\sum_{i=1}^m \theta_i f_i^{\eps}\left(x, t, u^{\varepsilon}\right) \leq K_\theta\left(1+\sum_{i=1}^m\left(u_i^{\varepsilon}\right)^r\right) \quad \forall\left(x, t, u^{\varepsilon}\right) \in \Omega \times \mathbb{R}_{+} \times \mathbb{R}_{+}^m

for some constant $K_\theta$ depending on $\theta, g_i$, and $K_3$ in (<ref>).

We now use the $L^p$-energy functions  to obtain the $L^p$-estimates of $u^{\eps}$.

Assume (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>),  (<ref>) and (<ref>) with $1\leq r<1+b+\frac{2}{N}$. Then for any $1 \leq p<\infty$ and any $T>0$, there exists a constant $C_{T, p}$ depending on $T, p$ and other parameters such that

\sup _{t \in(0, T)}\|u_i^{\eps}(t)\|_{L^p(\Omega)} \leq C_{T, p} \quad \forall i=1, \ldots, m .

Let $u^{\eps}$ solve (<ref>), and $\mathscr{L}_p(t):=\mathscr{L}_p\left[u^{\eps}\right](t)$ be defined in (<ref>). Then

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t)= & \int_{\Omega} \sum_{|\beta|=p-1}\left(\begin{array}{c}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u^{\eps}(x, t)^\beta \sum_{k=1}^m \theta_k^{2 \beta_k+1} \frac{\partial}{\partial t} u_k^{\eps}(x, t)\, d x \\
= & \int_{\Omega} \sum_{|\beta|=p-1}\left(\begin{array}{c}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u^{\eps}(x, t)^\beta \sum_{k=1}^m \theta_k^{2 \beta_k+1} \\
& \times\left[\nabla \cdot\left(D_k(x, t)\,\Phi(u^{\eps})\, \nabla u_k^{\eps}(x, t)\right)+f^{\eps}_k\left(u^{\eps}(x, t)\right)\right]\, d x\\
= & \int_{\Omega} \sum_{|\beta|=p-1}
\left(\begin{array}{c}
p \\
\beta
\end{array}\right)
\theta^{\beta^2} u^{\eps}(x, t)^\beta \sum_{k=1}^m \theta_k^{2 \beta_k+1}\nabla \cdot\left(D_k(x, t) \,\Phi(u^{\eps})\, \nabla u_k^{\eps}(x, t)\right)\, d x\\
&+ \int_{\Omega} \sum_{|\beta|=p-1}\left(\begin{array}{c}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u^{\eps}(x, t)^\beta \sum_{k=1}^m \theta_k^{2 \beta_k+1} f^{\eps}_k\left(u^{\eps}(x, t)\right)\, d x\\
\end{split}
\end{equation*}

For (I),  we apply Lemma <ref> and integration by parts, we have

\begin{equation*}
\begin{split}
(I):=& \int_{\Omega} \sum_{|\beta|=p-1}
\left(\begin{array}{c}
p \\
\beta
\end{array}\right)
\theta^{\beta^2} u^{\eps}(x, t)^\beta \sum_{k=1}^m \theta_k^{2 \beta_k+1}\nabla \cdot\left(D_k(x, t) \,\Phi(u^{\eps})\, \nabla u_k^{\eps}(x, t)\right)\, d x\\
=&-\int_{\Omega} \sum_{|\beta|=p-2}
\left(\begin{array}{l}
p \\
\beta
\end{array}\right)
\theta^{\beta^2} u^{\varepsilon}(x, t)^\beta \,\Phi(u^{\eps})\, \sum_{k=1}^m \sum_{l=1}^m C_{k, l}(\beta)\left(D_k(x,t) \nabla u_k^{\eps}\right) \cdot \nabla u_l^{\eps} \,dx
\end{split}
\end{equation*}

\begin{equation*}
C_{k, l}(\beta)=\left\{\begin{array}{cc}
\theta_k^{2 \beta_k+1} \theta_l^{2 \beta_l+1}, & k \neq l, \\
\theta_k^{4 \beta_{k}+4}, & k=l.
\end{array}\right.
\end{equation*}

For a given $\beta$ with $|\beta|=p-2$, create an $m N \times m N$ matrix $B(\beta)$ made up of $m^2$ blocks $B_{k, l}(\beta)$, each of size $N \times N$, where

B_{k, l}(\beta)=\frac{1}{2} C_{k, l}(\beta)\left(D_k+D_l\right).

Note that for each $k=1, \ldots, m$,

B_{k, k}(\beta)=\theta_k^{4 \beta_k+4} D_k.

I=-\int_{\Omega} \sum_{|\beta|=p-2}\left(\begin{array}{c}
p \\
\beta
\end{array}\right) \theta^{\beta^2} u^{\eps}(x, t)^{\beta} \,\Phi(u^{\eps})\, \nabla u^{\eps}(x, t)^T B(\beta) \nabla u^{\eps}(x, t)\, dx,

where $\nabla u^{\eps}(x, t)$ is a column vector of size $m N \times 1$, and for $j=1, \ldots, m$, entries $N(j-1)+1$ to $N j$ of $\nabla u^{\varepsilon}(x, t)$ are $\nabla u_j^{\eps}(x, t)$. We claim that if all of the entries in $\theta$ are sufficiently large, then $B(\beta)$ is positive definite. In fact, it is a simple matter to show it is positive definite if and only if the $m N \times m N$ matrix $\widetilde{B}(\beta)$ made up of $N \times N$ blocks

\widetilde{B}_{k, l}(\beta)=\left\{\begin{array}{cc}
\theta_k^2 D_k, & k=l \\
\frac{1}{2}\left(D_k+D_l\right), & k \neq l
\end{array}\right.

is positive definite. However, if we recall the uniform positive definiteness of the matrices $D_k$, we can show that if $\theta_i$ is sufficiently large for each $i$, then we have what we need. Consequently, returning to the above, we can show there exists $\alpha_p>0$ so that

\begin{equation}\label{eq-E}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\alpha_p \sum_{k=1}^m \int_{\Omega}u_k^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_k(x, t)\right|^2 \,dx \\
&\leq \int_{\Omega} \sum_{|\beta|=p-1}
\left(\begin{array}{l}
p \\
\beta
\end{array}\right)
\theta^{\beta^2} u^{\eps}(x, t)^\beta \sum_{k=1}^m \theta_k^{2 \beta_k+1}f_k^{\eps}\left(u^{\eps}(x, t)\right)\, dx,
\end{split}
\end{equation}

where we used (<ref>).

From (<ref>) and Lemma <ref>,  we choose the components of $\theta = (\theta_i)$ inductively so that $\theta_i$ are sufficiently large that the previous positive definiteness condition of $\widetilde{B}_{k, l}(\beta)$ is satisfied, and

\begin{equation}\label{star}
\begin{split}
\theta_i\geq g_i(\theta_{i+1}^{2p-1},\cdots,\theta_{m}^{2p-1}) \quad\text{for}\quad i=1,\cdots,m-1,
\end{split}
\end{equation}

where $g_i$ are functions constructed in Lemma <ref>. Note that $\theta_i \leq \theta_i^{2\beta_i + 1} \leq \theta_i^{2p - 1}$. Since $g_i$ is componentwise increasing,  the relation (<ref>) implies

\begin{equation*}
\theta_i^{2\beta_i + 1} \geq g_i\bra{\theta_{i+1}^{2\beta_{i+1}+1}, \ldots, \theta_m^{2\beta_m + 1}}, \quad \forall i=1,\ldots, m-1.
\end{equation*}

Now we can apply Lemma <ref>, to obtain some $K_{\widetilde{\theta}}$  so that for all $\beta\in\mathbb{Z}_+$ with $|\beta|=p-1$, we have

\begin{equation*}
\begin{split}
\sumi \theta_{i}^{2\beta_i+1}f^{\eps}_i(x,t,u^{\eps})\leq K_{\widetilde{\theta}}\bra{1+\sumi (u^{\eps}_i)^r}, \quad \forall (x,t,u)\in \Omega\times\mathbb{R}_+\times \mathbb{R}^m_+.
\end{split}
\end{equation*}

It follows that there exists $C_p>0$, such that (<ref>) implies

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\alpha_p \sum_{k=1}^m \int_{\Omega}u_k^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_k(x, t)\right|^2 \,dx \leq C\left(\intO\sumi u_i^{\eps}(x,t)^{p-1+r}\,dx+1\right).
\end{split}
\end{equation*}

Thus, we have

\begin{equation}\label{p}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\alpha_p \sum_{k=1}^m \int_{\Omega}\left(u_k^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_k(x, t)\right|^2 +u_i^{\eps}(x,t)^{b+p}\right)\,dx\\
& \leq C\left(\intO\sumi \left(u_i^{\eps}(x,t)^{p-1+r} + u_i^{\eps}(x,t)^{b+p}\right) \,dx+1\right).
\end{split}
\end{equation}

Combining Lemma <ref>, we get

\begin{equation}\label{eq-p}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\frac{\alpha_p}{2} \sum_{k=1}^m \int_{\Omega}\left(u_k^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_k(x, t)\right|^2 +u_i^{\eps}(x,t)^{b+p}\right)\,dx \leq C(T,p).
\end{split}
\end{equation}

This implies

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t) + \sigma\mathscr{L}_p(t)  \leq C(T,p).
\end{split}
\end{equation*}

Thus, we have

\begin{equation*}
\sup _{t \in(0, T)}\|u_i^{\eps}(t)\|_{L^p(\Omega)} \leq C_{T, p} \quad \forall i=1, \ldots, m .
\end{equation*}

This finishes the proof of Lemma <ref>.

Assume (<ref>), (<ref>),  (<ref>), (<ref>), (<ref>), (<ref>),  (<ref>), (<ref>) and (<ref>) with $1\leq r<1+b+\frac{2}{N}$. Then for any $T>0$, the solution of (<ref>) is bounded in $L^{\infty}$ in time, i.e.,

\begin{equation*}
\|u_i^{\eps}(t)\|_{L^{\infty}(Q_T)} \leq C_{T} \quad \forall i=1, \cdots, m
\end{equation*}

for some constant $C_T$ depending on $T$ and independent of $\eps>0$.

The proof is similar to the proof of Lemma 2.3 in Ref. [9]. For the convenience of reading, we give the specific proof process in the Appendix <ref>.

§.§ Passing to the limit - Global existence

Assume (<ref>). For any $k>0$, we have

\begin{equation*}
\int^T_0\intO \chi_{\{|u^{\eps}_i|\leq k\}} \Phi(u^{\eps}) |\nabla u^{\eps}_i|^2 \,dxdt \leq \frac{k}{\lambda} \left(\|u_0\|_{\LO 1}+\|f^{\eps}_i(u^{\eps})\|_{\LQ 1}\right).
\end{equation*}

Let $T_k(z)$ be defined

\begin{equation*}
-k  & if z\leq-k,\\
z   & if -k<z<k,\\
k   & if z\geq k.
\end{cases}
\end{equation*}

Define $S_k(z)=\int^z_0 T_k(\tau) \,d\tau$ and multiply (<ref>) by $T_k(u^{\eps}_i)$ we obtain

\begin{equation*}
\begin{split}
\intO S_k(u_i^{\eps})\,dx & + \int^T_0\intO D_i(x,t) \Phi(u^{\eps})  \nabla u_i^{\eps} \cdot \nabla(T_k(u^{\eps}_i))\, dxdt\\
&=\intO S_k(u_0^{\eps})\,dx + \int^T_0\intO f^{\eps}_i(u^{\eps})T_k(u^{\eps}_i)\,dxdt.
\end{split}
\end{equation*}

By applying the properties of $T_k$ and $S_k$, the right hand is bounded by

\begin{equation*}
\begin{split}
\intO S_k(u_0^{\eps})\,dx + \int^T_0\intO f^{\eps}_i(u^{\eps})T_k(u^{\eps}_i)\,dxdt\leq k \left(\|u_0\|_{\LO1}+\|f_i^{\eps}\|_{\LQ1} \right).
\end{split}
\end{equation*}

From (<ref>) and $ \nabla(T_k(u^{\eps}_i))=\chi_{\{|u^{\eps}_i|\leq k\}} \nabla u^{\eps}_i$, it follows that

\begin{equation*}
\begin{split}
\int^T_0\intO D_i(x,t) \Phi(u^{\eps}) \nabla u_i^{\eps} \cdot \nabla(T_k(u^{\eps}_i))\, dxdt &\geq \lambda \int^T_0\intO \chi_{\{|u^{\eps}_i|\leq k\}} \Phi(u^{\eps})  |\nabla u_i^{\eps}|^2 \, dxdt.
\end{split}
\end{equation*}

Finally, by using $S_k(u^{\eps}_i)\geq 0$ we obtain our desired estimate.

Assume (<ref>). For any $\beta>0$, there exists constants $C$ such that

\begin{equation}\label{eq-nabla-u}
\int^T_0\intO \frac{\Phi(u^{\eps}) |\nabla u^{\eps}_i|^2}{(1+|u^{\eps}_i|)^{1+\beta} } \,dxdt \leq C \left(\|u_0\|_{\LO 1} + \|f^{\eps}_i(u^{\eps})\|_{\LQ 1}\right).
\end{equation}

Let $M:=\|u_0\|_{\LO 1}+\|f^{\eps}_i(u^{\eps})\|_{\LQ 1}$, we can apply Lemma <ref> to obtain

\begin{equation*}
\begin{split}
\int^T_0\int_{\Omega} \frac{ \Phi(u^{\eps})  |\nabla u^{\eps}_i|^2}{(1+|u^{\eps}_i|)^{1+\beta} } \,dxdt &=\sum^{\infty}_{j=0} \int^T_0\intO \mathbf{1}_{\{2^j-1\leq |u^{\eps}_i|\leq 2^{j+1}-1\}}  \frac{ \Phi(u^{\eps}) |\nabla u^{\eps}_i|^2}{(1+|u^{\eps}_i|)^{1+\beta} } \,dxdt\\
&\leq \sum^{\infty}_{j=0}2^{-j(1+\beta)} \int^T_0\intO \mathbf{1}_{\{|u^{\eps}_i|\leq 2^{j+1}-1\}}  \Phi(u^{\eps})  |\nabla u^{\eps}_i|^2 \,dxdt\\
&\leq \frac{M}{\lambda} \sum^{\infty}_{j=0}2^{-j(1+\beta)}  2^{j+1}\\
&\leq \frac{2M}{\lambda} \sum^{\infty}_{j=0}(2^{-\beta})^{j}.
\end{split}
\end{equation*}

Thus (<ref>) holds for $C:=\frac{2}{\lambda} \sum^{\infty}_{j=0}(2^{-\beta})^{j}$, which is finite because of $\beta>0$.

From subsection <ref> we have the following bound

\begin{equation*}
\|u_i^{\eps}(t)\|_{L^{\infty}(Q_T)} \leq C_{T} \quad \forall i=1, \cdots, m.
\end{equation*}

Due to the polynomial growth $\eqref{A5}$, we can get

\begin{equation}\label{eq-f}
\|f_i(u_i^{\eps}(t))\|_{L^{\infty}(Q_T)} \leq C_{T} \quad \forall i=1, \cdots, m.
\end{equation}

By multiply (<ref>) by $u_i^{\eps}$ then integrating on $Q_T$ gives

\begin{equation*}
\begin{split}
\frac{1}{2}\|u^{\eps}_i(T)\|^2_{\LO 2}&+\int^T_0\intO D_i(x,t) \Phi(u^{\eps})  \nabla u_i^{\eps} \cdot \nabla u_i^{\eps} \,dxdt\\
&=\frac{1}{2}\|u^{\eps}_i(0)\|^2_{\LO 2}+\int^T_0\intO f_i^{\eps}(u^{\eps}) u_i^{\eps} \,dxdt\\
&\leq \frac{1}{2}\|u^{\eps}_i(0)\|^2_{\LO 2}+\|f_i^{\eps}(u^{\eps}) \|_{\LQ {\infty}}\|u_i^{\eps}\|_{\LQ 1}\\
&\leq C\left(T, |\Omega|, \|u_{i,0}\|_{\LO{\infty}}\right).
\end{split}
\end{equation*}

By (<ref>), (<ref>)  and Hölder inequality give

\begin{equation*}
\begin{split}
\int^T_0&\intO D_i(x,t) \Phi(u^{\eps})  \nabla u_i^{\eps} \cdot \nabla u_i^{\eps} \,dxdt \\
&\geq \lambda M_i \int^T_0\intO  (u_i^{\eps})^b |\nabla u_i^{\eps}|^2\,dxdt \\
&=\lambda M_i (\frac{2}{b+2})^2 \int^T_0\intO  |\nabla ((u_i^{\eps})^{\frac{b+2}{2}})|^2\,dxdt\\
&=\lambda  M_i (\frac{2}{b+2})^2 \|\nabla (( u_i^{\eps})^{\frac{b+2}{2}})\|_{\LQ 2}^2.
\end{split}
\end{equation*}

In particular

\begin{equation}\label{eq-p(1)-es}
\begin{split}
\|\nabla ((u_i^{\eps})^{\frac{b+2}{2}})\|_{\LQ 2}\leq C\left(T, |\Omega|, \|u_{i,0}\|_{\LO{\infty}}\right).
\end{split}
\end{equation}

By testing the equation (<ref>) with $\varphi\in L^2(0,T; H^1(\Omega))$,
we have

\begin{equation*}
\begin{split}
&\int^T_0\intO\pa_tu_i^{\eps} \varphi\,dxdt=-\int^T_0\intO D_i(x,t) \Phi(u^{\eps}) \nabla u_i^{\eps} \cdot\nabla \varphi \,dxdt+\int^T_0\intO f_i^{\eps}(u^{\eps}) \varphi \,dxdt\\
&=-\int^T_0\intO D_i(x,t) \Phi(u^{\eps})^{\frac{1}{2}} (1+u_i^{\eps})^{\frac{1+\beta}{2}} \frac{ \Phi(u^{\eps})^{\frac{1}{2}}\nabla u_i^{\eps}}{(1+u_i^{\eps})^{\frac{1+\beta}{2}}} \cdot\nabla \varphi \,dxdt +\int^T_0\intO f_i^{\eps}(u^{\eps}) \varphi \,dxdt\\
&\leq \|D_i(x,t)\|_{\LQ{\infty}} \|\Phi(u^{\eps})\|^{\frac{1}{2}}_{\LQ{\infty}} (1+\|u_i^{\eps}\|_{\LQ{\infty}})^{\frac{1+\beta}{2}} \|\frac{\Phi(u^{\eps})^{\frac{1}{2}}|\nabla u_i^{\eps}|}{(1+u_i^{\eps})^{\frac{1+\beta}{2}}}\|_{\LQ2} \|\nabla\varphi\|_{\LQ2}\\
& \quad +\|f_i^{\eps}(u^{\eps})\|_{\LQ{\infty}} (|\Omega| T)^{\frac 1 2}\|\varphi\|_{\LQ2}\\
&\leq C\left(T, |\Omega|,\|u_{i,0}\|_{\LO{\infty}}, \|D_i(x,t)\|_{\LQ{\infty}}\right) \|\varphi\|_{L^2(0,T; H^1(\Omega))},
\end{split}
\end{equation*}

where we use (<ref>), (<ref>), (<ref>) and (<ref>).
Thus we can get

\begin{equation}\label{eq-pa-t}
\|\pa_t u_i^{\eps}\|_{L^2(0,T; (H^1(\Omega))')}\leq C\left(T, |\Omega|, \|u_{i,0}\|_{\LO{\infty}}\right).
\end{equation}

From Lemma <ref>, (<ref>) and (<ref>), we can  apply a nonlinear version of the well known Aubin-Lions Lemma (see e.g. <cit.>) to ensure

\begin{equation}\label{eq-strong-conv}
u_i^{\eps}\stackrel{\eps\rightarrow 0}{\longrightarrow}u_i \quad \text{strongly in } \LQ2
\end{equation}

Thanks to the $L^{\infty}$ bound of Lemma <ref>, this convergence in fact holds in $\LQ p$ for any $1\leq p<\infty$.

Since the sequence $(u_i^\eps)^{\frac{b+2}{2}}$ is uniformly bounded in $L^2(0,T; W^{1,2}(\Omega))$, there exists a subsequence (still denoted by $(u_i^\eps)$) and $v\in \LQ 2$ such that

\begin{align*}
\nabla (u_i^\eps)^{\frac{b+2}{2}} \rightharpoonup v \quad \text{weekly in } \LQ 2.
\end{align*}

Now, we are ready to identify the weak limit $v$ by the following computation: for testing functions $\varphi \in C^{\infty}_0(Q_T)$ we have

\begin{align*}
\intQT v \cdot \varphi \,dxdt &= \lim_{\eps\rightarrow 0} \intQT \nabla (u_i^\eps)^{\frac{b+2}{2}} \cdot \varphi \,dxdt\\
& = -  \lim_{\eps\rightarrow 0} \intQT  (u_i^\eps)^{\frac{b+2}{2}} \cdot \nabla \varphi \,dxdt\\
& = -  \intQT  u_i^{\frac{b+2}{2}} \cdot \nabla \varphi \,dxdt.
\end{align*}

This implies that the function $u_i^{\frac{b+2}{2}}$ is weakly differentiable with weak derivative $v$, i.e., $\nabla u_i^{\frac{b+2}{2}} = v$.

\begin{align}\label{eq-weak-convergence-nabla}
\nabla (u_i^\eps)^{\frac{b+2}{2}} \rightharpoonup \nabla u_i^{\frac{b+2}{2}} \quad \text{weekly in } \LQ 2.
\end{align}

It remains to pass  to the limit $\eps\rightarrow 0$ in the weak formulation of (<ref>) with $\eta\in L^2(0,T; H^1(\Omega))$,

\begin{equation}\label{eq-ori-limit}
\begin{split}
\int^T_0\langle \pa_t u_i^{\eps}, \eta\rangle_{(H^1(\Omega))', H^1(\Omega)} \,dt &+\int^T_0\intO D_i(x,t) \Phi(u^{\eps})  \nabla u_i^{\eps}\cdot \nabla \eta \,dxdt\\
&=\int^T_0\intO f_i^{\eps}(u^{\eps}) \eta \,dxdt.
\end{split}
\end{equation}

The convergence of the first term on the left hand side and the term on the right hand side of (<ref>) is immediate.

The convergence

\begin{equation}\label{eq-conv-term2}
\begin{split}
\int^T_0\intO D_i(x,t)  \Phi(u^{\eps})  \nabla u_i^{\eps}\cdot \nabla \eta \,dxdt \stackrel{\eps\rightarrow 0}{\longrightarrow} \int^T_0\intO D_i(x,t) \Phi(u) \nabla u_i\cdot \nabla \eta \,dxdt \quad \text{ in } Q_T
\end{split}
\end{equation}

follows from $\frac{\Phi(u^{\eps})}{(u^{\eps}_i)^{\frac{b}{2}}}   \stackrel{\eps\rightarrow 0}{\longrightarrow} \frac{\Phi(u)}{u_i^{\frac{b}{2}} } $ a.e. in $Q_T$, Lemma <ref> and the equation (<ref>).

Passing to the limit in the weak formulation (<ref>), we have

\begin{equation}\label{eq-limit}
\begin{split}
\int^T_0 \intO \pa_t u_i \eta \,dxdt  &+\int^T_0\intO D_i(x,t) \Phi(u) \nabla u_i\cdot \nabla \eta \,dxdt=\int^T_0\intO f_i(u) \eta \,dxdt.
\end{split}
\end{equation}

We obtain that $u=\left(u_1,\cdots,u_m\right)$ is a global weak solution to (<ref>) and additionally

\begin{equation*}
\|u_i\|_{\LQ{\infty}}\leq C_T, \quad \forall i=1,\cdots,m.
\end{equation*}

§ UNIFORM-IN-TIME BOUNDNESS

Assume (<ref>), (<ref>) and (<ref>) with either $K_1<0$ or $K_1=K_2=0$. Then, there exists a constant M independent of time such that

\begin{equation}\label{eq-uni-l1}
\sup_{t\geq 0}\|u_i(t)\|_{\LO1}\leq M:=\|u_{i0}\|_{\LO1}, \quad \forall i=1,\cdots,m.
\end{equation}

Similar to the proof of Lemma <ref>, we have

\begin{equation*}
\frac{d}{dt}\sumi\intO c_iu_i \,dx \leq \intO \left(K_1 \sumi u_i^{\eps}(x,t) + K_2\right) \,dx.
\end{equation*}

Integrating over $(s,t)$ we have

\begin{equation*}
\sumi c_i\intO u_i(x,t) \,dx \leq \sumi c_i\intO u_i(x,s) \,dx +\int^t_s\intO K_1 \sumi u_i(x,r) \,dx dr+ K_2 |\Omega|(t-s)
\end{equation*}

for all $t>s\geq0$.

If $K_1=K_2=0$, the bound (<ref>) follows immediately.

If $K_1<0$, we get for some constant $\sigma>0$ and for all $t>s\geq0$,

\begin{equation}\label{eq-mod}
\begin{split}
\sumi c_i\intO u_i(x,t) \,dx &+ \sigma \int^t_s \left(\sumi c_i\intO  u_i(x,r) \,dx\right) dr \\
&\leq \sumi c_i\intO u_i(x,s) \,dx+ K_2 |\Omega|(t-s).
\end{split}
\end{equation}

\begin{equation*}
\psi(t)= \sumi c_i\intO u_i(x,t) \,dx ~\text{ and } ~\phi(s) =  \int^t_s \left(\sumi c_i\intO  u_i(x,r) \,dx\right) dr.
\end{equation*}

It follows from (<ref>) that

\begin{equation*}
\phi^{\prime}(s)=-\sum_{i=1}^m c_i \int_{\Omega} u_i(\cdot, s) d x \leq-\psi(t)-\sigma \phi(s)+K_2|\Omega|(t-s),
\end{equation*}

which leads to

\begin{equation*}
\left(e^{\sigma s} \phi(s)\right)^{\prime}+e^{\sigma s} \psi(t) \leq K_2|\Omega| e^{\sigma s}(t-s) .
\end{equation*}

Integrating with respect to $s$ on $(0, t)$, and using $\phi(t)=0$, we have

\begin{equation*}
-\phi(0)+\psi(t) \frac{e^{\sigma t}-1}{\sigma} \leq \frac{K_2|\Omega|}{\sigma}\left(-t+\frac{e^{\sigma t}-1}{\sigma}\right) .
\end{equation*}

Since $\sigma \phi(0)-K_2|\Omega| t \leq \sum_{i=1}^m c_i \int_{\Omega}u_{i 0}(x)\,d x$ (see (<ref>)) it follows that

\begin{equation*}
\psi(t) \leq\left(e^{\sigma t}-1\right)^{-1} \sum_{i=1}^m c_i \int_{\Omega} u_{i, 0}(x) d x+K_2|\Omega| \sigma^{-1},
\end{equation*}

which finishes the proof of  Lemma <ref>.

We first show that for any $1\leq p<\infty$, there exists a constant $C_p>0$ such that

\begin{equation}\label{eq-lp-uni}
\sup_{t\geq 0}\|u_i(t)\|_{\LO p}\leq C_p, \quad \forall i=1,\cdots,m.
\end{equation}

Indeed, by using the $L^p$-energy function $\mathscr{L}_p[u]$ defined in (<ref>) and similar to Lemma <ref>, we obtain

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\alpha_p \sum_{k=1}^m \int_{\Omega}u_k(x,t)^{b+p-2} \left|\nabla u_k(x, t)\right|^2 \,dx \leq C\left(\intO\sumi u_i(x,t)^{p-1+r}\,dx+1\right).
\end{split}
\end{equation*}

Combining Lemmas <ref> and <ref>, we get

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\frac{\alpha_p}{2} \sum_{k=1}^m \int_{\Omega}\left(u_k(x,t)^{b+p-2} \left|\nabla u_k(x, t)\right|^2 +u_k(x,t)^{b+p}\right)\,dx \leq C(p).
\end{split}
\end{equation*}

This implies

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t) + \sigma\mathscr{L}_p(t) \leq C(p).
\end{split}
\end{equation*}

Thus, we have

\begin{equation*}
\sup _{t\geq 0}\|u_i(t)\|_{L^p(\Omega)} \leq C_{p} \quad \forall i=1, \ldots, m.
\end{equation*}

Next we will show that the solution is bounded uniformly in time in sup norm. We use a smooth time-truncated function $\psi: \R\rightarrow [0,1]$ with

\begin{equation*}
\psi(s)=
\begin{cases}
0, &s\leq0\\
1, &s\geq 1
\end{cases}
\end{equation*}

and $0\leq \phi'\leq C$, and its shifted version $\psi_{\tau}(\cdot)=\psi(\cdot,-\tau)$ for any $\tau\in\mathbb{N}$.  Let $\tau\in\mathbb{N}$ be arbitrary. It is straightforward to show that since $u=(u_i)_{i=1,\cdots,m}$ is a weak solution to (<ref>), the function $\psi_{\tau}u=(\psi_{\tau}u_i)_{i=1,\cdots,m}$ is a weak solution to the following

\begin{equation*}
\begin{cases}
\pa_t (\psi_{\tau}u_i)-\nabla\cdot(D_i(x,t) \Phi(\psi_{\tau}u) \nabla (\psi_{\tau}u_i))=\psi'_{\tau}u_i+\psi_{\tau}f_i(u),  &x\in \Omega, t\in(\tau,\tau+2)\\
(D_i(x,t) \Phi(\psi_{\tau}u) \nabla (\psi_{\tau}u_i))\cdot \nu=0,   &x\in \pa \Omega, t\in(\tau,\tau+2)\\
(\psi_{\tau}u_i)(x,\tau)=0,   &x\in \Omega.
\end{cases}
\end{equation*}

Thanks to (<ref>) and the polynomial growth (<ref>), we have

\begin{equation*}
\psi_{\tau}' u_i+\psi_{\tau}f_i(u)\leq G_i(x,t,u):=C(1+\sum^m_{k=1} u_k^l)\in L^p(\Omega\times(\tau,\tau+2)),
\end{equation*}

where $1\leq p<\infty$.
Therefore, similar to the proof of Lemma <ref>, we get

\begin{equation*}
\|\psi_{\tau}u_i\|_{L^{\infty}(\Omega\times(\tau,\tau+2)}\leq C, \quad \forall i=1,\cdots,m,
\end{equation*}

where C is a constant independent of $\tau\in \mathbb{N}$. Thanks to $\psi_{\tau}\geq 0$ and $\psi|_{(\tau+1,\tau+2)}\equiv1$, we obtain finally the uniform-in-time bound

\begin{equation*}
\sup_{\tau\in \mathbb{N}}\|u_i\|_{L^\infty(\Omega\times(\tau,\tau+1))}\leq C,\quad \forall i=1,\cdots,m,
\end{equation*}

and the proof of Theorem <ref> is complete.

§ PROOF OF THEOREMS <REF>–<REF>

For the global existence, we only need to show that for any $1\leq p<\infty$ and any $T>0$, there exists $C_{T,p}>0$ such that

\begin{equation}\label{p-1}
\sup _{t \in(0, T)}\|u_i^{\eps}(t)\|_{L^p(\Omega)} \leq C_{T,p}, \quad \forall i=1,\cdots,m.
\end{equation}

The rest follows exactly as in Lemma <ref>. To show (<ref>), we use the $L^p$-energy functions $\mathscr{L}_p(t)$ constructed in Lemma <ref> until (<ref>) we end up with

\begin{equation}\label{e-p-1}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\alpha_p \sum_{i=1}^m \int_{\Omega}\left(u_i^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_i(x, t)\right|^2 +u_i^{\eps}(x,t)^{b+p}\right)\,dx\\
& \leq C\left(\sumi\intO \left(u_i^{\eps}(x,t)^{p-1+r} + u_i^{\eps}(x,t)^{b+p}\right) \,dx+1\right).
\end{split}
\end{equation}

Since (<ref>) and (<ref>), we apply Lemma <ref> to estimate

\begin{equation*}
\begin{split}
&\intO\left(u_i^{\eps}(x,t)^{p-1+r} + u_i^{\eps}(x,t)^{b+p}\right) \,dx\\
&\leq  \frac{\alpha_p}{2}(\int_{\Omega} u_i^{\eps}(x,t)^{p-2+b}|\nabla u_i^{\eps}(x,t)|^{2}\,dx+\int_{\Omega} u_i^{\eps}(x,t)^{b+p}\,dx)+C_{T},
\end{split}
\end{equation*}

where $C_{T}$ depends on $\mathscr{F}(T)$. Inserting this into (<ref>) yields

\begin{equation*}
\begin{split}
\mathscr{L}_p^{\prime}(t)+\sigma \mathscr{L}_p(t) \leq C(T,p),
\end{split}
\end{equation*}

which implies (<ref>).

For the uniform-in-time bounds, we use an argument similar to the proof of Theorem <ref>.

For the global existence, we only need to show that for any $1\leq p<\infty$ and any $T>0$, there exists $C_{T,p}>0$ such that

\begin{equation}\label{p-2}
\sup _{t \in(0, T)}\|u_i^{\eps}(t)\|_{L^p(\Omega)} \leq C_{T,p}, \quad \forall i=1,\cdots,m.
\end{equation}

The rest follows exactly as in Lemma <ref>. To show (<ref>), we use the $L^p$-energy functions $\mathscr{L}_p(t)$ constructed in Lemma <ref> until (<ref>)  we obtain

\begin{equation}\label{e-p-2}
\begin{split}
\mathscr{L}_p^{\prime}(t)+&\frac{\alpha_p}{2} \sum_{i=1}^m \int_{\Omega}\left(u_i^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_i(x, t)\right|^2 +u_i^{\eps}(x,t)^{b+p}\right)\,dx\\
& \leq C\left(\sumi\intO u_i^{\eps}(x,t)^{p-1+r} \,dx+1\right).
\end{split}
\end{equation}

We integrate (<ref>) in time to obtain

\begin{equation}\label{eq-t3-1}
\begin{split}
\sup_{t\in(0,T)}\mathscr{L}_p(t)+&\frac{\alpha_p}{2} \sum_{i=1}^m \int^T_0\int_{\Omega}\left(u_i^{\eps}(x,t)^{b+p-2} \left|\nabla u^{\eps}_i(x, t)\right|^2 +u_i^{\eps}(x,t)^{b+p}\right)\,dxdt\\
& \leq C\sumi\int^T_0\intO u_i^{\eps}(x,t)^{p-1+r}  \,dxdt+C_{p,T}+\mathscr{L}_p(0).
\end{split}
\end{equation}

Denote by $y_i:=u_i^{\eps}(x,t)^{\frac{b+p}{2}}$. The left-hand side (LHS) of (<ref>) can be estimated below by

\begin{equation}
\begin{split}
C\sumi \left( \|y_i\|^{\frac{2p}{b+p}}_{L^{\infty}(0,T; L^{\frac{2p}{b+p}}(\Omega))}+\|y_i\|^2_{L^2(0,T; H^1(\Omega))} \right)\leq \text{(LHS) of } \eqref{eq-t3-1}.
\end{split}\label{eq-t3-2}
\end{equation}

For the right-hand side of (<ref>), we first consider

\begin{equation}\label{eq-t3-3}
\begin{split}
\vt>p-1+r
\end{split}
\end{equation}

as a constant to be determined later. Of course we are only interested in the case when $p-1+r>q$, otherwise the right hand side of (<ref>) is bounded thanks to (<ref>). By the interpolation inequality we have

\begin{equation}\label{eq-t3-4}
\begin{split}
\intQT u^{\eps}_i(x,t)^{p-1+r}dxdt &= \|u^{\eps}_i(x,t)\|_{\LQ{p-1+r}}^{p-1+r}\\
&\leq \|u^{\eps}_i(x,t)\|_{\LQ{q}}^{\theta(p-1+r)}\|u^{\eps}_i(x,t)\|_{\LQ{\vt}}^{(1-\theta)(p-1+r)},
\end{split}
\end{equation}

where $\theta\in (0,1)$ satisfies

\begin{equation*}
\frac{1}{p-1+r} = \frac{\theta}{q} + \frac{1-\theta}{\vt},
\end{equation*}

which implies

\begin{equation*}
(1-\theta)(p-1+r) = \frac{\vt(p-1+r-q)}{\vt - q}.
\end{equation*}

Using this, and taking into account (<ref>), (<ref>) implies

\begin{equation}\label{eq-t3-5}
\begin{aligned}
\intQT u^{\eps}_i(x,t)^{p-1+r}dxdt &\leq \mathscr{F}(T)^{\theta(p-1+r)}\|u^{\eps}_i(x,t)\|_{\LQ{\vt}}^{\frac{\vt(p-1+r-q)}{\vt - q}}\\
&= \mathscr{F}(T)^{\theta(p-1+r)}\bra{\intQT u^{\eps}_i(x,t)^{\vt}dxdt}^{\frac{p-1+r-q}{\vt -q}}\\
&= C_T\bra{\intQT y_i^{\frac{2\vt}{b+p}}dxdt}^{\frac{p-1+r-q}{\vt - q}}.
\end{aligned}
\end{equation}

For $p$ large enough, we choose $\vt$ close enough to (but bigger than) $p-1+r$ so that $H^1(\Omega)\hookrightarrow \LO{\frac{2\vt}{b+p}}$. That means $\vt$ is arbitrary for $N\leq 2$ and

\begin{equation}\label{eq-t3-6}
\frac{2\vt}{b+p} \leq \frac{2N}{N-2} \Leftrightarrow \vt \leq \frac{(b+p) N}{N-2} \quad\text{ for }\quad N\geq 3.
\end{equation}

Thus, we can use the Gagliardo-Nirenberg's inequality to estimate

\begin{equation}\label{eq-t3-7}
\intO y_i^{\frac{2\vt}{b+p}}dx = \|y_i\|_{\LO{\frac{2\vt}{b+p}}}^{\frac{2\vt}{b+p}} \leq C\|y_i\|_{H^1(\Omega)}^{\alpha\cdot \frac{2\vt}{b+p}}\|y_i\|_{\LO{\frac{2p}{b+p}}}^{(1-\alpha)\cdot \frac{2\vt}{b+p}}
\end{equation}

where $\alpha\in (0,1)$ satisfies

\begin{equation*}
\frac{b+p}{2\vt} = \bra{\frac 12 - \frac 1N}\alpha + \frac{(1-\alpha)(b+p)}{2p}.
\end{equation*}

From this

\begin{equation*}
\alpha\cdot\frac{2\vt}{b+p} = \frac{2N(\vt-p)}{Nb+2p} \quad \text{ and }
\quad (1-\alpha)\cdot\frac{2\vt}{b+p} = \frac{2Np(b+p) - 2(N-2)p\vt}{(b+p)(Nb+2p)}.
\end{equation*}

Therefore, we obtain from (<ref>) that

\begin{equation*}
\intO y_i^{\frac{2\vt}{b+p}}dx \leq C\|y_i\|_{H^1(\Omega)}^{\frac{2N(\vt-p)}{Nb+2p}}
\|y_i\|_{\LO{\frac{2p}{b+p}}}^\frac{2Np(b+p) - 2(N-2)p\vt}{(b+p)(Nb+2p)}.
\end{equation*}

It follows that

\begin{align*}
\intQT y_i^{\frac{2\vt}{b+p}}dxdt &\leq C\int_0^T\|y_i\|_{H^1(\Omega)}^{\frac{2N(\vt-p)}{Nb+2p}}
\|y_i\|_{\LO{\frac{2p}{b+p}}}^{\frac{2Np(b+p) - 2(N-2)p\vt}{(b+p)(Nb+2p)}} \,dt\\
&\leq C\|y_i\|_{L^{\infty}(0,T;\LO{\frac{2p}{b+p}})}^{\frac{2Np(b+p) - 2(N-2)p\vt}{(b+p)(Nb+2p)}} \int_0^T\|y_i\|_{H^1(\Omega)}^{\frac{2N(\vt-p)}{Nb+2p}} \,dt.
\end{align*}

We choose

\begin{equation}\label{eq-t3-8}
\frac{2N(\vt-p)}{Nb+2p} \leq 2 \Leftrightarrow \vt \leq \frac{Np+Nb+2p}{N}=b+p+\frac{2p}{N},
\end{equation}

which is possible since $p - 1 + r < p + b+ \frac{2p}{N}$ for $p$ large enough. Thus, by Hölder's inequality,

\begin{equation*}
\intQT y_i^{\frac{2\vt}{b+p}}dxdt \leq C_T\|y_i\|_{L^{\infty}(0,T;\LO{\frac{2p}{b+p}})}^{\frac{2Np(b+p) - 2(N-2)p\vt}{(b+p)(Nb+2p)}}
\|y_i\|_{L^2(0,T;H^1(\Omega))}^{\frac{2N(\vt-p)}{Nb+2p}}.
\end{equation*}

Inserting this into (<ref>) yields

\begin{equation}\label{eq-t3-9}
\begin{split}
\intQT &u_i^{\eps}(x,t)^{p-1+r}dxdt \leq C_T\bra{\|y_i\|_{L^{\infty}(0,T;\LO{\frac{2p}{b+p}})}^{\frac{2Np(b+p) - 2(N-2)p\vt}{(b+p)(Nb+2p)}}
\|y_i\|_{L^2(0,T;H^1(\Omega))}^{\frac{2N(\vt-p)}{Nb+2p}}}^{\frac{p-1+r-q}{\vt - q}}
\end{split}
\end{equation}

where $\theta\in(0,1).$

Using Young's inequality of the form

\begin{equation*}
x^{\lambda_1}y^{\lambda_2} \leq \eps(x^{\frac{2p}{b+p}}+y^2) + C_{\eps} \quad \text{ for } \quad \lambda_1\frac{b+p}{2p}+\frac{\lambda_2}{2} < 1,
\end{equation*}

we can choose $\vt$ such that

\begin{equation}\label{eq-t3-11}
\begin{split}
\frac{N(b+p) - (N-2)\vt}{(Nb+2p)}\frac{p-1+r-q}{\vt - q} +
\frac{N(p-1+r-q)(\vt - p)}{(Nb+2p)(\vt - q)}<1.
\end{split}
\end{equation}

This is equivalent to

\begin{align}\label{eq-t3-12}
\vt > \frac{Nb(p-1+r)+2pq}{Nb+2(1+q-r)}.
\end{align}

We now check that we can choose $\vt$ which satisfies all the conditions (<ref>), (<ref>), (<ref>) and (<ref>). This is fulfilled provided

\begin{equation*}
\frac{Nb(p-1+r)+2pq}{Nb+2(1+q-r)}<b+p+\frac{2p}{N}, \quad\text{which is equivalent to}\quad r < 1 + \frac{N}{N+2}b +  \frac{2q}{N+2},
\end{equation*}

and this can be obtain from our assumption (<ref>).
We estimate (<ref>) further as

\begin{equation*}
\intQT u_i^{\eps}(x,t)^{p-1+r}dxdt \leq \eps\bra{\|y_i\|_{L^{\infty}(0,T;\LO{\frac{2p}{b+p}})}^{\frac{2p}{b+p}} + \|y_i\|_{L^2(0,T;H^1(\Omega))}^2} + C_{T,\eps}.
\end{equation*}

\begin{equation*}
\text{RHS of (\ref{eq-t3-1})} \leq \mathscr{L}_p(0) + \eps\sumi \bra{\|y_i\|_{L^{\infty}(0,T;\LO{\frac{2p}{b+p}})}^{\frac{2p}{b+p}} + \|y_i\|^2_{L^2(0,T;H^1(\Omega))}} + C_{T,\eps}.
\end{equation*}

Combining this with (<ref>) gives us the desired estimate (<ref>).

For the uniform-in-time bounds, We use an argument similar to the proof of Theorem <ref>.

We give a formal proof as its rigor can be easily obtained through approximation. Since $h_k(\cdot)$ is convex, we have

\begin{align*}
\nabla\cdot(D_i\Phi(u)\nabla(h_i(u_i))) &= \nabla\cdot(D_i\Phi(u)h_i^{\prime}(u_i)\nabla u_i)\\
&\geq h'(u_i)\nabla\cdot(D_i\Phi(u)\nabla u_i) + \lambda h_i''(u_i)\Phi(u)|\nabla u_i|^2\\
&\geq h'(u_i)\nabla \cdot(D_i\Phi(u)\nabla u_i).
\end{align*}

Therefore, by defining $v_i:= h_i(u_i)\geq 0$ we have

\begin{equation}\label{ee2}
\pa_t v_i - \nabla\cdot(D_i\Phi(u)\nabla v_i) \leq G_i(x,t,u):= h_i^{\prime}(u_i)f_i(x,t,u), \quad x\in\Omega, \; t>0,
\end{equation}

with initial data $v_i(x,0) = h_i(u_{i,0}(x))$ and homogeneous Dirichlet boundary condition $v_i(x,t) = 0$ for $x\in\pa\Omega$ and $t>0$. Thanks to (<ref>) and  we have

\begin{equation*}
\sumi G_i(x,t,u) \leq K_5\sumi v_i + K_6
\end{equation*}

\begin{equation*}
A\begin{pmatrix}G_i(x,t,u)\\ \cdots\\ G_m(x,t,u)\end{pmatrix} \leq K_7\overrightarrow{1}\bra{\sumi v_i + 1}^r.
\end{equation*}

We can now reapply the methods in the proof of Theorem <ref> to obtain $v_i \in L^\infty_{\text{loc}}(0,\infty;\LO{\infty})$, and in case $K_5<0$ or $K_5 = K_6 = 0$ in (<ref>),

\begin{equation*}
\text{ess}\sup_{t\geq 0}\|v_i(t)\|_{\LO{\infty}} <+\infty, \quad \forall i=1,\ldots, m.
\end{equation*}

Due to the assumption (<ref>), the global existence and boundedness of (<ref>) immediately hold.

A vector of non-negative concentrations $u = (u_1, \ldots, u_m)$ is called a weak solution to (<ref>) on $(0,T)$ if

\begin{equation*}
u_i\in C([0,T]; \LO{2}), ~~ \Phi(u)\nabla u_i\in L^2(0,T;L^2(\Omega)), \quad f_i(u)\in L^2(0,T;\LO{2}),
\end{equation*}

with $u_i(\cdot,0)= u_{i,0}(\cdot)$ for all $i=1,\ldots, m$, and for any test function $\varphi \in L^2(0,T;H^1(\Omega))$ with $\pa_t\varphi \in L^2(0,T;H^{-1}(\Omega))$, it holds that

\begin{equation}\label{weak_bc}
\begin{aligned}
\intO u_i(x,t)\varphi(x,t)dx\bigg|_{t=0}^{t=T} - \intQT u_i\pa_t\varphi  \,dxdt + \intQT D_i(x,t)\Phi(u)\nabla u_i \cdot \nabla \varphi \,dxdt\\
+\alpha_i\int_0^T\int_{\pa\Omega}u_i\varphi d\mathscr{H}^{n-1} \,dt = \intQT f_i(x,t,u)\varphi \,dxdt.
\end{aligned}
\end{equation}

The proof of this theorem is similar to that of Theorems <ref>, <ref> and <ref>, except for the fact that the $L^1$-norm can be obtained in a different, and easier, way. Since $\varphi \equiv 1$ is an admissible test function we get from (<ref>) and (<ref>) that

\begin{align*}
\sumi \intO c_iu_i(\cdot,t)  \,dx\bigg|_{t=0}^{t=T} + \sumi c_i\alpha_i\int_0^T\int_{\pa\Omega} u_id\mathscr{H}^{n-1} \,dt\\
\leq K_1\intQT u_i  \,dxdt + K_2|\Omega|T.
\end{align*}

The $L^1$-bound is uniform in time in case $K_1<0$ or $K_1 = K_2 =0$ (using a similar idea to Lemma <ref>).

§ APPLICATIONS

We show the application of our results to SEIR(-D) model considered in [1, 31, 30] and its variants. First, the SEIR model introduced in  [31, 30, 1] reads as

\begin{equation}\label{sys-model}
\begin{cases}
\pa_t s= \alpha n - (1 - A_0/n) \beta_i si - (1 - A_0/n) \beta_e se - \mu s + \nabla\cdot(n \nu_s \nabla s), &x\in\Omega\\
\pa_t e= (1 - A_0/n) \beta_i si + (1 - A_0/n) \beta_e se - \sigma e - \phi_e e - \mu e + \nabla\cdot(n \nu_e \nabla e), &x\in\Omega\\
\pa_t i = \sigma e - \phi_d ni - \phi_r i - \mu i + \nabla\cdot(n \nu_i \nabla i), &x\in\Omega\\
\pa_t r = \phi_r i + \phi_e e - \mu r + \nabla\cdot(n \nu_r \nabla r), &x\in\Omega,\\
\nabla s \cdot \eta=\nabla e \cdot \eta=\nabla i \cdot \eta=\nabla r \cdot \eta=0, &x\in\pa\Omega,
\end{cases}
\end{equation}

where the symbols $s$, $e$, $i$ and $r$ denote the susceptible-, exposed-, infected-, recovered populations, respectively, and $n:=s+e+i+r$ is the total living population, $\nu_s$, $\nu_e$, $\nu_i$, $\nu_r$, $\beta_i$, $\beta_e$, $\alpha$, $\mu$, $\phi$, $\sigma$, $\phi_d$, $\phi_r$ and $\phi_e$ are positive constants. The global existence of bounded weak solutions was shown in [1] where it was imposed that all diffusion rates are the same, i.e. $\nu_s = \nu_e = \nu_i = \nu_r$, the term $1-A_0/n$ is replaced by a non-singular function $A(n)$, for instance $A(n) = (1-A_0/n)_+$, and the term $-\phi_dni$ is replaced by $-\phi_d i$. These are needed to show the total population $n$ is pointwise bounded from below by a positive constant, making the diffusion operators non-degenerate[In fact, thanks to the lower bound of $n$, it was possible to replace all diffusion operators by a common function $\kappa(n)$, where $\kappa: (0,\infty)\to (0,\infty)$ is a continuous function.]. The method therein seems to break down when the diffusion are different. By applying our main results in Theorem <ref>, we show that these assumptions can be removed. We have the following result.

Assume all diffusion coefficients $\nu_s, \nu_e, \nu_i, \nu_r$, and other parameters in (<ref>) to be positive. Then, for any non-negative and bounded initial data $(s_0,e_0,i_0,r_0)\in L_+^{\infty}(\Omega)^4$, there is a global bounded weak solution to (<ref>), i.e. for any $T>0$,

\begin{equation}\label{model-ext}
\sup_{t\in (0,T)}\left(\|s(t)\|_{\LO{\infty}}+\|e(t)\|_{\LO{\infty}}+ \|i(t)\|_{\LO{\infty}}+ \|r(t)\|_{\LO{\infty}} \right) <\infty.
\end{equation}

Moreover, if  $\alpha\leq \mu$, then this solution is bounded uniformly in time, i.e.

\begin{equation}\label{model-uni}
\sup_{t>0}\big(\|s(t)\|_{\LO{\infty}}+ \|e(t)\|_{\LO{\infty}}+ \|i(t)\|_{\LO{\infty}}+ \|r(t)\|_{\LO{\infty}} \big) < +\infty.
\end{equation}

It is noted that the nonlinearities in (<ref>) are not locally Lipschitz continuous around $0$. We circumvent this problem by considering the following approximating system for each $\delta>0$,

\begin{equation}\label{sys-approximate-model}
\begin{cases}
\pa_t s^{\delta}= \alpha n^{\delta} - (1 - \frac{A_0}{n^{\delta} + \delta} ) \beta_i s^{\delta} i^{\delta} - (1 - \frac{A_0}{n^{\delta} + \delta} ) \beta_e s^{\delta} e^{\delta} - \mu s^{\delta} + \nabla\cdot(n^{\delta} \nu_s \nabla s^{\delta}), &x\in\Omega\\
\pa_t e^{\delta}= (1 - \frac{A_0}{n^{\delta} + \delta} ) \beta_i s^{\delta} i^{\delta} + (1 - \frac{A_0}{n^{\delta} + \delta} ) \beta_e s^{\delta} e^{\delta} - \sigma e^{\delta} - \phi_e e^{\delta} - \mu e^{\delta} + \nabla\cdot(n^{\delta} \nu_e \nabla e^{\delta}), &x\in\Omega\\
\pa_t i^{\delta} = \sigma e^{\delta} - \phi_d n^{\delta} i^{\delta} - \phi_r i^{\delta} - \mu i^{\delta} + \nabla\cdot(n^{\delta} \nu_i \nabla i^{\delta}), &x\in\Omega\\
\pa_t r^{\delta} = \phi_r i^{\delta} + \phi_e e^{\delta} - \mu r^{\delta} + \nabla\cdot(n^{\delta} \nu_r \nabla r^{\delta}), &x\in\Omega,\\
\nabla s^{\delta} \cdot\eta=\nabla e^{\delta} \cdot\eta=\nabla i^{\delta} \cdot\eta=\nabla r^{\delta} \cdot \eta=0, &x\in\pa\Omega,
\end{cases}
\end{equation}

where $u^\delta = (s^\delta, e^\delta, i^\delta, r^\delta)$, $n^\delta = s^\delta + e^\delta + i^\delta + r^\delta$, $f^{\delta}_s, f^{\delta}_e, f^{\delta}_i$ and $f^{\delta}_r$ the nonlinearities in the equations of $s^{\delta}, e^{\delta}, i^{\delta}$ and $r^{\delta}$, respectively.

\begin{equation*}
\begin{cases}
f^{\delta}_s(u^{\delta}) \leq \alpha n^{\delta}  + (A_0\beta_i + A_0\beta_e)s^{\delta} \\
f^{\delta}_s(u^{\delta}) + f^{\delta}_e(u^{\delta})\leq \alpha  n^{\delta}\\
f^{\delta}_s(u^{\delta}) + f^{\delta}_e(u^{\delta}) + f^{\delta}_i(u^{\delta}) \leq \alpha  n^{\delta}\\
f^{\delta}_s(u^{\delta}) + f^{\delta}_e(u^{\delta}) + f^{\delta}_i(u^{\delta}) + f^{\delta}_r(u^{\delta}) \leq (\alpha - \mu) n^{\delta}.
\end{cases}
\end{equation*}

Following from  Theorem <ref>, we can obtain the global existence of a weak solution for the system (<ref>) and

\begin{align}\label{model-n-infty}
\|s^{\delta}\|_{\LQ{\infty}}, \|e^{\delta}\|_{\LQ{\infty}}, \|i^{\delta}\|_{\LQ{\infty}}, \|r^{\delta}\|_{\LQ{\infty}}  \leq C_T,
\end{align}

where $C_T$ independent of $\delta$.
Due to the polynomial growth of $f^{\delta}_s, f^{\delta}_e, f^{\delta}_i$ and $f^{\delta}_r$, we can get

\begin{align*}
\|f^{\delta}_s (u^{\delta})\|_{\LQ{\infty}}, \|f^{\delta}_e (u^{\delta})\|_{\LQ{\infty}}, \|f^{\delta}_i (u^{\delta})\|_{\LQ{\infty}}, \|f^{\delta}_r (u^{\delta})\|_{\LQ{\infty}}  \leq C_T.
\end{align*}

By multiply the first equation of (<ref>) by $s^{\delta} $ then integrating on $\Omega_T$ gives

\begin{align}\label{b1}
\sup_{t\in(0,T)}\intO s^{\delta}(t)^2 + \nu_s \intQT n^{\delta}|\nabla s^{\delta}|^{2}dxdt \leq C_T.
\end{align}

Since $n^\delta \ge s^\delta$, it follows

\begin{align}\label{model-nabla-s}
C_T \ge \intQT s^\delta|\nabla s^\delta|^2dxdt \ge C \|\nabla  (s^{\delta})^{\frac 32} \|_{\LQ{2}}^2.
\end{align}

By testing the first equation of (<ref>) with $\phi \in L^2(0,T; H^1(\Omega))$, we have

\begin{align*}
&\intQT \pa_t s^{\delta} \phi \,dxdt = - \nu_s\intQT n^{\delta}  \nabla s^{\delta} \cdot\nabla \phi \,dxdt + \intQT f^{\delta}_s(u^{\delta}) \phi \,dxdt\\
&\le \nu_s\|\sqrt{n^\delta}\|_{\LQ{\infty}}\| \sqrt{n^\delta}|\nabla s^\delta|\|_{\LQ{2}}\|\nabla \phi\|_{\LQ{2}} + \|f_s^{\delta}(u^\delta)\|_{\LQ{2}}\|\phi\|_{\LQ{2}}\\
&\leq C_T \|\phi\|_{L^2(0,T; H^1(\Omega))},
\end{align*}

where we used (<ref>) and (<ref>). Thus,

\begin{align*}
\|\pa_t s^{\delta}\|_{L^2(0, T;(H^1(\Omega))')} \leq C_T.
\end{align*}

Applying a nonlinear version of Aubin-Lions Lemma, see e.g. <cit.> to ensure

\begin{equation}\label{eq-strong-conv-s^delta}
s^{\delta}\stackrel{\delta\rightarrow 0}{\longrightarrow} s \quad \text{strongly in } \LQ2
\end{equation}

Thanks to the $L^{\infty}$ bound of $s^{\delta}$, this convergence in fact holds in $\LQ p$ for any $1\leq p<\infty$. By the inequality (<ref>),
we obtain

\begin{align*}
\nabla s^\delta \rightharpoonup \nabla s, \quad \text{weakly in } \LQ 2.
\end{align*}

Similarly, for  $e^{\delta}, i^{\delta}$ and $r^{\delta}$, we can get the strong convergence in $\LQ{p}$ for any $1\leq p<\infty$ and weak convergence in $L^2(0, T; H^1(\Omega))$. Now, let $\delta\rightarrow0$ in the weak formulation of (<ref>) with $\psi\in L^2(0, T; (H^1(\Omega))')$ and using <cit.>, we get the weak solution $u= (s, e, i, r)$ of (<ref>).

Furthermore, if  $\alpha \leq \mu$, we have  the uniform-in-time boundedness first for solutions to (<ref>) and consequently for (<ref>).

Thanks to the robustness of our approach, we can also consider a variant (<ref>) where

* $\beta_i$ and $\beta_e$ are functions of $n$ or the combinations of $n$ and  species $i$ or $e$ (see the original model proposed in [31]);

* the diffusion rates are functions of $(x,t)$, which shows the possibly high heterogeneity of the environment; and

* and the other rates $\alpha, \mu, \sigma, \phi_e, \phi_d, \phi_r$ are non-negative functions of $(x,t)$,

and therefore the system reads as

\begin{equation}\label{sys-model-2}
\begin{cases}
\pa_t s= \alpha n - (1-\frac{A_0}{n}) \beta_i(i,n) si - (1-\frac{A_0}{n}) \beta_e(e,n) se - \mu s + \nabla\cdot(n \nu_s(x,t) \nabla s),\\
\pa_t e= (1-\frac{A_0}{n}) \beta_i(i,n) si + (1-\frac{A_0}{n}) \beta_e(e,n) se - \sigma e - \phi_e e - \mu e + \nabla\cdot(n \nu_e(x,t) \nabla e),\\
\pa_t i = \sigma e - \phi_dn i - \phi_r i - \mu i + \nabla\cdot(n \nu_i (x,t)\nabla i),\\
\pa_t r = \phi_r i + \phi_e e - \mu r + \nabla\cdot(n \nu_r(x,t) \nabla r), \\
\nabla s \cdot\eta=\nabla e \cdot\eta=\nabla i \cdot\eta=\nabla r \cdot\eta=0.
\end{cases}
\end{equation}

We assume the following:

*  $\beta_i(i,n), \beta_e(e,n)$ are non-negative continuous functions and satisfy

\begin{equation*}
\beta_i(i,n)  \leq c_i \bra{1 + i + n}, \quad \beta_e(e,n) \le c_e\bra{1 + e+n},
\end{equation*}

for some constants $c_i, c_e>0$,

*  for $z\in \{s,e,i,r\}$, $\nu_z: \Omega\times[0,\infty) \to \R^{n\times n}$ such that $\nu_z\in L^\infty(\Omega\times(0,T);\R^{n\times n})$ for each $T>0$, and there is a constant $\lambda_z>0$ such that

\begin{equation*}
\lambda_z|\xi|^2 \le \xi^\top \nu_z(x,t) \xi, \quad \forall (x,t)\in\Omega\times[0,\infty), \quad \forall \xi \R^n;
\end{equation*}

*  all non-negative functions $\alpha, \mu, \sigma, \phi_e, \phi_d, \phi_r: \Omega\times[0,\infty)$ are bounded by a common constant, i.e. $\exists M>0$ such that

\begin{equation*}
(\alpha + \mu + \sigma + \phi_e + \phi_d+\phi_r)(x,t) \le M, \quad \forall (x,t)\in \Omega\times [0,\infty).
\end{equation*}

The following theorem follows directly from Theorem <ref>.

Assume (<ref>)–(<ref>). Then, for any non-negative bounded initial data, there exists a  global bounded  weak solution to (<ref>), i.e. for any $T>0$,

\begin{equation*}
\sup_{t\in (0,T)}\big(\|s(t)\|_{\LO{\infty}}+\|e(t)\|_{\LO{\infty}}+ \|i(t)\|_{\LO{\infty}}+\|r(t)\|_{\LO{\infty}} \big) < +\infty.
\end{equation*}

In particular, if  $\alpha(x,t)\leq \mu(x,t)$ for all $(x,t)\in \Omega\times[0,\infty)$, then the solution is bounded uniformly in time, i.e.

\begin{equation*}
\sup_{t>0}\big(\|s(t)\|_{\LO{\infty}}+ \|e(t)\|_{\LO{\infty}}+ \|i(t)\|_{\LO{\infty}}+ \|r(t)\|_{\LO{\infty}}\big) < +\infty.
\end{equation*}

§ PROOF OF LEMMA <REF>

We first state the following interpolation inequality.

Let $\gamma\geq 1$ and $\alpha>0$. Then we have the following interpolation inequality

\begin{equation}\label{eq-key2-1}
\|u\|^{\lambda}_{\LQ{\sigma}}\leq C \left( \|u\|^{\gamma}_{L^{\infty}(0,T;L^{\gamma}(\Omega))} + \|u\|^{\alpha}_{L^{\alpha}(0,T;L^{\beta}(\Omega))} \right),
\end{equation}

\begin{equation}\label{eq-key2-2}
\begin{cases}
\lambda=\frac{\alpha\beta+\gamma(\beta-\alpha)}{2\beta-\alpha} \quad \text{and} \quad   \sigma=\frac{\alpha\beta+\gamma(\beta-\alpha)}{\beta}, &\text{ if } \quad 0<\beta<\infty,\\
\lambda=\frac{\alpha+\gamma}{2} \quad \text{and} \quad   \sigma=\alpha+\gamma,   &\text{ if }\quad \beta=\infty.
\end{cases}
\end{equation}

For $0<\beta<\infty$. We apply an interpolation inequality to get

\begin{equation}\label{eq-key2-3}
\|u\|_{\LQ{\sigma}}\leq \|u\|^{\theta}_{L^{\infty}(0,T;L^{\gamma}(\Omega))}  \|u\|^{1-\theta}_{L^{\alpha}(0,T;L^{\beta}(\Omega))}
\end{equation}

with $\theta\in(0,1)$ satisfying

\begin{equation}
\frac{1}{\sigma}=\frac{\theta}{\infty} + \frac{1-\theta}{\alpha}=\frac{\theta}{\gamma} + \frac{1-\theta}{\beta}.
\end{equation}

It follows that

\begin{equation}\label{eq-key2-4}
\theta=\frac{\gamma(\beta-\alpha)}{\alpha\beta+\gamma(\beta-\alpha)}, \quad
1-\theta=\frac{\alpha\beta}{\alpha\beta + \gamma(\beta-\alpha)} \quad \text{and} \quad
\sigma=\frac{\alpha\beta + \gamma(\beta-\alpha)}{\beta}.
\end{equation}

From (<ref>) and $\lambda>0$ to be chosen, using Young's inequality we obtain

\begin{equation*}
\begin{split}
\|u\|^{\lambda}_{\LQ{\sigma}} & \leq \|u\|^{\theta\lambda}_{L^{\infty}(0,T;L^{\gamma}(\Omega))} \|u\|^{(1-\theta)\lambda}_{L^{\alpha}(0,T;L^{\beta}(\Omega))}\\
& \leq \frac{\theta\lambda}{\gamma} \|u\|^{\gamma}_{L^{\infty}(0,T;L^{\gamma}(\Omega))} + \frac{\gamma-\theta\lambda}{\gamma} \|u\|^{(1-\theta)\lambda\frac{\gamma}{\gamma-\theta\lambda}}_{L^{\alpha}(0,T;L^{\beta}(\Omega))}.
\end{split}
\end{equation*}

Now, we choose $\lambda$ such that

$$\frac{(1-\theta) \lambda \gamma}{\gamma-\theta\lambda}=\alpha,$$

which is equivalent to

where we used (<ref>).

\begin{equation*}
\|u\|^{\lambda}_{\LQ{\sigma}}\leq C \left(\|u\|^{\gamma}_{L^{\infty}(0,T;L^{\gamma}(\Omega))} + \|u\|^{\alpha}_{L^{\alpha}(0,T;L^{\beta}(\Omega))} \right)
\end{equation*}

\begin{equation*}
\lambda=\frac{\alpha\beta+\gamma(\beta-\alpha)}{2\beta-\alpha} \quad \text{and} \quad   \sigma=\frac{\alpha\beta+\gamma(\beta-\alpha)}{\beta}.
\end{equation*}

For $\beta=\infty$, we can use a similar method to get the result. This completes the proof of Lemma <ref>.

Let $\left\{y_{n}\right\}_{n \geq 1}$ be a sequence of positive numbers which satisfies

\begin{equation*}
y_{n+1} \leq K B^{n}\left(y_{n}^{\gamma}+y_{n}^{\kappa}\right)
\end{equation*}

where $K, B>0$ and $\gamma, \kappa>1$ are independent of $n$. Then there exists $\varepsilon>0$ such that, if $y_{1} \leq \varepsilon$, then

\begin{equation*}
\lim _{n \rightarrow \infty} y_{n}=0 .
\end{equation*}

We are now presenting a proof of Lemma <ref> using the idea of Moser iteration.

From (<ref>),

\begin{equation*}
f_i^{\eps}(u^{\eps})\leq f_i(u^{\eps})\leq F_i(u^{\eps}):=K_4\left(1+\sum^m_{j=1}(u_j^{\eps})^l\right).
\end{equation*}

Using Lemma <ref>, we can get for any $1\leq p<\infty$ a constant $C_p>0$ exists depending on $p$ (but not on $T$) such that

\begin{equation}\label{eq-inb}
\|f^{\eps}_i(u^{\eps})\|_{\LQ p}\leq C_p, \quad \forall i=1,\cdots,m.
\end{equation}

Let $k \geq 1$ be a constant which will be specified later. For each $j \geq 0$, we define

v_j:=\left(u^{\eps}_i-k+\frac{k}{2^j}\right)_{+}=\max \left\{u^{\eps}_i-k+\frac{k}{2^j},  0\right\}

A_j:=\left\{(x, t) \in Q_T: u^{\eps}_i(x, t) \geq k-\frac{k}{2^j}\right\} .

The following simple observations will be helpful

\begin{equation}\label{eq-inf-1}
\begin{aligned}
&v_{j+1}(x, t) \leq v_j(x, t), &&   \forall (x, t) \in A_j, \\
&v_j(x, t)  \geq \frac{k}{2^{j+1}}, && \forall (x, t) \in A_{j+1} \subset A_j.
\end{aligned}
\end{equation}

By multiplying the equation $\partial_t u^{\eps}_i - \Delta\left(D_i(x,t) \,\Phi(u^{\eps})\, \nabla u^{\eps}_i\right)=f^{\eps}_i$ by $v_{j+1}$ and integrating on $Q_T$, we have

\begin{equation*}
\begin{aligned}
& \sup _{t \in(0, T)}\|v_{j+1}(t)\|^2_{\LO2} + 2 \lambda \int_0^T \int_{\Omega}  \Phi(u^{\eps}) \left|\nabla v_{j+1}\right|^2 \,d x d t \\
& \quad \leq\|v_{j+1}(0)\|^2_{\LO2} + 2 \int_0^T \int_{\Omega} f^{\eps}_i v_{j+1} \,d x d t,
\end{aligned}
\end{equation*}

where we used (<ref>).
Using (<ref>), we have

\begin{equation}\label{eq-inf-2}
\begin{aligned}
& \sup _{t \in(0, T)}\|v_{j+1}(t)\|^2_{\LO2} + 2 \lambda M_i \int_0^T \int_{\Omega}(u^{\eps}_i)^{b}\left|\nabla v_{j+1}\right|^2 \,d x d t \\
& \quad \leq\|v_{j+1}(0)\|^2_{\LO2}+2 \int_0^T \int_{\Omega} f^{\eps}_i v_{j+1} \,d x d t .
\end{aligned}
\end{equation}

Note that $u^{\eps}_i \geq k-\frac{k}{2^j} \geq \frac{k}{2}$ on $A_j$ ($j\geq1$), we have

\begin{equation*}
\begin{aligned}
2\lambda M_i \int_0^T \int_{\Omega}(u^{\eps}_i)^{b}\left|\nabla v_{j+1}\right|^2 \,d x d t & \geq 2\lambda M_i \iint_{A_j}(u^{\eps}_i)^{b}\left|\nabla v_{j+1}\right|^2 \,d x d t \\
& \geq 2\lambda M_i \left(\frac{k}{2}\right)^{b} \iint_{A_j}\left|\nabla v_{j+1}\right|^2 \,dx dt \\
& \geq \frac{\lambda M_i}{2^{b-1}} \int_0^T \int_{\Omega}\left|\nabla v_{j+1}\right|^2 \,dxdt,
\end{aligned}
\end{equation*}

thanks to $k\geq1$  and the fact that $v_{j+1} \equiv 0$ on $Q_T \backslash A_{j+1} \supset Q_T \backslash A_j$ since $A_{j+1} \subset A_j$.
By adding $\frac{\lambda M_i}{2^{b-1}} \int_0^T\|v_{j+1}\|^2_{\LO2} d t$ to both sides of (<ref>), we get

\begin{equation*}
\begin{aligned}
& \sup _{t \in(0, T)}\|v_{j+1}(t)\|^2_{\LO2} + \frac{\lambda M_i}{2^{b-1}} \int_0^T\|v_{j+1}\|_{H^1(\Omega)}^2 \, d t \\
& \leq \frac{\lambda M_i}{2^{b-1}}  \int_0^T\|v_{j+1}\|^2_{\LO2}\, d t+\|v_{j+1}(0)\|^2_{\LO2} + 2\int_0^T \int_{\Omega} f^{\eps}_i v_{j+1} \,d x d t,
\end{aligned}
\end{equation*}

which yields

\begin{equation}\label{eq-inf-3}
\|v_{j+1}\|_{W(0, T)}^2 \leq C\left(\int_0^T\|v_{j+1}\|^2_{\LO2} \, d t+\|v_{j+1}(0)\|^2_{\LO2} + \int_0^T \int_{\Omega} f^{\eps}_i v_{j+1} \,d x d t \right).
\end{equation}

By definition, when we choose $k \geq 2\|u^{\eps}_0\|_{L^{\infty}(\Omega)}$, we have

\begin{equation}\label{eq-inf-4}
\|v_{j+1}(0)\|^2_{\LO2}=\|\left(u^{\eps}_0-k+\frac{k}{2^{j+1}}\right)_{+}\|^2_{\LO2}=0.
\end{equation}

By using (<ref>), we have with $1 \leq \frac{2^{j+1}}{k} v_j$ on $A_{j+1}$

\begin{equation}\label{eq-inf-5}
\begin{aligned}
\int_0^T \int_{\Omega}\left|v_{j+1}\right|^2 d x d t & =\int_0^T \int_{\Omega} \mathbf{1}_{A_{j+1}}\left|v_{j+1}\right|^2 \,d x d t \\
& \leq \int_0^T \int_{\Omega} \mathbf{1}_{A_{j+1}}\left|v_j\right|^2 \,d x d t \\
& \leq\left(\frac{2^{j+1}}{k}\right)^{\frac{4}{N}} \int_0^T \int_{\Omega} \mathbf{1}_{A_{j+1}}\left|v_j\right|^{2+\frac{4}{N}} \,d x d t \\
& \leq C\left(2^{\frac{4}{N}}\right)^j\|v_j\|_{W(0, T)}^{2+\frac{4}{N}}.
\end{aligned}
\end{equation}

Choose $p>\frac{N+2}{2}$, we have

\begin{equation*}
\sigma:=\frac{p-1}{p}\left(2+\frac{4}{N}\right)>2.
\end{equation*}

\begin{equation*}
\frac{\sigma p}{p-1}=2+\frac{4}{N},
\end{equation*}

\begin{equation*}
\|v_j\|_{L^{\frac{\sigma p}{p-1}}\left(Q_T\right)} \leq C\|v_j\|_{W(0, T)}.
\end{equation*}

We now can use Hölder's inequality to estimate with (<ref>)

\begin{equation}\label{eq-inf-6}
\begin{aligned}
\int_0^T \int_{\Omega} f^{\eps}_i v_{j+1}  \,d x d t & \leq \int_0^T \int_{\Omega} f^{\eps}_i v_{j+1}\left(\frac{2^{j+1}}{k}\right)^{\sigma-1} v_j^{\sigma-1} \,d x d t \\
& \leq\left(\frac{2^{j+1}}{k}\right)^{\sigma-1} \int_0^T \int_{\Omega}\left|f^{\eps}_i \| v_j\right|^\sigma \,d x d t \\
& \leq C\left(2^{\sigma-1}\right)^j\|f^{\eps}_i\|_{L^p\left(Q_T\right)}\|v_j\|^\sigma_{L^{\frac{\sigma p}{p-1}} \left(Q_T\right)} \\
& \leq C\left(2^{\sigma-1}\right)^j\|f^{\eps}_i\|_{L^p\left(Q_T\right)}\|v_j\|_{W(0, T)}^\sigma .
\end{aligned}
\end{equation}

Inserting (<ref>), (<ref>) and (<ref>) into (<ref>) leads to

\begin{equation}\label{eq-inf-7}
\|v_{j+1}\|_{W(0, T)}^2 \leq C\left(1+\|f^{\eps}_i\|_{L^p\left(Q_T\right)}\right) B^j\left(\|v_j\|_{W(0, T)}^{2+\frac{4}{N}}+\|v_j\|_{W(0, T)}^\sigma\right)
\end{equation}

for all $j \geq 0$, where $B=\max \left\{2^{\frac{4 }{N}} , 2^{\sigma-1}\right\}$.

By setting $Y_j=\|v_j\|_{W(0, T)}^2$, we obtain a sequence $\left\{Y_n\right\}_{n \geq 1}$ satisfying the property in Lemma <ref>. It remains to show that $Y_1$ is small enough.

We show now that for any $\eta>0$, there exists $k \geq \max \left\{1, 2\|u^{\eps}_0\|_{L^{\infty}(\Omega)}\right\}$ large enough such that

\begin{equation}\label{eq-inf-8}
Y_1=\|v_1\|^2_{W(0, T)} \leq \eta.
\end{equation}

Let $p> 1$. By multiplying (<ref>) by $p |u_i^{\eps}|^{p-1}$ and integrating over $\Omega$, we obtain

\begin{equation}\label{eq-e1}
\begin{split}
\frac{d}{dt}\|u_i^{\eps}\|^{p}_{\LO{p}} -p\int_{\Omega}\nabla\cdot(D(x,t) \Phi(u^{\eps}) \nabla u_i^{\eps}) |u_i^{\eps}|^{p-1} \,dx =p\int_{\Omega}f_i^{\eps} |u_i^{\eps}|^{p-1} \,dx.
\end{split}
\end{equation}

Integration by parts and the homogeneous Neumann boundary condition $(D(x,t) \Phi(u^{\eps}) \nabla u_i^{\eps})\cdot \nu= 0$ lead to

\begin{equation*}
\begin{aligned}
-p\int_{\Omega}&\nabla\cdot(D(x,t) \Phi(u^{\eps})  \nabla u_i^{\eps}) |u_i^{\eps}|^{p-1}\,dx\\
&\geq \lambda M_ip(p-1)\int_{\Omega}|u_i^{\eps}|^{b+p-2}|\nabla u_i^{\eps}|^2 \,dx\\
&=\lambda M_ip(p-1)\frac{4}{(b+p)^2}\int_{\Omega}|\nabla \left((u_i^{\eps})^{\frac{b+p}{2}}\right)|^2 \,dx\\
&=:C(p,\lambda, M_i)\int_{\Omega}|\nabla \left((u_i^{\eps})^{\frac{b+p}{2}}\right)|^2 \,dx,
\end{aligned}
\end{equation*}

where we used (<ref>) and (<ref>).
By Hölder's inequality

\begin{equation*}
\left|p\int_{\Omega}f_i^{\eps} |u_i^{\eps}|^{p-1} \,dx\right| \leq p\|f_i^{\eps}\|_{\LO p}\|u_i^{\eps}\|_{\LO{p}}^{p-1}.
\end{equation*}

Therefore, it follows from (<ref>) that

\begin{equation}\label{eq-mu-ori}
\frac{d}{dt}\|u_i^{\eps}\|_{\LO{p}}^{p} + C(p,\lambda, M_i) \int_{\Omega}\left|\nabla\left((u_i^{\eps})^{\frac{b+p}{2}}\right)\right|^2dx \leq p\|f_i^{\eps}\|_{\LO {p}}\|u_i^{\eps}\|_{\LO{p}}^{p-1}.
\end{equation}

By applying for $r<1$ the elementary inequality

\begin{equation*}
y' \leq \alpha(t)y^{1-r} \quad \Longrightarrow \quad y(T) \leq \left[y(0)^{r}+ r \int_0^T\alpha(t) \,dt\right]^{1/r},
\end{equation*}

to (<ref>) with $r=1/p$ and $y(t) = \|u_i^{\eps}(t)\|_{\LO p}^{p}$, we obtain

\begin{equation}\label{eq-CT}
\begin{split}
\|u_i^{\eps}(T)\|_{\LO p}^{p} &\leq \left[\|u^{\eps}_0\|_{\LO p} + \int_0^{T}\|f_i^{\eps}\|_{\LO p} \,dt \right]^{ p} \\
&\leq \left[\|u^{\eps}_0\|_{\LO p} + \|f_i^{\eps}\|_{\LQ{p}}T^{\frac{(p-1)}{p}} \right]^{p} =: C_{T}.
\end{split}
\end{equation}

That means

\begin{equation}\label{eq-p-1}
u_i^{\eps} \in L^{\infty}(0,T;L^{p}(\Omega)) \quad \text{ and } \quad \|u_i^{\eps}(T)\|_{\LO p}^{p} \leq C_{T}
\end{equation}

where $C_{T}$ defined in (<ref>) grows at most polynomially in $T$. By integrating (<ref>) with respect to $t$ on $(0,T)$ and by using Young's inequality and the convention $r := b+p\geq p>1$, we get

\begin{equation*}
\begin{aligned}
C(p,\lambda, M_i)\int^T_0\intO\left|\nabla\left((u_i^{\eps})^{\frac{r}{2}}\right)\right|^2 \,dxdt &\leq \|u_0^{\eps}\|_{\LO p}^{p} +
p\int_0^T\|f_i^{\eps}\|_{\LO p}\|u_i^{\eps}\|_{\LO p}^{p-1} \,dt\\
&\leq \|u_0^{\eps}\|_{\LO p}^{p} + p\|f_i^{\eps}\|_{L^{p}(Q_T)}\|u_i^{\eps}\|_{L^{p}(Q_T)}^{p-1}.
\end{aligned}
\end{equation*}

By adding $C(p,\lambda, M_i)\int^T_0\intO (u_i^{\eps})^{r}  \,d x d t$ to both sides, we have

\begin{equation}\label{eq-h0}
\begin{aligned}
C(p,\lambda, M_i)&\int_0^T\|(u_i^{\eps})^{\frac{r}{2}}\|_{H^1(\Omega)}^2 \,dt\\
& = C(p,\lambda, M_i)\int_0^T\left[
\int_{\Omega}\left|\nabla\left((u_i^{\eps})^{\frac{r}{2}}\right)\right|^2 \,dx + \int_{\Omega}\left| (u_i^{\eps})^{\frac{r}{2}} \right|^2  \,dx\right]dt\\
&\leq \|u_0^{\eps}\|_{\LO p}^{p} + p\|f_i^{\eps}\|_{L^{p}(Q_T)}\|u_i^{\eps}\|_{L^{p}(Q_T)}^{p-1} + C(p,\lambda, M_i)\int_0^T\|u_i^{\eps}\|_{\LO {r}}^{r} \,dt.
\end{aligned}
\end{equation}

By the Sobolev's embedding, we have

\begin{equation}\label{eq-h1}
C(p,\lambda, M_i)\int_0^T\|(u_i^{\eps})^{\frac{r}{2}}\|^2_{H^1(\Omega)}\,dt \geq C(p,\lambda, M_i) \,C_S^2\int_0^T
\|u_i^{\eps}\|_{\LO s}^{r} \,dt
\end{equation}

\begin{equation*}
s =
\begin{cases}
\frac{r N}{N-2} &\text{ if } N \geq 3,\\
r < s < \infty \text{ arbitrary } &\text{ if } N = 1,2.
\end{cases}
\end{equation*}

On the other hand, by using the bound $\|u_i^{\eps}(t)\|_{\LO p}^{p} \leq C_{T}$ in (<ref>) and the interpolation inequality

\begin{equation*}
\|u_i^{\eps}\|_{\LO {r}} \leq \|u_i^{\eps}\|_{\LO p}^{\alpha}\|u_i^{\eps}\|_{\LO s}^{1-\alpha}
\leq C_{T}^{\frac{\alpha}{p}}\|u_i^{\eps}\|_{\LO s}^{1-\alpha}
\end{equation*}

\begin{equation*}
\frac{1}{r} = \frac{\alpha}{p} + \frac{1-\alpha}{s} \quad \text{and} \quad \alpha = \frac{2p}{2p+bN}\in(0,1],
\end{equation*}

we estimate in the cases $b>0$ for which $\alpha<1$

\begin{equation}\label{eq-h2}
\begin{split}
C(p,\lambda, M_i)\int_0^T\|u_i^{\eps}\|_{\LO {r}}^{r} \,dt &\leq C(p,\lambda, M_i )\int_{0}^{T}C_{T}^{\frac{\alpha r}{p}} \|u_i^{\eps}\|_{\LO s}^{(1-\alpha)r} \,dt \\
&\leq \frac{C(p,\lambda, M_i)\,C_S^2}{2}\int_0^T\|u_i^{\eps}\|_{\LO s}^{r} \,dt + CC_{T}^{\frac{r}{p}}T,
\end{split}
\end{equation}

where we have used Young's inequality (with the exponents $1=(1-\alpha) + \alpha$) in the last step. Note that if $b=0$, the bound (<ref>) still holds true without the first term and with $\frac{r}{p}=1$.
Inserting (<ref>) and (<ref>) into (<ref>) leads to

\begin{equation}\label{eq-DT}
\begin{aligned}
\int_0^T\|u_i^{\eps}\|_{\LO s}^{r}dt &\leq \frac{2}{C(p,\lambda, M_i)\,C_S^2}\left[\|u_0^{\eps}\|_{\LO p}^{p} + p \|f_i^{\eps}\|_{L^{p}(Q_T)}\|u_i^{\eps}\|_{L^{p}(Q_T)}^{p-1}  + CC_{T,0}^{\frac{r}{p}}T\right]\\
&\leq \frac{2}{C(p,\lambda, M_i)\,C_S^2}\left[\|u_0^{\eps}\|_{\LO p}^{p} + p\|f_i^{\eps}\|_{L^{p}(Q_T)} C_{T}^{\frac{p-1}{p}}  + CC_{T}^{\frac{r}{p}}T\right] \\
&=: D_{T} \quad (\text{use (\ref{eq-p-1}))}.
\end{aligned}
\end{equation}

It follows that

\begin{equation*}
u^{\eps}_i \in L^{r}(0,T;L^{s}(\Omega))
\end{equation*}

\begin{equation*}
\begin{cases}
s = \frac{r N}{N-2} &\text{ if } N\geq 3,\\
r< s <\infty \text{ arbitrary } &\text{ if } N=1,2
\end{cases}
\end{equation*}

\begin{equation*}
\int_0^T\|u_i^{\eps}\|_{\LO s}^{r} \,dt \leq D_{T}
\end{equation*}

with $D_{T}$ defined in (<ref>).

Thus, we have

\begin{equation*}
\|u^{\eps}_i\|_{L^{\infty}\left(0, T ; L^p(\Omega)\right)}+\|u^{\eps}_i\|_{L^{r}\left(0, T ; L^s(\Omega)\right)} \leq C_T,
\end{equation*}

where $r=b+p \geq p$ and $s=\frac{r N}{N-2}$ if $N \geq 3$ and $r<s<+\infty$ arbitrary if $N \leq 2$.

Using Lemma <ref>, we have

\begin{equation*}
\|u^{\eps}_i\|_{L^{\tau}\left(Q_T\right)} \leq C_T \quad \text { with } \quad \tau=
\begin{cases}
\frac{N r +2 p}{N} & \text { if } N \geq 3, \\
<r +p \text { arbitrary } & \text { if } N \leq 2 .
\end{cases}
\end{equation*}

Direct calculations show that $\tau>2+\frac{4}{N}$ if $N \geq 2$ and $\tau>3$ if $N=1$. In particular,

\begin{equation}\label{eq-inf-9}
\|u^{\eps}_i\|_{L^{2+\frac{4}{N}}\left(Q_T\right)} \leq C_T \text { for } d \geq 2 \text { and }\|u^{\eps}_i\|_{L^3\left(Q_T\right)} \leq C_T \text { for } d=1 .
\end{equation}

From (<ref>),

\begin{equation}\label{eq-inf-10}
\|v_1\|_{W(0, T)}^2 \leq C\left(\int_0^T\|v_1(t)\|^2_{\LO2} \,d t+\|v_1(0)\|^2_{\LO2} + \int_0^T \int_{\Omega} f^{\eps}_i v_1 \,dxdt\right).
\end{equation}

Since $k \geq 2\|u^{\eps}_0\|_{L^{\infty}(\Omega)},\|v_1(0)\|^2_{\LO2}=\|\left(u^{\eps}_0-\frac{k}{2}\right)_+\|^2_{\LO2}=0$.

Consider now the case $N \geq 2$. By using (<ref>), it yields

\begin{equation}\label{eq-inf-11}
\begin{aligned}
\int_0^T \int_{\Omega}\left|v_1\right|^2 \,d x d t & =\int_0^T \int_{\Omega} \mathbf{1}_{A_1}\left|v_1\right|^2 \,d x d t \leq\left(\frac{4}{k}\right)^{\frac{4}{N}} \int_0^T \int_{\Omega}\left|v_0\right|^{2+\frac{4}{N}} \,d x d t \\
& \leq\left(\frac{4}{k}\right)^{\frac{4}{N}}\|u^{\eps}_i\|_{L^{2+\frac{4}{N}}(Q_T)}^{2+\frac{4}{N}} \leq\left(\frac{4}{k}\right)^{\frac{4}{N}} C_T,
\end{aligned}
\end{equation}

recalling that $v_0=(u^{\eps}_i)_{+}$. Similarly to (<ref>), we get

\begin{equation}\label{eq-inf-12}
\int_0^T \int_{\Omega} f^{\eps}_i v_1 \,d x d t \leq\left(\frac{4}{k}\right)^{\sigma-1}\|f^{\eps}_i\|_{L^p\left(Q_T\right)}\|u^{\eps}_i \|_{L^{2+\frac{4}{N}}\left(Q_T\right)}^\sigma \leq\left(\frac{4}{k}\right)^{\sigma-1} C_T .
\end{equation}

From (<ref>), (<ref>) and (<ref>), we get (<ref>) if

k=4 \max \left\{\left(\frac{C_T}{\eta}\right)^{\frac{N}{4}}, \left(\frac{C_T}{\eta}\right)^{\frac{1}{\sigma-1}}\right\}.

Thus, with this choice of $k$, it follows that

0=\lim _{j \rightarrow \infty} Y_j=\|(u^{\eps}_i-k)_{+}\|^2

and hence,

\|u^{\eps}_i\|_{L^{\infty}\left(Q_T\right)} \leq k=4 \max \left\{\left(\frac{C_T}{\eta}\right)^{\frac{N}{4}}, \left(\frac{C_T}{\eta}\right)^{\frac{1}{\sigma-1}}\right\}

which is our desired estimate.

The proof for the case $N=1$ is very similar using

\int_0^T \int_{\Omega}\left|v_1\right|^2 \,d x d t \leq \frac{4}{k} \int_0^T \int_{\Omega}\left|v_0\right|^3 \,d x d t \leq \frac{4}{k} C_T

\int_0^T \int_{\Omega} f^{\eps}_i v_1 \,d x d t \leq\left(\frac{4}{k}\right)^{\frac{4 \xi}{1+2 \xi}}\|f^{\eps}_i\|_{L^p\left(Q_T\right)}\|u^{\eps}_i\|_{L^3(Q_T)}^{1+\frac{4 \xi}{1+2 \xi}} \leq\left(\frac{4}{k}\right)^{\frac{4 \xi}{1+2 \xi}} C_T

where $\xi=\frac{1}{2}(2p-3)>0$. We therefore omit the details.

Thus, we  completed the proof of Lemma <ref>.

Acknowledgement. The authors would like to thank Prof. Pierluigi Colli, Prof. Gabriela Marinoschi, and Prof. Elisabetta Rocca for their inspired and fruitful discussion.

B.Q. Tang received funding from the FWF project “Quasi-steady-state approximation for PDE”, number I-5213. J. Yang are partially supported by NSFC Grants
No. 12271227 and China Scholarship Council (Contract No. 202206180025).

[1]
Ferdinando Auricchio, Pierluigi Colli, Gianni Gilardi, Alessandro Reali, and
Elisabetta Rocca.
Well-posedness for a diffusion–reaction compartmental model
simulating the spread of covid-19.
Mathematical Methods in the Applied Sciences, 46:12529–12548,

[2]
Dieter Bothe, André Fischer, Michel Pierre, and Guillaume Rolland.
Global wellposedness for a class of
reaction-advection-anisotropic-diffusion systems.
J. Evol. Equ., 17(1):101–130, 2017.

[3]
M Cristina Caputo, Thierry Goudon, and Alexis F Vasseur.
Solutions of the 4-species quadratic reaction-diffusion system are
bounded and ${C^{\infty}}$-smooth, in any space dimension.
Analysis & PDE, 12(7):1773–1804, 2019.

[4]
José A Cañizo, Laurent Desvillettes, and Klemens Fellner.
Improved duality estimates and applications to reaction-diffusion
Communications in Partial Differential Equations,
39(6):1185–1204, 2014.

[5]
M Cristina Caputo and Alexis Vasseur.
Global regularity of solutions to systems of reaction–diffusion with
sub-quadratic growth in any dimension.
Communications in Partial Differential Equations,
34(10):1228–1250, 2009.

[6]
Laurent Desvillettes, Klemens Fellner, Michel Pierre, and Julien Vovelle.
Global existence for quadratic systems of reaction-diffusion.
Advanced Nonlinear Studies, 7(3):491–511, 2007.

[7]
Julian Fischer, Katharina Hopf, Michael Kniely, and Alexander Mielke.
Global existence analysis of energy-reaction-diffusion systems.
SIAM Journal on Mathematical Analysis, 54(1):220–267, 2022.

[8]
Klemens Fellner, Evangelos Latos, and Takashi Suzuki.
Global classical solutions for mass-conserving, (super)-quadratic
reaction-diffusion systems in three and higher space dimensions.
Discrete Contin. Dyn. Syst. Ser. B, 21(10):3441–3462, 2016.

[9]
Klemens Fellner, Evangelos Latos, and Bao Quoc Tang.
Global regularity and convergence to equilibrium of
reaction-diffusion systems with nonlinear diffusion.
J. Evol. Equ., 20(3):957–1003, 2020.

[10]
Klemens Fellner, Jeff Morgan, and Bao Quoc Tang.
Global classical solutions to quadratic systems with mass control in
arbitrary dimensions.
Annales de l'Institut Henri Poincaré C, Analyse non
linéaire, 37(2):281–307, 2020.

[11]
Klemens Fellner, Jeff Morgan, and Bao Quoc Tang.
Uniform-in-time bounds for quadratic reaction-diffusion systems with
mass dissipation in higher dimensions.
Discrete and Continuous Dynamical Systems – Series S,
14(2):635–651, 2021.

[12]
William E. Fitzgibbon, Jeffrey J. Morgan, Bao Q. Tang, and Hong-Ming Yin.
Reaction-diffusion-advection systems with discontinuous diffusion and
mass control.
SIAM J. Math. Anal., 53(6):6771–6803, 2021.

[13]
Selwyn L. Hollis, Robert H. Martin, Jr., and Michel Pierre.
Global existence and boundedness in reaction-diffusion systems.
SIAM J. Math. Anal., 18(3):744–761, 1987.

[14]
Said Kouachi.
Existence of global solutions to reaction-diffusion systems via a
Lyapunov functional.
Electron. J. Differential Equations, pages No. 68, 10, 2001.

[15]
El-Haj Laamri.
Global existence of classical solutions for a class of
reaction-diffusion systems.
Acta Applicandae Mathematicae, 115(2):153–165, 2011.

[16]
Philippe Laurençot and Bogdan-Vasile Matioc.
Bounded weak solutions to the thin film muskat problem via an
infinite family of liapunov functionals.
Transactions of the American Mathematical Society,
375(08):5963–5986, 2022.

[17]
El Haj Laamri and Michel Pierre.
Global existence for reaction-diffusion systems with nonlinear
diffusion and control of mass.
Ann. Inst. H. Poincaré C Anal. Non Linéaire,
34(3):571–591, 2017.

[18]
Jeff Morgan.
Global existence for semilinear parabolic systems.
SIAM journal on mathematical analysis, 20(5):1128–1144, 1989.

[19]
Jeff Morgan.
Boundedness and decay results for reaction-diffusion systems.
SIAM Journal on Mathematical Analysis, 21(5):1172–1189, 1990.

[20]
A. Moussa.
Some variants of the classical Aubin-Lions lemma.
J. Evol. Equ., 16(1):65–93, 2016.

[21]
Jeff Morgan and Bao Quoc Tang.
Boundedness for reaction–diffusion systems with lyapunov functions
and intermediate sum conditions.
Nonlinearity, 33(7):3105, 2020.

[22]
Jeff Morgan and Bao Quoc Tang.
Global well-posedness for volume-surface reaction-diffusion systems.
Commun. Contemp. Math., 25(4):Paper No. 2250002, 63, 2023.

[23]
Simon Malham and Jack X. Xin.
Global solutions to a reactive Boussinesq system with front data on
an infinite domain.
Comm. Math. Phys., 193(2):287–316, 1998.

[24]
Michel Pierre.
Global existence in reaction-diffusion systems with control of mass:
a survey.
Milan J. Math., 78(2):417–455, 2010.

[25]
Michel Pierre and Didier Schmitt.
Blowup in reaction-diffusion systems with dissipation of mass.
SIAM Review, 42(1):93–106, 2000.

[26]
Michel Pierre and Didier Schmitt.
Examples of finite time blow up in mass dissipative
reaction-diffusion systems with superquadratic growth.
Discrete Contin. Dyn. Syst., 43(3-4):1686–1701, 2023.

[27]
Vandana Sharma.
Global existence and uniform estimates of solutions to reaction
diffusion systems with mass transport type boundary conditions.
Commun. Pure Appl. Anal., 20(3):955–974, 2021.

[28]
Philippe Souplet.
Global existence for reaction–diffusion systems with dissipation of
mass and quadratic growth.
Journal of Evolution Equations, 18(4):1713–1720, 2018.

[29]
Bao Quoc Tang.
Global classical solutions to reaction-diffusion systems in one and
two dimensions.
Communications in Mathematical Sciences, 16(2):411–423, 2018.

[30]
Alex Viguerie, Guillermo Lorenzo, Ferdinando Auricchio, Davide Baroli, Thomas
J. R. Hughes, Alessia Patton, Alessandro Reali, Thomas E. Yankeelov, and
Alessandro Veneziani.
Simulating the spread of COVID-19 via a spatially-resolved
susceptible-exposed-infected-recovered-deceased (SEIRD) model with
heterogeneous diffusion.
Appl. Math. Lett., 111:Paper No. 106617, 9, 2021.

[31]
Alex Viguerie, Alessandro Veneziani, Guillermo Lorenzo, Davide Baroli, Nicole
Aretz-Nellesen, Alessia Patton, Thomas E. Yankeelov, Alessandro Reali, Thomas
J. R. Hughes, and Ferdinando Auricchio.
Diffusion-reaction compartmental models formulated in a continuum
mechanics framework: application to COVID-19, mathematical analysis, and
numerical study.
Comput. Mech., 66(5):1131–1152, 2020.