# Analysing spatial point patterns in digital pathology: immune cells in high-
grade serous ovarian carcinomas

Jonatan A. González
Computer, Electrical and Mathematical Sciences and Engineering (CEMSE),
King Abdullah University of Science and Technology, Jeddah, Saudi Arabia
and
Julia Wrobel
Department of Biostatistics and Informatics,
Colorado School of Public Health, University of Colorado,
Anschutz Medical Campus, Aurora, Colorado, USA
and
Simon Vandekar
Department of Biostatistics, Vanderbilt University Medical Center,
Nashvill, Tennessee, USA and
Paula Moraga
Computer, Electrical and Mathematical Sciences and Engineering (CEMSE),
King Abdullah University of Science and Technology, Jeddah, Saudi Arabia

###### Abstract

Multiplex immunofluorescence (mIF) imaging technology facilitates the study of
the tumour microenvironment in cancer patients. Due to the capabilities of
this emerging bioimaging technique, it is possible to statistically analyse,
for example, the co-varying location and functions of multiple different types
of immune cells. Complex spatial relationships between different immune cells
have been shown to correlate with patient outcomes and may reveal new pathways
for targeted immunotherapy treatments.

This tutorial reviews methods and procedures relating to spatial point
patterns for complex data analysis. We consider tissue cells as a realisation
of a spatial point process for each patient. We focus on proper functional
descriptors for each observation and techniques that allow us to obtain
information about inter-patient variation.

Ovarian cancer is the deadliest gynaecological malignancy and can resist
chemotherapy treatment effective in cancers. We use a dataset of high-grade
serous ovarian cancer samples from 51 patients. We examine the immune cell
composition (T cells, B cells, macrophages) within tumours and additional
information such as cell classification (tumour or stroma) and other patient
clinical characteristics. Our analyses, supported by reproducible software,
apply to other digital pathology datasets.

Keywords: Multiplex single-cell imaging; Multitype point patterns; Second-
order descriptors; Spatial point patterns; Tumour microenvironment

##  1 Introduction

Novel spatial omics assays, including multiplex immunofluorescence imaging,
imaging mass cytometry, and spatial transcriptomics, offer the ability to
provide simultaneous high-resolution spatial and proteomic or transcriptomic
cell-level information. These spatial assays are expected to provide valuable
spatial information about interactions between cells and cellular
neighbourhoods that inform scientific understanding of how the immune system
responds within the tumour microenvironment (TME) across tumour types and
stages of cancer development 6; 29; 21; 53. Tumour microenvironments are
complex and dynamic ecosystems with various cell types, including cancer and
immune cells. The presence and function of immune cells in the tumour
microenvironment have been shown to play a critical role in tumour growth,
progression, and response to therapy. Image processing techniques combined
with spatial statistical methods can answer complex questions about the
interactions of cell subtypes within the TME 54. For example, recent analyses
of multiplex imaging data have been applied to ovarian cancer and indicate
that spatial interactions of specific types of immune cells are prognostic of
patient survival 50.

Novel analytic methods have been developed to address analysis goals in
characterising spatial relationships between cells. For example, Kim et al.32
describe a novel method for analysing single-cell transcriptomics data using a
combination of deep learning and graph theory. This approach allows for
identifying cell types and exploring their gene expression patterns, which can
be especially useful in spatial analysis applications. Similarly, Zhenghao et
al. 7 introduce Spatial-LDA, a topic modelling algorithm that considers the
spatial organisation of cells in a tissue sample. This method can infer gene
expression patterns and identify cell types with a spatial organisation,
providing insight into the spatial heterogeneity of the sample. Zhenzhen et
al. 58 present a new algorithm for solving the $K$-nearest neighbour problem
using locality-sensitive hashing. This algorithm has potential applications in
spatial analysis tasks such as computer vision and recommendation systems.
Together, these papers highlight the importance of spatial analysis methods in
modern biology research.

Many spatial omics methods face theoretical and practical limitations that may
affect their application due to implicit simplifying assumptions that are
often ignored. The theoretical limitations 2 include the (1) assumption of
complete spatial randomness (CSR), one of the most common simplifying
assumptions, whereby the cell locations occur within the sampled tissue in a
completely random fashion; (2) homogeneity or constant intensity, which means
that the expected value of the number of cells in every subset of the sampled
tissue remains constant; (3) the assumption of no interaction among cells,
i.e., no attraction or repulsion; and (4) the independence between immunity
markers (or any other characteristic associated with each cell) and cell
locations. Most spatial omics methods do not explicitly define their
assumptions, limiting their transparency and reproducibility. In addition,
many methods do not permit inference across slides (patients), restricting
their utility. These practical limitations could impede the results’ accuracy,
reliability, and replicability.

Spatial point processes are convenient and widely used mathematical concepts
for analysing event distributions that can be considered as points in a
geometric space 39. A variety of methods and models allow us to understand how
points are located in some space and how they interact 28; 19; 2. These
methods study the possible independence between points and consider from
complete randomness to aggregation or repulsion conceived from various point
generation mechanisms 23. Furthermore, marked point process theory, a
particular topic of spatial point processes, can be used to study differences
in spatial characteristics between samples with more realistic data models.
While these methods were not developed for the analysis of features of the
tumour microenvironment, they can be used to provide rigorous statistical
inference. These methods can be applied in spatial omics analysis to create
accurate models of intercellular processes and their relationship with various
pathologies 51 to understand the spatial organisation of various cells in
tissues within the tumour microenvironment (TME) 56; 55.

In this tutorial, we use spatial point process methodologies to analyse the
spatial distribution and interactions of immune cells in ovarian cancer tissue
samples from several patients collected using multiplex immunofluorescence
imaging. We focus specifically on B-cells, macrophages, CD8 T-cells, and CD4
T-cells, immune cell subtypes identified as essential players in the tumour
microenvironment of ovarian cancer. Through this investigation, we hope to
gain a deeper understanding of the spatial organisation of these immune cell
subtypes and their association with patient-level outcomes. Taken as a whole,
our tutorial explicitly outlines and tests underlying point process
assumptions using a large open-source ovarian cancer dataset. We provide a
principled and reproducible set of analyses for single-cell spatial omics data
intended to guide future work in this area. To aid in ease of use, we provide
code as supplementary material.

##  2 Data

###  2.1 Spatial omics data structure

Multiplex imaging might have more than 40 channels of markers that bind to
specific proteins and can be used to establish cell identification and
function. Many datasets, including the one we employ in this work, have fewer
markers. After image acquisition, which varies substantially based on the
imaging platform, image processing steps are applied before data can be
analysed as a point pattern, including autofluorescence adjustment, single-
cell segmentation, and cell phenotyping. Segmentation identifies individual
cells, nuclei, and tissue areas (i.e. tumour and stroma regions) so that
analyses can be performed at the single-cell level. Then cells are phenotyped
or given a cell-type label based on the expression of different immune and
other markers using clustering or gating marker channels 24; 49; 59. Finally,
cell marker intensities must be normalised to remove non-biological batch
effects across slides or other experimental units 20; 56; 25.

###  2.2 Ovarian tumour microenvironment data

Data acquisition and processing details can be found in Jordan et al.30, and
Steinhart et al. 50. The data are available through a public repository on
Bioconductor, 57. Briefly, 128 cancerous tissue samples from 51 human ovaries
were obtained through 5-micron slices of the tissue microarray (TME) stained
with specific antibodies for CD8 (T cells, C8/144B, Agilent Technologies),
CD68 (macrophages, KP1, Agilent Technologies), cytokeratin (CK, tumour cells,
AE1/AE3, Agilent Technologies), CD3 (T cells, LN10, Leica), and CD19 (B cells,
BT51E, Leica). Immune cells were classified according to the following
criteria: B-cell (CD19+), macrophage (CD68+), CD8 T-cell (CD3+ and CD8+), CD4
T-cell (CD3+ and CD8-). Throughout this article, we consider those samples
with at least eight cells with each of the four categories of immunity markers
in order to have a sufficiently balanced sample to perform the analyses.
Figure 1 shows a composite image from our dataset (a) and point patterns of
immune cells for four randomly selected subjects (b).

(a)

(b)

Figure 1: (a): Image of an ovarian tumour sample stained with antibodies
specific for CD4+ (T-cells), CD8+(T-cells), CD68+ (macrophages), and
cytokeratin (tumour cells) from the VectraPolarisData ovarian cancer dataset
on Bioconductor. (b): Locations and immune classification of ovarian tissue
cells of four out of fifty-one randomly selected patients with high-grade
serous ovarian cancer.

Each immune cell can belong to the tumour or stroma compartment. A tumour is
an abnormal growth of cells that can occur in any part of the body. Stroma, on
the other hand, refers to the supportive tissue surrounding the cells of an
organ or tissue. Additionally, we have several clinical covariates for each
patient, including:

  * •

Whether the tumour is primary from the initial diagnosis. This is a factor
with two categories: Yes or no (binary)

  * •

Whether the tumour has undergone chemotherapy before image acquisition
(binary)

  * •

Whether the patient has a BRCA (breast cancer gene) mutation (binary)

  * •

Whether the patient received a PARPi inhibitor (binary)

  * •

