\section{Introduction}

\label{sec1}

The amount of theoretical and observational results strongly suggest that our universe had undergone one accelerated phase of expansion during its very early phase of evolution, known as inflation \cite{guth,linde}, and presently it is again undergoing an accelerating phase of expansion, known as late time cosmic acceleration,  dominated either by some hypothetical dark energy fluid adjusted in the Einstein's General Relativity (GR) \cite{copeland} or by some extra geometrical terms appearing due to  modifications of GR or due to new gravitational theories \cite{Clifton1,Nojiri:2017ncd}. So, one accelerating phase (inflation) occurred in the high energy regime and the second acceleration is occurring in the low energy regime and the difference of the energy scales between the early and late accelerating phases of the universe is significantly high. Therefore, understanding the evolution of the universe in both these regimes has been one of the longstanding issues in modern cosmology. Despite of many investigations performed in these sectors, the actual nature is still unrevealed. The unification of the gravitational theory
both at classical and quantum levels has therefore remained as an open problem in modern cosmology.

Concerning the early evolution of the universe, considerable attention on the nature of the gravitational theory at the quantum level has been paid by many investigators. The existing literature demands that the theory of quantum gravity \cite{MB01, Niem01, Kowa01, Niem01b, Kempf01, KN01, EGKS01, Amjad1, Amjad2} played a very crucial role to understand the very early evolution of the universe. However, the approach to reveal the quantum nature of the gravity is not unique although. Various approaches were proposed over the years that include  string theory (ST) \cite{Mukhi:2011zz},
doubly special relativity (DSR) \cite{KowalskiGlikman:2004qa,AmelinoCamelia:2010pd,Ghosh:2006bx,Pramanik:2012fj}, black hole physics (BHP) \cite{Bekenstein1, Bekenstein2, Bekenstein3, Bekenstein4} etc. All of them  reported  the existence of a minimum length scale of the order of the Planck length $l_{\rm pl}$ (equivalently, there exists a maximum energy scale in nature). This actually motivated to generalize  the
Heisenberg's Uncertainty Principle to some Generalized Uncertainty principle (GUP) in quantum gravity  \cite{Maggiore}.

The Generalized Uncertainty principle has been found to be very effective to explain several issues related to the dynamics of the universe. For instance, the origin of the magnetic fields with microgauss
strength \cite{SFW86,K94}  which is observed in the intergalactic regimes of the universe acquires some explanation if GUP is taken into account  \cite{Amjad04}. The potential application of GUP can be found in the context of black hole mechanics. Following  the Bekenstein's
entropy relation \cite{Bekenstein1, Bekenstein2, Bekenstein3, Bekenstein4}
and the Hawking temperature \cite{Hawking74, BD82}, the small black holes in the universe should emit black body radiation. Due to this radiation, these small black holes must loss their mass and consequently becomes hotter.  This emission of radiation  continues until they are completely evaporated. However, the above formulation was based under the assumption of a classical black hole metric together with the consideration that the emitting radiation
of the black hole was much much smaller (can be ignored) compared to the rest mass energy of the black hole. Now, during this emission process, the black holes naturally become smaller and lighter, and certainly within this length scale,  the assumption of classical metric and the ignorance of the
radiation seem to be invalid. Thus,  whether the black holes will completely evaporate to either photons, ordinary matter particles or vacuum, or whatever it be (some remnant for instance) is highly questionable due to the lack of a definite quantum gravitational theory (see the discussions \cite{Bunch81, CMP88, York, Parikh99, Susskind95}). Here,  GUP finds an answer towards this discrepancy \cite{Adler}. As argued in Ref. \cite{Adler} a small black
hole having temperature greater than the ambient temperature will keep  radiating  photons and ordinary particles until it reaches the Planck scale. When the black hole will reach the Planck size, it will stop radiating and its entropy reaches zero while  its effective temperature will reach the maximum limit, and finally, it reduces to an inert
remnant having no radiating power, but only gravitational interactions may exist. It is not necessary that the created remnant will enjoy the horizon structure as in classical black hole \cite{York}. Such remnants might be consider as the dark matter candidate \cite{Cline}. This is a very surprising result because GUP indicates the origin of one of the heavy fluids of the universe. Naturally, the effects of GUP in various cosmological theories can be studied in order to see whether such effects give rise to some interesting directions related the dynamical history of the universe.
Since  dark energy is another heavy resource of the universe playing the leading role behind its accelerated expansion, so how GUP modified dark energy models behave, will be an interesting direction of research. As the minimum length effects are quite universal at any stage, so, there is no reason to exclude the GUP modified models, rather, one can generalize such cosmic theories.
Thus, in this work shall focus on the cosmological models that are modified by the GUP. The work has been organized in the following way.

In Section \ref{sec2} provide with a brief description on the GUP and the corresponding algebra. In Section \ref{sec3} we discuss the quintessence scalar field model modified by the GUP. Then in Section \ref{sec4} we discuss the gravitational field equations of the modified quintessence scalar field in the Friedmann-Lema\^{i}tre-Robertson-Walker universe.
In Section \ref{sec5} we define the dimensionless dynamical variables and perform the stability analysis. After that in Section \ref{sec6} we introduce the kinematical quantities and compare with the present model. Finally, in Section \ref{sec7} we close the work with a brief summary of all the findings.

\section{Generalized Uncertainty Principle}

\label{sec2}

The GUP is related with the existence of a minimum measurable length, such
that the Heisenberg uncertainty is to be modified as
\begin{equation}
\Delta X_{i}\Delta P_{j}\geqslant\frac{\hbar}{2}[\delta_{ij}(1+\beta
P^{2})+2\beta P_{i}P_{j}]. \label{GUP}%
\end{equation}
where $\beta$ is the deformation parameter defined as $\beta={\beta}_{0}%
\ell_{Pl}^{2}/2\hbar^{2}$, which is also known as quadratic generalized GUP.
From (\ref{GUP}) the generalized deformation of the Heisenberg algebra in a
four-dimensional Minkowski spacetime of signature
\cite{Quesne2006,Vagenas,Kemph1,Kemph2} follows%
\begin{equation}
\lbrack X_{\mu},P_{\nu}]=-i\hbar\lbrack(1-\beta(\eta^{\alpha\gamma}P_{\alpha
}P_{\gamma}))\eta_{\mu\nu}-2\beta P_{\mu}P_{\nu}]. \label{con11}
\end{equation}
From the algebra (\ref{con11}) it leads to the deformation of the coordinate
representation of the momentum operator \cite{ss1}.

By using the latter commutation relation we can write the coordinate
representation of the operators $X_{\mu},P_{\nu}$. We select $X_{\mu}$
undeformed, that is $X_{\mu}=x_{\mu}$, and the momentum operator to modified
as \cite{Moayedi}
\begin{equation}
P_{\mu}=p_{\mu}(1-\beta(\eta^{\alpha\gamma}p_{\alpha}p_{\gamma}))~,
\label{xp-form}%
\end{equation}
in which~$p^{\mu}$ are defined as usual, that is, $p^{\mu}=i\hbar
\frac{\partial}{\partial x_{\mu}},$ and $[x_{\mu},p_{\nu}]=-i\hbar\eta_{\mu
\nu}$.~

In the context of the GUP, the Klein-Gordon equation which describes a spin-0
particle, described by the wave function $\Psi$, it is now a fourth-order
partial differential equation as given
\begin{equation}
\Delta\Psi+2\beta\hbar\Delta\left(  \Delta\Psi\right)  +\left(  \frac
{mc}{\hbar}\right)  ^{2}\Psi+O\left(  \beta^{2}\right)  =0, \label{kg.01}%
\end{equation}
where $\Delta$ is the Laplace operator defined as $\Delta=\square$ when the
underlying space is $\eta_{\mu\nu}$, or in general $\Delta=\frac{1}{\sqrt{-g}%
}\partial_{\mu}\left(  g^{\mu\nu}\sqrt{-g}\partial_{\mu}\right)  $ for any
Riemannian space of Lorentzian signature with metric $g_{\mu\nu}$. Equation
(\ref{kg.01}) is a singular perturbative equation.

An important observation that helps us to generalize the quintessence scalar
field model in GUP, is that the fourth-order differential equation
(\ref{kg.01}) follows from the variation of the action integral%
\begin{equation}
S_{SF}=\int dx^{4}\sqrt{-g}L\left(  \Psi,\mathcal{D}_{\sigma}\Psi\right),
\label{kg.02}%
\end{equation}
where $L\left(  \Psi,\mathcal{D}_{\sigma}\Psi\right)  $ is the usual
Lagrangian for the Klein-Gordon equation%
\begin{equation}
L\left(  \Psi,\mathcal{D}_{\sigma}\Psi\right)  =\frac{1}{2}g^{\mu\nu
}\mathcal{D}_{\mu}\Psi\mathcal{D}_{\nu}\Psi-\frac{1}{2}\left(  \frac{mc} %
{\hbar}\right)  ^{2}\Psi^{2}, \label{kg.03} %
\end{equation}
where the new operator $\mathcal{D}_{\mu}$ is defined as%
\begin{equation}
\mathcal{D}_{\mu}=\nabla_{\mu}+\beta\hbar^{2}\nabla_{\mu}\left(
\Delta\right)  , \label{kg.04} %
\end{equation}
where $\nabla_{\mu}$ is the covariant derivative for the metric tensor
$g_{\mu\nu}$.

We introduce the new field $\Phi=\Delta\Psi$, and the Lagrange multiplier
$\lambda.$ Variation with respect to the Lagrange multiplier gives,
$\frac{\delta S}{\delta\lambda}=0$ from where it follows $\lambda=-2\beta
\hbar^{2}\Phi$, such that the action integral (\ref{kg.02}) to be written as%
\begin{equation}
S=\int dx^{4}\sqrt{-g}\left(  \frac{1}{2}g^{\mu\nu}\nabla_{\mu}\Psi\nabla
_{\nu}\Psi+2\beta\hbar^{2}g^{\mu\nu}\nabla_{\mu}\Psi\nabla_{\nu}\Phi
+\beta\hbar^{2}\Phi^{2}-\frac{1}{2}V_{0}\Psi^{2}\right)  . \label{pd.66}%
\end{equation}

Therefore, the new Lagrangian is%
\begin{equation}
L\left(  \Psi,\Psi_{;\mu},\Phi,\Phi_{;\mu}\right)  =\sqrt{-g}\left(  \frac
{1}{2}g^{\mu\nu}\Psi_{;\mu}\Psi_{;\nu}+2\beta\hbar^{2}g^{\mu\nu}\Psi_{;\mu
}\Phi_{;\nu}\right)  -\sqrt{-g}\left(  \frac{1}{2}V_{0}\Psi^{2}-\beta\hbar
^{2}\Phi^{2}\right)  \label{pd.67} %
\end{equation}
Hence, the Euler-Lagrange equations
for Lagrangian (\ref{pd.67})
are%
\begin{equation}
g^{\mu\nu}\Psi_{,\mu\nu}-\Gamma^{\mu}\Psi_{,\mu}-\Phi=0 \label{pd.69}%
\end{equation}%
\begin{equation}
2\beta\hbar^{2}\left(  g^{\mu\nu}\Phi_{,\mu\nu}-\Gamma^{\mu}\Phi_{,\mu
}\right)  +\left(  V_{0}\Psi+\Phi\right)  =0 \label{pd.70}%
\end{equation}
where equation (\ref{pd.69}) is the constraint $\Phi=\Delta_{g}\Psi$, and
system (\ref{pd.69}), (\ref{pd.70}) is equivalent with the fourth-order
partial differential equation (\ref{kg.01}).

\section{Quintessence modified by the GUP}

\label{sec3}

Quintessence is one of the simple dark energy models. Specifically, a canonical
scalar field $\phi\left(  x^{\mu}\right)  $ is responsible for the
acceleration era of the universe. The late time acceleration or the early
acceleration phase, where the scalar field plays the role of the inflaton
\cite{ratra1}.

In the case of a four-dimensional Riemannian manifold described by the metric
tensor $g_{\mu\nu}$, with Ricci scalar $R$, the Action integral which describes
the gravitational theory for the quintessence model is%
\begin{equation}
S_{Q}=\int dx^{4}\sqrt{-g}\left(  -\frac{R}{2}+\frac{1}{2}g^{\mu\nu}%
\nabla_{\mu}\phi\nabla_{\nu}\phi-V(\phi)\right)  . \label{kg.05}%
\end{equation}
However, by using the operator $\mathcal{D}_{\mu}$ defined in (\ref{kg.04}),
we can generalize the latter Action integral as follows%
\begin{equation}
S_{Q}^{GUP}=\int dx^{4}\sqrt{-g}\left(  -\frac{R}{2}+\frac{1}{2}g^{\mu\nu
}\mathcal{D}_{\mu}\phi\mathcal{D}_{\nu}\phi-V(\phi)\right)  , \label{kg.06}%
\end{equation}
which reduces to (\ref{kg.05}) when $\beta=0$.

However, we can always define a new field $\psi=\Delta\phi$, as before and
write the Action Integral (\ref{kg.06}) like that of a second-order theory,
that is,%
\begin{equation}
S_{Q}^{GUP}=\int dx^{4}\sqrt{-g}\left(  -\frac{R}{2}+\frac{1}{2}g^{\mu\nu
}\nabla_{\mu}\phi\nabla_{\nu}\phi-V(\phi)+\beta\hbar^{2}\left(  2g^{\mu\nu
}\nabla_{\mu}\phi\nabla_{\nu}\psi+\psi^{2}\right)  \right)  . \label{kg.09}%
\end{equation}

