\section{Introduction}
Neuroimaging, conventional machine learning, and deep learning have offered unprecedented opportunities to understand the underlying mechanism of brain disorders over the past decades (\cite{ml_neuro}), and pave the road towards individualized precision medicine (\cite{precision}). A large body of case-control studies leverage mass univariate group comparisons to derive structural or functional neuroimaging signatures (\cite{biomarker_ad1}, \cite{biomarker_ad2}, \cite{biomarker_ad3}). However, these studies suffer from under-powered statistical inferences since the heterogeneous nature of neurological diseases usually violate the assumption that each group population are pathologically homogeneous.\\
A developing body of methodologies aim to disentangle this heterogeneity with various machine learning methods. \cite{lda_brain} used the Bayesian Latent Dirichlet Allocation (LDA) model to identify latent atrophy patterns from voxel-wise gray matter (GM) density maps. However, due to the nature of the LDA model, the density maps need to be first discretized, and the model hypothesized to exclusively capture brain atrophy while ignoring the potential regional enlargement. \cite{sustain} proposed another method, Sustain, to uncover temporal and phenotypic heterogeneity by inferring both subtypes and stages. However, the model only handles around 10-20 large brain regions without more detailed information. Notably, both methods applied directly in the PT domain, confront main limitations in avoiding potential disease-unrelated brain variations. In contrast, semi-supervised clustering methods (\cite{HYDRA}, \cite{Chimera}, \cite{MAGIC}, \cite{SmileGAN}) were proposed to cluster patients via the patterns or transformations between the reference group (CN) and the target patient group (PT). A recent proposed deep learning-based semi-supervised clustering method, termed Smile-GAN (\cite{SmileGAN}), achieved better clustering performance by learning multiple mappings from CN to PT with an inverse clustering function for both regularization and clustering. Albeit these methods demonstrated potential in clinical applications, they have a major limitation. In particular, these methods model disease heterogeneity as a dichotomous process and derive a discrete output, i.e., subtype, ignoring the fact that brain disorders progress along a continuum. Moreover, variations among subtypes were contributed by both spatial and temporal differences in disease patterns, thus confounding further clinical analysis.(Fig. \ref{fig1}a) \\
To address the aforementioned limitations, we propose a novel method, Surreal-GAN (Semi-Supervised Representation Learning via GAN), for deriving heterogeneity-refined imaging signatures. Surreal-GAN is a significant extension of Smile-GAN under the same principle of semi-supervised learning by transforming data from CN domain ${\mathbb{X}}$ to PT domain ${\mathbb{Y}}$. Herein, the key extension is to represent complex disease-related heterogeneity with low dimensional representations with each dimension indicating the severity of one relatively homogeneous imaging patterns (Fig. \ref{fig1}b). We refer to these low dimensional representations as R-indices (${\mathbf{r}}_i$, where $i$ represent the $i_{th}$ dimension) of disease pattern. The first key point of the method is modelling disease as a continuous process and learning infinite transformation directions from CN to PT, with each direction capturing a specific combination of patterns and severity. The idea is realized by learning one transformation function which takes both normal data and a continuous latent variable as inputs and output synthesized-PT data whose distribution is indistinguishable from that of real PT data. The second key point is controlling the monotonicity (from light to severe with increasing index value) of disease patterns through one monotonicity loss and double sampling procedure. The third key point is boosting the separation of captured disease patterns through one orthogonality regularization. The fourth key point is introducing an inverse mapping function consisting of a decomposition part and a reconstruction part which further ensure that patterns synthesized through transformation are exactly captured for inferring representations of disease patterns, the R-indices. Besides the above mentioned regularization, function continuity, transformation sparsity and inverse consistency further guide the model to capture meaningful imaging patterns.\\
To support our claims, we first validated Surreal-GAN on semi-synthetic data with known ground truth of disease patterns and severity. We compared performance of Surreal-GAN to NMF (\cite{nmf}), opNMF (\cite{OPNMF}), Discriminant-NMF (\cite{DNMF}), Factor Analysis and LDA (\cite{lda}) models, and subsequently studied the robustness of model under various conditions. Finally, we applied the model to one real dataset for Alzheimer's disease (AD), and demonstrated significant correlations between the proposed imaging signatures and clinical variables.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\linewidth]{figure/figure1.png}
\end{center}
\caption{Model Ideas and Schematic Diagram. (a) Results from Semi-supervised clustering for AD disease: model clusters PT into discrete subtypes based on both spatial and temporal variations: P1 and P4 show mild and severe stages and are mixture of two patterns; (b) Improvements from Surreal-GAN: two patterns are captured by two indices with values indicating severity; (c) General architecture of networks; (d) Schematic diagram of Surreal-GAN.}
\label{fig1}
\end{figure}

