\section{Motivation}
After first being discovered in the late 1960s \cite{AllanNature}, radio
emission of air showers has again received increasing interest in the past
years. This is mainly due to the achievements in information technology,
which open the possibility to use phased arrays of radio antennas in
combination with digital beam-forming. A first proof of this technology has been
provided by the {\sc Lopes} experiment \cite{Falcke:2002tp,Falcke:2005tc}. Its main
advantage are low costs and simple design of the radio antennas as
compared to classical air shower detector elements (such as scintillators or
photomultipliers) and the very large field of view provided by the phased array
(usually close to $2\pi$). This technology is now not only being developed as an
extension at many existing air shower experiments ({\sc Kascade, Auger,
Tunka,...}), but also forms the key building block of
{\sc Lofar} -- a multi-site, multi-purpose radio detector spanning half of
Europe \cite{Lofar:2005wt}.\\
In this paper, we will first outline the science case for an extension of the
\textsc{IceCube}\xspace observatory with an extended $\O(10\un{km^2})$ radio array, then describe
the current status and finally outline a roadmap towards a full detector.

\section{Science Case}
The \textsc{IceCube}\xspace observatory \cite{Berghaus:2008bk} at the South Pole was originally designed and intended
for detection
of $\un{TeV}-\un{PeV}$ neutrinos penetrating Earth, while at the same time using
it as a shield against the abundant atmospheric muon flux. Beyond the original
anticipation, \textsc{IceCube}\xspace has by now become a $4\pi$ observatory. In
particular, the \textsc{IceTop}\xspace air-shower array \cite{Stanev:2009ce} on the surface,
expands the physics reach of \textsc{IceCube}\xspace in a threefold manner, as detailed in the
sections below.
\onefig{muons}
{Detectable muons in \textsc{IceCube}\xspace vs. electrons in \textsc{IceTop}\xspace for different
energies and primaries. The dashed line shows an analytic calculation from
\citet{Halzen:1996ik}.\label{fig:compo}}
\subsection{Composition\label{sec:composition}}
Two different observables are currently employed to determine primary composition
of cosmic rays:
\begin{itemize}
\item the depth of the {\bf shower maximum $X_{max}$}, where on
average protons of the same energy will penetrate deeper in the atmosphere and
show larger variations than heavier nuclei. This is used to measure composition
e.g. with the fluorescence detectors in \textsc{Auger}\xspace.  \item the {\bf electron over
muon ratio}, where heavier primaries will
cause higher muon multiplicities (c.f. fig.\ref{fig:compo}), as exploited e.g. by
\textsc{Kascade}\xspace.
\end{itemize}
Significant statistics of coincident \textsc{InIce}\xspace and \textsc{IceTop}\xspace events have been
collected and are being analyzed. Due to the strong correlation of the muon and
electron flux with the total primary energy and inclination of the shower,
very good control of systematic effects in the reconstruction of the events is
required \cite{Klepser:2008pm}.

For air showers that are coincidentally detected by an additional future large
radio array, both -- the muon and the electron component of the shower -- will
be measured by more than one type of detector, allowing an overconstrained
determination of both components. Furthermore, simulations show that the steepness
of the lateral distribution of the geosynchrotron signal will depend on the
depth of the shower maximum \cite{Huege:2005rd,deVries:2010ti} -- an effect that will be
enhanced for higher observation levels \cite{Huege:2006kd}. Radio detection may
thus provide a view of the shower development at high altitudes. It should be noted,
though, that the simulation of the geosynchrotron emission is still a developing
field. The reconstruction of the radio signals will thus have to be validated
against the \textsc{IceTop}\xspace measurements first.

The benefit will go beyond the independent systematic
errors of the three detector types on the shower energy and directional
reconstruction. With the increased statistics provided by a large array, the
energy range that is covered with good statistics will be extended to the
$\un{EeV}$ regime, overlapping and complementing not only the measurements
provided by the \textsc{Kascade}\xspace but also of the \textsc{Auger}\xspace experiments. Using the
additional information on the shower maximum, a combined measurement will
eventually maximally exploit all handles on primary particle composition that
are currently available to earth-bound detectors. \\

\subsection{Neutrinos}
At energies up to $\O(\un{PeV})$, the neutrino field of view of \textsc{IceCube}\xspace is
limited to the $2\pi\un{sr}$ of the northern hemisphere due to the abundant flux
of atmospheric muons from air showers from above. This flux shows a much steeper
spectrum than the expected signal flux. At energies above several
$\O(100)\un{TeV}$ the neutrinos themselves get absorbed in Earth, so that
the search must be extended to the southern hemisphere.  Still, even at these
higher energies, the dominant background is from unresolved
muon bundles induced by air showers, which can not be easily distinguished from
single muons in \textsc{IceCube}\xspace.
While \textsc{IceTop}\xspace in principle offers the possibility to detect the air
showers from which these muons stem, its vetoing power is limited by its small
size compared to the aperture available to \textsc{IceCube}\xspace to around $0.25\pi\un{sr}$. A radio array with
sufficient lateral extent can significantly increase the field of view in which
air shower induced muons can be vetoed. An estimate shows that the
sensitivity to the GZK neutrino flux may be increased by a factor of more than
three for an array extending to $3\un{km}$ around the center of \textsc{IceTop}\xspace
\cite{Jan:2010}.

\subsection{Ultra-high energy gamma rays}
Photon induced air showers will show a muon flux that is a factor of $\approx 100$
lower than for hadron induced showers \cite{Drees:1988wu} (c.f. fig.\ref{fig:compo}). At the same time, the flux of
muons induced by hadron showers is energetic enough that \textsc{IceCube}\xspace will have
a nearly $100\%$ detection efficiency for showers of energies above
$10^{15}\un{eV}$ if the shower axis intersects the geometric volume. This can
be exploited to discriminate hadronic air showers against electromagnetic
showers, turning the observatory in a $\un{km^2}$-scale $\un{PeV}$ photon
detector.

While at $\un{PeV}$ energies, the UHE gamma horizon is limited to our own galaxy
due to absorption on the cosmic microwave background, both the \textsc{H.E.S.S.}\xspace and \textsc{Magic}\xspace
experiment report unbroken $E^{-2}$ spectra for some sources within our galaxy
\cite{Abdo:2007ad,Hess:2007}.  Photons have been detected from the same sources with energies
up to energies of $90\un{TeV}$ \cite{Abdo:2007ad}. If point sources of gamma
rays with unbroken spectra up to $\un{PeV}$ exist, they should generate $\O(10)
\un{evt}\un{yr^{-1}}$ in an array as described above.  Non-excluding limits on
the diffuse photon fraction in cosmic rays have been set previously at
$<10^{-4}$ for the energy range up $5\cdot10^{16}\un{eV}$ by the \textsc{Casa-Mia}\xspace
collaboration\cite{Chantell:1997gs} and at $<10^{-1.6}$ for energies above
$2\cdot10^{18}\un{eV}$ by \textsc{Auger}\xspace\cite{Auger:2009qb}.  Ongoing
analyses \cite{BuitinkPriv} with \textsc{IceTop}\xspace show that a comparably high sensitivity can be
reached, but only for the limited section of the galactic plane covered by the
\textsc{IceTop}\xspace-\textsc{InIce}\xspace field-of-view.

Again a future radio array will enhance the physics reach for UHE gamma
detection. As the extent of a future radio array will be much larger than the
extent of \textsc{IceTop}\xspace, sources at higher declination -- including the entire
galactic plane -- can be probed. At the same time, the larger collection area
will increase event statistics, and will thus bridge the gap in energy
between \textsc{Casa-Mia}\xspace and \textsc{Auger}\xspace. Finally, as is discussed in section
\ref{sec:composition}, the
radio signal can additionally provide a handle on the depth of the shower
maximum, hence further improve the photon identification in this search.

\onefig[t]{antenna2}
  {Custom {\em Fat Wire-Dipole} (FWD) antenna\label{fig:antenna}}
\section{Current status and experimental setup}
A set of dedicated measurements has been performed to explore the possibility of
such a radio array in Antarctic conditions~\cite{Jan:2010}. A {\em fat wire-dipole} antenna
depicted in Figure \ref{fig:antenna} has been developed with the aim of a high
and uniform broadband response in the $\un{MHz}$-regime while limiting antenna
dimensions and cost.

Two of these antennas were deployed at South Pole in the polar season 2008/2009
around the SPASE building and were read out using two channels of the
pre-existing RICE DAQ \cite{Kravchenko:2007fy}. The sensitivity of the deployed instrumentation was severely
limited by cable losses as only two $600\un{m}$ and $500\un{m}$ pieces of legacy
$75\un{\Omega}$ coaxial cable were available for use. In the frequency range
between 50 and 100 MHz signal transmission losses vary from 20-30 dB.
In absence of a dedicated trigger for the surface antennas the $> 2
\un{\mu{}s}$ cable delays also resulted in signals that were often truncated by the
limited buffer length of the RICE waveform capture.
Nevertheless, the essential pre-conditions for establishing an enlarged
footprint at South Pole have been achieved via two fundamental performance
milestones.

\onefig[t]{spectrum1}
  {$8\un{\mu{}s}$-sweep noise spectrum from fat wire-dipole\label{fig:noisespectrum}}
1) Analysis of the spectrum shows the noise floor being limited by galactic
and thermal contributions. Figure \ref{fig:noisespectrum} shows the background spectrum
from a single $8\un{\mu{}s}$ sweep corrected for the attenuation effects of
cable, filter, and amplifier. The strong increase below
$25\un{MHz}$ is an artefact of the high-pass filter being used. As in a previous short-term
measurement during the 2007/2008 season, the spectrum lacks the typical
monochromatic features from TV stations and other anthropogenic sources observed
in other more inhabitated areas.

2) To test the ability of surface antennas to identify sources at South Pole, a set of
calibration pulser runs were taken in January 2009 using a RICE transmitter
located at the $500\un{m}$ distant MAPO
building. Taking advantage of the nanosecond time resolution of the RICE DAQ,
the source position has been reconstructed using only information from the two
fat wire-dipoles plus two auxilliary dipole antennas at
the roof of MAPO building. The vast majority of the events reconstruct to within
$\sim 2\un{m}$ of the actual pulse location.

\section{Roadmap}
Here we outline a roadmap to assess the viability of an extended air-shower
array for geosynchrotron radiation as described above. Considering the
restricted hardware deployment cycle due to the remote location, a
staged proposal spanning three years has been made.
The setup in the first year will be aimed towards unambiguously establishing the
detection of air-showers above the background radio noise and thus provide a
proof-of-viability for this detection technique in the Antarctic
environment. Using a conservative estimate based on REAS2 simulations
\cite{Jan:2010}, air-showers with a primary energy of $E_{prim} > 10^{17}
\un{eV}$ should be detectable at $> 5\sigma$ above the ambient noise level at a
distance of $125\un{m}$.
Antennas will be deployed in pairs in orthogonal
polarization. Four such pairs at a baseline of $55\un{m}$ and requiring a four-fold coincidence
will provide an effective area of $3\cdot10^{4}\un{m^2}$. Using the charged cosmic ray flux as measured by the
Kascade experiment \cite{Ulrich:2002}, each of these clusters will detect $\sim 8
\un{events/day}$. Even allowing for considerable inefficiencies (e.g. due to
anthropogenic backgrounds), some $10^{3}$ events will be accumulated per year.

