Similarly to how learning happens in humans, multi-task architectures aim at simultaneously learning multiple tasks that are related to each other. In principle, learning two related tasks generates mutually beneficial signals that lead to more general and robust representations than traditional or multimodal learning. Figure 1 a) and b) better illustrate the original concept proposed by Caruana (1997). Suppose that a complex model, e.g. a CNN, is trained on the main task M. In Figures 1 (a) and (c) the optimization objective of M has two local minima, represented as the set {a, b}. in Figure 1 (a), the auxiliary task A is related to the main task, with which it shares the local minimum in a. The joint optimization of M and A is thus likely to identify the shared local minima a as the optimal solution (Caruana, 1997). The search is biased by the extra information given by task A towards representations that lay at the intersection of what could be learned individually for each task. In Figure 1 (b), the auxiliary task is totally unrelated to the main task. No local minima are shared in this case and a negative transfer may happen without positive improvements to the performance.

Multi-task architectures divide into two families depending on the hard or soft sharing of the parameters. In architectures with hard parameter sharing such as the one proposed in this paper, multiple supervised tasks share the same input and some intermediate representation. The parameters learned up to this intermediate point are called generic parameters since they are shared across all tasks. In soft parameter sharing, the weight updates are not shared among the tasks and the parameters are task-specific, introducing only a soft constraint on the training process (Duong et al., 2015).

As explained by Caruana (1997), multi-task learning leads to various benefits if the tasks are linked by a valid relationship. For instance, the generalization error bounds are improved and the risk of overfitting is reduced(Baxter, 1995). The speed of convergence is also increased since fewer training samples are required per task (Baxter, 2000). Because of this, multi-task learning has been successful in various applications, such as natural language processing (Subramanian et al., 2018), computer vision (Kokkinos, 2017), autonomous driving (Leang et al., 2020), radiology (Andrearczyk et al., 2021) and histology (Gamper et al., 2020b; Marini et al., 2022). The preliminary work by Gamper et al. (2020b), in particular, shows a decrease in the loss variance as an effect of multi-task for oral cancer, suggesting that this work may have a high potential for histology applications. Depending on the applications and on the loss functions used to represent the multiple tasks, multiple strategies exist for weighting each task contribution in the objective function Gong et al. (2019). Alternatively to uniform weighting, for example, dynamical task re-weighting during training was proposed by Leang et al. (2020) and Kendall et al. (2018). The latter, in particular, exploits uncertainty estimates to weight each task, making the model convergence more robust to distribution shifts and unseen samples.

The existing frameworks do not consider yet the combination of multitask learning with adversarial tasks. In Figure 1 (c), we illustrate how our contribution in this paper differs from existing studies. For instance, we extend the multi-task framework to introduce an adversarial task C that is learned by adversarial training Ganin et al. (2016). In this case, the main task M has local minima in {a, b, c}, but the minimum in c is also a solution to the adversarial task C. We hypothesize that, by being adversarial to C, the optimization will be likely to prefer solutions that satisfy M and A, while avoiding solutions that satisfy the adversarial task C. This is because the main and auxiliary tasks (M and A) work in a cooperative setting with the feature encoder to learn a representation that will minimize both of their losses (Caruana, 1997). On the other hand, the adversarial task (i.e. C) plays an adversarial game with the feature encoder as it is described by Ganin et al. The encoder will optimize the representations so as to maximize the loss on the adversarial task, reducing the possibilities of recovering the information about the undesired target no matter how hard the adversarial branch tries. The architecture in this work implements the extension of multitask learning to include adversarial tasks, and this is one of the main contributions of our work.