Variation with respect to the metric tensor of (\ref{kg.09}) leads to the
gravitational field equations, which are,
\begin{equation}
G_{\mu\nu}=T_{\mu\nu}, \label{kg.10}%
\end{equation}
where $G_{\mu\nu}$ is the Einstein tensor and~$T_{\mu\nu}$ is the energy
momentum tensor for the two fields defined as
\begin{align}
T_{\mu\nu}  &  =\left(  \frac{1}{2}\nabla_{\mu}\phi\nabla_{\nu}\phi
+2\beta\hbar^{2}\nabla_{\mu}\phi\nabla_{\nu}\psi\right)  +\nonumber\\
&  -g_{\mu\nu}\left(  \frac{1}{2}g^{\alpha\beta}\nabla_{\alpha}\phi
\nabla_{\beta}\phi+2\beta\hbar^{2}g^{\alpha\beta}\nabla_{\alpha}\phi
\nabla_{\beta}\psi-V(\phi)+\beta\hbar^{2}\psi^{2}\right)  , \label{kg.11}%
\end{align}

Furthermore, variation of (\ref{kg.09}) with respect to the scalar fields $\phi,~\psi$ leads to
the equation of motions for the two fields%
\begin{equation}
\Delta\phi-\psi=0, \label{kg.12}%
\end{equation}%
\begin{equation}
\beta\hbar^{2}\Delta\psi+\frac{1}{2}\left(  \psi+V_{,\phi}\right)  =0.
\label{kg.13}%
\end{equation}

\section{FLRW background space}

\label{sec4}

In the large scales our universe is almost homogeneous and isotropic and
such description of the universe is well described by the Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) line element with zero spatial curvature given by,
\begin{equation}
ds^{2}=dt^{2}-a^{2}\left(  t\right)  \left(  dx^{2}+dy^{2}+dz^{2}\right)  ,
\label{frw.01}%
\end{equation}
where $a\left(  t\right)  $ is the expansion scale factor of the universe, and $H=\dot{a}/a$ is
the Hubble function.

For the line element (\ref{frw.01}) the gravitational field equations
(\ref{kg.10}) are%
\begin{equation}
3H^{2}=\left(  \frac{1}{2}\dot{\phi}^{2}+V\left(  \phi\right)  \right)
+\beta\hbar^{2}\left(  2\dot{\phi}\dot{\psi}-\psi^{2}\right)  +\rho_{m}
\label{frw.02}%
\end{equation}%
\begin{equation}
2\dot{H}+3H^{2}=-\left(  \frac{1}{2}\dot{\phi}^{2}-V\left(  \phi\right)
+\beta\hbar^{2}\left(  2\dot{\phi}\dot{\psi}+\psi^{2}\right)  \right)  ,
\label{frw.03}%
\end{equation}
where $\rho_{m}$ is the energy density for the dust fluid source, i.e.
$T_{\mu\nu}^{(m)}=\rho_{m}u_{\mu}u_{\nu}$, in which $u^{\mu}=\delta_{t}^{\mu}$.

For the fields $\phi,~\psi$ the equation of motions are%
\begin{align}
\ddot{\phi}+\frac{3}{a}\dot{a}\dot{\phi}-\psi &  =0,\label{frw.04}\\
\varepsilon\left(  \ddot{\psi}+\frac{3}{a}\dot{a}\dot{\psi}\right)  +\frac
{1}{2}\left(  \psi+V_{,\phi}\right)   &  =0, \label{frw.05}%
\end{align}
while for the dust fluid source, the Bianchi identity $T_{~\ ~~~~;\nu
}^{(m)~\mu\nu}=0$ gives $\dot{\rho}_{m}+3H\rho_{m}=0,$ that is $\rho_{m}%
=\rho_{m0}a^{-3}$, where $\rho_{m0}$ is the present value of $\rho_m$.

An equivalent way to write the field equations is to define the new variables%
\begin{align}
&  3H^{2}=\rho_{GUP}+\rho_{m},\\
&  2\dot{H}+3H^{2}=-p_{GUP},
\end{align}
where
\[
\rho_{GUP}=\rho_{\phi}+\rho_{\psi},~p_{GUP}=p_{\phi}+p_{\psi}.
\]

Parameters $\rho_{\phi}$, $~p_{\phi}$ are the energy density and pressure for
the quintessence
\begin{equation}
\rho_{\phi}=\frac{1}{2}\dot{\phi}^{2}+V\left(  \phi\right)  ~,~p_{\phi}%
=\frac{1}{2}\dot{\phi}^{2}-V\left(  \phi\right)  ,
\end{equation}
and $\rho_{\psi}$, $p_{\psi}$ are the energy and pressure for the second
interacting field $\psi$, that is,
\begin{equation}
\rho_{\psi}=\beta\hbar^{2}\left(  2\dot{\phi}\dot{\psi}-\psi^{2}\right)
~,~p_{\psi}=\beta\hbar^{2}\left(  2\dot{\phi}\dot{\psi}+\psi^{2}\right)  .
\end{equation}
At this point we remark that the second field $\psi$ is not a physical field,
but it is introduced by the quantum corrections of GUP and it describes the
new geometrodynamic degrees of freedom of GUP.

As far as concerns the parameter for the equation of state of the scalar
field, that is calculated
\begin{equation}
w_{GUP}=\frac{p_{GUP}}{\rho_{GUP}}=\frac{\left(  \frac{1}{2}\dot{\phi}%
^{2}-V\left(  \phi\right)  \right)  +\beta\hbar^{2}\left(  2\dot{\phi}%
\dot{\psi}+\psi^{2}\right)  }{\left(  \frac{1}{2}\dot{\phi}^{2}+V\left(
\phi\right)  \right)  +\beta\hbar^{2}\left(  2\dot{\phi}\dot{\psi}-\psi
^{2}\right)  }. \label{ac.20}%
\end{equation}
where expanding around~$\beta\hbar^{2}\rightarrow0$, it follows
\begin{equation}
w_{GUP}\simeq w_{\phi}+\frac{1}{\rho_{\phi}}\left(  p_{\psi}-w_{\phi}%
\rho_{\psi}\right)  +O\left(  \left(  \beta\hbar^{2}\right)  ^{2}\right)
\label{ac.20c}%
\end{equation}
while when the scalar field potential $V\left(  \phi\right)  $, dominates, that
is $\dot{\phi}^{2}\ll V$ it follows%
\begin{equation}
w_{GUP}\simeq-1+\frac{4\beta\hbar^{2}}{V}\dot{\phi}\dot{\psi}%
\end{equation}
which can cross the phantom divide line when $~\dot{\phi}\dot{\psi}\leq0$. An
important observation here is that the field equations form a singular
pertubation system, also known as slow-fast dynamical system. That means that
it is possible$~\left\vert \dot{\psi}\right\vert $ to be large enough such
that the term $\left\vert \beta\hbar^{2}\dot{\psi}\right\vert \ $\ will not be
negligible. That is a main characteristic of the slow-fast dynamical systems.

In the following sections we study the global evolution of the dynamics, in
order to understand the effects of GUP in the quintessence field. This work
extends and generalize the previous work on the specific theory \cite{Paliathanasis:2015cza}.

\section{Dynamical systems formulation}

\label{sec5}

We continue our analysis by define the new dimensionless variables%
\begin{equation}
x_{1}=\frac{\dot{\phi}}{\sqrt{6}H}~,~y_{1}=\sqrt{\frac{V}{3H^{2}}}%
~,~x_{2}=\beta\hbar^{2}\frac{2\sqrt{2}\dot{\psi}}{\sqrt{3}H}~,~y_{2}%
=\frac{\beta\hbar^{2}\psi^{2}}{3H^{2}},~\Omega_{m}=\frac{\rho_{m}}{3H^{2}%
},\label{gp.01}%
\end{equation}
such that the gravitational field equations to be written as%
\begin{equation}
\frac{dx_{1}}{d\tau}=\frac{1}{4}\left(  6x_{1}\left(  x_{1}\left(  x_{1}%
+x_{2}\right)  -y_{1}^{2}-1\right)  +y_{2}\left(  6x_{1}-\sqrt{6}\mu\right)
\right)  ,\label{gp.02}%
\end{equation}%
\begin{equation}
\frac{dy_{1}}{d\tau}=\frac{1}{2}y_{1}\left(  3\left(  1-y_{1}^{2}%
+y_{2}\right)  +x_{1}\left(  3\left(  x_{1}+x_{2}\right)  -\sqrt{6}%
\lambda\right)  \right)  ,\label{gp.03}%
\end{equation}%
\begin{equation}
\frac{dx_{2}}{d\tau}=\frac{1}{2}\left(  3x_{1}x_{2}^{2}-3x_{2}\left(
1-x_{1}^{2}+y_{1}^{2}-y_{2}\right)  +\sqrt{6}\left(  2\lambda y_{1}^{2}+\mu
y_{2}\right)  \right)  \label{gp.04}%
\end{equation}%
\begin{equation}
\frac{dy_{2}}{d\tau}=\frac{1}{4}y_{2}\left(  12x_{1}\left(  x_{1}%
+x_{2}\right)  +12\left(  1-y_{1}^{2}-y_{2}\right)  -\sqrt{6}x_{2}\mu\right)
\label{gp.05}%
\end{equation}%
\begin{equation}
\frac{d\mu}{d\tau}=\frac{1}{4}\sqrt{\frac{3}{2}}x_{2}\mu^{2}~,~\frac{d\lambda
}{d\tau}=\sqrt{6}x_{1}\lambda^{2}\left(  \Gamma\left(  \lambda\right)
-1\right)  \label{gp.06}%
\end{equation}
where the new independent variable $\tau=\ln a$, and variables $\lambda,\mu$
are defined as%
\begin{equation}
\lambda=-\frac{V_{,\phi}}{V}~,~\beta\hbar^{2}\mu=-\frac{2}{\psi}\label{gp.07}%
\end{equation}
while $\Gamma\left(  \lambda\right)  =\frac{V_{,\phi\phi}V}{V_{,\phi}^{2}}$.
Moreover, the constraint equation is written as%
\begin{equation}
\Omega_{m}\left(  \mathbf{x,y}\right)  =1-x_{1}\left(  x_{1}+x_{2}\right)
-y_{1}^{2}+y_{2},\label{gp.08}%
\end{equation}
from where it follows that the paramters are constraint as $\,0\leq\Omega
_{m}\leq1.$

The parameter for the equation of state for the effective fluid is derived to
be%
\begin{equation}
w_{tot}\left(  \mathbf{x,y}\right)  =x_{1}\left(  x_{1}+x_{2}\right)
-y_{1}^{2}+y_{2}, \label{gp.09}%
\end{equation}
while we define the physical variables $\Omega_{\phi},~\Omega_{\psi}$ such as,
$\Omega_{\phi}\left(  \mathbf{x,y}\right)  =x_{1}^{2}+y_{1}^{2}~,~\Omega
_{\psi}\left(  \mathbf{x,y}\right)  =x_{1}x_{2}-y_{2},~$where the constraint
equation is written as $\Omega_{m}\left(  \mathbf{x,y}\right)  =1-\Omega
_{\phi}\left(  \mathbf{x,y}\right)  -\Omega_{\psi}\left(  \mathbf{x,y}\right)
$.

Finally, the parameter for the equation of state for the scalar field is
expressed in terms of the new variables as
\begin{equation}
w_{\phi}\left(  \mathbf{x,y}\right)  =\frac{\left(  x_{1}^{2}-y_{1}%
^{2}\right)  +x_{1}x_{2}+y_{2}}{\left(  x_{1}^{2}+y_{1}^{2}\right)
+x_{1}x_{2}-y_{2}}. \label{gp.10}%
\end{equation}

\subsection{Equilibrium points for exponential potential}

Consider now that the scalar field potential is exponential, that is $V\left(
\phi\right)  =V_{0}e^{-\lambda\phi}$, such that $\Gamma\left(  \lambda\right)
=1$ and $\frac{d\lambda}{d\tau}=0$. For that potential the dynamical system
(\ref{gp.02})-(\ref{gp.06}) is reduced by one-dimension. Therefore, the
equilibrium points of the dynamical system are the points with coordinates
$P=\left(  x_{1}\left(  P\right), x_{2}\left(  P\right), y_{1}\left(P\right), y_{2}\left(  P\right), \mu\left(  P\right)  \right) $, where the
right hand side of (\ref{gp.02})-(\ref{gp.06}) are zero.

The equilibrium points are calculated to be%
\begin{equation}
P_{0}=\left(  0,0,0,0,0,\mu\right)  ~,~P_{1}^{\pm}=\left(  \pm1,0,0,0,\mu
\right)  ~,
\end{equation}%
\begin{equation}
~P_{2}=\left(  x_{1},\frac{1}{x_{1}}-x_{1},0,0,0\right)  ~,~P_{3}=\left(
0,\sqrt{\frac{2}{3}}\lambda y_{1}^{2},y_{1},y_{1}^{2}-1,0\right)  .
\end{equation}
In the following lines we discuss the physical properties of the points as
also their stability.