Cancer stage, a categorical factor which can take the values Stage I, Stage
II, Stage III, Stage IV

  * •

Patient age (continuous)

  * •

Patient survival status (time-to-event)

###  2.3 Scientific objectives

The cells in each patient image are considered a multitype marked point
pattern. Roughly speaking, a multitype marked point pattern is a collection of
points in a spatial domain, where each point is associated with multiple
attributes or “marks”, and the points are categorised into multiple types
based on some additional characteristics or attributes. At a high level, our
goal is to characterise each patient point pattern and extract features that
can be used to analyse patient-level outcomes (e.g., survival). The specific
scientific aims we address with this analysis are detailed below.

Our first scientific aim is to explore whether the number of immune cells of
different subtypes relates to patient outcomes. This analysis does not
incorporate spatial or point process information but is commonly the first
analysis step for cell imaging data when several tissue samples (several
patients) are available.

Our next scientific aim is to use intensity functions of the point processes,
which are first-order point pattern summary statistics described in detail in
Section 3.2.1, to explore and visualise the density and distribution of cells
and continuous-valued functional antibody markers in single images. One
quantity we obtain from this analysis is a spatial relative risk of different
cell types (e.g., the relative risk of B-cells compared to macrophages in a
given tissue sample), and within an image, hypothesis tests can be performed
to test the significance of relative risk at a particular spatial location.
Intensity functions of continuous markers, for example, pSTAT3, show how these
markers are distributed throughout the tissue.

The third scientific objective is to use second-order summary statistics of
point processes, detailed in Section 3.3, to explain relationships between
cell types within the tumour microenvironment. Multitype second-order
descriptors such as Ripley’s $K$-function explore pairwise relationships
between cell types, and we can conduct hypothesis tests to determine if
spatial clustering is significant within an image.

Finally, we discuss how second-order summary statistics of spatial
interactions between cells can be associated with patient-level outcomes
(i.e., compared across images).

##  3 Methods

###  3.1 Point process fundamentals

We introduce the following notation: $\mathbf{X}$ will denote a point process
and $X$ a point pattern, i.e., a realisation of $\mathbf{X}$ 18; 2. The points
in a point pattern refer to the spatial locations of a set of individual
events or objects that are of interest in a particular study or analysis. In
this specific context, points represent cell location within a TME. Therefore,
we indistinctly discuss points, spatial or cell locations throughout this
paper.

$N(B)$ will denote the number of points of $\mathbf{X}$ in a planar region
$B\subset W\subset\mathbb{R}^{2}$. Since spatial locations are recorded with
$m$ different labels, one for each cell type, we are interested in multitype
point processes analysis. We then divide the process $\mathbf{X}$ into subsets
denoted by $\mathbf{X}^{(1)},\ldots\mathbf{X}^{(M)}$, where each
$\mathbf{X}^{(m)}$, $m=1,\ldots,M$, is the point process with points of type
$m$. In this case, we will denote by $N_{m}(B)$ the points of type $m$ in the
region $B$. When the point process is considered regardless of its labels or
marks, it is denoted by $\mathbf{X}^{(\bullet)}$. Note that the terms “marks”
and “markers” arise from different contexts. The former is related to the
characteristics associated with each cell location, for example, cell type or
cell size, and it comes from point process analysis, whereas the latter comes
from the biological context of antibody markers, for example, B-Cell or
macrophages. In order to avoid very complex notation, the subscript of the
descriptors will be understood according to the context. The subscripts $j$
and $m$ usually represent marks (B-Cell, macrophages, etc.), and the
subscripts $i$ and $k$ usually represent point patterns (tissue samples or
patients).

In practice, the observations window $W$, which describes the area where cells
can occur, is usually attached to the data points. However, when the window is
not given, we need to be able to propose a suitable window for the data. This
can be done in several ways, for example, by putting all the points inside the
minimum rectangle that contains the data or inside the convex hull of the
sample. However, these methods have the disadvantage that many data points lie
on the edge of the window, which can introduce unwanted edge effects. A
frequently used method that avoids this difficulty is Ripley and Rasson’s 47,
which defines the window as a dilation of the convex hull (centred at the
centroid of the convex hull that contains the data) by a factor of
$1/\sqrt{1-\omega/n}$, where $n$ is the number of data points, and $\omega$ is
the number of vertices of the convex hull.

###  3.2 Analysing first-order descriptors

In univariate statistics, moments are fundamental quantities that provide a
valuable mean of describing probability distributions of random variables. The
spatial distributions of points have analogous quantities. For example, the
numerical means and variances of random variables are replaced by moment
measures 12; 8. Particularly, first-order measures are those associated with
the expected values of a variable, which in the case of spatial point
patterns, this random variable represents the number of points in each subset
of the observation region.

####  3.2.1 Intensity function

The intensity function is the first-order descriptor of a point process, i.e.,
the intensity function describes the expected value of points of some kind
anywhere in the observation window; it is given by

$\lambda_{m}(\mathbf{u})=\lim_{|\mathrm{d}\mathbf{u}|\rightarrow{0}}\frac{\mathbb{E}\left(N_{m}(\mathrm{d}\mathbf{u})\right)}{|\mathrm{d}\mathbf{u}|},\qquad\mathbf{u}\in
W,m=1,\ldots,M.$ (1)

When $\lambda$ assumes a constant value, the process is called stationary or
homogeneous. This assumption is violated in the ovarian cancer data set, as
seen by the uneven spatial dispersion of cells in each sample (Figure 1). The
unmarked point process $\mathbf{X}^{(\bullet)}$, which does not distinguish
between cell types, has an intensity given by

$\lambda_{\bullet}(\mathbf{u})=\sum_{i=1}^{M}\lambda_{i}(\mathbf{u}).$ (2)

Many techniques are available to estimate first-order intensity functions,
including adaptive methods 17; 5; 2; 14. Of all these techniques, non-
parametric estimation by kernels is the most common, as well as most
appropriate for spatial omics data. Kernel estimation depends on a bandwidth
$(\epsilon)$ that controls the amount of smoothing and is given by the
standard deviation of the kernel function. Many methods can be used to
estimate the bandwidth: cross-validation 17; 3, likelihood cross-validation
36, Scott’s rule of thumb 48, Cronie and van Lieshout’s criterion 10 among
others 2.

Adaptive bandwidth selection is preferred over fixed bandwidth in intensity
estimation because it allows a more accurate estimation of the intensity
function in different data regions. Depending on the underlying data
distribution, fixed bandwidth selection can lead to oversmoothing or
undersmoothing of the intensity estimate. On the other hand, adaptive
bandwidth selection adjusts the bandwidth locally based on the density of the
data points near the location of interest.

In practice, many of these bandwidth estimation methods do not always agree,
which often causes the intensity not to be appropriately estimated, especially
when some regions of the point pattern are particularly crowded, and some
other areas seem to be empty or at least very sparsely populated by points 14.
This is true in the samples in Figure 1(b). In these cases, assigning a large
bandwidth will produce oversmoothing in regions with many points and
undersmoothing in areas with few. Adaptive kernel estimators were developed to
address this limitation of fixed bandwidth estimators. We employ an adaptive
kernel estimator to accommodate the varying density of cells observed in
spatial omics data (Figure 1), defined as follows,

$\hat{\lambda}_{m}\left(\mathbf{u}\right)=\frac{1}{e\left(\mathbf{u}\right)}\sum_{i=1}^{n_{m}}{\kappa_{\epsilon(\mathbf{u}_{i})}\left(\mathbf{u}-\mathbf{u}_{i}\right)},\qquad\mathbf{u}\in
W,\mathbf{u}_{i}\in X^{(m)},m=1,\ldots M,$

where $n_{m}$ is the number of points of $X^{(m)},\kappa()$ is the Gaussian
kernel and $e()$ represents an edge correction factor given by

$e\left(\mathbf{u}\right)=\int_{W}\kappa_{\epsilon(\mathbf{u^{\prime}})}\left(\mathbf{u}-\mathbf{u^{\prime}}\right)d\mathbf{u^{\prime}}.$
(3)

The bandwidth function $\epsilon()$ is defined as

$\epsilon(\mathbf{u})=\frac{\epsilon^{\star}}{\gamma}\sqrt{\frac{n_{m}}{\lambda^{\star}_{m}(\mathbf{u})}},$

where $\epsilon^{\star}$ is a fixed global bandwidth that can be estimated
through classical methods such as maximal smoothing principle 52;
$\lambda^{\star}_{j}(\mathbf{u})$ is a pilot estimate of the intensity of
$X^{(m)}$ usually computed through classical kernel estimation with fixed
bandwidth $(\epsilon^{\star})$, and $\gamma$ is the geometric mean term for
the pilot intensity evaluated in the points of the point pattern $X^{(m)}$.
The geometric mean is intended to free the bandwidth factor from dependence on
the data scale 15.

We calculate the global bandwidth using Scott’s isotropic rule 48. The pilot
estimate is the usually fixed bandwidth $(\epsilon^{\star})$ kernel estimate.

####  3.2.2 Segregation