Proposed by Ganin et al. (2016), adversarial learning introduced a novel approach to solving the so-called problem of domain adaptation, namely the minimization of the domain shift in the distributions of the training (also called source distribution) and testing data (i.e. target). Typically treated as either an instance re-weighting operation (Gong et al., 2013) or as an alignment problem (Long et al., 2015), domain adaptation is handled by adversarial learning as the optimization of a domain confusion loss. A domain classifier discriminates between the source and the target domains during training and its parameters are optimized to minimize the error when discriminating the domain labels. This can be extended to more than two domains by a multi-class domain classifier. The adversarial learning of domain-related features is obtained by a gradient reversal operation on the branch learning to discriminate the domains. Because of this operation, the network parameters are optimized to maximize the loss of the domain classifier, thus making multiple domains impossible to distinguish one from another in the internal network representation. This causes competition between the main task and the domain branch during training which is referred to as a min-max optimization framework. As a downside, the optimization of adversarial losses may be complicated, with the min-max operation affecting the stability of the training (Ganin et al., 2016). Convergence can be promoted, however, by activating and de-activating the gradient reversal branch according to a training schedule as in Lafarge et al. (2017).

Adversarial learning is one building block of our architecture, since we incorporate the adversarial task as an additional branch of our multi-task architecture. Differently from the main work by Ganin et al. (2016), we introduce an additional hyper-parameter that controls the activation of the gradient reversal operator, so that the gradient flow is reversed for the adversarial tasks and kept unchanged for the rest of the auxiliary tasks in the architecture.

One important feature of our architecture is that it introduces the learning of interpretable, high-level features as additional tasks to regularize the training process. This introduces transparency to the architecture, since the additional tasks introduce interpretable inductive biases during training. We revise, in the following, the concept of interpretability of deep learning, clarifying how this constitutes an important building block for our contributions. Interpretability of deep learning has become increasingly important over the past decade, which saw the definition of a plethora of methods and approaches with discording terminologies (Graziani et al., 2022). In the context of digital pathology, interpretability analyses mainly focus on achieving post-hoc explanations rather than direct interpretable modelling (H.M. et al., 2022). Linear models, in particular, were proposed as an inherently interpretable approach to probe the internal activations of the network after training. Regression Concept Vectors (RCVs) and Concept Activation Vectors (CAVs) were shown to generate insightful explanations in terms of diagnostic measures (Kim et al., 2018; Graziani et al., 2018, 2019b, 2019c; Yeche et al., 2019). RCVs solve a linear regression problem in the space of the internal activations of a CNN layer. They were used to evaluate the relevance of a concept, e.g. a clinical feature such as nuclei area or tumour extension, in the model output generation. The performance of the linear regressor indicated how well the concept was learned. Both CAVs and RCVs constitute a baseline of linear interpretability of CNNs, formalized for applications in the medical domain as concept attribution (Graziani et al., 2020).

No possibility is given by the existing interpretability methods to act on the training process and modify the learning of a concept. It is not possible, for example, to discourage the learning of a confounding concept, e.g. domain, staining, watermarks. Similarly, the learning of discriminant concepts cannot be further encouraged. With this paper, we aim at addressing this gap.

The architecture for guiding the training of CNNs is described in this section for a general application with pre-defined features. This general framework can be applied to multiple tasks. The diagnosis of cancerous tissue in breast microscopy images is proposed in this paper as an application for which the implementation details are described in Sections 3.5 and 3.6.