While this primary setup will heavily rely on
commercially available components and existing infrastructure, the setup in the second season
should employ all of the key technologies that are required for a
several-$\un{km}^2$ array, hence will also allow to demonstrate the scalability
of the approach. In particular these include each antenna pair to be read out by
an independent, custom DAQ board deployed in close proximity. Connecting the
stations in a peer-to-peer grid, in which each of the boards generates threshold
triggers and communicates them to neighboring antennas will allow to establish
local coincidences. Local storage of the waveform information with a buffer
depth of several $100\un{ms}$ will be required to read out the array following a
global trigger either generated internally from the radio
array or externally from \textsc{IceCube}\xspace or \textsc{IceTop}\xspace. With a significantly larger number
of antennas on an enlarged footprint, this second stage will collect a
data sample large enough to allow detailed verification of the array
performance and shower reconstruction.

The third year will be dedicated to investigation of long-term systematics. With
data accumulated over two years, potential environmental changes (such as
compactification of the snow around the antennas) or degradation of the antenna
system will be studied in detail. Any needed in-situ recalibration of the
antennas and comparison to previous calibrations would be conducted in this
campaign. For this purpose a radio-frequency transmitter can be deployed at
different locations in the array.

\section{Summary}
\textsc{IceCube}\xspace{}'s core mission includes several cosmic ray initiatives. While combined analysis of
\textsc{IceTop}\xspace and \textsc{InIce}\xspace data is being used to target cosmic-ray composition, there will continue to be
uncertainties driven by the ability to completely characterize individual events.
An extended air-shower radio array using geosynchrotron radiation will significantly enhance the
acceptance of \textsc{IceCube}\xspace for the veto and gamma-ray missions. It will also provide complementary
information for shower reconstruction. The combined analysis of
\textsc{InIce}\xspace{}/\textsc{IceTop}\xspace{}/RASTA events
will permit cross-calibration of all detection techniques and improved composition measurements.
The proposed three year development plan will allow to establish and test all
the key ingredients required towards reaching these goals.

\bibliographystyle{model1a-num-names}


\section{Introduction}
Cancer progression involves a complex interaction of genetic and
genomic factors that jointly subvert normal cell development.  The~genomic component, which encompasses gene expression and regulation,
is substantially impacted by the biochemical composition of the local
environment in which a cell grows. So-called \textit{micro-environmental}
parameters, including levels of oxidation, lactate, acidity,
nutrients of various kinds and other factors affecting physical
interactions between cells, are increasingly studied for their
potential to improve our understanding of cancer biology, and for
their promise to lead to new therapeutic strategies. Changes in such
parameters can impact gene transcription, which in turn impacts
protein production.  Variation in these fundamental parameters can
therefore induce a cascade of effects, producing disruptions of normal
cellular processes in downstream biological pathways [\citet{hanahan}].
For example, changes to the pH level in the cellular environment may
effect glycolysis, thus impacting on numerous genes involved in the
glycolysis pathway.  Some of these genes may also play roles in the
regulation of cell growth, and their suppression may engender
tumorigenesis and promote the aggressive advance of existing cancerous
states.  Microarray gene expression assays can be used to generate
data on the transcriptional response of cancer cells to controlled
manipulations of environmental factors such as pH.  This data is
useful for characterizing these \textit{micro-environmental response pathways}.

Our study concerns changes to cellular pH levels, and the resulting
\textit{neutralization} and \textit{lactic acidosis response pathways}.  Section
\ref{signaturesection} describes the application of sparse Bayesian
regression models [\citet{lucas06};\break \citet{seowest07}] to microarray data
generated through a series of laboratory experiments on cultured
breast tumor cells in which cellular pH levels were manipulated in a
controlled manner.  These analyses yield statistical \textit{expression
  signatures} of the cellular responses to various interventions on
the pH level.  The~main challenge lies in relating these signatures,
and the biological pathways they characterize, to variation in gene
expression across large samples of human breast tumors.  This
integration of in vitro and in vivo data sets is the
driving focus of this and related studies. In addition to comprising a
detailed study of new data and experimental results, through which are
generated several directions for biomedical research, this work
exemplifies an overall strategy for cross-study, integrative analysis
of gene expression data for exploring and relating pathway-related
experimental findings to clinical contexts and patient outcomes.

When considering variability in expression patterns of genes in
observational tumor data, we face questions of differences due to the
differing contexts. It is to be expected that a tumor in vivo
evidences far more complex and heterogeneous biological variation than
in the controlled in vitro setting, and this will be manifest in
measures of gene expression. Normal cell processes held in quiescence
in cell cultures may when active co-regulate the expression of
relevant signature genes in in vivo, confounding the pattern of
expression that was evident in vitro.  Hence, when aiming to
translate experimental findings to tumor populations, thereby
providing a mapping of an in vitro signature to its in vivo counterpart, we require statistical models capable of
discovering and representing the additional complexity surrounding and
interacting with the original response signature.  Section
\ref{statsection} describes our analysis of a large and heterogeneous
breast cancer data set using sparse latent factor models
[\citet{west2003bayes7}; \citet{carvalho2008}] that satisfy these
desiderata. This analysis includes a targeted factor search that
facilitates estimation of statistical factors associated with an
initial set of genes underlying the in vitro experimental
signatures.  The~factors discovered in this way represent a modular
decomposition of the biological patterns evident in the in vivo
breast cancer data, while retaining connections to the experimental
signatures.

Section \ref{biosurvsection} discusses aspects of the biological and
clinical interpretations of these estimated factors, which can be
viewed as a refined in vivo set of summary biomarkers of
variation in the \textit{neutralization} and \textit{lactic acidosis
  response pathways} of these breast cancers. In survival analyses, we
find that these factor model derived biomarkers have substantial
prognostic value in connection with long-term survival and, hence, the
sets of genes comprising these factors warrant further study.  We
present predictive validation of this key finding in analyses of two
separate breast cancer data sets.  We then provide biological
interpretation of one key factor that emerged from the evolutionary
factor search, which plays a key role as a predictive variable in the
survival analyses.  It turns out that this factor is a single component
of a specific biological pathway that has previously been noted as a
risk biomarker in cancer, but not, to date, connected at all into
response pathways linked with variation in cellular pH.  This finding
has generated follow-on biological research and initiated a new line
of experimentation on the role of this pathway in connection with
cancer cell micro-environmental influences.

\section{Neutralization experiments and analysis}\label{signaturesection}

\subsection{Biological and experimental context}
Investigating the effects of\break changes in the micro-environment in which
cells grow is of increasing interest in cancer research.  The~tumor
micro-environment is typically characterized by oxygen depletion, high
lactate and extracellular acidosis coupled with vascular leakage,
glucose and energy deprivation. These and other micro-environmental
features vary widely across tumors and generally exhibit substantial
temporal and spatial differences in a tumor. Micro-environmental
stresses trigger biochemical changes in cancer cells that directly
modulate physiological, metabolic and ultimately clinical
phenotypes. Improved understanding of the molecular mechanisms of such
tumor responses holds promise for immediate translational impact and
clinical care, as relevant therapies can be brought to bear to modify
the micro-environment.

Currently, with the exception of hypoxia, very little is understood
about how each individual stress affects cellular phenotypes and tumor
progression.  To examine how cancer cells respond to increased acidity
or pH neutralization at different time points, MCF7 cell cultures (a
commonly-used breast tumor cell line) were grown in neutral media and
then exposed to varying interventions in several assays in
parallel. For some cells, lactic acid was added to the medium (25 mM
lactic acid at pH 6.7) for 1 and 4 hours; others cells experienced
strong lactic acidosis conditions (25 nM lactic acid at pH 5.5) for 4
hours.  Similarly, the effects of neutralization were assayed by
shifting the MCF7 cultures from overnight lactic acidosis conditions
at pH 6.7 to neutral regular media at pH 7.4 for 1 and 4 hours.
Control cells were grown in each starting condition (neutral
conditions and lactic acidosis conditions).  The~complete set of
experiments is summarized in Table \ref{exptable}.  The~mRNA extracted
from each of the resulting $n=27$ batches of MCF7 cultures was
purified using Ambion miRVana RNA purification kits and standard
microarray assays were performed using Affymetrix U133 Plus 2 Genechip
platforms. All raw microarray data were preprocessed using RMA
[\citet{rma}], the log (base 2) scale output of which were used in all ensuing
statistical analyses.

\begin{table
\caption{Summary of neutralization/acidosis experiments.
Cell entries indicate the number of replicates per experimental group}\label{exptable}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lccccc@{}}
 \hline
  & \multicolumn{5}{c@{}}{\textbf{Exposure condition}}\\[-6pt]
 & \multicolumn{5}{c@{}}{\hrulefill}\\
 & \multicolumn{2}{c}{\textbf{pH 7.4}}& \multicolumn{2}{c}{\textbf{pH 6.7}} & \multicolumn{1}{c@{}}{\textbf{pH 5.5}} \\[-6pt]
 & \multicolumn{2}{c}{\hrulefill}& \multicolumn{2}{c}{\hrulefill} & \multicolumn{1}{c@{}}{\hrulefill} \\
 \textbf{Growth condition} & \textbf{1hr} & \textbf{4hr} & \textbf{1hr} & \textbf{4hr} & \textbf{4hr}\\
\hline
pH 7.4 & 3$\times $ & 3$\times $ & 3$\times $ & 3$\times $ & 3$\times $\\
pH 6.7 & 3$\times $ & 3$\times $ & 3$\times $ & 3$\times $ &\\
\hline
\end{tabular*}
\end{table}