In the context of multitype point processes, segregation refers to the
tendency of points of different types to spatially separate from one another.
A point process is a mathematical model that describes the spatial
distribution of points in a region of interest, where each point can be
assigned a type or category. Segregation can occur when points of one type
tend to cluster together and are spatially separated from points of other
types.

In a multitype point process with a spatially varying probability distribution
of types, the intensity function $\lambda_{m}(\mathbf{u})$ can be used to
calculate the conditional probability,

$p(m|\mathbf{u})=\frac{\lambda_{m}(\mathbf{u})}{\lambda_{\bullet}(\mathbf{u})},$

of a point being of type $m$ given its location $\mathbf{u}$. When all the
types share a common baseline intensity, i.e., when
$\lambda_{m}(\mathbf{u})=a_{m}\beta(\mathbf{u})$, $p(m|\mathbf{u})$ becomes
constant, meaning that spatial variation in the probability distribution of
types is equivalent to segregation of types. To obtain non-parametric
estimates of $p(m|\mathbf{u})$, kernel smoothing estimates of the numerator
and the denominator can be substituted into the formula.

We can use a non-parametric Monte Carlo test for the spatial segregation
between different cell types in a multitype point process 16. This test is
based on a randomised version of the statistic given by

$T=\sum_{k=1}^{n}\sum_{m=1}^{M}\left[\hat{p}_{m}(m|\mathbf{u}_{k})-n_{m}/n\right]^{2}$

In this case, $\hat{p}_{m}(m|\mathbf{u}_{k})$ is the leave-one-out kernel
estimate of the probability that the $k-$th point has the label $m$.

If the Monte Carlo test indicates that the observed spatial pattern of cells
is significantly different from what would be expected under complete spatial
randomness, then this suggests spatial segregation between the different cell
types. This could imply that the cells are interacting with each other or with
the surrounding environment in a non-random manner, leading to the observed
spatial organisation.

####  3.2.3 Spatial relative risk

The concept of spatial relative risk dates back to Bithell (1990)5 and began
as a way of viewing and interpreting the relative abundance of cases of a
disease with respect to the distribution of the population at risk in a given
geographic region. Relative risk is a measure of the strength of association
between two densities. It can be used to analyse spatial cell data to compare
two cell type densities or, for example, compare a cell type density with the
complete density, including all the other types.

A ratio between two intensity functions is known as relative risk31. In
general, the spatial log relative risk, defined as a ratio of two intensity
functions, is considered in the literature instead of the relative risk and is
estimated by

$\rho_{ij}(\mathbf{u})=\log\left\\{\frac{\hat{\lambda}_{i}(\mathbf{u})}{\hat{\lambda}_{j}(\mathbf{u})}\right\\}+\log\left\\{\frac{n_{j}}{n_{i}}\right\\},\quad
i,j=1,\ldots,M,$

Fixed or variable bandwidth kernels usually perform the numerator and
denominator estimators (Section 3.2.1). When the bandwidth is fixed,
calculating the numerator and denominator bandwidths individually does not
usually give good results 15; it tends to under-smooth the relative risk
surface. This means that, in general, one would have to choose a considerably
larger bandwidth than those used for individual estimates. In the adaptive
case, selecting a joint global bandwidth is convenient as doing the pilot
estimates separately, i.e., we set a shared adaptive bandwidth in numerator
and denominator. Terrel’s maximal principle 52 is often used for estimating
the global bandwidth 15.

By examining the relative risks across different regions within the ovaries,
we can identify areas with higher or lower chances of a given cell type and
potentially identify patterns or clusters of high-risk areas. However, it is
essential to remember that the relative risk alone cannot establish causality
and that other factors, such as confounding or bias, may also influence the
results.

####  3.2.4 First-order analysis of continuous marks

Some marker channels for spatial omics data can be continuously expressed
within the tissue. We call “marks” these characteristics. For example, among
the many marks of the TME of a given tissue sample, our dataset has three: ck
marker, ki67 marker, and pSTAT3. Specifically, pSTAT3 is a functional marker
used to detect a specific type of cell signalling.

We can consider a spatially weighted average of the marks at each point in the
observation window for these markers. This measure can be calculated using a
Nadaraya Watson estimator 44,

$\tilde{m}(\mathbf{u})=\frac{\sum_{i}m_{i}\kappa_{\epsilon}(\mathbf{u}-\mathbf{u}_{i})e^{-1}(\mathbf{u}_{i})}{\sum_{i}\kappa_{\epsilon}(\mathbf{u}-\mathbf{u}_{i})e^{-1}(\mathbf{u}_{i})},\qquad\forall\mathbf{u}\in
W,$

where $m_{i}$ are the mark values, $\kappa_{\epsilon}()$ is a kernel function
with bandwidth $\epsilon$, and $e(\mathbf{u}_{i})$ is Eq. (3) evaluated in
each data point, this factor is known as Diggle’s edge correction 17; 2.

Through this technique, we can understand the distribution of specific
continuous variables within the observation region. Variables sampled across
cell locations can be interpolated to the entire observation region, in this
case, the tissue sample of a given patient. Note that this technique allows us
to interpolate a continuous mark in all window locations, even if no cells
exist. Therefore, we predict in a specific location where it might be
impossible for a cell to locate itself.

###  3.3 Analysing second-order descriptors

Just as the first-order intensity is the statistic analogous to the mean in a
univariate distribution, second-order statistics are analogous to the variance
or covariance of a marked point process. In this case, this measure can be
associated with the count of pairs of points and is an index of the
statistical association between points 2. This can rigorously quantify
relationships between cell types or markers, such as the spatial association
between B-cells and Macrophages in ovarian cancer samples. Several statistics
allow for analysing of these second-order statistics. We employ some non-
parametric descriptors below.

####  3.3.1 Multivariate second-order characteristics

In the context of spatial point processes, it is common to use second-order
summary descriptors, such as the pair correlation function or the
$K$-function, to describe the spatial arrangement of points in a point
pattern. Homogeneous versions of these descriptors assume that the intensity
of points is constant throughout each sample, which may not hold true in the
case of spatial omics data. Therefore, it is necessary to use inhomogeneous
versions of these descriptors that consider the variations in the point
process intensity across different types of points and locations. This
approach allows for a more accurate description of the spatial patterns and
relationships within the point pattern data.

One of the most popular second-order descriptors for analysing point patterns
is Ripley’s $K$-function 45. Ripley’s $K$-function in its inhomogeneous
multivariate version can be written by

${K}_{ij}(r)=\mathbb{E}\left[\left.\sum_{\mathbf{u}_{k}\in
X^{(j)}}\frac{1}{\lambda(\mathbf{u}_{k})}\mathbf{1}\\{||\mathbf{u}-\mathbf{u}_{k}||\leq
r\\}\right|\mathbf{u}\in\mathbf{X}^{(i)}\right],$ (4)

as long as this value does not depend on the choice of the location
$\mathbf{u}$ 39; 2 and where $r$ belongs to a suitably chosen range
$T=(0,r_{0}]$, 26. An estimator for the $K$-function is given by

$\hat{K}_{ij}(r)=\frac{1}{|W|}\sum_{x_{\ell}\in X^{(i)}}\sum_{x_{k}\in
X^{(j)}}\frac{\mathbf{1}\\{||\mathbf{u}_{k}-\mathbf{u}_{\ell}||\leq
r\\}e(\mathbf{u}_{k},\mathbf{u}_{\ell};r)}{\hat{\lambda}_{i}(\mathbf{u}_{\ell})\hat{\lambda}_{j}(\mathbf{u}_{k})},$
(5)

where $|\cdot|$ denotes, in this particular case, the area of the window $W$,
and $e(\mathbf{u},\mathbf{v};r)$ is an edge-correction weight (see, e.g. 46;
2). The $L$-function is a popular transformation of the $K$-function and it is
given by4

$L(r):=\sqrt{\frac{K(r)}{\pi}}.$

This transformation initially intends to convert the theoretical Poisson
(Complete Spatial Randomness) $K$-function given by $K_{CSR}(r)=\pi r^{2}$, to
the straight line $L_{CSR}(r)=r$ making a visual inspection of the plot more
straightforward. In the multitype context, this benchmark represents the
independence or no association between types $i$ and $j$ without the
assumption of CSR.

When there are many types of cells, comparing a particular type with the
remaining cell types is convenient35. For this, we can define a second-order
descriptor in the sense of the one presented in Eq. (2). The dot $K$-function,
$K_{i\bullet}(r)$, is the expected number of cells of any type lying within a
distance $r$ of a point of type $i$, standardised by dividing by the intensity
of the unmarked point pattern, i.e.,

${K}_{i\bullet}(r)=\mathbb{E}\left[\left.\sum_{\mathbf{u}_{k}\in
X^{(\bullet)}}\frac{1}{\lambda_{\bullet}(\mathbf{u}_{k})}\mathbf{1}\\{||\mathbf{u}-\mathbf{u}_{k}||\leq
r\\}\right|\mathbf{u}\in\mathbf{X}^{(i)}\right],\quad r\geq 0.$ (6)

The corresponding $L_{i\bullet}(r)$ is defined similarly.

Now, consider the cross-type pair correlation function defined as the
derivative of the $K$ in the following way28; 8,