Point $P_{0}\,\ $describes a universe dominated by the dust fluid source,
where $\Omega_{m}\left(  P_{0}\right)  =1$, while the parameter for the
equation of state $w_{tot}$ is derived to be $w_{tot}\left(  P_{0}\right)
=0$. The eigenvalues of the linearized system are calculated to be
$e_{1}\left(  P_{0}\right)  =3~,~e_{2}\left(  P_{0}\right)  =-\frac{3}%
{2}~,~e_{3}\left(  P_{0}\right)  =-\frac{3}{2},~e_{4}\left(  P_{0}\right)
=\frac{3}{2}$ and $e_{5}\left(  P_{0}\right)  =0$, from where we infer that
the point is a saddle point and the exact solution at the point is always unstable.

Poins $P_{1}^{\pm}$ describe universes dominated by the kinetic term of the
scalar field, i.e. $\Omega_{m}\left(  P_{1}^{\pm}\right)  =0$, and
$w_{tot}\left(  P_{1}^{\pm}\right)  =1.$ The eigenvalues are calculated to be
$e_{1}\left(  P_{1}^{\pm}\right)  =6~,~e_{2}\left(  P_{1}^{\pm}\right)
=3~,~e_{3}\left(  P_{1}^{\pm}\right)  =\frac{1}{2}\left(  6\mp\sqrt{6}%
\lambda\right)  ,~e_{4}\left(  P_{1}^{\pm}\right)  =0$ and $e_{5}\left(
P_{1}^{\pm}\right)  =0$, from where we infer that the exact solutions at the
equilibrium points are always unstable. Point $P_{1}^{+}$ is a saddle point
when $\lambda>\sqrt{6}$ while point $P_{1}^{-}$ is a saddle point when
$\lambda<-\sqrt{6}$.

Point $P_{2}$ is actually family of points that lie on a two-dimensional surface.
The physical properties of this points are similar with that of points
$P_{1}^{\pm}$. The eigenvalues of the linearized system are derived
$e_{1}\left(  P_{2}\right)  =6~,~e_{2}\left(  P_{2}\right)  =3~,~e_{3}\left(
P_{2}\right)  =\frac{1}{2}\left(  6-\sqrt{6}x_{1}\lambda\right)
,~e_{4}\left(  P_{2}\right)  =0$ and $e_{5}\left(  P_{2}\right)  =0$, from
where we infer that the exact solutions at the family of points are always unstable.

The family of points on the surface with coordinates $P_{3}$ describe de
Sitter universes where $\Omega_{m}\left(  P_{3}\right)  =0$ and $w_{tot}%
\left(  P_{3}\right)  =-1$. The eigevalues of the linearized system are
determined to be $e_{1}\left(  P_{3}\right)  =-3~,~e_{2}\left(  P_{3}\right)
=-3~,~e_{3}\left(  P_{3}\right)  =-3,~e_{4}\left(  P_{3}\right)  =0$ and
$e_{5}\left(  P_{2}\right)  =0$, hence the Center Manifold Theorem (CMT) should be applied to find the
manifold where the de Sitter universes are attractors.

We observe that no scaling solutions or a tracking solutions exist in this
specific model like in the quintessence theory. However, the critical points
which describe the de Sitter solution do not exist in the case of quintessence
for the exponential potential; these are the new equilibrium points provided by
the new terms given by GUP.

We solve numerically the field equations (\ref{gp.02})-(\ref{gp.06}) for three
different sets of initial conditions, and we present the evolution of the
physical parameters $\Omega_{m},~\Omega_{\phi}=1-\Omega_{m}$ and
$w_{tot},~w_{\phi}$ if Fig. \ref{fig0}. For the three different sets of the
initial conditons the final attractor is the de Sitter universe, while for the
parameters of the equation of state we observe that they can cross the phantom
divide line. In order to understand that behaviour we present the phase space
diagram for the dynamical system in the plane $\left\{  x_{1}-y_{1}\right\}  $
from where it is clear that $P_{3}$ can be a local attractor.

To analyze $P_{3}$ in the following we proceed with
the application of the CMT.

\begin{figure}[ptb]
\centering\includegraphics[width=0.8\textwidth]{gup1.pdf} \caption{Qualitative
evolution of the physical variables $\Omega_{m},~\Omega_{\phi}$ (Left
figs.)\ and $w_{tot},~w_{\phi}$ (Right figs) for three sets of initial
conditions. The plots of the first row is for initial conditions near to the
point $P_{0}$, the second row is for initial condition near to the point
$P_{1}^{+}$, while the plots of the third row are for initial conditions near
to the point $P_{1}^{-}$. For the energy density we observe that at late times
the scalar field dominates $\Omega_{\phi}\rightarrow1$ and $\Omega
_{m}\rightarrow0$, while the parameter for the equation of state $w_{tot}$
have the limits $w_{tot}\rightarrow-1$. Thus, we observe that $w_{\phi}$ and
$w_{tot}$ can cross the phantom divide line and take values smaller than minus
one. }%
\label{fig0}%
\end{figure}

\begin{figure}[ptb]
\centering\includegraphics[width=0.5\textwidth]{gup2.pdf} \caption{Phase space
diagram in the plane $\left\{  x_{1}-y_{1}\right\}  $ where for the rest of
the parameters we selected $\lambda=-1$, $y_{2}=0.01$, $x_{2}=0.01$ and
$\mu=0.1$.~From the diagram we observe the attractor $P_{3}.$}%
\label{fig1}%
\end{figure}

\subsubsection{Center manifold theorem  for line of  points $P_{3}$}

Assuming $y_{1c}\notin\left\{  0, \frac{\sqrt{2}}{2}, 1\right\}, \lambda \neq 0$, and introducing the new variables
\begin{align}
  u_1 &= \frac{1}{3} {y_{1c}} \left(\left({y_{1c}}^2-1\right) \left({y_{1c}} \left(\sqrt{6} \lambda  {x_1}+\lambda  \mu  \left({y_{1c}}^2-1\right)+3\right)-6 {y_1}\right)+3 {y_{1c}}
   {y_2}\right),\\
  u_2 &= -\frac{1}{2} \lambda  \mu  {y_{1c}}^2 \left(2 {y_{1c}}^4-3 {y_{1c}}^2+1\right),\\
   &v_1= \frac{1}{18} \Big\{6 \lambda  {y_{1c}}^2 \left(-2 \lambda  {x_1}
   \left({y_{1c}}^2-1\right)-\sqrt{6} \left({y_{1c}}^2+{y_2}\right)\right)+18 {x_2}+12 \sqrt{6} \lambda  {y_1} \left({y_{1c}}^2-1\right) {y_{1c}} \nonumber \\
   & -\sqrt{6} \mu  \left({y_{1c}}^2-1\right) \left(\lambda ^2
   {y_{1c}}^2 \left(3 {y_{1c}}^2-2\right)+3\right)\Big\},\\
  v_2 &= \frac{1}{6} \lambda  {y_{1c}}^2 \left(-\lambda  \left(6 {x_1} \left({y_{1c}}^2-2\right)+\sqrt{6} \mu  \left({y_{1c}}^2-1\right)^2\right)-3
   \sqrt{6} \left(-2 {y_1} {y_{1c}}+{y_{1c}}^2+{y_2}+1\right)\right),\\
   v_3 & = -\frac{1}{3} \left({y_{1c}}^2-1\right) \left({y_{1c}}^2 \left(\sqrt{6} \lambda  {x_1}+\lambda  \mu
   \left({y_{1c}}^2-1\right)+3\right)-6 {y_1} {y_{1c}}+3 ({y_2}+1)\right),
\end{align}
a particular  point $P_{3c}=\left(
0,\sqrt{\frac{2}{3}}\lambda y_{1c}^{2},y_{1c},y_{1c}^{2}-1,0\right)$ on line of points $P_3$ is translated to the origin and the linearization of the system \eqref{gp.01}, \eqref{gp.02}, \eqref{gp.03}, \eqref{gp.04}, \eqref{gp.05} is transformed to its real Jordan canonical form:
\begin{equation}
 \left(
\begin{array}{c}
 v_1'\\
 v_2'\\
v_3'\\
u_1'\\
u_2'\\
\end{array}
\right)   =  \left(
\begin{array}{ccccc}
 -3 & 0 & 0 & 0 & 0 \\
 0 & -3 & 1 & 0 & 0 \\
 0 & 0 & -3 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 \\
 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)  \left(
\begin{array}{c}
 v_1\\
 v_2\\
v_3\\
u_1\\
u_2\\
\end{array}
\right).
\end{equation}

Assuming $y_{1c}\notin\left\{ 0, \frac{1}{\sqrt{2}}, 1\right\}, \lambda \neq 0$, and applying the CMT we obtain that the center manifold is given locally (up to second order) by the graph
\begin{align}
   \Bigg\{(u_1,u_2, v_1, v_2, v_3)\in\mathbb{R}^5: \quad   & v_1 = a_1 u_1^2 + a_2 u_1 u_2 + a_3 u_2^2 + \mathcal{O}(3), \nonumber \\
    & v_2 = b_1 u_1^2 + b_2 u_1 u_2 + b_3 u_2^2 + \mathcal{O}(3), \nonumber \\
    & v_3 = c_1 u_1^2 + c_2 u_1 u_2 + c_3 u_2^2 + \mathcal{O}(3)\Bigg\},
\end{align}
where
\begin{align}
& a_1= -\frac{\lambda  \left({y_{1c}}^2-1\right)}{2 \sqrt{6} {y_{1c}}^2}, a_2= \frac{\lambda ^2 \left(-11 {y_{1c}}^4+10 {y_{1c}}^2-1\right)-6}{3 \sqrt{6} \lambda  {y_{1c}}^2 \left(2 {y_{1c}}^4-3
   {y_{1c}}^2+1\right)}, a_3=\frac{6 \left(4 {y_{1c}}^2-3\right)+\lambda ^2 \left(40 {y_{1c}}^6-57 {y_{1c}}^4+20 {y_{1c}}^2-1\right)}{18 \sqrt{6} \lambda  \left({y_{1c}}^2-1\right) \left({y_{1c}}-2
   {y_{1c}}^3\right)^2},\\
   & b_1=-\frac{1}{4} \sqrt{\frac{3}{2}} \lambda , b_2=\frac{5 \lambda }{\sqrt{6} \left(2-4 {y_{1c}}^2\right)}, b_3= \frac{\lambda ^2 {y_{1c}}^2 \left(12 {y_{1c}}^2-5\right)+6}{12 \sqrt{6} \lambda  \left({y_{1c}}-2
   {y_{1c}}^3\right)^2},\\
   & c_1= \frac{1}{4} \left(\frac{1}{{y_{1c}}^2}-1\right), c_2=\frac{1}{2-4 {y_{1c}}^2}-\frac{1}{6 {y_{1c}}^2}, c_3= \frac{{y_{1c}}^2 \left(\lambda ^2 \left(12 {y_{1c}}^4-9 {y_{1c}}^2+1\right)+6\right)-6}{36
   {y_{1c}}^4 \left(\lambda -2 \lambda  {y_{1c}}^2\right)^2}.
\end{align}
Hence, the dynamics on the center manifold is given locally (up to second order) by
\begin{align}
   &u_1'=\frac{{u_1} {u_2} \left(2 {y_{1c}}^6+7 {y_{1c}}^4-8 {y_{1c}}^2+1\right)}{4 {y_{1c}}^6-6 {y_{1c}}^4+2 {y_{1c}}^2}+\frac{{u_2}^2 \left(\lambda ^2 \left(2 {y_{1c}}^4-3
   {y_{1c}}^2+1\right)-6\right)}{6 \lambda ^2 {y_{1c}}^2 \left(1-2 {y_{1c}}^2\right)^2}+{u_2}  + \mathcal{O}(3),  \label{center1a}\\
   & u_2'=-\frac{{u_1} {u_2}^2}{4 {y_{1c}}^6-6 {y_{1c}}^4+2 {y_{1c}}^2}+\frac{{u_2}^3 \left(\lambda ^2
   {y_{1c}}^4+3\right)}{6 {y_{1c}}^4 \left({y_{1c}}^2-1\right) \left(\lambda -2 \lambda  {y_{1c}}^2\right)^2}-\frac{{u_2}^2}{4 {y_{1c}}^4-6 {y_{1c}}^2+2} + \mathcal{O}(3) \label{center1b}.
\end{align}

\begin{figure}[ptb]
\centering\includegraphics[width=0.8\textwidth]{center1.pdf} \caption{Dynamics on the center manifold is given locally (up to second order) by \eqref{center1a}, \eqref{center1b} for several choices of the parameters. }%
\label{fig3}%
\end{figure}
Defining the function
\begin{equation}
    U(u_1, u_2)=\frac{{u_1}^2}{2}+{u_2}^2,