In the following, we clarify the notation used to describe the model. We assume that a set of $N$ observations, i.e. the input images, is drawn from an unknown underlying distribution and split into a training subset $\{\bm{x}_{i}\}_{i=1}^{n}$ and a test subset $\{\bm{x}_{i}\}_{i=n+1}^{N}$. The main task, namely the one for which we aim at improving the generalization, is the prediction of the image labels $\bm{y}=\{y_{i}\}_{i=1}^{n}$, for which ground truth annotations are available. A CNN of arbitrary structure is used as a feature encoder, of which the features are then passed through a stack of dense layers. The model parameters up to this point are defined as $\bm{\theta}_{f}$. The parameters of the label prediction output layers are identified by $\bm{\theta}_{y}$. The structure described up to this point replicates a standard CNN with a single main task branch that is addressing the classification. The remaining parameters of the architecture implement (i) the learning of auxiliary tasks by multi-task learning (Caruana, 1997) and (ii) the adversarial learning of detrimental features to induce invariance in the representations, as in the domain adversarial approach by Ganin et al. (2016). We combine these two approaches by introducing $K$ extra targets representing desired and undesired tasks that must be introduced to the learning of the representations. The targets are modeled as the prediction of the feature values $\{c_{k,i}\}_{i=1}^{N}$, where $k\in{1,\dots,K}$ is an index representing the extra task being considered. The feature values may be either continuous or categorical. Additional parameters $\bm{\theta}_{k}$ are trained in parallel to $\bm{\theta_{y}}$ for the $K$ extra targets. We refer to all model outputs for all inputs $\bm{x}$ as $f(\bm{x})\in\mathcal{R}^{K+1}$ .

The architecture is illustrated in Figure 2 and consists of two blocks. The first block is used to extract features from the input images. A state-of-the-art CNN of arbitrary choice without the decision layer is used as a feature encoder generating a set of feature maps. The feature maps are passed through a Global Average Pooling (GAP) operation that is performed to spatially aggregate the responses and connect them to a stack of dense layers. For this specific architecture, we use a stack of three dense layers of 1024, 512, and 256 nodes respectively. The second block comprises one branch per task, taking as input the output of the first block. The main task branch consists of the prediction of the labels $\mathbf{y}$ and has as many dense nodes as there are of unique classes in $\mathbf{y}$. For binary classification tasks, e.g. discrimination of tumorous against non-tumorous inputs, the main task branch has a single node with a sigmoid activation function. $K$ branches are added to model the extra targets. We refer to extra tasks for all the additional targets to the main task whether desired or undesired. Auxiliary tasks refer to the modeling of the desired targets, while adversarial tasks refer to that of undesired targets. The extra tasks are modeled by linear models as in Graziani et al. (2018). For continuous-valued targets, the extra branch consists of a single node with a linear activation function. For categorical targets, the extra branch has multiple nodes followed by a softmax activation function. A gradient reversal operation (Ganin et al., 2016) is performed on the branches of the undesired targets to discourage the learning of these features.

The objective function of the proposed architecture balances the losses of the main task and the extra tasks for the desired and undesired targets. This is obtained by a combination of multi-task and adversarial learning. The main task loss is $\mathcal{L}_{y}^{i}(\bm{\theta_{f}},\bm{\theta_{y}})=\mathcal{L}_{y}(\bm{x}_{i},y_{i};\bm{\theta_{f}},\bm{\theta_{y}})$, where $\bm{\theta}_{f}$ are the parameters of the first block (namely of the CNN encoder and the dense layers) in Figure 2 and $\bm{\theta}_{y}$ those of the main task branch in the second block of the same figure. The extra parameters $\bm{\theta}_{k}$ ($k\in{1,\dots,K}$) are trained for the branches of the desired and undesired target predictions, with the loss being $\mathcal{L}_{k}^{i}(\bm{\theta}_{f},\bm{\theta}_{k})=\mathcal{L}_{k}(\bm{x}_{i},c_{k,i};\bm{\theta}_{f},\bm{\theta}_{k})$.

Training the model on $n$ training and $(N-n)$ testing samples consists of optimizing the function:

The gradient update is:

where $\lambda_{m}$ and $\lambda_{k}$ are positive scalar hyperparameters to tune the trade-off between the losses. For each extra branch, the hyperparameter $\alpha_{k}\in\{-1,1\}$ is used to specify whether the update is adversarial or not. A value of $\alpha_{k}=-1$ activates the gradient reversal operation and starts an adversarial competition between the feature extraction and the corresponding $k$^th extra branch. The main task is only trained on the training data, since $\mathcal{L}_{y}^{i}=0$ for $i>n$ in Eq. (2) and (3). Following the work in Ganin et al. (2016), the additional branches can be trained on slightly different dataset splits. If additional task labels are available for the test set, they can be included in the training of the additional branches. Moreover, if only partially labeled data are available, they can also be introduced in the training of the additional branches. For instance, the gradient updates can be kept only for the labeled portion of the data, setting the loss to zero for the unlabeled inputs.