$g_{ij}(r):=\frac{1}{2\pi r}\frac{\text{d}}{\text{d}r}K_{ij}(r).$

This function considers the contributions of points lying at distances equal
to $r$, and its benchmark under CSR (in the univariate case) is 1. In the
multivariate case, the benchmark is also one, but it represents independence
when the indexes $i$ and $j$ are different. Specifically, $g_{ii}(r)$ is
consistent with CSR, and $g_{ij}(r),i\neq j,$ is consistent with a lack of
correlation between types $i$ and $j$2. A pragmatic and useful interpretation
is conceived through probabilities: the probability $p(r)$ of finding two
points of types $i$ and $j$ at locations $\mathbf{u}$ and $\mathbf{v}$
separated by a distance $r$ is

$\lambda_{i}(\mathbf{u})\lambda_{j}(\mathbf{v})g_{ij}(r)\mathrm{d}\mathbf{u}\mathrm{d}\mathbf{v},\quad\mathbf{u},\mathbf{v}\in
W,$

where $\mathrm{d}\mathbf{u}$ and $\mathbf{v}$ stand for the area of two
infinitesimal regions around $\mathbf{u}$ and $\mathbf{v}$, respectively.

Multivariate second-order descriptors are statistical tools that can analyse
spatial patterns in multiplex imaging data. Specifically, they can be used to
investigate the spatial relationships between several different types of
molecules or cellular structures in the same tumour microenvironment by
measuring the extent to which two types of cells are spatially clustered or
dispersed relative to each other. The functions generate plots showing the
expected distance between cell pairs of each type as a function of the
distance between them, allowing researchers to visualise the degree of
clustering or dispersion.

####  3.3.2 Shortest distance

As Baddeley et al. 2 pointed out, there is sometimes a confusing ambivalence
between counting points in a point pattern in an observation region and
measuring the shortest distances. For instance, we can say that there are four
cells (on average) in a given $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ tissue
sample or that one cell appears every $0.25\text{\,}\mathrm{\SIUnitSymbolMicro
m}$. This duality impacts the mathematics behind point pattern analysis;
therefore, a complete statistical analysis should include the study of
shortest distances in addition to correlation.

Let $\mathbf{X}$ be a spatial stationary multitype point process and
$\mathbf{u}\in\mathbb{R}^{2}$ a fixed location. Then if we consider the
distance
$d(\mathbf{u},\mathbf{X}^{(m)}):=\min\left\\{||\mathbf{u}-\mathbf{u}_{k}||:\mathbf{u}_{k}\in\mathbf{X}^{(m)}\right\\}$
(called the empty-space distance), the empty-space function is the cumulative
distribution function of the empty-space distance, i.e.,

$F_{m}(r):=\mathbb{P}\left\\{d(\mathbf{u},\mathbf{X}^{(m)})\leq r\right\\},$

for all the interpoint distances $r$. The nearest neighbour distance
distribution function $G_{ij}(r)$ is the cumulative distribution function of
the nearest neighbour distance from a point of type $i$ to the nearest point
of type $j$, which can be written by as
$d(\mathbf{u},\mathbf{X}^{(j)}\setminus\mathbf{u})$, i.e., the shortest
distance from a point $\mathbf{u}$ to the process
$\mathbf{X}^{(j)}\setminus\mathbf{u}$. Thus,

$G_{ij}(r):=\mathbb{P}\left\\{d(\mathbf{u},\mathbf{X}^{(j)}\setminus\mathbf{u})\leq
r|\text{ }\mathbf{u}\text{ is a typical point of }\mathbf{X}^{(i)}\right\\},$

for any location $\mathbf{u}$. The theoretical values of the univariate
versions of these functions for CSR processes are
$F_{\text{CSR}}(r)=G_{\text{CSR}}(r)=1-\exp\\{-\lambda\pi r^{2}\\}$, where
$\lambda$ is the intensity of the point process; it means that for completely
random patterns, the empty space and the nearest-neighbour distance have the
same distribution 18; 39; 2. In the multitype case, $G_{ij}$ measures the
association between types $i$ and $j$, when there is independence,
$G_{ij}(r)=F_{j}(r)$.

In practice, it is convenient to compare these two functions ($F$ and $G$)
because although in the case of independence (or CSR in the univariate case),
they coincide, in the other cases, they tend to behave oppositely; that is,
while one increases, the other decreases. We can then consider a new function
as the quotient between these two, simultaneously summarising both behaviours.
This quotients is known as the $J$-function, defined as 34; 11

$J_{ij}(r):=\frac{1-G_{ij}(r)}{1-F_{j}(r)},\quad\forall r\geq 0,\text{ such
that }F_{j}(r)<1.$

Under independence, $J_{ij}(r)=1,\forall r$, values of $J$ above 1 represent
repulsion and values below 1 represent attraction.

In the multivariate non-stationary case, the functions $F_{j}$, $G_{ij}$ and
$J_{ij}$ require advanced mathematical treatment to be defined and
estimated11. However, its interpretation remains as simple as in the
stationary bivariate case. This is why we recommend its use as a complement to
other descriptors such as $K$ and $G$. In the same vein, it is possible to
define dot versions, $G_{i\bullet}$ and $J_{i\bullet}$ by considering
distances to points of any other type in a multitype point process.

####  3.3.3 Complete Spatial Randomness and Independence

When the data points are of a single type, we ignore cell type; the null
reference model is complete spatial randomness 19; 2. Although this model is
unrealistic for most real-world phenomena, it functions as a dividing
hypothesis9 or benchmark from which spatial aggregation or inhibition
(regularity) is inferred.

For point patterns of various types, there is an analogous reference model to
which complete spatial independence is added to complete spatial randomness.
For multitype point patterns, there are two possible choices for a benchmark:
random labelling, where type labels are randomly assigned to points, and
independence, where points of different types are independent. When both
conditions hold together, the point process is known as having Complete
Spatial Randomness and Independence (CSRI).

In the homogeneous case, the theoretical values of all second-order
descriptors such as $K_{ij}$, $L_{ij}$, $g_{ij}$, $K_{ii}$, $L_{ii}$,
$g_{ii}$, $K_{i\bullet}$, $L_{i\bullet}$, $g_{i\bullet}$, $G_{ij}$, $J_{ij}$,
$G_{ii}$, $J_{ii}$ and $F_{i}$ are consistent with CSRI. Therefore, performing
statistical tests to detect discrepancies is easy. However, this does not
necessarily occur in the inhomogeneous case. For this case, the $J$-function
can be used to test independence and random labelling.

The traditional approach to test independence between two types of points $i$
and $j$ consists of selecting one of the point patterns, e.g., $X^{(i)}$, and
making random displacements of its points as many times as desired37. Each
shifting breaks the possible spatial dependency between types without
affecting the dependency within the type. After each displacement, the desired
second-order descriptor is calculated and stored in a vector of observations
on which a global envelope test can be performed later (see Section 3.3.4).
Each random shifting will produce some points that lie outside the reference
window. There are some methods for dealing with these points 40; in our case,
we use the erosion technique40. This method uses an eroded window $W_{c}$ such
that $W_{c}\setminus\mathbf{u}\subset W$ for all possible shift vectors
$\mathbf{u}$. Then the statistics are computed with the restricted point
patterns to the new window, i.e., $X^{(i)}|_{W_{c}}$. Since some points are
lost when eroding, this method often loses power in the tests. Shifting a
point pattern $X$ could affect its distribution in the inhomogeneous case
since the intensity changes with the observation region. The intensity of the
pattern should be shifted with each random displacement to fix this problem.

In the case of random labelling, the null hypothesis is that each cell
immunity marker (the label) is determined randomly, independently of other
cells, with constant probabilities. A Monte Carlo test may be deliberated in
the following way: To generate different observations of the wanted second-
order statistic, it is only necessary to randomly permute the labels through
the cells of the tissue samples without changing their locations. If the null
hypothesis is true, e.g., if the immunity markers are random, then the point
patterns that result from the relabelling are statistically equivalent to the
original point pattern. The “dot” functions are the most useful for evaluating
random labelling for envelope tests2; 11.

####  3.3.4 Global envelope test

The statistical treatment of point patterns is usually complex, especially
when testing hypotheses about the interaction (covariance) between points.
Summary statistics are often used to depict the characteristics of the
observations in a mathematical fashion; for example, various statistics
capture spatial interaction. These functions depend on interpoint distances
$r$; we denote this type of function as $T(r)$. We described classic
descriptors here: Ripley’s $K$ function, the $L$ function, the nearest
neighbour distance distribution function $G(r)$, the empty space function
$F(r)$ and the $J$ function 39; 28; 8; 18; 2.

Choosing a proper second-order functional descriptor involves assessing the
research question and the characteristics of the point pattern data. For
example, the cross-type $K$-function helps to identify the clustering or
regularity of points of different types relative to each other, while the
$J$-function is useful in detecting the association between different marks.
It is essential to select the appropriate second-order functional descriptor
for the research question and data properties, considering the types of marks,
spatial scales of interest, and spatial dependence 2.