\end{equation}
we have
\begin{align}
   & U'(\tau)=\nabla U \cdot (u_1', u_2')=u_1 u_2 + \frac{{u_1}^2 {u_2} \left(2 {y_{1c}}^6+7 {y_{1c}}^4-8 {y_{1c}}^2+1\right)}{4 {y_{1c}}^6-6 {y_{1c}}^4+2 {y_{1c}}^2} \nonumber \\
    &+{u_1} \left(\frac{{u_2}^2 \left(\lambda ^2 \left(2 {y_{1c}}^4-3
   {y_{1c}}^2+1\right)-6\right)}{6 \lambda ^2 {y_{1c}}^2 \left(1-2 {y_{1c}}^2\right)^2}-\frac{{u_2}^3}{2 {y_{1c}}^6-3 {y_{1c}}^4+{y_{1c}}^2}\right)\nonumber \\
   & +\frac{{u_2}^4 \left(2 \lambda ^2 {y_{1c}}^4+6\right)}{6
   {y_{1c}}^4 \left({y_{1c}}^2-1\right) \left(\lambda -2 \lambda  {y_{1c}}^2\right)^2}+\frac{{u_2}^3}{-2 {y_{1c}}^4+3 {y_{1c}}^2-1}\sim u_1 u_2 + \mathcal{O}(3).
\end{align}
Hence, for small $u_1, u_2$ such as $u_1 u_2 < 0$ the origin is locally stable, whereas,  for small $u_1, u_2$ such as $u_1 u_2 >0$ it is unstable.
\begin{figure}[ptb]
\centering\includegraphics[width=0.8\textwidth]{center2.pdf} \caption{$U'(\tau)$ for several choices of the parameters, showing that for small $u_1, u_2$ such as $u_1 u_2 < 0$ the origin is locally stable, whereas,  for small $u_1, u_2$ such as $u_1 u_2 >0$ it is unstable.}%
\label{fig4}%
\end{figure}

In the figure \ref{fig3} is it presented the dynamics on the center manifold is given locally (up to second order) by \eqref{center1a}, \eqref{center1b} for several choices of the parameters.
In figure \ref{fig4} it is presented  $U'(\tau)$ for the same choices of the parameters, showing that for small $u_1, u_2$ such as $u_1 u_2 < 0$ the origin is locally stable, whereas,  for small $u_1, u_2$ such as $u_1 u_2 >0$ it is unstable.

\textbf{Special case $y_{1c}=\frac{\sqrt{2}}{2}$}. Using the transformation
\begin{align}
 & {u_1}= \mu , \\
 & {u_2}= \frac{1}{24} \left(\lambda  \mu -2 \sqrt{6} \lambda  {x_1}+12 \sqrt{2} {y_1}+12 {y_2}-6\right),\\
 & {v_1}= \frac{1}{144} \left(-\sqrt{6} \left(\lambda ^2-12\right) \mu -12
   \lambda  \left(-2 \lambda  {x_1}+4 \sqrt{3} {y_1}+2 \sqrt{6} {y_2}+\sqrt{6}\right)\right)+{x_2},\\
 & {v_2}= -\frac{1}{48} \lambda  \left(\sqrt{6} (\lambda  \mu +18)-36 \lambda  {x_1}-24 \sqrt{3}
   {y_1}+12 \sqrt{6} {y_2}\right),\\
 & {v_3}= \frac{1}{24} \left(-\lambda  \mu +2 \sqrt{6} \lambda  {x_1}-12 \sqrt{2} {y_1}+12 {y_2}+18\right),
\end{align}
the linearization  transforms to
\begin{equation}
 \left(
\begin{array}{c}
 v_1'\\
 v_2'\\
v_3'\\
u_1'\\
u_2'\\
\end{array}
\right) =   \left(
\begin{array}{ccccc}
 -3 & 0 & 0 & 0 & 0 \\
 0 & -3 & 1 & 0 & 0 \\
 0 & 0 & -3 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)\left(
\begin{array}{c}
 v_1\\
 v_2\\
v_3\\
u_1\\
u_2\\
\end{array}
\right).
\end{equation}
Applying the CMT we obtain that the center manifold is given locally (up to second order) by the graph
\begin{align}
   \Bigg\{(u_1,u_2, v_1, v_2, v_3)\in\mathbb{R}^5: \quad   & v_1 = a_1 u_1^2 + a_2 u_1 u_2 + a_3 u_2^2 + \mathcal{O}(3), \nonumber \\
    & v_2 = b_1 u_1^2 + b_2 u_1 u_2 + b_3 u_2^2 + \mathcal{O}(3), \nonumber \\
    & v_3 = c_1 u_1^2 + c_2 u_1 u_2 + c_3 u_2^2 + \mathcal{O}(3)\Bigg\},
\end{align}
where
\begin{align}
& a_1= 0, a_2= \frac{1}{\sqrt{6}}, a_3= 0, b_1=0, b_2= 0, b_3= 0, c_1=-\frac{1}{192}, c_2= 0, c_3= \frac{1}{4}.
\end{align}
Hence, the dynamics on the center manifold is given locally (up to second order) by
\begin{align}
   &u_1'= \frac{1}{192} {u_1}^2 \left({u_1} \left(-\lambda ^2+24 {u_2}-12\right)+24 (\lambda +2 \lambda  {u_2})\right) + \mathcal{O}(3),  \label{center3a}\\
   & u_2'=-\frac{1}{32} {u_1} ({u_1}-8 \lambda  {u_2}) + \mathcal{O}(3) \label{center3b}.
\end{align}
According to Figure \ref{center3}, the origin is a saddle.
\begin{figure}[ptb]
\centering
\includegraphics[width=0.3\textwidth]{center3.pdf}
\includegraphics[width=0.3\textwidth]{center3b.pdf} \caption{Dynamics on the center manifold is given locally (up to second order) by \eqref{center3a}, \eqref{center3b} for $\lambda=-1$, (resp. $\lambda=1$). The dynamics is qualitative the same for other values of $\lambda<0$ (resp. $\lambda>0$).}%
\label{center3}%
\end{figure}

\textbf{Special case $y_{1c}=1$}. Using the transformation
\begin{equation}
    u_1= \mu ,u_2= y_2, v_1= \frac{1}{2} \left(-2 \sqrt{6} \lambda +2 \lambda ^2 x_1+2 \sqrt{6} \lambda  y_1-\sqrt{6} \lambda  y_2\right), v_2= \frac{1}{2} (2
   y_1-y_2-2), v_3= \frac{1}{3} \left(-\sqrt{6} \lambda +3 x_2-\sqrt{6} \lambda  y_2\right),
\end{equation}
the linearization  transforms to
\begin{equation}
 \left(
\begin{array}{c}
 v_1'\\
 v_2'\\
v_3'\\
u_1'\\
u_2'\\
\end{array}
\right)= \left(
\begin{array}{ccccc}
 -3 & 1 & 0 & 0 & 0 \\
 0 & -3 & 0 & 0 & 0 \\
 0 & 0 & -3 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)\left(
\begin{array}{c}
 v_1\\
 v_2\\
v_3\\
u_1\\
u_2\\
\end{array}
\right).
\end{equation}
Applying the CMT we obtain that the center manifold is given locally (up to second order) by the graph
\begin{align}
   \Bigg\{(u_1,u_2, v_1, v_2, v_3)\in\mathbb{R}^5: \quad   & v_1 = a_1 u_1^2 + a_2 u_1 u_2 + a_3 u_2^2 + \mathcal{O}(3), \nonumber \\
    & v_2 = b_1 u_1^2 + b_2 u_1 u_2 + b_3 u_2^2 + \mathcal{O}(3), \nonumber \\
    & v_3 = c_1 u_1^2 + c_2 u_1 u_2 + c_3 u_2^2 + \mathcal{O}(3)\Bigg\},
\end{align}
where
\begin{align}
& a_1=0, a_2=0, a_3=-\frac{1}{4} \sqrt{\frac{3}{2}} \lambda , b_1=0, b_2= \frac{\lambda }{12}, b_3= -\frac{1}{8}, c_1=0, c_2=\frac{\lambda ^2+3}{3 \sqrt{6}}, c_3=0.
\end{align}
Hence, the dynamics on the center manifold is given locally (up to second order) by
\begin{align}
   &u_1'=\frac{1}{24} u_1^2 \left(\left(\lambda ^2+3\right) u_1 u_2+6 \lambda  (u_2+1)\right)  + \mathcal{O}(3),  \label{center4a}\\
   & u_2'=-\frac{1}{2} \lambda  u_1 u_2 (3 u_2+1) + \mathcal{O}(3) \label{center4b}.
\end{align}
According to Figure \ref{center4}, the origin is a saddle.
\begin{figure}[ptb]
\centering
\includegraphics[width=0.3\textwidth]{center4b.pdf}
\includegraphics[width=0.3\textwidth]{center4.pdf} \caption{Dynamics on the center manifold is given locally (up to second order) by \eqref{center4a}, \eqref{center4b} for $\lambda=-1$, (resp. $\lambda=1$). The dynamics is qualitative the same for other values of $\lambda<0$ (resp. $\lambda>0$).}%
\label{center4}%
\end{figure}

\subsection{Equilibrium points for arbitrary potential}

Consider now an arbitrary potential $V\left(  \phi\right)  $, such that
$\Gamma\left(  \lambda\right)  $ to be arbitrary. In that consideration, the
dynamical system of our consideration has dimension six, where the equilibrium
points have coordinates $\bar{P}=\left(  x_{1}\left(  \bar{P}\right)
,x_{2}\left(  \bar{P}\right)  ,y_{1}\left(  \bar{P}\right)  ,y_{2}\left(
\bar{P}\right)  ,\mu\left(  \bar{P}\right)  ,\lambda\left(  \bar{P}\right)
\right)  $.

In particular the equilibrium points are $\bar{P}_{A}=\left(  x_{1}\left(
P_{A}\right)  ,x_{2}\left(  P_{A}\right)  ,y_{1}\left(  P_{A}\right)
,y_{2}\left(  P_{A}\right)  ,\mu\left(  P_{A}\right)  ,\lambda\left(
P_{A}\right)  \right)  $, where $A=1,2,3$ and $\lambda\left(  P_{A}\right)  $
solves the algebraic equation $\lambda^{2}\left(  \Gamma\left(  \lambda
\right)  -1\right)  =0$. While there are more equilibrium points because of the
different values of $\lambda\left(  P_{A}\right)  $, the physical physical
properties of the equilibrium points are exactly the same with the exponential
potential. That is an interesting result, since any kind of potential has a
similar physical behaviour.

Since the physical properties on the additional equilibrium points do not change, we omit the presentation of the stability analysis.

\section{Evolution of the Cosmographic parameters}
\label{sec6}

There are usually two known ways through which one tries to understand the dynamical history of our universe. One is the model independent approach in which one deals with the kinematical quantities that are
extracted  directly from the space-time metric and secondly, one needs to fix either the matter sector of the universe or the underlying gravitational sector and then explore the dynamics of the universe.
Indeed, the kinematical quantities are more flexible and novel compared to the model dependent approach since within the kinematical approach the dynamics of the universe is solely dependent on the geometrical quantities and hence this is more robust.  The study of the kinematical quantities derived from a homogeneous and isotropic universe is called the cosmography \cite{Weinberg-GR}. Using the variation of the kinematical quantities, one can measure the viability of a proposed cosmological scenario.

Under the background of the homogeneous and isotropic universe characterized by the FLRW universe, one can define the kinematical quantities as  \cite{Visser:2003vq,Bolotin:2018xtq}

\begin{eqnarray}
	&& H(t)  = \frac1a\frac{da}{dt}, \nonumber \\
	&&q(t) =  -\frac1a\frac{d^2a}{dt^2}\left[\frac1a\frac{da}{dt}\right]^{-2}, \nonumber \\
	&&j(t) =  \frac1a\frac{d^3a}{dt^3}\left[\frac1a\frac{da}{dt}\right]^{-3},\nonumber \\
	&&s(t) =  \frac1a\frac{d^4a}{dt^4}\left[\frac1a\frac{da}{dt}\right]^{-4},\nonumber
\end{eqnarray}
where $H$, $q$, $j$, $s$  are respectively called the Hubble, deceleration, jerk, snap  parameters and except for the Hubble parameter, the last three parameters are dimensionless. With these terminology, the expansion scale factor $a(t)$ can be written as

\begin{eqnarray}
    a(t) =  a_0 \left[1 + H_0 (t-t_0) - \frac{1}{2} q_0^2 (t-t_0)^2 + \frac{1}{3!} j_0 (t-t_0)^3 + \frac{1}{4!} s_0 (t-t_0)^4 + \mathcal{O} [(t-t_0)^5] \right],
\end{eqnarray}
in which any sub-index attached to any quantity refers to its present value.
We note that one can extend this series by allowing higher order derivatives of the scale factor beyond order four. Here we are interested to investigate the deceleration parameter, and its next two hierarchies, namely, jerk and the snap parameter which in terms of the Hubble parameter can alternatively be rewritten as

\begin{align}
q  & =-1-\frac{\dot{H}}{H^{2}}\\
j  & =\frac{\ddot{H}}{H^{3}}-3q-2\\
s  & =\frac{\dddot{H}}{H^{4}}+4j+3q\left(  q+4\right)  +6,
\end{align}
where an overhead dot represents the cosmic time derivative; therefore, two or three overhead dots represents the two- or three-times derivative with respect to the cosmic time.
Thus, from the evolution of $q, j, s$, one can understand the expansion of the universe (accelerating or decelerating), the rate of acceleration ($j$) and its next derivative.
In Fig. \ref{fig3} we present the qualitative evolution of the cosmographic  parameters.

\begin{figure}[ptb]
\centering\includegraphics[width=0.8\textwidth]{gup3.pdf} \caption{Qualitative
evolution of the cosmographic parameters $q,~j$ and $s$, for the same initial  conditions as used for Fig. \ref{fig1}. Right Figs are subplots of the Left Figs in the  range near the de Sitter attractor. }%
\label{fig3}%
\end{figure}

\section{Conclusions}
\label{sec7}

The dynamics of the universe both at its very early and late evolution are two mysterious issues in cosmology.  Various approaches have been proposed to understand both the phases of the universe but none of them are found to be extremely satisfactory.  The classical theory of gravity does not work out during the very early evolution of the universe but the quantum gravity does a bit, according to the available quantum gravitational theories. The quantum gravity theory predicts
the existence of a minimum length scale of the order of the Planck length $l_{\rm pl}$ and consequently, this motivated to generalize  the
Heisenberg's Uncertainty Principle to some Generalized Uncertainty principle (GUP) in quantum gravity.
The Generalized Uncertainty principle has been found to be very effective to explain several important issues that are directly related to our universe and its evolution. In particular, in the context of black hole evaporation, GUP has a significant role. Moreover, GUP has a explanation towards the existence of the dark matter fluid in the universe's sector. Thus, it will be interesting whether the GUP modified cosmological models could give rise to some interesting features of the universe.
Hence, in the present work we study the modified quintessence model and investigate the stability of the model by performing its dynamical analysis.

We first investigate the modified quintessence model for the exponential potential which gives rise to various possibilities. For various critical points,  we get different solutions such as the matter dominated universe and the universe dominated by the kinetic term only, however, all of them are found to be unstable in nature. Only the family of points on the surface with coordinates $P_3$ describe the de Sitter universes and they all are attractors. It is interesting to note that although no scaling solutions or a tracking solution exist in this specific model contrary to the quintessence
theory, however, the critical points describing the de Sitter  solution do not exist in the case of a quintessence model with an
the exponential potential.  This is a new result in this work. In order to understand the behavior of the family of critical points, $P_3$, we applied the Center Manifold Theorem.

Furthemrore, we considered an arbitrary potential and performed the dynamical system analysis similar to the exponential potential. We found that the  physical properties of the equilibrium points for the general potential remain exactly same as in the case with the exponential potential. That is a very striking and exciting result in this work, because this clearly indicates that the investigations with the exponential potential is enough for the modified quintessence models with GUP.

Finally, we observe that the parameter for the equation of state for the effective cosmological fluid can cross the phantom divine line and as it was observed by the numerical simulations without the appearance of ghosts, while the late time attractor is the de Sitter universe.

\begin{acknowledgments}
AG was funded by FONDECYT grant number 1200293. The research of AP and GL was funded by Agencia Nacional de Investigaci\'{o}n y
Desarrollo - ANID through the program FONDECYT Iniciaci\'{o}n grant no.
11180126. Additionally, by Vicerrector\'{\i}a de Investigaci\'{o}n y
Desarrollo Tecnol\'{o}gico at Universidad Catolica del Norte. SP was supported
by the  Mathematical Research Impact-Centric Support Scheme (MATRICS), File No.
MTR/2018/000940, awarded by the Science and Engineering Research Board (SERB),  Govt. of India.
\end{acknowledgments}


\section*{Acknowledgement}
This work was supported in part by The National Key Research and Development Program of China under grant 2018YFB1800804, the National Nature Science Foundation of China under U1936217,  61971267, 61972223, 61941117, 61861136003, Beijing Natural Science Foundation under L182038, Beijing National Research Center for Information Science and Technology under 20031887521,  research fund of Tsinghua University - Tencent Joint Laboratory for Internet Innovation Technology, and The Information Core Technology Center at Institute for Interdisciplinary.

\bibliographystyle{ACM-Reference-Format}

\section{Introduction}
\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\textwidth]{figures/overview.png}
    \caption{The overview of reinforced epidemic control system.}
    \label{fig:Overview}
    \vspace{-0.5cm}
