Screening cell-cell communication in spatial transcriptomics via collective optimal transport

Spatial transcriptomic technologies and spatially annotated single cell RNA-sequencing (scRNA-seq) datasets provide unprecedented opportunities to dissect cell-cell communication (CCC). How to incorporate the spatial information and complex biochemical processes in reconstructing CCC remains a major challenge. Here we present COMMOT to infer CCC in spatial transcriptomics, which accounts for the competition among different ligand and receptor species as well as spatial distances between cells. A novel collective optimal transport method is developed to handle complex molecular interactions and spatial constraints. We introduce downstream analysis tools on spatial directionality of signalings and genes regulated by such signalings using machine learning models. We apply COMMOT to simulation data and eight spatial datasets acquired with five different technologies, showing its effectiveness and robustness in identifying spatial CCC in data with varying spatial resolutions and gene coverages. Finally, COMMOT reveals new CCCs during skin morphogenesis in a case study of human epidermal development. Both the method and the computational package have broad applications in inferring cell-cell interactions within spatial genomics datasets.

to estimate CCC activities from scRNA-seq data 2,3 often using knowledge databases of 43 signaling 4-6 . Most of these methods rely on the expression levels of ligand and receptor 44 pairs and explicitly defined functions. In particular, the products of ligand and receptor 45 gene expression 5,7 or dedicated nonlinear models like Hill function-based formulas 6 are 46 used. In addition, these methods emphasize different aspects of CCC. To name a few, 47 partial differential equation models and outperforms three related optimal transport 120 methods. We then apply COMMOT to analyze scRNA-seq data spatially annotated 121 using paired spatial datasets and five types of spatial transcriptomics data which have 122 different properties in terms of spatial resolution and gene coverage. Finally, we 123 examine a specific system of human epidermal development and reveal connections 124 between CCC and skin development. 125 126

Overview of COMMOT 128
In ligand-receptor interactions, multiple ligand species can bind to the same receptor 129 species, or a receptor species can accommodate multiple ligand species in space (Fig.  130 1a). To incorporate this into an inference model, we develop a method named collective 131 optimal transport (Fig. 1b). This method has three important properties: 1) taking 132 unnormalized distributions and controlling the marginals of the transport plan by the 133 original distributions to maintain the species units for more accurate estimate of the 134 ligand-receptor coupling, 2) allowing infinity entries in the cost matrix to enforce spatial 135 distance limits of CCC to avoid connecting cells that too far apart spatially, and 3) 136 transporting multi-species distributions (ligands) to multi-species distributions (receptors) 137 to account for the multi-species interactions among ligands and receptors (Fig. 1c). 138 spaces (Fig. 2a). 189

190
As an example, in the case of two ligand species and one receptor species (Fig. 2a, b), 191 the two ligands could be Wnt2 and Wnt5a and the receptor could be Fzd5 35 . Then [ ] 192 and [ ] describe the spatial distribution of free Wnt2 and Wnt5a whose production 193 rates are observed from spatial data. Similarly, [ ] describes the spatial distribution of 194 unoccupied Fzd5 whose initial distribution is taken from spatial data. inferred spatial distribution of bound ligand-receptor complexes is compared with the 202 simulated one. In the 1D physical model, the two ligand species are located on the two 203 sides of the location of receptor, respectively and will compete to bind to the receptor 204 with a limit spatial range of interactions. The effects of symmetric ligand locations, 205 biased ligand locations, and different ligand diffusivities are correctly captured by 206 collective optimal transport (Fig. 2a). In the 2D physical model, five interacting species 207 are modeled, and the interaction patterns are also recapitulated (Fig. 2b). In addition to 208 revealing ligand-receptor complexes, collective optimal transport also annotates the 209 source of ligands that contributed to the ligand-receptor complexes (Fig. 2a,b).
Further quantitative comparison is carried out using simulations with randomly 212 generated spatial expression of ligands and receptors in various cases of ligand-213 receptor binding (Fig. 2c). To account to different ligand-receptor binding situations, ten 214 different cases are considered with different numbers of ligand and receptor species 215 and different levels of competition (the ratio of interacting pairs over total number of 216 species) among the species (Fig. 2d). For each case, ten random initializations are 217 generated. The levels of bound ligand-receptor complex inferred by collective optimal 218 transport are compared to the PDE simulation results using the root mean square error 219 and Spearman's correlation coefficient. Our results show that collective optimal 220 transport consistently derives accurate predictions across the different cases and 221 outperforms the results by evaluating each ligand-receptor pair independently using 222 traditional optimal transport. The improvement over the results from the pairwise 223 approach is especially significant when competition among species is high (Fig. 2d, 224 Supplementary Figs. 1-4), for example, cases 4, 5 and 10. A comparison with simulated 225 data demonstrates that collective optimal transport used by COMMOT accurately 226 reconstructs CCC in space. The ability to handle non-probability distribution is 227 especially useful when the total levels of ligand and receptor expressions are 228 significantly different ( Supplementary Fig. 5). Also, incorporating spatial information 229 helps prioritizing interactions among cells that are spatially close ( Supplementary Fig. 6). 230 Moreover, our method outperforms two other variants of optimal transport that also do 231 not require the input data to be probability distributions: unbalanced optimal transport 232 and partial optimal transport (Supplementary Figs. 7-8). In unbalanced optimal transport, Next, we used COMMOT to examine spatial signaling activities among epidermal cells 258 by considering ligand-receptor pairs annotated in the CellChatDB database. Our 259 computational analysis predicts that molecular interactions between ligands GAS6 and 260 PROS1 with their receptor TYRO3 (TYRO3-GAS6 and TYRO3-PROS1) are significant 261 in granular cells and moderately present in basal cells (Fig. 3b). This prediction was 262 confirmed by both immunostaining for proteins (Fig. 3d) and using RNAscope to stain 263 for RNA (Fig. 3e).

