title: Modality Attention and Sampling Enables Deep Learning with Heterogeneous Marker Combinations in Fluorescence Microscopy
authors: Alvaro Gomariz Tiziano Portenier Patrick M. Helbling Stephan Isringhausen Ute Suessbier César Nombela-Arrieta Orcun Goksel
text: 
  
  
      Modality Attention and Sampling Enables Deep Learning with Heterogeneous Marker Combinations in Fluorescence Microscopy
  
      
  
      Alvaro Gomariz^1,2,*, Tiziano Portenier^1, Patrick M. Helbling^2, Stephan Isringhausen^2, Ute Suessbier^2, 
   César Nombela-Arrieta^2,†, Orcun Goksel^1,3,† 
   
      
  
      ^1 Computer-assisted Applications in Medicine, Computer Vision Lab, ETH Zurich, Switzerland
  
      ^2 Department of Medical Oncology and Hematology, University Hospital and University of Zurich, Switzerland
  
      ^3 Department of Information Technology, Uppsala University, Sweden
      ^* Corresponding author e-mail: alvaroeg@ethz.ch
      ^† These two last authors have equal contributions
  
  
  
  
  
  
  Fluorescence microscopy allows for a detailed inspection of cells, cellular networks, and anatomical landmarks by staining with a variety of carefully-selected markers visualized as color channels. 
  Quantitative characterization of structures in acquired images often relies on automatic image analysis methods. 
  Despite the success of deep learning methods in other vision applications, their potential for fluorescence image analysis remains underexploited. 
  One reason lies in the considerable workload required to train accurate models, which are normally specific for a given combination of markers, and therefore applicable to a very restricted number of experimental settings.
  We herein propose Marker Sampling and Excite — a neural network approach with a modality sampling strategy and a novel attention module that together enable (i) flexible training with heterogeneous datasets with combinations of markers and (ii) successful utility of learned models on arbitrary subsets of markers prospectively. 
  We show that our single neural network solution performs comparably to an upper bound scenario where an ensemble of many networks is naïvely trained for each possible marker combination separately. 
  In addition, we demonstrate the feasibility of this framework in high-throughput biological analysis by revising a recent quantitative characterization of bone marrow vasculature in 3D confocal microscopy datasets and further confirm the validity of our approach on an additional, significantly different dataset of microvessels in fetal liver tissues. 
  Not only can our work substantially ameliorate the use of deep learning in fluorescence microscopy analysis, but it can also be utilized in other fields with incomplete data acquisitions and missing modalities.
  
  
  
  
  § INTRODUCTION
  
  Deep neural networks have been largely successful in many computer vision problems <cit.>, by learning network model parameters in layers to produce feature maps called activations, to arrive at a desired output often given as ground truth during a training phase.
  The learned parameters can then be deployed in an inference phase to make predictions on new input data. 
  Countless advances in this field have occurred in a relatively short timeframe, especially in the use of supervised segmentation, where the goal is the semantic partitioning of an input image with the ground truth typically consisting of pixels annotated interactively by experts.
  For instance, UNet <cit.> is a well-known deep Convolutional Neural Network (CNN) architecture with proven success on semantic segmentation in various biomedical domains. 
  Nevertheless, some aspects of biological images still pose several practical challenges in the application of deep CNN architectures (hereafter also called models).
  
  Fluorescence-based microscopy (FM) is  a mainstay technology for the study of living or fixed tissues in biomedical research. It operates by detecting microscopic signals emanating from inorganic molecules or genetically encoded proteins. 
  Fluorescent dyes are often coupled to antibodies, which target structures or cells of interest within complex samples in a highly specific fashion, a process known as immunostaining. 
  Fluorescent signals are registered and separately encoded as independent image channels due to their distinct spectral properties, thereby allowing the visualization of stained anatomical landmarks of interest. 
  Herein we refer to these channels as markers (also called labels in the literature), which are analogous to the acquisition of modalities in other imaging fields, such as the specific imaging sequences for quantifying different tissue properties with Ultrasound or Magnetic Resonance Imaging. 
  
  The inherent nature of markers in bioimaging studies poses some additional limitations in the creation of datasets that can be processed by typical CNN frameworks. 
  First, the number of markers that can be simultaneously imaged is limited, due to the overlapping spectral profiles of different fluorochromes, which preclude their reliable separation in individual channels.
  Therefore, any detailed characterization of tissues and their pathological perturbations often requires the use of different permutations of a restricted number of markers, which in turn can only provide a limited level of insight into the biological structures studied. 
  Moreover, sample availability is typically a limiting factor, and processing, immunostaining and image acquisition are laborious and time consuming tasks, especially for whole-organ or 3D imaging techniques. 
  Thus, it is not always technically feasible to increase the number of markers, although such additional sources of information would simplify image processing techniques.
  Finally, the process of immunostaining does not always work  consistently, leading to cells and structures stained with variable intensity despite using the same markers.
  Altogether, these issues hinder the generation of a large number of datasets of images stained consistently with all combinations of possible markers. 
  FM datasets thereby often consist of heterogeneous combinations of markers, and with each combination often being limited in number of samples, applications of deep learning become strongly limited. 
  Furthermore, typical supervised segmentation algorithms allow a trained model for only limited future applicability, i.e. when the exact same marker combination is used as in training. 
  This limitation leads to the tradeoff that either a separate specific model is trained each time a new combination is desired and data is available, or a small set of intersection of markers is found in the data; either way neglecting large amounts of precious data and any possibility of using the models later with alternative marker combinations. 
  
  In a general image analysis framework, the problem settings above can be referred to as missing modalities, and are somewhat related to Multi-Task Learning, a field that studies whether information should be learnt jointly or separately <cit.>, and to Domain Adaptation, which aims to bring datasets from different sources into a common space to improve generalization performance <cit.>. 
  It is agreed in both these fields that using a unique model that shares certain amount of information is advantageous. 
  Despite many advances in these fields, the presented missing marker or modality problem is, however, largely unexplored. 
  Recently, synthesis approaches for completing missing data have been proposed for both markers <cit.> and time-sequences <cit.> in FM, as well as for modalities in magnetic resonance imaging <cit.>. 
  Such modality synthesis is cumbersome and potentially sub-optimal when the segmentation model could instead encode information across modalities with shared features. 
  Methods combining different modalities in shared feature spaces were proposed in <cit.> as Hetero-Modal Image Segmentation framework (HeMIS) as well as in <cit.>.
  One would reasonably expect a multi-modal network model to behave differently in the existence or absence of a particular modality.
  Such processes of conditioning the models explicitly are known as attention mechanisms.
  For example, soft attention mechanisms transform the activations of a model conditioned on the activations themselves <cit.>.  
  Notably, Squeeze and Excitation (SE) <cit.> and similar modules <cit.> have been very successful and since been integrated in several different network architectures to improve their performance and the interpretability of extracted features. 
  
  Building on these ideas, we herein devise a method that addresses the fundamental problem of multi-modality heterogeneous sets and evaluate this on a 3D microscopy image dataset of bone marrow vascular network <cit.>, where the annotations are divided into two vascular types, namely sinusoids and arteries. 
  The variable size and morphology of the vasculature had hindered precise segmentation and thus a reliable vasculature characterization in earlier works <cit.>.
  In this experimental setting with 5 FM imaging markers, we first evaluate multiple baseline conventional scenarios to analyze the effect of each possible marker combination on the performance of semantic (not to be confused with instance) vessel segmentation, which also serves as an upper bound for assessing our proposed methods. 
  Then, we show that a Marker Sampling strategy enables a single CNN model to successfully perform in the presence of any marker combination, while outperforming the current state-of-the-art HeMIS. 
  Next, we present our novel Marker Excite soft attention module, which learns how to recalibrate the network activations as a function of available markers, showcasing this by training a single model that performs comparably to the upper bound scenario with an ensemble of 31 separate models individually trained and specialized for each of the marker combinations. 
  We further demonstrate that our model can even outperform the upper bound, by leveraging information shared across markers when the training dataset contains such practical variations. 
  Next, we present a case study  on the problem setting of <cit.> to show the application of our method on an existing practical research question. 
  Finally, we demonstrate the widespread applicability of our methods by applying them to a largely different and independently generated fetal liver FM dataset, thereby providing the first quantitative characterization of microvascular networks in this organ across different embryonic stages.  
  
  
  
  § RESULTS
  
  
  
   §.§ Segmenting fluorescence microscopy samples stained with multiple markers
  
  
  
  
  Markers in FM label specific biological structures, and their efficient combination enables the visualization of distinct cellular/subcellular components or networks thereof.
  Manual annotation of these structures as classes is usually possible when the available markers accurately portray them. 
  Meanwhile, image segmentation algorithms can target these classes by employing arbitrary combinations of markers, but their performance will largely vary depending on the combination employed. 
  To be able to achieve such segmentation with different markers, we designed a neural network-based image processing pipeline applicable to tissue-wide FM imaging (illustrated in Fig.  <ref>a with details in the Methods section) that processes markers as channels. 
  In this work we employ a dataset with 8 large samples divided in a number of 2D patches (Supplementary Table <ref>). 
  Each patch consists of two ground truth classes (sinusoids and arteries) and five markers (DAPI, endomucin, endoglin, collagen, and CXCL12-GFP) denoted as m_G,   G ⊆{1, …, 5}. 
  Hereafter, we use multiple subscripts successively to indicate combinations of these markers. 
  All following CNNs results were validated with 4-fold cross-validation, and  claims are made only when the respective null hypothesis can be rejected with a p-value≤0.05. 
  More details are given in the Methods Section. 
  
  We build an extensive semantic segmentation baseline by training a distinct UNet separately for each of the 31 possible combinations of 5 markers.
  These also serve in the following sections as upper bound performance given this architecture.
  We present results for models trained separately for two classes of vessels, since empirically (measured with F1-score in Fig.  <ref>b,f) this performed  superior to training them simultaneously. 
  To analyze marker importance, we rank in Fig.  <ref>c-e&g-i the marker combinations according to their mean F1-score.
  Although the segmentation accuracy is seen overall to increase with more markers, it is seen to also highly depend on specific markers: 
  Sinusoid segmentation in Fig.  1c-e has higher accuracy when either m_2 or m_3, which specifically label these structures, are present (blue). 
  Without them, m_4 helps with segmentation (red), compared to the least helpful m_1 or m_5 (green). 
  For arteries in Fig.  1g-i, having the arterial-specific m_4 marker seems essential (blue), and results can be further improved by adding m_2 (red). 
  These observations are used later herein to interpret outcomes, e.g. when some suitable markers are missing.
  
  
  
  
   §.§ Marker Sampling for segmenting with missing markers
  
  
  
  
  The strategy described above allows to  determine the best combination of markers to segment a given class.
  Nevertheless, the need for a fixed marker combination would be restraining in different applications due to the practical limitation on number of markers to use simultaneously in sample preparation. 
  It is also not fail-safe if individual markers fail during acquisition. 
  Moreover, such a naïve approach of training a distinct CNN for each combination would require 2^K - 1 models for K markers (31 models in our study), an exponential growth, which is prohibitive for training CNNs in reasonable time frames.
  Finally, retraining this many models becomes highly impractical as new samples are added to the training dataset. 
  To address this challenge we propose a Marker Sampling (MS) approach (illustrated in Fig.  <ref>a). 
  A key component of MS is a sampling layer at the input of a segmentation network, herein UNet. 
  During training, we provide all the available markers to the network, while this sampling layer randomly selects a subset to be processed by the proceeding network. 
  This single-model framework consequently can learn to generalize to any subset of markers at inference.
  
  Since the accuracy for each marker combination may vary widely (also seen in Fig.  <ref>c,g), we herein report comparative improvements in a paired test manner; i.e. 
  when we test a hypothesis with respect to an alternative method (indicated in the figures as “Rel. to”), we calculate the relative metric difference separately for each marker combination experiment and report the distribution and statistical significance of such differences (with details in the Methods Section).
  Fig. <ref>b shows a comparison of our MS method with the state-of-the-art HeMIS.
  In addition we train our single CNN without MS, i.e., the network always uses all markers during training without sampling. 
  When performing inference on a subset of markers, we simply set the missing input to zero. 
  We refer to this simplistic baseline as Marker Zero (MZ).
  The results show that our MS model vastly outperforms both HeMIS and MZ, indicating that training a network simply with randomly sampled marker subsets generalizes across the possible input combinations better than the other two approaches. 
  Surprisingly, MZ does not perform significantly different than HeMIS, suggesting that the latter does not provide any advantage for this task compared to a standard UNet architecture. 
  Based on this, learning shared features can be considered comparable to marker-specific representations for this purpose.
  For completeness we also evaluate incorporating the MS strategy in the HeMIS architecture (Fig.  <ref>c) and find that there is no improvement compared to UNet-based MS, therefore concluding that any advantage originates from our MS strategy per se.
  
  It could be beneficial to normalize the output of the sampling layer across channels to keep its signal magnitude constant. 
  In fact, MS can be interpreted as a variant of Dropout regularization of neural networks <cit.>, where it is applied only on the input layer and, not only during training but also during inference (deterministically via available markers). 
  We investigate such relation by training two variants of our model: Marker Sampling with Dropout (MS-DR), which scales the intensities of the available images by a constant Dropout ratio, and Marker Sampling with Variable Ratio (MS-VR), which scales the intensities by the ratio of available markers in each sample. 
  Among these, MS is found to be the best performing model (Fig.  <ref>c), indicating that our improvements cannot be matched by such handcrafted normalization in the sampling layer. 
  
  In summary, while HeMIS, MZ, and MS all perform well when all markers are present (m_12345), our proposed MS generalizes significantly better when some markers are missing at the time of segmentation. 
  Indeed, for instance, with m_4 alone our proposed MS can already segment satisfactorily (Fig.  <ref>d), which the baselines are unable to.
  
  
  
   §.§ Marker Excite for learning attention to markers
  
  
  
  
  As demonstrated above, MS can effectively generalize across marker combinations using a traditional neural network backbone such as UNet to learn shared features for distinct markers. 
  We herein additionally propose Marker Excite (ME) — a novel attention module for a further boost by learning a weighting of deep features as a function of available markers.  
  Note that existing attention modules such as SE learn such weighting as a function of layer activations.
  Instead, our ME module learns attention to marker availability, provided explicitly as an additional input to the network in the form of a one-hot encoded symbolic vector (Fig.  <ref>a). 
  We integrate ME modules at different layers of a UNet as shown in Fig.  <ref>b.
  Integrating our ME approach with the previously presented MS strategy of prepending a sampling layer, we attain our ultimate model MS-ME.
  We compare this below also using the conventional attention by replacing ME with SE, a baseline referred hereafter as MS-SE.
  
  Fig.  <ref>c shows that MS-ME yields improved overall accuracy with respect to MS as well as MS-SE, especially for the sinusoids. 
  A major advantage of attention modules is that they increase the model complexity only marginally (MS-ME having merely 0.64% more parameters than the UNet architecture it was based on). 
  To demonstrate that the presented improvement does not originate from inflation of model size, we conducted an additional experiment with a much larger UNet baseline (MS+) of over 20% more parameters and show that the clear improvement from our proposed approach persists (Fig.  <ref>d). 
  We also performed an ablation study by placing ME attention at different network layers (Fig.  <ref>f,g), concluding that placing ME after every convolutional block (outlined in orange in Fig.  <ref>b) is the optimal configuration (namely MS-ME) for our task.
  An analysis of the influence of the proposed ME modules reveals that the recalibration effect is stronger for activations of higher resolution (Supplementary Fig.  <ref>).
  
  We also compare our proposed MS-ME model with an upper-bound (UB) presented in Fig.  <ref>c,g.
  Note that UB consists of 31 individual models each of which separately trained and specialized to the availability of a specific marker combination, whereas MS-ME is a single model aiming to achieve well on whatever marker combination may be available at a given time.
  Fig.  <ref>e shows that across all markers combinations our MS-ME is not significantly different than UB, and even slightly superior for sinusoid segmentation potentially thanks to leveraging additional information shared across combinations. 
  A sample qualitative result is provided in Fig.  <ref>h for a marker combination (m_15) known to be suboptimal (based on Fig.  <ref>c,g) with a specialized training of a dedicated UB model.
  Even for this case, MS-ME is seen to still perform somewhat satisfactorily, while UB fails completely.
  
  
  
  
   §.§ Training with heterogeneous panels of markers
  
  We next study the scenario with an incomplete training set, i.e. with acquisitions of heterogeneous combinations and different number of markers — a typical setting in the field in practice.
  For this purpose, we used subsets of our fully-annotated dataset.
  Training a separate network of all combinatorial test settings would be computationally prohibitive, thus we artificially ablate the data to create a number of case studies (Fig.  <ref>) and thereby emulate various practical or extreme scenarios. 
  
  
  
  In Cases 1 to 4 we studied settings in which training data contain a fixed number of markers per sample, each with a different marker combination. 
  Case 5 simulates a common scenario where two different sets of samples are prepared with two different staining protocols.
  Case 6 studies a practical scenario with samples available from different studies, which have so far been largely unusable for machine learning due to their heterogeneity.
  
  While our proposed CNN models are by construction applicable for the above scenarios, any baseline segmentation methods such as UNet require an overlapping set of markers. These baseline models are therefore trainable only for a few intersections of marker combinations, and discard any potentially useful information outside such common intersections.
  For instance, the marker combination m_23 in Case 3 is only trainable using samples s1 or s2, and in Case 1 not even trainable.
  
  With our proposed models, in all cases any marker combination is trainable. 
  Therefore, for numerous test combinations, our proposed models are inherently superior by design, since the baselines cannot even accommodate test combos unseen during training.
  For many other tests, the baselines would be using a small intersection subset of samples, again at a major disadvantage and also forcing us to retrain a network for each of these subsets.
  We herein compared our methods (Fig.  <ref>) to the UB models only for marker combinations achievable by baseline in training. 
  As this clearly biases results towards UB, we also report herein the ratio of marker combinations subject to such evaluation, i.e. available in the training set.
  Small ratios in most cases demonstrate that our contribution is not only superior in F1-score, but also in the ability to include and utilize datasets hitherto impossible.
  
  The superiority of MS over MZ and of MS-ME over MS is emphasized with increasing number of markers available for individual samples, since there are more combinations from which information can be leveraged through sampling and attention. 
  Thus, MS and MS-ME do not add any advantage only when all samples have a single marker (Case 1). 
  Remarkably, when considering both sinusoids and arteries, MS-ME outperforms UB in Cases 3, 4, and 6, where there is more heterogeneity of markers across samples. 
  Thus, best relative gains of MS-ME are observed with more heterogeneity and number of markers across samples. 
  We see that MS-ME not only requires a single model and is applicable for several marker combinations previously unattainable, but it also presents a superior segmentation performance, even over the UB for several cases.
  
  
  
   §.§ A practical application on segmenting bone marrow vasculature
  
  As a showcase of the proposed methods in an application scenario, we herein revise our segmentation results from <cit.>, where the use of traditional Morphological Image Processing (MIP) algorithms had restricted our analyses both in terms of the number of samples included as well as in the accuracy of quantification. 
  As indicated in Table <ref>, the herein proposed MS-ME (i) permits full automation of the analysis, (ii) substantially increases F1-score by 47.3±9.6%, and (iii) almost quadruples the number of samples that can be successfully processed, allowing the inclusion of 35 samples which could not be previously analyzed due to marker heterogeneity (Fig.  <ref>b) or  insufficient image quality for MIP. 
  Segmentation examples with MS-ME can be seen in Supplementary Fig.  <ref>.
  
  
  
  We separately quantify the vascular volumes in two distinct regions of the bone marrow, namely diaphysis and metaphysis, shown in Fig.  <ref>a. 
  Our revision of the previous  analysis <cit.> indicates that the volume occupied by sinusoids is significantly lower than we previously reported (Supplementary Fig.  <ref>).
  The results in Fig.  <ref>c also confirm the previously observed trend (with more statistical power) that they are more abundant in the metaphysis (15.65 ± 2.55 %) than in the diaphysis (10.69 ± 2.16 %). 
  The differences can be explained by the higher accuracy of our new method. 
  In addition, our methods allow for an accurate segmentation of much rarer arterial networks, which to the best of our knowledge have only been quantified in very limited tissue regions in previous studies <cit.>. 
  We find that these vessels occupy 0.63 ± 0.40 % of the volume in the diaphysis and 0.52 ± 0.32 % in the metaphysis.
  
  
  
  
  
  
   §.§ Marker Sampling and Excite for the characterization of fetal liver vasculature
  
  
  
  
  To further demonstrate the ability to generalize across samples and the applicability of our proposed methods on significantly different biological tissues in which different marker combinations are available, we provide a comprehensive quantitative analysis of a murine fetal liver FM dataset.
  The segmentation label (vessels) and the 5 different markers employed for this dataset, largely divert from those employed in bone marrow and are illustrated in Fig.  <ref>a, denoted as m_G,   G ⊆{a, …, e} (details in Methods section - Segmentation of fetal liver vasculature). 
  
  Evaluating the UB for each of the possible marker combinations shows that a superior F1-score is always achieved when m_b is available (Fig.  <ref>b).
  If all markers are available at training time, our models behave similarly as with the bone marrow dataset; i.e., MS is superior to MZ and MS-ME is superior to MS in terms of F1-score evaluated across all possible marker combinations at test time (Fig.  <ref>c). 
  We further evaluate the segmentation performance when the training samples have different combinations of markers by artificially assigning a random marker combination to each of the samples. 
  The results in Fig.  <ref>d show again the superiority in F1-score of MS-ME relative to MS and MZ.
  In this case, MZ has a comparable F1-score to that of MS, probably because the ablation strategy employed in this setting already creates enough training examples with the different combinations of markers, rendering sampling unnecessary. The F1-score of all three models is comparable to that of UB in this setting, corroborating results from some training settings in Fig.  <ref> for the bone marrow dataset.
  
  Next, we follow the analysis strategy of the previous section to demonstrate the potential of our framework in the characterization of different biological samples.
  We employ MS-ME to segment the liver vasculature in 15 fetal liver samples at two different embryonic stages denominated E13.5 and E17.5 as illustrated in Fig.  <ref>e, each having different marker combinations (Fig.  <ref>f).
  The results in Fig.  <ref>g reveal isometric vascular growth during these embryonic stages, as no significant differences in vascular occupancy between timepoints is observed, with a vascular ratio of 19.4±3.60% across all samples.  
  
  
  
  § DISCUSSION
  
  In this work we formalize the widespread challenges concerning segmentation of structures in FM datasets with different markers, and address them with unprecedented accuracy and efficiency.
  We first studied the contribution of all marker combinations by naïvely training a distinct CNN model for each combination. 
  The segmentation results were ranked for each of the possible markers employed, of which even the worst combinations could produce meaningful results.
  Albeit useful as an upper bound performance, training this number of models is infeasible in practice. 
  We hence proposed a number of methods that allow for a practicable workflow with a single model that can operate with any combination of markers both during inference (application-time) and training.
  
  During our preliminary experiments, the fact that a simple UNet adaptation MZ could perform as well as HeMIS, the state-of-the-art deep learning approach for missing modalities, indicated to us that the solution to this challenge was potentially not in the design of novel network architectures, but rather in adapting suitable training strategies that can inherently address such problem. 
  Accordingly, we devised MS as a sampling strategy, the results of which confirmed our hypothesis, providing a drastic improvement at inference for segmenting different structures with any marker combination.
  We also tested normalization strategies for the sampled markers, but none helped: MS-DR may have failed because it scales image intensities by a constant factor at training time, which makes sense in the original Dropout idea that samples weights as a regularization technique, but sampling in our missing marker setting intrinsically occurs at inference too. 
  MS-VR, which scales intensities by the ratio of available markers, did not improve the results either. The reason could be that such scaling is reasonable as regularization in the case of marker-independent features, whereas the employed network may learn different features for each marker.
  
  Results with our novel ME method show that attention mechanisms do not have to be limited to layer activations, but using other sources of information as weak labels can also boost outcomes substantially with only a slight increase in model complexity.
  When our two contributions were combined as MS-ME the best segmentation results were obtained, performing as well as an ensemble of 31 specialized UB models.
  In addition, when training samples had different markers, MS-ME was applicable to combinations inaccessible with typical CNNs such as UNet or HeMIS, while our MS-ME segmentation accuracy surpassed even an upper bound setting, especially when marker heterogeneity was high.
  This surprising result, as a far-reaching consequence, shows that the solution we originally devised to accommodate missing markers can easily help to improve outcomes even in traditional learning scenarios, potentially by leveraging and incorporating complementary information from multiple sources. 
  
  When applied to the bone marrow vasculature dataset, our approach based on MS-ME showed a dramatic increase in accuracy, speed, and number of samples which could be processed, indicating the importance of chosen image analysis techniques in scientific studies.
  We therewith increased the statistical power of the characterization reported in <cit.>, allowing us to confirm the previously hinted tendency of sinusoids occupying larger volumes of the metaphysis than the diaphysis.
  Application of the same analysis pipeline to a dataset depicting fetal liver vasculature with a different marker set, further allowed to quantitatively describe the vascular occupancy in different embryonic stages. 
  Our proposed techniques will be instrumental in imaging studies to accurately study the morphometric and structural features of microvascular networks in hematopoiteic organs, to gain key deeper insight on their described fundamental roles in the pathophysiology of blood-forming tissues.
  Beyond this, vascular segmentation in large multidimensional datasets remains a widespread and largely unresolved challenge for many groups using advanced microscopy technologies. The methods presented may reveal as extremely powerful tools for studies aimed at uncovering novel mechanisms regulating tissue and cellular dynamics in almost any biological tissue, from animal models, human samples, to model organoids grown in vitro <cit.>.
  
  Despite the superiority of MS-ME in segmentation results without compromises in number of parameter or inference speed (Supplementary Table <ref>), some limitations exist. 
  For instance, although MS-ME can be applied to previously unseen combinations of markers, quantitative quality assessment is only possible if labeled examples exist featuring the respective combination.
  In addition, although not observed in our study, our MS strategy may amplify existing imbalances of markers in the training set, increasing the danger of overfitting to specific combinations.
  
  These challenges open a number of interesting questions for future work. 
  For instance, novel uncertainty estimation methods for deep neural networks <cit.> can be employed as proxy for assessment of segmentation quality upon previously unseen combinations. 
  In addition, although overfitting to specific markers can be naïvely addressed with a validation set containing as many markers as possible, investigating how to adapt the dropout rate of MS to different experimental settings may lead to solutions with superior results.
  Lastly, our MS-ME model may enable novel possibilities in the context of transfer learning <cit.> that enables seamless application of trained networks to distinct biological tissues with new markers. 
  For example, fine-tuning only the few parameters in the ME module as opposed to expensive, classical fine-tuning strategies may yield competitive results if the core network already learned suitable general purpose features.
  
  
  A major challenge and limitation for utilizing deep learning in FM has been the difficulty in establishing standardized staining protocols that would enable more homogeneous marker combinations to train supervised models.
  With our methods proposed herein, a single model is shown to perform comparably or superior to a number of individual problem-specific models that would be infeasible in practice due to the exponential growth in model parameters and training time with an increasing number of markers. 
  In addition, the versatility of our methods enables them to be easily applied to different network architectures for tasks beyond semantic segmentation, such as classification <cit.>, detection <cit.>, or instance segmentation <cit.>.
  These contributions can facilitate the sharing and exchange of trained CNNs across labs in the field as well as a faster adoption of neural network solutions in routine lab work at, e.g., microscopy facilities. 
  
  
  
  
  § METHODS
  
  
  
   §.§ Dataset composition and notation
  
  The dataset S under evaluation consists of 8 samples s^i,   i ∈{1, …, 8} and each sample s^i is composed of J patches  p^i_j,   j ∈{1, …, J}. 
  Samples are prepared with a set of different markers k ∈{1, …, K}, where K is the total number of markers. Marker combinations are denoted as m_G, where G is a non-empty subset of available markers, i.e. G ⊆{1, …, K}, denoted in the subscript as a successive sequence, e.g. m_12345 indicating the combination of all markers for K=5 in our study. 
  Each sample s_i was manually annotated for C classes c ∈{1, …, C}. Thus, each patch p^i_j consists of | G | input images x_k ∈ℝ^h_in× w_in with k ∈ G and two output annotations y_c ∈ℝ^h_out× w_out as follows: 
  
      p^i_j = {x_k, y_c} ∀ k ∈ G,  ∀ c ∈{1, …, C}
  
  All patches are of the same size, i.e. h_in = w_in = 572 and h_out = w_out = 388.
  
  The complete bone marrow dataset S (Supplementary Table <ref>) was prepared using the following 5 markers: DAPI stains DNA in all nuclei, endomucin and endoglin have both been reported to have a high specificity for bone marrow sinusoids, collagen mostly stains vessel walls (including sinusoids and arteries), as well as extracellular matrix, and GFP (in the CXCL12-GFP mouse) is a genetically encoded marker which stains reticular mesenchymal stromal cells. 
  We use C=2 annotated classes: sinusoids and arteries (the latter accounting for both central larger arteries and endosteally located smaller arterioles). 
  The immunostaining and imaging protocols are detailed extensively in <cit.>.
  Note that such a dataset with a fixed combination of markers and classes is challenging to obtain in practice. 
  We artificially ablate parts of the data as described in the different experiments to simulate different realistic scenarios of marker and class combinations. 
  The fetal liver dataset is separately described in the Segmentation of fetal liver vasculature subsection. 
  
  
  
  
   §.§ Marker combination strategies
  
  We study herein two fundamentally different settings. First, we study a relatively simpler scenario where training data is acquired with a consistent, fixed set of available markers m_1, …, K; while at test time segmentation is requested for samples only with a subset of trained markers. 
  In this setting we ablate subsets of markers to generate test samples with missing markers, i.e.: m_G  ∀ G ⊆{1, …, K}.
  Note that such training data with many markers are scarce as they require extensive workforce, effort, and budget; not only for the sample preparation and acquisitions, but also for their annotations.
  
  Next, we study the more common (and more challenging) scenario where a training set with a fixed set of markers is unavailable.
  Since it is computationally not feasible to evaluate each possible set of markers during training, we propose to evaluate a number of illustrative cases (Fig.  <ref>). 
  To this end, we predefine a set of marker combinations M^𝑡𝑟𝑎𝑖𝑛 so that specific combinations m_G^i ∈ M^𝑡𝑟𝑎𝑖𝑛 are assigned to each sample s^i in the training set.
  The test set is constructed with all possible combinations of markers as above, i.e. m_G  ∀ G ⊆{1, …, K}. 
  To avoid overfitting to a specific marker combination, we use a validation set with samples that contain all the markers, i.e. m_1, …, K. 
  In addition,  M^𝑡𝑟𝑎𝑖𝑛 is shuffled for each cross-validation step so that different m_G^i are applied to different samples.
  
  The resulting segmentation performance understandably varies largely across different marker sets m_G. 
  Such variation in any metric would be difficult to interpret and to use to compare the models. 
  We therefore report relative changes, i.e. the score Q achieved by a model ϕ on a marker combination m_G relative to a reference model ϕ_ref as Q_ϕ, ϕ_ref(m_G) = Q_ϕ(m_G) - Q_ϕ_ref(m_G). 
  Consequently, a resulting vector encodes metric differences of a method with respect to a reference, on a given marker combination m_G.
  To put such results in context, we also report results referenced to an ideal Upper Bound model (UB) as Q_ϕ, ϕ_UB(m_G) = Q_ϕ(m_G) - Q_ϕ_UB^m_G, where ϕ_UB^m_G is a UNet model trained exclusively on m_G.
  Note that such UB reference results Q_ϕ_UB^m_G can only be computed for marker combinations explicitly available in at least a sample, i.e. for marker combinations
  
      m_G ∈ M^UB,   where  
           M^UB = {m_G   |   G ⊆ m^i,  ∀ m^i ∈  M^𝑡𝑟𝑎𝑖𝑛} .
  
  In contrast, Q_ϕ with the proposed models can be computed on any marker combination, i.e.
  
      m_G ∈ M^total,   where  
          M^total = {m_G   |   G ⊆⋃_m^i ∈ M^𝑡𝑟𝑎𝑖𝑛 m^i }
  
  In other words, the proposed methods can produce results even for combinations m_G that never occur in any training sample, as long as the individual markers are present in at least one sample.
  Thus, the ratio |M^UB| / |M^total| is also reported to indicate the maximum possible marker combinations achievable by conventional networks. 
  Lower ratios then emphasize that, besides requiring comprehensive training sets, standard neural network models for segmentation can only be utilized for few marker combinations, while our proposed model can accommodate all combinations.
  
  
  
  
   §.§ Network architectures
  
  The network architectures employed in this work are described below, first the two baseline models adapted from previous work, followed by our proposed models:
      
          
    * UNet <cit.> is one of the most commonly used segmentation models of biomedical images. It is a CNN based on an encoder-decoder architecture with skip connections, targeted to extract features at different resolution levels. We used its standard settings, except for decreasing the number of parameters in each layer by a factor of two, having empirically found that this produces superior results on our dataset. 
          
          
    * HeMIS <cit.> is the current state-of-the-art for image segmentation with missing modalities. Each marker is processed by separate models (backends) that are subsequently combined with mean and standard deviation (abstraction layer). Such feature aggregation using statistical moments allows to seamlessly apply backpropagation regardless of the missing modalities. The merged output is then processed with additional convolutional layers (frontend) to obtain the final segmentation. 
  
          
    * Marker Zero (MZ) is a UNet where missing markers are padded with zeros, i.e., for any marker combination m_G, an input image
          
      x_k ←
                      x_k    if  k ∈ G
      
                      0    otherwise .
  
          
          
    * Marker Sampling (MS) is a UNet prepended with a sampling layer that randomly deletes markers at training time with probability r_drop∈ (0, 1), i.e. x_k ←𝐵𝑒𝑟𝑛(r_drop) x_k, where 𝐵𝑒𝑟𝑛(r) denotes random sampling from a Bernoulli distribution with probability r. 
          r_drop=0 would be equivalent to MZ, and r_drop=1 is excluded because that would create an input image without markers. 
          Note that no such sampling occurs at inference, but it is implicit to the application strategy in that samples with different markers are expected. 
          We set r_drop = 0.5 in all our experiments, hence all marker combinations are sampled with equal probability.
          In-depth analysis revealed that this choice leads to the best overall accuracy (Supplementary Fig.  <ref>a).
          Note that in this study we consider the scenario where all marker combinations are equally likely to occur at inference, which justifies the choice of setting r_drop=0.5. 
          However, if combinations of specific number of markers are expected to occur more frequently at inference, r_drop can be adjusted accordingly (Supplementary Fig.  <ref>b), potentially also incorporating marker-specific preferences.
          
          
    * Marker Sampling - Dropout (MS-DR) is similarly to MS, where the sampling layer is replaced with a Dropout layer, i.e. x_k ←r^-1_drop𝐵𝑒𝑟𝑛(r_drop) x_k only at training time.
          
          
    * Marker Sampling - Variable Ratio (MS-VR) is a modification of MS where the sampling layer normalizes the output according to the number of sampled markers, i.e. x_k ← p^-1𝐵𝑒𝑟𝑛(r_drop) x_k at training time, with p being the ratio of markers chosen by the sampling layer.
          In this case, the normalization p is also applied at inference time, i.e. x_k ← p^-1 x_k, since it it is a function of the number of utilized markers, instead of an expected probability as in MS-DR.
          
          
    * HeMIS - Marker Sampling (HeMIS-MS) is a HeMIS network that is prepended with our MS sampling layer. In contrast to zeroing out missing markers as in the UNet counterpart, we herein ablate any backend corresponding to a missing marker, in order to avoid statistical moment calculations from empty feature maps.
      
          
    * Marker Sampling - Squeeze And Excitation (MS-SE) is an MS model with a Squeeze and Excitation (SE) module as proposed in <cit.> in each UNet block (see Fig.  <ref>b). 
          These modules each learn a weighting function δ: ℝ^1 × F→ℝ^1 × F of a feature map X with F activations, where each activation layer X_f∈{1, …, F} is weighted as a function of its layer-wise mean value ⟨ X_f ⟩, i.e. X_f ←δ_f(⟨ X_f ⟩) X_f.
          δ is parametrized using two fully-connected network layers, respectively, with F/2 and F nodes and no biases. The first layer is followed by a rectified linear unit (ReLU) and the second one by a sigmoid function.
          
          
    * Marker Sampling - Marker Excite (MS-ME) instead weights the activations using our proposed Marker Excite (ME) module (illustrated in Fig.  <ref>a) that learns a weighting from a binary vector V encoding which of the K markers are presented to the network as input. 
          Our module thus learns a function ω: {0, 1}^K→ℝ^1 × F to weight the activations as X_f ←ω_f(V) X_f. Similar to SE we parametrize ω using two fully-connected layers. 
          The first layer has as many nodes as possible marker combinations (2^K-1) and we employ biases in both layers, which is seen to increase F1-score (Supplementary Fig.  <ref>). 
          Similarly to MS-SE, we add ME modules after each UNet block (Fig.  <ref>b).
      
      
  We do not apply Batch Normalization in any of the models as we empirically found this to decrease the F1-score, presumably due to the relatively small batch sizes. 
  Supplementary Table <ref> shows that the HeMIS architecture has less parameters than UNet, but a higher memory footprint and lower speed. 
  Meanwhile, our MS and ME methods have little to no additional burden, neither in memory space nor in speed. 
  
  
  
   §.§ Training and evaluation details
  
  We implement all CNN models in TensorFlow 2.1 <cit.> and train them on an NVIDIA GeForce GTX TITAN X GPU with 12 GB of VRAM, with 3 Intel Xeon E5-2640 v3 CPU cores and 40GB of host RAM. 
  Due to large image sizes and GPU memory limitations, we use a rather small batch size of 2. 
  Using the Adam optimizer <cit.>, we minimize the following weighted cross-entropy loss:
  
      L = ∑_c ∈ C W_c ∑_y_j^c ∈ Y^c y_j^clog(y_j^c),
  
  where y_j^c is the ground truth annotation for class c, y_j^c is its corresponding network prediction, and W_c is a class-specific weight to account for class imbalance. 
  We found that using a weight that is linear in the class cardinality leads to training instabilities, especially with large class imbalance. 
  Therefore we instead use a logarithmic weighting as follows:
  
      W_c = log(ϵ + 1/C |Y^c|∑_i ∈{1, …, C} |Y^i| ),
  
  where |Y^c| is the total number of annotated pixels for class c and ϵ is a small number for stability in logarithm calculation.
  
  The results are evaluated using F1-score, individually for each of the classes, as follows:
  
      F1 = 2TP/2TP + FP + FN
  
  All networks are trained for 200 epochs, and we choose for evaluation the one that yields the highest mean F1-score across all classes on the validation set.
  The evolution of the loss function and the F1-score is shown in Supplementary Fig.  <ref>.
  Note that there exists a considerable class imbalance (see Supplementary Table <ref>), which can explain the inferior F1-scores when classes are targeted simultaneously (Fig.  <ref>b,f) as compared to targeting them individually, which is why we chose the latter approach to conduct our experiments.
  There has been substantial research in class imbalance <cit.>, but this is not within the scope of this paper, where we focus on effects of markers, i.e. input image modalities.
  
  The 230 image patches were split according to which of the 8 FM samples in Supplementary Table <ref> they belonged to, with 5 samples for training, 1 for validation, and 2 for testing. 
  For all experiments we employ 4-fold cross-validation, by ensuring that all samples are once in test set and that the same sample folds are used for different methods in each experimental setting such that their results can be directly compared. 
  
  
  
   §.§ Tiling strategy for handling large datasets
  
  Wide-tissue FM images typically have very large pixel sizes, exceeding memory capacity of typical GPUs used for deep learning. In addition, FM samples are often acquired with different resolutions. 
  To address these problems, we herein adopt a pipeline similar to those included in recent UNet variants <cit.> by decomposing each sample s^i in J patches p^i_j,   j ∈{1, …, J} with constant size and resolution to segment them individually and to subsequently stitch them together.
  These patches are constructed as follows (illustrated in Supplementary Fig.  <ref>):
  
      
    * We resample the resolution of all FM samples to 1μm pixel size.
      
    * Since convolution operations decrease the spatial extent of the image, we first zero-pad the complete sample with a margin of 92 pixels (half of the difference between input and output patch sizes) required to preserve the original size when stitching the output patches. 
      
      
    * The input patches are taken with an overlap of 92 pixels (same as the padding) between them so that there is no overlap in the output patches. In this way, we limit padding artefacts to the border of the samples instead of introducing them for each individual patch used in the CNNs.
      
      
    * To normalize appearance differences across samples, for each patch we apply Gaussian standardization (zero mean, unit variance) using the mean and standard deviation of the respective sample that the patch comes from.
      
      
    * Output patches are neither overlapping nor padded. Since the sample size does not have to be divisible by the patch size, the last patch of each dimension was taken with as much overlap with the previous one as needed to cover the whole sample. At inference, the result for the overlapping region is averaged for the involved patches. Finally, the slices reconstructed by 2D tiling are inserted in their corresponding axial position to form 3D images. 
   
  
  
  
   §.§ Bone marrow quantification pipeline
  
  For the results reported on the characterization of bone marrow vasculature we employed additional unlabeled samples from the dataset described in <cit.>. 
  In that earlier work, several samples had to be discarded because their image quality was not sufficient for the simplistic image processing tools (MIP) employed therein. 
  
  The MS-ME model proposed herein allows for processing of samples with diverse marker combinations. 
  However, as shown in Fig.  <ref>c-e&g-i, not all markers are beneficial for achieving a precise segmentation. 
  In order to provide an accurate quantification without sacrificing many samples, we employ in our analysis only samples that are stained with the marker sets in M_2^s for sinusoids and in M_2^a for arteries, as shown respectively in Fig.  <ref>c and Fig.  <ref>g. 
  Quantified on the samples stained with the best marker combination, the above marker sets guarantee an F1-score of 75% or higher than that achievable given the best marker combination. 
  With this criterion, 47 samples are employed for the quantification of sinusoids and 29 for arteries. 
  
  All employed samples were visually inspected to qualitatively confirm that the segmentation was satisfactory (examples in Supplementary Fig.  <ref>).
  The F1-score reported for MS-ME in Table <ref> is calculated from the results reported in Fig.  <ref>e aggregated over the markers in M_2^s for sinusoids and in M_2^a for arteries. 
  The MIP method was evaluated on the same annotations as the other methods proposed in this work (Supplementary Table <ref>), although not all of these samples were employed for the quantification in <cit.>.
  
  FM images often contain out of tissue regions which must be discarded for analysis purposes. To this end, we account in this analysis for an extra class denoted as tissue.
  In several works, segmentation of this class is easily achieved by thresholding of the marker m_1 (DAPI) and some simple morphological image processing <cit.>. 
  Given the simplicity of the task, we train a separate UNet model that uses only m_1 (available in all the samples) and achieves an F1-score of 88.7±6.8 with the same evaluation strategy applied for other models in this work. 
  The reported fraction of bone marrow volume for sinusoids and arteries is calculated as the ratio of their respective foreground pixels within the newly defined tissue class to the total number of foreground tissue pixels. 
  
  
  
   §.§ Segmentation of fetal liver vasculature
  
  For the fetal liver dataset presented in Results section Marker Sampling and Excite for the quantification of fetal liver vasculature, we employed the same strategies previously described for the earlier bone marrow dataset.
  We herein report the dataset description and technical details that differ from the bone marrow images. 
  
  We denote the 5 available markers with letters instead of numbers (i.e. the complete marker set is denoted as m_abcde): DAPI (a), lyve1 (b), Hlf (c), Evi1 or α-catulin (d), and smooth muscle actin (e). 
  The class to segment in this dataset is denoted as vessels.
  The number of samples, patches, and percentage of annotated pixels can be seen in Supplementary Table <ref>.
  
  The vasculature quantification for different embryonic stages denoted as E13.5 and E17.5 displayed in Fig.  <ref>g is performed on 15 unlabeled samples: 8 samples of E13.5 and 7 samples of E17.5.
  A manually annotated 3D tissue mask is included in these images that defines the total tissue volume, which is employed to disregard out-of-tissue regions in the analysis of vasculature occupation. 
  Note that that the annotated fetal liver vasculature data employed for training and evaluation of segmentation CNNs is a subset of these samples. 
  
  
  
   §.§ Statistical tests
  
  Unless otherwise specified, the two-sided Wilcoxon signed-rank test was employed to assess differences between paired test results, which is non-parametric to avoid any assumption of normal distribution. 
  Since the aim of the models was to perform best across different classes, the employed tests aggregated the data for both sinusoids and arteries in a paired-wise manner.  
  When the data under evaluation was unpaired, we used two-sided Mann-Whitney U test, which is also non-parametric. 
  Unless otherwise stated, measurements were always taken from distinct samples.
  Significance was established with a p-value≤0.05.
  Boxplots employed in the different figures consist of a box that represents the data quartiles and whiskers that indicate the extent of data points within 1.5 times the interquartile range. 
  
  
  
   §.§ Animal studies
  
  Mice C57BL/6J were purchased from Charlesriver. 
  For bone marrow studies, mice were analysed at 12–20 weeks of age. For fetal liver analyses, embryos were extracted from the previously euthanized pregnant dams at the developmental stages indicated, through a small abdominal incision and further dissected with forceps under a stereo microscope, where single lobes were separated using a surgical scalpel.
  The lobes were fixed in 2% paraformaldehyde (diluted in PBS) (6h, 4°C), washed twice in ice-cold PBS for five minutes and subsequently dehydrated in 30% sucrose in PBS (24–48h, 4°C). The liver lobes were placed into base molds (15×15×5 mm) and completely covered in OCT medium, snap-frozen using liquid nitrogen and stored at -80°C until use. Slices were stained following the same protocols described in <cit.> for marrow tissues.
  Experimental animals were not randomized and experiments were performed in a non-blinded fashion. 
  Animals were maintained at the animal facility of the University Hospital Zürich and treated in accordance with the guidelines of the Swiss Federal Veterinary Office. 
  Experiments and procedures were approved by the Veterinäramt des Kantons Zürich, Switzerland.
  
  
  
  § DATA AVAILABILITY
  
  The labeled dataset employed for training and evaluation of the models described is included as a single HDF5 file within a CodeOcean capsule in <https://codeocean.com/capsule/8424915/tree/v1><cit.>. 
  
  
  
  § CODE AVAILABILITY
  
  The code employed for training the models described in this paper is publicly available on the CodeOcean platform as <https://codeocean.com/capsule/8424915/tree/v1><cit.>. 
  This capsule also includes the trained models employed for the different presented figures. 
  MS-ME is also implemented within MiNTiF, our Fiji plugin for ImageJ for  user-friendly training and deployment of CNNs by non-experts of deep learning: <https://github.com/CAiM-lab/MiNTiF>.
  
  
  
  
  § ACKNOWLEDGMENTS
  
  We thank Takashi Nagasawa for the Cxcl12-GFP reporter mice, Tomomasa Yokomizo for the Hlf-tdTomato mouse strain and Mineo Kurokawa for the Evi1-GFP mouse strain. 
  This work was enabled by funding from the Hasler Foundation (A.G., O.G.).
  We also acknowledge the support from the Swiss National Science Foundation (SNSF) 179116 (O.G.) & 310030_185171 (C.N.-A.), the Swiss Cancer Research Foundation KFS-3986-08-2016 (C.N.-A.), the Clinical Research Priority Program “ImmunoCure” of the University of Zurich to (C.N.-A.).
  We extend our thanks to NVIDIA for their GPU support.
  
  
  
  § AUTHOR CONTRIBUTIONS
  
  A.G. designed and performed the experiments and drafted the manuscript.
  T.P. and O.G. provided key methodological ideas and input.
  A.G., P.M.H., S.I. and U.S. conducted data acquisition, manual labeling of the images, and helped interpret bone marrow analysis in the biological context.
  T.P., O.G., and C.N.-A. critically revised the text. 
  A.G., C.N.-A., and O.G. conceived this project.
  This work was jointly directed by C.N.-A. for the biological and data acquisition aspects and O.G. for computer scientific and methodological aspects. 
  All authors discussed the results and provided feedback to the writing. 
  
  
  10
  url<#>1urlprefixURL 
  
  
  lecun2015deep
  authorLeCun, Y., authorBengio, Y. &
    authorHinton, G.
  titleDeep learning.
  journalNature volume521,
    pages436–444 (year2015).
  
  litjens2017survey
  authorLitjens, G. et al.
  titleA survey on deep learning in medical image analysis.
  journalMedical image analysis
    volume42, pages60–88
    (year2017).
  
  ronneberger2015u
  authorRonneberger, O., authorFischer, P. &
    authorBrox, T.
  titleU-net: Convolutional networks for biomedical image
    segmentation.
  In booktitleInternational Conference on Medical
    image computing and computer-assisted intervention,
    pages234–241 (organizationSpringer,
    year2015).
  
  zhang2018overview
  authorZhang, Y. & authorYang, Q.
  titleAn overview of multi-task learning.
  journalNational Science Review
    volume5, pages30–43 (year2018).
  
  wang2018deep
  authorWang, M. & authorDeng, W.
  titleDeep visual domain adaptation: A survey.
  journalNeurocomputing
    volume312, pages135–153
    (year2018).
  
  guo2020revealing
  authorGuo, S.-M. et al.
  titleRevealing architectural order with quantitative
    label-free imaging and deep learning.
  journalElife volume9,
    pagese55502 (year2020).
  
  ounkomol2018label
  authorOunkomol, C., authorSeshamani, S.,
    authorMaleckar, M. M., authorCollman, F. &
    authorJohnson, G. R.
  titleLabel-free prediction of three-dimensional
    fluorescence images from transmitted-light microscopy.
  journalNature methods
    volume15, pages917–920
    (year2018).
  
  christiansen2018silico
  authorChristiansen, E. M. et al.
  titleIn silico labeling: predicting fluorescent labels in
    unlabeled images.
  journalCell volume173,
    pages792–803 (year2018).
  
  ouyang2018deep
  authorOuyang, W., authorAristov, A.,
    authorLelek, M., authorHao, X. &
    authorZimmer, C.
  titleDeep learning massively accelerates super-resolution
    localization microscopy.
  journalNature biotechnology
    volume36, pages460–468
    (year2018).
  
  li2014deep
  authorLi, R. et al.
  titleDeep learning based imaging data completion for
    improved brain disease diagnosis.
  In booktitleInternational Conference on Medical
    Image Computing and Computer-Assisted Intervention,
    pages305–312 (organizationSpringer,
    year2014).
  
  iglesias2013synthesizing
  authorIglesias, J. E. et al.
  titleIs synthesizing mri contrast useful for
    inter-modality analysis?
  In booktitleInternational Conference on Medical
    Image Computing and Computer-Assisted Intervention,
    pages631–638 (organizationSpringer,
    year2013).
  
  chartsias2017multimodal
  authorChartsias, A., authorJoyce, T.,
    authorGiuffrida, M. V. & authorTsaftaris, S. A.
  titleMultimodal mr synthesis via modality-invariant latent
    representation.
  journalIEEE transactions on medical imaging
    volume37, pages803–814
    (year2017).
  
  lee2019collagan
  authorLee, D., authorKim, J., authorMoon,
    W.-J. & authorYe, J. C.
  titleCollagan: Collaborative gan for missing image data
    imputation.
  In booktitleProceedings of the IEEE Conference on
    Computer Vision and Pattern Recognition, pages2487–2496
    (year2019).
  
  havaei2016hemis
  authorHavaei, M., authorGuizard, N.,
    authorChapados, N. & authorBengio, Y.
  titleHemis: Hetero-modal image segmentation.
  In booktitleInternational Conference on Medical
    Image Computing and Computer-Assisted Intervention,
    pages469–477 (organizationSpringer,
    year2016).
  
  dorent2019hetero
  authorDorent, R., authorJoutard, S.,
    authorModat, M., authorOurselin, S. &
    authorVercauteren, T.
  titleHetero-modal variational encoder-decoder for joint
    modality completion and segmentation.
  In booktitleInternational Conference on Medical
    Image Computing and Computer-Assisted Intervention, pages74–82
    (organizationSpringer, year2019).
  
  varsavsky2018pimms
  authorVarsavsky, T., authorEaton-Rosen, Z.,
    authorSudre, C. H., authorNachev, P. &
    authorCardoso, M. J.
  titlePimms: permutation invariant multi-modal
    segmentation.
  In booktitleDeep Learning in Medical Image Analysis
    and Multimodal Learning for Clinical Decision Support,
    pages201–209 (publisherSpringer,
    year2018).
  
  jaderberg2015spatial
  authorJaderberg, M., authorSimonyan, K.,
    authorZisserman, A. et al.
  titleSpatial transformer networks.
  In booktitleAdvances in neural information
    processing systems, pages2017–2025 (year2015).
  
  wang2017residual
  authorWang, F. et al.
  titleResidual attention network for image classification.
  In booktitleProceedings of the IEEE Conference on
    Computer Vision and Pattern Recognition, pages3156–3164
    (year2017).
  
  hu2018squeeze
  authorHu, J., authorShen, L. & authorSun,
    G.
  titleSqueeze-and-excitation networks.
  In booktitleProceedings of the IEEE conference on
    computer vision and pattern recognition, pages7132–7141
    (year2018).
  
  roy2020squeeze
  authorRoy, A. G., authorSiddiqui, S.,
    authorPölsterl, S., authorNavab, N. &
    authorWachinger, C.
  title‘squeeze & excite’guided few-shot segmentation
    of volumetric images.
  journalMedical image analysis
    volume59, pages101587
    (year2020).
  
  roy2018recalibrating
  authorRoy, A. G., authorNavab, N. &
    authorWachinger, C.
  titleRecalibrating fully convolutional networks with
    spatial and channel “squeeze and excitation” blocks.
  journalIEEE transactions on medical imaging
    volume38, pages540–549
    (year2018).
  
  rickmann2019project
  authorRickmann, A.-M., authorRoy, A. G.,
    authorSarasua, I., authorNavab, N. &
    authorWachinger, C.
  title‘project & excite’modules for segmentation of
    volumetric medical scans.
  In booktitleInternational Conference on Medical
    Image Computing and Computer-Assisted Intervention, pages39–47
    (organizationSpringer, year2019).
  
  wang2019towards
  authorWang, X., authorCai, Z., authorGao,
    D. & authorVasconcelos, N.
  titleTowards universal object detection by domain
    attention.
  In booktitleProceedings of the IEEE Conference on
    Computer Vision and Pattern Recognition, pages7289–7298
    (year2019).
  
  gomariz2018quantitative
  authorGomariz, A. et al.
  titleQuantitative spatial analysis of
    haematopoiesis-regulating stromal cells in the bone marrow microenvironment
    by 3d microscopy.
  journalNature communications
    volume9, pages1–15 (year2018).
  
  gomariz2019imaging
  authorGomariz, A., authorIsringhausen, S.,
    authorHelbling, P. M. & authorNombela-Arrieta, C.
  titleImaging and spatial analysis of hematopoietic stem
    cell niches.
  journalAnnals of the New York Academy of Sciences
    (year2019).
  
  acar2015deep
  authorAcar, M. et al.
  titleDeep imaging of bone marrow shows non-dividing stem
    cells are mainly perisinusoidal.
  journalNature volume526,
    pages126–130 (year2015).
  
  coutu2017three
  authorCoutu, D. L., authorKokkaliaris, K. D.,
    authorKunz, L. & authorSchroeder, T.
  titleThree-dimensional map of nonhematopoietic bone and
    bone-marrow cells and molecules.
  journalNature biotechnology
    volume35, pages1202 (year2017).
  
  coutu2018multicolor
  authorCoutu, D. L., authorKokkaliaris, K. D.,
    authorKunz, L. & authorSchroeder, T.
  titleMulticolor quantitative confocal imaging cytometry.
  journalNature methods
    volume15, pages39 (year2018).
  
  christodoulou2020live
  authorChristodoulou, C. et al.
  titleLive-animal imaging of native haematopoietic stem and
    progenitor cells.
  journalNature pages1–6
    (year2020).
  
  srivastava2014dropout
  authorSrivastava, N., authorHinton, G.,
    authorKrizhevsky, A., authorSutskever, I. &
    authorSalakhutdinov, R.
  titleDropout: a simple way to prevent neural networks from
    overfitting.
  journalThe journal of machine learning research
    volume15, pages1929–1958
    (year2014).
  
  Kunisaki_2013
  authorKunisaki, Y. et al.
  titleArteriolar niches maintain haematopoietic stem cell
    quiescence.
  journalNature volume502,
    pages637–643 (year2013).
  
  gal2015cnn
  authorGal, Y. & authorGhahramani, Z.
  titleBayesian convolutional neural networks with bernoulli
    approximate variational inference.
  In booktitleICLR workshop track
    (year2016).
  
  kendall2017uncertainties
  authorKendall, A. & authorGal, Y.
  titleWhat uncertainties do we need in bayesian deep
    learning for computer vision?
  In booktitleAdvances in Neural Information
    Processing Systems 30, pages5574–5584 (year2017).
  
  transfusion_bengion19
  authorRaghu, M., authorZhang, C.,
    authorKleinberg, J. & authorBengio, S.
  titleTransfusion: Understanding transfer learning for
    medical imaging.
  In booktitleAdvances in Neural Information
    Processing Systems, vol. volume32, pages3347–3357
    (year2019).
  
  AlexNet
  authorKrizhevsky, A., authorSutskever, I. &
    authorHinton, G. E.
  titleImagenet classification with deep convolutional
    neural networks.
  In editorPereira, F., editorBurges, C. J. C.,
    editorBottou, L. & editorWeinberger, K. Q. (eds.)
    booktitleAdvances in Neural Information Processing Systems
    25, pages1097–1105 (year2012).
  
  He2016
  authorHe, K., authorZhang, X., authorRen,
    S. & authorSun, J.
  titleDeep residual learning for image recognition.
  In booktitle2016 IEEE Conference on Computer
    Vision and Pattern Recognition (CVPR) (publisherIEEE,
    year2016).
  
  Xie2018
  authorXie, Y. et al.
  titleEfficient and robust cell detection: A structured
    regression approach.
  journalMedical Image Analysis
    volume44, pages245–254
    (year2018).
  
  Ren2017
  authorRen, S., authorHe, K., authorGirshick,
    R. & authorSun, J.
  titleFaster r-CNN: Towards real-time object detection
    with region proposal networks.
  journalIEEE Transactions on Pattern Analysis and
    Machine Intelligence volume39,
    pages1137–1149 (year2017).
  
  He2017
  authorHe, K., authorGkioxari, G.,
    authorDollar, P. & authorGirshick, R.
  titleMask r-CNN.
  In booktitle2017 IEEE International Conference on
    Computer Vision (ICCV) (publisherIEEE,
    year2017).
  
  Yang2020
  authorYang, L. et al.
  titleNuSeT: A deep learning tool for reliably separating
    and analyzing crowded cells.
  journalPLOS Computational Biology
    volume16, pagese1008193
    (year2020).
  
  falk2019u
  authorFalk, T. et al.
  titleU-net: deep learning for cell counting, detection,
    and morphometry.
  journalNature methods
    volume16, pages67–70
    (year2019).
  
  tensorflow2015-whitepaper
  authorAbadi, M. et al.
  titleTensorFlow: Large-scale machine learning on
    heterogeneous systems (year2015).
  <https://www.tensorflow.org/>.
  noteSoftware available from tensorflow.org.
  
  kingma2014adam
  authorKingma, D. P. & authorBa, J.
  titleAdam: A method for stochastic optimization.
  journalarXiv preprint arXiv:1412.6980
    (year2014).
  
  li2019overfitting
  authorLi, Z., authorKamnitsas, K. &
    authorGlocker, B.
  titleOverfitting of neural nets under class imbalance:
    Analysis and improvements for segmentation.
  In booktitleInternational Conference on Medical
    Image Computing and Computer-Assisted Intervention,
    pages402–410 (organizationSpringer,
    year2019).
  
  gomez2019deepimagej
  authorGómez-de Mariscal, E. et al.
  titleDeepimagej: A user-friendly plugin to run deep
    learning models in imagej.
  journalbioRxiv pages799270
    (year2019).
    
  MSME_CO
  authorGomariz, A. et al.
  titleMarker Sampling and Excite: Deep Learning with Heterogeneous Marker Combinations in Fluorescence Microscopy.
  journalCodeOcean (year2021)
  <https://codeocean.com/capsule/8424915/tree/v1>.
  
  
  
  
  
  
  
  § SUPPLEMENTARY MATERIALS
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      Example images for the qualitative assessment of the segmentation of bone marrow images employed for the quantification of vasculature. Input images contain different combinations of markers shown as an overlay of different colors. The white rectangle within the input images represents the size of the output image when processing with a neural network. White arrows depict inference with the MS-ME model. Different colors are employed in the output images to show the different predicted classes. Since the tissue class overlaps with the other two, its contour is used instead for visual purposes. 
      
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      Comparison of bone marrow volume ratio occupied by sinusoids in both diaphysis (DIA) and metaphysis (META) when segmenting them with the morphological image processing (MIP) algorithm previously proposed (n=12 for both DIA and META) and with our MS-ME method proposed here (n=61 for DIA, n=24 for META).
      
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      Illustration of the image tiling pipeline employed to create suitable patches for CNNs. 
      An example of a slice within the 3D image frame is shown in the upper part using Imaris (Bitplane AG). 
      That slice is decomposed in patches as illustrated in the lower part. Each output patch (red) is smaller than their corresponding input patch (cyan) due to the convolutional operations in CNNs. We position the output patches next to each other without overlap in order to avoid padding artifacts in the application of CNNs. Instead, zero padding is only applied along the borders of the whole slice (area with white stripes). When an overlap between output patches cannot be avoided to fill the slice (area with green stripes), the average of the different patches in that region is used. 
      
      
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      
      Effect of bias term on ME module. F1-score of MS-ME model where the bias terms for all ME modules have been removed, relative to the proposed MS-ME model with bias across all marker combinations and cross-validation steps. 
      
      
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      
      Effect of marker dropout rate r_drop in MS. 
      F1-score of MS models with different r_drop evaluated on the sinusoids (left) and arteries (right) relative to the proposed MS with r_drop=0.5.
      (a) Evaluation for all 31 possible marker combinations (n=124).
      Whereas some r_drop≠ 0.5 produce a slightly superior F1-score for sinusoids, r_drop=0.5 is the best option for arteries and overall. 
      (b) Median relative F1-score for models evaluated on combinations of specific numbers of markers, each represented by a different color for the different r_drop (n=K!/(K-k)!k!, where K is the number of markers available, and k the number of markers considered for each evaluation).
      Smaller r_drop are shown to be beneficial for combinations of more markers, and vice-versa.
      However, this trend becomes noisier for r_drop>0.5, as illustrated with the gray dashed line. 
      This effect can be due to the decrease in markers observed over time, although it is a question worth of further investigation in future work. 
      
      
  
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      
      Analysis of attention parameters in ME modules for sinusoids (left) and arteries (right).
      We estimate recalibration strength by calculating cosine distances between the ME attention subnetwork outputs obtained for each of the possible input marker combinations. 
      Results are represented as the mean of all such pair-wise distances between all possible marker combinations, at a given layer where ME is placed, with the bars depicting the standard deviations of these distances. 
      Using the colored network layers shown in Fig.  <ref>b, Encoder layers correspond to the network layers in blue, Decoder to the layers in green, and Bottleneck to the yellow layers.
      The numbers next to each layer indicate the resolution level, where 0 corresponds to the highest (original resolution) and 3 to the lowest (i.e., right before and after the bottleneck, for the encoder and decoder, respectively). 
      It can be seen in this representation that recalibration strength is higher in layers with higher resolution, especially near the input of the network.
      This observation may indicate that high resolution layers of the network focus on effectively combining features from available markers, and in this way create shared abstract features that are common across markers for subsequent processing in lower resolution layers.
      
      
  
  
  [hbt!]
      
      < g r a p h i c s >
  
      
      Evolution of the weighted cross-entropy loss (top) and the F1-score (bottom) across epochs for the training (blue) and validation (orange) sets with our proposed MS-ME model, both for models trained for segmentation of sinusoids (left) and arteries (right). 
      The red dashed line marks the epoch at which we choose the model, based on the highest validation F1-score.
      
      
  
  
  
  
  [hbt!]
      0.8
      
      sample id (s^i)     # patches (J^i)     sinusoids (%)     arteries (%) 
  
      1     54     9.71     0.29 
  
      2     12     7.87     0.99 
  
      3     30     11.97     0.21 
  
      4     12     11.20     0.98 
  
      5     12     10.93     0.58 
  
      6     48     13.60     0.24 
  
      7     54     11.23     0.24 
  
      8     8     15.19     0.27 
  
      TOTAL     230     11.41     0.34 
  
      
      
      Summary of the annotated bone marrow dataset employed for training and evaluation of the proposed methods. There are 8 different samples composed of different annotated patches. The annotated classes are either sinusoids or arteries. The remaining pixels are considered background. There is a notable class imbalance, with most of the pixels being background and arteries being a minority class.  
      
  
  
  [hbt!]
      0.8
      
      sample id (s^i)     # patches (J^i)     vessels (%)
  
      0     4     20.33 
  
      1     8     20.81 
  
      2     4     18.88 
  
      3     4     13.25 
  
      4     4     20.54 
  
      5     8     21.11 
  
      6     4     25.43 
  
      7     4     15.41 
  
      8     4     27.20 
  
      9     4     10.04 
  
      10     8     12.78 
  
      11     4     25.04 
  
      12     4     17.90 
  
      13     4     21.67 
  
      14     4     31.79 
  
      15     4     5.32 
  
      16     8     6.97 
  
      17     4     23.02 
  
      18     9     18.69 
  
      19     4     18.39 
  
      20     4     20.79 
  
      21     4     27.87 
  
      22     4     22.87 
  
      23     8     12.27 
  
      24     4     1.41 
  
      25     17     16.97 
  
      TOTAL     142     17.72 
  
      
      
      
      Summary of the annotated dataset employed for the segmentation of fetal liver vasculature, where the only class is named vessels and the remaining pixels are considered background. 
      
      
  
  
  [hbt!]
      0.8
      
      Models     Training time 
  (ms/batch)     Inference time 
  (ms/batch)     GPU memory 
  footprint (GB)     # parameters
  (x10^6) 
  
      MZ     0.28     0.20     0.95     7.76 
  
      HeMIS     1.07     0.17     3.8     0.35 
  
      MS     0.32     0.19     0.95     7.76 
  
      HeMIS-MS     1.00     0.20     3.9     0.35 
  
      MS-DR     0.34     0.17     0.95     7.76 
  
      MS-VR     0.32     0.18     0.95     7.76 
  
      MS-SE     0.35     0.15     1.12     8.2 
  
      MS-ME     0.34     0.2     1.16     7.81 
  
      MS+     0.25     0.18     1.05     9.35 
  
      UB     7.57     0.2     0.95     240.56 
  
      
      
      
      Computational complexity of the different models employed. 
      Time is calculated as the median across batches. 
      GPU memory footprint refers to the maximum usage recorded during training. 
      Note that UB is an ensemble of 31 models, and hence the increase in parameters only affects the training time.
      A single model is employed at inference for the pertinent marker combination, making the inference time equivalent to that of a single model. 
      
      
  
  
  