\end{figure}

The epidemic has always been a threat to human society by exposing us in front of a dilemma between saving lives or economy. The virus infects gathering people and spreads through daily commute \cite{poletto2012heterogeneous, balcan2010modeling, wesolowski2012quantifying}. Controlling the spread of the virus must cut off daily mobility, which is a pillar of the modern economy. For instance, the recent outbreak of COVID-19 has caused millions of infections and hundreds of thousands of death tolls. The epidemic forces many municipal governments to issue  a stay-at-home order, which is a Fully LockDown (FLD) policy. FLD in most cities lasts for weeks thus deeply hurts the economy~\cite{barua2020understanding}. Some municipalities try to only quarantine symptomatic people and their close contacts at the early stage of the epidemic. However, this infected-individual-quarantine policy would be only implementable when governments are able to accurately and comprehensively trace risky people. It is also unreliable when there exist many asymptomatic infected people. Current computer-science explorations pursue using smartphone data to infer and trace highly-risky people \cite{ferretti2020quantifying, oliver2015mobile}. However, fully tracing individual mobility and contacts requires full coverage of smartphones and further raises the concern of threatening privacy~\cite{privacy}. According to an investigation by the University of Maryland and The Washington Post, around 60\% of respondents either prefer not sharing their private information or do not own a smart phone~\cite{timberg2020most}.
In summary, the vast amount of complex individual mobility and asymptomatic infected people prevent current epidemic-control policies from cutting off virus spread without hurting the economy when private information cannot be fully captured.

We in this research demonstrate that a smart epidemic control policy is still available even if private mobility information is unavailable. We develop a dynamic control framework to avoid an epidemic outbreak by limiting the probability of risky mobility's occurrence. Instead of targeting and limiting risky individual's mobility according to private data, our framework estimates each urban region's risk of having a high infected population and uses the estimation to control inter-regional mobility. Highly-risky inter-regional mobility will be limited to suppress the probability of infected people's movement. Because the infected people are a small proportion of the population even in a seriously infected city, only a small number of mobility must be restricted. It is possible to avoid an epidemic outbreak by heterogeneously limiting little inter-regional mobility. Furthermore, the estimation is based on the regional aggregate demand for mobility and the regional epidemic statistics. Thus, private data is dispensable.

However, there exist three specific challenges causing the complexity of estimating and controlling inter-regional mobility for suppressing the infection and protecting the economy.
First, urban mobility is vast and temporally varying, making it hard to target the really risky mobility.
Further, the requirements of the policy's practicality sophisticate the design of the epidemic-control policy. An implementable control policy cannot continuously quarantine the same urban region for too long.
Last but not the least, the search for policy is difficult. Due to the exponentially-increasing nature of epidemics, the number of future infections is a highly non-convex function of each previous decision, making it hard to explore the policy space. Furthermore, the dual objectives cause the policy exploration often end up stuck in local optimums, which is also exacerbated by the non-convexity of infections.

With the consideration of the above challenges, we develop a DUal-objective Reinforcement-Learning Epidemic Control Agent (DURLECA) framework by combining Graph Neural Network (GNN) and Reinforcement Learning (RL) approach, to search out an effective mobility-control policy.
DURLECA hires a GNN to estimate the virus-transmission risk induced by urban mobility, which is a dynamic flow on a graph. Based on the estimated risk, the RL agent periodically determines the extent of the restriction on each inter-regional mobility.
The GNN of DURLECA is developed with a novel architecture, namely Flow-GNN, to fit the virus spread process on mobility flows, which existing GNN architectures are incompatible to characterize.
We also carefully construct a reward function for the RL agent to precisely capture the natural trade-off relation between epidemic control and urban-mobility retaining. The reward function also considers the difference between continuous and intermittent restrictions on the same region.
Furthermore, we develop two RL exploration strategies that appropriately incorporate epidemic expert knowledge for guiding and stabling policy exploration.