Once we have chosen an appropriate statistic that summarises a point pattern,
testing against a null model is helpful to determine whether the cell types or
markers exhibit interesting spatial patterns. We want to know if our estimated
statistic from the data is incompatible with its distribution under some null
model. In practice, the null distribution is often unknown, so using Monte
Carlo techniques to generate the null distribution is typical in spatial
statistics. Once the necessary simulations are done, the envelope tests are
performed.

In the context of multiplex imaging, functional envelope tests can be applied
to investigate the distribution of any functional descriptor or biomarker
calculated from the TME data. For example, researchers may be interested in
studying the correlation between the expression of two proteins in the same
tissue sample and how this correlation varies across different tissue regions.
Envelope tests can provide a way to quantify the degree of similarity or
difference between functions computed from the samples and functions computed
from synthetic data.

Normally, we can simulate point patterns from some null model. More precisely,
we can generate simulations under $H_{0}:T(r)=T_{\text{obs}}(r)$. Then we can
apply the following steps to test the null hypothesis.

  1. 1.

Compute the test statistic $T_{0}(r):=T_{\text{obs}}(r)$, for example Ripley’s
$K$-function.

  2. 2.

Generate $s$ simulated patterns $X_{i},i=1,\ldots,s$ from the null model, for
example, Poisson (no interaction).

  3. 3.

For each $i=1,\ldots,s,$ compute the test statistic $T_{i}(r)$ based on the
point pattern $X_{i}$.

We can use a non-parametric global envelope test based on a measure called the
extreme rank length (ERL) 43; 42. The ERL measure can compare the functions
$T_{i}$ without multiple testing problems 41. The test may also be interpreted
graphically as it points out the distances where the data contradicts the null
hypothesis.

The functions $T_{i}$ should be evaluated on a fixed number of distances
$r_{1},\ldots,r_{d}\in(0,r_{0}]$, so every function in finitely discretised by
$T_{i}=T_{i}(r_{1}),\ldots,T_{i}(r_{d})$. The ERL measure ranks the $T_{i}$
among each other. If $R_{0j},R_{1j},\ldots,R_{sj}$ are the ranks of
$T_{0}(r_{j}),T_{1}(r_{j}),\ldots,T_{s}(r_{j})$, such that the largest
$T_{i}(r_{j})$ has rank 1. The vector of pointwise ranks
$(R_{i1},R_{i2},\ldots,R_{id})$ is associated with each $T_{i}$. We define
$\mathbf{R}_{i}:=(R_{i[1]},R_{i[2]},\ldots,R_{i[d]})$ as these pointwise ranks
ordered from smallest to largest, i.e. $R_{i[j]}\leq R_{i[j^{\prime}]}$ if and
only if $j\leq j^{\prime}$. The ERL measure is then defined by using the
lexicographic ordering $\prec$ of the $\mathbf{R}_{i}$ 43, and the ERL measure
is given by

$E_{i}=\frac{1}{s+1}\sum_{j=1}^{s}\mathbf{1}\left\\{\mathbf{R}_{j}\prec\mathbf{R}_{i}\right\\}.$

Then, the $p$-value associated with the Monte Carlo test is

$p=\frac{1}{s+1}\sum_{i=0}^{s}\mathbf{1}\left\\{E_{i}\leq E_{1}\right\\}.$

In the case of composite null hypotheses, where some parameters of the null
model must be estimated, the tends to be conservative. There are some double-
stage Monte Carlo testing approaches to tackle this issue 13; 1.

##  4 Results

###  4.1 Associations between cell counts and phenotypic variables

Before getting into any spatial analyses, it is important to understand how
the counts of the immune cells are related to continuous and categorical
regressors. To analyse possible differences in cell populations between
patients, we use a linear log-Poisson model for the expected immune cell
counts. We ungroup each patient count into several cell types corresponding to
the number of B-cells, CD4 T-cells, CD8 T-cells and Macrophages within the two
tissue domains of the ovary samples: stroma or tumour. In this particular
context, let $n_{ijk}$ denote the number of immune cells of the patient
$i,(i=1,\ldots,g)$, for cell type $j,(j=1\ldots,m)$ of type of tissue
$k,(k=1,\ldots,l)$. In our case, $g=51,m=4$ and $l=2$. A Poisson log-linear
model for the expected cell counts may be used to analyse possible differences
in counts, considering the possible overdispersion of the data 38.

However, we opt in this context for a Generalised Estimation Equation (GEE)
model 33. These statistical models are used to analyse longitudinal data,
which involve multiple measurements of the same individuals or patients. In
longitudinal data, observations within each patient might be correlated, and
traditional statistical models that assume independent observations are
inappropriate. GEE models extend the generalised linear model (GLM) framework
to account for such correlation and allow for consistent estimation of model
parameters, even when the covariance structure is unknown or misspecified.

We consider the number of points of each tissue sample over the sampling areas
as an offset; it is included in the model to account for differences in the
tissue sample size and the number of points. We also model the possible
effects of design covariates $\mathbf{X}$ on cell counts. The model can be
expressed as

$\mathbb{E}\left(n_{ijk}|\mathbf{X}\right)=|W_{i}|\exp{\left(\mathbf{X}^{\top}\cdot\beta\right)},\quad\text{and}\quad\mathbb{V}\text{ar}\left(n_{ijk}|\mathbf{X}\right)=\phi\mathbb{E}\left(n_{ijk}|\mathbf{X}\right),$
(7)

where $\beta$ is the regression coefficient vector. $\phi$ accounts for
potential overdispersion. Table 1 shows the summary of the fitted model.

Regressor |  | Estimate | Naïve S.E. | Robust S.E |  $p$-value |
---|---|---|---|---|---|---
Intercept $(\beta_{0})$ |  | $0.03$ | $0.83$ | $1.58$ | $0.9856$ |
Immune |  |  |  |  |  |
| CD4 T-cells | $-0.07$ | $0.33$ | $0.24$ | $0.7785$ |
| CD8 T-cells | $1.56$ | $0.22$ | $0.16$ |  $<$$2\text{\times}{10}^{-23}$ | ${***}$
| Macrophages | $1.43$ | $0.24$ | $0.18$ |  $<$$2\text{\times}{10}^{-15}$ | ${***}$
Tissue |  |  |  |  |  |
| Tumour | $1.38$ | $0.16$ | $0.13$ |  $<$$2\text{\times}{10}^{-25}$ | ${***}$
Primary tumour |  |  |  |  |  |
| Yes | $-0.06$ | $0.19$ | $0.29$ | $0.8391$ |
Prior chemo |  |  |  |  |  |
| Yes | $-0.58$ | $0.20$ | $0.33$ | $0.0769$ | ${*}$
Cancer stage |  |  |  |  |  |
| Stage II | $-3.18$ | $1.10$ | $2.57$ | $0.2189$ |
| Stage III | $-0.68$ | $0.47$ | $0.45$ | $0.1304$ |
| Stage IV | $0.07$ | $0.49$ | $0.47$ | $0.8853$ |
BRCA mutation |  |  |  |  |  |
| Yes | $0.49$ | $0.19$ | $0.36$ | $0.1741$ |
PARPi inhibitor |  |  |  |  |  |
| Yes | $0.16$ | $0.17$ | $0.41$ | $0.6990$ |
Death |  |  |  |  |  |
| Yes | $-0.18$ | $0.16$ | $0.29$ | $0.5396$ |
Age at diagnosis |  | $-0.01$ | $0.01$ | $0.02$ | $0.8003$ |
|  |  |  ${}^{***}p<0.001$; ${}^{**}p<0.01$; ${}^{*}p<0.05$

Table 1: Estimates with associated naïve and robust standard errors for the
regression coefficients of generalised estimation equation (gee) model.

The estimated scale parameter is $\hat{\phi}=441.7\gg 1$, indicating an
extreme extra-Poisson variation of the cell counts within each combination of
the factors considered. This high additional variance could be fundamentally
due to a non-Poisson variation (departures from complete spatial randomness
(CSR)) within each sample, for example, due to the attraction or repulsion
between cells or to inter-patient variation in mean cell count, possibly
related to covariates that could not be measured in the study or random
effects.

In Table 1, we see significant increments in the counts of CD8 T-cells, and
Macrophages relative to the reference level, B-cells. We also appreciate that
there is more abundance of cells in tumour tissue than in the stroma. A common
side effect of chemotherapy is a decrease in immune cell counts, evidenced by
chemotherapy patients, who experienced a mildly significant decrease in immune
cell counts. The other regressors are not statistically significant, meaning
there is insufficient evidence for an influence on the counts.

###  4.2 Analysing expected counts of immune cells within TMEs

The first-order descriptors refer to the rate at which immune cells are
distributed in the tissue. This function describes the density of immune cells
at any given point in the tissue and can be used to understand the spatial
arrangements of these cells.

The observation window $W_{i}$, where $i$ denotes the patient’s number, is
defined through Ripley and Rasson’s technique 47. In this case, we have
different spatial windows for different point patterns, i.e., we have
different sampling regions for different patients, as every patient has an
associated point pattern. In a point process, the intensity function
represents the first-order properties, which describes how the mean number of
points in a dataset varies with spatial coordinates. This measurement is
fundamental in point processes because it shows how the mass or the total
number of points is distributed throughout the study region.