The proposed architecture requires the combination of multiple objectives in the same loss function. The vanilla formulation in Eq. 1 simply performs a weighted linear sum of the losses for each task. This is the predominant approach used in prior work with multi-objective losses (Gong et al., 2019) and adversarial updates (Ganin et al., 2016; Lafarge et al., 2017). The appropriate choice of weighting of the different task losses is a major challenge of this setting. The tuning of the hyperparameters may reveal tedious and non-trivial due to the combination of classification and regression tasks with different ranges of the loss function values (e.g. combining the bounded binary cross-entropy loss in [0,1] with the unbounded mean squared error loss).

An optimal weighting approach may be learned simultaneously with the other tasks by adding network parameters for the loss weights $\lambda_{m}$ and $\lambda_{k}$. The direct learning of $\lambda_{m}$ and $\lambda_{k}$ would just result in weight values quickly converging to zero. Kendall et al. (2018) proposed a Bayesian approach that makes use of the homoscedastic uncertainty of each task to learn the optimal weighting combination. In loose words, homoscedastic uncertainty reflects a task-dependent confidence in the prediction. The main assumption to obtain an uncertainty-based loss weighting strategy is that the likelihood of the task output can be modeled as a Gaussian distribution with the mean given by the model output and a scalar observation noise $\sigma$:

This assumption is also applied to the outputs of the additional tasks. The loss weights $\lambda_{m}$ and $\lambda_{k}$ are then learned by optimizing the minimization objective given by the negative log likelihood of the joint probability of the task outputs given the model predictions. To clarify this concept, let us focus on a simplified architecture with the main task being the logistic regression of binary labels (e.g. tumor v.s. non-tumor) with noise $\sigma_{1}$ and one auxiliary task consisting of the linear regression of feature values $\bm{c}=\{c_{i}\}_{i=1}^{N}$, with noise $\sigma_{2}$. The minimization objective for this multi-task model is:

By minimizing Eq. 6 with respect to $\sigma_{1}$ and $\sigma_{2}$, the optimal weighting combination is learned adaptively based on the data (Kendall et al., 2018). As $\sigma_{1}$ increases, the weight for its corresponding loss decreases, and vice-versa. The last term $\log\sigma_{1}+\log\sigma_{2}$, besides, acts as a regularizer discouraging each noise to increase unreasonably. This construction can be extended easily to multiple regression outputs and the derivation for classification outputs is given in Kendall et al. (2018).