Supported by a Susceptible-Infected-Hospitalized-Recovered (SIHR) epidemic simulation environment developed from the traditional SIR model~\cite{SIR}, DURLECA is able to successfully search out a mobility-control policy that suppresses the epidemic and retains most of the mobility. Our experiments on a real-world mobility dataset collected in Beijing demonstrate the effectiveness of DURLECA. Even if the city starts to suppress an epidemic\footnote{The $R0$ is $0.19 \sim 1.08$ for SARS, $1.4 \sim 2.8$ for Influenza, $1.94 \sim 5.7$ for COVID-19, according to \url{https://en.wikipedia.org/wiki/Basic_reproduction_number}.} whose $R0=2.1$ after 20 days of discovering the first patient, DURLECA still finds out a policy where:
\begin{itemize}

    \item
    The peak demand for hospitalization is under 1.3\textperthousand\footnote{The hospital bed density is 2.9 \textperthousand \text{ in} U.S., 4.2\textperthousand  \text{ in} China, and 13.4\textperthousand  \text{ in} Japan, according to \url{https://www.indexmundi.com/g/r.aspx?v=2227&l=en}.} of the whole population. The average demand for hospitalization is controlled under 0.4\textperthousand.
    \item 76\% of the total mobility is retained. In more than $70\%$ intervened days, two-thirds regions retain over $70\%$ mobility. No region ever experiences a stringent day, i.e., daily retained mobility lower than 20\%.
\end{itemize}

In summary, the contribution of this paper is in three-folds:
\begin{itemize}
   \item
   We bypass the privacy concern for smart epidemic control. Instead of directly tracing and quarantining risky individuals, we suppress the risk of an epidemic outbreak by estimating and restricting risky inter-regional aggregate mobility.
   \item
   We develop DURLECA to dynamically generate customized control actions for inter-regional mobility, which allows a smart solution for the life-or-economy dilemma of epidemic control.
  \item
    We propose innovative approaches to guarantee DURLECA's capability. We design a novel GNN architecture that can fit the epidemic transmission dynamics. Our RL reward function captures the nature of the trade-off relation between epidemic suppressing and mobility retaining, and reflects practical requirements. We also develop two RL exploration strategies that appropriately incorporate epidemic expert knowledge for guiding and stabling policy exploration.
\end{itemize}

\section{Preliminary and Problem Formulation}
During the stay-home order of COVID-19, governments distribute mobility quotas per day to each household for retaining the basic economic activities, such as people's procurement for food. According to the current quota regulation, we develop a new policy environment. We assume that the government periodically predicts or collects aggregate demands for inter-regional mobility of every Origin-Destination (OD) pair. The government also collects information about the number and location of current discovered patients.
Those pieces of information are used to determine the quotas for each inter-regional mobility.
The quota-distribution aims at minimizing the risk of epidemic break out in the foreseeable future periods while maximizing the mobility demands. We in this section present the modeling of mobility and epidemic that supports quota allocation.

\subsection{Mobility Modeling}
We model a city's urban-mobility demand at time step $\tau$ as a mobility matrix $M_d^\tau$, whose element $M_{i,j}^\tau$ represents the inter-regional mobility demand, i.e., the number of people who demand to move, from $i$ to $j$. According to $M_d^\tau$ and the epidemic information, the city government determines a mobility quota matrix $p^\tau$ at $\tau$, whose element $p_{i,j}^\tau$ is the quota rate distributing to the mobility demand from $i$ to $j$. Therefore, the allowed inter-regional mobility denoted by $M_{p,i,j}^\tau$ is calculated according to the following equations.
\begin{align}
    & M_{p,i,j}^\tau = p_{i,j}^\tau M_{d,i,j}^\tau ,  \\
    & M_p^\tau = \mathcal{T}(M_d^\tau, p^\tau) = M_{d}^\tau \odot p^\tau,
\end{align}
where $\mathcal{T}$ refers to the mobility control function and  $\odot$ denotes for element-wise multiplication. Note that $M_d^\tau$, $M_p^\tau$, and $p^\tau$ are $K \times K$ matrices, where $K$ is the number of regions in the studied city. We summarize the mobility-related notations in Appendix.

\subsection{Epidemic Modeling}
The main challenge for urban epidemic control comes from infected people who are infectious but asymptomatic. Therefore, we develop a new epidemic model to capture the difference between asymptomatic people and symptomatic people. Our model is based on the traditional Susceptible-Infected-Recovered ($SIR$) model in public-health literature \cite{SIR}. We introduce a new state beyond $SIR$ and denote it by Hospitalized ($H$). People in state $H$ are infected with symptoms and thus will be quarantined or hospitalized. They will not participate in urban mobility and will not contribute to new infections.
We refer our model as $SIHR$ model\footnote{This modeling is different from the SEIR model \cite{SEIR} which assumes the asymptomatic people are not infectious.}.

We use our $SIHR$ model to capture the dynamic process of infection spread over urban mobility. We denote region $i$'s epidemic state by $E^\tau_i = \{S^\tau_i,I^\tau_i,H^\tau_i,R^\tau_i\}$, whose each element respectively denotes the susceptible, infected, hospitalized, and recovered population of $i$ at $\tau$. We use $E_{v,i}^\tau =\{S^\tau_i+I^\tau_i,H^\tau_i,R^\tau_i\}$ to represent the visible state of $i$ at $\tau$, where the healthy people cannot be differentiated from infectious asymptomatic people. We denote the total population of $i$ at $\tau$ by $N_i^\tau$.

The epidemic state $E_i^\tau$ is updated in each time step. For each time step $\tau$, we separate $\tau$ into two sub-steps: mobility happens and infection occurs. At the mobility-happening sub-step, people accomplish their moves between regions. We use $E_i^{s,\tau}$ to represent the epidemic state of the staying people while $E_i^{m,\tau}$ represents the new arrival's. The overall epidemic state at the mobility-happening sub-step, denoted as $\hat{E}_i^{\tau}$, is calculated as follows:
\begin{align}
    & E_i^{s, \tau} = E_i^{\tau} - \sum_j \frac{M_{p,i,j}^{\tau}}{N_i^{\tau}}E_i^{\tau}, \label{eq:flow_stay}\\
    & E_i^{m, \tau} =  \sum_j \frac{M_{p,j,i}^{\tau}}{N_j^{\tau}}E_j^{\tau}, \label{eq:flow_movein}\\
    & \hat{E}_i^{\tau} = E_i^{s, \tau} + E_i^{m, \tau}. \label{eq:flow_end}
\end{align}

At the infection-occurring sub-step, people that stay at $i$ infect each other. Simultaneously, new arrivals at $i$ infect each other. Therefore, the epidemic state is updated as follows:
\begin{align}

    & S_i^{\tau+1} = \hat{S}_i^{\tau} -
    \frac{\beta_i^{s,\tau} S_i^{s,\tau} I_i^{s,\tau}}{N_i^{s,\tau}} -
    \frac{\beta_i^{m,\tau} S_i^{m,\tau} I_i^{m,\tau}}{N_i^{m,\tau}} , \\

    & I_i^{\tau+1} = \hat{I}_i^{\tau} +
    \frac{\beta_i^{s,\tau} S_i^{s,\tau} I_i^{s,\tau}}{N_i^{s,\tau}} +
    \frac{\beta_i^{m,\tau} S_i^{m,\tau} I_i^{m,\tau}}{N_i^{m,\tau}} - \gamma \hat{I}_i^{\tau}, \\
    %
    & H_i^{\tau+1} = H_i^{\tau} + \gamma_i^\tau \hat{I_i}^{\tau} \label{eq:hospital} - \theta_i^\tau H_i^{\tau},\\

    & R_i^{\tau+1} = \hat{R_i^{\tau}} + \theta_i^\tau H_i^{\tau}, \label{eq:R}
\end{align}
where $\{\hat{S}_i^{\tau},\hat{I}_i^{\tau},\hat{H}_i^{\tau},\hat{R}_i^{\tau}\}$ are elements of $\hat{E}_i^{\tau}$. $\{ \beta^{s,\tau}_i, \beta^{m,\tau}_i \}$ are the epidemic's transmission rate for the staying people and the moving people respectively. $\gamma^\tau_i$ is the hospitalized rate and $\theta^\tau_i$ is the recover rate. We use one set of $\{\beta^{s}, \beta^{m}, \gamma,\theta \}$ for all regions at all time steps for simplification. We introduce how we estimate $R0$ in Appendix.

\subsection{Multi-Objective Sequential Control Problem Formulation}
The above mobility and epidemic modeling allow us to formulate the dynamic inter-regional mobility control problem for minimizing infections and maximizing mobility retaining, shown in Equation (\ref{eq:trans_model})$\sim$(\ref{eq:reward}):
\begin{align}
    & M_p^\tau = \mathcal{T}(M_d^\tau, p^\tau) \label{eq:trans_model}, \quad E_p^{\tau+1} = \mathcal{E}(M_p^\tau, {E}_p^\tau) \\
    & P^{t,T} = \mathrm{arg}\,\underset{P}{\mathrm{max}} \sum_{\tau=t}^T \mathcal{O}(M_P^\tau, E_P^\tau), \label{eq:reward}
\end{align}
where $\mathcal{O}$ is the objective function, satisfying $\partial^2 \mathcal{O}/\partial M_p \partial E_p < 0$ because of the trade-off nature between epidemic control and mobility retaining. $\mathcal{O}$ should also meet some practical requirements.
In the next section, we detail the design of the objective function and use it as the reward function of the RL module of DURLECA. Besides, we particularly consider the fact that the frequency of government interventions is lower than the frequency of mobility. Therefore, the mobility and infection updates per hour while the government determines mobility quotas per four hours.

\section{DURLECA}
DUal-objective Reinforcement-Learning Epidemic Control Agent (DURLECA) is a GNN-enhanced RL agent to estimate regional infection risk and determine mobility quota. An overview of DURLE-CA is shown in Figure~\ref{fig:system}.  At each time step $\tau$, DURLECA acquires an observation $E_v^\tau$ from the environment. According to $E_v^\tau$ and the demand mobility $\{M_d^\tau,M_d^{\tau+1},M_d^{\tau+2},M_d^{\tau+3}\}$, our RL agent gives a control action $p^\tau$ for the optimization problem in Equation (\ref{eq:reward}).
In the rest of this section, we provide the details of DURLECA.

\begin{figure}[t]
\vspace{-0.2cm}
    \centering
    \includegraphics[width=0.5\textwidth]{figures/system.png}
    \caption{The details of the proposed DURLECA.}
    \label{fig:system}
    \vspace{-0.4cm}
\end{figure}

\subsection{Reinforcement Learning}

We now re-formulate the multi-objective sequential control problem using the basic factors of RL, i.e., state, action, reward, and learning algorithm.

\textbf{State:} We take the visible epidemic state $E_v^\tau$, its temporal one-order derivatives $\nabla_\tau E_v^\tau$, the mobility demand $\{M_d^\tau,M_d^{\tau+1},M_d^{\tau+2},M_d^{\tau+3}\}$, and the historical mobility loss $L^\tau$ (defined later) as the state for RL.

\textbf{Action:} The action of RL is defined as the mobility restriction $p^\tau$ determining the quota rate for each inter-regional mobility at $\tau$. Each element of $p^\tau$ is a real number between 0 and 1.

\textbf{Reward:} The reward function is designed to reflect the objective $\mathcal{O}$ of the optimization problem in Equation (\ref{eq:reward}). It includes two terms: an infection-spread-cost term and a mobility-restriction-cost term. In order to guide the RL agent to effectively find an effective and practical mobility-control policy. We design the reward function to satisfy the following three requirements:
\begin{itemize}
    \item Reflecting the trade-off relation between infection control and mobility retaining.

    \item Capturing the exponential growth of the social cost caused by infection spread. The social cost is low when the infected population is small. However, the social cost will skyrocket once the infection population exceeds the capacity of the city's healthcare system.
    \item Penalizing continuous mobility restrictions in the same region. People's tolerance for mobility restrictions is limited. Thus, the reward function has to include a growing penalty for continuously restricting the same region.
\end{itemize}
According to the above three requirements, we separately design the infection-spread-cost term and mobility-restriction-cost term. We denote the infection-spread-cost by $R_h^\tau$ and model it as follows:
    \begin{align}
        R_h^\tau = k_h \text{exp}(\frac{\frac{1}{K}\sum_i H_i^\tau}{H_0}), \label{eq:infection}
    \end{align}
where $k_h$ is a hyper-parameter determining the start-up social cost of a city having the first patient while $H_0$ is a hyper-parameter determining how the social cost increase along with the number of patients. The mobility-restriction-cost term is denoted by $R_m^\tau$ and defined below:
\begin{align}
        & L_i^\tau = \sum_{t=0}^{\tau-1} \lambda^{\tau - t} \frac{M_{d,i}^{t}- M_{p,i}^{t}}{ \overline{M_{d,i}}}, \\
        & R_m^\tau = \frac{1}{K} \sum_i \text{exp}(\frac{L_i^\tau}{L_0}) \frac{M_{d,i}^{\tau}- M_{p,i}^{\tau}}{\overline{M_{d,i}}}. \label{eq:restriction}
    \end{align}
Here, $L_i^\tau$, the historical mobility loss, is the amount of mobility restricted in history and induces an exponentially-growing penalty $R_m^\tau$ on the current restriction. The hyper-parameter $\lambda$ determines the discount rate of historical restrictions' impacts. The hyper-parameter $L_0$ determines how large the penalty is for continuously limiting the same region. Finally, we develop the RL reward function $\mathcal{R}$ as follows:
    \begin{align}
        \mathcal{R}(M_p^{\tau}, E_p^{\tau}) = -(R_m^\tau + R_h^\tau).
    \end{align}
Note that our design enables the reward function to reflect that the infection-spread-cost booms once the whole city's hospitalized population exceeds the city's healthcare system capacity while the mobility-restriction-cost skyrockets if any single region is continuously restricted for multiple periods.

\textbf{Learning Algorithm:}
DURLECA employs a Deep Deterministic Policy Gradient (DDPG) \cite{DDPG} agent to search for mobility-control policy because the action space is continuous. The DDPG agent is composed of a critic network and an actor network. The critic network aims to estimate the expected reward gained by a control action. The actor network searches for the best action, which gives quota rates for all inter-regional mobility, by maximizing the critic network's output. We use Parameter Noise \cite{param_noise} to improve exploration during RL training.

\subsection{Flow-GNN}
Both critic and actor networks have to well capture the graph nature of urban mobility, where regions are nodes connected by OD flows.
Therefore, we adopt GNN to develop both of them. We design a novel GNN architecture so that GNN can characterize the epidemic transmission process driven by regional infection aggregation upon inter-regional mobility.
We refer our proposed GNN as Flow-GNN, which is developed on the basis of GraphSage~\cite{graphsage}.

In particular, we design Flow-GNN fit for the low-frequency mobility control associated with high-frequency mobility dynamics. Considering that we determine mobility quota per four hours, we include 4 Flow-GNN layers in our network and input edge information chronologically. The edge-input information for the $k$-th layer is $M_*^{\tau+k-1}$. We use $f_i^{k}$ to denote the feature of region $i$ outputted by the $k$-th GNN layer and calculate it according to the following equations:
\begin{align}
    & f_{stay,i}^{k-1} = (1-\sum_j \frac{M_{*,i,j}^{\tau+k-1}}{N_i^{\tau+k-1}})f_i^{k-1}, \label{eq:stay_feature}\\
    & f_{in,i}^{k-1} = \sum_j \frac{M_{*,j,i}^{\tau+k-1}}{N_j^{\tau+k-1}}f_j^{k-1}, \label{eq:in_feature} \\
    & f_i^k = \sigma (W^k (f_{in,i}^{k-1}, f_{stay,i}^{k-1}) + B^k). \label{eq:all_feature}
\end{align}
Here $(f_{in,i}^{k-1}, f_{stay,i}^{k-1})$ denotes for concatenation, $\sigma$ is a non-linear activation function, and $W^k, B^k$ are trainable parameters. Specifically, we input the first layer with $f_i^0 =\{E_{v,i}^\tau, \nabla_\tau E_{v,i}^\tau \}$.

The above equations correspond to our modeling of epidemic transmission in Section 2.2, where we separate each time step $\tau$ into mobility-happened sub-step and infection-occurred sub-step. Equation (\ref{eq:stay_feature}) describes the epidemic feature of staying population at mobility-happened sub-step while Equation (\ref{eq:in_feature}) represents the new-arrival population's epidemic feature in the same sub-step. Equation (\ref{eq:all_feature}) characterizes the epidemic transmission in the staying population and the new-arrival population.

\subsection{Exploration Strategies}
The exponentially-increasing nature of epidemics and our dual objectives cause difficulties for RL exploration and increase the agent's risk of falling into local optimums.
We design two RL exploration strategies to address this problem. The first strategy is to incorporate pseudo-expert knowledge to improve RL searching efficiency. The second is to protect the agent from falling into local optimums by stopping it from exploring apparently unreasonable policies.

\textbf{Generating-and-Incorporating Pseudo Expert:} We can generate simple but dynamic policies according to current epidemic-management experience, which can be a good start point for RL exploration. For instance, most cities currently restrict the mobility of regions with a large symptomatic population while a region has urgent reopening demand if it has been continuously locked down for a long time. Thus, we design a pseudo expert, which control $p_{i,j}^\tau$ as follows:
\begin{align}
    p_{i,j}^\tau =
    \begin{cases}
        0 & \quad \text{if } H_i^\tau > X_h \text{ and } L_i^\tau < X_l  \\
        1 & \quad \text{else} .
    \end{cases} \label{eq:expert}
\end{align}
The expert will lock a region down based on two conditions: 1) the number of hospitalized, or symptomatic patients in this region, exceeds the threshold $X_h$; 2) this region has not been restricted very much in history, reflected by that $L_i^\tau$ does not exceed the threshold $X_l$. During testing, this expert is also used as a comparing baseline.

We let the agent first explore with expert's guidance and then gradually learn to explore by itself to outperform the expert. The idea is inspired by the approach adopted to develop AlphaGO~\cite{go}. Specifically, we set an adaptive probability for the agent to directly choose the expert action instead of taking an action by itself during training. This design enables the agent to compare the pseudo-export strategy with its own, which avoids the agent to move towards inefficient directions at the initial stage of training. The adaptive probability decreases along with training steps, which enables the agent to broadly explore and outperform the expert at the later stage of training.

\textbf{Avoiding Extreme Points:} The wide exploration might lead the agent to fall into some extreme points. The training might be unstable due to a sudden large loss caused by a poor control action.
Meanwhile, the strong incentive of avoiding the large loss will force the agent to fall into local optimal control policies, such as a forever fully-lockdown.
To avoid such extreme points, we set two rules:
\begin{itemize}
    \item The infection threshold $I_{t}$: If the agent explores into a state where the regional mean number of infected people exceeds $I_{t}$, it will end the episode and receive a large penalty.
    \item The lockdown threshold  $L_{t}$: If the agent explores into a state where there exists a region $i$ that $L_i^\tau$ exceeds $L_{t}$, it will end the episode and receive a large penalty.
\end{itemize}
The two rules are straightforward but effective to help the RL agent avoid potential local optimums.

\section{Experimental Evaluations}
In this section, we conduct extensive experiments to answer the following research questions:

\textbf{RQ1:} Can DURLECA resolve the life-or-economy dilemma?

\textbf{RQ2:} Can DURLECA adapt to both early intervention and late intervention?

\textbf{RQ3:} Can DURLECA be generalized to different cities and different diseases?

Besides, we conduct ablation studies in Appendix to evaluate the effectiveness of our proposed Flow-GNN and RL exploration strategies.

\subsection{Dataset}
We use a real-world OD dataset collected by a mobile operator in Beijing to evaluate DURLECA. The dataset divides Beijing into $17 \times 19$ regions and covers 544,623 residents. Averagely, each region has $1686$ observed residents. The dataset covers 24-hour OD-flows for the whole month of January 2019. We repeat the one-month data 24 times and get a prolonged dataset of 24 months so that we have a sufficiently long period for discussing epidemic control. We list other details in Appendix.

\subsection{Metrics and Settings}
We design six metrics to evaluate the performance of DURLECA on resolving the life-or-economy dilemma. We introduce the metrics in the following and summarize them in Table \ref{tab:metrics}.

\begin{table}[ht]
\vspace{-0.3cm}
\centering
\begin{tabular}{l|l|l}
\hline
\hline
Metric & Value & Physical Meaning \\
\hline
Mean/Max $H$ & 0-1686 & Temporal mean/max of $H$\\
Total $R$ & 0-1686 & Total $R$ after the epidemic \\ \hline
$Q$ & 0-1 &  Total quota rate \\
$T_{20\%}^c$ & 0-744 & The city 20\%-mobility duration\\
$T_{20\%}^r$  & 0-744 & The region 20\%-mobility duration \\
\hline
\hline
\end{tabular}
\caption{The summary of metrics and related value ranges.}
\label{tab:metrics}
\vspace{-0.7cm}
\end{table}

We select three metrics to assess the epidemic-suppressing performance of an epidemic control policy, including the total number of infected people that is equal to $R$ at the end of the epidemic period, the mean number of hospitalized people whose value is the mean of $H$ over time, and the peak demand for hospitalization capacity that is equal to the max value of $H$ over time. Total $R$ determines the total social medical costs while both Mean/Max of $H$ reflect the sustained and peak pressure on the healthcare system. We also select three metrics to assess the mobility-retaining performance of an epidemic control policy, including the total ratio of retained mobility $Q$, the duration of stringent mobility restrictions on the whole city $T_{20\%}^c$, and the duration of stringent mobility restrictions on the most restricted region $T_{20\%}^r$.

\textbf{Epidemic Settings:} Without the loss of generality, we set $\beta_s = \frac{0.3}{24}, \beta_m=\frac{3}{24}, \gamma=\frac{0.3}{24}, \theta = \frac{0.3}{24}$ in most of our experiments. The estimated basic reproduction number $R_0$ is 2.1.

\textbf{Intervention time:} We define $t_{start}$ as the time when the policymakers discover the epidemic and start to intervene. In our experiments, we compare results with $t_{start}=0,10,20$.

For more details about our experiment settings, please refer to Appendix.
Our implementation is available online at \url{https://github.com/anyleopeace/DURLECA}.

\subsection{Performance Comparison}

\textbf{Baselines:}
We set four different expert baselines to simulate different real-world expert policies and compare them with DURLECA on resolving the life-or-economy dilemma.
\begin{itemize}
    \item \textbf{EP-Fixed:} In the real world, a simple but inflexible control is to restrict all mobility in the city. For simulation, we design EP-Fixed to give a fixed quota rate $X_q$ to all inter-regional mobility during the whole epidemic period. We set  $X_q=\{0.15,0.2\}$ in our experiments, as we find them at the boundary of successfully controlling the epidemic.

    \item \textbf{EP-Soft:} We design an expert baseline following Equation (\ref{eq:expert}), which softly depends on the historical mobility loss $L_i^\tau$ and the current hospitalized population $H_i^\tau$ to determine whether to lock down a region. We set $X_h=0,X_l=168$ in our experiments. $X_h=0$ guarantees the expert receive equivalent information compared with DURLECA. $X_l=168$ corresponds to the real-world control policy in some countries: a continuous 7-day (168-hour) lockdown.

    \item \textbf{EP-Hard:} Without softly depending on the historical mobility loss, an expert can reopen a region if it has been locked down for successive $X_t$ days. This expert, namely EP-Hard, gives daily quota as follows:
    \begin{align}
    p_{i,j}^\tau &=
    \begin{cases}
        0 & \quad \text{if } H_i^\tau > X_h \text{ and } \sum_{t=1}^{X_t} M_{p,i}^{\tau-t} > 0  \\
        1 & \quad \text{else} .
    \end{cases} \label{eq:expert_hard}
    \end{align}
    We set $X_h=0,X_t=7$ for a similar reason of EP-Soft.

    \item \textbf{EP-Lockdown:} The most robust and conservative policy is to lock down the whole city until the epidemic ends. To simulate it, we design an expert following Equation (\ref{eq:expert}) but with $X_h=0,X_l=inf$. It can lock down a region for an any-long time until the hospitalized population becomes zero.
\end{itemize}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.46\textwidth]{figures/sim.png}
    \caption{The simulation visualization of DURLECA and selected baselines when $t_{start}=20$. Main Figure: The number of hospitalized people along with time. Upper Right: The histogram for $Q_i^\tau$. Best viewed in color.}
    \label{fig:results}

\end{figure}

\begin{table}[ht]
\scalebox{0.95}{
    \centering
        \begin{tabular}{l|lllll}
        \hline
        \hline
          & Mean/Max $H$ & Total $R$ & $Q$ & $T_{20\%}^c$  & $T_{20\%}^r$  \\
        \hline

        No Intervention   & 27.21/157.55 & 1069.29 & 1 & 0 & 0 \\ \hline
        EP-Fixed $20\%$   & 4.03/17.96 & 877.32  & 0.20 & 724 & 724 \\
        EP-Fixed $15\%$    & 0.44/1.55 & 10.18  & 0.15 & 724 & 724 \\

        EP-Soft & 4.66/53.80 & 1040.68 & 0.57 & 9 & 19\\

        EP-Hard  & 0.45/1.78 & 8.31 & 0.13 & 36 & 36 \\
        EP-Lockdown  & 0.41/1.42  & 5.75   & 0 & 27 & 27 \\ \hline

        DURLECA  & 0.60/2.28 & 19.07 & 0.76 & 0 & 0 \\
        \hline
        \hline
        \end{tabular}
        \vspace{-0.7cm}
    }
    \caption{The simulation results of DURLECA and all baselines when $t_{start}=20$.}
\label{tab:results}
\vspace{-0.5cm}
\end{table}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.53\textwidth]{figures/20_results_quotas.png}

    \caption{The spatial distribution and histogram of $Q_i^\tau$ given by DURLECA. We select four periods from a "high-risk day" and a "low-risk day". Each grid in one of the four 17$\times$19 maps represents a region in Beijing. Each histogram summarizes the distribution of quota rates in the respective period.}
    \label{fig:quota}
    \vspace{-0.3cm}
\end{figure}

\textbf{Results and Analysis:}
We compare DURLECA with all baselines when $t_{start} = 20$ in Table~\ref{tab:results}. We also visualize three selected results in Figure~\ref{fig:results}. Expert baselines can achieve only one goal in the life-or-economy dilemma, while DURLECA can achieve both.

EP-Soft can retain 57\% of the total mobility. However, it leads to an epidemic outbreak, reflected by the super-high value of Mean/Max $H$ and Total $R$. The healthcare system will break down.
EP-Fixed ($X_q=15\%$), EP-Hard and EP-Lockdown can keep Mean/Max $H$ at a low level so that the healthcare system will not be overwhelmed. However, the low value of $Q$ indicates that all of them fail to retain mobility. The large value of $T_{20\%}^c$ and $T_{20\%}^r$  demonstrates that some regions and the city have to experience long-term lockdown, which is an unacceptable damage to the economy. Besides, the differences of EP-Fixed (15\%) and EP-Fixed (20\%) in Mean/Max $H$ and Total $R$ also indicate that the expert control is very vulnerable to mobility perturbation. The above results also manifest that all those expert policies fail at resolving the life-or-economy dilemma of epidemic control.

Compared with those baselines, DURLECA simultaneously suppresses the epidemic and retains a large amount of mobility. DURLE-CA achieves low values of Mean/Max $H$, which guarantee the demands for hospitalization will not exceed the capacity of most countries. DURLECA also suppresses the total infected population at a low level, about $1\%$ of the total population. The red curve in Figure \ref{fig:results} presents the performance of DURLECA in epidemic suppression.

DURLECA also retains the most mobility. 76\% of the total mobility during the intervened period is retained. Furthermore, no regions will be fully locked down. DURLECA retains 70-100\% mobility for most regions in the city. The economic loss due to epidemic control can be significantly reduced. In all, DURLECA successfully resolves the life-or-economy dilemma.

DURLECA's control is highly customized and dynamic, which is hard to be mimicked by human experts. In Figure \ref{fig:quota}, we visualize the spatial distribution of quota rates and the associated histogram in four selected periods. Figure \ref{fig:quota} manifests that DURLECA's smartness in distributing quotas according to both epidemic risks and mobility patterns. The agent tends to give more quotas in a low-risk and low-mobility period and give fewer quotas in either a high-risk or a high-mobility period.

\subsection{Comparison of the Scenarios of  Early/Late Intervention}

To examine whether DURLECA is still effective if the government's intervention is later than the discovery of the first patient, we compare DURLECA's performance in three scenarios.
We have discussed the very-late-intervention scenario where the government starts to act 20 days after discovering the first patient ($t_{start}=20$) in Section 4.3. Here, we compare the early-intervention scenario ($t_{start}=0$) and the late-intervention scenario ($t_{start}=10$). The results are shown in Table \ref{tab:early_late}.

We find that:
1) EP-Soft can control the epidemic in the early-intervention scenario. Because the virus has not widely spread, restricting the few infected areas is enough for epidemic suppressing. However, it fails to avoid an epidemic outbreak in the late-intervention scenario.
2) EP-Hard and EP-Lockdown can control the epidemic under both scenarios. However, it will lock all risky regions down and cut off most mobility.
3) DURLECA successfully suppresses the epidemic while retains the majority of urban mobility in both scenarios.