In Figure 2, we show, as an example, a first-order analysis of one of the
patients in the sample. We first show the point pattern associated with the
patient with the delimitation of the window $(W_{6})$ and without
differentiating between cell types, that is, without assigning any mark
neither discrete nor continuous (Figure 2 (a)). As each cell can have several
characteristics, we mainly highlight two of them, the type of immunity and the
type of tissue; when we choose one or the other, we will have two different
coloured point patterns (Figure 2 (b) and (c)).

####  4.2.1 Intensity function

Specifically, the first-order intensity function provides information about
the average number of immune cells per unit area and can help researchers
identify regions of the tissue where immune cells are more or less
concentrated. Understanding the first-order intensity function of a point
pattern in this context can provide insight into the immune response to
ovarian cancer and potentially inform strategies for improving immunotherapy
treatments.

Intensity estimates for each immune marker for a particular patient (patient
six) are displayed in Figure 2 (d). Note that the estimates can be performed
taking into account another set of labels such as tissue region (Figure 2
(c)).

Figure 2: Example of kernel estimates of the intensities of patient six. This
conveys information on the abundance of the immune cells and the estimated
probability of each immunity marker. (a) Representation of the immune cells
without considering marks. (b) Representation considering immune markers. (c)
Representation considering tissue type. (d) Kernel estimates of intensity for
each immune marker. (e) Spatial probability of each immune marker with respect
to the other types. (f) Most likely immunity markers at each location of the
tissue sample of patient six.

####  4.2.2 Segregation

We can speak about the probability distribution of different cell types. The
probability that any cell belongs to type $m$ is

$p_{m}(m|\mathbf{u})=\frac{\lambda_{m}(\mathbf{u})}{\lambda_{\bullet}(\mathbf{u})},$

as long as $\lambda_{\bullet}(\mathbf{u})\neq 0,$ where
$\lambda_{\bullet}(\cdot)$ is given in Eq. (2). Estimates of the varying
spatial proportion of each immune marker are displayed in Figure 2 (e),
suggesting that immune markers are segregated; for example, Macrophages appear
clustered more in the top of the sample, whereas B-cells are shown in higher
proportions in the right. Spatial segregation occurs when certain types of
points, in particular subregions, predominate rather than randomly mix within
some observation window. A texture plot can be used to summarise which immune
marker of four has the highest probability at this location, i.e., we divide
the tissue sample into regions where different immunity markers predominate
(Figure 2 (f)). We can appreciate how macrophages and B-cells predominate
throughout patient six’s sample.

We can apply this test using, for example, 999 simulations for each region of
interest in our dataset. It is typical to adjust the $p$-values to account for
multiple tests across the $g$ patients using Bonferroni correction 27. We
obtain mildly significant $p$-values (about 0.051) for the segregation of
immune markers.

####  4.2.3 Relative risk

Understanding the relative risk between different pairs of immune markers in
cancerous tissues such as ovaries may be relevant for the development of
effective diagnostic and treatment strategies. While each immune cell type
plays a distinct role in the immune response against cancer cells, their
functions might be intricately connected and can influence one another. By
understanding the relative risk associated with each pair of immune markers in
cancerous tissues, researchers can identify the most effective combinations of
immune markers for targeting cancer cells and develop therapies that enhance
the effectiveness of the immune response while minimising the potential
negative impact of immune cells that may promote tumour growth. Estimates of
spatial varying relative risks of immunity markers are displayed in Figure 3.

Figure 3: Adaptive bandwidth spatial log-relative risk surfaces of patient six
data for each combination of levels of Immune markers, with asymptotic
tolerance contours.

We can appreciate the northeast region of the tissue sample, particularly with
a significantly increased risk (up to 1.5 times) of having B-cells given the
counts of CD4 and CD8 T-cells; this risk is also slightly with respect to
Macrophages. The southwestern region, on the other hand, shows a significant
decrease in the first two cases.

####  4.2.4 Smoothing continuous markers

A marker for phosphorylated signal transducer and activator of
transcription-3, pSTAT3, is available for all cells in each sample. It is
present in multiple types of cells, including tumour cells, immune cells, and
stromal cells. pSTAT3 regulates various cellular processes, such as cell
growth, proliferation, and survival, and has been implicated in the
development and progression of numerous types of cancer. Due to its widespread
presence and functional significance, pSTAT3 has become a popular target for
cancer therapy research. pSTAT3 is known to be upregulated in ovarian cancer
tissue compared to normal ovarian tissue 22. pSTAT3 exhibits significant
spatial heterogeneity that differs across patients (Figure 4).

Figure 4: Spatially varying average pSTAT3 of four tissue samples (patients).

###  4.3 Analysing interactions amongst cells: correlation and spacing

At the microscopic scale, neighbour cell-cell interactions in cancer tissue
may involve a complex interplay between different cell types. Cancer cells can
interact with neighbouring cells through direct contact or by releasing
signalling molecules that bind to specific receptors on the surface of
neighbouring cells. These interactions can promote cancer cell proliferation,
migration, and invasion and modulate the behaviour of surrounding cells, such
as immune cells and stromal cells. We approach the complex nature of neighbour
cell-cell interactions at the microscopic scale through second-order
statistics that depend on a suitable range of distances $r\in T=(0,r_{0}]$
(see Section 3.3.1). In our case, we set a sensible default that depends on
the geometry of the ovaries; indeed,
$r_{0}=$350\text{\,}\mathrm{\SIUnitSymbolMicro m}$$.

####  4.3.1 Multivariate second-order descriptors

We use multivariate Ripley’s $K$ and pair correlation to quantify associations
between cell types across all samples in the dataset. The negative
discrepancies between the estimates and $L(r)=0$ (the benchmark) suggest that
the cell type pairs show a spatially repulsive behaviour. We can arrange the
centred $L$-functions by cell type pairs to have a simple judgement about the
interaction (Figure 5).

Figure 5: The centred cross-type $L$-functions $L_{ij}(r)-r$ for each pair of
types $i$ and $j$ of immune cells of all patients. Every colour and panel
represents a combination $ij$ of immune types of cells.

In this case, where the point pattern represents immune cells of two different
immunity markers, these values falling below the independence benchmark on
average might suggest that the two types of immune cells are more closely
associated with one another than would be expected by chance, i.e., they would
interact to each other by inhibiting. This might indicate a situation where
the immune system is overwhelmed by the cancer cells, and the immune cells
cannot work together efficiently. The inhibitory effect could also be due to
the competition between the two types of immune cells for resources or space
within the tumour microenvironment. However, this independence hypothesis
should be formally tested, and we present that analysis in the next section.

####  4.3.2 Complete Spatial Randomness and Independence

To test the independency, we use the bivariate inhomogeneous $J$-function and
the shifting technique described in Section 3.3.3. For illustration purposes,
we show the test results of a random patient in Figure 6; the same procedure
can be applied to each one in the sample.

Figure 6: $J$-functions $J_{ij}(r)$ for patient 3 where $i$ and $j$ are all
the cell types. The envelopes of 2999 simulations were generated by random
shifts (grey shading with extreme rank length). Results suggest the different
components are independent.

We obtain non-significant $p$-values for every combination of immunity markers
indicating that we are in the presence of independence. In the context of
immune cells with different immunity markers, the independence of the two
types would suggest that the two immunity markers do not have any significant
spatial association with each other. This could indicate that the two markers
behave independently within the tumour microenvironment. However, it is worth
noting that independence does not necessarily mean that the two types of
immune cells do not interact. There could still be functional interactions
occurring between the two types of immune cells that are not reflected in
their spatial distribution.

In the case of random labelling, we use the centred version of the dot
$J$-function and a set of 2999 permutations of the labels of the immunity
markers. The results for patient 3 are displayed in Figure 7.

Figure 7: Test of random labelling for patient 3. The summary function is
$\hat{J}_{i\bullet}(r)-\hat{J}(r)$, where each level $i$ is each panel title.
The solid lines are the row means of 2999 random labellings, and the grey
shading represents a central region based on the extreme rank length. Results
suggest that the cell types are not random labels.

We obtain enough statistical evidence to reject the random labelling
hypothesis in this case. If the immune cells in a cancerous tissue sample have
non-randomly assigned immunity markers, it might suggest that there is a
selective pressure on the immune cells to express specific markers that enable
them to recognise and target cancer cells more effectively, for instance. This
selective pressure could be due to the cancer cells, which may have developed
mechanisms to evade the immune system and require a more specific immune
response. Alternatively, it could be due to external factors such as
treatment, which may have induced a selection pressure on the immune cells to
express markers more effectively in combating the cancer cells.

###  4.4 Testing differences between patients groups

Finally, we can also test the hypothesis that two (or more) observed point
patterns are realisations of the same spatial point process model to
investigate whether all types of immune cells interact similarly across
patient groups. For example, we would want to compare whether the patients
with four different cancer stages have the same spatial distribution of immune
cells. We use the inhomogeneous $L$-function in this case, but we could use
any other second-order descriptor. The null hypothesis that we would like to
test can be formally formulated as,

$H_{0}:\bar{L}_{\text{stage I}}(r)=\bar{L}_{\text{stage
II}}(r)=\bar{L}_{\text{stage III}}(r)=\bar{L}_{\text{stage
IV}}(r),\quad\forall r\in T.$ (8)

