[1]Saul Fuster

[1]Department of Electrical Engineering and Computer Science, University of
Stavanger, Stavanger, 4021, Norway

2]Department of Pathology, Stavanger University Hospital, Stavanger, 4011,
Norway

3]Department of Chemistry, Stavanger University Hospital, Stavanger, 4021,
Norway

# NMGrad: Advancing Histopathological Bladder Cancer Grading with Weakly
Supervised Deep Learning

<EMAIL_ADDRESS>Umay Kiraz<EMAIL_ADDRESS>Trygve Eftestøl
<EMAIL_ADDRESS>Emiel A.M. Janssen
<EMAIL_ADDRESS>Kjersti Engan<EMAIL_ADDRESS>*  [
[

###### Abstract

The most prevalent form of bladder cancer is urothelial carcinoma,
characterized by a high recurrence rate and substantial lifetime treatment
costs for patients. Grading is a prime factor for patient risk stratification,
although it suffers from inconsistencies and variations among pathologists.
Moreover, absence of annotations in medical imaging difficults training deep
learning models. To address these challenges, we introduce a pipeline designed
for bladder cancer grading using histological slides. First, it extracts
urothelium tissue tiles at different magnification levels, employing a
convolutional neural network for processing for feature extraction. Then, it
engages in the slide-level prediction process. It employs a nested multiple
instance learning approach with attention to predict the grade. To distinguish
different levels of malignancy within specific regions of the slide, we
include the origins of the tiles in our analysis. The attention scores at
region level is shown to correlate with verified high-grade regions, giving
some explainability to the model. Clinical evaluations demonstrate that our
model consistently outperforms previous state-of-the-art methods.

###### keywords:

Computational Pathology, Deep Learning, Grading, Multiscale, Urothelial
Carcinoma, Weakly Supervised Learning

##  1 Introduction

Bladder cancer, a prevalent urological malignancy, poses significant clinical
challenges in terms of both diagnosis and prognosis [1]. Non-muscle-invasive
bladder cancer (NMIBC) accounts for approximately 75% of the newly diagnosed
cases of urothelial carcinoma. NMIBC is particularly known for its variable
outcomes, necessitating accurate and consistent grading for optimal patient
management [2]. The 2022 edition of the European Association of Urology
guidelines on NMIBC recommends a stratification of patients into risk groups
based on the risk of progression to muscle-invasive disease [3]. Grade, stage,
and various other factors contribute to the risk. Precise risk assessment is
vital in the management of NMIBC, since treatment strategies not only rely on
the presence of muscle invasion.

Grading is based on assessment of the cellular morphology abnormalities of
urothelial tissue. In 2004, the WHO introduced a grading classification system
(WHO04) for NMIBC based on histological features. WHO04 encompasses three
categories: papillary urothelial neoplasm of low malignant potential (PUNLMP),
non-invasive papillary carcinoma low grade (LG), and high grade (HG), ranging
from lower to higher malignancy respectively [4]. HG is related to lower
differentiation, loss of polarity, and pleomorphic nuclei, among others. The
intricate evaluation of heterogeneous scenarios contributes to significant
inter- and intraobserver variability. Disparities potentially lead to
misclassification, and consequently, to inappropriate treatment decisions [5].
WHO04 grading system was subsequently retained in the updated 2016/2022 WHO
classifications [6]. However, according to several multi-institutional
analyses of individual patient data, the proportion of tumors classified as
PUNLMP (WHO 2004/2016) has markedly decreased to very low levels in the last
decade. This trend resulted in suggestions to reassess PUNLMP tumors as LG [7,
8]. Therefore, the upcoming grading system will reasonably undergo a
modification, shifting towards the inclusion of only LG and HG categories. In
recent years, the integration of deep learning techniques within the field of
computational pathology (CPATH) has offered promising avenues for enhancing
the precision of computer-aided diagnosis (CAD) systems and elucidating
discrepancies among pathologists [9, 10]. Consequently, the synergy between
histological expertise and CAD technologies is vital for accurate grading
assessments.

CPATH is the field of pathology that leverages the potential of CAD systems to
thoroughly analyze high-resolution digital images known as Whole Slide Images
(WSIs) for diverse diagnostic and prognostic purposes [11]. WSIs are produced
by slide scanners and pre-stored at different magnification levels, emulating
the functionality of physical microscopes. Lower magnification is suited for
tissue-level morphology examination, while higher magnification for cell-level
scrutiny [12]. WSIs are characterized by their substantial size, which can
introduce adversarial noise. Bladder cancer WSIs present unique challenges due
to their disorganized nature and the presence of diagnostically non-relevant
tissue. These slides often include artifacts, such as cauterized or stretched
tissue [13]. Moreover, tissue such as blood, muscle and stroma are less
informative for grading a tumor. Therefore, the absence of annotations
presents a significant challenge for identifying regions of interest (ROIs)
[9]. It is crucial to distinguish between region-based labels (e.g., tissue
type, grade) and WSI-based labels, including follow-up information and overall
patient grade. Grade exemplifies a label that encompasses both perspectives
[14, 15]. While clinical reports assign the worst grade observed to a patient,
a WSI may exhibit diverse urothelium regions with normal, LG, and HG. This
dual nature underscores the complexity of label interpretation when
considering both the medical, WSI-focused perspective and the more technical
perspective involving regional data analysis and processing.

CPATH is amid a transformative era, aiming to reshape the landscape of digital
pathology as we know it [16]. Among the diverse practices in the field,
imaging methodologies rooted in convolutional neural networks (CNNs) have
emerged as the foundation of feature extraction from histological images [11].
These deep learning networks possess a remarkable capacity to automatically
discern morphological and cellular patterns within WSIs [17, 18, 19, 20, 21,
22]. Ultimately, CNNs contribute to more precise and timely clinical
decisions. However, training deep learning models in CPATH presents challenges
when only WSI-level labels are available, lacking region-based annotations [9,
23, 24, 25]. To address these, weakly supervised learning techniques, like
attention-based methods and multiple instance learning (MIL), are employed.
However, MIL methods can be susceptible to individual instances dominating the
weighted aggregation of the WSI representation [26, 27, 28]. In the context of
WSIs, the tissue is distributed across the slide, for which regions typically
present similar features within and pathologists pinpoint ROIs with crucial
information [29]. Specifically while grading, situations may arise where
multiple instances in close proximity exhibit HG characteristics, while other
regions may concurrently display LG attributes. As a result, constraining
instances to specific regions enhances our understanding of the diverse
features within WSIs. Consequently, a conventional MIL architecture approach
might not be appropriate, because there is a susceptibility to information
leakage between regions. A model accommodating for the nested structure of
WSIs, wherein tissues are part of a region, and regions, in turn, belong to a
WSI, may more effectively capture the clinical WSI-level grading label [30,
31, 32].

In this study, we introduce a novel pipeline for grading NMIBC using
histological slides, referred to as nested multiple grading (NMGrad). The
proposed solution starts by tissue segmentation of the WSI, separating
urothelium from other tissue types. The next step categorizes extracted tiles
of urothelium areas into location-dependent regions for predicting the
patient’s WHO04 grade. We implemented a weakly supervised learning framework
using attention mechanisms, and a nested aggregation architecture for ROI
differentiation. Our method offers an innovative approach for generating
diagnostic suggestions, with generated heatmaps for highlighting tiles and
ROIs independently.

##  2 Related Work

Numerous studies in the domain of computational pathology for bladder cancer
diagnostics have emerged in recent years [33]. In Wetteland et al. [34], a
pipeline for grading NMIBC is introduced. This pipeline identifies relevant
areas in the WSIs and predicts cancer grade by considering individual tile
predictions and applying a decision threshold to determine the overall patient
prediction. Their results demonstrate promising performance, with potential
benefits for patient care. In Zheng et al. [35], the authors focus on the
development of deep learning-based models for bladder cancer diagnosis and
predicting overall survival in muscle-invasive bladder cancer patients. They
introduce two deep learning models for diagnosis and prognosis respectively.
They show that their presented algorithm outperformed junior pathologists. In
Jansen et al. [36], the authors propose a fully automated detection and
grading network based on deep learning to enhance NMIBC grading
reproducibility. The study employs a U-Net-based segmentation network to
automatically detect urothelium, followed by a VGG16 CNN network for
classification. Their findings demonstrate that the automated classification
achieves moderate agreement with consensus and individual grading from a group
of three senior uropathologists. Spyridonos et al. [37] investigate the
effectiveness of support vector machines and probabilistic neural networks for
urinary bladder tumor grading. The results indicate that both SVM and PNN
models achieve a relatively high overall accuracy, with nuclear size and
chromatin cluster patterns playing key roles in optimizing classification
performance.

Zhang et al. [38] address a common limitation of interpretability in CAD
methods. To tackle this, they introduce MDNet, a novel approach that
establishes a direct multimodal mapping between medical images and diagnostic
reports. This framework consists of an image model and a language model.
Through experiments on pathology bladder cancer images and diagnostic reports,
MDNet demonstrates superior performance compared to comparative baselines.
Zhang et al. [39] propose a method that leverages deep learning to automate
the diagnostic reasoning process through interpretable predictions. Using a
dataset of NMIBC WSIs, the study demonstrates that their method achieves
diagnostic performance comparable to that of 17 uropathologists.

Two critical challenges we identified include summarizing information from
local image features into a WSI representation, and the scarcity of annotated
datasets. Effectively translating detailed local information to the WSI level
is complex, particularly in tasks like grading NMIBC. Moreover, the limited
availability of well-annotated datasets hinders the development and evaluation
of robust models. To tackle these issues, weakly supervised methods have
emerged as a standout tool in CPATH [9]. While weakly supervised methods are
widespread, some studies still rely on annotations and supervised learning.
However, there is a growing consensus for the future of CPATH to predominantly
embrace weakly supervised approaches. This shift is driven by the
impracticality of obtaining detailed annotations for large datasets covering
various cancers and tasks. Among the various weakly supervised methods,
attention-based MIL (AbMIL), a popular instance aggregation method, exploits
attention mechanisms in order to mitigate the uncertainty from individual
instances and enhance interpretability [40, 41]. AbMIL bridges the gap between
limited supervision and the spatial details necessary for accurate analysis
and explainability. An evolution of MIL model architectures relies on the
arrangement of the data within bags, where instances are further subdivided
into finer groups. This concept is referred to as nesting [30, 31]. Nested
architectures preserve a sense of localization or categorization by
selectively processing data instances within individual subgroups.
Subsequently, they aggregate summarized information from the subgroups into a
final bag representation.

In our work, we aim to bridge the gap between non-annotated datasets, weakly
supervised methods and the intrinsic categorization of WSI data. Therefore, we
leverage the nested MIL with attention mechanisms (NMIA) model architecture we
proposed in Fuster et. al. [32], for accurate and interpretable NMIBC grading.
Finally, in order to overcome the lack of annotations for defining the tissue
of interest, a tissue segmentation algorithm TRI-25x-100x-400x was proposed by
our research group in [42]. More recent works from our research group on
tissue segmentation have found adoption within the scientific literature [43,
44]. The utilization of this segmentation algorithm offers the opportunity to
extract tiles specifically from the urothelium, contributing to a refined and
targeted extraction process.

##  3 Data Material

The dataset comprises a total of 300 digital whole-slide images (WSIs) derived
from 300 patients diagnosed with NMIBC, from the Department of Pathology,
Stavanger University Hospital (SUH) [14, 45]. The glass slides were digitized
using a Leica SCN400 slide scanner, saved in the vendor-specific SCN file
format. Collected over the period spanning 2002 to 2011, this dataset
encompasses all risk group cases of non-muscle invasive bladder cancers. The
biopsies were processed through formalin-fixation and paraffin-embedding, and
subsequently, 4µm thick sections were prepared and stained using Hematoxylin,
Eosin, and Saffron (HES). Furthermore, all WSIs underwent meticulous manual
quality checks, ensuring the inclusion of only high-quality slides with
minimal or no blur. Due to the cauterization process used in the removal of
NMIBC, some slides may exhibit areas with burned and damaged tissues. All WSIs
originate from the same laboratory, resulting in relatively consistent
staining color across the dataset.

All WSIs were graded by an expert uropathologist in accordance with the WHO04
classification system, as either LG or HG, thus providing slide-level
diagnostic information. However, the dataset lacks region-based annotations
pinpointing the precise areas of LG or HG regions within the WSI.
Consequently, the dataset is considered weakly labeled. For WSIs labeled as
LG, at least one LG region is expected, with the possibility of presenting
non-cancerous tissue in other regions. As for HG slides, at least one region
should display HG tissue, while other regions may exhibit a LG appearance or
non-cancerous tissue. Given the absence of alternative gold standards, we are
compelled to continue utilizing a grading assessment that may have limitations
for training and evaluating our algorithms. The dataset employed in this study
was divided into three subsets: 220 WSI/patients for training, 30 for
validation, and 50 for test. The split employed ensures that each subset
maintains the same proportional representation of diagnostic outcomes. This
stratification encompassed factors such as WHO04 grading, cancer stage,
recurrence, and disease progression to best mirror the diversity of the
original data material. The distribution of LG and HG WSIs within each dataset
is detailed in Table 1 for reference.

Within a subset of the test set, denoted as Test${}_{\text{ANNO}}\in$ Test, 14
WSIs contain either one or two annotated regions of confirmed LG or HG tissue,
verified by an expert uropathologist. It is noteworthy that not all regions
were annotated. The labels of these regions correspond to the associated weak
label of the WSI.

Table 1: Number of WSIs and corresponding grade per subset.

| Low-grade | High-grade
---|---|---
Train | 124 (0) | 96 (0)
Validation | 17 (0) | 13 (0)
Test | 28 (7) | 22 (7)

00footnotetext: The number between parenthesis corresponds to slides
containing some annotations giving region-based labels.

##  4 Methods

We propose NMGrad, a pipeline that begins with a tissue segmentation algorithm
for extracting urothelium tissue. Subsequently, the urothelium is divided into
localized regions. Thereafter we employ a weakly supervised learning method to
predict tumor grade from the segmented urothelium regions. We exploit the
sense of region locations by adopting a nested architecture with attention,
NMIA [32]. The rationale behind employing this structured data arrangement
analysis is to identify relevant instances and regions within the WSI. The
attention mechanism and the nested bags/regions also contributes to a more
precise and insightful analysis of the data. An overview of NMGrad can be
visualized in Fig. 1.

Figure 1: NMGrad pipeline. Initially, we apply a tissue segmentation algorithm
for ROI extraction. Then, we pinpoint diagnostically significant urothelium
areas within WSIs. Subsequently, we split the urothelium mask into regions,
based on proximity and size, and extract tile triplets. In a hierarchical
fashion, we further transform these triplets within their corresponding
regions into region feature embeddings using an attention-based aggregation
method. All the region representations are then consolidated into a
comprehensive WSI-level representation through a weight-independent attention
module. Finally, this WSI feature embedding is input into the WHO04 grading
classifier in order to produce accurate WSI grade predictions.

###  4.1 Automatic Tissue Segmentation & Region Definition

We utilize the tissue segmentation algorithm introduced by Wetteland et al.
[42] to automatically generate tissue type masks, facilitating the subsequent
extraction of tiles. We define triplets $\mathcal{T}$ of tissue, which consist
of a set of three tiles at various magnifications levels, namely 25x, 100x and
400x. An example is shown in Fig. 2. We used a tile size of $128\times 128$
for all magnification levels. Triplets are formed to maintain consistency,
ensuring that the center pixel in every tile accurately represents the same
physical point. The tissue segmentation algorithm works at tile level, and
classifies all triplets $\mathcal{T}$ in the WSI as $y\in{\cal{Y}}=$
{urothelium, lamina propria, muscle, blood, damage, background}. As grading
relies on urothelium alone, we utilize the urothelium mask for defining valid
areas for tile extraction, as described in [46]. In this work, various
magnifications are explored for defining the model’s input, either using mono-
scale MONO (400x), di-scale DI (400x, 100x) or tri-scale TRI (400x, 100x,
25x). We employed 400x magnification to establish a tight grid of tiles for
data extraction purposes. These sets of tiles are later fed to the grading
models, where each magnification tile is processed by its respective weight-
independent CNN. To preserve the sense of location within the image, we define
regions. This results in the following stratified division of data: all
urothelium in the WSI, scattered regions of urothelium, and finally,
individual tiles of urothelium.

Figure 2: We obtain comprising sets of three tiles at different magnification
levels named triplets $\mathcal{T}$, enabling detailed examination. Tile
triplets demonstrate regions associated with low- and high-grade features.

####  4.1.1 Region Definition

For defining regions out of the extracted urothelium tiles, we define blobs of
tiles $\text{URO}\textsubscript{BLOB}\subseteq\text{URO}$. UROBLOB is formed
when tiles are 8-connected, and this joint set of tiles is the representation
of a region
$\text{URO}\textsubscript{BLOB}=\\{\mathcal{T}_{1},\mathcal{T}_{2}....\\}$. A
region is eligible for inclusion if the number of tiles $N_{B}$ is higher than
a threshold number $T_{\text{LOWER}}$. Any blob with $N_{B}<T_{\text{LOWER}}$
tiles is discarded, along with the tiles within. As NMIBC WSI can contain
large tissue bundles, resulting in sizable blobs, we also define an upper
limit threshold $T_{\text{UPPER}}$. For blobs where $N_{B}\geq
T_{\text{UPPER}}$, we split the region into several sub regions for more
detailed analysis. We apply KMeans clustering over the coordinates of the
tiles $x,y$ within the blob for location sense, defining the number of
clusters as $N_{C}=\lceil N_{B}/T_{\text{UPPER}}\rceil$. This results in joint
regions within a bundle of tissue of consistent size, as observed in regions
5-8 in Fig. 3.

Figure 3: Region definition. Urothelial tissue within a WSI is eligible for
tile extraction. Blobs of tiles are formed, and blobs smaller than a threshold
$T_{\text{LOWER}}$ are discarded. From the remaining blobs, any smaller than
$T_{\text{UPPER}}$ are kept and defined as a region. For blobs bigger than
$T_{\text{UPPER}}$, the blobs is subdivided into smaller pieces using the
location of individual tiles within and KMeans clustering. The obtained
clusters are designated as regions.

###  4.2 Multiple Instance Learning in a WSI context

Multiple instance learning (MIL) is a weakly supervised learning method where
unlabeled instances are grouped into bags with known labels [26, 27]. A
dataset $\mathcal{X},\mathcal{Y}=\left\\{(\mathbf{X}^{i},y^{i}),\forall
i=1,...,N\right\\}$ is formed of pairs of sample sets $\mathbf{X}$ and their
corresponding labels $y$, where $i$ denotes a bag index. In the context of
WSI, the bag can be one patient or one WSI or one region. In a conventional
MIL data arrangement, we consider the bag $\mathbf{X}$ to be one WSI
consisting of instances $\mathbf{x}_{l}$:

$\mathbf{X}=\left\\{\mathbf{x}_{l},\forall l=1,...,L\right\\}$ (1)

where $L$ is the number of instances in the bag. In our study, $\mathbf{x}$
refers to individual tiles in a set of extracted tiles from a patient slide
$\mathbf{X}$. A feature extractor $G_{\theta}:{\mathcal{X}}\to\mathcal{H}$
transforms image tiles, $\mathbf{x}_{l}$, into low-dimensional feature
embeddings, $\mathbf{h}_{l}$. At this point, the bag structure previously
formed remains intact, as instances have been simply transformed. Given a
label $y$ for a WSI, the training objective of the model is to predict the
grade observed in the WSI. However, to deduce the specific region(s) within a
WSI that leads to the patient’s diagnosis of either LG or HG is of utmost
importance. This entails the model’s capability to discern and highlight the
critical areas within the WSI that play a pivotal role in the diagnosis. In
order to accomplish this goal, we adopted attention-based multiple instance
learning (AbMIL) as our MIL framework, using attention-based aggregation as
shown in Fig. 1. An attention score $a_{i}$ for a feature embedding
$\mathbf{h}_{i}$ can be calculated as:

$a_{i}=\frac{\mathrm{exp}\\{\mathbf{w}^{\top}(\mathrm{tanh}(\mathbf{Vh_{i}^{\top}})\odot\mathrm{sigm}(\mathbf{Uh_{i}^{\top})})\\}}{\sum_{l=1}^{L}\mathrm{exp}\\{\mathbf{w}^{\top}(\mathrm{tanh}(\mathbf{Vh_{l}^{\top}})\odot\mathrm{sigm}(\mathbf{Uh_{l}^{\top})})\\}}$
(2)

where $\mathbf{w}\in\mathbb{R}^{L\times 1}$, $\mathbf{V}\in\mathbb{R}^{L\times
M}$ and $\mathbf{U}\in\mathbb{R}^{L\times M}$ are trainable parameters and
$\odot$ is an element-wise multiplication. Furthermore, the hyperbolic tangent
$\mathrm{tanh}(\cdot)$ and sigmoid $\mathrm{sigm}(\cdot)$ are included to
introduce non-linearity for learning complex applications. The benefit of
attention modules extends beyond interpretability for understanding the
model’s decision-making process, as it also grants enhanced predictive
capabilities by prioritizing salient features. This is because attention
scores directly influence the forward propagation of the model, allowing it to
focus on the most relevant and informative regions within the data. Once the
attention scores $\mathbf{A}$ are obtained, we obtain the patient prediction
$\hat{y}$ using a patient classifier $\Psi_{\rho}$ as:

$\hat{y}=\Psi_{\rho}(\mathbf{H}^{a})=\Psi_{\rho}(\mathbf{A}\cdot\mathbf{H})=\Psi_{\rho}(\mathbf{A}\cdot
G_{\theta}(\mathbf{X}))$ (3)

####  4.2.1 Nested Multiple Instance Architecture

An evolution of the conventional AbMIL architecture defines levels of bags
within bags where only the innermost bags contain instances. This is referred
to as nested multiple instance with attention (NMIA) [32]. A bag-of-bags for a
WSI $\mathbf{H}_{\text{WSI}}$ contains a set of inner-bags, or regions,
$\mathbf{H}_{\text{REG},k}$:

$\mathbf{H}_{\text{WSI}}=\left\\{\mathbf{H}_{\text{REG},k},\forall
k=1,...,K\right\\}$ (4)

where the number of inner-bags $K$ varies between different WSI. Ultimately,
$\mathbf{H}_{\text{REG},k}$ contains instance-level representations
$\mathbf{h}_{\text{TILE},l}$ of tiles located within the physical region. This
serves to further stratify into clusters or regions, and accurately represent
the arrangement of the scattered data, where tiles belong to particular tissue
areas and these themselves to the WSI.

$\mathbf{h}_{\text{REG}}=\mathbf{A}_{\text{TILE}}\cdot\mathbf{H}_{\text{TILE}}=\mathbf{A}_{\text{TILE}}\cdot
G_{\theta}(\mathbf{X})$ (5)

Finally, the ultimate WSI representation $h_{\text{WSI}}$ is fed to the
classifier for obtaining the grade prediction leveraging the region
representations $\mathbf{H}_{\text{REG}}$ and the attention scores
$\mathbf{A}_{\text{REG}}$ obtained from those same representations:

$\hat{y}=\Psi_{\rho}(h_{\text{WSI}})=\Psi_{\rho}(\mathbf{A}_{\text{REG}}\cdot\mathbf{H}_{\text{REG}})$
(6)

##  5 Experiments

Within the scope of our experimental investigation, we systematically assessed
and compared the impact of diverse magnification levels and combinations of
them, namely MONO, DI and TRI-scale models. We further investigate the impact
of several weakly supervised aggregation techniques in the performance of our
deep learning model. These aggregation techniques vary in their ability to
distill valuable diagnostic insights from the data, namely mean, max and
attention-based. We put special emphasis in the comparison between a standard
AbMIL architecture and the nested model proposed in NMGrad. Finally, we
compare our solution to current state-of-the-art methods. The code is
available at https://github.com/Biomedical-Data-Analysis-Laboratory/GradeMIL.

We list the details of design choice of the models during training. VGG16 is
used as CNN backbone, with ImageNet pretrained weights [47]. In preliminary
experiments, we explored various architectures and found that VGG16 exhibited
favorable performance across multiple tasks on NMIBC WSIs. Stochastic gradient
descent (SGD) was set as optimizer with a learning rate of 1e-4, a batch size
of 128, and a total of 200 epochs, with 30 epoch for early stopping based on
the AUC score on the validation set. A total number of 5000 tiles per WSIs
were preemptively randomly sampled to be further sub-sampled during training.
Focal Tversky loss (FTL) is employed [48]. The Tversky index (TI) leverages
false predictions, emphasizing on recall in case of large class imbalance
tuning parameters $\alpha$ and $\beta$. TI is defined as:

$\text{TI}_{c}(\hat{y},y)=\frac{1+\sum_{i=1}^{N}\hat{y}_{i,c}y_{i,c}+\epsilon}{1+\sum_{i=1}^{N}\hat{y}_{i,c}y_{i,c}+\alpha\sum_{i=1}^{N}\hat{y}_{i,\hat{c}}y_{i,c}+\beta\sum_{i=1}^{N}\hat{y}_{i,c}y_{i,\hat{c}}+\epsilon}$
(7)

where $\hat{y}_{i,\hat{c}}=1-\hat{y}_{i,c}$ and $y_{i,\hat{c}}=1-y_{i,c}$ are
the probability that sample $i$ is not of class $c\in\mathcal{C}$. $\epsilon$
is used for numerical stability, preventing zero division operations. FTL
employs another parameter $\gamma$ for leveraging training examples hardship:

$\text{FTL}_{c}(\hat{y},y)=\sum(1-\text{TI}_{c}(\hat{y},y))^{1/\gamma}$ (8)

##  6 Results & Discussion

The utilization of our proposed pipeline NMGrad using TRI-scale input emerged
as a standout performer, as of Table 2. TRI-scale models showcased the ability
to capture relevant patterns across different magnifications, substantially
enhancing the overall grading accuracy. Comparing the effects of scales on the
model performance shows a larger gap in performance between MONO and TRI
models than structuring the processing of data using NMGrad or a standard
AbMIL architecture. Furthermore, our exploration of aggregation techniques
extended to mean and max aggregation methods, which do not include in-built
attention mechanisms, yielded less promising outcomes. The absence of
attention mechanisms rendered these techniques less effective in capturing
nuanced features, underscoring the significance of attention mechanisms.

Table 2: Test performance for various aggregation techniques in weakly supervised learning. We provide the average of five runs, with the standard deviation shown in parentheses. The table presents the different approaches employed in the field, including aggregation techniques that involve considering spatial separation of the instances. We also explore the use of multiple magnification levels, considering 400x as the foundation for all magnification analysis. We also shown the results from other bladder cancer grading works, although in other datasets. Model | Accuracy | Precision | Recall | F1 Score | $\kappa$ | AUC
---|---|---|---|---|---|---
AbMILMONO | 0.68(0.07) | 0.71(0.07) | 0.68(0.06) | 0.67(0.07) | 0.36(0.12) | 0.81(0.07)
AbMILDI | 0.79(0.09) | 0.80(0.10) | 0.78(0.09) | 0.78(0.09) | 0.57(0.18) | 0.85(0.13)
AbMILTRI | 0.82(0.07) | 0.82(0.07) | 0.82(0.07) | 0.82(0.07) | 0.64(0.14) | 0.91(0.04)
MEANTRI | 0.81(0.03) | 0.83(0.03) | 0.80(0.03) | 0.80(0.03) | 0.61(0.05) | 0.92(0.03)
MAXTRI | 0.80(0.06) | 0.80(0.06) | 0.78(0.06) | 0.79(0.07) | 0.58(0.13) | 0.85(0.06)
NMGradMONO | 0.68(0.09) | 0.71(0.08) | 0.69(0.08) | 0.68(0.09) | 0.37(0.16) | 0.80(0.06)
NMGradDI | 0.83(0.03) | 0.85(0.03) | 0.82(0.03) | 0.82(0.03) | 0.65(0.06) | 0.91(0.04)
NMGradTRI | 0.86(0.03) | 0.87(0.02) | 0.85(0.04) | 0.85(0.03) | 0.71(0.06) | 0.94(0.01)
Wetteland et al. [34] | 0.90(-) | 0.87(-) | 0.80(-) | 0.83(-) | - | -
Jansen et al. [36] | 0.74(-) | - | 0.71(-) | - | 0.48(0.14) | -
Zhang et al. [39] | 0.95(-) | - | - | - | - | 0.95(-)

NMIA architecture embedded in NMGrad, which employs attention mechanisms for
both tile and region aggregation, marked a substantial leap in performance in
comparison to mean and max aggregation. This strategy provided the strengths
of attention-based aggregation and ROI localization via a nested architecture.
The incorporation of attention mechanisms allowed the model to pinpoint and
emphasize critical visual cues within WSIs related to urothelial cell
differentiation, ultimately resulting in a notable enhancement in predictive
accuracy for bladder cancer grading. Finally, in a direct comparison to the
previous best performing model proposed by Wetteland [34], they implemented
weakly supervised learning in a naive manner, where all patches were assigned
a weak label. Predictions are made at the patch level, and the determination
of WSI-level prediction relies somewhat arbitrarily on the summation of patch
predictions, neglecting consideration of localized regions. Using our proposed
solution NMGradTRI, the solution aligns more closely with clinical
expectations, such as the presence of one or more regions indicative of HG if
the WSI is classified as HG, rather than scattered patches. Furthermore,
NMGradTRI capacity to learn attention scores offers interpretability, as
opposed to relying on ad-hoc rules for post-processing patches into a final
prediction. Ultimately, we obtained slightly better F1 score, with a trade-off
on accuracy. We have also adhered the works of Jansen [36] and Zhang [39] for
comparing the performance of state-of-the-art grading algorithms, although the
results correspond to their in-house cohorts, different from ours. The
development of our deep learning model NMGradTRI for predicting the grading of
bladder cancer represents a significant advancement in the realm of accurate
grading of bladder tumors.

In order to enhance the fidelity of binary decisions, we opted to introduce an
uncertainty spectrum, thereby introducing a third class. Given the output
predictions of test set WSIs, we define the uncertainty spectrum as
$[\mu_{\hat{y}(y=LG)}+\sigma_{\hat{y}(y=LG)},\mu_{\hat{y}(y=HG)}-\sigma_{\hat{y}(y=HG)}]$.
A plot illustrating the concept and WSI predictions is shown in Fig. 4. It was
observed that by excluding predictions falling within the uncertainty
spectrum, the overall F1 score increased to 0.89. This underscores the
potential utility of considering skepticism regarding non-confident
predictions for robust clinical interpretation.

Figure 4: Plot displaying the WSI predictions of the test set, with green
shading representing the LG confidence interval, red for HG, and gray denoting
the uncertainty interval. Additionally, a blue line depicts the regression
line fitting the predictions.

Due to the attention scores generated at the inference stage, we are able to
visualize a heatmap, as in Fig. 5. NMGradTRI demonstrated effectiveness in
correctly assessing individual tiles and ROIs, despite being trained solely on
patient-level weak labels. We consulted the generated heatmaps with
experienced pathologists for qualitative analysis. Results on the annotated
set of regions, Test${}_{\text{ANNO}}$, exhibited competence to discern LG and
HG regions, as of Table 3. Region attention scores were considered for direct
comparison to region prediction scores utilizing the classifier $\Psi_{\rho}$,
as values for both scores are restrained between 0 and 1. These two scores
were evaluated individually against the annotated areas. It was corroborated
that higher attention scores are associated with HG areas in high-grade WSIs.
However, the same does not apply for low-grade. For low-grade WSIs, we
observed a wide range of possibilities, where high-attention is spread across
regions. Ideally, one would anticipate a direct correlation between HG and
elevated attention, and conversely, a correlation between LG and reduced
attention. However, this correlation is primarily observable in the positive
class (HG), aligning with the inherent design of MIL, which is tailored for
identifying positive instances. In contrast to attention, we observed a high
degree of correspondence between the label of annotated ROIs and the output
region predictions.

Figure 5: Region-level attention score heatmaps. Example regions of annotations of low- and high-grade ROIs annotated by an uropathologist are compared to the output attention provided by the proposed model NMGrad, left to right respectively. The choice of annotated ROIs correspond to highest attention scores, for red and blue correspond to low and high attention correspondingly. We have included the WSI-level prediction score for reference. Table 3: Region-level prediction and attention correspondence on annotated areas, using NMGradTRI. We individually compare the degree of consensus of annotations with both the highest attributed attention within the WSI and the output region prediction of the classifier $\Psi_{\rho}$. | Accuracy | Precision | Recall | F1 Score
---|---|---|---|---
Attention | 0.76 | 0.81 | 0.69 | 0.75
Prediction | 0.89 | 0.83 | 0.91 | 0.87

We further investigated the correlation between region attention scores and
the output region predictions, as displayed in Fig. 6. We observe a
generalized reciprocity between low-grade having smaller attention scores and
prediction outputs, and vice versa. Moreover, the high-grade range of values
is more limited to a lower range compared to the low-grade. For instance, low-
grade areas with predictions rounding zero show the broadest range of values.
This observation aligns with the earlier statement, wherein the positive class
typically exhibits a more focused distribution of attention scores,
predominantly linked with positive HG instances. In contrast, the negative
class disperses attention across various regions within the WSI despite all
presenting similar LG features. In regards to missclassified WSIs, we noted
that LG WSIs manifest both high attention and prediction scores, whereas HG
slides display a broad range of values. When examining the regression lines of
TP and FP, they exhibit similar trends, as do TN and FN respectively.
Essentially, misspredicted WSIs exhibit characteristics contrary to their
assigned class.

Figure 6: Correlation between region attention scores and output prediction of
region embeddings on the test set of WSIs. Regions from accurately predicted
WSIs (TN, TP) are denoted by squares, while those from incorrectly predicted
WSIs (FN, FP) are marked with crosses. A discernible pattern emerges, where
low attention scores align with diminished predictions, and conversely, higher
attention scores correlate with elevated predictions. Additionally, an
observable trend indicates that mispredictions tend to manifest on the
opposite end of the spectrum, with low-grade instances concentrating high
attention and prediction scores, and vice versa. The trend is represented with
a polynomial regression line (RL).

To augment the evaluation process, we integrate correlation calculations with
follow-up information, thereby ensuring a more thorough assessment of our
model’s performance. We employed Cramér’s V correlation coefficient
$\varphi_{c}$ for calculating the intercorrelation between grading and the
event of recurrence and progression, for the test set [49]. We observed a lack
of correlation between grading and recurrence for either the manual or
automatic grading, aligning with expectations. However, for progression,
NMGrad exhibited a higher correlation than the uropathologist (0.32 vs. 0.26,
respectively). In accordance with [50], a strong correlation between grading
and progression is observed. Notably, these correlations suggest that NMGrad’s
grade may hold greater predictive value for assessing the likelihood of
progression in the context of NMIBC.

##  7 Conclusion

Accurate grading of NMIBC is paramount for patient risk stratification, but it
has long suffered from inconsistencies and variations among pathologists.
Furthermore, the pathological workload is increasing, as well as its expenses.
In response to this challenge, we introduced the NMGrad pipeline, a pioneering
approach in bladder cancer grading using WSIs. NMGrad starts by using a tissue
segmentation algorithm, finding areas of urothelium in the slides. Thereafter,
it leverages a nested AbMIL model architecture to precisely identify
diagnostically relevant regions within WSIs and collectively predict tumor
grade. Moreover, through a multiscale CNN model, NMGrad processes urothelium
tissue tiles at multiple magnification levels. We observed that in high-grade
patients, attention scores pinpoint specific ROIs, while in low-grade
patients, attention is more dispersed, deviating from the expected MIL
pattern. Our clinical evaluations demonstrate that NMGrad consistently
outperforms previous state-of-the-art methods, achieving 0.94 AUC score. This
achievement represents a significant advancement in the field of bladder
cancer diagnosis, with the potential to improve patient outcomes, reduce
economic burdens, and enhance the quality of care in the management of this
challenging disease.

## Declarations

Ethical Approval This work involved human subjects or animals in its research.
Approval of all ethical and experimental procedures and protocols was granted
by the Regional Committees for Medical and Health Research Ethics (REC) under
Application No. 2011/1539, and performed in line with the Norwegian Health
Research Act.

Funding This research has received funding from the European Union’s Horizon
2020 research and innovation program under grant agreements 860627 (CLARIFY).

Conflict of Interest The authors declare no competing interests.

## References

  *   * [1] Teoh, J. Y.-C. _et al._ (2020).  Global trends of bladder cancer incidence and mortality, and their associations with tobacco use and gross domestic product per capita.  _European urology_ 78:893–906.
  * [2] Burger, M. _et al._ (2013).  Epidemiology and risk factors of urothelial bladder cancer.  _European urology_ 63:234–241.
  * [3] Babjuk, M. _et al._ (2022).  European association of urology guidelines on non–muscle-invasive bladder cancer (ta, t1, and carcinoma in situ).  _European urology_ 81:75–94.
  * [4] Eble, J. (2004).  World health organization classification of tumours.  _Pathology and genetics of tumours of the urinary system and male genital organs_ 68–69.
  * [5] Tosoni, I. _et al._ (2000).  Clinical significance of interobserver differences in the staging and grading of superficial bladder cancer.  _BJU International_ 85\.
  * [6] Netto, G. J. _et al._ (2022).  The 2022 world health organization classification of tumors of the urinary system and male genital organs—part b: prostate and urinary tract tumors.  _European urology_ .
  * [7] Hentschel, A. E. _et al._ (2020), Papillary urothelial neoplasm of low malignant potential (pun-lmp): Still a meaningful histo-pathological grade category for ta, noninvasive bladder tumors in 2019?  _Urologic Oncology: Seminars and Original Investigations_ , Vol. 38, 440–448.
  * [8] Jones, T. D. & Cheng, L. (2020).  Reappraisal of the papillary urothelial neoplasm of low malignant potential (punlmp).  _Histopathology_ 77:525–535.
  * [9] Van der Laak, J., Litjens, G. & Ciompi, F. (2021).  Deep learning in histopathology: the path to the clinic.  _Nature medicine_ 27:775–784.
  * [10] Komura, D. & Ishikawa, S. (2018).  Machine learning methods for histopathological image analysis.  _Computational and structural biotechnology journal_ 16:34–42.
  * [11] Cui, M. & Zhang, D. Y. (2021).  Artificial intelligence and computational pathology.  _Laboratory Investigation_ 101:412–422.
  * [12] Azam, A. S. _et al._ (2020).  Diagnostic concordance and discordance in digital pathology: a systematic review and meta-analysis.  _Journal of Clinical Pathology_ .
  * [13] Kanwal, N., Pérez-Bueno, F., Schmidt, A., Engan, K. & Molina, R. (2022).  The devil is in the details: Whole slide image acquisition and processing for artifacts detection, color variation, and data augmentation: A review.  _IEEE Access_ 10:58821–58844.
  * [14] Kvikstad, V. _et al._ (2019).  Prognostic value and reproducibility of different microscopic characteristics in the who grading systems for pta and pt1 urinary bladder urothelial carcinomas.  _Diagnostic pathology_ 14:1–8.
  * [15] van der Kwast, T. _et al._ (2022).  International society of urological pathology expert opinion on grading of urothelial carcinoma.  _European Urology Focus_ 8:438–446.
  * [16] Berbís, M. A. _et al._ (2023).  Computational pathology in 2030: a delphi study forecasting the role of ai in pathology within the next decade.  _EBioMedicine_ 88\.
  * [17] Ciompi, F. _et al._ (2017), The importance of stain normalization in colorectal tissue classification with convolutional networks.  _2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI)_ , 160–163.
  * [18] Fuster, S. _et al._ (2022), Invasive cancerous area detection in non-muscle invasive bladder cancer whole slide images.  _2022 IEEE 14th Image, Video, and Multidimensional Signal Processing Workshop (IVMSP)_ , 1–5.
  * [19] Cruz-Roa, A. _et al._ (2018).  High-throughput adaptive sampling for whole-slide histopathology image analysis (hashi) via convolutional neural networks: Application to invasive breast cancer detection.  _PloS one_ 13:e0196828.
  * [20] Litjens, G. _et al._ (2016).  Deep learning as a tool for increased accuracy and efficiency of histopathological diagnosis.  _Scientific reports_ 6:26286.
  * [21] Reis, S. _et al._ (2017).  Automated classification of breast cancer stroma maturity from histological images.  _IEEE Transactions on Biomedical Engineering_ 64:2344–2352.
  * [22] Dundar, M. M. _et al._ (2011).  Computerized classification of intraductal breast lesions using histopathological images.  _IEEE Transactions on Biomedical Engineering_ 58:1977–1984.
  * [23] Andreassen, C., Fuster, S., Hardardottir, H., Janssen, E. A. & Engan, K. (2023), Deep learning for predicting metastasis on melanoma wsis.  _2023 IEEE 20th International Symposium on Biomedical Imaging (ISBI)_ , 1–5.
  * [24] Tabatabaei, Z., Colomer, A., Moll, J. O. & Naranjo, V. (2023).  Toward more transparent and accurate cancer diagnosis with an unsupervised cae approach.  _IEEE Access_ 11:143387–143401.
  * [25] Deng, R. _et al._ (2023).  Omni-seg: A scale-aware dynamic network for renal pathological image segmentation.  _IEEE Transactions on Biomedical Engineering_ 70:2636–2644.
  * [26] Dietterich, T. G., Lathrop, R. H. & Lozano-Pérez, T. (1997).  Solving the multiple instance problem with axis-parallel rectangles.  _Artificial intelligence_ 89:31–71.
  * [27] Maron, O. & Lozano-Pérez, T. (1997).  A framework for multiple-instance learning.  _Advances in neural information processing systems_ 10:570–576.
  * [28] Carbonneau, M.-A., Cheplygina, V., Granger, E. & Gagnon, G. (2018).  Multiple instance learning: A survey of problem characteristics and applications.  _Pattern Recognition_ 77:329–353.
  * [29] Mercan, E. _et al._ (2016).  Localization of diagnostically relevant regions of interest in whole slide images: a comparative study.  _Journal of digital imaging_ 29:496–506.
  * [30] Tibo, A., Frasconi, P. & Jaeger, M. (2017), A network architecture for multi-multi-instance learning.  _Machine Learning and Knowledge Discovery in Databases: European Conference (ECML PKDD)_ , 737–752.
  * [31] Tibo, A., Jaeger, M. & Frasconi, P. (2020).  Learning and interpreting multi-multi-instance learning networks.  _The Journal of Machine Learning Research_ 21:7890–7949.
  * [32] Fuster, S., Eftestøl, T. & Engan, K. (2022), Nested multiple instance learning with attention mechanisms.  _2022 21st IEEE International Conference on Machine Learning and Applications (ICMLA)_ , 220–225.
  * [33] Khoraminia, F. _et al._ (2023).  Artificial intelligence in digital pathology for bladder cancer: Hype or hope? a systematic review.  _Cancers_ 15:4518.
  * [34] Wetteland, R. _et al._ (2021).  Automatic diagnostic tool for predicting cancer grade in bladder cancer patients using deep learning.  _IEEE Access_ 9:115813–115825.
  * [35] Zheng, Q. _et al._ (2022).  Accurate diagnosis and survival prediction of bladder cancer using deep learning on histological slides.  _Cancers_ 14:5807.
  * [36] Jansen, I. _et al._ (2020).  Automated detection and grading of non–muscle-invasive urothelial cell carcinoma of the bladder.  _The American journal of pathology_ 190:1483–1490.
  * [37] Spyridonos, P. _et al._ (2006).  Comparative evaluation of support vector machines and probabilistic neural networks in superficial bladder cancer classification.  _Journal of Computational Methods in Sciences and Engineering_ 6:283–292.
  * [38] Zhang, Z., Xie, Y., Xing, F., McGough, M. & Yang, L. (2017), Mdnet: A semantically and visually interpretable medical image diagnosis network.  _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 6428–6436.
  * [39] Zhang, Z. _et al._ (2019).  Pathologist-level interpretable whole-slide cancer diagnosis with deep learning.  _Nature Machine Intelligence_ 1:236–245.
  * [40] Ilse, M., Tomczak, J. & Welling, M. (2018), Attention-based deep multiple instance learning.  _International conference on machine learning_ , 2127–2136.
  * [41] Campanella, G. _et al._ (2019).  Clinical-grade computational pathology using weakly supervised deep learning on whole slide images.  _Nature medicine_ 25:1301–1309.
  * [42] Wetteland, R., Engan, K., Eftestøl, T., Kvikstad, V. & Janssen, E. A. (2020).  A multiscale approach for whole-slide image segmentation of five tissue classes in urothelial carcinoma slides.  _Technology in Cancer Research & Treatment_ 19:1533033820946787.
  * [43] Dalheim, O. N., Wetteland, R., Kvikstad, V., Janssen, E. A. M. & Engan, K. (2020), Semi-supervised tissue segmentation of histological images.  _Colour and Visual Computing Symposium_ , Vol. 2688, 1–15.
  * [44] Fuster, S., Khoraminia, F., Eftestøl, T., Zuiverloon, T. C. & Engan, K. (2023), Active learning based domain adaptation for tissue segmentation of histopathological images.  _2023 31st European Signal Processing Conference (EUSIPCO)_ , 1045–1049.
  * [45] Kvikstad, V. (2022).  _Better prognostic markers for nonmuscle invasive papillary urothelial carcinomas_.  Ph.D. thesis, Universitetet i Stavanger.
  * [46] Wetteland, R., Engan, K. & Eftestøl, T. (2021), Parameterized extraction of tiles in multilevel gigapixel images.  _2021 12th International Symposium on Image and Signal Processing and Analysis (ISPA)_ , 78–83.
  * [47] Krizhevsky, A., Sutskever, I. & Hinton, G. E. (2012).  Imagenet classification with deep convolutional neural networks.  _Advances in neural information processing systems_ 25\.
  * [48] Abraham, N. & Khan, N. M. (2019), A novel focal tversky loss function with improved attention u-net for lesion segmentation.  _2019 IEEE 16th international symposium on biomedical imaging (ISBI)_ , 683–687.
  * [49] Bergsma, W. (2013).  A bias-correction for cramér’s v and tschuprow’s t.  _Journal of the Korean Statistical Society_ 42:323–328.
  * [50] Akoglu, H. (2018).  User’s guide to correlation coefficients.  _Turkish journal of emergency medicine_ 18:91–93.

.

Saul Fuster is currently pursuing a Ph.D.in Science and Technology at the
Department of Electrical and Computer Engineering, University of Stavanger.
His main research interest is on the application of deep learning algorithms
for medical image analysis.

Umay Kiraz is a Ph.D. candidate at the Department of Pathology at Stavanger
University Hospital, Norway. Her field of study is improving diagnosis and
prognosis in breast cancer with artificial intelligence, including molecular
and digital pathology.

Trygve Eftestøl is currently a Professor at the Department of Electrical and
Computer Engineering, University of Stavanger (UiS). He is deputy leader of
the Biomedical Data Analysis Laboratory and in the steering comittee of the
Cognitive and Behavioral Neuroscience Lab at UiS. His main research interests
include biomedical signal analysis and applied machine learning.

Emiel A.M. Janssen is currently director of the Department of Pathology and
head of the Pathology Research Group at Stavanger University Hospital, Norway.
He is also a professor in Biomedicine at the University of Stavanger, Norway,
and adjunct professor at Griffith University, Australia.

Kjersti Engan is currently a Professor at the Department of Electrical and
Computer Engineering, University of Stavanger (UiS). She is the leader of the
Biomedical Data Analysis Laboratory. Her research interests include signal
processing and machine learning, emphasizing on medical applications.


# Optimal land conservation decisions for multiple species

Cassidy K. Buhler  Department of Decision Sciences and MIS, Bennett S. LeBow
College of Business, Drexel University, Philadelphia, PA 19104
<EMAIL_ADDRESS>and  Hande Y. Benson  Department of Decision Sciences and
MIS, Bennett S. LeBow College of Business, Drexel University, Philadelphia, PA
19104<EMAIL_ADDRESS>
###### Abstract.

Given an allotment of land divided into parcels, government decision-makers,
private developers, and conservation biologists can collaborate to select
which parcels to protect, in order to accomplish sustainable ecological goals
with various constraints. In this paper, we propose a mixed-integer
optimization model that considers the presence of multiple species on these
parcels, subject to predator-prey relationships and crowding effects.

##  1\. Introduction

###  1.1. Motivation

Climate change poses a massive threat to the health of humanity [8]. The
ramifications of the climate crisis extend beyond increased global
temperatures, as necessities like food and water will become more and more
scarce. Humans are not the only species that will be affected. At this rate,
it’s been estimated that one in six species could face extinction [10]. In an
interconnected world, it is, therefore, important to incorporate
sustainability with a focus on biodiversity into every level of decision
making.

Government decision-makers at the municipal, county, state, and federal levels
frequently work with private companies, including engineering firms and
developers, to determine where development should take place and where
“protected areas” must be established. That is, spaces where species are
protected from human interference. In this paper, we will be focusing on the
design of protected areas.

###  1.2. Designing protected areas

Given an allotment of land divided into parcels, our task is to select which
parcels to protect, in order to accomplish ecological goals subject to various
financial and economic constraints. These decisions are aided by optimization
models. [1, 2]

Protected areas can take on many configurations. Recent work in this area is
developing models to ensure particular spatial properties. For instance, some
researchers have focused on selecting parcels such that the protected area is
connected [3, 6, 4] or contiguous [11, 7]. However, many models that include
these additional constraints focus on a single species; incorporating multiple
species adds more dimensions to an already computationally expensive problem
to solve.

##  2\. Model and Methodologies

In this section, we present the optimization models for making land
preservation decisions in the presence of multiple species and budget
constraints. We start by introducing our notation, then present a model where
the species do not interact and one where predator-prey relationships are
present.

###  2.1. Notation

The notation used in the optimization models and in the rest of the paper are
as follows:

Sets:

$P$ | the set of parcels
---|---
$S$ | the set of species

Parameters:

$w_{i}$ |  Weight to prioritize species $i$
---|---
$N_{i}(p)$ |  Number of individuals of species $i$ observed at parcel $p$
$N_{i}$ |  Total population of species $i$ across all parcels
$\tilde{N}_{i}(p)$ |  Number of individuals of species $i$ that are simulated to be at parcel $p$ in the future
$c_{p}$ |  The cost of preserving parcel $p$
$B$ |  Budget

Decision Variables:

$x_{p}$ | Binary variable denoting whether or not a parcel is preserved
---|---

###  2.2. Model without Species Interaction

Our baseline model assumes that land preservation decisions can be made by
taking into account only the number of individuals of each species as
currently observed. This means we assume the different species do not interact
with each other and that crowding effects do not occur. Another interpretation
is that the model takes only the present conditions into consideration,
instead of focusing on sustainability.

This baseline model (1) is a knapsack problem, where the objective is to save
a weighted combination of species’ populations subject to a budget constraint.
The weights, $w_{i}$, can be chosen to prioritize endangered species or
reflect other conservation concerns.

(1) $\begin{array}[]{ll}\text{maximize }&\displaystyle\sum_{i\in
S}w_{i}\sum_{p\in P}N_{i}(p)x_{p}\\\ \text{subject to
}&\displaystyle\sum_{p\in P}c_{p}x_{p}\leq B\\\ &x_{p}\in\\{0,1\\}\quad\forall
p\in P\end{array}$

###  2.3. Model with Species Interaction

The main model we present is the model with species interaction, such as
predator-prey relationships. The model itself has the same overall form as (1)
with one critical difference. The species’ populations $\tilde{N}_{i}(p)$ are
calculated using a simulation that models the numbers of individuals present
in each parcel after $T$ time periods. In our numerical testing, we have taken
$T=2000$, which represents steady-state populations.

(2) $\begin{array}[]{ll}\text{maximize
}&\displaystyle\sum_{i}w_{i}\sum_{p}\tilde{N}_{i}(p)x_{p}\\\ \text{subject to
}&\displaystyle\sum_{p}c_{p}x_{p}\leq B\\\ &x_{p}\in\\{0,1\\}\quad\forall p\in
P\end{array}$

$\tilde{N}_{i}(p)$ is obtained from $N_{i}(p)$ using the Gavina et.al.’s model
from [5], which adapts the classical Lotka-Volterra equations describing
predator-prey relationships to multiple species and takes into account
crowding effects. The simulation is described in Algorithm 1.

input : $N_{i}(p)$: Number of individuals of species $i$ observed at parcel
$p$ for each $i$ and $p$

output : $\tilde{N}_{i}(p)$: Number of individuals of species $i$ that are
simulated to be at parcel $p$ in the future for each $i$ and $p$

foreach  _parcel $p$_ do

Initialize the species counts in parcel $p$, i.e. let
$\tilde{N}_{i}(p)=N_{i}(p)$ for each $i$;

for  _$t\leftarrow 1$ to $T$_ do

Let $\Delta N=$ (Birth Rate of species $i$)($\tilde{N}_{i}(p)$) - (Competition
Effect of all species on $i$) - (Crowding Effect for species $i$);

Let $\tilde{N}_{i}(p)=\tilde{N}_{i}(p)+\Delta N$.

Algorithm 1 Lotka-Volterra competition with crowding effects

##  3\. Data

The data was generated using the code provided in Gupta et.al.’s paper [6].
The process is outlined below.

A landscape is an $n\times n$ grid of parcels, where each parcel is a piece of
land that can be protected. Each parcel has a value between $[0,1]$ where a
higher value represents a worse habitat for that species. A landscape where
better habitat parcels are clustered with each other (and the worse habitat
parcels are clustered with each other) has a low habitat fragmentation. A high
habitat fragmentation is an absence of these clusters, as the better habitat
parcels are dispersed among worse habitat parcels, and vice versa. Based on
these landscapes, density of each species are simulated. Given a species
population $N_{i}$, these individuals are distributed among the landscape
using an inhomogeneous point process. This distribution is $N_{i}(p)$.

Following the model given in [5], we obtained $\tilde{N}_{i}(p)$ by inputting
$N_{i}(p)$ into Algorithm 1. Due to the scaling in this code, the values
outputted are non-integers. To address this, we rounded the output,
$\tilde{N}_{i}(p)$, to the nearest integer.

Habitat preferences were addressed by assigning landscapes and populations to
species, which yields unique distributions of each species on a grid. We
generated 10,000 $10\times 10$ landscapes with random habitat fragmentation
levels and the two most and two least fragmented landscapes were selected. For
each of these four landscapes, we explored two population sizes, $N_{i}=100$
and $N_{i}=250$. This yields 8 different distributions that represent 8
species. Additional details can be found in Table 1.

##  4\. Numerical Testing

With the 8 species generated, we grouped them into sets of size $2$ and $5$ in
order to implement a 2-species reserve and a 5-species reserve. For the
2-species reserves, the sets are
$\\{S_{0},S_{1}\\},\\{S_{2},S_{3}\\},\\{S_{4},S_{5}\\},$ and
$\\{S_{6},S_{7}\\}$. The 5-species reserve, the sets are
$\\{S_{0},S_{1},S_{2},S_{3},S_{4}\\}$ and
$\\{S_{5},S_{6},S_{7},S_{0},S_{1}\\}$. These give 6 scenarios total: 4 for
2-species reserves and 2 for 5-species reserves.

For each scenario, we varied the budget $B$ from 0 to 100 with a stepsize of 5
and solved (1) with $N_{i}(p)$ and (2) with $\tilde{N}_{i}(p)$. We tracked the
similarity of each model’s solution by counting the number of parcels that had
the same protection status, and taking the minimum, mean, and median over all
budgets except where budget is 0 and 100. We omitted those values because the
solution will always be to protect nothing or protect everything, which is
uninteresting for comparison.

The numerical results for all cases can be found in Table 2. In addition, one
solution for a 2-species reserve and 5-species reserve are included in Figure
3 and Figure 2. The reserves found using the two models show a high degree of
similarity when the objective function weights for the two species are the
same. However, when these weights differed, the solutions showed more
variation among lower budgets, which is illustrated by Figure 4. A comparison
between Figures 3 and 1 show a specific example of the configuration change.

##  5\. Discussion

A drawback to Algorithm 1, is that $\tilde{N}$ needs to be rounded. It is
worth exploring alternatives or variations in order to obtain integer
solutions. Not only that, but it would be interesting to modify the parameters
in Algorithm 1 to explore dynamics that would yield a larger difference
between $N$ and $\tilde{N}$.

For future directions, we hope to explore extensions to make the model and
numerical testing more realistic. This includes increasing the grid size,
expanding methods to obtain $N$, and investigating different parameters to use
for the species weights and parcel costs. A larger grid size would be valuable
to pursue because real-world landscapes are typically larger than $10\times
10$. Also, this would allow more possible reserve solutions, thus making the
problem more interesting. With regards to the methods to obtain $N$,
estimating the location and movement of a species using the spatial capture-
recapture (SCR) model [9] and using these as inputs into our model, as done in
Gupta et.al’s paper, would provide a more accurate depiction of animal
behavior.

## References

  * [1] Diogo Alagador and Jorge Orestes Cerdeira.  Operations research applicability in spatial conservation planning.  Journal of Environmental Management, 315:115172, 2022.
  * [2] Alain Billionnet.  Mathematical optimization ideas for biodiversity conservation.  European Journal of Operational Research, 231(3):514–534, 2013\.
  * [3] Jennifer K Costanza, James Watling, Ron Sutherland, Curtis Belyea, Bistra Dilkina, Heather Cayton, David Bucklin, Stephanie S Romañach, and Nick M Haddad.  Preserving connectivity under climate and land-use change: No one-size-fits-all approach for focal species in similar habitats.  Biological Conservation, 248:108678, 2020.
  * [4] Bistra Dilkina, Rachel Houtman, Carla P Gomes, Claire A Montgomery, Kevin S McKelvey, Katherine Kendall, Tabitha A Graves, Richard Bernstein, and Michael K Schwartz.  Trade-offs and efficiencies in optimal budget-constrained multispecies corridor networks.  Conservation Biology, 31(1):192–202, 2017.
  * [5] Maica Krizna A Gavina, Takeru Tahara, Kei-ichi Tainaka, Hiromu Ito, Satoru Morita, Genki Ichinose, Takuya Okabe, Tatsuya Togashi, Takashi Nagatani, and Jin Yoshimura.  Multi-species coexistence in lotka-volterra competitive systems with crowding effects.  Scientific reports, 8(1):1–8, 2018.
  * [6] Amrita Gupta, Bistra Dilkina, Dana J Morin, Angela K Fuller, J Andrew Royle, Christopher Sutherland, and Carla P Gomes.  Reserve design to optimize functional connectivity and animal density.  Conservation Biology, 33(5):1023–1034, 2019.
  * [7] Hayri Önal, Yicheng Wang, Sahan TM Dissanayake, and James D Westervelt.  Optimal design of compact and functionally contiguous conservation management areas.  European Journal of Operational Research, 251(3):957–968, 2016\.
  * [8] Annette Prüss-Üstün, Jennyfer Wolf, Carlos Corvalán, Robert Bos, and Maria Neira.  Preventing disease through healthy environments: a global assessment of the burden of disease from environmental risks.  World Health Organization, 2016.
  * [9] Chris Sutherland, Angela K Fuller, and J Andrew Royle.  Modelling non-euclidean movement and landscape connectivity in highly structured ecological networks.  Methods in Ecology and Evolution, 6(2):169–177, 2015.
  * [10] Mark C Urban.  Accelerating extinction risk from climate change.  Science, 348(6234):571–573, 2015.
  * [11] Yicheng Wang, Qiaoling Fang, Sahan TM Dissanayake, and Hayri Önal.  Optimizing conservation planning for multiple cohabiting species.  PloS one, 15(6):e0234968, 2020.

##  6\. Appendix

Species | Fragmentation level | $N_{i}$
---|---|---
$S_{0}$ | highest | 100
$S_{1}$ | 2nd highest | 100
$S_{2}$ | highest | 250
$S_{3}$ | 2nd highest | 250
$S_{4}$ | lowest | 100
$S_{5}$ | 2nd lowest | 100
$S_{6}$ | lowest | 250
$S_{7}$ | 2nd lowest | 250
Table 1. Details for each species’ distribution. Fragmentation level describes the density of habitat quality on a landscape, and $N_{i}$ is the total population of a species on that landscape.  Case | Min | Average | Median
---|---|---|---
1 | 98 | 99.79 | 100
2 | 92 | 97.26 | 98
3 | 92 | 95.89 | 96
4 | 92 | 95.89 | 96
5 | 94 | 96.95 | 96
6 | 94 | 97.16 | 98
Table 2. These values represent the similarity of solutions given from (1) and
(2) over budgets varying from 5 to 95 with stepsizes of 5. For the 100
parcels, the number of parcels with the same protection status is recorded.
The minimum, average, and median of this is computed and displayed. Case 1:
$\\{S_{0},S_{1}\\}$, Case 2: $\\{S_{2},S_{3}\\}$, Case 3: $\\{S_{4},S_{5}\\}$,
Case 4: $\\{S_{6},S_{7}\\}$, Case 5: $\\{S_{0},S_{1},S_{2},S_{3},S_{4}\\}$,
and Case 6: $\\{S_{5},S_{6},S_{7},S_{0},S_{1}\\}$ were all conducted.  Figure
1. This graph depicts the same scenario as Figure 3, but rather than setting
$w_{2}=w_{3}$, we used $w_{2}=0.9$ and $w_{3}=0.1$. In this example, $90/100$
parcels have the same protection status. Figure 2. Displaying the solution to
(1) and (2) for a 5-species reserve using
$\\{S_{0},S_{1},S_{2},S_{3},S_{4}\\}$ when the budget is 55. Each parcel is
annotated with its corresponding $N_{i}(p)$ or $\tilde{N}_{i}(p)$ for
$i=0,\ldots,4$. The order of the values in each parcel are top left, top
right, middle, bottom left, bottom right, for species
$S_{0},S_{1},S_{2},S_{3},S_{4}$, respectively. The color of each parcel
represents the preservation decision: green if preserved, orange if not. In
this example, $96/100$ parcels have the same protection status. Figure 3.
Displaying solutions to (1) and (2) for a 2-species reserve using
$\\{S_{2},S_{3}\\}$ and $B=55$. Each parcel is annotated with its
corresponding $N_{i}(p)$ or $\tilde{N}_{i}(p)$ for $i=2,3$. For each parcel,
the value for $S_{2}$ is on top and $S_{3}$ is on bottom. The color of each
parcel represents the preservation decision: green if preserved, orange if
not. In this example, $92/100$ parcels have the same protection status. Figure
4. The similarity of solutions from (1) and (2) vary based on the objective
function weights, $w_{i}$. In this example, we used $S_{2}$ and $S_{3}$ for a
2-species reserve, and solved the models where $w_{2}=w_{3}$ (unweighted) and
where $w_{2}=0.9$ and $w_{3}=0.1$ (weighted). For lower budgets, the
difference is more apparent.