\subsection{Cellular response signatures}\label{sparsereg}
Quantitative summaries of the cellular responses to lactic acidosis
and neutralization treatments were obtained using a standard sparse
multivariate regression model [\citet{lucas06}; \citet{seowest07}].  We analyzed
$19{,}375$ genes (technically, probe-sets from the Affymetrix array; we
will use ``gene'' and ``probe'' interchangeably) whose median
expression level is at least 5.5 and whose expression ranges more than
$0.5$-fold across the $n=27$ experimental samples.  Let $X^{\mathrm{exp}}$
denote the $19{,}375 \times 27$ matrix of expression values.  Rows
represent genes and columns correspond to three replicate samples for
each of the following experimental groups: (i) control (pH~$7.4\rightarrow 7.4$) at 1~hour; (ii) control at 4 hours; (iii)
lactic acidosis (pH $7.4\rightarrow 6.7$) at 1~hour;  (iv) lactic acidosis at 4 hours; (v) neutralization
(pH~$6.7\rightarrow 7.4$) at 1~hour; (vi) neutralization at 4 hours;
(vii) acidic growth (constant pH of 6.7) at 1~hour; (viii)
acidic growth at 4 hours; (ix) strong lactic acidosis (pH
$7.4\rightarrow 5.5$) at 4 hours.  Let $H^{\mathrm{exp}}$ denote the $11 \times
27$ design matrix where the first 8 rows contain binary indicators for
effects associated with differential expression relative to the 1hr
control group: 1hr lactic acidosis effect, 1hr neutralization effect,
1hr acidic growth effect, 4hr control effect, 4hr lactic acidosis
effect (relative to 4hr control), 4hr neutralization (relative to 4hr
control), 4hr acidic growth effect (relative to 4hr control) and 4hr
strong lactic acidosis effect (relative to 4hr control).  The~last
three rows contain artefact control factors derived from the first
three principle components of the expression levels associated with
the AFFX series control genes included on the Affymetrix microarrays.
These control genes are not variably expressed in humans, and so
patterns of variation across samples manifest in control genes
represents systematic errors arising from different experimental
conditions.  Use of these artefact control factors provides
opportunity for sample-specific correction of artefactual effects on
genes that may otherwise result in false-discovery or obscure
meaningful biological variation
[following \citet{lucas06} and \citet{carvalho2008}].  After deriving the
artefact control factors, rows corresponding to Affymetrix control
genes are removed from subsequent analyses.

The~model for the
expression of gene $g$ in sample $i$ is
\[
x_{g,i}^{\mathrm{exp}}  =  \mu_g + \sum_{k=1}^{11}\beta_{g,k}h_{k,i}^{\mathrm{exp}} +
\nu_{g,i}
\]
or in matrix from
\[
X^{\mathrm{exp}}  =  \mu 1 + BH + N,
\]
where  $\mu_g$ denotes the mean expression of
gene $g$ in the 1hr control samples,
each $\beta_{g,k}$ is the change in
expression of gene $g$ due to design factor  $k$, and the
$\nu_{g,i}$ are independent,  normally distributed idiosyncratic noise terms representing
residual biological variation, experimental and measurement errors with individual
variances $\psi_g.$
Sparsity is induced via prior distributions that place
positive probability on $\beta_{g,k}=0$ for each $g,k$ pair, and
resulting posterior analysis allows investigation of posterior
sparsity patterns via probabilities  $\pi_{g,k}^* = \operatorname{Pr}(\beta_{g,k}\ne 0|X^{\mathrm{exp}}).$
Full details follow \citet{lucas06} and prior specifications, including priors
for the~$\mu_k,$ variance parameters and all hyper-parameters, are given in
 \ref{appendix}.  Posterior inference via MCMC is achieved using the BFRM software [\citet{bfrm}].

Figure \ref{sigskeleton} broadly illustrates genes uniquely associated
with individual treatment effects as well as those involved in
multiple responses.  This gives some indication of the degree of
intersection of the cellular pathways being queried by the different
treatments.  Across the 8 treatments, the sparsity, as measured by the
percent of genes for which $\pi_{g,k}^*>0.99$, ranges from $29\%$ (4~hour neutralization) to $46\%$ (4~hour lactic acidosis).  The~fold-change associated with the involved genes ($2^{|\beta_{g,k}|}$
for $g$ such that $\pi_{g,k}^*>0.99$) ranges from 1.06$\times $ to 13$\times $, with a
mean of 1.4$\times $.