Signaling analysis in high spatial resolution ST data 297
We first explore CCC in ST data with high spatial resolution using the CellChatDB 298 database 6 . We analyzed MERFISH data of the mouse hypothalamic preoptic region ligand-receptor pairs and differential signaling activities were identified for the relatively 339 rare cell types (oligodendrocytes, astrocytes, endothelial, microglia, and mural) (Fig. 5c). 340 Among the two most active signaling pathways, WNT signaling is colocalized in local 341 regions, while noncanonical WNT (ncWNT) signaling is widely distributed across the 342 sample (Fig. 5d). COMMOT predicted significant WNT signaling in neurons 343 ( Supplementary Fig. 21), correlating well with critical roles of WNT signaling in neuronal 344 migration and activity in the somatosensory cortex 53 . In addition to cluster 2 and 4 in the differentially expressed genes that match the signaling patterns of each of the CCC-348 induced clusters (Fig. 5f). This analysis reveals known signaling components of the 349 relevant pathways and may represent novel regulators of each pathway. For example, 350 the positive differentially expressed genes associated with cluster 0 (WNT signal 351 senders) include known WNT ligands Wnt5b, Wnt10a, and Wnt2b, and the differentially 352 expressed genes in cluster 3 (WNT signal receivers) include known target genes of 353 WNT signaling pathway such as Gja1 and Acsf2, and the known corresponding 354 intracellular signaling transductors Lrp5 and Lrp6 (Fig. 5f). 355

356
We further jointly analyzed mouse cortex datasets generated with three different 357 technologies including Visium, seqFISH+, and STARmap. Across these three datasets, 358 AGT signaling was identified in neurons. Spatially, neurons in the L2-3 region were 359 identified as strong receivers of AGT ligands across the three datasets. Interestingly, a 360 striped signaling pattern was observed, wherein strong signals within individual layers 361 form stripes, while weak signals form inter-stripe regions. Strong AGT signaling activity 362 among oligodendrocytes was also identified in both STARmap and seqFISH+ datasets. 363 ( Supplementary Fig. 22). In both Visium and seqFISH+ cortex datasets, we inferred 364 WNT signaling to be active across different cortical layers. In both datasets, we 365 identified WNT signaling to be relatively low in layer 5, compared to other layers. This Finally, we applied COMMOT to signaling analysis with Visium 16 spatial transcriptomics 383 data, where each spatial spot consists of multiple cells. By analyzing the breast cancer 384 data with 3798 spots and 36601 genes, we found clear spatial signaling directions of 385 Midkine (MK) signaling which was identified to be the most active (Fig. 6a), and the 386 regions receiving such signals (Fig. 6b). To identify the genes that may be regulated by 387 or regulate CCC, we used the tradeSeq 57 to perform a differential expression test for 388 which the amount of MK reception was used as the cofactor, analogous to a temporal 389 differential expression test where pseudotime is used as the cofactor (Fig. 6c,d). 390 expressed gene with its own unique spatial pattern (Fig. 6c). Furthermore, the level of 393 COL1A1 expression increases while the level of S100G expression decreases as MK 394 signals increase (Fig. 6d). Similar to the temporal DE gene analysis in scRNA-seq data, 395 the signaling DE gene analysis in ST data identifies relation between gene expression 396 and signaling activity, for example, the connection between COL1A1 expression and 397 MK signaling. Therefore, ST data with good coverage of genes and large number of 398 cells or spots is preferred for such analysis. expressed genes. Using this model, we observe that COL1A1 and S100G are distinctly 410 impacted by various MK ligand-receptor pairs (Fig. 6e). Generally, this analysis can be 411 carried out for any ligand-receptor pair that is expressed in the data, for example, the 412 PD1 signaling pathway related to T cell functions ( Supplementary Fig. 25). 413 genes (Fig. 6f). We found significant PSAP signaling activity widely across the tissue 416 ( Fig. 6g), where broad protective roles of PSAP in the nervous system have been 417 discovered 60 . In addition, FGF signaling activity was identified on the border of 418 cerebellar cortex (Fig. 6g), consistent with its known role in patterning of the cerebellum 419 during development 61 . 420 421