We can establish analogous hypotheses for the patients grouped by other
factors such as prior chemotherapy or death.

Before using a specific statistic for testing $H_{0}$, we should evaluate the
assumption of homoscedasticity, i.e., that the variances are equivalent across
the different patient groups. We can evaluate the equality of variances of the
$L$-functions calculated for two or more groups through a Levene’s-style test
42. This equality would be convenient since some statistical procedures assume
that the populations’ variances are equal. We obtain $p$-values of
$0.121,0.520$ and $0.067$ for cancer stage, prior chemo and death,
respectively. It means that the $L$-functions show equal variances across
groups. The graphical output for the test for cancer stage as grouping factor
is shown in Figure 8(a).

To test the hypothesis $H_{0}$ shown in Eq. (8), we can use the original
$L$-functions or a suitable rescaled statistic in case of heteroscedasticity
42. We can apply a rank envelope test by permuting the $L$-functions in a non-
parametric one-way ANOVA fashion, i.e., permuting the functional descriptors
across the different patient groups. We set 50,000 random permutations and
obtain $p$-values of $0.829,0.652$ and $0.225$ for cancer stage, prior chemo
and death, respectively. We show the graphical output for the test in the case
of cancer stage in Figure 8(b) for illustration.

(a)

(b)

Figure 8: (a) The test for equality of variances of the centred $L$-functions
in the four cancer stages groups. (b) The graphical functional ANOVA based on
functional descriptors. Test for equality of means of the centred
$L$-functions in the four groups according to cancer status using the group
means.

This analysis suggests that the interaction of different immune cells remains
consistent across various cancer stages. The presence or absence of
chemotherapy did not significantly affect the interactions between these
immune cells. Additionally, the analysis found that patients who died during
the sampling did not display any statistical differences in cell interaction
compared to those who survived. These findings suggest that the cancer stage
does not strongly influence the interactions between immune cells, the use of
chemotherapy, or the patient’s survival.

##  5 Conclusion

This paper introduces non-stationary spatial point processes to model
interactions between cell types in ovarian cancer samples. We have made a tour
of the methods of point processes related to the descriptive analysis of
observations. These observations or realisations can be analysed separately or
together through replicated point process methodologies. We have seen how both
first-order characteristics, related to the spatial distribution of the number
of points, and second-order characteristics, related to the covariance
structure between points, can be estimated. These methods can reveal
biologically interesting features of ovarian cancer samples that can inform
our understanding of the immune response in ovarian cancer.

The techniques for analysing replicated patterns help study cell arrangements
in multiplex imaging because more than one sample is available, and inference
across samples is of primary interest in understanding how spatial features
relate to patient-level phenotypes.

Although the number of statistical procedures related to point patterns is
vast28; 19; 2, the concepts and contexts do not usually fit the reality of
cell patterns, especially when the tissue that has been observed is not
complete, that is when we observe only a piece of the tissue. This is the case
of the ovarian cancer samples analysed here, which present irregular geometric
shapes that change from one patient to another. In technical terms, these cell
characteristics translate into inhomogeneity and heteroskedasticity in the
functional descriptors8 since the functional form of the variances usually
implies window geometry and intensity. The methods we proposed for multiplex
image analysis here are appropriate for inhomogeneous and heteroskedastic
point processes, so they are ideal for modelling spatial interactions of cells
in these datasets. However, this problem of non-stationarity, far from being a
limitation, constitutes a marvellous open field of investigation.

Spatial statistics research in digital pathology has the potential to enhance
our understanding of disease progression and the effectiveness of treatments.
However, analysing large-scale digitised tissue images requires powerful
computational resources and efficient algorithms to handle the vast amount of
data generated. Rapid and scalable computation is essential in this context.
It enables researchers to analyse massive datasets in real-time and identify
patterns and trends that would otherwise be missed. This is particularly
important in digital pathology, where the timely analysis of tissue images can
significantly impact patient outcomes. Possible directions for spatial
statistics research in digital pathology include developing novel algorithms
designed explicitly for high-dimensional data and exploring distributed
computing architectures that can leverage cloud-based resources for efficient
data processing. Ultimately, the ability to rapidly and efficiently analyse
spatial data in digital pathology can greatly enhance our ability to diagnose
and treat diseases, leading to better patient outcomes and improved public
health.