The experiments are run on three publicly available datasets, namely Camelyon 16, Camelyon 17 (Litjens et al., 2018a) and the breast subset of PanNuke (Gamper et al., 2019, 2020b) (The data are available for download at the links https://camelyon17.grand-challenge.org and https://warwick.ac.uk/fac/cross_fac/tia/data/pannuke. The Camelyon challenge collections of 2016 and 2017 contain respectively 270 and 899 WSIs. All training slides of both challenges contain annotations of metastasis type (i.e. negative, macro-metastases, micro-metastases, isolated tumor cells), and 320 images contain manual segmentations of tumor regions. The analysis also includes the breast tissue scans of the PanNuke dataset, for which multiple nuclei types were annotated by the semi-automatic instance segmentation tool described in Gamper et al. (2019). Labels of neoplastic, inflammatory, connective, epithelial and dead nuclei are given together with the images by the dataset creators. Training, validation and test splits are built on these two datasets as reported in Table 1. The pre-existing PanNuke folds are used, two of which (i.e. fold 1 and fold 2) are used in the training set. For external validation of our model we create the two splits described in Table 2, namely the CamExt and the PanExt splits. CamExt is built by leaving out from training the images of Camelyon17 that were acquired at one specific acquisition center, namely center 4 in the original collection. The lesions of patients in this center have a much larger incidence of macro lesions than the patients in the internal validation set, with 26 slides out of 100 reporting a tumor diameter larger than 2.0 mm. Because of the expanded tumor size, true positives may be relatively easy to obtain on this set since the patches are less likely to be sampled from the tumor borders, and thus do not contain normal and tumorous tissue at the same time. The PanExt dataset is further used to test model performance, and it comprises the data in the testing split of the PanNuke collection, namely the images in Fold 3 of the original dataset. We extract 775 image patches of which 480 contain tumor and 295 represent healthy tissue without tumorous cells. The target labels for the additional tasks were not computed for this dataset, hence these images were never seen by the model during training at the patch, slide and patient levels. All WSIs are pre-processed in the same way. Patches of $224\times 224$ pixels are extracted at the highest magnification level and the staining variability is reduced by Reinhard normalization Reinhard et al. (2001); Otálora et al. (2022). The PanNuke images do not depict an entire WSI but only a small portion, from which we also extract image patches of $224\times 224$ pixels. During training we perform oversampling with slightly overlapping patches by extracting patches located in the center, upper left, upper right, bottom left and bottom right corners of the image. This fully covers the WSI portion depicted in the original PanNuke image content and it is applied to balance the domain under-representation of these input images.

The main task that we address is the binary classification of input images that include tumor tissue from those without tumor. Inception V3 pretrained on ImageNet (Szegedy et al., 2016) is used as the backbone CNN for feature encoding. Following the observations in(Graziani et al., 2019a), the parameters up to the last convolutional layer are kept frozen to avoid overfitting to the pathology images¹¹1 A total of 7,804,427 trainable parameters are used in the first block for the feature extraction with InceptionV3 as the convolutional backbone and 6,427,392 with ResNet50. In the second block there are 257 to 2,827 trainable parameters, depending on the number of tasks being trained jointly. . The output of the CNN is passed through the GAP and the three fully connected layers as illustrated in Figure 2. The fully connected layers have respectively 2048, 512 and 256 units. A dropout probability of 0.8 and L2 regularization are added to these three fully connected layers to avoid overfitting. This configuration was obtained by optimizing performance only on the main task and by searching the hyperparameter space to identify the optimal dense block width, regularization strength and dropout probability. The main task is the detection of patches containing tumor as a binary classification task. The branch consists of a single node with sigmoid activation function connected to the output of the third dense layer. The architecture as described up to here, hence without extra branches, is used as the baseline for the experiments. The additional tasks consist of either the linear regression or the linear classification of continuous or categorical labels respectively. For linear regression, the additional branch is a single node with linear activation function. The Mean Squared Error (MSE) between the predicted value and the label is added to the optimization function in Eq. 1. For the linear classification, the extra branch has a number of dense nodes equal to the number of classes to predict and a softmax activation function, also connected to the third dense layer. The Categorical Cross-Entropy (CCE) loss is added to the optimization in Eq. 1. Further details about the extra branches used for the experiments are given in Section 3.6.

The architecture is trained end-to-end with mini-batch Stochastic Gradient Descent (SGD) with standard parameters (learning rate of $10^{-4}$ and Nesterov momentum of 0.9). The main task is learned by optimizing the class-Weighted Binary Cross Entropy (WBCE) loss, where positive samples are given a higher weight than the negative ones to counter-balance their under-representation in the training set. The weights are set so that they sum to one. The weight for the positive class corresponds to the number of samples missing to reach one, namely one minus the ratio of positive samples, i.e. 0.82. Similarly, the weight for the negative class is set to 0.18.

We evaluate the convergence of the network by early stopping on the total validation loss with patience of 5 epochs. The Area Under the ROC Curve (AUC) is used to evaluate model performance. For each experiment, we perform five runs with multiple initialization seeds to evaluate the performance variation due to initialization. The splits are kept unchanged for the multiple seed variations. To evaluate the performance on multiple test splits, we perform bootstrapping of the test sets. A number of 50 test sets of 7589 images (the total number of test images in the two sets) are obtained by sampling with replacement from the total pool of testing images. This method evaluates the variance of the test set without prior assumption on the data distribution and it shows the performance difference due to variation of the sampling of the population.

The experiments focus on the integration of four desired and one undesired targets with multiple combinations. The auxiliary targets relate to the main task, being important diagnostic features. We expect that learning of these desired features will improve the solution robustness and generalization of the model. Discarding the undesired targets may improve the invariance of the learned features to confounding factors. The Nottingham Histologic Grading (NHG) of breast tissue identifies the key diagnostic features for breast cancer (Bloom and Richardson, 1957). By analyzing this we derived the desired and undesired features that are illustrated in Figure 4. From this set, we retain cancer indicators at the nuclear level, since the input images are at the highest magnification. We model the variations of the nuclei size, appearance (e.g. irregular, heterogeneous texture) and density shown in Figure 4 as real-valued variables. Because of the heterogeneity of the data, we also guide the network training to discard information about the image acquisition center, which is modeled as an undesired target. The staining variability constitutes part of the shift, as differences are visible by the naked eye as shown in Figure 3. The staining variability, however, is only one of the components that contribute to creating the distributional shift among centers. Within the domain shift, other sources of variability further enhance the heterogeneity of the data, e.g. tissue fixation times, processing temperature and scanner resolution. These factors vary across institutions and jointly contribute to the distributional shifts observed among centers²²2For this reason, the staining normalization during preprocessing does not interfere with detecting the acquisition site..

Hand-crafted features representing the variations in the nuclei size and appearance are automatically extracted either from the images of from the nuclear contours. The nuclear contours are available in the form of manual annotations only for the PanNuke data. The manual delineation of nuclei contours may be cumbersome, costly and not entirely reflect the standard clinical routine procedure for cancer diagnosis. For this reason, we include nuclei segmentation labels that are obtained by an automatic segmentation model. For the images in Camelyon, automated contours of the nuclei are obtained by the multi-instance deep segmentation model in Otálora et al. (2020). This model is a Mask R-CNN model (He et al., 2017), fine-tuned from ImageNet weights on the Kumar dataset for the nuclei segmentation task (Kumar et al., 2017). The R-CNN identifies nuclei entities in the individual patches obtained from each WSI, and it then generates pixel-level masks by optimizing the Dice score. From the pre-existing labeled dataset of 30 annotated WSIs, we generate weak labels for over a thousands slides in the Camelyon datasets. ResNet50 (He et al., 2017) is used for the convolutional backbone as in (Otálora et al., 2020). The network is optimized by SGD with standard parameters (learning rate of 0.001 and momentum of 0.9).

Figure 5 shows the distribution of the feature values considered as additional targets for this study. Nuclei density in Figure 5(a) is estimated by counting the nuclei in each patch. The number of pixels inside nuclear contours is averaged for each input patch to represent variations of the nuclei area (b), referred to as area in the experiments.

Haralick descriptors of texture correlation and contrast (Haralick, 1979) are also extracted from the patches as in Graziani et al. (2018), for which the value distributions are shown in Figures 5(c) and 5(d) respectively. Being continuous and unbounded measures, the values for these features are normalized to have zero mean and unitary standard deviation before training the model. In the paper, we refer to these features as area, density, contrast and correlation. The values of these features are used as prediction labels for the auxiliary target branches, that are also named as the feature that they should predict. These auxiliary branches perform a linear regression task, trying to minimize the Mean Squared Error between the predicted value of the feature and the extracted values used as labels.

The center that performed the data acquisition may act as a confounder on both the main task and the additional, desired targets. Figure 6, for example, illustrates how the acquisition center impacts the nuclei density, i.e. the observed values of nuclei count. For this reason, we use information about the center that performed the data acquisition to model an adversarial task, using the center value that is present in the dataset as metadata. We model it as a categorical variable that may take values from 0 to 7, namely one for each known center in the data.

Since there is no specific information on acquisition centers in Camelyon16 and PanNuke, these have been modeled as two distinct acquisition centers in addition to the five known centers of Camelyon17. This information is partly inaccurate, since we know that in both datasets more than a single acquisition center was involved (Litjens et al., 2018b; Gamper et al., 2020a). The noise introduced by this information may limit the benefits introduced by the adversarial branch but it should not affect the performance negatively. In the future, unsupervised domain alignment methods may also be explored. The prediction of this variable is added to the architecture as an undesired target branch, referred to as center in the experiments.

We validate the predictive power of the features that we selected as auxiliary targets by a simple test. We train multi-layer perceptron models (MLP) with 20 hidden layers to predict the main task from the auxiliary label values that are, in this case, passed as input features and not as additional, multi-task targets. The models were trained with BCE loss and with the Adam optimizer with standard parameters such as learning rate of 0.01, and early stopping patience of 15. We evaluate the predictive AUC on the internal test set, which is reported in Table 3. As the results show, the main task predictions are better than random guessing for all the models that we trained, suggesting that the auxiliary labels contain information about the main task. Using multiple features at the same time (as in model ID MLP4), moreover, improves the performance of the classifier, showing that each auxiliary target contributes to the discrimination of the classes in the main task.

In the baseline model, only the main task branch is trained and no extra tasks are used. The results for the baseline models are shown in the first and second rows of Table 4, that are identified by unique IDs, namely model-ID 1 for the model with Inception V3 and model-ID R1 for the Res-Net model. The experiments are performed on the dataset configurations in Tables 1 and 2, where internal and external test sets are used to evaluate model generalization. Two columns are used to report the results on the internal (int.) and external (Cam.Ext.) test sets. Table 5 reports the results on PanExt. Where not stated otherwise, the average AUC (avg. AUC) over ten repetitions with multiple initialization seeds is used for the evaluation.

An ablation study is performed by adding a single additional task at a time to the baseline model. This study aims at identifying the benefits of encouraging each task individually. The desired and undesired targets described in Section 3.6 are added as additional branches to the architecture detailed in Section 3.5. The gradient reversal operation is only active for the center branch. The losses of each task are combined by two strategies, namely a vanilla and the uncertainty-based approach in Kendall et al. (2018). In the vanilla configuration, the loss weight values are set to 1 for all branches. The results for single task combinations are reported in Table 4 with unique IDs ranging from 2 to 5. Model-ID 2, for example, is given by the combination of the main task branch with the additional task area, namely of predicting the area of the nuclei in the images. A single auxiliary branch already outperforms the baseline (internal avg AUC $0.819\pm 0.001$, external avg. AUC $0,868\pm 0.005$, int. avg. F1 $0.783\pm 0.002$, ext. avg. F1 $0.315\pm 0.001$), as for example in model-ID 3 by encouraging nuclei count (internal avg AUC $0.836\pm 0.005$, external avg. AUC $0,890\pm 0.009$, internal avg. F1 $0.768\pm 0.001$, external avg. F1 $0.660\pm 0.003$) The avg. F1 score of model-ID3 on the two test sets is at $0.746\pm 0.009$, a significantly higher value than the $0.701\pm 0.001$ of the baseline model. On the external test set, the best generalization is achieved by adding count as a desired target (ext. avg. AUC $0.890\pm 0.009$).

The most promising branches are then combined to further improve performance. The following combinations of multiple auxiliary and adversarial tasks are tested in the experiments: center + density, center + area, center + density + area. The combination of all the branches leads to the best performance on the internal test for both the models with the Inception V3 and ResNet backbones, namely in model-ID 8 and R8. These models yield an increase from the baseline of 0.055 AUC points for model-ID 8 and 0.091 AUC points for R8. The highest performance on the internal test set is given by the R8 model with an AUC of $0.893\pm 0.001$ Model-ID 6 reports comparable performance on the external test set to model-ID 3. The addition of the center adversarial branch in model-ID 6 leads to the best Inception V3 based model overall with average AUC on both internal and external sets at $0.824\pm 0.006$ and avg. F1 score $0.755\pm 0.005$ for the uncertainty trained model. This represents a significant improvement compared to the overall average AUC $0.79\pm 0.001$ and avg. F1 score $0.701\pm 0.004$ of the baseline model, with $p-value<0.001$. The statistical significance of the results is evaluated by the non-parametric Wilcoxon test (two-sided) applied on the bootstrapping of the test set as described in Sec. 3.5. The performance on PanExt is reported for the baselines and the best performing models in Table 5.

To confirm the benefit of the added related tasks, we compare these results with those obtained with random noise as additional targets. This experiment is performed as a sanity check, where an auxiliary task is trained to predict random values. As expected, the overall, internal and external avg. AUCs are lower for this experiment and have larger standard deviations (overall avg. AUC $0.819\pm 0.04$, int. test AUC $0.834\pm 0.001$ and ext. avg. AUC $0.879\pm 0.03$). This shows that the selected tasks are more relevant to the main task than the regression of random values.

At this point, one may ask whether it is possible to interpret the internal representation of the model to verify that the additional tasks were learned by the guided architectures. Table 6 evaluates how well the representations of nuclei area, count and contrast are learned by the baseline and the proposed CNNs. For model-ID3 (trained with the uncertainty-based weighting strategy), the prediction of the nuclei count values has average determination coefficient $R^{2}=0.81\pm 0.05$, showing that the concept was learned during training, passing from an initial Mean Squared Error (MSE) of the prediction of $0.46$ to $0.17$ at the end of training. Similar results apply to the other model-IDs 2 to 4 when only a single branch is added. Table 6 compares the performance on the extra-tasks to learning the concepts directly on the baseline model activations, where the network parameters are not optimized to learn the extra tasks. The classification of the center in model-ID 5 reduces in accuracy as the gradient reversal is used during training. The centers of the validation sets are predicted with accuracy $0.29\pm 0.01$ at the end of the training (starting from an initial accuracy of $0.53\pm 0.01$). When more additional tasks are optimized together the performance on the side tasks is affected, with Model-IDs 6, 7 and 8 not reporting high $R^{2}$ values. The average $R^{2}$ of nuclei count for model-ID 6, for example, decreases from $-2.25\pm 0.05$ and plateaus at around $-0.63\pm 0.05$.

Figure 7 shows the dimensionality reduction of the internal representations learned by the baseline and model-ID 3. The visualization is obtained by applying the Uniform Manifold Approximation and Projection (UMAP) method by McInnes et al. (2018) (the hyper-parameters for the visualization were kept to the default values of 15 neighbors, 0.1 minimum distance and local connectivity of 1). The model-ID 3 selected for visualization was trained with the uncertainty-based weighting strategy. In the representation, the two classes are represented with different colors, whereas the size of the points in the plot is indicative of the values of nuclei counts in the images. The top row shows the projection of the internal representation of the last convolutional layer (known as mixed10 in the standard InceptionV3 implementation) of the two models. The bottom row shows the projection of the first fully connected layer after the GAP operation. Since the nuclei count values were normalized to zero mean and unit variance, these are represented in the plot as ranging between a minimum of -2 and a maximum of 2. For clarity of the representation, the image shows the UMAP of a random sampling of 4000 input images.