Robust identification of spatial signaling direction and downstream targets 422
To further evaluate COMMOT, we assess the method robustness, study the correlation 423 between inferred CCC and the expression of known downstream genes, and compare 424 with two existing methods, CellChat 6 and Giotto 23 that are designed for scRNA-seq and 425 spatial transcriptomics data respectively. with both good depth and spatial organization. We first integrated them using 433 SpaOTsc 31 to generate two imputed datasets, the imputed scRNA-seq with predicted 434 spatial coordinates of cells and the imputed spatial data with predicted expression of the 435 8925 genes from the scRNA-seq data. Using subsampling as a test for robustness, we 436 randomly subsampled the data and compared them to the results from the full dataset. and Jaccard distance (for cluster level signaling) were found to remain relatively small 439 showing that both signaling directions and cluster level signaling are robustly captured 440 in the subsampled data ( Supplementary Fig. 26; See Methods -Evaluation metrics for 441 detailed calculation of the cosine distance and Jaccard distance). The differentially 442 expressed signaling genes were then identified by using the inferred amount of received 443 signal as the co-factor, analogous to the approach for pseudotime 57 . The DE genes 444 identified from subsampled datasets were found to be comparable to the ones identified 445 from the full dataset ( Supplementary Fig. 26). Finally, using Visium human breast cancer and mouse brain datasets, we studied the 476 difference between COMMOT and the three other methods. When assessing the 477 distributions of inter-cluster spatial distance of the identified significant CCCs, the 478 significant CellChat interactions have a relative uniform distribution of the spatial 479 distances. This is expected since CellChat was designed for scRNA-seq data and 480 spatial organization is not considered in the method. Both COMMOT and Giotto 481 prioritize CCC between spatially close clusters while COMMOT identifies more CCC 482 links that are beyond the neighboring cells because Giotto uses cell-contact graph that might limit the range of inferred CCC. CellPhoneDB v3 identifies more CCC links than 484 COMMOT potentially due to its consideration of inter-cluster distances rather than 485

distance among individual cells (Supplementary Figs. 31-36). In the comparison to 486
CellChat, COMMOT uniquely detects CCC that are not differentially expressed in 487 This tool is based on a novel computational method, collective optimal transport, which 511 incorporates both competing marginal distributions and constrained transport plans, two 512 important properties that cannot be dealt with using the current variants of optimal 513 transport. 514

515
We have studied a wide range of data types of different spatial resolutions and gene 516 coverage: in silico spatial transcriptomics data obtained by integrating scRNA-seq and 517 spatial staining data, Visium, Slide-seq, and seqFISH+ spatial transcriptomics. 518 COMMOT can consistently capture the CCC activities from the simulation data and the 519 real spatial data. For instance, dorsal-to-ventral Dpp signaling and anterior-to-posterior 520 Wg signaling were observed in the Drosophila embryo, consistent with previous in vivo 521 data 66, 67 . Our simulations point to a significant TGF-β signaling activity in the 522 hippocampal region and WNT signaling activity in the cerebellar cortex of the mouse 523 brain, coinciding with the known signaling roles of the above pathways in those 524 regions 68,69 . COMMOT is a useful predictive tool to study spatial positions of 525 differentially expressed gene products. In human skin, COMMOT identified that higher 526 WNT signaling increases the expression of BCAM, POSTN, and STMN1. As WNT 527 signaling has known roles in stem cell proliferation 41 , our analysis suggested that the 528 gene products of WNT signaling would likely be expressed in the basal stem cells, which we confirmed by immunofluorescence staining the three genes above of human 530