## References

  * Baddeley et al. 2017 Baddeley, A., A. Hardegen, T. Lawrence, R. K. Milne, G. Nair, and S. Rakshit (2017).  On two-stage Monte Carlo tests of composite hypotheses.  Computational Statistics & Data Analysis 114, 75–87.
  * Baddeley et al. 2015 Baddeley, A., E. Rubak, and R. Turner (2015).  Spatial Point Patterns: Methodology and Applications with R.  Chapman & Hall Interdisciplinary Statistics Series. CRC Press, Boca Raton, Florida.
  * Berman and Diggle 1989 Berman, M. and P. J. Diggle (1989).  Estimating weighted integrals of the second-order intensity of a spatial point process.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 51(1), 81–92.
  * Besag 1977 Besag, J. (1977).  Contribution to the discussion of dr Ripley’s paper.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 39, 193–195.
  * Bithell 1990 Bithell, J. F. (1990).  An application of density estimation to geographical epidemiology.  Statistics in Medicine 9(6), 691–701.
  * Chen et al. 2021 Chen, B., S. Cherie’R, E. T. McKinley, A. J. Simmons, M. A. Ramirez-Solano, X. Zhu, N. O. Markham, C. N. Heiser, P. N. Vega, A. Rolong, et al. (2021).  Differential pre-malignant programs and microenvironment chart distinct paths to malignancy in human colorectal polyps.  Cell 184(26), 6262–6280.
  * Chen et al. 2020 Chen, Z., I. Soifer, H. Hilton, L. Keren, and V. Jojic (2020).  Modeling multiplexed images with spatial-lda reveals novel tissue microenvironments.  Journal of Computational Biology 27(8), 1204–1218.
  * Chiu et al. 2013 Chiu, S. N., D. Stoyan, W. S. Kendall, and J. Mecke (2013).  Stochastic Geometry and Its Applications (Third ed.).  Wiley Series in Probability and Statistics. John Wiley & Sons, Chichester.
  * Cox 1977 Cox, D. R. (1977).  The role of significance tests.  Scandinavian Journal of Statistics 4(2), 49–70.
  * Cronie and van Lieshout 2018 Cronie, O. and M. van Lieshout (2018, 02).  A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions.  Biometrika 105(2), 455–462.
  * Cronie and van Lieshout 2015 Cronie, O. and M. N. M. van Lieshout (2015).  Summary statistics for inhomogeneous marked point processes.  Annals of the Institute of Statistical Mathematics, 1–24.
  * Daley and Vere-Jones 2007 Daley, D. and D. Vere-Jones (2007).  An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure (Second ed.).  Springer-Verlag, New York.
  * Dao and Genton 2014 Dao, N. A. and M. G. Genton (2014).  A Monte Carlo-adjusted goodness-of-fit test for parametric models describing spatial point patterns.  Journal of Computational and Graphical Statistics 23(2), 497–517.
  * Davies and Baddeley 2018 Davies, T. M. and A. Baddeley (2018).  Fast computation of spatially adaptive kernel estimates.  Statistics and Computing 28(4), 937–956.
  * Davies and Hazelton 2010 Davies, T. M. and M. L. Hazelton (2010).  Adaptive kernel estimation of spatial relative risk.  Statistics in Medicine 29(23), 2423–2437.
  * Diggle et al. 2005 Diggle, P., P. Zheng, and P. Durr (2005).  Nonparametric estimation of spatial segregation in a multivariate point process: Bovine tuberculosis in Cornwall, UK.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 54(3), 645–658.
  * Diggle 1985 Diggle, P. J. (1985).  A kernel method for smoothing point process data.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 34(2), 138–147.
  * Diggle 2003 Diggle, P. J. (2003).  Statistical Analysis of Spatial Point Patterns (Second ed.).  Chapman & Hall Monographs on Statistics & Applied Probability. CRC Press, Boca Raton, Florida.
  * Diggle 2013 Diggle, P. J. (2013).  Statistical Analysis of Spatial and Spatio-Temporal Point Patterns (Third ed.).  Chapman & Hall Monographs on Statistics & Applied Probability. CRC Press, Boca Raton, Florida.
  * Dries et al. 2021 Dries, R., J. Chen, N. Del Rossi, M. M. Khan, A. Sistig, and G.-C. Yuan (2021).  Advances in spatial transcriptomic data analysis.  Genome Research 31(10), 1706–1718.
  * Elhanani et al. 2023 Elhanani, O., R. Ben-Uri, and L. Keren (2023).  Spatial profiling technologies illuminate the tumor microenvironment.  Cancer Cell.
  * Gao et al. 2021 Gao, S., W. Zhang, N. Yan, M. Li, X. Mu, H. Yin, and J. Wang (2021).  The impact of stat3 and phospho-stat3 expression on the prognosis and clinicopathology of ovarian cancer: a systematic review and meta-analysis.  Journal of Ovarian Research 14(1), 1–18.
  * Gelfand et al. 2010 Gelfand, A. E., P. Diggle, P. Guttorp, and M. Fuentes (Eds.) (2010).  Handbook of spatial statistics.  Chapman & Hall Handbooks of Modern Statistical Methods. CRC Press, Boca Raton, Florida.
  * Geuenich et al. 2021 Geuenich, M. J., J. Hou, S. Lee, S. Ayub, H. W. Jackson, and K. R. Campbell (2021).  Automated assignment of cell identity from single-cell multiplexed imaging and proteomic data.  Cell Systems 12(12), 1173–1186.
  * Harris et al. 2022 Harris, C. R., E. T. McKinley, J. T. Roland, Q. Liu, M. J. Shrubsole, K. S. Lau, R. J. Coffey, J. Wrobel, and S. N. Vandekar (2022).  Quantifying and correcting slide-to-slide variation in multiplexed immunofluorescence images.  Bioinformatics 38(6), 1700–1707.
  * Ho and Chiu 2006 Ho, L. P. and S. N. Chiu (2006).  Testing the complete spatial randomness by Diggle’s test without an arbitrary upper limit.  Journal of Statistical Computation and Simulation 76(7), 585–591.
  * Hochberg 1988 Hochberg, Y. (1988, 12).  A sharper Bonferroni procedure for multiple tests of significance.  Biometrika 75(4), 800–802.
  * Illian et al. 2008 Illian, J., P. Penttinen, H. Stoyan, and D. Stoyan (2008).  Statistical Analysis and Modelling of Spatial Point Patterns.  Statistics in Practice. Wiley.
  * Johnson et al. 2021 Johnson, A. M., J. M. Boland, J. Wrobel, E. K. Klezcko, M. Weiser-Evans, K. Hopp, L. Heasley, E. T. Clambey, K. Jordan, R. A. Nemenoff, et al. (2021).  Cancer cell-specific major histocompatibility complex ii expression as a determinant of the immune infiltrate organization and function in the nsclc tumor microenvironment.  Journal of Thoracic Oncology 16(10), 1694–1704.
  * Jordan et al. 2020 Jordan, K. R., M. J. Sikora, J. E. Slansky, A. Minic, J. K. Richer, M. R. Moroney, J. Hu, R. J. Wolsky, Z. L. Watson, T. M. Yamamoto, J. C. Costello, A. Clauset, K. Behbakht, T. R. Kumar, and B. G. Bitler (2020).  The capacity of the ovarian cancer tumor microenvironment to integrate inflammation signaling conveys a shorter disease-free interval.  Clinical Cancer Research 26(23), 6362–6373.
  * Kelsall and Diggle 1995 Kelsall, J. E. and P. J. Diggle (1995).  Kernel estimation of relative risk.  Bernoulli 1(1/2), 3–16.
  * Kim et al. 2022 Kim, J., S. Rustam, J. M. Mosquera, S. H. Randell, R. Shaykhiev, A. F. Rendeiro, and O. Elemento (2022).  Unsupervised discovery of tissue architecture in multiplexed imaging.  Nature methods, 1–9.
  * Liang and Zeger 1986 Liang, K.-Y. and S. L. Zeger (1986).  Longitudinal data analysis using generalized linear models.  Biometrika 73(1), 13–22.
  * van Lieshout and Baddeley 1996 van Lieshout, M. N. M. and A. J. Baddeley (1996).  A nonparametric measure of spatial interaction in point patterns.  Statistica Neerlandica 50(3), 344–361.
  * van Lieshout and Baddeley 1999 van Lieshout, M. N. M. and A. J. Baddeley (1999).  Indices of dependence between types in multivariate point patterns.  Scandinavian Journal of Statistics 26(4), 511–532.
  * Loader 1999 Loader, C. (1999).  Local Regression and Likelihood.  Springer, New York.
  * Lotwick and Silverman 1982 Lotwick, H. W. and B. W. Silverman (1982).  Methods for analysing spatial processes of several types of points.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 44(3), 406–413.
  * McCullagh and Nelder 1989 McCullagh, P. and J. A. Nelder (1989).  Generalized Linear Models (Second ed.).  Chapman & Hall Monographs on Statistics and Applied Probability. CRC Press, Boca Raton, Florida.
  * Møller and Waagepetersen 2004 Møller, J. and R. P. Waagepetersen (2004).  Statistical Inference and Simulation for Spatial Point Processes.  Chapman & Hall Monographs on Statistics & Applied Probability. CRC Press, Boca Raton, Florida.
  * Mrkvička et al. 2021 Mrkvička, T., J. Dvořák, J. A. González, and J. Mateu (2021).  Revisiting the random shift approach for testing in spatial statistics.  Spatial Statistics 42, 100430.
  * Mrkvička et al. 2017 Mrkvička, T., M. Myllymäki, and U. Hahn (2017).  Multiple Monte Carlo testing, with applications in spatial point processes.  Statistics and Computing 27(5), 1239–1255.
  * Mrkvička et al. 2020 Mrkvička, T., M. Myllymäki, M. Jílek, and U. Hahn (2020).  A one-way ANOVA test for functional data with graphical interpretation.  Kybernetika 56(3), 432–458.
  * Myllymäki et al. 2017 Myllymäki, M., T. Mrkvička, P. Grabarnik, H. Seijo, and U. Hahn (2017).  Global envelope tests for spatial processes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79(2), 381–404.
  * Nadaraya 1989 Nadaraya, E. A. (1989).  Nonparametric estimation of probability densities and regression curves.  Kluwer, Dordrecht.
  * Ripley 1977 Ripley, B. D. (1977).  Modelling spatial patterns (with discussion).  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 39(2), 172–212.
  * Ripley 1988 Ripley, B. D. (1988).  Statistical Inference for Spatial Processes.  Cambridge University Press. Cambridge.
  * Ripley and Rasson 1977 Ripley, B. D. and J. P. Rasson (1977).  Finding the edge of a Poisson forest.  Journal of Applied Probability 14(3), 483–491.
  * Scott 2015 Scott, D. W. (2015).  Multivariate density estimation: theory, practice, and visualization (Second ed.).  Wiley series in probability and statistics. John Wiley & Sons.
  * Seal et al. 2022 Seal, S., T. Vu, T. Ghosh, J. Wrobel, and D. Ghosh (2022).  DenVar: density-based variation analysis of multiplex imaging data.  Bioinformatics Advances 2(1), vbac039.
  * Steinhart et al. 2021 Steinhart, B., K. R. Jordan, J. Bapat, M. D. Post, L. W. Brubaker, B. G. Bitler, and J. Wrobel (2021).  The spatial context of tumor-infiltrating immune cells associates with improved ovarian cancer survival.  Molecular Cancer Research 19(12), 1973–1979.
  * Sucaet and Waelput 2014 Sucaet, Y. and W. Waelput (2014).  Digital Pathology (1 ed.).  Springer Briefs in Computer Science. Springer, New York.
  * Terrell 1990 Terrell, G. R. (1990).  The maximal smoothing principle in density estimation.  Journal of the American Statistical Association 85(410), 470–477.
  * van Dam et al. 2022 van Dam, S., M. J. Baars, and Y. Vercoulen (2022).  Multiplex tissue imaging: Spatial revelations in the tumor microenvironment.  Cancers 14(13), 3170.
  * Vu et al. 2022 Vu, T., J. Wrobel, B. G. Bitler, E. L. Schenk, K. R. Jordan, and D. Ghosh (2022).  SPF: A spatial and functional data analytic approach to cell imaging data.  PLOS Computational Biology 18(6), 1–21.
  * Wilson et al. 2022 Wilson, C., A. C. Soupir, R. Thapa, J. Creed, J. Nguyen, C. M. Segura, T. Gerke, J. M. Schildkraut, L. C. Peres, and B. L. Fridley (2022).  Tumor immune cell clustering and its association with survival in african american women with ovarian cancer.  PLoS Computational Biology 18(3), e1009900.
  * Wilson et al. 2021 Wilson, C. M., O. E. Ospina, M. K. Townsend, J. Nguyen, C. Moran Segura, J. M. Schildkraut, S. S. Tworoger, L. C. Peres, and B. L. Fridley (2021).  Challenges and opportunities in the statistical analysis of multiplex immunofluorescence data.  Cancers 13(12), 3031.
  * Wrobel and Ghosh 2022 Wrobel, J. and T. Ghosh (2022).  VectraPolarisData: Vectra Polaris and Vectra 3 multiplex single-cell imaging data.  Bioconductor R package version 1.0.
  * Xun et al. 2023 Xun, Z., X. Ding, Y. Zhang, B. Zhang, S. Lai, D. Zou, J. Zheng, G. Chen, B. Su, L. Han, et al. (2023).  Reconstruction of the tumor spatial microenvironment along the malignant-boundary-nonmalignant axis.  Nature Communications 14(1), 933.
  * Zhang et al. 2022 Zhang, W., I. Li, N. E. Reticker-Flynn, Z. Good, S. Chang, N. Samusik, S. Saumyaa, Y. Li, X. Zhou, R. Liang, et al. (2022).  Identification of cell types in multiplexed in situ images by combining protein expression and spatial information using celesta.  Nature Methods 19(6), 759–769.