# paper.pdf ## Metadata - PDFFormatVersion=1.7 - IsLinearized=true - IsAcroFormPresent=false - IsXFAPresent=false - IsCollectionPresent=false - IsSignaturesPresent=false - Author=Sanjay Patel - CreationDate=D:20260625150117-04'00' - Creator=Adobe Acrobat (64-bit) 26.1.21662 - ModDate=D:20260625150217-04'00' - Producer=Adobe Acrobat (64-bit) 26.1.21662 - xmp:modifydate=2026-06-25T15:02:17-04:00 - xmp:createdate=2026-06-25T15:01:17-04:00 - xmp:metadatadate=2026-06-25T15:02:17-04:00 - xmp:creatortool=Adobe Acrobat (64-bit) 26.1.21662 - dc:format=application/pdf - dc:creator=Sanjay Patel - xmpmm:documentid=uuid:df3900e8-82d8-48aa-ab38-1a65fffc2ea7 - xmpmm:instanceid=uuid:2fc83b3a-7c9a-44c8-8cd9-b074dd3cc48a - pdf:producer=Adobe Acrobat (64-bit) 26.1.21662 ## Contents ### Page 1 SUPPLEMENTARY MATERIALS Spatial remodeling of bone marrow architecture defines tissue-state signatures of disease activity and therapeutic response in myelodysplastic neoplasms Ryan Nachman,1 Aleksandra Kopacz,1 Caitlin Unkenholz,1,4 Jian Chai, 2 Arvin Ruiz,3 Itzel Valencia,4 Jeanne Jiang,4 Fabio Socciarelli,4 Jiwoon Park,5 Christopher Mason,(5-8) Ling Zhang,9 David Sallman,10 Gail J. Roboz,11 Pinkal Desai,11 Justin Kaner,11 Joshua Fein,11 Monica Guzman,11 Neal Lindeman,12 Amy Chadburn,2 Madhu Ouseph,2,12 Paul Simonson,2 Julia Geyer,2 Giorgio Inghirami,2,3 Shahin Rafii,1 David Redmond,1* Sanjay Patel2,4* 1Division of Regenerative Medicine, Hartman Institute for Therapeutic Organ Regeneration, Ansary Stem Cell Institute, Department of Medicine, Weill Cornell Medicine, New York, NY, USA 2Division of Hematopathology, Department of Pathology and Laboratory Medicine, Weill Cornell Medicine/NewYork-Presbyterian Hospital, New York, NY, USA 3Immunopathology Laboratory, Weill Cornell Medical College/NewYork-Presbyterian Hospital, New York, NY, USA 4Multiparametric In Situ Imaging (MISI) Laboratory, Department of Pathology and Laboratory Medicine, Weill Cornell Medicine, New York, NY, USA 5Department of Physiology and Biophysics, Weill Cornell Medicine, New York, NY, USA 6The Feil Family Brain and Mind Research Institute, Weill Cornell Medicine, New York, NY 7WorldQuant Initiative for Quantitative Prediction, Weill Cornell Medicine, New York, NY 8The HRH Prince Alwaleed Bin Talal Bin Abdulaziz Alsaud Institute for Computational Biomedicine, Weill Cornell Medicine, New York, NY 9Department of Hematopathology and Laboratory Medicine, Moffitt Cancer Center, Tampa, FL, USA 10Department of Malignant Hematology, Moffitt Cancer Center, Tampa, FL, USA 11Division of Hematology and Medical Oncology, Department of Medicine, Weill Cornell Medicine/NewYork-Presbyterian Hospital, New York, USA 12Division of Molecular Genetic Pathology, Department of Pathology and Laboratory Medicine, Weill Cornell Medicine/NewYork-Presbyterian Hospital, New York *Correspondence: Sanjay Patel, MD () & David Redmond, PhD () | ### Page 2 METHODS Clinicopathologic and Laboratory Data Electronic medical records were reviewed for all patients to extract clinical data concurrent with the evaluated biopsy, including time of collection, body mass index (BMI), complete blood count parameters (hemoglobin, mean corpuscular volume, absolute neutrophil count, platelet count), and marrow blast percentage. Diagnostic hematopathology reports were re-reviewed by board-certified hematopathologists (SSP, JTG) to confirm diagnostic classification (NBM, CCUS, MDS) and to extract aspirate differential counts, including myeloid-to-erythroid (M:E) ratio and lineage-specific proportions. International Prognostic Scoring System–Molecular (IPSS-M) scores were calculated for all treatment- naïve MDS samples. Cytogenetic and molecular findings were incorporated into diagnostic categorization as described in the main Methods. Clinical and laboratory parameters, flow cytometry and NGS data, as well as longitudinal clinical/treatment response data are presented in Supplemental Data Tables 1-5. Study Cohort and Clinical Definitions The cohort comprised 49 unique patients with myelodysplastic neoplasms (MDS; n=36) or clonal cytopenia(s) of undetermined significance (n=13), contributing 77 trephine biopsy samples. Twenty-nine serial biopsies from 14 MDS patients (range 2–7 per patient) were obtained after the diagnostic biopsy. Twenty-one normal bone marrow (NBM) controls from a previously published cohort were included. A tissue microarray (TMA) incorporating 7 additional TP53-mutated MDS cases and 5 NBMs was included in analyses corresponding to main text Figures 2–3. Clinical definitions were prespecified and adjudicated prior to development of the MDS Microarchitectural Perturbation Score (MDS-MAPS), and therefore without knowledge of MAPS results. Diagnostic (treatment-naïve) samples were those establishing diagnosis prior to therapy. Serial samples were those obtained after diagnosis. Disease-modifying therapies primarily included hypomethylating agents (azacytidine or decitabine) with or without venetoclax (Supplemental Data Tables 1 and 3). ### Page 3 Remission versus Active Disease Classification of Individual Biopsies Remission versus active disease status for serial biopsies was determined through integration of 2023 International Working Group (IWG) response criteria,1 contemporaneous pathology reports, and concurrent cytogenetic and molecular findings. Biopsies were classified as remission if they met IWG criteria for complete remission (CR) or CR with incomplete hematologic recovery (CRi), including marrow blasts <5%, and were interpreted morphologically as remission without persistent dysplasia or excess blasts. Cases were classified as active disease if driver mutations and/or clonal cytogenetic abnormalities consistent with residual MDS persisted, even when blasts were <5%. Molecular remission was defined as absence of detectable mutations to a sensitivity of 5% variant allele frequency (VAF) using clinically validated non–minimal residual disease targeted sequencing. Blast progression was defined as any increase in marrow blast percentage relative to baseline. Clinical laboratory data were abstracted from medical records. Clinical Laboratory Characterization All diagnostic samples underwent standard clinical evaluation including histomorphologic review, multiparameter flow cytometry (MFC), cytogenetic analysis, and targeted next-generation sequencing (NGS). MFC was performed using a 10-color EuroFlow2 panel with acquisition of ≥200,000 events per sample and analyzed using BD FACSDiva. Blast populations were defined as CD34+/CD117+ cells. Flow-derived proportions were used for orthogonal validation of multiplex immunofluorescence–derived cell types. Conventional cytogenetic analysis was performed on 20 metaphases per case and reported according to International System for Human Cytogenetic Nomenclature guidelines. FISH was performed when clinically indicated using probes targeting recurrent MDS-associated abnormalities. Targeted NGS was performed using either a 45-gene panel (analytical sensitivity ~2% VAF for SNVs, 1% for indels) or the TruSight Oncology 500 assay (~3% VAF sensitivity) on Illumina platforms, with analysis using clinically validated pipelines. Mutation status was derived from diagnostic reports. ### Page 4 Multiplex Immunofluorescence Tissue Staining Four-micron-thick Bouin-fixed, decalcified, paraffin-embedded trephine sections were prepared and stained using the Opal multiplex immunofluorescence system (Akoya Biosciences) on a Bond RX automated stainer (Leica Biosystems).3 Slides were deparaffinized and subjected to EDTA-based antigen retrieval (Leica ER2, 20 minutes). Cyclical horseradish peroxidase–mediated tyramide signal amplification was performed, with antibody stripping between cycles using citrate-based retrieval (Leica ER1). The primary antibody panel (applied sequentially) included: - 1. CD71 (Opal 480; clone 10F11; 1:80) - 2. CD61 (Opal 520; clone 2F2; ready-to-use) - 3. CD117 (Opal 570; clone D3W6Y; 1:100) - 4. CD38 (Opal 620; clone 38C03; 1:50) - 5. CD34 (Opal 690; clone QBEND/10; 1:100) - 6. CD15 (Opal 780; clone MMA; 1:50) For tissue microarray (TMA) analyses, CD38 was replaced with p53 (clone DO-7) [examples below]. | Representative single channel and merge immunofluorescence images from validation of p53-containingpanel applied to TMA. Scale bars shown are 100μm (tonsil, 11.2x magnification), and both 100μm (top row,11.2x magnification) and 20μm (bottom row, 48.2x magnification) for all other case examples. | | - | | | ### Page 5 | | | - | | ![A screenshot of a grid of squares AI-generated content may be incorrect.]() | For both the stained whole trephine biopsy sections and the TMA slides, whole-slide images (WSIs) were acquired at 20× magnification using the Vectra Polaris / PhenoImager platform (Akoya Biosciences). Spectral ### Page 6 unmixing was performed in InForm (v2.4.8), tiles were fused in HALO (v3.3.2541.231), and multispectral image stacks were exported as multi-layer TIFF files for computational analysis. Image Processing and Segmentation Pipeline Image analysis was performed in Python (v3.11) using open-source libraries. Whole-cell and nuclear segmentation was performed using Mesmer (DeepCell v0.12.9),4 a convolutional neural network–based segmentation framework optimized for multiplex immunofluorescence images. The Mesmer architecture consists of a feature pyramid network coupled to a ResNet50 backbone, with dual semantic heads that predict inner-distance transforms and pixel-wise classification maps to delineate individual cells. The model was pretrained on an augmented TissueNet dataset and subsequently fine-tuned using manually annotated 256 × 256 pixel bone marrow image patches derived from this cohort. Fine-tuning employed the default Mesmer parameter configuration with a four-channel semantic head output (inner distance, outer distance, foreground–background classification, and pixel-wise transformations). Segmentation employed a watershed-based postprocessing pipeline with an interior probability threshold of 0.3 and a centroid threshold of 0.1. A deep watershed transform was applied to extract cell masks, using parameters maxima_threshold = 1, interior_threshold = 0.3, and interior_smoothing = 1. The minimum nuclear and whole-cell area was set to 15 pixels, with no imposed maximum area constraint. Edge-touching cells were retained and not excluded from downstream analyses. Whole-slide segmentation was performed by tiling images into overlapping 256 × 256 pixel patches. Segmentation masks were subsequently merged across tile boundaries using a Dask-image–based patch stitching algorithm to generate continuous whole-slide masks. Segmentation accuracy was validated against manually annotated masks, with performance metrics reported in Supplementary Table 5. ### Page 7 Cell Phenotyping and Classification Single-cell membrane marker classification was performed using a ResNet18 convolutional neural network initialized with ImageNet pretrained weights. The neural network was fine-tuned on a dataset of expert hematopathologist annotations, including examples of all 6 membrane markers, which is detailed in Supplementary Table 5. For each segmented cell, the classifier received a three-channel cropped image centered on the cell of interest, consisting of the binary whole-cell segmentation mask, the DAPI nuclear channel, and the individual membrane marker channel being evaluated. Marker positivity classifications were validated against expert hematopathologist annotations on a held-out set of image examples, with performance metrics reported in Supplemental Table 5. The classifier was applied separately to each membrane marker. . ### Page 8 Each classifier generated a continuous probability score (range 0–1) representing the likelihood of marker positivity. Marker intensities were normalized using a slide-wise normalization approach prior to classification. Positivity thresholds were determined using manual gating aligned with expert hematopathologist review. Unless otherwise specified, cells were considered positive for a given marker at a probability threshold of >0.1. For CD34 and CD61, more stringent thresholds were applied (>0.25 probability) to minimize false- positive assignments in structurally complex regions. Following marker-level classification, cell identities were assigned using a vector-based approach. For each cell, the predicted marker probability vector was compared to predefined lineage signature vectors using cosine similarity, and the phenotype corresponding to the closest matching signature was assigned as the final cell identity. The lineage signature vectors were constructed based on canonical immunophenotypic profiles established in clinical hematopathology, with each vector specifying a binary pattern of expected marker positivity and negativity for a given cell type (e.g., myeloblast [CD34⁺/CD117⁺]: [1,1,0,0,0,0]; proerythroblast [CD117⁺/CD71⁺]: [0,1,0,1,0,0]). Cosine similarity was selected for vector matching because it emphasizes the pattern of marker expression rather than absolute classifier confidence magnitude, providing robustness to inter- slide staining intensity variation inherent in multiplex immunofluorescence imaging. | | | - | ### Page 9 | | | - | Discrimination of CD34+ hematopoietic cells from CD34+ endothelial cells Endothelial cells and CD34⁺ hematopoietic cells are distinguished in our pipeline through two complementary approaches. First, sinusoidal vasculature is identified using a dedicated structural segmentation mask that detects vessel-like structures based on morphologic features (elongated, branching architecture) distinct from round hematopoietic cells. Second, CD34 mean fluorescence intensity (MFI) is used as a discriminating feature: endothelial cells exhibit higher CD34 MFI than hematopoietic progenitors. Specifically, CD34⁺ nuclei with >50% spatial overlap with the vascular segmentation mask are classified as endothelial-associated and excluded from hematopoietic cell quantification; remaining CD34⁺ cells are available for hematopoietic cell type classification. This dual criterion — morphology-based structural masking combined with intensity-based and spatial overlap thresholding — minimizes misclassification between the two CD34⁺ populations. p53 Classification (TMA) ### Page 10 To enable in situ identification of putative TP53-mutated cells, a separate convolutional neural network classifier was developed to categorize nuclear p53 protein expression levels. A ResNet18 architecture initialized with ImageNet-pretrained weights was trained to classify individual nuclei into three expression categories: no expression, low expression, and moderate/high expression. Cells classified as moderate/high were operationally designated as p53-positive for downstream spatial analyses. A total of 76 manually annotated single-cell nuclear p53 images were used for model development. Annotations were performed by expert review and stratified into the three expression categories described above. Images were randomly partitioned into training (80%) and validation (20%) sets, with 25 annotated single-cell images used for validation. On the held-out validation set, the classifier achieved an overall accuracy of 0.95, precision of 1.00, recall of 0.86, and an F1 score of 0.92 for identification of moderate/high p53 expression. These performance metrics supported the use of the model for in situ detection of p53-overexpressing cells within multiplex immunofluorescence datasets. A comprehensive summary of model architecture, preprocessing steps, training parameters, and validation procedures is provided in Supplementary Table 5. Structural Masks Segmentation CNN-based segmentation was used to generate binary masks for endothelium (CD34+ elongated structures), bone trabeculae, and adipose tissue, as described in Sarachakov et al. Cell-to-structure distances were calculated ### Page 11 by subtracting estimated cell radius from centroid-to-boundary distance using OpenCV functions (findContours, pointPolygonTest). Morphologic Feature Extraction Morphologic features were computed directly from the whole-cell and nuclear segmentation masks generated during image processing. For each segmented cell, total cell area was calculated by converting the pixel area of the whole-cell mask into square microns (μm²) using the known image resolution. Nuclear area was similarly computed from the corresponding nuclear mask. Cell shape was quantified using eccentricity, defined as the ratio of the major axis to the minor axis of an ellipse fitted to the segmentation mask, providing a measure of deviation from circularity. For megakaryocytes, multinucleation was quantified by counting the number of discrete nuclear mask objects contained within the corresponding whole-cell boundary. All morphologic computations, including area measurements, ellipse fitting, and nuclear enumeration, were performed using the OpenCV Python library. Spatial Statistics Implementation Permutation Testing Spatial enrichment analyses were performed using a permutation-based framework to assess whether specific cell phenotypes were preferentially localized relative to structural elements. For each whole-slide image, cell phenotype labels were randomly permuted 1,000 times while preserving the original spatial coordinates, thereby generating a null distribution for enrichment statistics. Distances between cells and structural elements were calculated using the Euclidean distance metric in two-dimensional space. Enrichment was evaluated within predefined radii of 15 μm for vascular proximity and 30 μm for trabecular proximity, selected based on the spatial scale of these structural components within bone marrow tissue. ### Page 12 P-values were calculated as the proportion of permuted datasets demonstrating equal or greater enrichment than observed. No Benjamini–Hochberg false discovery rate correction was applied to permutation- derived p-values in these analyses. Ripley’s K (Rolling Signal) Spatial clustering of erythroid normoblasts was quantified using a rolling-window implementation of Ripley’s K statistic.5 Whole-slide images were subdivided into overlapping windows measuring 256 × 256 μm, with 50% overlap between adjacent windows to ensure continuous spatial coverage while preserving local resolution. ### Page 13 Within each window, Ripley’s K was computed across a range of radii spanning 8.5 to 12.5 μm in 0.5 μm increments, capturing clustering behavior at distances relevant to erythroid island organization. A border edge-correction method was applied to account for truncated neighborhoods at window boundaries. Whole-slide K values were calculated as the mean of the window-level K estimates, providing a summary measure of erythroid clustering across the entire tissue section. HDBSCAN Erythroid clustering used HDBSCAN, with min_cluster_size = 5, min_samples = 3, cluster detection method = eom, and Euclidean values as the distance metric. Cluster density was computed using convex hull area. Cluster density was computed using convex hull area. Atypical Localization of Immature Precursors (ALIP) Analysis Atypical localization of immature precursors (ALIP) clusters were defined as ≥2 myeloblasts and/or promyelocytes within 12 μm of one another and located >150 μm from trabeculae and >40 μm from vasculature. MDS-MAPS Detailed Implementation Patch Definition WSIs were subdivided into non-overlapping 256 × 256 μm patches. Exclusion criteria included patches with <100 cells or >95% non-hematopoietic elements. Feature Scaling Features were scaled using z-score normalization derived exclusively from diagnostic samples. Scaling parameters were locked prior to longitudinal analysis. Weight Derivation ### Page 14 Weights were derived using diagnostic samples only and locked prior to remission modeling. For each 256×256 μm patch, 82 features spanning cell-type proportions, cell morphology, spatial proximity metrics, and microenvironmental localization were computed and z-score normalized across the diagnostic cohort. Feature weights were derived in two steps. First, features were ranked by their correlation with disease severity through unsupervised sample stratification, identifying disease-associated perturbations receiving positive weights (e.g., myeloblast proportion, ALIP clusters) and normal bone marrow features receiving negative weights (e.g., erythroid spatial density, fat-to-cell ratio). Second, weights were adjusted through expert curation by a board- certified hematopathologist (SSP) according to established MDS pathological criteria, organized into four biological domains: blast burden, erythroid dysplasia, megakaryocyte dysplasia, and microenvironmental features. Biologically-motivated interaction terms were additionally incorporated as products of z-scored feature pairs reflecting known pathological co-occurrences, including myeloblast–megakaryocyte proportion (multi- lineage dysplasia), myeloblast proportion–trabecular distance (ALIP pattern), and MMC proportion–fat elements (preserved normal hematopoiesis). This amounted to 32 total weighted features and 4 interaction terms (Supplementary Table 6). The patch-level MAPS was then computed as: MAP = Σ (Z(i) × W(i)) + Σ (Z(j) × Z(k) × W(jk)) where Z(i) is the z-score of feature i, W(i) is its assigned weight, and Z(j) × Z(k) × W(jk) represents weighted interaction terms between feature pairs j and k. ### Page 15 CCUS/MDS-LB Discrimination and Remission Classification Modeling To evaluate the ability of MAPS to discriminate CCUS from MDS-LB and separately, remission from active disease, logistic regression models were implemented using leave-one-out cross-validation (LOOCV) for CCUS/MDS-LB discrimination or a leave-one-patient-out cross-validation (LOPO-CV) framework for remission classification. L2-regularized logistic regression was applied with a regularization strength of C = 1.0, using the lbfgs solver. To prevent data leakage, feature standardization was performed independently within each training fold, and the derived scaling parameters were applied to the corresponding held-out patient samples (or individual samples in LOOCV analysis). Statistical significance of incremental discrimination (ΔAUC) was assessed using permutation testing with 1,000 iterations, in which MAPS values were randomly permuted and the LOOCV or LOPO-CV procedures repeated to generate a null distribution. Bootstrap confidence intervals were not calculated. Specifically for remission classification, sensitivity analyses were performed through stratified evaluation separating stem cell transplantation (SCT)–associated remission samples from non-SCT remission samples, as presented in the corresponding figure. ### Page 16 Blast Progression Modeling To evaluate the relative contributions of clinical and spatial features to blast progression, analyses were performed in a subset of 16 patients for whom follow-up data were available, including 8 patients who experienced an increase in marrow blast percentage and 8 who remained stable. Two L2-regularized logistic regression models were constructed with a regularization strength of C = 0.5 to mitigate overfitting in this limited cohort. The clinical model incorporated established clinical and laboratory variables, including IPSS-M score, mean corpuscular volume (MCV), age at diagnosis, hemoglobin level, absolute neutrophil count (ANC), and platelet count. The spatial model included mean MAPS along with the top five multiplex immunofluorescence (MxIF)–derived features ranked by univariate non-parametric AUC. Feature importance was evaluated using SHAP (SHapley Additive exPlanations) values, calculated as the product of each model coefficient and the corresponding standardized feature value. SHAP values were sign-corrected based on coefficient directionality to facilitate interpretation, such that positive values indicated association with blast increase and negative values indicated association with stability. Bootstrap confidence intervals were generated using 1,000 resamples with replacement, refitting the model in each iteration. Model discrimination was quantified using a non-parametric AUC computed as the Mann–Whitney U statistic normalized by group sizes, representing the fraction of correctly ranked stable versus blast-increase patient pairs. Mixed-Effects Modeling To assess whether changes in MAPS were associated with remission status independent of blast burden while accounting for repeated measurements within patients, a linear mixed-effects model was fitted. The model specified MAPS as the continuous outcome and included fixed effects for remission status and blast percentage, with a random intercept for each patient to account for within-patient correlation across serial biopsies. An identity link function was used, as the outcome variable (MAPS) was continuous. Only random intercepts were included; random slopes were not modeled due to the limited number of serial samples per ### Page 17 patient. Statistical significance of fixed effects was assessed using Wald tests. Formal residual diagnostics were not performed. Linear mixed-effects model: MAPS ~ Remission + Blast% + (1|Patient) Public Dataset Re-Analysis To evaluate transcriptional correlates of spatial niche remodeling observed in our imaging cohort, we re- analyzed publicly available bulk and single-cell transcriptomic datasets of CD34⁺ hematopoietic progenitor populations from patients with MDS and healthy controls. Datasets included GSE136816, GSE111085, and EGAS00001007568. In addition, protein-level expression data were examined from the mass cytometry dataset reported by Behbehani et al. For bulk RNA sequencing datasets (GSE136816 and GSE111085), raw gene-level count matrices were obtained when available. Expression values were analyzed in count space and normalized using the trimmed mean of M-values (TMM) method. Differential expression testing was performed using generalized linear modeling with likelihood ratio tests. Multiple hypothesis testing correction was applied using the Benjamini– Hochberg procedure, and genes were considered significantly differentially expressed at a false discovery rate (FDR) < 0.05 (adjusted P < 0.05). For the single-cell RNA sequencing dataset (EGAS00001007568), normalized expression matrices were used as provided by the original study authors. Cell-type annotations were retained from the published metadata. Gene expression comparisons between MDS and healthy CD34⁺ progenitor populations were performed using rank-based or model-based testing as appropriate to the dataset structure, with multiple testing correction applied using the Benjamini–Hochberg method (FDR < 0.05). Mass cytometry data from Behbehani et al. were analyzed at the sample level to quantify relative CXCR4 protein expression across hematopoietic cell populations. Protein abundance values were used as ### Page 18 reported in the original publication, and comparisons between diagnostic categories were performed using nonparametric statistical testing. Across all datasets, statistical significance was defined as an adjusted P value < 0.05. Code and Data Availability Custom scripts were written in Python (v3.11) using Scikit-learn, HDBSCAN, NetworkX, Scipy, Statsmodels, Scimap, and SHAP. Derived feature matrices are available upon request. Whole-slide images are available under institutional data use agreement due to protected health information. Abbreviations Used in Manuscript | Term | Abbreviation | | - | - | | Artifical Intelligence | AI | | Atypical Localization of Immature Precursors | ALIP | | Clonal Cytopenia of Undetermined Significance | CCUS | | Clonal Hematopoiesis of Indeterminate Potential | CHIP | | C-X-C motif chemokine ligand 12 | CXCL12 | | C-X-C motif chemokine receptor 4 | CXCR4 | | Disease-Modifying Therapy | DMT | | Hematopoietic Stem Cell Transplantation | HSCT | | Hematopoietic Stem and Progenitor Cells | HSPCs | | International Consensus Classification | ICC | | International Prognostic Scoring System-Molecular | IPSS-M | | Myelodysplastic Neoplasms | MDS | | MDS with Increased Blasts | MDS-IB | | MDS with Low Blasts | MDS-LB | | MDS Severity Score | MDSSS | | Multiparameter Flow Cytometry | MFC | | Maturing Myeloid Cells | MMCs | | Multiplex Immunofluorescence | MxIF | | Non-Hematopoietic Elements | NHEs | | Normal Bone Marrow | NBM | | t-distributed Stochastic Neighbor Embedding | t-SNE | | Weill Cornell Medicine/NewYork-Presbyterian Hospital | WCM/NYP | | World Health Organization | WHO | | | ### Page 19 | | | - | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | Whole Slide Images | WSIs | ### Page 20 SUPPLEMENTAL REFERENCES 1 Zeidan AM, Platzbecker U, Bewersdorf JP, Stahl M, Adès L, Borate U et al. Consensus proposal for revised International Working Group 2023 response criteria for higher-risk myelodysplastic syndromes. Blood 2023;141: 2047–2061. 2 van Dongen JJM, Lhermitte L, Böttcher S, Almeida J, van der Velden VHJ, Flores-Montero J et al.EuroFlow antibody panels for standardized n-dimensional flow cytometric immunophenotyping of normal, reactive and malignant leukocytes. Leukemia 2012; 26: 1908–1975. 3 Sarachakov A, Varlamova A, Svekolkin V, Polyakova M, Valencia I, Unkenholz C et al. Spatial mapping of human hematopoiesis at single-cell resolution reveals aging-associated topographic remodeling. Blood 2023;142: 2282–2295. 4 Greenwald NF, Miller G, Moen E, Kong A, Kagel A, Dougherty T et al. Whole-cell segmentation of tissue images with human-level performance using large-scale data annotation and deep learning. Nat Biotechnol2022; 40: 555–565. 5 Healy CP, Adler FR, Deans TL. Rolling Signal-based Ripley’s K: A new algorithm to identify spatial patterns in histological specimens. 2020; : 2020.05.21.109314. ### Page 21 ### Page 22 ### Page 23 ### Page 24 ### Page 25 ### Page 26 ### Page 27 ### Page 28 ### Page 29 ### Page 30 | Supplemental Table 5. Summary of model architectures, training, and validation | | - | | Domain | Network Architecture | Prediction Method | Data Set | Validation Details | Metrics | | Cell Detection Segmentaion | Mesmer: ResNet50 backbone with a multi-resolution feature pyramid. Two semantic heads for pixelwise and inner distance transforms. Model weights pretrained on TissueNet. | 256x256 overlapping patches MxIF channels | 1036 annotated 256x256 patches amounting to around 100,000 individual cell segmentations derived from training set of 6 BM WSI samples. Split 80%/20% into training and validation sets. | 340 annotated 256x256 image patches amounting to around 28,500 individual cell segmentations derived from hold- out set of 3 BM WSI samples. | Jaccard IoU: .71 Recall: .77 Precision: .91 F1-score: .83 | | Membrane Marker Classificaton | Resnet18 with global average pooling layer connected to single neuron binary prediction output. Model weights pretrained on ImageNet. | 128x128 cropped patches centered on the cell of interest. Three channel images: cell segmentation, DAPI, and membrane marker channel of interest. | 872 single cell membrane images annotated as positive or negative expression. Split into 75%/25% training and validation sets. | 155 annotated single cell membrane images. | Accuracy: .97 Precision: 1. Recall: .94 F1-score: .97 | | TP53 Classification | Resnet18 with global average pooling layer connected to three neuron classification output. Model weights pretrained on Membrane Marker Classification dataset and then fine- tuned for nuclear expression classification. | 128x128 cropped patches centered on the cell of interest. Three-channel images: cell segmentation, DAPI, and TP53 channel. | 76 single cell TP53 nuclear expression images annotated as No Expression, Low Expression, Moderate/High Expression. Split into 80%/20% training and validation sets. | 25 annotated singe cell TP53 nuclear expression images. | Accuracy: .95 Precision: 1. Recall: .86 F1-score: .92 | ### Page 31 ### Page 32 SUPPLEMENTARY DATA FIGURES | | | - | | Supplementary Fig. 1. Somatic mutation profiles and cytogenetic abnormalities in CCUS and MDS patients.Stacked heatmap of patient metadata as each column. Top row of heatmap details patient age, sex, IPSS-M score (if applicable), sample group designated for each patient sample according to color legend on right-hand side. Patient sample groups listed are, in order: CCUS, MDS-LB and MDS-LB-SF3B1m, MDS-IB, MDS-IB-TP53m. Below patient metadata are results of targeted next- generation sequencing (only genes mutated in ≥10% samples are shown) with quantification of clonal mutation burden colored by variant allele frequency according to color bar on right-hand side. Lowest segment of heatmap are blue squares denoting cytogenetic abnormalities. | ### Page 33 | | | - | | Supplementary Fig. 2. Single-cell segmentation, phenotyping, and quantitative validation of the multiplex immunofluorescence analysis pipeline.(a) Schematic of the artificial intelligence–driven digital image analysis workflow. Shown are representative raw multiplex immunofluorescence (MxIF) image tiles (top left), patch-level Mesmer-based segmentation masks (top right), stitched whole-slide segmentation masks generated from overlapping 256 × 256 μm patches (bottom left), and convolutional neural network (CNN)– based cell phenotype classification masks used for downstream spatial analyses (bottom right).(b) Ridge plots showing distribution of mean fluorescence intensity (MFI) values across cell types classified by the CNN-based binary marker labeling algorithm, demonstrating concordance between predicted phenotypes and marker signal intensity distributions.(c) Correlation heatmap comparing MxIF-derived cell-type proportions with orthogonal measurements from manual differential counts and multiparameter flow cytometry. Spearman correlation coefficients (ρ) and corresponding P values are shown for each comparison. | ### Page 34 | | | - | | Supplementary Fig. 3. Cell-type frequencies and morphologic characteristics across disease and genotype-defined groups.(a) Box plots showing proportions of hematopoietic stem and progenitor cells (HSPCs), myeloblasts, and megakaryocytes across genotype-defined sample groups. Statistical comparisons were performed using the Mann–Whitney U test (*P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001).(b) Box plots showing proportions of CD38⁺ hematolymphoid cells (non-plasma cells), mast cells, endothelial cells, and non- hematopoietic elements across sample groups (Mann–Whitney U test; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001).(c) Top: Stacked bar plots showing cell-type proportion distributions in tissue microarray clot section samples. Bottom: Proportion of total cells exhibiting high p53 expression (putative TP53-mutated cells) within each sample.(d) Left: Box plots comparing megakaryocyte cell area across sample groups stratified by CD34 expression status. Right: Proportion of CD34⁺ megakaryocytes (mean fluorescence intensity >60) by sample group (Mann–Whitney U test; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001).(e) Box plots comparing erythroid normoblast cell area across sample groups (Mann–Whitney U test; *P ≤ 0.05, **P ≤ 0.01).(f) Scatter plot demonstrating significant inverse correlations between proerythroblast and promyelocyte cell area and IPSS-M risk score. TP53-mutated MDS-IB samples are highlighted (red circles). | ### Page 35 | | | | - | | | ![A screenshot of a computer screen AI-generated content may be incorrect.]() | | - | - | ### Page 36 | | | | - | | Supplementary Fig. 4. Spatial proximity of hematopoietic cell types to vasculature and trabecular bone.(a) Heatmap of permutation-based P values showing enrichment of selected cell types within 15 μm of endothelial structures. Lower P values (yellow) indicate significant spatial enrichment relative to phenotype-label permutations. Individual samples are annotated by diagnostic group (NBM, CCUS, MDS-LB, MDS-IB); IPSS-M scores are shown for MDS samples. Bold column labels represent group-level combined P values calculated using Stouffer’s method.(b) Box plots showing relative distance to vasculature for selected cell types, defined as the mean cell-type distance to endothelium normalized to the mean distance of all other cell types within the same sample. Higher values indicate relative displacement (reduced enrichment). Statistical comparisons were performed using the Mann–Whitney U test (*P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001).(c) Box plots showing G-function based spatial enrichment test to vasculature for selected cell types. The value plotted is the AUC of G function curve representing the net enrichment of each cell type to vasculature, quantified by normalizing the proportion of each cell type within radius of vasculature to the expected proportion by random label shuffle and integrating this value over the range from 0 to 25 μm. Higher values indicate increased spatial enrichment. Statistical comparisons were performed using the Mann–Whitney U test (*P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001)(d) Heatmap of permutation-based P values showing enrichment of selected cell types within 30 μm of trabecular bone. Lower P values (yellow) indicate significant spatial enrichment.(e) Box plots showing relative distance to trabeculae, defined as mean cell-type distance to trabecular mask normalized to the mean distance of all other cell types within the same sample (Mann–Whitney U test; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001).(f) Heatmap of permutation-based P values for enrichment within 15 μm of endothelium in the tissue microarray cohort containing TP53-mutated MDS and TP53–wild-type NBM samples. Columns represent individual samples; rows represent indicated cell types. Bold column labels denote group-level combined P values (Stouffer’s method). HSPCs remain relatively enriched within 15 μm of endothelium in clot section material from NBMs represented on the tissue microarray. | | | - | - | - | | | | - | ### Page 37 | | | - | | Supplementary Fig. 5. Reduced CXCR4 expression in CD34+ progenitors across publicly available MDS transcriptomic and proteomic datasets.(a) Volcano plots showing differential gene expression between healthy and MDS CD34+ progenitors in four independent datasets from three separate publications. CXCR4 is highlighted as significantly down-regulated (FDR-adjusted P < 0.05 & abs(Log2 Fold Change) > .5) in Dataset 1 (GSE136816; bulk RNA-seq of GMPs and CMPs), Dataset 2 (GSE111085; bulk RNA-seq of CD34+ HSPCs, high-grade MDS), Dataset 3 (EGAS00001007568; single-cell RNA-seq of CD34+ HSPCs), and Dataset 1 single-cell RNA- Seq of HSPCs.(b) UMAP embedding of CD34+ HSPCs from Dataset 1 single-cell RNA-seq, annotated by sample condition, with feature plots showing CD34 and CXCR4 expression. | ### Page 38 | | | - | | (c) CXCR4 expression (FPKM) in bulk RNA-seq datasets comparing healthy versus MDS CD34+ progenitors: healthy vs MDS GMPs (adjusted P = 0.008), healthy vs MDS CMPs (adjusted P = 0.02), healthy vs MDS–low-grade HSPCs (adjusted P = 0.17), and healthy vs MDS–high-grade HSPCs (adjusted P = 0.008).(d) Integration of differential expression results across four datasets (two bulk RNA-seq and two single-cell RNA-seq) identified CXCR4 as the only gene significantly and consistently down-regulated (FDR-adjusted P < 0.05) in all datasets.(e) Violin plots from Dataset 3 single-cell RNA-seq showing CXCR4 expression across annotated CD34+ progenitor subsets in healthy versus MDS samples.(f) Venn diagram illustrating overlap of significantly differentially expressed genes across datasets comparing healthy versus MDS CD34+ progenitors; 27 genes were consistently overexpressed in MDS, with selected genes labeled.(g) Mass spectrometry–based quantification of CXCR4 protein abundance (Dataset 4) across defined progenitor populations. Columns represent individual samples ordered from healthy to high-risk MDS/AML (left to right); rows represent cell types. Higher values (yellow) indicate greater CXCR4 protein expression; lower values (gray/blue) indicate reduced expression. | | | | - | ### Page 39 | | | - | | Supplementary Fig. 6. Quantitative assessment of erythroid island architecture in normal and MDS bone marrow.(a) Representative whole-slide images displayed as merged DAPI and CD71 channels (left) and corresponding erythroid normoblast segmentation masks (right) for CCUS, MDS-LB, and MDS-IB samples. Insets show magnified regions of interest (ROIs) highlighting differences in erythroid spatial clustering.(b) Stacked bar plots showing the distribution of erythroid normoblast cluster sizes per sample, derived using HDBSCAN unsupervised clustering. Bars represent the proportion of total erythroid normoblasts within each cluster size category. Cells not assigned to clusters or belonging to clusters of <4 cells are shown in gray; increasing cluster sizes are shown in red. Samples are | ### Page 40 | | | - | | grouped as NBM (ordered by age), CCUS, and MDS (ordered by disease severity), as indicated by the color bar beneath the x-axis.TP53 and SF3B1 variant allele frequencies are annotated for each sample. For each patient, mean cluster density and clustering coefficient (calculated across clusters of >4 cells) are indicated.(c) Box plots comparing mean erythroid cluster density across sample groups (Mann–Whitney U test; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001). Lower cluster density reflects dispersed, poorly formed erythroid aggregates; higher density indicates compact, well-formed erythroid islands.(d) Scatter plot showing positive correlation between mean erythroid cluster density and peripheral hemoglobin concentration in MDS-LB and MDS-IB samples. TP53-mutated cases are denoted by black circles and SF3B1-mutated cases by yellow circles. | ### Page 41 | | | - | | Supplementary Fig. 7. Robustness analysis of atypical localization of immature precursors (ALIP).(a) Left: Bar plot showing ALIP cluster frequency per MDS sample (x-axis), normalized to the total number of promyelocytes and myeloblasts in each biopsy (y-axis). Right: Box plots comparing normalized ALIP cluster frequency across genotype-defined groups (Mann–Whitney U test; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001).(b) Sensitivity analysis of ALIP definition parameters. Bar plot showing ALIP cluster frequency per sample under eight parameter configurations, varying intercellular distance thresholds (12, 16, 18, and 20 μm) and minimum cluster size (≥2 or ≥3 cells). Black bars denote clusters of ≥2 cells; gray bars denote clusters of ≥3 cells. Colored dots indicate the corresponding distance threshold.(c) Box plots comparing genotype-defined groups using six alternative ALIP parameter configurations (specified in subplot titles). Statistical comparisons were performed using Mann–Whitney U tests (*P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001). | | | | - | | | ### Page 42 | | | - | | Supplementary Fig. 8. Structure and quality control of the MDS Microarchitectural Perturbation Score (MDS-MAPS). | ### Page 43 | | | - | | (a) Top: t-distributed stochastic neighbor embedding (t-SNE) of sample-level features (corresponding to Fig. 6B), with individual samples labeled. Bottom: t-SNE plot annotated with TP53 and SF3B1 variant allele frequencies, highlighting a distinct MDS-IB-TP53 cluster characterized by low t-SNE2 values.(b) Top: t-SNE of 256 × 256 μm patch-level feature data. MDS-derived patches segregate from NBM and clonal hematopoiesis patches in patterns reflecting disease severity. Bottom: patch-level MDS-MAPS values (derived from additive weighting of features correlated with disease severity) overlaid as a color-coded feature map on the patch-level t-SNE projection.(c) Bar plots showing features most strongly correlated with the t-SNE1 (left) and t-SNE2 (right) axes at the sample level. Red bars indicate features with strong positive correlation. Features positively correlated with t-SNE1 trend toward the MDS-IB cluster; features positively correlated with t-SNE2 trend toward the MDS-IB-TP53 cluster.(d) Bar plot of adjusted feature weights used in MAPS score computation. Correlative features from t-SNE sample stratification were curated by experts according to known MDS biology.(e) Box plots showing within-sample distributions of patch-level MDS-MAPS values. Samples are grouped by diagnosis (NBM/CCUS in gray; MDS color-coded by IPSS-M risk score).(f) Quality-control metrics per sample, including total patch count, total cell count, total analyzed tissue area (mm²), and distribution of cells per patch (minimum 100 cells per patch threshold applied). | | | | - | | | ### Page 44 | | | - | | Supplementary Fig. 9. MDS-MAPS discriminates CCUS and low-risk MDS trephine biopsies.(a) Leave-one-out cross-validation (LOOCV) for MDS-LB vs. CCUS classification (n = 33 samples). Receiver operating characteristic (ROC) and precision-recall analyses comparing MAPS against mutation burden features: Max VAF AUC = 0.706; Sum VAF AUC = 0.685; MAPS AUC = 0.815. Max VAF AP = 0.757; Sum VAF AP = 0.741; MAPS AP = 0.890. MAPS demonstrated superior discrimination relative to VAF-based mutation burden features across ROC and precision-recall metrics.(b) Top: Box plot comparing baseline mean MAPS between CCUS and MDS-LB groups (CCUS median = -4.8; increase median = -1.1; Mann–Whitney U AUC = 0.86, P = 0.001; effect size r = 0.72). Bottom: Correlation between baseline MAPS and Max VAF from next-generation sequencing (n = 33 patients; Spearman ρ = 0.25, P = 0.17).(c) Feature ranking by univariate discrimination capacity. Left: top 20 spatial features ranked by oriented nonparametric AUC.Right: clinical features ranked by oriented nonparametric AUC. Bar color indicates direction of association with MDS-LB classification. | ### Page 45 | | | - | ### Page 46 | Supplementary Fig. 10. Exploratory modeling of baseline predictors of blast count increase.(a) SHAP (Shapley additive explanation) forest plots showing feature contributions to prediction of blast count increase from baseline biopsies (n = 16 patients; 8 increase, 8 stable). Two L2-regularized logistic regression models (C = 0.5) were evaluated: clinical features (purple; AUC = 0.78) and spatial/MxIF features (black; AUC = 0.94). Each dot represents a patient-level SHAP value (red = blast increase; blue = stable). Diamonds indicate sign-corrected mean SHAP values for the increase group; horizontal lines represent 95% bootstrap confidence intervals (1,000 resamples). Positive values indicate association with blast increase; negative values indicate association with stability. Features are ranked by absolute mean SHAP magnitude within each model.(b) Top: Correlation between baseline MAPS and IPSS-M risk score (n = 13 patients; Spearman ρ = 0.65, P = 0.016). Bottom: Box plot comparing baseline mean MAPS between stable and blast-increase groups (stable median = 5.5; increase median = 7.5; Mann– Whitney U AUC = 0.61, P = 0.505; effect size r = 0.22).(c) Box plots comparing z-scored feature values between stable and blast-increase groups for top discriminating spatial features. Nonparametric AUC was calculated as the normalized Mann–Whitney U statistic (U / [n₀ × n₁]), representing the fraction of correctly ranked stable/increase patient pairs.(d) Feature ranking by univariate discrimination capacity. Left: top 20 spatial features ranked by oriented nonparametric AUC.Right: clinical features ranked by oriented nonparametric AUC. Bar color indicates direction of association with blast increase. | | - | ### Page 47 | | | - | | Supplementary Fig. 11. Multimodal longitudinal assessment in patients with transient remission. | ### Page 48 | | | - | | (a) Longitudinal stacked line plots showing clinical parameters and MAPS over treatment course for two patients with transient remission: WCM44 (left; 7 biopsies) and WCM13 (right; 3 biopsies). Red bars denote active disease; light green bars denote remission. Treatment at biopsy timepoints is annotated (Aza, azacitidine; Dec, decitabine).(b) t-SNE trajectories for WCM44 (left) and WCM13 (right), showing sample-level embeddings overlaid on the full cohort embedding of NBM and MDS cases. Post-treatment samples are highlighted with black circles. Regions enriched in Non-MDS samples are indicated. For WCM13, corresponding whole-slide MxIF images and MDS-MAPS heatmaps at each timepoint are shown, with arrows linking each biopsy to its t-SNE embedding position.(c) Longitudinal assessment of HSPC spatial enrichment to endothelium in WCM44, shown as permutation testing p-values over time (low p-value indicates significant enrichment; high p-value indicates displacement).(d) Box plots comparing sample-level quantifications of selected features between active disease and remission biopsies in WCM44 and WCM13. Individual points are color-coded by patient. |