\section{Methods}
The schematic diagram of Surreal-GAN is shown in Fig. \ref{fig1}d. The model is applied on regional volume data (volumes of different graymatter (GM) and whitematter (WM) regions) derived from structural MRI. The essential idea is to learn one transformation function $f:{\mathbb{X}}*{\mathbb{Z}} \rightarrow {\mathbb{Y}}$ which transforms CN data ${\mathbf{x}}$ to different synthesized PT data ${\mathbf{y}}'=f({\mathbf{x}},{\mathbf{z}})$, with latent variable ${\mathbf{z}}$ specifying distinct combinations of patterns and severity. Here, the Latent (LAT) domain, ${\mathbb{Z}}=\{{\mathbf{z}}:0\leq {\mathbf{z}}_i\leq 1\ \forall 1\leq i\leq M\}$, is a class of vectors with dimension $M$ (predefined number of patterns). We denote four different data distributions as ${\mathbf{x}}\sim p_{cn}({\mathbf{x}})$, ${\mathbf{y}}\sim p_{pt}({\mathbf{y}})$, ${\mathbf{y}}'\sim p_{syn}({\mathbf{y}}')$ and ${\mathbf{z}}\sim p_{lat}({\mathbf{z}})$, respectively, where ${\mathbf{z}}\sim p_{lat}({\mathbf{z}})$ is sampled from a multivariate uniform distribution $U[0,1]^M$, rather than a categorical distribution as in Smile-GAN. A continuous z variable lay the foundation for learning continuous representations. In addition to the transformation function, an adversarial discriminator $D$ is introduced to distinguish between real PT data ${\mathbf{y}}$ and synthesized PT data ${\mathbf{y}}'$, thereby ensuring that the synthesized PT data are indistinguishable from real PT data (i.e. minimizing the distance between real PT data distribution and synthesized PT data distribution).\\
A continuous latent variable does not necessarily lead to desired properties of CN to PT transformation (e.g. separation/monotonicity of captured patterns). That is, there are a number of functions potentially achieving equality in distributions, thus making it hard to guarantee that the transformation function learned by the model is closely related to the underlying pathology progression process. Moreover, during the training procedure, the transformation function backboned by the neural network tends to ignore the latent variable ${\mathbf{z}}$. Even when the Latent variable ${\mathbf{z}}$ is emphasized through regularization, different components of ${\mathbf{z}}$ variable do not naturally lead to separate disease patterns and are not positively correlated with pattern severity. Therefore, with the assumption that there is one true underlying function for real PT data ${\mathbf{y}}=h({\mathbf{x}},{\mathbf{z}})$, Surreal-GAN aims to boost the transformation function $f$ to approximate the true underlying function $h$, while having the latent variable, ${\mathbf{z}}$, contributing to recognizable imaging patterns rather than being random "activators" of the GAN. For this purpose, we constrain the function class of $f$ via different types of regularization as follows: 1) encourage sparse transformations; 2) enforce Lipschitz continuity of functions; 3) introduce one inverse function $g:{\mathbb{Y}}\rightarrow{\mathbb{Z}}$ serving to decompose and reconstruct; 4) boost orthogonality among synthesized/captured patterns; 5) enforce monotonicity of synthesized patterns and their positive correlations with components of ${\mathbf{z}}$. The training procedure of Surreal-GAN is thus a adversarial learning process constrained by these regularization.
\subsection{Adversarial Loss}
The adversarial loss (\cite{GAN}) is used for iterative training of discriminator $D$ and the transformation function $f$, which can be written as:
\begin{align}
L_{\text{GAN}}(D,f)&=\mathbb{E}_{y\sim p_{pt}({\mathbf{y}})}[\log(D({\mathbf{y}}))]+\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim P_{lat}({\mathbf{z}})}[\log(1-D(f({\mathbf{x}},{\mathbf{z}})))]\\
&=\mathbb{E}_{{\mathbf{y}}\sim p_{pt}({\mathbf{y}})}[\log(D({\mathbf{y}}))]+\mathbb{E}_{{\mathbf{y}}'\sim p_{syn}({\mathbf{y}}')}[\log(1-D({\mathbf{y}}'))]
\end{align}
The transformation function $f$ attempts to transform CN data to synthesized PT data so that they follow similar distributions as real PT data. The discriminator, providing a probability that ${\mathbf{y}}$ comes from the real data rather than transformation function, tries to identify the synthesized PT data and distinguish it from real PT data. Therefore, the discriminator $D$ attempts to maximize the adversarial loss function while the transformation function $f$ attempts to minimize against it. The training process can be denoted as:
\begin{align}
    \min_f\max_DL_{\text{GAN}}(D,f)=\mathbb{E}_{{\mathbf{y}}\sim p_{pt}({\mathbf{y}})}[\log(D({\mathbf{y}}))]+\mathbb{E}_{{\mathbf{y}}'\sim p_{syn}({\mathbf{y}}')}[\log(1-D({\mathbf{y}}'))]
\end{align}
\subsection{Regularization}
\subsubsection{Sparse Transformation}
We assumed that disease process will not change brain anatomy dramatically and primarily only affect certain regions throughout most of disease stages, which means that the true underlying transformation only moderately changes some regions while keeping the rest mildly affected with most of ${\mathbf{z}}$ variable values. Therefore, to control sparsity and distance of transformations, we defined one change loss to be the $l_1$ distance between the synthesized PT data and the original CN data:
\begin{align}
    L_{\text{change}}(f)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim p_{lat}({\mathbf{z}})}[||f({\mathbf{x}},{\mathbf{z}})-{\mathbf{x}}||_1]
\end{align}

\subsubsection{Lipschitiz Continuity and Inverse Mapping}
First, by enforcing the transformation function $f$ to be $K_1$-Lipschitz continuous, we can derive that, for fixed latent variable ${\mathbf{z}}={\bm{a}}$ and $\forall {\bm{x}}_1,{\bm{x}}_2\in {\mathbb{X}}$, $||f({\bm{x}}_1,{\bm{a}})-f({\bm{x}}_2,{\bm{a}})||_2\leq K_1||{\bm{x}}_1-{\bm{x}}_2||_2$. Thus, by controlling the constant $K_1$, we constrained the transformation function to preserve the original distances among CN data instead of scattering them into the PT domain if a same combination of disease patterns is imposed to them. \\
Second, to prevent latent variable ${\mathbf{z}}$ from being ignored, we introduced an inverse mapping function from PT domain to Latent domain, $g:{\mathbb{Y}}\rightarrow{\mathbb{Z}}$. By Lemma1 (section A.1), with function $g$ being $K_2$-Lipschitz continuous and $d(\cdot,\cdot)$ being any distance metric satisfying triangular inequality, we can derive that, $\forall {\bm{z}}_1,{\bm{z}}_2\sim p_{lat}({\mathbf{z}}), {\bm{z}}_1\neq {\bm{z}}_2$, and $\bar{{\bm{x}}}\sim p_{cn}({\mathbf{x}})$,  $d(f(\bar{{\bm{x}}},{\bm{z}}_1),f(\bar{{\bm{x}}},{\bm{z}}_2))$ is lower bounded by $(\frac{d({\bm{z}}_1,{\bm{z}}_2)}{K_2}-\frac{1}{K_2}(d(g(f(\bar{{\bm{x}}},{\bm{z}}_1)),{\bm{z}}_1)+d(g(f(\bar{{\bm{x}}},{\bm{z}}_2)),{\bm{z}}_2)))$. Therefore, by minimizing the distance between sampled latent variable ${\mathbf{z}}$ and reconstructed latent variable $g(f({\mathbf{x}},{\mathbf{z}}))$, we can control differences between synthesized patterns to be positively related to distances between latent variables, and thus to be non-trivial (i.e., same CN data is transformed to significantly different PT data with different ${\mathbf{z}}$ variables). Therefore, we define a reconstruction loss as the $l_2$ distance between $g(f({\mathbf{x}},{\mathbf{z}}))$ and ${\mathbf{z}}$:
\begin{align}
    L_{\text{recons}}(f,g)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim p_{lat}({\mathbf{z}})}[||g(f({\mathbf{x}},{\mathbf{z}}))-{\mathbf{z}}||_2]
\end{align}

\subsubsection{Pattern Decomposition}
Nevertheless, a simple function $g$ directly reconstructing sampled ${\mathbf{z}}$ from synthesized PT ${\mathbf{y}}'$ do not reconstruct each component ${\mathbf{z}}_i$ completely based on changes led by ${\mathbf{z}}_i$ in transformation. We are more interested in pattern representations of real PT data than generating synthesized PT, and inverse function $g$ will be used for inferring R-indices as introduced in section 2.4. Thus, it is very important for the inverse function to accurately capture patterns synthesized by intensively regularized transformation function $f$. For this purpose, we further separated the function $g:{\mathbb{Y}}\rightarrow {\mathbb{Z}}$ into $g_1:{\mathbb{Y}}\rightarrow \mathbb{R}^{M*S}$ (where $M$ is number of patterns and $S$ is number of input ROIs) and $g_2:\mathbb{R}^S\rightarrow \mathbb{R}$. Decomposer $g_1$ serves to reconstruct changes synthesized by each component $z_i$ in transformation process: ${\mathbf{q}}_i=f({\mathbf{x}},{\mathbf{a}}^i)-{\mathbf{x}}$, where ${\mathbf{a}}^i$ is a vector with $i_{th}$ component ${\mathbf{a}}^i_i={\mathbf{z}}_i$ and ${\mathbf{a}}^i_j=0$, $\forall i\neq j$. Let $\hat{{\mathbf{q}}}_{f(x,z)}=[{\mathbf{q}}_1^T,\cdots,{\mathbf{q}}_M^T]^T$ be concatenation of all synthesized changes ${\mathbf{q}}_i$, we defined the decomposition loss as:
\begin{align}
    L_{\text{decom}}(f,g_1)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim p_{lat}({\mathbf{z}})}[||g_1(f({\mathbf{x}},{\mathbf{z}}))-\hat{{\mathbf{q}}}_{f({\mathbf{x}},{\mathbf{z}})}||_2]
\end{align}
$g_2$ serves to further reconstruct each component of sampled ${\mathbf{z}}$ variable from $g_1(f({\mathbf{x}},{\mathbf{z}}))$  independently, so that the reconstruction loss can be rewritten as:
\begin{align}
    L_{\text{recons}}(f,g) = L_{\text{recons}}(f,g_1,g_2)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim p_{lat}({\mathbf{z}})}[||\hat{{\mathbf{l}}}_{g_2(g_1(f({\mathbf{x}},{\mathbf{z}}))}-{\mathbf{z}}||_2]
\end{align}
where $\hat{{\mathbf{l}}}_{g_2(g_1(f({\mathbf{x}},{\mathbf{z}}))}=g(f({\mathbf{x}},{\mathbf{z}}))=[g_2(g_1(f({\mathbf{x}},{\mathbf{z}}))_{0:S}),\cdots,g_2(g_1(f({\mathbf{x}},{\mathbf{z}}))_{S*(M-1):S*M})]^T$.\\ Function $g$ here can be also considered as the approximation of expectation of posterior distribution. In this sense, minimization of $l_2$ distances can also be interpreted as maximization of mutual information between latent variable ${\mathbf{z}}$ and synthesized PT data ${\mathbf{y}}'=f({\mathbf{x}},{\mathbf{z}})$ as explained in Remark1 in section A.2 (\cite{infogan}). From this perspective, the transformation function is forced to best utilize information from latent variable ${\mathbf{z}}$, while also keeping mutual information between synthesized PT data ${\mathbf{y}}'$ and original CN data ${\mathbf{x}}$ by controlling transformation distances.

\subsubsection{Orthogonality of Patterns}
Inclusion of latent variable ${\mathbf{z}}$ in the transformation function does not guarantee that different components of the latent variable, ${\mathbf{z}}_i$, contributes to relatively different patterns in synthesized PT data. Instead, during training procedure, different components, ${\mathbf{z}}_i$, are more likely to synthesize patterns in same regions and lead to accumulated severity. Therefore, we define another orthogonal loss to encourage separation among patterns synthesized by different components. With changes led by each component, ${\mathbf{q}}_i$ defined in section 2.2.3, we constructed a matrix ${\mathbf{A}}_{f({\mathbf{x}},{\mathbf{z}})}$ with the $i_{th}$ column ${\mathbf{A}}_{{f({\mathbf{x}},{\mathbf{z}})}_{:,i}}=|{\mathbf{q}}_i|/||{\mathbf{q}}_i||_2$. To encourage separation of changes led by different components, we boosted the matrix ${\mathbf{A}}_{f({\mathbf{x}},{\mathbf{z}})}$ to be relatively orthogonal by minimizing the following orthogonality loss:
\begin{align}
    L_{\text{ortho}}(f)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim p_{lat}({\mathbf{z}})}[||{\mathbf{A}}_{f({\mathbf{x}},{\mathbf{z}})}^T{\mathbf{A}}_{f({\mathbf{x}},{\mathbf{z}})}-{\bm{I}}||_F]
\end{align}
\subsubsection{Monotonicity and Positive Correlation}
Besides separation of synthesized patterns, a positive correlation between values of each component ${\mathbf{z}}_i$ (from 0 to 1) and severity of synthesized patterns can not be acquired spontaneously either. Regarding severity of the disease pattern, we assumed that, as a pattern is becoming more severe, absolute values of changes in regional volumes can only increase monotonically or remain constant. In another word, there is no oscillation in regional volumes as disease pattern is becoming more severe, while more regions can be included gradually with absolute changes switched from 0 to a positive value. To satisfy this requirement, we sampled another latent variable ${\mathbf{z}}'\sim p_{sev}({\mathbf{z}}'|{\mathbf{z}})$, conditioned on previous sampled ${\mathbf{z}}$ variable, such that ${\mathbf{z}}'_i\geq {\mathbf{z}}_i$, $\forall 1\leq i\leq M$. With these double sampled latent variables, we defined the monotonicity loss to be:
\begin{align}
    L_{\text{mono}}(f)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}\sim p_{lat}({\mathbf{z}}),{\mathbf{z}}'\sim p_{sev}({\mathbf{z}}'|{\mathbf{z}})}[||\max(|f({\mathbf{x}},{\mathbf{z}})-{\mathbf{x}}|-|f({\mathbf{x}},{\mathbf{z}}')-{\mathbf{x}}|,\mathbf{0})||_2]
\end{align}
Minimization of this monotonicity loss penalize the case that volume changes induced by ${\mathbf{z}}$ is larger than that induced by ${\mathbf{z}}'$, while the other direction is permitted. To further ensure a positive correlation (small ${\mathbf{z}}$ lead to mild patterns; large ${\mathbf{z}}$ lead to severe patterns), we introduced a cn loss, which serves to more strictly constrain the distance between synthesized PT data ${\mathbf{y}}'$ and CN data ${\mathbf{x}}$, when latent variable ${\mathbf{z}}$ is close to $\mathbf{0}$ vector. By letting $p_{cn}({\mathbf{z}})=U(0,0.05)^M$ to be a multivariate uniform distribution, the cn loss is defined as:
\begin{align}
    L_{\text{cn}}(f)=\mathbb{E}_{{\mathbf{x}}\sim p_{cn}({\mathbf{x}}),{\mathbf{z}}^{cn}\sim p_{cn}({\mathbf{z}})}[||f({\mathbf{x}},{\mathbf{z}}^{cn})-{\mathbf{x}}||_1]
\end{align}

\subsection{Full Objective}
With the adversarial loss and all other loss functions defined for regularization, we can write the full objective as:
\begin{align}
    L(D,f,g_1,g_2)&=L_{\text{GAN}}(D,f)+\gamma L_{\text{change}}(f)+\kappa L_{\text{decom}}(f,g_1)+\zeta L_{\text{recon}}(f,g_1,g_2)\notag\\&+ \lambda L_{\text{ortho}}(f) +\mu L_{\text{mono}}(f) + \eta L_{\text{cn}}(f)
\end{align}
with $\gamma$, $\kappa$, $\zeta$, $\lambda$, $\mu$ and $\eta$ be parameters controlling the relative importance of different loss functions. Through training process, we want to derive parametrized functions $f$,$g_1$ and $g_2$ satisfying:
\begin{align}
    f,g_1,g_2 = \arg\min_{f,g}\max_{D}L(D,f,g_1,g_2)
\end{align}

\subsection{Representation-Indices Inference}
With the assumption that $p_{syn}(f({\mathbf{x}},{\mathbf{z}}))\approx p_{pt}({\mathbf{y}})$ and $f$ satisfying all constraints, we consider learned transformation function $f$ to be a good approximation of the true underlying function $h$, such that $f({\mathbf{x}},{\mathbf{z}})\approx h({\mathbf{x}},\sigma({\mathbf{z}}))$, where $\sigma\in\Omega$ and $\Omega$ is a class of permutation functions which change orders of elements in vector ${\mathbf{z}}$. Since the order of indices in derived representation is not important and we can always reorder them to find the best matching, we simply rewrite it as $f({\mathbf{x}},{\mathbf{z}})\approx h({\mathbf{x}},{\mathbf{z}})$ without loss of generality. For any real PT data, $\bar{{\mathbf{y}}} = h(\bar{{\mathbf{x}}},{\mathbf{r}})\sim p_{pt}({\mathbf{y}})$, we can derive that the ground truth R-indices ${\mathbf{r}}\approx g(f(\bar{{\mathbf{x}}},{\mathbf{r}}))\approx g(h(\bar{{\mathbf{x}}},{\mathbf{r}}))\approx g(\bar{{\mathbf{y}}})$. Therefore, reconstruction function $g$ can be directly applied to infer R-indices of disease patterns, ${\mathbf{r}}$, of any unseen PT data.
\section{Implementation Details}
\subsection{Model Architecture}
The general architecture of networks can be understood from Fig. \ref{fig1}c. The transformation function $f$ utilizes a encoding-decoding structure. Latent variable ${\mathbf{z}}$ and CN data ${\mathbf{x}}$ are first mapped to two vectors with dimension 34 and the element-wise multiplication of them are then decoded to construct the synthesized PT data with dimension 139. The inverse function $g$ first decodes real/synthesized PT data ${\mathbf{y}}/{\mathbf{y}}'$ to larger vector with dimension 139*M and then maps them into a vector with dimension M. The discriminator $D$ utilizes the encoding structure which maps $y/y'$ to a vector with dimension 2. More details of model architectures can be found in Appendix Table \ref{tab3} and \ref{tab4}.
\subsection{Lipschitz Continuity}
In implementation, we performed weight clipping to ensure Lipschitz continuity of transformation function $f$ and inverse function $g$ (\cite{WGAN}). With $\Theta_f$ and $\Theta_g$ denoting the weight spaces of function $f$ and $g$, the compactness of them imply the Lispchitz continuity of two functions, with Lipschitz constants $K_1$ and $K_2$ only depends on $\Theta_f$ and $\Theta_g$. The compactness of $\Theta_f$ and $\Theta_g$ was guaranteed by clapping weight spaces into two closed box, $\Theta_f = [-c_f,c_f]^d$ and $\Theta_g = [-c_g,c_g]^d$. In implementation, both $c_f$ and $c_g$ were set to be 0.5 and these bounds can be further relaxed.  Weight clipping was claimed not to be the most satisfactory method to ensure Lipschitz continuity because of difficulty in choosing the appropriate clipping bounds. However, in our case, we performed weight clipping for function $f$ and $g$ rather than discriminator $D$ for a different purpose and weight clipping contributed to good performances for deriving representations.

\subsection{Training Details}
Regarding optimization procedure, ADAM optimizer(\cite{ADAM}) was used with a learning rate (lr) $4*10^{-5}$ for Discriminator and $2*10^{-4}$ for transformation function $f$ and clustering function $g$. $\beta_1$ and $\beta_2$ are set to be 0.5 and 0.999 respectively. For hyper-parameters, we set $\gamma=6$, $\kappa=80$, $\zeta=80$, $\mu=500$, $\eta=6$. However, performance of model on different tasks are robust to varying hyper-parameters as shown in section B.2. Parameter $\lambda$ determines the relative importance of orthogonality loss and controls the degree of overlapping among synthesized patterns. Without prior knowledge on overlapping degree of underlying disease patterns, we set $\lambda$ to different values in the following experiments and used agreements among multiple trained models to determine the optimal one. Moreover, for all experiments, the batch size was set to be 1/8 of the PT data sample sizes. The model was trained for at least 100000 iterations and saved until the reconstruction loss and the monotonicity loss are smaller than 0.003 and $6*10^{-4}$ respectively. Higher lr might improve convergence speed in some cases. Detailed training procedure is revealed by Appendix Algorithm1.
\section{Experiments}
\subsection{Data Preprocessing}
For semi-synthetic and real data experiments, we used T1-weighted (T1w) MR imaging (MRI) and cognition from the iSAGING (Imaging-based coordinate SysTem for AGIng and NeurodeGenerative diseases) consortium for cognitive impairment and dementia. All baseline MRIs were corrected for intensity inhomogeneities (\cite{data_preprocess1}). Then, 139 regions of interest (ROIs) of brain volumes were extracted using a multi‐atlas label fusion method (\cite{MUSE}), harmonized (\cite{Harmonization}) to remove site effect, and then utilized as input features for Surreal-GAN. (Fig. \ref{supple_fig3})
\subsection{Semi-synthetic experiments}
\textbf{Semi-synthetic data construction} To set the ground truth of disease patterns and severity a prior, we generated semi-synthetic data by imposing simulated disease-related patterns to real CN data, thereby retaining the variation from the real CN data. Specifically, we split 1392 cognitively normal (age $<$ 70 years) into one real CN group with 492 subjects and one Pseudo-PT group with 900 subjects. For the $i_{th}$ Pseudo-PT subject, we sampled a three-dimensional vector as pattern severity from a multivariate uniform distribution, ${\textnormal{s}}_i\sim U[0,1]^3$. Three types of atrophy patterns were introduced to 900 Pseudo-PT subjects. Different regions were included in different patterns as shown in Appendix Table \ref{tab5}. Five different datasets were constructed for testing the robustness of model under different conditions: (1) \textbf{Basic dataset}: For the $i_{th}$ Pseudo-PT subject and the $j_{th}$ ROI included in pattern $k$, we decreased the volume $v_{ij}$ depending on sampled severity: $v_{ij} = v_{ij}-v_{ij}*s_{ik}*N(1,0.05)*0.3$. Each pattern consists of changes in 14 ROIs and has four of them shared with other patterns; (2) \textbf{Large overlapping dataset}: Same as the Basic dataset, except that each pattern shares 8 ROIs with others; (3)
\textbf{Scarce regions dataset}: Each pattern consists of changes in only 4 ROIs with other parts the same as the Basic dataset; (4)\textbf{Within-pattern noise dataset}: Regional changes incorporated in the same pattern have larger within-pattern variances with $v_{ij} = v_{ij}-v_{ij}*s_{ik}*N(1,\mathbf{0.2})*0.3$; (5) \textbf{Mild atrophy dataset}: Compared to the basic case, smaller changes were imposed by letting $v_{ij} = v_{ij}-v_{ij}*s_{ik}*N(1,0.05)*\mathbf{0.2}$.\\
\textbf{Evaluation and Agreement Metric} In our experiments, Concordance Index (c-index) was used as the evaluation metric. Specifically, for each pattern, we calculated a c-index between the inferred R-index and the ground truth. The average of M derived C-indices, refereed as pattern-c-index, was used for model evaluations. Moreover, to quantify the pair-wise agreement between two trained models, we proposed another Pattern-agr-index which equals pattern-c-index between R-indices derived by two different models.
\\
\textbf{Experiment Procedure}
We first validated the relationship between model agreements and their accuracy. On each dataset, we ran the model 10 times with parameter $\lambda=0.1,0.2,0.4,0.6,0.8$ respectively. With each $\lambda$, the average pair-wise pattern-agr-index among 10 models was compared with their average pattern-c-indices. Further, the $\lambda$ value leading to the best agreement was considered as the best choice, and results from 10 corresponding models were reported as model performances for each task. Moreover, we fixed the optimal $\lambda$ value for each task and changed other parameters to 0 to understand contributions from each regularization term (section B.3). Lastly, we compared our model with five other models. NMF and Factor Analysis (FA) are well-known for deriving lower-dimensional representation of complex dataset. The variant opNMF (\cite{OPNMF}) has shown promise in parsing complex brain imaging data with extra orthogonal constraint. Discriminant-NMF (DNMF) (\cite{DNMF}) was proposed to learn representations which best classify subgroups. LDA, though known for extracting topics from documents, have been recently applied to uncover heterogeneity of AD disease from neuroimaging data (\cite{lda_brain}). Data preprocessing for these models were introduced in section B.4.

\subsection{Real Data Experiments}
\textbf{Data Selection}
We applied Surreal-GAN to an AD dataset. For AD, we defined the CN group (N=850) to be subjects with Mini-mental state examination (MMSE) scores above 29, and the PT group (N=2204) as subjects diagnosed as mild cognitive impairment (MCI) or AD at baseline.\\
\textbf{Experiments Procedure}
All CN and PT subjects were first residualized to rule out the covariate effects estimated in the CN group using linear regression. Then, adjusted features were standardized with respect to CN group. Without ground truth, we selected both the optimal number of patterns, $M$, and $\lambda$ parameter by measuring agreements among repetitively trained models. For each combination of $M=2,3,4$ and $\lambda=0.1,0.2,0.4,0.6,0.8$, we ran the model 10 times respectively. The $M$ and $\lambda$ values leading to the highest agreement was considered optimal. Among the 10 corresponding models, the one having the highest mean pair-wise agreements with the other models was used to derive R-indices for all PT subjects. To visualize patterns corresponding to the $i_{th}$ dimension of R-indices, ${\mathbf{r}}_i$, we grouped subjects into three different subgroups: ${\mathbf{r}}_i<0.4$, $0.4<{\mathbf{r}}_i<0.7$ and ${\mathbf{r}}_i>0.7$, with ${\mathbf{r}}_j<0.4$ for all $j\neq i$. Voxel-wise group comparisons were performed between each subgroup and the CN group via AFNI-3dttest (\cite{afni}) and GM-tissue map (\cite{RAVEN}) which encodes the volumetric changes in GM observed during the registration. Pattern representations, R-indices, of PT subjects were further compared with other clinical variables to test their significance. Details of clinical variables and statistical tests can be found in section B.9.

\section{Results}
\subsection{Results on Semi-synthetic Data}
\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{figure/figure2.png}
\end{center}
\caption{Results on Semi-synthetic data (results presented as mean value ± standard error). (a) Representation accuracies on different tasks (measured by Pattern-c-index); (b) Agreements among repetitively trained models on different tasks (measured by Pattern-agr-index).}
\label{fig2}
\end{figure}
\begin{table}[t]
\caption{Model Comparison Results}\label{tab1}
\begin{center}
\begin{tabular}{l@{\hspace{0.02cm}}l@{\hspace{0.02cm}}l@{\hspace{0.02cm}}l@{\hspace{0.02cm}}l@{\hspace{0.02cm}}l@{\hspace{0.02cm}}l}
\multicolumn{1}{c}{}&\multicolumn{1}{c}{\bf Surreal-GAN} &\multicolumn{1}{c}{\bf NMF} &\multicolumn{1}{c}{\bf LDA}&\multicolumn{1}{c}{\bf FA}&\multicolumn{1}{c}{\bf OpNMF}&\multicolumn{1}{c}{\bf DNMF}\\
\hline
\begin{tabular}{c}Basic Case \end{tabular} & \begin{tabular}{c}\bf{0.838}\\ $\pm$ \bf{0.011}\end{tabular}& \begin{tabular}{c}0.680\\ $\pm$ 0.000\end{tabular}& \begin{tabular}{c}0.605\\ $\pm$ 0.001\end{tabular}&\begin{tabular}{c}0.631\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.648\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.562\\ $\pm$ 0.017\end{tabular}\\
\begin{tabular}{c}Large\\ Overlapping\end{tabular} & \begin{tabular}{c}\bf{0.837}\\ $\pm$ \bf{0.008}\end{tabular} &\begin{tabular}{c}0.692\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.591\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.612\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.670\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.566\\ $\pm$ 0.022\end{tabular}\\
\begin{tabular}{c}Scarce Regions \end{tabular}& \begin{tabular}{c}\bf{0.783}\\ $\pm$ \bf{0.007}\end{tabular} &\begin{tabular}{c}0.625\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.582\\ $\pm$ 0.001\end{tabular}&\begin{tabular}{c}0.526\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.608\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.525\\ $\pm$ 0.006\end{tabular}\\
\begin{tabular}{c}Within-pattern\\ Noise\end{tabular}& \begin{tabular}{c}\bf{0.818}\\ $\pm$ \bf{0.013}\end{tabular} &\begin{tabular}{c}0.675\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.626\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.619\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.651\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.550\\ $\pm$ 0.013\end{tabular}\\
\begin{tabular}{c}Mild Changes \end{tabular}& \begin{tabular}{c}\bf{0.762}\\ $\pm$ \bf{0.023}\end{tabular} &\begin{tabular}{c}0.611\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.574\\ $\pm$ 0.004\end{tabular}&\begin{tabular}{c}0.526\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.594\\ $\pm$ 0.000\end{tabular}&\begin{tabular}{c}0.539\\ $\pm$ 0.013\end{tabular}\\
\end{tabular}
\end{center}
\end{table}
We first validated that agreements among multiple trained models (Fig. \ref{fig2}b) are good indicators of model accuracy (Fig. \ref{fig2}a) for parameter selections. Besides, from Fig. \ref{fig2}a, we proved the model's ability in capturing ground truth of patterns and severity under different cases including pattern overlapping and large within-pattern noise. However, the model experienced declining performances when each pattern includes only scarce disease-related regions, and showed even worse performances when very mild atrophies are led by disease process. This is expected since, for the mild atrophy case, the simulated atrophy rate, 0-20\% of the regional volume, is even smaller than the standard deviation of most regional volumes among the normal population. Finally, Surreal-GAN significantly outperformed all five methods on five data sets we constructed (Table \ref{tab1}). However, none of these methods optimally utilize the CN data as a reference distribution, which is one of important limitations preventing them from precisely learning disease-specific pattern representations rather than capturing disease-unrelated confounding variations. The family of DNMF methods capture discriminant representation by minimizing and maximizing within-class and among-class variations respectively, thus did not capture heterogeneity within the PT class accurately either. Designed for parsing disease heterogeneity, semi-supervised methods focused on clustering patients into hard categorical subtypes (\cite{HYDRA}, \cite{Chimera}, \cite{MAGIC}, \cite{SmileGAN}). With different end points, these methods (including Smile-GAN) can not be directly compared with Surreal-GAN. To the best of our knowledge, our model is the first semi-supervised approach for learning disease-related representations, using a healthy reference group as a means for comparison. Therefore, to better understand the model performance, we leveraged the simulated ground truth information (not available in real experiments), derived "practical upper bounds" for compared unsupervised models, and also computed the gap between Surreal-GAN and supervised methods. (section B.5 and B.6)

\subsection{Results on Real Data}
\begin{figure}[t]
\centering
   \includegraphics[width=1\linewidth]{figure/figure3.png}
\caption{Voxel-wise statistical comparison results between selected subgroups and CN. (FDR method with p-value threshold of 0.05 was used for multiple comparison). In each direction guided by arrows, we increased values of one R-index and keep the other relatively fixed. Therefore, we chose out groups of subjects who have much stronger expression of one pattern for visualization.}
\label{fig3}
\end{figure}

On the AD dataset, models show highest agreement when the number of patterns $M=2$. From voxel-based comparison results (Fig. \ref{fig3}), we can clearly visualize two patterns related to AD. Each dimension of R-indices, ${\mathbf{r}}_i$, shows a clear positive correlation with the corresponding pattern severity. Pattern1, correlated with ${\mathbf{r}}_1$, shows diffuse cortical atrophies (e.g., bilateral insula, orbitofrontal and frontal poles), while pattern2 reveals focal atrophy in medial temporal lobes. Reproducibility of these two patterns were validated in section B.7. Furthermore, the 2D R-indices were effective and interpretable signatures for disease diagnosis and prognosis, achieving similar performance in AD classification compared to original 139 ROIs (section B.8).\\
Continuity and monotonicity of inferred R-indices further enable us to easily test correlation among them and other clinical variables (Table \ref{tab2}). Both ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$ are highly correlated with CSF-Abeta value and memory dysfunction measured by ADNI-MEM. However, ${\mathbf{r}}_1$ shows stronger correlations with white matter lesion volumes, presence of hypertension, executive function and language functioning, while ${\mathbf{r}}_2$ is more correlated with CSF-Tau, CSF-pTau and presence of APOE-E4 alleles.

\section{Conclusion}
In this study, we have proposed a novel method, Surreal-GAN, for learning representations of underlying disease-related imaging patterns. This model has overcome limitations in previous published semi-supervised clustering methods and shown great performance on semi-synthetic data sets. On AD dataset, the model parse heterogeneous PT data into two concise and clinically informative indices with preserved discriminant power, further proving its huge potential in uncovering disease heterogeneity from MRI data.
\newpage


\section{Introduction and statement of the results}
One of the most important problem in conformal geometry is the problem of finding conformal metrics with a prescribed curvature quantity. An example of curvature quantity which has received a lot of attention in the last decades is the Branson's \;$Q$-curvature. It is a Riemannian scalar invariant introduced by Branson-Oersted\cite{bo} (see also Branson\cite{bran1}) for closed four-dimensional Riemannian manifolds.
\vspace{8pt}

\noindent
Given \;$(M, g)$\; a four-dimensional closed Riemannian manifold with Ricci tensor \;$Ric_{g}$, scalar curvature \;$R_{g}$, and Laplace-Beltrami operator \;$\D_g$, the $Q$-curvature of \;$(M, g)$\; is defined by
\begin{equation}\label{eq:P}
Q_{g}=-\frac{1}{12}(\D_{g}R_{g}-R_{g}^{2}+3|Ric_{g}|^{2}).
\end{equation}
Under the conformal change of metric $g_u=e^{2u}g$\; with \;$u$\; a smooth function on \;$M$, the $Q$-curvature transforms in the following way
\begin{equation}\label{eq:P1}
P_gu+2Q_g=2Q_{g_{u}}e^{4u},
\end{equation}
where \;$P_g$ \;is the Paneitz operator introduced by Paneitz\cite{p1} and is defined by the following formula
\begin{eqnarray}\label{eq:P2}
P_{g}\varphi=\D_{g}^{2}\varphi+div_{g}\left((\frac{2}{3}R_{g}g-2Ric_{g})\n_g\varphi\right),
\end{eqnarray}
where \;$\varphi$\; is any smooth function on \;$M$, $div_g$\; is the divergence of with respect to \;$g$, and \;$\n_g$\; denotes the covariant derivative with respect to \;$g$. When, one changes conformally \;$g$\; as before, namely by \;$g_u=e^{2u}g$\; with \;$u$\; a smooth function on \;$M$, $P_g$\; obeys the following simple transformation law
\begin{equation}\label{eq:tranlaw}
P_{g_u}=e^{-4u}P_{g}.
\end{equation}
The equation \eqref{eq:P1} and the formula \eqref{eq:tranlaw} are analogous to classical ones which hold on closed Riemannian surfaces. Indeed, given a closed Riemannian surface $(\Sig, g)$ and \;$g_u=e^{2u}g$\; a conformal change of \;$g$\; with \;$u$\; a smooth function on \;$\Sig$, it is well know that
\begin{eqnarray}\label{eq:g}
\Delta_{g_u}=e^{-2u}\Delta_{g},\;\;\;\;\;\;\;\;\;-\Delta_{g}u+K_{g}=K_{g_u}e^{2u},
\end{eqnarray}
where for a background metric \;$\tilde g$\; on \;$\Sigma$ , \;$\Delta_{\tilde g}$\; and \;$K_{\tilde g}$\; are respectively the Laplace-Beltrami operator and the Gauss curvature of \;($\Sigma, \tilde g$). In addition to these, we have an analogy with the classical Gauss-Bonnet formula
$$
\int_{\Sigma}K_{g}dV_{g}=2\pi\chi(\Sigma),
$$
where\;$\chi(\Sigma)$\; is the Euler characteristic of \;$\Sigma$\; and \;$dV_g$\; is the volume form of \;$\Sig$\; with respect to \;$g$. In fact, we have the Chern-Gauss-Bonnet formula
$$
\int_{M}(Q_{g}+\frac{|W_{g}|^{2}}{8})dV_{g}=4\pi^{2}\chi(M),
$$
where \;$W_{g}$\; denotes the Weyl tensor of \;($M,g$) and \;$\chi(M)$\; is the Euler characteristic of  \;$M$. Hence, from the pointwise conformal invariance of  \;$|W_{g}|^{2}dV_{g}$, it follows that \;$\int_{M}Q_{g}dV_g$\; is also conformally invariant and will be denoted by \;$\kappa_g$, namely
\begin{equation}\label{eq:kappa}
\kappa_{g}:=\int_{M}Q_{g}dV_g.
\end{equation}
When \;$(M, g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$\; is the \;$4$-dimensional standard sphere, we have
\begin{equation}\label{kgsphere}
\kappa_g=\kappa_{g_{\mathbb{S}^4}}=8\pi^2.
\end{equation}
\vspace{10pt}

\noindent
Of particular importance in Conformal Geometry is the following  Kazdan-Warner type problem. Given a smooth positive function \;$K$ defined on a closed \;$4$-dimensional Riemannian manifold $(M, g)$, under which conditions on \;$K$\; there exists a  Riemannian metric conformal to \;$g$ \;with $Q$-curvature equal to \;$K$ . Thanks to \eqref{eq:P1}, the problem is equivalent to finding a smooth solution of the fourth-order nonlinear partial differential equation
\begin{eqnarray}\label{eq:qequation}
P_{g}u+2Q_{g}=2Ke^{4u}\;\;\;\;\;\;in\;\;M.
\end{eqnarray}
Equation \eqref{eq:qequation} is usually refereed to as the prescribed \;$Q$-curvatre equation and has been studied in the framework of Calculus of Variations, Critical Points Theory, Morse Theory and Dynamical Systems, see \cite{cy}, \cite{dm}, \cite{mal1}, \cite{ms}, \cite{nd5}, \cite{nd6}, \cite{nd7}, and the references therein.
\vspace{8pt}

\noindent
In this paper, we investigate equation \eqref{eq:qequation} in the context of Optimal Control Theory. Precisely, we study the following optimal control problem for the Paneitz obstacle problem
\begin{equation}\label{stateop}
\text{Finding} \;\;\;u_{min}\in H^2_Q(M)\;\;\text{such that}\;\;\; I(u_{min})=\min _{u\in H^2_Q(M)}I(u),
\end{equation}
where
 $$
 I(u)= \, \left<u,u\right>_{g}  - \kappa_{g} \log\left(\int_{M}Ke^{4T_g(u)} dV_g\right),\;\;\;\;u\in H^{2}_Q(M)
$$
with
$$
 <u,u>_{g}=\int_M\D_gu\D_gvdV_g+\frac{2}{3}\int_{M}R_g\n_g u\cdot\n_gvdV_g-\int_M2Ric_g(\n_g u, \n_gv)dV_g
$$
 $$
 T_g(u)=\arg{\min_{v\in H^2_Q(M),\\\;\;v\geq u}\left<v, v \right>_g}
 $$
 and
 $$H^2_Q(M):=\{u\in H^2(M):\;\; \int_MQ_gudV_g=0\}$$
 with \;$H^2(M)$\; denoting the space of functions on \;$M$\; which are of class \;$L^2$, together with their first and second derivatives. Moreover, the symbol $$\arg{\min_{v\in H^2_Q(M),\\\;\;v\geq u}\left<v, v \right>_g}$$ denotes the unique solution to the minimization problem \;$$\min_{v\in H^2_Q(M),\\\;\;v\geq u}\left<v, v \right>_g,$$ see Lemma \ref{obstacleq}. We remark that for \;$u$\; smooth,
 $$
 \left<u,u\right>_{g}=\left<P_gu, u\right>_{L^2(M)},
 $$
where \;$\left<\cdot, \cdot\right>_{L^2(M)}$\;denotes the \;$L^2$\; scalar product.
\vspace{6pt}

\noindent
 In the subcritical case, namely \;$0<\kappa_g<8\pi^2$, we  prove the following result.
\begin{thm}\label{positivevarmass}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq\mathbb{R}$, and \;$0<\kappa_g <8\pi^2$, then \;there exists \;$$u_{min}\in C^{\infty}(M)\cap H^2_Q(M)$$\; such that
$$
I(u_{min})=\min _{v\in H^2_Q(M)}I(v)\;\;\;\;{and}\;\;\;\;u_{min}= T_g(u_{min}).
$$
Moreover, setting $$
u_c=u_{min}-\frac{1}{4}\log\int_MKe^{4u_{min}}+\frac{1}{4}\log \kappa_g\;\;\;\text{and}\;\;\;\;g_c=e^{2u_c}g,
$$ we have
$$
Q_{g_c}=K.
$$
\end{thm}
\vspace{10pt}

\noindent
To state our existence result in the critical case, i.e $\kappa_g=8\pi^2$, we first set some notations.
We define \;$\mathcal{F}_K: M\longrightarrow \mathbb{R}$\; as follows
\begin{equation}\label{eq:limitfs}
\mathcal{F}_K(a):=2\left(H(a, a)+\frac{1}{2}\log(K(a))\right),\;\;\;a\in M
\end{equation}
where \;$H$\; is the regular part of the Green's function \;$G$\; of $P_g(\cdot)+2Q_g$ satisfying the  normalization $\int_M Q_g(x) G(\cdot, x)dV_g(x)=0$, see Section \ref{eq:notpre}. Furthermore, we define
\begin{equation}\label{eq:critfk}
Crit(\mathcal{F}_K):=\{a\in M: \;\;a\;\;\text{is critical point of} \;\;\mathcal{F}_K\}.
\end{equation}
Moreover, for \;$a\in M$\; we set
\begin{equation}\label{eq:partiallimit}
\mathcal{F}^a(x):=e^{(H(a, x)++\frac{1}{4}\log(K(x))},\;\;\;x\in M
\end{equation}
and define
\begin{equation}\label{eq:defindexa}
\mathcal{L}_K(a):=- \mathcal{F}^{a}(a)L_g(\mathcal{F}^{a})(a),
\end{equation}
where $$L_g:=-\D_g+\frac{1}{6}R_g$$ is the conformal Laplacian associated to \;$g$.  We set also
\begin{equation}\label{eq:critsett}
\mathcal{F}_{\infty}^+:=\{a\in Crit(\mathcal{F}_K):\:\;\mathcal{L}_K(a)>0\}.
\end{equation}
With this notation, our existence result in the critical case reads as follows:
\begin{thm}\label{positivevarmassi}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq\mathbb{R}$,  \;$\kappa_g =8\pi^2$, and \;$\mathcal{F}_{\infty}^+=Crit(\mathcal{F}_K)$, then \;there exists \;$$u_{min}\in C^{\infty}(M)\cap H^2_Q(M)$$\; such that
$$
I(u_{min})=\min _{v\in H^2_Q(M)}I(v)\;\;\;\;{and}\;\;\;\;u_{min}= T_g(u_{min}).
$$
Moreover, setting $$
u_c=u_{min}-\frac{1}{4}\log\int_MKe^{4u_{min}}+\frac{1}{4}\log \kappa_g\;\;\text{and}\;\;\;\;g_{c}=e^{2u_c}g ,
$$
we have
$$
Q_{g_c}=K.
$$
\end{thm}
\begin{rem}
\begin{itemize}
\noindent

\item The relation \;$u_{min}=T_g(u_{min})$\; in the above theorems is an additional information with respect to the existence results based on Calculus of Variations, Critical Points Theory, Morse Theory, and Dynamical Systems. It provides the inequality
\begin{equation}\label{obsinq}
\left<u_{min}, u_{min}\right>_g\leq \left<u, u\right>_g, \;\;\;\;\;\forall \;u_{min}\leq  u\in H^2_Q(M).
\end{equation}
\item
We remark that the nonlocal character of \;$e^{4T_g(u)}$\; in the definition of  \;$I$ with respect to  \;$e^{4u}$\; appearing in the definition of  \;$J$\; defined by
$$
J(u):= \, <u,u>_{g}  - \kappa_{g} \log\left(\int_{M}Ke^{u} dV_g\right),\;\;\;\;u\in H^{2}_Q(M)
$$
used in the existence approaches of \eqref{eq:qequation} via Calculus of Variations, Critical Points Theory, Morse Theory, and Dynamical Systems. The trade of the local character to non-local is in contrast with the traditional approach in the study of Differential Equations, but have the advantage of providing automatically the variational inequality \eqref{obsinq}.
\item
The \;$Q$-curvature functional \;$J$\; is invariant by translation by constants, while the \;$Q$-optimal control functional \;$I$\; is not. The functional \;$J$ is weakly lower semicontinuous, but the functional \;$I$\; is not. This makes it difficult to apply the Direct Methods in Calculus of Variations to study \eqref{stateop}.
\item
We expect formula \eqref{obsinq} to be useful to deal with the case \;$\kappa_g=8\pi^2$\; by helping to track down the loss of coercivity in the Variational Analysis of equation \eqref{eq:qequation}.
\end{itemize}
\end{rem}
\vspace{10pt}

\noindent
As a byproduct of our existence argument, we have the following regularity result for solutions of the optimal control problem \eqref{stateop}.
\begin{thm}\label{regularity}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq \mathbb{R}$,  \;$0<\kappa_g \leq 8\pi^2$, and \;$u\in H^2_Q(M)$\; is a minimizer of \;$I$\; on \;$H^2_Q(M)$, then
$$
u\in  C^{\infty}(M).
$$
\end{thm}
\vspace{10pt}

\noindent
An other consequence of our existence argument is the following compactness theorems for the set of minimizers of \;$I$\; on \;$H^2_Q(M)$.  We start with the subcritical case.
\begin{thm}\label{compsub}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq \mathbb{R}$, and \;$0<\kappa_g <8\pi^2$, then \;$\forall m\in \mathbb{N}$\; there exists \;$C_m>0$\; such that \;$\forall u \in C^{\infty}(M)\cap H^2_Q(M)$\; minimizer of \;$I$\;  on \;$H^2_Q(M)$, we have
$$
||u||_{C^k(M)}\leq C_m.
$$
\end{thm}
\vspace{8pt}

\noindent
For the critical case, setting
\begin{equation}\label{eq:critsett0}
\mathcal{F}_{\infty}^0:=\{a\in Crit(\mathcal{F}_K):\:\;\mathcal{L}_K(a)\neq 0\},
\end{equation}
we have:
\begin{thm}\label{compcri}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq\mathbb{R}$,  \;$\kappa_g =8\pi^2$, and \;$\mathcal{F}_{\infty}^0=Crit(\mathcal{F}_K)$, then \;$\forall m\in \mathbb{N}$\; there exists \;$C_m>0$\; such that \;$\forall u \in C^{\infty}(M)\cap H^2_Q(M)$\; minimizer of \;$I$\;  on \;$H^2_Q(M)$,  we have
$$
||u||_{C^m(M)}\leq C_m.
$$
\end{thm}
\vspace{10pt}

\noindent
We prove also some results in the particular case of the \;$4$-dimensional standard sphere, see Theorem \ref{sharpomt} and Corollary \ref{sfixedt} in Section \ref {casesphere}.
\vspace{10pt}

\noindent
To prove Theorem \ref{positivevarmass}-Theorem \ref{compcri}, we first use the variational characterization of the solution of  Paneitz obstacle problem \;$T_g(u)$ (see Lemma \ref{obstacleq})\;  to show that the Paneitz obstacle solution map \;$T_g$\; is  idempotent, i.e \;$T^2_g=T_g$, see Proposition \ref{nilpotent}. Next, using the idempotent property of \;$T_g$, we establish some monotonicity formulas, see Lemma \ref{decreasingformula}, Lemma \ref{minimaltunnel}, and Lemma \ref{decreasingformulaop}. Using the later monotonicity formulas, we show that  any minimizer of \;$J$\; or any solution of the optimal control problem \eqref{stateop} is a fixed point of \;$T_g$, see Corollary \ref{rigidityminimizer} and Corollary \ref{rigidityminimizeri}. This allows us to show that the \;$Q$-curvature functional \;$J$\; and the \;$Q$-optimal functional have the same minimizers on \;$H^2_Q(M)$, see Proposition \ref{sameminimizer} . With this at hand, Theorem \ref{positivevarmass} follows from the work of Chang-Yang\cite{cy} in the subcritical case, while Theorem \ref{positivevarmassi} follows from our work in the critical case in \cite{nd6}. Moreover, Theorem \ref{regularity} follows from the regularity result of Uhlenbeck-Viaclosky\cite{uv}. Furthermore, Theorem \ref{compsub} follows from the compactness result of Malchiodi\cite{mal} and Druet-Robert-\cite{dr}, while Theorem \ref{compcri} follows our compactness theorem in \cite{nd6}.
\vspace{10pt}

\noindent
The structure of the paper is as follows. In Section \ref{eq:notpre}, we collect some preliminaries and fix some notations. In Section \ref{opp}, we discuss the Paneitz obstacle problem and some monotonicity formulas involving the \;$Q$-curvature functional \;$J$. We also present some consequences of the latter monotonicity formulas. In Section \ref{socp}, we establish some monotonicity formulas for the \;$Q$-optimal functional \;$I$\; and their consequences as well. In Section \ref{pots}, we present the proof of Theorem \ref{positivevarmass}-Theorem \ref{compcri}. Finally, in Section \ref{casesphere}, we discuss the particular case of the \;$4$-dimensional standard sphere.
\section{Notations and Preliminaries}\label{eq:notpre}
 In this brief section, we fix our notations and give some preliminaries. First of all, from now until the end of the paper, \;$(M, g)$\;  and $\;K: M\longrightarrow \mathbb{R}_+$\; are respectively the given underlying closed four-dimensional Riemannian manifold and the smooth positive function to prescribe.
 \vspace{6pt}

  \noindent
We recall the function \;$J$\; used in other approaches to study \eqref{eq:qequation}.
\begin{equation}\label{eq:defj1}
J(u):= \, \left<u,u\right>_g \, + \, 4 \int_M Q_g u dV_g \,  - \,\kappa_{g} \log\left(\int_{M}Ke^{4u} dV_g\right),\;\;\;\;u\in H^2(M).
\end{equation}
Moreover, we recall the perturbed functional \;$J_t$ \;($0<t\leq 1$) which plays also an important role in the study of  minimizers of \;$J$.
\begin{equation}\label{eq:defjt1}
J_t(u):= \, \left<u,u\right>_g \, + \, 4t \int_M Q_g u dV_g \,  - \,t\kappa_{g} \log\left(\int_{M}Ke^{4u} dV_g\right),\;\;\;\;u\in H^2(M).
\end{equation}
We observe that
$$
J=J_1.
$$
Moreover, we define
$$
\overline{(u)}_Q=\frac{1}{\kappa_g}\int_MQ_gudV_g, \;\;\;\;u\in H^2(M),
$$
so that
$$
H^2_Q(M)=\{u\in H^2(M):\;\;\;\overline{(u)}_Q=0\}.
$$
For \;$a\in M$, we let \;$G (a, \cdot)$\; be the unique solution of the following system
\begin{equation}\label{eq:defG4}
\begin{cases}
P_g G(a, \cdot)+2Q_g(\cdot)=16\pi^2\delta_a(\cdot)\;\;\text{in}\;\;\;\;M\\
\int_M Q_g(x) G(a, x)dV_g(x)=0.
\end{cases}
\end{equation}
It is a well know fact that \;$G(\cdot, \cdot)$\; has a logarithmic singularity. In fact \;$G(\cdot, \cdot)$\; decomposes as follows
\begin{equation}\label{eq:decompG4}
 G(a,x)=\log \left(\frac{1}{\chi^2(d_{g}(a, x))}\right)+H(a, x),  \;\;\;x\neq a \in M.
\end{equation}
where \;$H(\cdot, \cdot)$ is the regular par of \;$G(\cdot, \cdot)$ and  \;$\chi$\; is some smooth cut-off function, see for example \cite{zw}.
\vspace{8pt}

\noindent
The decomposition of the Green's function \;$G$\; and the arguments of the proof of the Moser-Trudinger's inequality of  Chang-Yang\cite{cy} imply the following Moser-Trudinger type inequality.
\begin{pro}\label{cmt}
Assuming that  \;$P_g \geq 0$, $\ker P_g=\mathbb{R}$,  then \;there exists \; $C=C(M, g)>0$\; such that
$$
\log \int_Me^{4u}dV_g\leq C+\frac{1}{8\pi^2}\left<u, u\right>_g, \;\;\;\forall u\in H^2_Q(M).
$$
\end{pro}
\vspace{6pt}

\noindent
When \;$(M, g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$, we say \;$v$\; is a standard bubble if

\begin{equation}\label{standard}
P_{g_{\mathbb{S}^4}} v+6=6e^{4v} \;\;\text{on}\;\;\;\mathbb{S}^4.
\end{equation}
By the result of Chang-Yang \cite{cy1}, \;$v$\; satisfies
$$
e^{2v}g_{\mathbb{S}^4}=\varphi^*(g_{\mathbb{S}^4}),.
$$
for some \;$\varphi$\; conformal transformation of \;$\mathbb{S}^4$. It is well-known that the standard bubbles are related to the classical Moser-Trudinger-Onofri inequality. Indeed, we have:
\begin{pro}\label{cmto}
Assuming that \;$(M, g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$\; and \;$K=1$, then
\begin{equation}\label{icmto}
J(u)\geq 0,\;\;\;\;\forall u\in H^2(M).
\end{equation}
Moreover, equality in \eqref{icmto} holds if and only if \;$$v:=u-\frac{1}{4}\log \int _Me^{4u}+\frac{1}{4}\log \frac{\kappa_g}{3}$$ is a standard bubble.
\end{pro}
\vspace{8pt}

\noindent
To end this section, we say \;$w$\; is a \;$Q$-normalized standard bubble, if
\begin{equation}\label{qnormalized}
w=v-\overline{(v)}_Q,
\end{equation}
with \;$v$\;a standard bubble.
\section{Obstacle problem for the Paneitz operator}\label{opp}
In this section, we study the obstacle problem for the Paneitz operator. Indeed in analogy to the classical obstacle problem for the Laplacian, given \;$u\in H^2_Q(M)$, we look for a solution to the minimization problem
\begin{equation}\label{obspan}
\min_{v\in H^2_Q(M),\\\;\;v\geq u}\left<v, v \right>_g.\end{equation} We start with the following lemma providing the existence and unicity of solution for the obstacle problem for the Paneitz operator \eqref{obspan}.
\begin{lem}\label{obstacleq}
Assuming that \;$P_g \geq 0$\;  and \;$\ker P_g \simeq \mathbb{R}$, then \;$\forall u\in H^2_Q(M)$,
there exists a unique $T_g(u)\in H^2_Q(M)$ such that
\begin{equation}\label{tg}
\left<T_g(u),T_g(u) \right>_g=\min_{v\in H^2_Q(M),\\\;\;v\geq u}\left<v, v \right>_g
\end{equation}
\end{lem}
\begin{pf}
Since \;$P_g$\; is self-adjoint, $P_g\geq 0$\; and \;$\ker P_g\simeq\mathbb{R}$, then \;$<\cdot, \cdot>_g$\; defines a scalar product on \;$H^2_Q(M)$\; inducing a norm equivalent to the standard \;$H^2(M)$-norm on \;$H^2_Q(M)$. Hence, as in the classical obstacle problem for the Laplacian, the lemma follows from Direct Methods in the Calculus of Variations.
\end{pf}
\vspace{6pt}

\noindent
We study now some properties of the obstacle solution map \;$T_g\; : H^2_Q(M)\longrightarrow  H^2_Q(M)$. We start with the following algebraic one.
\begin{pro}\label{nilpotent}
Assuming that \;$P_g \geq 0$, \;$\ker P_g \simeq \mathbb{R}$, then the obstacle solution map \;$T_g\; : H^2_Q(M)\longrightarrow  H^2_Q(M)$\; is idempotent, i.e $$T^2_g=T_g.$$
\end{pro}

\begin{pf}
Let $v\in H^2_Q(M)$ such that $v\geq T_g(u)$. Then $T_g(u)\geq u$ implies
$v\geq u$. Thus by minimality, we obtain \;$$\left<v, v \right>_g \geq \left<T_g(u),T_g(u)\right>_g.$$ Hence, since \;$T_g(u)\geq T_g(u)$\; then by unicity we have
$$
T_g(T_g(u))=T_g(u),
$$
thereby ending the proof.
\end{pf}

\vspace{8pt}

\noindent
Next, we discuss some monotonicity formulas. We start with the following one.
\begin{lem}\label{decreasingformula}
Assuming that \;$P_g \geq 0$, \;$\ker P_g \simeq \mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then $$
J_t(u)-J_t(T_g(u))\geq \left<u, u\right>_g-\left<T_g(u)),T_g(u) \right>_g\geq 0, \;\;\;\forall u\in H^2_Q(M).
$$
\end{lem}
\begin{pf}
Using the definition of \;$J_t$ (see \eqref{eq:defjt1}), we have
\begin{equation}
J_t(u)-J_t(T_g(u))= \left<u, u\right>_g-\left<T_g(u),T_g(u) \right>_g-t\kappa_g\left(\log\frac{ \int_M Ke^{4u}dV_g}{ \int_M Ke^{4T_g(u)}dV_g}\right).
\end{equation}
Hence the result follows from \;$K>0$,  \;$T_g(u)\geq u$, and Lemma \ref{obstacleq} .
\end{pf}
\vspace{6pt}

\noindent
Lemma \ref{decreasingformula} imply the following rigidity result.
\begin{cor}\label{rigidity}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq\mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then \;\;$\forall u\in H^2_Q(M)$,
\begin{equation}\label{ineq}
J_t(T_g(u)) \leq  J_t(u)
\end{equation}
and
\begin{equation}\label{eql}
J_t(u)=J_t(T_g(u))\implies u = T_g(u).
\end{equation}
\end{cor}
\begin{pf}
Using lemma \ref{decreasingformula}, we have
\begin{equation}\label{eq1}
J_t(u)-J_t(T_g(u))\geq \left<u, u\right>_g-\left<T_g(u),T_g(u) \right>_g\geq 0.
\end{equation}
Thus, \eqref{ineq} follows from \eqref{eq1}. If \;$J_t(u)=J_t(T_g(u))$, then \eqref{eq1} implies
$$
\left<u, u\right>_g=\left<T_g(u),T_g(u) \right>_g.
$$
Hence, since \;$u\geq u$, then the unicity part in Lemma \ref{obstacleq} implies
$$
u=T_g(u),
$$
thereby ending the proof of the corollary.
\end{pf}

\vspace{8pt}

\noindent
Corollary \ref{rigidity} implies that minimizers of \;$J_t$\; on \;$H^2_Q(M)$\; are fixed points of  the obstacle solution map \;$T_g$. Indeed, we have:
\begin{cor}\label{rigidityminimizer}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq\mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then $$u\in H^2_Q(M)\;\;\;\text{is a minimizer of}\;\;\; J_t\implies u = T_g(u).$$
\end{cor}
\begin{pf}
$u\in H^2_Q(M)$\; is a minimizer of  \;$J_t$\; on \;$H^2_Q(M)$\; implies
\begin{equation}\label{eq3}
J_t(u)\leq J_t(T_g(u)).
\end{equation}
Thus combining \eqref{ineq} and \eqref{eq3}, we get
\begin{equation}\label{eq4}
J_t(u)= J_t(T_g(u)).
\end{equation}
Hence, combining \eqref{eql} and \eqref{eq4}, we obtain
$$
u=T_g(u).
$$
\end{pf}

\begin{rem}
Under the assumption of Corollary \ref{rigidity}, we have Proposition \ref{nilpotent} and Corollary \ref{rigidity} imply that we can assume without loss of generality that any minimizing sequence \;$(u_l)_{l\geq 1}$\; of \;$J_t$\; on  \;$H^2_Q(M)$\;satisfies
$$
u_l=T_g(u_l), \;\;\;\;\forall l\ge 1.
$$
\end{rem}

\section{Optimal control for the Paneitz operator}\label{socp}
In this section, we study a natural optimal control problem associated to the obstacle problem for the Paneitz operator . Indeed, we look for solutions of
$$\min_{u\in H^2_Q(M)}I(u),$$ where \;$I$\; is the $Q$-optimal control functional defined by

\begin{equation}\label{eq:defi12}
I(u):= \, \left<u,u\right>_g  - \kappa_{g} \log\left(\int_{M}Ke^{4T_g(u)} dV_g\right),\;\;\;\;u\in H^{2}_Q(M).
\end{equation}
Similarly to the \;$Q$-curvature functional $J$, for $0<t\leq 1$\; we define \;$I_t$\; by
\begin{equation}\label{eq:defi12t}
I_t(u):= \, \left<u,u\right>_g  - t\kappa_{g} \log\left(\int_{M}Ke^{4T_g(u)} dV_g\right),\;\;\;\;u\in H^{2}_Q(M).
\end{equation}
\vspace{6pt}

\noindent
We start with the following comparison result.
\begin{lem}\label{minimaltunnel}
Assuming that \;$P_g \geq 0$, $\ker P_g\simeq\mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then
$$
I_t\leq J_t\;\;\;\;\;\;\text{on}\;\;\;\;\; H^2_Q(M)\;\;\;\;\;
\text{and}\;\;\;\;\;\;
J_t\circ T_g=I_t\circ T_g\;\;\;\;\text{on}\;\;\; H^2_Q(M).
$$
\end{lem}
\begin{pf}
By definition of \;$J_t$\; and \;$I_t$ (see \eqref{eq:defjt1} and \eqref{eq:defi12t}), we have
$$
J_t(u)-I_t(u)=t\kappa_g\log \left(\frac{\int_{M}Ke^{4T_g(u)}}{\int_{M}Ke^{4u} } \right).
$$
Thus \;$I_t(u)\leq J_t(u)$\; follows from \;$T_g(u)\geq u$\; and \;$K>0$. Moreover, we have
$$
J_t(T_g(u))-I_t(T_g(u))=t\kappa_g\log \left(\frac{\int_{M}Ke^{4T^2_g(u)}}{\int_{M}Ke^{4T_g(u)} }\right) .
$$
Hence, \;$T_g^2=T_g$\; (see Lemma \ref{nilpotent}) implies \;$$J_t(T_g(u))=I_t(T_g(u)).$$
\end{pf}
\vspace{6pt}

\noindent
We have the following monotonicity formula for the \;$Q$-optimal control functional \;$I_t$.
\begin{lem}\label{decreasingformulaop}
Assuming that \;$P_g \geq 0$, \;$\ker P_g \simeq \mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then \;$\forall u\in H^2_Q(M)$,
$$
I_t(u)-I_t(T_g(u))=\left<u, u\right>_g-\left<T_g(u),T_g(u)\right>_g\geq 0.
$$
\end{lem}
\begin{pf}
By definition of \;$I_t$ (see \eqref{eq:defi12t}), we have
$$
I_t(u)-I_t(T_g(u))=\left<u, u\right>_g-\left<T_g(u),T_g(u) \right>_g-t\kappa_g\log \left(\frac{\int_{M}Ke^{4T_g(u)}}{\int_{M}Ke^{4T_g^2(u)} } \right).
$$
Using \;$T_g^2(u)=T_g(u)$\; and the definition of \;$T_g$\; (see Lemma \ref{obstacleq}), we  get
$$I_t(u)-I_t(T_g(u))=\left<u, u\right>_g-\left<T_g(u),T_g(u) \right>_g\geq 0.$$
\end{pf}
\vspace{6pt}

\noindent
Lemma \ref{obstacleq} and Lemma \ref{decreasingformulaop} imply that minimizers of \;$I_t$\; are fixed points of \;$T_g$.
\begin{cor}\label{rigidityminimizeri}
Assuming that \;$P_g \geq 0$, \;$\ker P_g=\mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then  $$u\in H^2_Q(M)\;\;\;\;\;\text{is a minimizer of }\;\;\;\;\; \;I_t\implies u = T_g(u).$$
\end{cor}
\begin{pf}
$u\in H^2_Q(M)$\;is a minimizer of \;$I_t$\; implies
$$I_t(u)\leq I_t(T_g(u)).$$
Thus Lemma \ref{decreasingformulaop} gives
$$
\left<u, u\right>_g=\left<T_g(u),T_g(u) \right>_g.
$$
Hence, by unicity we have
$$u=T_g(u).$$
\end{pf}
\vspace{6pt}

\noindent
\begin{rem}
Under the assumptions of Corollary \ref{decreasingformulaop}, we have that Proposition \ref{nilpotent} and Corollary \ref{decreasingformulaop} imply that for a minimizing sequence \;$(u_l)_{l\geq 1} $\; of \;$I_t$\; on  \;$H^2_Q(M)$, we can  assume without loss of generality that
$$
u_l=T_g(u_l), \;\;\;\;\forall l\ge 1.
$$
\end{rem}
\vspace{10pt}

\noindent
We have the following proposition showing that \;$I_t$\; and \;$J_t$\; have the same minimizers on \;$H^2_Q(M)$.
\begin{pro}\label{sameminimizer}
Assuming that \;$P_g \geq 0$, \;$\ker P_g\simeq\mathbb{R}$, \;$0<t\leq 1$\; and \;$0<\kappa_g\leq 8\pi^2$, then $$u\in H^2_Q(M)\;\;\text{is a minimizer of}\;\;J_t\;\;\text{is equivalent to}\;\;u\in H^2_Q(M)\;\;\text{is a minimizer of }\;\;I_t.$$
\end{pro}
\begin{pf}
Suppose $u\in H^2_Q(M)$\; is a minimizer of \;$J_t$. Then Corollary \ref{rigidityminimizer} implies
$$
u=T_g(u).
$$
Thus using Lemma \ref{minimaltunnel} we have
$$
I_t(u)=J_t(u)
$$
For $v\in H^2_Q(M)$, we have Lemma \ref{minimaltunnel}, Lemma \ref{decreasingformulaop}, and $u\in H^2_Q(M)$\; is a minimizer of \;$J_t$\; imply
$$
I_t(v)\geq I_t(T_g(v))=J_t(T_g(v))\geq J_t(u)=I_t(u).
$$
Hence $u\in H^2_Q(M)$\; is a minimizer of \;$I_t$\;on \;$H^2_Q(M)$. Similarly, suppose  $u\in H^2_Q(M)$\; is a minimizer of \;$I_t$. Then Corollary \ref{rigidityminimizeri} implies
$$
u=T_g(u).
$$
Thus using again Lemma \ref{minimaltunnel} we have
$$
I_t(u)=J_t(u).
$$
For $v\in H^2_Q(M)$, we have Lemma \ref{minimaltunnel}, Lemma \ref{decreasingformula}, and \;$u\in H^2_Q(M)$\; is a minimizer of \;$I_t$\; imply
$$
J_t(v)\geq J_t(T_g(v))=I_t(T_g(v))\geq I_t(u)=J_t(u).
$$
Hence \;$u\in H^2_Q(M)$\; is a minimizer of \;$J_t$\; on \;$H^2_Q(M)$.
\end{pf}
\vspace{8pt}

\noindent

\section{Proof of Theorem \ref{positivevarmass} -Theorem \ref{compcri}}\label{pots}
In this section, we present the proof of Theorem \ref{positivevarmass} -Theorem \ref{compcri}. As already mentioned in the introduction, the proofs are based on Proposition \ref{sameminimizer} and some contributions of Chang-Yang\cite{cy}, Druet-Robert\cite{dr}, Malchiodi\cite{mal}, the author\cite{nd6} and Uhlenbeck-Viaclovsky\cite{uv} in the the study of the fourth-order nonlinear partial differential equation \eqref{eq:qequation}.
\vspace{8pt}

\noindent
\begin{pfn}{ of Theorem \ref{positivevarmass}}\\
Since \;$P_g \geq 0$, $\ker P_g=\mathbb{R}$, and \;$0<\kappa_g <8\pi^2$, then the works of Chang-Yang\cite{cy} and Uhlenbeck-Viaclosvky\cite{uv} imply the existence of \;$u_0\in C^{\infty}(M)$ such that
$$
J(u_0)=\min_{u\in H^2(M)} J(u).
$$
Since \;$J$\; is translation invariant, then setting $$u_{min}=u_0-\overline{(u_0)}_Q,$$ we have $$u_{min}\in C^{\infty}(M)\cap H^2_Q(M)$$ and
$$
J(u_{min})=\min_{u\in H^2_Q(M)} J(u).
$$
Using Proposition \ref{sameminimizer}, we get
$$
I(u_{min})=\min_{u\in H^2_Q(M)} I(u)
$$
Thus Corollary \ref{rigidityminimizeri} implies
$$
u_{min}=T_g(u_{min}).
$$
Recalling that $$J(u_{min})=J(u_0)=\min_{u\in H^2(M)} J(u),$$ we have
$$
P_gu_{min}+2Q_g=2\kappa_g\frac{Ke^{4u_{min}}}{\int_MKe^{4u_{min}}}.
$$
Thus, setting
$$
u_c=u_{min}-\frac{1}{4}\log\int_MKe^{4u_{min}}+\frac{1}{4}\log \kappa_g,
$$ we have
$$
P_gu_{c}+2Q_g=2Ke^{4u_{c}}.
$$
Hence, setting
$$
g_{u_c}=e^{2u_c}g,
$$
we obtain
$$
Q_{g_{u_c}}=K.
$$
thereby ending the proof.
\end{pfn}

\begin{pfn} {of Theorem \ref{positivevarmassi}}\\
Let \;$\varepsilon_l\in (0, 1)$\; with \;$\varepsilon_l\rightarrow 0$. For \;$l\geq 1$, we define
$$
J_l:=J_{1-\varepsilon_l}\;\;\;\text{and}\;\;\;I_l:=I_{1-\varepsilon_l}
$$
As in the proof of Theorem \ref{positivevarmass}, for \;$l\geq 1$\; the works of Chang-Yang\cite{cy} and Uhlenbeck-Viaclosvky\cite{uv} give the existence of $$u^{l}_{min}\in C^{\infty}(M)\cap H^2_Q(M)$$ such that
\begin{equation}\label{eq:minjl}
J_{l}(u^{l}_{min})=\min_{u\in H^2(M)} J_{l}(u).
\end{equation}
Thus, using Proposition \ref{sameminimizer}, we get
\begin{equation}\label{minoptep}
I_{l}(u^{l}_{min})=\min_{u\in H^2_Q(M)} I_{l}(u).
\end{equation}
Clearly \eqref{eq:minjl} imply,
\begin{equation}\label{presqlimit}
P_g u^{l}_{min}+2Q_g(1-\varepsilon_l)=2\kappa_g(1-\varepsilon_l)\frac{Ke^{4u^{l}_{min}}}{\int_MKe^{4u^{l}_{min}}}.
\end{equation}
Hence, setting
\begin{equation}\label{ul}
u_c^{l}=u_{min}^{l}-\frac{1}{4}\log\int_MKe^{4u_{min}^{l}}+\frac{1}{4}\log \kappa_g,
\end{equation}
we obtain
\begin{equation}\label{eulerl1}
P_gu_{c}^{l}+2Q_g(1-\varepsilon_l)=2K(1-\varepsilon_l)\e^{4u_{c}^{l}}.
\end{equation}
Thus our bubbling rate formula  in \cite{nd6} and the assumption \;$\mathcal{F}_{\infty}^+=Crit(\mathcal{F}_K)$\;  prevents the sequence \;$u_{c}^{l}$\; from bubbling. Hence we have
\begin{equation}\label{smoothc}
u_{c}^{l}\longrightarrow u_c \;\;\;\;\text{smoothly, as } \;\;l\longrightarrow \infty.
\end{equation}
Thus \eqref{eulerl1} gives
 \begin{equation}\label{euler1}
 P_gu_{c}+2Q_g=2Ke^{4u_{c}}.
 \end{equation}
Recalling \;$u_{min}^l\in H^2_Q(M)$, we have  \eqref{ul} and \eqref{smoothc} imply
\begin{equation}\label{smoothconv}
u_{min}^{l}\longrightarrow u_{min} \;\;\;\;\text{smoothly}.
\end{equation}
and
$$
u_c=u_{min}-\frac{1}{4}\log\int_MKe^{4u_{min}}+\frac{1}{4}\log \kappa_g.
$$
Clearly \eqref{smoothconv} and \eqref{minoptep} imply
$$
I(u_{min})=\min_{u\in H^2_Q(M)} I(u).
$$
Hence Corollary \ref{rigidityminimizeri} and \eqref{euler1} imply
$$
u_{min}=T_g(u_{min}).
$$
and
$$
Q_{g_{u_c}}=K.
$$
\end{pfn}
\vspace{6pt}

\noindent
\begin{pfn}{ of Theorem \ref{regularity}}\\
It follows  directly from Proposition \ref{sameminimizer}, the translation invariant property of \;$J$\; and the regularity result of Uhlenbeck-Viaclovsky\cite{uv}.
\end{pfn}
\vspace{6pt}

\noindent
\begin{pfn}{ of Theorem \ref{compsub}}\\
Let $u\in C^{\infty}(M)\cap H^2_Q(M)$ be a minimizer of \;$I$\; on \;$H^2_Q(M)$. Then  the translation invariance property of \;$J$\; and Proposition \ref{sameminimizer} imply \;$u$\; is a minimizer of \;$J$\; on \;$H^2(M)$. Hence \;$u$\; satisfies
$$
P_gu+2Q_g=2\kappa_g \frac{Ke^{4u}}{\int_MKe^{4u}}.
$$
Then, setting
\begin{equation}\label{defiv}
v=u-\frac{1}{4}\log \int_{M}Ke^{4u}+\frac{1}{4}\log \kappa_g,
\end{equation}
we get
$$
P_gv+2Q_g=2 Ke^{4v}
$$
Thus, since $0<\kappa_g<8\pi^2$, then the compactness result of Malchiodi\cite{mal} and Druet-Robert\cite{dr} imply \;$\forall m\in \mathbb{N}$, there exists \;$\tilde C_m>0$\; such that
$$
||v||_{C^m(M)}\leq \tilde C_m.
$$
Hence, \;$u\in H^2_Q(M)$\; and \; \eqref{defiv} give the existence of  \;$C_m>0$\; such that
$$
||u||_{C^m(M)}\leq C_m,
$$
thereby ending the proof.
\end{pfn}
\vspace{6pt}

\noindent
\begin{pfn}{ of Theorem \ref{compcri}}\\
The proof is a small modification of the one of Theorem \ref{compsub}. For the sake of completeness, we repeat all the steps. Let $u\in C^{\infty}(M)\cap H^2_Q(M)$ be a minimizer of \;$I$\; on \;$H^2_Q(M)$. Then  as in the proof of Theorem \ref{compsub}, \;$u$ is a minimizer of \;$J$\; on \;$H^2(M)$. Hence \;$u$\; satisfies
$$
P_gu+2Q_g=2\kappa_g \frac{Ke^{4u}}{\int_MKe^{4u}}.
$$
Then, setting
$$
v=u-\frac{1}{4}\log \int_{M}Ke^{4u}+\frac{1}{4}\log \kappa_g,
$$
we get
$$
P_gv+2Q_g=2 Ke^{4v}.
$$
Thus since \;$\mathcal{F}_{\infty}^0=Crit(\mathcal{F}_K)$, then our compactness theorem in \cite{nd6} imply that \;$\forall m\in \mathbb{N}$\; there exists \;$\tilde C_m>0$ such that
$$
||v||_{C^m(M)}\leq \tilde C_m.
$$
Hence recalling that \;$u\in H^2_Q(M)$, we have there exists \;$C_m>0$\; such that
$$
||u||_{C^m(M)}\leq C_m.
$$
\end{pfn}
\section{Obstacle problem and Moser-Trudinger type inequality}\label{casesphere}
In this section, we discuss some Moser-Trudinger type inequalities related to the Paneitz obstacle problem. In particular, we specialize to the case of the \;$4$-dimensional standard sphere \;$(\mathbb{S}^4, g_{\mathbb{S}^4})$.
\vspace{6pt}

\noindent
We have the following obstacle Moser-Trudinger type inequality.
\begin{pro}
Assuming that  \;$P_g \geq 0$, $\ker P_g=\mathbb{R}$,  then there exists \;$C=C(M, g)>0$\; such that
$$
\log \int_Me^{4T_g(u)}dV_g\leq C+\frac{1}{8\pi^2}\left<u, u\right>_g, \;\;\;\forall u\in H^2_Q(M).
$$
\end{pro}
\begin{pf}
Clearly \;$u\leq T_g(u)$\; gives
\begin{equation}\label{imp0}
\log \int_Me^{4u}dV_g\leq \log \int_Me^{4T_g(u)}dV_g.
\end{equation}
 Since \;$P_g \geq 0$\; and \;$\ker P_g=\mathbb{R}$, then \;the classical Moser-Trudinger inequality in Proposition \ref{cmt} implies the existence of \; $C=C(M, g)>0$\; such that
\begin{equation}\label{imp1}
\log \int_Me^{4T_g(u)}dV_g\leq C+\frac{1}{8\pi^2}\left<T_g(u), T_g(u)\right>_g.
\end{equation}
Using the definition of \;$T_g$, we get
\begin{equation}\label{imp2}
\left<T_g(u), T_g(u)\right>_g\leq \left<u, u\right>_g.
\end{equation}
Hence combining \eqref{imp0}-\eqref{imp2}, we get
$$
\log \int_Me^{4T_g(u)}dV_g\leq C+\frac{1}{8\pi^2}\left<u, u\right>_g.
$$
\end{pf}
\vspace{8pt}

\noindent
When \;$(M, g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$\; and  \;$K=1$, we have the following sharp obstacle Moser-Trudinger type inequality.
\begin{thm}\label{sharpomt}
Assuming that \;$(M,  g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$\; and \;$K=1$, then
$$I\geq 0 \;\;\;\text{on}\;\;\; H^2_Q(M),$$ i.e
\begin{equation}\label{opmt}
\log \int_Me^{4T_g(u)}dV_g\leq \frac{1}{8\pi^2}\left<P_gu, u\right>, \;\;\;\forall u\in H^2_Q(M).
\end{equation}
Moreover equality in \eqref{opmt}  holds if and only if \;$$v:=u-\frac{1}{4}\log \int _Me^{4u}+\frac{1}{4}\log \frac{\kappa_g}{3}$$ is a standard bubble, see \eqref{standard} for its definition.
\end{thm}
\begin{pf}
Since \;$(M, g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$\; and \;$K=1$, then by the classical Moser-Trudinger-Onofiri ineqality in Proposition \ref{cmto}, we have
\begin{equation}\label{imp3}
J\geq 0 \;\;\;\text{on}\;\;\; \;H^2(M)
\end{equation}
and \;\begin{equation}\label{impes}
J(u)=0\;\; \;\text{is equivalent to }\;\;\;v:=u-\frac{1}{4}\log \int _Me^{4u}+\frac{1}{4}\log \frac{\kappa_g}{3}\;\; \;\text{is a standard bubble}.
\end{equation}
 Using Lemma \ref{decreasingformulaop}, we get
\begin{equation}\label{imp4}
I\geq I\circ T_g\; \;\;\text{on}\;\;\; H^2_Q(M).
\end{equation}
Thus, using Lemma \ref{minimaltunnel} and \eqref{imp4}, we have
\begin{equation}\label{imp5}
I\geq J\circ T_g \;\;\;\text{on}\;\;\; H^2_Q(M).
\end{equation}
So, combining \eqref{imp3} and \eqref{imp5}, we get
\begin{equation}\label{ipositive}
I\geq 0 \;\;\;\text{on}\;\;\; H^2_Q(M).
\end{equation}
Hence, recalling the definition of \;$I$ (see \eqref{eq:defi12}) and \eqref{kgsphere}, we  have \eqref{ipositive} is equivalent to
$$
\log \int_Me^{4T_g(u)}dV_g\leq \frac{1}{8\pi^2}\left<u, u\right>_g, \;\;\;\forall u\in H^2_Q(M).
$$
Suppose \;$$v:=u-\frac{1}{4}\log \int _Me^{4u}+\frac{1}{4}\log \frac{\kappa_g}{3}$$is a standard bubble with\;$u\in H^2_Q(M)$. Then \eqref{impes} implies
\begin{equation}\label{impes1}
J(u)=0
\end{equation}
Thus  \eqref{impes1}, Lemma \ref{minimaltunnel} and the first part (namely \eqref{ipositive})  imply
$$
I(u)=0.
$$
Hence we have the equality case in \eqref{opmt} . Suppose we have the equality case in \eqref{opmt} with \;$u\in H^2_Q(M)$. Then
\begin{equation}\label{izero}
I(u)=0.
\end{equation}
Thus, using \eqref{ipositive} and \eqref{izero}  we get
\begin{equation}\label{minimum0}
I(u)=\min_{v\in H^2_Q(M)}I(v).
\end{equation}
Using \eqref{minimum0} and Corollary \ref{rigidityminimizeri}, we obtain
\begin{equation}\label{fixedpt}
u=T_g(u).
\end{equation}
So Lemma \ref{minimaltunnel} , \eqref{izero} and \eqref{fixedpt} imply
\begin{equation}\label{jzero}
J(u)=0.
\end{equation}
Hence using \eqref{impes} and \eqref{jzero}, we have \;$v:=u-\frac{1}{4}\log \int _Me^{4u}+\frac{1}{4}\log \frac{\kappa_g}{3}$\; is a standard bubble.
\end{pf}
\vspace{8pt}

\noindent
Theorem \ref{sharpomt} implies the following corollary stating that $Q$-normalized standard bubbles (see \eqref{qnormalized} for their definitions) are fixed points of the obstacle solution map \;$T_g$.
\begin{cor}\label{sfixedt}
Assuming that \;$(M,  g)=(\mathbb{S}^4, g_{\mathbb{S}^4})$\; and \;$w$\; is a \;$Q$-normalized standard bubble (see \eqref{qnormalized} for its definition), then
$$
T_g(w)=w.
$$
\end{cor}
\begin{pf}
Since \;$w$\; is a \;$Q$-normalized standard bubble, then
  $$w:=v-\overline{(v)}_Q$$ with \;$v$\; is a standard bubble. Thus, Lemma \ref{minimaltunnel}, Theorem \ref{sharpomt}, and the translation invariant property of \;$J$\; imply
$$
0\leq I(w)\leq J(w)=J(v)=0.
$$
Using again Theorem \ref{sharpomt}, we obtain
$$
I(w)=\min_{v\in H^2_Q(M)}I(v)
$$
Hence using Corollary \ref{rigidityminimizeri}, we get
$$
w=T_g(w).
$$
\end{pf}