\begin{figure

\includegraphics{261f01.eps}

     \caption{Neutralization signature skeleton: black indicates genes $g$ (rows)
    with posterior probability $\pi_{g,k}^*>0.99$ for each experimental
    group $k$ (columns). Genes are ordered to emphasize which genes are unique
    to each successive experiment relative to the previous.}\label{sigskeleton}
\end{figure}

\begin{figure

\includegraphics{261f02.eps}

  \caption{Neutralization treatment
        signature scores ($s_{k,i}$) for each sample in the original
        study. Separate treatment groups are color coded.} \label{sigcohorts}
\end{figure}

The~cellular response to each treatment, also called the \textit{signature} of the treatment, is characterized by estimated effects
$\beta_{g,k}^* = E(\beta_{g,k} | \beta_{g,k}\ne 0,X^{\mathrm{exp}})$ together
with the $\pi_{g,k}^*.$ The~ability of each signature to uniquely
identify the treatment it reflects can be further explored using
summary \textit{signature scores} as defined in \citet{lucas09}.  Based on
posterior means $\beta_{g,k}^*$ and $\psi_g^*,$ let
\[
s_{k,i} = \sum_{g=1}^{19{,}375}
\beta_{g,k}^*x_{g,i}^{\mathrm{exp}}/{\psi_g^*}
\]
define the score for treatment signature $k$ on sample $i.$ This
expression is derived from the data-driven component of the Bayes
factor that weighs the evidence in favor of the given signature
describing the variation in a sample ($p(x_i | h_{k,i}=1)/p(x_i |
h_{k,i}=0)$).  Figure \ref{sigcohorts} shows the values of the scores
associated with 7 treatment signatures plotted across samples.  As
expected, the highest scoring samples for each signature are those
upon which that signature is based, but important connections between
signatures can be identified on the basis of other high- or
low-scoring treatment groups.  For example, there is a inverse
relationship between the 4hr acidosis score and the 4hr neutralization
score.  Also evident is the similarity between the 1hr and 4hr acidic
growth signatures, which can also be inferred through the large
intersection of the genes defining the two signatures (Figure
\ref{sigskeleton}).

\section{Latent factor analysis of breast tumor gene expression}\label{statsection}

\subsection{In vivo  breast cancer data}\label{millersection}
The~primary goals of this study are to uncover shared structures in
the cell response signatures defined above, and to quantify the extent
to which these structures can be used to predict clinical phenotypes
in real human cancers.  Here we make use of the gene expression data
for a collection of 251 surgically removed breast tumors as reported
in \citet{miller}.  Affymetrix 133A and 133B
GeneChip microarrays were generated for each tumor sample, and
relevant clinico-pathological variables were collected for each
patient.  This included age at diagnosis, tumor size, lymph node
status (an indicator of metastatic cancer) and Elston histological
grade (a categorical rating of malignancy as deemed by pathologists).
Molecular assays to identify the presence of absence of mutations in
the estrogen receptor (ER), progesterone receptor (PgR) and P53 genes
were also performed. These data are representative of a variety of
different presentations of human breast cancer on these clinical
measures.

We first evaluate the signature score as defined above for each tumor.  These
scores are then standardized across samples so that each vector of 27
scores for a particular signature has mean and variance equal to the
mean gene-specific expression and mean gene-specific variance.  This
transformation places the signature scores on the same scale as gene
expression in the tumor data set, thus enabling a ``metagene''
interpretation of a vector of scores [\citet{west01}; \citet{pittman04}]; see
Figure \ref{sigscores}.
The~relationships between the tumor signature scores bear some
similarities to those observed in the cell line study.  There is once
again evidence of correlation between the two acidic growth signatures
and the 1hr neutralization signature.  These three signatures, in
turn, display patterns opposite that of the 4hr neutralization
signature.  The~patterns are less prominent, however, than was evident
in the cell culture data.  Although the variation in these scores
presumably relates, in part, to underlying biological variation in the
activity of the lactic acidosis and neutralization response pathways
within these tumors, as mentioned above, the set of genes
characterizing the in vivo effects of lactic acidosis and
neutralization may differ substantially from those characterizing the
in vitro responses as a result of the more complex interactions
with other cellular processes.

\begin{figure

\includegraphics{261f03.eps}

     \caption{Initial evaluation of neutralization signature levels
      across tumor samples.  Samples are ordered by first principle
      component to emphasize dominant signature
      gradients.}\label{sigscores}
      \end{figure}

We thus aim to refine our evaluation of the response pathway activity
levels in the tumor data by using the signature scores as initial
``anchors'' in an analysis using sparse latent factor models.  The~main idea is to define statistical factors on sets of genes related to
these initial scores, and to link in other genes that may connect with
the different response pathways active in vivo.  This is
accomplished as follows.

\subsection{Sparse factor model specification}\label{sparsefac}
Sparse latent factor models represent common patterns in gene
expression via latent factors in which the factor-gene relationships
are sparse; this notion of statistical sparsity\break is key for
representing the intersecting subsets of genes potentially\break related to
underlying networks of biological pathways
[\citet{west2003bayes7};\break \citet{seowest07}; \citet{lucas08};\break \citet{carvalho2008}].  The~form of
the statistical model is an extension of the sparse regression
model. A key part of our analysis strategy stems from augmenting the
$44{,}592 \times 251$ matrix of gene expression data for the tumor data
with the 7 values of the projected treatment signature scores.  Let
$p=44{,}592+7=44{,}599$ and $n=251$, and let $X^{\mathrm{obs}}$ denote the $p
\times n$ matrix in which the first 7 rows are the projected scores
across tumor samples, and rows $8-p$ are the gene expression values.
Here we will make use of $K=4$ artefact control factors derived from
the first four principal components of the control genes of the breast
tumor microarrays.  A latent factor model consisting of $L$ latent
factors is therefore
\[
x_{g,i}^{\mathrm{obs}}  =  \mu_g + \sum_{k=1}^K \alpha_{g,k}\lambda_{k,i} + \sum_{l=K+1}^{K+L}\alpha_{g,l}\lambda_{l,i} + \nu_{g,i}
\]
or, in matrix form,
\[
X^{\mathrm{obs}}  =  \mu 1 + A \Lambda + N,
\]
where: (i) the first $K$ rows of the $(K+L) \times n$ matrix
$\Lambda$ are the known artefact controls; (ii) the remaining
$L$ rows contain latent factor scores; (iii) the first $K$
columns of the $p\times (K+L)$ matrix $A$ are regression parameters on
the artefact controls (changing notation from the earlier $\beta$ to
$\alpha$ for notational convenience here); (iv) the remaining
$L$ columns of $A$ are factor loadings parameters relating factors to
genes and to the projected scores; and (v) $A$ is sparse, with
sparsity pattern to be inferred along with estimation of nonzero
values.  The~model is completed by assigning sparsity priors over
columns of $A,$ precisely as was done for $B$ in the sparse regression
model; prior specification for $A,$ variance components and other
hyper-parameters follows default recommendations for the BFRM
framework (see \ref{appendix}).

Flexibility in representing potentially complicated patterns
underlying expression is achieved using nonparametric Bayesian
Dirichlet process models for the factor scores.  The~$L$-vectors $
(\lambda_{K+1,i},\ldots,\lambda_{K+L,i})'$, representing the latent
factor values on tumor sample $i$, are modeled as draws from an
unknown latent factor distribution subject to a Dirichlet process
prior with a multivariate normal base measure.  This standard
nonparametric mixture model allows great flexibility in adapting to
nonnormal structures commonly manifest in factor scores
[\citet{carvalho2008}; \citet{bfrm}].

Ensuring the identifiability of latent factors requires the use of a
modified prior on $A$ such that the leading L rows have an upper
triangle of zeros and positive upper diagonal elements; that is, for
$g=1\dvtx L,$ we have $\alpha_{g,g+K}>0$ and $\alpha_{g,l}=0$ for $l>g+K$.
The~first $L$ variables in $X^{\mathrm{obs}}$ then represent ``founders''
of the $L$ latent factors, with variable $g$ associated with a
$\alpha_{g,g}$-fold change in expression due to factor $g,$
$(g=1,\ldots,L).$ It also defines an hierarchical dependence on the
factors, namely,
\begin{eqnarray*}
x_{1,i}^{\mathrm{obs}} & = & \cdots + \alpha_{1,K+1}\lambda_{K+1,i} + \nu_{1,i}, \nonumber\\
x_{2,i}^{\mathrm{obs}} & = & \cdots + \alpha_{2,K+1}\lambda_{K+1,i} + \alpha_{2,K+2}\lambda_{K+2,i}+\nu_{2,i},\nonumber
\label{facreg}\\
x_{3,i}^{\mathrm{obs}} & = & \cdots + \alpha_{3,K+1}\lambda_{K+1,i} + \alpha_{3,K+2}\lambda_{K+2,i}+\alpha_{3,K+3}\lambda_{K+3,i}+\nu_{3,i}\nonumber
\end{eqnarray*}
and so on.  This structure aids the interpretation of the latent
factor loadings as representing interconnected components of a complex
biological process.  The~latent factor scores $\lambda_{i,l}$ quantify
variation across tumors for these expanding levels of complexity, with
each additional factor accounting for variation in observed gene
expression unaccounted for by the previous set of factors.  With our
use of projected in vitro signature scores here as the first 7
variables, the first 7 factors will now represent patterns underlying
co-variation in expression of sets of genes that link indirectly to
these treatment signatures. Additional factors then reflect other
dimensions of common variation in the set of genes analyzed.

\subsection{Targeted factor search}\label{facmodsearch}
Decomposition of the patterns of variation evident in the tumor gene
expression data into latent factors proceeds through evolutionary
model search, full details of which appear in
\citet{carvalho2008} and \citet{bfrm}.  The~evolutionary model search provides a
computationally efficient approximation to the computationally
prohibitive full factor analysis on the entire set of genes, and
produces full posterior results for the final set of factors and
genes.  A key novelty of this approach is that we exploit the
sensitivity of the model search procedure to its initial configuration
in order to explore the space of factor models surrounding an initial
model containing $7$ latent factors and representing only the 7
response metagenes.  By construction, these initial factors are each
defined, or ``founded,'' by the neutralization/lactic acidosis
treatment scores, thereby ensuring that the model search is primarily
concerned with patterns of variation related to these particular
response pathways.

Evolution of this initial model proceeds as follows.  Samples from the
joint posterior distribution of model parameters are obtained through
MCMC.  Based on these fitted values, we impute inferences for all
genes $g>7$ that are not currently included in the model, as described
in \citet{carvalho2008}.  Thus, after fitting the initial factor model
which considers only the signature scores, we examine expression
levels of the full set of $44{,}000+$ genes for evidence of association
with the current factors.  The~imputation process generates
approximate probabilities $\pi_{g,l}^* = \operatorname{Pr}(\alpha_{g,l}\ne 0
| X^{\mathrm{obs}})$ for all such genes $g.$ Genes are ranked on the basis of
these probabilities, and the model is then expanded to include a small
number of the genes with largest values of the projected
$\pi_{g,l}^*.$ The~model is then refitted to this expanded sample, and
if appropriate, the number of factors is increased in order to adapt
to additional common patterns of expression variation now evident in
the increased set of variables being modeled.  This process is
repeated until no new genes or factors can be added, or until the
model reaches a designated maximum size.  More details on the search
strategy, including control parameters governing model expansion, are
given in \ref{appendix}.

\begin{figure

\includegraphics{261f04.eps}

     \caption{Skeleton of fitted factor loadings for  tumor data.  Black indicates
      variable-factor loadings with $\pi_{g,l}^*>0.99.$ The~first 7 variables are the
       projected neutralization scores, followed by 493 genes
       reordered for a clear visual presentation of the sparsity
       structure of, and cross-talk in, gene-factor loadings.}\label{facskel}
\end{figure}

The~initial 7-gene, 7-factor model evolved under this process to reach
a terminal size of 500 variables (the designated maximum)
incorporating 30 latent factors.  Figure \ref{facskel} shows the
skeleton of the factor structure, in terms of major patterns of
gene-factor relationships.  The~ordering of the factors is determined
by the model search procedure, and represents the incremental
improvement to model fit provided by each subsequent factor.  In this
sense, each subsequent factor builds upon the complexity modeled by
the previous factors.  The~leading 7 factors
correspond to the following signatures, respectively: 4hr lactic
acidosis, 1hr lactic acidosis, 1hr neutralization, 4hr strong lactic
acidosis, 4hr acidic growth, 1hr acidic growth, and 4hr
neutralization.  Like their in vitro signature counterparts, the
in vivo factors loadings contain a great deal of sparsity.  Of the 493
genes included in the final model, only 333 are among those identified
in the in vitro signature analysis.  Factor 1, founded by the
the 4hr lactic acidosis signature score, has 173 genes with nonzero
loadings at the $0.99$ probability threshold, compared to $8909$ in
the in vitro signature.

\begin{figure

\includegraphics{261f05.eps}

 \caption{Heat map showing the magnitudes of the fitted factor loadings $\alpha_{g,l}^*$
for the first $L=30$ rows of $A$. The~founder gene for each factor is designated by its
U133$+$ probe ID.  The~terms in each row determine a linear combination
of latent factors that predict the observed expression levels of the founder gene. }\label{crosstalk}
  \end{figure}

Posterior estimates of the factor loadings ($\alpha_{g,l}^* =
E(\alpha_{g,l}|\alpha_{g,l}\ne 0,X^{\mathrm{obs}})$) aid in generating further
insights. In particular, the upper portion of the estimated loadings
matrix sheds light on the structure of connections between latent
factors; see Figure \ref{crosstalk}.  As described in Section
\ref{sparsefac}, one interpretation of a row $A$ is as a set of
coefficients determining a linear combination of factor scores that
predict the gene expression vector for the corresponding variable.
The~inset of Figure \ref{crosstalk} shows that the fitted values of
all 7 of signature scores involve positive contributions from factor
1, the factor version of the 4hr lactic acidosis signature.  Thus, the
pattern of 4hr lactic acidosis signature activity across samples
describes a fundamental pattern of pathway activation that underlies
the activity patterns of the other 6 signatures.  The~seventh factor
(i.e., the factor representation of the 4hr neutralization signature)
sits atop this hierarchy of pathway complexity, represented as a
linear combination of factors 1 (4hr neutralization), 3 (1hr
neutralization), 4 (4hr strong lactic acidosis) and 5 (4hr acidic
growth), plus the additional pattern of expression unique to this
pathway.

\section{Biomedical connections of factor profiles}\label{biosurvsection}

\subsection{Factor-based prediction  of long term survival}\label{survivalsection}
The~in vivo latent factors linked to neutralization pathways
represent complexity in the patterns of expression, and therefore in
the levels of underlying biological pathway activation evident across
the tumor samples.  For this reason, latent factors can be regarded as
candidate biomarkers of physiological states that link to these
pathways. Our study explores this using the posterior mean factor
scores $\lambda_{l,i}^*$ as candidate predictors in a survival analysis
of the breast cancer patient data.

We use Weibull regression models of patient survival that draw on the
30 estimated neutralization/lactic acidosis pathway factors, the 7
original projected signature scores and the clinical covariates
available for this data set [\citet{miller}].  The~latter include
histologic grade, ER mutation status, node status, P53 mutation
status, PgR mutation status, tumor size and age at diagnosis.  This
analysis allows both integration and comparison of the prognostic
value of these traditional markers with specific pathway-related
signature scores, and their latent factor representations---an
integrative clinico-genomic analysis.  Let $t_i$ denote the survival
time of patient $i.$ The~Weibull density function is $p(t_i|a,\gamma)
= at_i^{a-1}\exp(\eta_i - t_i^ae^{\eta_i}),$ where $a>0$ is the index
parameter and $\eta_i = \gamma'y_i$ the linear predictor based on a
covariate vector $y_i.$ We explore subsets of covariates and
regression model uncertainty using Bayesian shotgun stochastic search
[\citet{hans2007}; \citet{sss}]. This generates a list of regression covariate
subsets and the corresponding posterior regression model probabilities
for use in Bayesian model averaging for survival prediction and in
exploring relevance of variables. Details of model and prior
specification follow defaults in the SSS software [\citet{sss}] as noted
in \ref{appendix}.

Figure \ref{incprobs} shows posterior probabilities for each of the
46 candidate covariates.  Nodal status emerges as the leading
predictor of long term survival, followed by latent factor 30 and then
tumor size.  Note that none of the original signature scores, and no
other clinical variables, receive appreciable probability.  That nodal
status provides the best predictor of survival is to be expected.
Previous studies [\citet{west01}; \citet{pittman04}; \citet{dressman06}] have shown that
nodal status is not well predicted by gene expression and that
combined use of nodal status with gene expression predictors can
improve survival prognosis. Hence, it seems that the
information content of the nodal status predictor is unlikely to
overlap with that of any factor score.  This can also be clearly seen
through the pairwise inclusion probabilities in Figure
\ref{pairprobs}.

\begin{figure}[t]

\includegraphics{261f06.eps}

     \caption{Posterior inclusion probabilities for the 46 candidate covariates in the Weibull
      models. The~candidate covariates include the 30 estimated
      latent factors, followed by the 7 original signature scores,
      followed by 10 traditional clinical covariates
         (Elston grades 1, 2 and 3, ER, node
      status, P53, PgR status, tumor size, age at diagnosis).}\label{incprobs}
  \end{figure}

\begin{figure}[b]

\includegraphics{261f07.eps}

     \caption{Pairwise  inclusion probabilities for the top
      6 predictor variables in breast cancer survival analysis.
      Darker colored tiles indicate higher
      probabilities.}\label{pairprobs}
  \end{figure}

The~pairwise inclusion probability of factor 30 and nodal status is
close to the marginal probability of
factor 30; however, the pairwise inclusion probability of factor 30
with any of the other factors is far less than any of the marginal
inclusions probabilities of those factors; thus, factor 30 is clearly a dominant and
preferred expression-based biomarker of survival risk.

\begin{figure}[t]

\includegraphics{261f08.eps}

     \caption{\textup{(a)} Kaplan--Meier curves demonstrating stratification
     of Miller data into high- and low- risk groups, based on the fitted Weibull mixture.
     \textup{(b)} Estimated survival curve associated with varying levels of
     factor 30, holding all other predictors at their median values.
     \textup{(c)} Kaplan--Meier curves demonstrating stratification of
     Miller data into high- and low- risk groups, based solely on the value of factor 30.}\label{survplots}
     \end{figure}

Posterior and predictive inferences are formally based on a mixture of
1000 Weibull survival models, mixed with respect to their posterior
probabilities.  For each patient $i,$ we can compute the implied
predicted survival function at her covariate vector $y_i$ and identify
the predicted median survival $m_i;$ this is the predicted median
survival time for a \textit{future} patient who has the same covariates
as patient $i.$ Figure \ref{survplots}(a) shows a Kaplan--Meier
survival plot of the Miller et al. data simply stratified by
$m_i\le m$ or $m_i>m,$ where $m = \operatorname{median}\{ m_i,
i=1,\ldots,n\}.$ There is approximately a 30\% difference in the
empirical 10-year survival probability between patients cohorts
stratified crudely on this basis, as a simple visual of the relevance
of the included covariates.  By way of focusing on factor 30, we plot
the model-averaged survival curve for a hypothetical patient whose
covariates are held constant at their median values in the data set,
save for variation in the factor~30~score; see Figure
\ref{survplots}(b), where factor 30 is set at its $5$th, $50$th
and $95$th percentiles in the data set, all other covariates
remaining fixed.  The~estimated effect of variation in factor 30 alone
accounts for approximately 20\% of the difference in 10-year patient
survival between the high-risk and low-risk subgroups.  This
prediction is confirmed by considering the Kaplan--Meier curves formed
by stratifying the patients on the basis of the patient-specific
factor 30 value compared to the median across samples; see Figure
\ref{survplots}(c).  The~pattern of gene expression comprising
the loading associated with factor 30 warrants further investigation,
to which we will return in Section \ref{biof30}.

\begin{figure}[b]

\includegraphics{261f09.eps}

 \caption{Kaplan--Meier curves demonstrating stratification of
\textup{(a)} the \protect\citet{pawitan05} patient samples, and
\textup{(b)} the \protect\citet{sotiriou06} patient samples (right)
  into high- and low- risk groups based on imputed values of factor
  30, as identified in the latent factor analysis of the Miller data.}
  \label{bcsurvpredict}
  \end{figure}

\subsection{Out-of-sample factor projection}\label{facproj}
It is critical to evaluate whether or not the above results can be
confirmed through out-of-sample prediction.  We do this with two
additional breast cancer data sets: that of \citet{pawitan05}, consisting of 159 primary breast tumors assayed on
Affymetrix U133A and U133B chips, and that of \citet{sotiriou06}, consisting of 189 primary breast tumors assayed on
U133A chips.

Fixing all factor model parameters at their posterior means, we can
directly predict values of the latent factors for each new patient;
see\break \ref{appendix} and [\citet{lucas08}].  Note that this calculation is purely
predictive; no model fitting nor additional analysis of the two
validation data sets was performed.  Using the predicted latent factor
vectors, we can produce the same survival plots for these data,
stratifying each of the two new patient cohorts on the basis of their
factor 30 scores as above; see Figure \ref{bcsurvpredict}, as compared
to Figure \ref{survplots}(c).
The~association between low factor 30 scores and good prognosis
remains evident in these out-of-sample predictions that draw on
different patient populations.  Further, the differences between high
and low risk groups is comparable across all three samples.

The~robustness of factor 30 as a prognostic biomarker provides strong
support for the view that it reflects a biologically meaningful module
of gene expression.  By evaluating the predicted factor scores in the
original experimental data, we are able to establish that factor 30,
despite its relatively late incorporation to the model, is linked to
the 4-hour lactic acidosis pathway.  Figure \ref{f1f30reverse}
compares the predicted factor scores for factors 1 and 30 as evaluated
in the experimental data.  Factor 1, which is founded by the 4-hour
lactic acidosis signature, is in fact comparable to the original
signature score as depicted in Figure \ref{sigcohorts}.  Factor 30 has
its highest values in the original 4-hour lactic acidosis samples, but
shows a different pattern of activity across the other sample cohorts.
In particular, factor 30 appears to be repressed in the samples
associated with the 4-hour neutralization and acidic growth
treatments.  This implies that factor 30 may characterize some
critical intersection between these pathways that is itself related to
tumor aggressiveness.

\begin{figure}[t]

\includegraphics{261f10.eps}

\caption{Predicted values of factors 1 and 30 for each sample in the original experimental study.
Separate treatment groups are color coded.}\label{f1f30reverse}
\end{figure}

\eject
\subsection{Biological evaluation of prognostic factor 30}\label{biof30}

Having established the clinical relevance of factor 30, the task
remains to ascribe to it biological meaning. The~loading of factor 30
is comprised of only 6 gene probe sets for which
$\pi_{g,l}^*>0.99$. Four of these, including the founder gene, are
related to the phosphoglycerate kinase 1 (PGK1) gene, while the other
two are related to a neuronal cell death-related protein and the CEGP1
protein.  The~factor is characterized by overexpression
($\beta_{g,k}>0$) of the PGK1 and neuronal cell death proteins and
suppression ($\beta_{g,k}<0$) of CEGP1.

A literature search generates detailed biological information on PGK1,
and its role in the glycolysis pathway where it is fundamental to
cell\break
growth and metabolism.  PGK1 catalyzes the reversible conversion
of\break
1,3-diphosphoglycerate to 3-phosphoglycerate with the generation of
one molecule of ATP and this represents an important step in
glycolysis pathways. In addition, PGK1 has been reported to induce
other processes related to cancer progression, such as conferring a
multi-drug resistant (MDR) phenotype [\citet{duan}] and affecting tumor
angiogenesis through affecting secreted plasmin [\citet{lay}]. Previous
studies have also shown that elevated levels of PGK1 predict \textit{poor} survival outcomes in lung cancers [\citet{chen03}], and that
PGK1 can often be expressed at high levels in pancreatic [\citet{hwang}]
and renal [\citet{unwin}] cancers.  The~association between high factor
30 levels and poor prognosis indicates a similar relationship between
PGK1 and survival may exist for breast cancers.

Since PGK1 is an important component of glycolysis pathways, our
findings here may implicate glycolysis activities in poor patient
survival. This is supported by previous findings that expression of
glycolysis pathways and PGK1 are repressed by lactic acidosis
[\citet{chen08}]. Factor 30 links the neutralization pathway response
signatures to a clear \textit{PGK1 factor} that may now serve as a
biomarker of one key aspect of tumor responses to changes in pH with
the potential to aid in predicting follow-on changes in tumor
metabolism via glycolysis pathway activation.  Further evaluation of
this chain of relationships is now initiated and will be explored
using independent methods such as tumor tissue microarrays
[\citet{chen03}]. Since PGK1 and glycolysis pathways are also controlled
by hypoxia [\citet{chi}], these results also highlight their potential
roles as integral mediators of multiple micro-environmental factors
affecting tumor progression and clinical outcomes.

\section*{Acknowledgments}
The~authors are grateful to an editor of AOAS and an anonymous
reviewer for constructive comments on this manuscript.

\eject
\begin{supplement}[id=supp:code+etc]
\sname{Supplement A}
\stitle{Software and Data}
\slink[doi]{10.1214/09-AOAS261SUPPA}
\sdatatype{url}
\slink[url]{http://ftp.stat.duke.edu/WorkingPapers/08-34.html}
\sdescription{
  This site contains all materials needed to reproduce the reported
  analyses.  This includes all data files, control files for the BFRM
  and SSS software, and MATLAB functions for producing graphical summaries.}
\end{supplement}

\begin{supplement}[id=appendix]
\sname{Supplement B}
\stitle{Appendix}
\slink[doi]{10.1214/09-AOAS261SUPPB}
\slink[url]{http://lib.stat.cmu.edu/aoas/261/supplement.pdf}
\sdatatype{.pdf}
\sdescription{The appendix Merl et al. (\citeyear{merl09}) provides further details
on prior specifications in the sparse regression and sparse latent factor models.
The~appendix also contains details on the control parameters for the evolutionary
factor search and shotgun stochastic search, and describes the procedure for imputing factor scores in new samples.\looseness=-1}
\end{supplement}
\iffalse

\fi


\section{Introduction}
 Let $G=(V,E)$ be a simple graph of order $n=|V|$.
Let $u,v\in V$ be two different vertices of $G$, the distance $d(u,v)$ between vertices $u$ and $v$ is the length
of the shortest path between $u$ and $v$. Given a set of vertices $S=\{v_1,v_2,...,v_k\}$ of $G$,
the {\it metric representation} of a vertex $v\in V$ with respect to $S$ is the vector $r(v|S)=(d(v,v_1),d(v,v_2),...,d(v,v_k))$.
We say that $S$ is a {\it resolving set} for $G$ if for every pair of different vertices $u,v\in V$, $r(u|S)\ne r(v|S)$.
The {\it metric dimension} of $G$ is the minimum cardinality of any resolving set for $G$ and it is denoted by $dim(G)$.
The concept of metric dimension was introduced first independently by Harary and Melter \cite{harary} and Slater \cite{leaves-trees}, respectively.

Let $G$ and $H$ be two graphs of order $n$ and $m$, respectively. The corona
product $G\odot H$ is defined as the graph obtained from $G$ and $H$ by taking one copy of $G$ and $n$ copies of $H$ and
then joining by edges all the vertices from the $i^{th}$-copy of $H$ with the $i^{th}$-vertex of $G$.

Given the graphs $G$ and $H$ with set of vertices
$V_1=\{v_1,v_2,...,v_{n}\}$ and $V_2=\{u_1,u_2,...,u_{m}\}$,
respectively, the Cartesian product of $G$ and $H$ is the graph
$G\times H$ formed by the vertices $V=\{(v_i,u_j)\;:\;1\le i\le
n,\,1\le j\le m\}$ and two vertices $(v_i,u_j)$ and $(v_k,u_l)$ are
adjacent in $G\times H$ if and only if ($v_i=v_k$ and $u_j\sim u_l$)
or ($v_i\sim v_k$ and $u_j=u_l$). The metric dimension of Cartesian
product graph is studied in \cite{pelayo1}.

 The following results related to the metric
dimension of the graph $(P_n\times P_m)\odot K_1$ and $(K_n\times P_m)\odot K_1$ were published in \cite{baskoro}.

\begin{theorem}{\em \cite{baskoro}} For $n\ge 1$ and $1\le
m\le 2$, $dim((P_n\times P_m)\odot K_1)=2$.
\end{theorem}

\begin{theorem}{\em \cite{baskoro}} For $n\ge 3$,
$$dim((K_n\times P_m)\odot K_1)=\left\{\begin{array}{ll}
                                        n-1, & m=1, \\
                                        n, & m=2.
                                      \end{array}
\right.$$
\end{theorem}

In this paper we show that the above results
are not correct for the case $m=2$ and $n\ge 3$. We also solve the
general case $m\ge 2$.

\section{Results}

\begin{theorem}
If $n\ge 3$ and $m\ge 2$, then $dim((P_{n}\times P_{m})\odot
K_1)=3$.
\end{theorem}

\begin{proof}
Let $\{v_1,v_2,...,v_n\}$ and $\{u_1,u_2,...,u_m\}$ be the set of vertices of the graphs $P_{n}$ and $P_{m}$,
respectively.  The vertices of $P_{n}\times P_{m}$ will be denoted by  $v_{ij}=(v_i,u_j)$ and the pendant vertex of
$v_{ij}$ in  $(P_{n}\times P_{m})\odot K_1$ will be denoted by $u_{ij}$.
We will show that $S=\{v_{11},v_{1m},v_{nm}\}$ is a
resolving set for $(P_{n}\times P_{m})\odot K_1$. The representations
of vertices of $(P_{n}\times P_{m})\odot K_1$ with respect
to $S$ are given by the following expressions,
\begin{align*}
r(v_{ij}|S)&=(d(v_{ij},v_{11}),d(v_{ij},v_{1m}),d(v_{ij},v_{nm}))\\
&=(i+j-2,m+i-j-1,m+n-i-j),
\end{align*}
\begin{align*}
r(u_{ij}|S)&=(d(u_{ij},v_{11}),d(u_{ij},v_{1m}),d(u_{ij},v_{nm}))\\
&=(i+j-1,m+i-j,m+n-i-j+1).
\end{align*}
Now, let us suppose there exist two different vertices $x,y$ of
$(P_{n}\times P_{m})\odot K_1$ such that $r(x|S)=r(y|S)$. If
$x=v_{ij}$ and $y=v_{kl}$, then $i\ne k$ or $j\ne l$ and we obtain
that
$$(i+j-2,m+i-j-1,m+n-i-j)=(k+l-2,m+k-l-1,m+n-k-l).$$
Which leads to $i=k$ and $j=l$, a contradiction. Analogously we
obtain a contradiction if $x=u_{ij}$ and $y=u_{kl}$. On the other
hand, if $x=v_{ij}$ and $y=u_{kl}$, then we have
$$(i+j-2,m+i-j-1,m+n-i-j)=(k+l-1,m+k-l,m+n-k-l+1),$$ which leads to
 $1=-1$, a contradiction. So,
for every different vertices $x,y$ of $(P_{n}\times P_{m})\odot
K_1$, we have $r(x|S)\ne r(y|S)$. Therefore, $dim((P_{n}\times
P_{m})\odot K_1)\le 3$.

On the other hand, since $(P_{n}\times
P_{m})\odot K_1$ is not a path, $dim((P_{n}\times
P_{m})\odot K_1)\ge 2$. Now let us suppose $S'=\{a,b\}$ is a resolving set for
$(P_n\times P_{m})\odot K_1$. If there exist two different paths of
length $d(a,b)$ between $a$ and $b$, then there exist two different
vertices $c,d$ of $(P_n\times P_{m})\odot K_1$ such that
$d(c,a)=d(d,a)$ and $d(c,b)=d(d,b)$, a contradiction. Let us suppose there is
only one path $Q$, of length $d(a,b)$, between $a$ and $b$. Thus,
all the vertices of $Q$, except possibly $a$ or $b$ which could be
pendant vertices, belong either to a copy of $P_n$ or to a copy of
$P_m$. We consider the following cases.

Case 1: If every vertex belonging to the path $Q$ has degree less or
equal than three, then $m=2$ and $S'\subset
\{u_{11},v_{11},u_{21},v_{21}\}$ or $S'\subset
\{u_{1n},v_{1n},u_{2n},v_{2n}\}$. Let us suppose $S'\subset
\{u_{11},v_{11},u_{21},v_{21}\}$. Now, for the vertices
$u_{1i},v_{1,i+1}$, $2\le i\le n-1$ we have that
\begin{align*}
d(u_{i1},a)&=d(u_{i1},v_{11})+d(v_{11},a)\\
&=d(v_{i+1,1},v_{11})+d(v_{11},a)\\
&=d(v_{i+1,1},a),
\end{align*}
\begin{align*}
d(u_{i1},b)&=d(u_{i1},v_{11})+d(v_{11},b)\\
&=d(v_{i+1,1},v_{11})+d(v_{11},b)\\
&=d(v_{i+1,1},b).
\end{align*}
Thus, $r(u_{i1}| S')=r(v_{i+1,1}|S')$, a contradiction.
On the contrary, if $S'\subset
\{u_{1n},v_{1n},u_{2n},v_{2n}\}$, then for the vertices
$u_{i1},v_{i-1,1}$, $2\le i\le n-1$ we have
\begin{align*}
d(u_{i1},a)&=d(u_{i1},v_{1n})+d(v_{1n},a)\\
&=d(v_{i-1,1},v_{1n})+d(v_{1n},a)\\
&=d(v_{i-1,1},a),
\end{align*}
\begin{align*}
d(u_{i1},b)&=d(u_{i1},v_{1n})+d(v_{1n},b)\\
&=d(v_{i-1,1},v_{1n})+d(v_{1n},b)\\
&=d(v_{i-1,1},b).
\end{align*}
Thus, $r(u_{i1}| S')=r(v_{i-1,1}|S')$, a contradiction.

Case 2: There exists a vertex $v$ of degree four belonging to the path $Q$. So, $v$ has two neighbors $c,d$ not
belonging to $Q$, such that $d(c,a)=1+d(v,a)=d(d,a)$ and
$d(c,b)=1+d(v,b)=d(d,b)$. Thus, $r(c| S')=r(d|S')$, a contradiction.
Hence,
$dim((P_{n}\times P_{m})\odot
K_1)\ge 3$. Therefore, the result follows.
\end{proof}

The
following lemmas are useful to obtain
the next result.

\begin{lemma}{\em \cite{pelayo1}}\label{remark-kn-pm}$\,$
If $n\ge 3$ then $dim(K_n\times P_m)=n-1$.
\end{lemma}

\begin{lemma}{\em \cite{buczkowski}}\label{one-pendant-edge}
If $G_1$ is a graph obtained by adding a pendant edge to a
nontrivial connected graph $G$, then
$$dim(G)\le dim(G_1)\le dim(G) + 1.$$
\end{lemma}

\begin{theorem}
If $m\ge 2$, then
$$dim((K_{n}\times P_{m})\odot K_1)=\left\{\begin{array}{cc}
    n-1, & \textrm{for $n\ge 4$,} \\
    3, & \textrm{for $n=3$.}
  \end{array}\right.
$$
\end{theorem}

\begin{proof}
Similarly to the above proof, let $v_{ij}=(v_i,u_j)$ be the set of
vertices of $K_{n}\times P_{m}$, where $v_i$, $1\le i\le n$ and
$u_j$, $1\le j\le m$ are vertices of the graphs $K_{n}$ and $P_{m}$,
respectively. Let us  denote by $u_{ij}$ the pendant vertex of
$v_{ij}$. Assume that $n=3$. We will
show that $S=\{v_{11},v_{21},v_{3m}\}$ is a resolving set for $(K_3\times P_m)\odot K_1$. Let
us consider two different vertices $x,y$ of $(K_3\times P_{m})\odot
K_1$. We have the following cases.

Case 1: $x=v_{ij}$ and $y=v_{kl}$. If $j=l$, then $i\ne k$ and
either $i\ne 3$ or $k\ne 3$, say $i\ne 3$. So, for $v_{i1}\in S$ we
have $d(x,v_{i1})=j-1<j=d(y,v_{i1})$. On the contrary, say $j<l$. If
$i\ne 3$ or $k\ne 3$, for instance, $i\ne 3$, then for $v_{i1}\in S$
we have $d(x,v_{i1})=j-1<l-1\le d(y,v_{i1})$. Now, if $i=k=3$, then
$d(x,v_{3m})=m-j>m-l=d(y,v_{3m})$.

Case 2: $x=u_{ij}$ and $y=u_{kl}$. Is analogous to the above case.

Case 3: $x=v_{ij}$ and $y=u_{kl}$. If $j=l$ and $i=k=3$, then we
have $d(x,v_{3m})=m-j<m-j+1=d(y,v_{3m})$. Also, if  $j=l$ and ($i\ne
3$ or $k\ne 3$), say $i\ne 3$, then for $v_{i1}\in S$ we have
$d(x,v_{i1})=j-1<j\le d(y,v_{i1})$. On the other hand, if $j\ne l$,
we consider the following subcases.

Subcase 3.1: $i=k$ and $i\ne 3$. If $j=l+1$, then we have that
$d(x,v_{3m})=m-j+1=m-l<m-l+2=d(y,v_{3m})$. On the other hand, if
$j\ne l+1$, then for $v_{i1}\in S$ we have $d(x,v_{i1})=j-1\ne
l=d(y,v_{i1})$.

Subcase 3.2: $i=k=3$. If $j=l-1$, then there exists $v_{r1}\in S$,
$r\ne 3$ such that $d(x,v_{r1})=j=l-1<l+1=d(y,v_{r1})$. On the other
hand, if $j\ne l-1$, then we have that $d(x,v_{3m})=m-j\ne
m-l+1=d(y,v_{3m})$.

Subcase 3.3: $i\ne k$. Hence, we have either $i\ne 3$ or $k\ne 3$,
for instance $i\ne 3$. If $d(x,v_{i1})=j-1=d(y,v_{i1})$, then there
exist $v_{r1}\in S-\{v_{i1}\}$, $r\ne 3$, such that
$d(x,v_{r1})=j>j-1\ge d(y,v_{r1})$.

Therefore, $dim((K_3\times P_{m})\odot K_1)\le 3$.

On the other hand, let $S'=\{a,b\}$ be a resolving set for  $(K_3\times
P_{m})\odot K_1$. If
there exist two different paths of length $d(a,b)$ between $a$ and
$b$, then there exist two different vertices $c,d$ of $(K_3\times
P_{m})\odot K_1$ such that $d(c,a)=d(d,a)$ and $d(c,b)=d(d,b)$.  Hence, $r(c| S')=r(d|S')$, a contradiction. Moreover,
if there is only one path $Q$, of length $d(a,b)$,
between $a$ and $b$, then
there exists a vertex $v$ of degree four belonging to the path $Q$. So, $v$ has two neighbors $c,d$ not
belonging to $Q$, such that $d(c,a)=1+d(v,a)=d(d,a)$ and
$d(c,b)=1+d(v,b)=d(d,b)$. Thus, $r(c| S')=r(d|S')$, a contradiction.
Thus, $dim((K_3\times P_{m})\odot K_1)\ge 3$. Therefore, for $n=3$, the
result follows.

Now, let $n\ge 4$.  We will show that $S=\{v_{1m},v_{31},v_{41},...,v_{n1}\}$ is a
resolving set for $(K_{n}\times P_{m})\odot K_1$. Let us consider two
different vertices $x,y$ of $(K_{n}\times P_{m})\odot K_1$. We have
the following cases.

Case 1: $x=v_{ij}$ and $y=v_{kl}$. If $j=l$, then $i\ne k$. Let us
suppose $i=1$ and $k=2$. Hence for $v_{1,m}\in S$ we have
$d(x,v_{1m})=m-j<m-j+1=d(y,v_{1m})$. Now, if $i\notin \{1,2\}$ or
$k\notin \{1,2\}$, then we have $v_{i1}\in S$ or $v_{k1}\in S$, say
$v_{i1}\in S$. Thus, we have $d(x,v_{i1})=j-1<j=l=d(y,v_{i1})$.

On the other hand, if $j\ne l$, say $j<l$, then there exists
$v_{t1}\in S$, $t\in \{3,...,n\}$, $t\ne k$, such that
\begin{align*}
d(x,v_{t1})&=d(x,v_{i1})+d(v_{i1},v_{t1})\\
&\le j-1+d(v_{k1},v_{t1})\\
&< l-1+d(v_{k1},v_{t1})\\
&=d(y,v_{k1})+d(v_{k1},v_{t1})\\
&=d(y,v_{t1}).
\end{align*}

Case 2: $x=u_{ij}$ and $y=u_{kl}$. Since $d(u_{ij},v)=d(v_{ij},v)+1$
for every $v\in S$, we proceed analogously to the above case and we
obtain that $r(u_{ij}|S)\ne r(u_{kl}|S)$.

Case 3: $x=v_{ij}$ and $y=u_{kl}$.   If $j\le l$, then for every $v_{t1}\in S$ we have
\begin{align*}
d(x,v_{t1})&=d(x,v_{i1})+d(v_{i1},v_{t1})\\
&=j-1+d(v_{i1},v_{t1})\\
&<l-1+d(v_{i1},v_{t1})\\
&\le l+d(v_{k1},v_{t1})\\
&=d(y,v_{k1})+d(v_{k1},v_{t1})\\
&=d(y,v_{t1}).
\end{align*}
Now, if $j>l$, then we have
\begin{align*}
d(x,v_{1m})&=d(x,v_{im})+d(v_{im},v_{1m})\\
&=m-j+d(v_{im},v_{1m})\\
&<m-l+d(v_{im},v_{1m})\\
&\le m-l+1+d(v_{km},v_{1m})\\
&=d(y,v_{km})+d(v_{km},v_{1m})\\
&=d(y,v_{1m}).
\end{align*}
Therefore, for every two different vertices $x,y$ of $(K_n\times
P_{m})\odot K_1$ we have, $r(x|S)\ne r(y|S)$ and, as a consequence, $S$ is a resolving set for $(K_n\times P_{m})\odot K_1$ of cardinality $n-1$.

On the other hand, by Lemma \ref{remark-kn-pm} and Lemma
\ref{one-pendant-edge} we have $dim((K_{n}\times P_{m})\odot K_1)\ge
n-1$. Hence, for $n\ge 4$, the result follows.
\end{proof}

\section*{Acknowledgements}
This work was partially supported   by the
Spanish Ministry of Education through projects TSI2007-65406-C03-01
 ``E-AEGIS" and Consolider Ingenio 2010 CSD2007-00004
``ARES''.


\section{INTRODUCTION}
To understand the event structure observed in different laboratories
i.e. $SLAC$, $HERA$, $DESY$ etc, and especially the one would be
expected in $LHC$ ($CERN$), the theoretical formalisms which
describe the small $x$ ($x$ is $Bejorken$ variable) region are
vital. The main unknown parameters  in these models are the
unintegrated parton distribution functions ($UPDF$)
\cite{3a,3b,3c,3d}. The $UPDF$ are two-scales dependent
distributions which are functions of $x$ (longitudinal momentum
fraction of the parent hadron) and the scales $k_{t}^2$ and
$\mu^{2}$, the squared transverse momentum of the parton and the
factorization scale, respectively. As we pointed out these
distributions are the  essential ingredients for the less exclusive
phenomenological computations in the high energy collisions of
particle physics.

It is well known that in the region of high energy and moderate
momentum transfer i.e. small $x$, the collinear factorization
theorem i.e.  $Dokshitzer$-$Gribov$-$Lipatov$-$Altarelli$-$Parisi$
($DGLAP$) \cite{5a,5b,5c,5d} evolution, breaks down. This happens
because of  the large increase of the phase space available for the
gluon emissions (i.e. a rapid rise in the gluon density), which
makes the quantum chromodynamics ($QCD$) perturbative expansions
unjustified and one can not obtain the $UPDF$. On the other hand, at
above high energy limit,  the cross section can be predicted by
using the $k_{t}$ factorization and the
$Balitsky$-$Fadin$-$Kuraev$-$Liptov$ ($BFKL$)
\cite{bfkl1,bfkl2,bfkl3} evolution. But the precision of $k_{t}$
factorization is not good e.g. the next-to-leading order ($NLO$)
corrections to $BFKL$ are very large \cite{bfkl4,bfkl5,bfkl6,bfkl7}.
Another approach to derive the $UPDF$ is the
$Ciafaloni$-$Catani$-$Fiorani$-$Marchesini$ ($CCFM$) equations
\cite{1a,1b,1c,1d,1e}. Although the $CCFM$ equations describe the
evolution of the $UPDF$ correctly, but working in this framework is
a complicated task, so practically they are used only in the Monte
Carlo event generators \cite{2b,2c,2d,2e,2f}. On the other hand, up
to now, there is not a complete quark version for these kind of
  equations \cite{1a,1b,1c,1d,1e,webber}, since   the
enhanced terms that are resumed by $CCFM$ come from gluon evolution.
However, to over come this problem, it has been shown that the
$CCFM$ equation can be reformulated (the linked dipole chain model)
by reducing the division between the initial and the final state
radiation diagrams using the colour dipole cascade model
\cite{LDCM1,LDCM2,Artem}.

The $Kimber$, $Martin$ and $Ryskin$ ($KMR$) \cite{4} approach is an
alternative prescription for producing the $UPDF$ which is based on
the standard $DGLAP$ equations \cite{5a,5b,5c,5d},
\begin{eqnarray}
\frac{\partial a(x,\mu^{2})}{\partial \ln(\mu^{2})}=\sum_{a'=q,g}
P_{aa'}\otimes a'(y,\mu^{2}),
 \label{1}
\end{eqnarray}
where $a(x,\mu^{2})=xq(x,\mu^{2})$ or $ xg(x,\mu^{2})$ and
$P_{aa'}(z)$ are the conventional (integrated) parton distribution
functions (PDF) and the well known $DGLAP$ splitting functions,
respectively. In equation (\ref{1}) the symbol $\otimes$ denotes a
convolution as,
\begin{eqnarray}
 f\otimes g=\int_x^{1}\frac{dy}{y}f(\frac{x}{y})g(y).
 \label{2}
\end{eqnarray}
In this approach under the certain approximation the $UPDF$ are
obtained from the PDF by introducing the scale $\mu$ only in the
last step of evolution with the inclusion of angular ordering
constraint ($AOC$). It has been shown that the $KMR$ prescription
gives the same results both for the $DGLAP$ and the unified
$BFKL$-$DGLAP$ equations \cite{4p} and the $AOC$ is applicable to
all orders as in the $CCFM$ formalism, i.e.
 all the loops contributions via the chain of evolution which are restricted by $AOC$, are resumed.

In this work, along  the lines of our recent calculations
\cite{8a,8b}, we intend to use the $KMR$ prescription with $MSTW
2008$ \cite{9} $PDF$ to produce three dimensional plots of $UPDF$ at
different energies ($\mu$) and discussed the various behavior of
$UPDF$ i.e. $f(x,k_t,\mu)$ . So the paper is organized as follows:
In section $II$ we briefly introduce the $KMR$ formalism and
finally, section $III$ is devoted to the results and the discussions
concerning the three dimensional (3D) graphs of the $UPDF$ produced
via this approach.
\section{The $KMR$ formalism \cite{4}}
The $KMR$ prescription \cite{4} works as a machine that by taking a
defined PDF as inputs, generates $UPDF$, as outputs. Using the
leading order (LO) splitting functions, $P_{aa'}$, the $DGLAP$
equations can be written in a modified form as \cite{4},
\begin{eqnarray}
\frac{\partial a( x,\mu^{2})}{\partial
\ln(\mu^{2})}=\frac{\alpha_s}{2\pi}\left[\int_x^{1-\Delta}P_{aa'}(z)\,a'\left(\frac{x}{z}
, \mu^{2}\right)dz-a(x ,
\mu^{2})\sum_{a'}\int_0^{1-\Delta}dz'P_{a'a}(z')\right],
 \label{3}
\end{eqnarray}
where $\Delta$ is a cutoff to prevent $z=1$ singularities in the
splitting functions arising from the soft gluon emission. In the
conventional $DGLAP$ formalism, $\Delta=0$ and the singularities are
canceled by the virtual terms. The value of $\Delta$ can be
determined by imposing an appropriate dynamical condition which is
replaced  by the angular ordering  constraint  arising  from the
coherency of the gluon emissions \cite{6a,6b},
\begin{eqnarray}
 ...>\theta_{n}>\theta_{n-1}>\theta_{n-2}> ...,
 \label{4}
\end{eqnarray}
where $\theta$'s are the radiation angels. This condition, at the
final step of  evolution, leads to \cite{1a,1b,1c,1d,4p},
\begin{eqnarray}
\mu>\frac{zk_t}{1-z}\Rightarrow
\Delta=1-z_{max}=\frac{k_t}{\mu+k_t}.
 \label{5}
\end{eqnarray}
The first part of the equation (\ref{3}), shows the contribution of
real emissions, that can change the transverse momentum $k_t$. The
second term expresses the evolutions due to the virtual effects
without changing the $k_t$. The latter can be re-summed, to obtain a
survival probability factor,
\begin{eqnarray}
T_a(k_t,\mu)=\exp\left[-\int_{k_t^2}^{\mu^2}\frac{\alpha_s({k'_t}^2)}{2\pi}
\frac{{dk'_t}^{2}}{{k'_t}^{2}}\sum_{a'}\int_0^{1-\Delta}dz'P_{a'a}(z')\right].
 \label{6}
\end{eqnarray}
Now, similar to the $Sudakov$ form factor, the above survival
probability, equation ({\ref{6}),  is imposed  into the equation
(\ref{1}), and by using equation (\ref{2}),  we find the equation
which describes the $UPDF$,
\begin{eqnarray}
f_{a}(x,k_{t}^2,\mu^{2})&=&T_a(k_t,\mu)\left[\left.
\frac{\partial\,a(x,\mu^{2})}{\partial\,ln(\mu^{2})}\right|_{\,\mu^2=k_t^2}\right]_{real}\nonumber\\
&=&T_a(k_t,\mu)\frac{\alpha_s({k_t}^2)}{2\pi}\int_x^{1-\Delta}P_{aa'}(z)\,a'\left(\frac{x}{z}
, {k_t}^2 \right)dz.
 \label{7}
\end{eqnarray}
More explicit forms of the above equation for the gluon $g$ and the
different quark flavors $q=u, d, s, ...$ are as follows,
\begin{eqnarray}
f_{q}(x,k_{t}^2,\mu^{2})&=&T_q(k_t,\mu)\frac{\alpha_s({k_t}^2)}{2\pi}
\nonumber\\&\times&
\int_x^{1-\Delta}dz\left[P_{qq}(z)\frac{x}{z}\,q\left(\frac{x}{z}
, {k_t}^2 \right) + P_{qg}(z)\frac{x}{z}\,g\left(\frac{x}{z} ,
{k_t}^2 \right)\right],
 \label{8}
\end{eqnarray}
and
\begin{eqnarray}
f_{g}(x,k_{t}^2,\mu^{2})&=&T_g(k_t,\mu)\frac{\alpha_s({k_t}^2)}{2\pi}
\nonumber\\&\times& \int_x^{1-\Delta}dz\left[\sum_q
P_{gq}(z)\frac{x}{z}\,q\left(\frac{x}{z} , {k_t}^2 \right) +
P_{gg}(z)\frac{x}{z}\,g\left(\frac{x}{z} , {k_t}^2 \right)\right].
 \label{9}
\end{eqnarray}
The key observation here is the dependency on the scale $\mu^2$,
which appears at the last step of the evolution. Another point is
that, the $Sudakov$ form factor which arises from the resumption of
virtual effects, can be used at every order of approximation.
Although the splitting functions must be used at the $NLO$ level,
but as it is shown in \cite{37},  the $NLO$ corrections to the
splitting functions, are relatively small in comparison to the $LO$
contributions. However, as stated above, only the $LO$ splitting
functions are used. On the other hand, although the definition of
$Sudakov$ form factor (like the $PDF$ themselves) has been started
intuitively from a probabilistic interpretation, but its role in the
mathematical description of the evolution remains in the equations.

The primary computations based on this kind of approach to evaluate
the $UPDF$, show very good agreement with the experimental data for
$F_2$ \cite{4}. Also, in recent years, the $KMR$ prescription have
been widely used for phenomenological calculations (see \cite{8a}
and the references therein). Recently the stability and the
reliability of the $KMR$ \ $UPDF$ have been investigated in
\cite{8a,8b}.

Finally, we should mention here that, the key property of the $CCFM$
approach (as given in their publications
\cite{1a,1b,1c,1d,1e,2b,2c}) is the $AOC$, which in turn has root in
the coherency of gluon radiation along the evolution chain, that is
valid for whole range of $x$ values. In the conventional $DGLAP$
formalism, the strong ordering constraint on the transverse momenta,
restricts the domain of study to the large and moderate values of
$x$:
$$
\hat{ \sigma}(\gamma^\ast q\rightarrow
qg)=\int_{p_{t_{min}}}^{p_{t_{max}}}dp^2_t{d\hat{\sigma}\over
{dp^2_t}},
$$
where
$$
p^2_{t_{min}}=\lambda^2,
$$
and
$$
p^2_{t_{max}}=p_t^2|_{\sin^2 \theta=1}=k^{\prime^2}={\hat{ s}\over
4}=Q^2{1-x\over 4x}.
$$
So to obtain the $DGLAP$ equations with $\ln({\hat{ s}\over
4})\simeq \ln (Q^2)$, $x$ should not be very low. In the $KMR$
prescription the $AOC$ property of the $CCFM$ formalism is applied
to modified $DGLAP$ evolution as a cut off on the integrals.
Therefore, the results of these modifications show that the effect
of application of $AOC$ is even more important than the inclusion of
the conventional low $x$ effects in the $BFKL$ approach \cite{4}.

\section{Results and discussion}
 As we stated in the section $II$, by using the equations (\ref{8}) and (\ref{9}), the $UPDF$ are
generated via the $KMR$ procedure. For the input $PDF$, the $MSTW
2008$ \cite{9} set of partons at the $NLO$ level are used
\footnote[1]{We use the $MSTW 2008$ code that is accessible from
\emph{http://durpdg.dur.ac.uk} and \emph{http://www.hepforge.org}.}.
Since the generated $UPDF$ ($f_{a}(x,k_{t}^2,\mu^{2})$) are three
variable functions, by fixing the scale $\mu^2$, their values versus
$x$ and $k_{t}^2$ are plotted in the various panels of figures 1, 2,
3 and 4 for the gluons, the up, the strange and the bottom  quarks,
respectively. For the better comparison, the values of the $\mu^2$
are chosen in a wide range $\mu^2=10, 10^{2}, 10^{4}, 10^{8}$
$GeV^2$ which is up to the $LHC$ working scales. The three typical
quark flavors, the $u$ quarks consists of the valence and the sea
contributions $u=u_v+u_{sea}$ and the $s$ and the $b$ quarks which
are completely sea distributions, are presented. The main feature of
these figures is exhibiting the general behavior of the $UPDF$ with
respect to the coupled contributions of $x$ and $k_t^2$. For
example, the most probable value of $k_t^2$ ($x$) at every $x$
($k_t^2$) for any kind of partons can be checked. As it can be seen,
by increasing the scale $\mu^2$ the graphs are shifted to the higher
$k_t^2$. This is expected, since the probability of finding partons
with larger $k_t^2$ is more probable at higher scales. The growth of
the values of the distributions by increasing $\mu^2$ and decreasing
$x$ and also the phenomenon of converging the quark distributions to
a unique value at small $x$ are known characteristics of the parton
distributions which are the heritage of their parent $PDF$. The
different behaviors of up and strange quarks at large $x$ have root
in the valence contribution in the case of up quark. The pronounced
peaks   become wider with respect to $k_t^2$, and move to higher
values of $k_t^2$. This behavior is much effective  for the up, the
strange and the bottom quarks. The peaks come from the concept of
distributions  and they are results of the dynamical evolution of
partons. The figures show that at given values of  hard scale and
$x$, at which $k_t$, it is more probable to detect the out going
partons. So based on the final partons, we can predict the dynamical
properties of the produced jets and their components, and on the
other hand it can inform us about the precision of the current
theoretical formalisms itself. The input PDF of $MSTW 2008$ are also
given in the figure 5, for comparison. With good approximation by
integrating over $UPDF$, we can get the input   $MSTW 2008$ PDF
($a(x,\mu^2)=\int^{\mu^2} {{dk_t^2\over k_t^2}}f_a(x,k_t^2,\mu^2)$).
For example for gluons, at $x=0.01$ and $\mu^2=100$ we get $6.7$
whereas $MSTW 2008$ gives the value of $6.5$ i.e. 3\% off.
Situations are the same for other points and parton distributions.
It is worth to say that in the original $KMR$ work, they get 25\%
discrepancies \cite{4} for above comparison. This is also evident by
comparison of figure 5 with those of 1 to 4 i.e. the $UPDF$ are
decreasing by increasing $x$. On the other hand, as have been
discussed in the $KMR$ and other related works, because of the
imposition of angular ordering, the $UPDF$ have values for
$k_t^2\geq \mu^2$ as $x$ decreases. But this will not affect the
above integration too much. The figures 1 to 4 also show that, for
low scales ($\mu^2\simeq 10$ $GeV^2$) the $UPDF$ become negative
when $x$ becomes close to one. This reflects the negative values of
$MSTW 2008$ gluon distributions at the $NLO$ level and beyond that.
So the negative values of $UPDF$ have root in the parent integrated
gluon distributions which in turn are the result of $MSTW 2008$
assumptions \cite{9}.  As it was pointed, in the $MSTW2008$
\cite{9}, for better data fitting it is allowed that, the gluon
distribution takes negative values, because there is no theorem that
imposes positivity condition on $PDF$ beyond the  $LO$
approximation. So they become negative  in order to fit the data (in
other words they can be traced to the slow evolution of $F_2$ at
small $x$ and  $Q^2$ i.e. a positive gluon would give too rapid
evolution to fit the $dF_2/d\ln(Q^2)$ data. Then in the $KMR$
integrals, the evaluation of input $g(x,k_t^2)$ at small $x$ and
$k_t$ (as a scale, instead of $Q^2$  in $g(x,Q^2)$) leads to the
negative values for the output $UPDF$. Finally, $(i)$ the comparison
of $UPDF$ produced from different  $PDF$ sets have been made in our
former works \cite{8a,8b} . The different parameterizations
procedures lead to different $PDF$, and a discussion about these
procedures is presented in \cite{8a,8b} and references therein.
$(ii)$ The differences between  the $LO$ and the $NLO$ $PDF$ are
parameterizations dependent. In the $MSTW2008$ this is noticeable,
but in some other parameterizations sets based on different
assumptions and procedures it can be less (e.g $GRV$ sets
\cite{8a,8b}), but as we have showed in \cite{8a,8b} (by
investigating the ratios of $KMR$ $UPDF$ compared to the
corresponding ratios of input $PDF$) the relative differences are
less in the output $UPDF$ and the $KMR$ prescription suppresses
these discrepancies. To show this point more transparently,  in
figure 6 we have plotted the gluon $UPDF$ with three different input
$PDF$, namely the original $KMR$ \cite{4} with $MRST99$ \cite{99}
$PDF$, our recent works \cite{8a,8b} with $GJR08$ $PDF$ \cite{GRV}
and present calculation ($MRST2008$) at $\mu^2=100$ $GeV^2$ and
$x=0.1, 0.01, 0.001$ and $0.0001$ in terms of $k_t^2$. It is clearly
seen that different input $PDF$ give very similar $UPDF$. $(iii)$ In
fact a complete prescription for producing the $NLO$ $UPDF$ needs to
include both the $PDF$ and the splitting functions at the $NLO$
level. This prescription is presented in \cite{37}, but as it is
shown in this reference \cite{37}, inclusion of the $NLO$ splitting
functions have very low effect comparing to the contribution of the
$NLO$ $PDF$. Therefore, ignoring the corrections due to the $NLO$
splitting functions do not affect our analysis of the general
behavior of the $NLO$ $UPDF$. $(iv)$ There is no restriction on the
$k_t$ dependency. As the orders of the approximation are in terms of
orders of $\alpha_s(k_t^2)$, the $NLO$ accuracy is contained in the
NLO PDF and splitting functions that discussed in the former
comments. Hence at scales $k_t^2\geq Q_0^2$ ,where $Q_0^2$ is the
scale that upper than it, the perturbative $QCD$ is still
applicable, these results are valid.

\begin{acknowledgments}
We would like to acknowledge  the Research Council of University of
Tehran and Institute for Research and Planning in Higher Education
for the grants provided for us. $HH$ would like to thank Professors
$Hannes$ $Jung$, $Francesco$ $Hautmann$, $Artem$ $Lipatov$, $Gavin$
$Salam$, $James$ $Stirling$, $Graeme$ $Watt$ and $Bryan$ $Webber$
for their helpful communications.
\end{acknowledgments}