\begin{table}[ht]
    \centering
    \scalebox{0.95}{
        \begin{tabular}{l|l|lll}
        \hline
        \hline
         & $t_{start}$ & Mean/Max $H$ & Total $R$ & $Q$ \\
        \hline
        EP-Soft    & 0  & 0.02/0.03 & 0.08 & 0.18 \\
        EP-Hard    & 0  & 0.03/0.04 & 0.10 & 0.13 \\
        EP-Lockdown    & 0  & 0.01/0.03 & 0.09 & 0 \\
        DURLECA    & 0  & 0.03/0.04 & 0.27 & 0.73 \\ \hline

        EP-Soft    & 10 &  4.67/62.10 & 1041.64 & 0.57 \\
        EP-Hard    & 10  & 0.08/0.14 & 0.64 & 0.11 \\
         EP-Lockdown    & 10 &  0.08/0.14 & 0.55 & 0 \\
         DURLECA    & 10 &  0.07/0.16 & 1.31 & 0.74 \\
        \hline
        \hline
        \end{tabular}
    }
    \caption{The simulation results of DURLECA and three baselines when $t_{start}=0,10$.}
\label{tab:early_late}
\vspace{-0.6cm}
\end{table}

\subsection{Generalization Ability}
We examine the generalization ability of DURLECA under different urban settings and diseases. Cities have different capacities for hospitalization treatment. Heterogeneous economic structures also cause cities' divergent tolerance for mobility restrictions. We vary the setting of $\{H_0, L_0\}$, which represents the change of urban features, and examine DURLECA's performance. The results are shown in Table~\ref{tab:h0l0}. The results demonstrate that DURLECA can find out different policies responding to the change of urban settings. For instance, we find that a higher $H_0$ leads to more mobility and more hospitalizations, which suggests that cities with higher hospitalization capacities can take more patients and retain more mobility. We also examine DURLECA's adaptiveness to various diseases. We vary the setting of $\{\beta_m$, $\beta_s$, $\gamma$, $\theta\}$ to simulate different diseases with different $R0$.
We find that DURELCA is also able to adjust epidemic-control policy to adapt to different diseases (Table \ref{tab:r0}). For instance, DURLECA provides loose mobility restrictions on low-$R0$ diseases but stringent mobility restrictions on high-$R0$ ones. DURLECA's adaptiveness to urban-setting and disease-setting changes not only demonstrates its generalization ability but also its smartness.

\begin{table}[ht]
    \centering
    \scalebox{1}{
        \begin{tabular}{ll|llll}
        \hline
        \hline
         $H_0$ & $L_0$ & Mean/Max $H$ & Total $R$ & $Q$ \\
        \hline
         1 & 72  & 0.54/1.74 & 10.02 & 0.38 \\
         3 & 72  & 0.60/2.28 & 19.07 & 0.76 \\
         10 & 72  & 2.79/6.34 & 223.29 & 0.90 \\ \hline
         3 & 48  & 1.69/4.60 & 153.54 & 0.88 \\
         3 & 72  & 0.60/2.28 & 19.07 & 0.76 \\
         3 & 168 & 0.45/1.58 & 16.69 & 0.71 \\
        \hline
        \hline
        \end{tabular}
    }
    \caption{The simulation results of DURLECA with different $H_0,L_0$ when $t_{start}=20$.}