epithelia. 531
Having said that, we acknowledge that false positives in our inferred CCC are inherently 533 possible because spatial transcriptomics data does not directly represent protein 534 abundancy and our method cannot capture protein-specific modifications, such as 535 protein phosphorylation, glycosylation, proteolytic cleavage into fragments, dimerization, 536 that certainly affect their signaling functions and, thus, CCC mechanisms, that 537 COMMOT aims to infer. The reliability of CCC predictions is expected to significantly 538 improve once the emerging spatial proteomics approaches will mature. 539

540
The spatial distance constraint in the model is used to capture the effect of ligand 541 diffusivity, which is determined by several factors, including protein weight and tortuosity 542 of extra-cellular space 70 . When screening all ligand-receptor pairs, it is difficult to 543 estimate this parameter for each pair accurately. In our model, the local short-range 544 interactions are emphasized even when the spatial distance range is increased 545 ( Supplementary Fig. 38). Thus, when screening many ligand-receptor pairs, a uniform 546 and relatively large spatial distance limit may be used in the model to avoid missing 547 important interactions. Once the important interactions are identified, using an accurate 548 estimation of this parameter would further refine the prediction to remove some false 549 positive CCC links. 550 Like other existing CCC inference methods using scRNA-seq data, COMMOT outputs 552 static CCC activity networks among cells or clusters using spatial data. With a 553 foreseeable deluge of spatial omics data, the temporal sequences of spatial 554 transcriptomics data will soon further enable the revealing of CCC dynamics 71 . When 555 paired with dynamic formulation of optimal transport, the collective optimal transport will 556 be a powerful tool to incorporate multi-omics data for spatiotemporal dynamics of CCC. 557 558 Our PDE model used to generate benchmark simulation data is a flexible and 559 customizable platform for modeling CCC in space. As one possible extension, gene 560 regulatory networks, modeled by a system of ordinary differential equations, could be 561 assigned to each cell, or coupled to each spot of cells to refine the PDE model for more 562 accurate simulation datasets. Building such complex multiscale models is important and 563 tractable as more multi-omics datasets become available, and the collective optimal 564 transport can be integrated with the PDE models to better uncover detailed dynamics of 565

CCC. 566 567
While traditional optimal transport is powerful at integrating a pair of datasets and 568 multimarginal optimal transport 72 integrates multiple datasets, the collective optimal 569 transport is able to effectively deal with competing species. Moreover, collective optimal 570 transport provides a new way to preserve units and control coupled mass. The 571 collective optimal transport in this work will be useful for a broad range of problems 572 beyond the CCC inference. by solving a global optimization problem that accounts for the higher-order interactions 583 among the multiple ligand and receptor species. To this end, we introduce collective 584 optimal transport that seeks a collection of optimal transport plans for all pairs of 585 species that can be coupled simultaneously. As a result, the coupling between a 586 species pair will affect other couplings and vice versa, which cannot be realized in 587 traditional optimal transport 34 . The collective optimal transport results in a large-scale 588 optimization problem for which new algorithms are needed, and thus we presented one 589 based on the efficient Sinkhorn iteration 73 . 590 591 For a ST data of spatial locations and a set of ligand species and receptor 592 species, a collective optimal transport problem is formulated as follows: 593 where , is the expression level of ligand on spot , , is the expression level of 595 receptor on spot , and penalizes the untransported mass and . The coupling 596 matrix , , , scores the signaling strength from spot to spot through the pair of 597 ligands and receptor for , ∈ where is the index set of ligand and receptor 598 species that can bind. The cost matrix , is based on the thresholded distance matrix 599 such that its -th entry equals , ) if , > , and infinity otherwise where is the 600 Euclidean distance matrix among the spots, , is the spatial limit of signaling through 601 the pair of ligand and receptor , and is a scaling function such as square or 602 exponential. When the ligands or receptors contain heteromeric units, the minimum of 603 units is used by default to represent the amount of ligand or receptor. 604 605 Collective optimal transport algorithm 606 To solve the collective optimal transport problem described above, we rewrite the 607 original problem as 608

Evaluation metrics 648
The spatial signaling direction is described by a vector field defined on grid points 649 discretizing the tissue space, and is represented by an array ∈ ℝ × . Cosine distance 650 is used to compare the vector field from subsampled data to the one from the full 651 data and is defined as In the CellChat analysis, the spatial data was treated as a scRNA-seq data, and the 690 count matrix was first normalized using the normalizeData function. The data was then 691      data. a-f CCC analysis of a seqFISH+ data mouse secondary somatosensory cortex. a