\label{tab:h0l0}
\vspace{-0.4cm}
\end{table}

\begin{table}[ht]
    \centering
    \scalebox{1}{
        \begin{tabular}{l|lll}
        \hline
        \hline
         $R_0$ & Mean/Max $H$ & Total $R$ & $Q$ \\
        \hline
         1.4  & 0.35/1.20 & 10.96 & 0.86 \\

         2.1  & 0.60/2.28 & 19.07 & 0.76 \\

         3.5  & 0.81/4.08 & 23.18 & 0.46\\
        \hline
        \hline
        \end{tabular}
    }
    \caption{The simulation results of DURLECA under epidemics with different $R0$.}
\label{tab:r0}
\vspace{-0.8cm}
\end{table}
\section{Related work}
\textbf{Epidemic Modeling:}
The SIR model is a widely used mathematical model in epidemiology, which divides the population into three states: susceptible, infected and recovered~\cite{SIR}. Based on the SIR model, Ogren and Martin used an embedded Newton algorithm to help find an optimal control strategy \cite{ogren2000optimal}. The distributed delay and discrete delay of SIR was also studied \cite{mccluskey2010complete}. Considering a more practical epidemic scenario, the SEIR model added an Exposed state to deal with the incubation period \cite{SEIR}. Others also strengthened the differential equations considering vaccination consequences for a measles epidemic \cite{allen10introduction}.
Later works also tried to incorporate human spatial patterns into the epidemic model. Sattenspiel et al. presented how contacts occur between individuals from different regions and how they influence epidemic spreads~\cite{sattenspiel1995structured}.
Balcan et al. presented the GLobal Epidemic and Mobility model, which integrated sociodemographic and population mobility data in a spatially structured stochastic approach~\cite{balcan2010modeling}.
Different from previous works, we distinguish visible and invisible infections and model epidemic transmission upon traffic flows, so that to support exploring mobility-control policies for epidemic control.

\textbf{Graph Neural Network for OD-flows:}
The problem of estimating, predicting and controlling human flows between regions has been addressed using neural networks since \cite{lorenzo2013od}. Especially due to the reason that most OD flows are modeled based on graphs,
Graph Neural Network (GNN) shows great importance and was first suggested in \cite{scarselli2009the}. Later GNNs were used to predict future mobility flows~\cite{chai2018bike, geng2019spatiotemporal, wang2019origin}. Besides, \cite{epidemic_gcn} borrowed knowledge from epidemic models to design GNN for node prediction in documents.
However, existing GNN architectures lack the ability to model the virus-spreading flow. Our designed Flow-GNN allows our model to characterize the virus-spreading flow and guarantees DURLECA's capability.

\textbf{Deep Reinforcement Learning:}
Deep Reinforcement Learning (DRL) has been proved to be effective for control problems that have a large action space \cite{DQN, DDPG, double_q}. DQN \cite{DQN} and DDPG \cite{DDPG} are two representative DRL algorithms, proposed for discrete control problems and continuous control problems, respectively.
To enable the agent to find an optimal solution, later works proposed to enhance exploration \cite{param_noise, explore_1, explore_2}. Imitation learning is another area of RL, where the goal is to enable the agent to behave like a human expert \cite{imitation}. AlphaGo proposed to start from imitation but further explore to outperform expert \cite{go}. In \cite{graph_protection}, DQN was also used for node protection against epidemic under a single objective. Compared with it, both our control action and objectives are more complex and practical, and thus our RL training are more challenging. We design two strategies to address the exploration challenge.

\section{Future research and Conclusion}

A series of problems ask for future study on smart-and-privacy-protected epidemic control while DURLECA can be the framework. In this research, we do not consider the uncertainty of mobility and epidemic information when DURLECA explores epidemic policies. It asks for future work to explore the algorithm for searching a robust policy when the information is uncertain. A practical policy has to be robust even if there exist errors in the input data.

To conclude, this research demonstrates a sequence of important facts, which broaden the vision of human society for epidemic control and are listed below:
\begin{itemize}
    \item Private data is dispensable because restricting the aggregated inter-regional mobility sufficiently lowers the probability of infectious people's movement and thus suppresses the risk of epidemic transmission.
    \item Resolving the life-or-economy dilemma of epidemic control must allow dynamic and customized regional policies.
    \item The powerfulness of our GNN-enhanced RL in epidemic control manifests that field knowledge is critical for AI-system architecture and valuable for neural network training.
\end{itemize}
In all, smart governance empowered by AI will protect future society from the loss of lives due to epidemics and the economic risk caused by epidemic control.

\section{Notation Summary}

\begin{table}[ht]
\centering
\scalebox{1}{
\begin{tabular}{l|l}
\hline
\hline
Term/Notation & Definition \\
\hline
superscript $\tau$  & At time step $\tau$. \\
subscript $d$  & The original demand without restrictions. \\
subscript $p$  & With restriction $p$. \\
subscript $i$, $j$ & Region index \\ \hline

$M_*^\tau$  & The mobility. A matrix.\\
$M_{*,i,j}^\tau$  & The OD flow from $i$ to $j$. A scalar.  \\
$M_{*,i}^\tau$  & $\sum_j M_{*,i,j}^\tau$. The out-flow from $i$. A scalar.  \\
$\overline{M_{*,i}}$ &  $\frac{1}{T}\sum_\tau M_{*,i}^\tau$. The mean out-flow from $i$. A scalar. \\ \hline
$Q_i^\tau$ & $\frac{M_{p,i}^\tau}{M_{d,i}^\tau}$. The region quota rate. A scalar. \\
$Q^\tau$ & $\frac{\sum_i M_{p,i}^\tau}{\sum_i M_{d,i}^\tau}$. The city quota rate. A scalar. \\
$Q$ & $\frac{\sum_\tau \sum_i M_{p,i}^\tau}{\sum_\tau \sum_i M_{d,i}^\tau}$. The total quota rate. A scalar. \\
\hline
\hline
\end{tabular}
}
\caption{The summary of mobility-related notations. Subscript $*$ can be either $d$ or $p$.}
\label{tab:definitions}
\end{table}

\section{Experiment Settings and Reproducibility}

\subsection{Dataset}
We list the dataset details in Table \ref{tab:dataset}, where $P_{move}$ counts the mean probability for an individual to move in one hour.
\begin{table}[ht]
\centering
\begin{tabular}{lllll}
\hline
\hline
City & Regions & Mean Population & $P_{move}$ & Duration\\
\hline
Beijing & $17\times19$ & 1686 & 0.18 & 744 Days\\
\hline
\hline
\end{tabular}
\caption{The summary of the prolonged dataset.}
\label{tab:dataset}
\end{table}

\textbf{Privacy and ethical concerns:}
We have taken the following procedures to address privacy and ethical concerns. First, all of the researchers have been authorized by the data provider to utilize the data for research purposes only. Second, the data is completely anonymized. Third, we store all the data in a secured off-line server.

\subsection{Implementation Details}
Without the loss of generality, we set the moving transmission rate $\beta^m = \frac{3}{24}$, the staying transmission rate $\beta^s = \frac{0.1}{24}$, the hospitalized rate $\gamma=\frac{0.3}{24}$ and the recover rate $\theta=\frac{0.3}{24}$. Without intervention, the estimated basic reproduction rate $R_0$ is 2.1 at the initial stage of the epidemic.
For the reward, we mainly set $\lambda = 0.99, L_0 = 72, H_0 = 3, k_h = 1$. For the pseudo expert, we set $X_{h} = 1, X_{l} = 168$. For the extreme point threshold, we set $I_{t} = 100, L_{t} = 336$.

During training, we randomly initialize an epidemic state at the start of each episode. We train DURLECA for 400,000 steps, using Adam optimizer with the learning rate as 0.0001.
During testing, we fix one epidemic-initialization setting and compare different baselines. Considering the randomness of training, we train DURLECA with different random seeds 5 times for each set of configurations, and choose the one that achieves the best episode reward to report as the result.

We mainly implement DURLECA based on Keras-RL \cite{keras_rl} with our modifications.

\subsection{Disease $R0$}
In a classical SIR model, the basic reproduction rate $R0$ is calculated as $R0=\frac{\beta}{\gamma}$ \cite{R0}. In our model, as the infection has been divided into two parts, we estimate an averaged $\overline{\beta}$ over $\beta_m$ and $\beta_s$ according to their corresponding population size,
\begin{align}
    \overline{\beta} = P_{move}\beta_m + (1-P_{move})\beta_s.
\end{align}
Then we estimate $R0=\frac{\overline{\beta}}{\gamma}$.

In Section 4.5, we vary the setting of $\{\beta_m, \beta_s\}$ and keep $\{\gamma, \theta\}$ the same. To make a fair comparison, we also vary $t_{start}$ for each simulated disease to make sure the city has nearly the same number of hospitalized people when we start the intervention. For $R0=1.4$, we set $\{\beta_m  = 1.9, \beta_s = 0.1, t_{start} = 45\}$. For $R0=2.1$, we set $\{\beta_m  = 3, \beta_s = 0.1, t_{start} = 20\}$. For $R0=3.5$, we set $\{\beta_m  = 5, \beta_s = 0.2, t_{start} = 10\}$.

\section{Ablation Study}

To evaluate the effectiveness of our proposed Flow-GNN and RL exploration strategies, we conduct ablation studies in this section.

\textbf{GNN-Baselines}: To evaluate the effectiveness of our proposed Flow-GNN, we use the well known GraphSageConv layer \cite{graphsage} and our modified GraphSageConv layer to replace the proposed Flow-GNN layer in the actor network and the critic network. We name the two baselines as \textbf{GNN-Mean} and \textbf{GNN-Softmax}.

The layer calculation of GNN-Mean follows Equation~(\ref{eq:GNN-mean}):
\begin{align}
    & f_i^k = \sigma (W^k (f_{i}^{k-1}, \frac{1}{N(i)}\sum_{j\in N(i)} f_{j}^{k-1}) + B^k), \label{eq:GNN-mean}
\end{align}
where $N(i)$ denotes the connected regions of $i$.

The layer calculation of GNN-Softmax follows Equation~(\ref{eq:GNN-softmax}):
\begin{align}
    & w_{j,i} = \frac{exp(M_{*,j,i}^{\tau+k-1})}{\sum_{j\in N(i)}exp(M_{*,j,i}^{\tau+k-1})}, \notag \\
    & f_i^k = \sigma (W^k (f_{i}^{k-1}, \sum_{j\in N(i)} w_{j,i} f_{j}^{k-1}) + B^k). \label{eq:GNN-softmax}
\end{align}

\textbf{RL-Baselines}: To evaluate the effectiveness of our RL exploration strategies, we remove the pseudo-expert strategy and the avoiding-extreme-points strategy, respectively. We refer to the two baselines as \textbf{RL-NoEP} and \textbf{RL-NoThre}.

\begin{table}[ht]
    \centering
    \scalebox{0.9}{
        \begin{tabular}{l|lllll}
        \hline
        \hline
          & Mean/Max $H$ & Total $R$ & $Q$ & $T_{20\%}^c$  & $T_{20\%}^r$ \\
        \hline

        No Intervention  & 27.21/157.55 & 1069.29 & 1 & 0 & 0 \\ \hline

        GNN-Mean & -/- & - & - & - & - \\
        GNN-Softmax  & 0.53/1.88 & 7.79 & 0.06 & 26 & 28 \\

        RL-NoEP  & 0.41/1.45 & 5.87 & 0.00 & 27 & 27 \\
        RL-NoThre  & 1.43/3.68 & 86.82 & 0.75 & 0 & 9 \\ \hline

        DURLECA & 0.60/2.28 & 19.07 & 0.76 & 0 & 0 \\
        \hline
        \hline
        \end{tabular}
    }
    \caption{Ablation study when $t_{start}=20$.}
\label{tab:ablation}
\end{table}

\textbf{Results and Analysis:}
As shown in Table \ref{tab:ablation}, without Flow-GNN or the proposed RL exploration strategies, the agent fails to learn a good policy.

The failure of GNN-Mean comes from its inability to learn weighted edge information, i.e., how many people move from one region to another. With considering weighted edge information, GNN-Softmax still fails to retain mobility, as it can not describe traffic flows and the epidemic transmission upon it. These prove the effectiveness of our proposed Flow-GNN.

RL-NoEP gives a long-term lockdown to the whole city, which is a typical local optimum. As for RL-NoThre, the agent successfully finds one policy that achieves relatively low hospitalizations and high mobility. However, this solution is worse than DURLECA. Besides, we find that the success of RL-NoThre highly relies on luck. During our five repeating experiments, the agent was stuck in local optimums for four times, giving a long-term lockdown to the whole city. These, as discussed earlier in Section 1, are due to the difficulty of exploration. The agent is easy to encounter extreme points during exploration, and the extreme points force the agent to adopt conservative policies, i.e., lock the whole city down. Compared with RL-NoEP and RL-NoThre, DURLECA is guided by a pseudo-expert and is designed to avoid extreme points. Thus, DURLECA can find much better solutions.