CDS: Counterfactual Directionality Score for Structured Interventions in Spatial Graphs

# CDS: Counterfactual Directionality Score for Structured Interventions in Spatial Graphs
Source: [https://arxiv.org/html/2607.13508](https://arxiv.org/html/2607.13508)
Humaira Anzum Md Ishtyaq Mahmud Jagan Mohan Reddy Dwarampudi Tania Banerjee Department of Electrical and Computer Engineering University of Houston, Houston, TX, USA humaira\.snigdha\.22@gmail\.com, ishtyaqtushar@gmail\.com jdwaramp@cougarnet\.uh\.edu, tbanerjee@uh\.edu

###### Abstract

Quantifying directional influence between node populations is a fundamental problem in graph\-based modeling, particularly in spatial biological systems where cell–cell interactions shape functional outcomes\. Existing approaches based on attention, attribution, or correlation capture associations but do not provide a principled framework for evaluating directional effects under controlled perturbations\. We introduce a framework for structured counterfactual interventions in graph\-based models to estimate directional influence between node types\. Our approach trains a Neighbor Influence Model \(NIM\) to predict node states from local neighborhoods and applies constrained interventions that modify neighborhood composition while preserving key spatial and structural properties\. We define the Counterfactual Directionality Score \(CDS\), which measures the change in predicted node state induced by targeted perturbations, and provide a theoretical interpretation of CDS as a finite\-difference measure of local intervention sensitivity\. To obtain valid uncertainty estimates, we introduce a core\-level bootstrap procedure that accounts for dependencies within spatial samples\. Experiments on synthetic spatial graphs with known directional structure show that CDS recovers directional influence, remains well calibrated under null conditions, and is robust to confounding signals, while preliminary results on spatial transcriptomics data reveal biologically plausible and consistent interactions across tissue cores\.

## 1Introduction

Spatially\-resolved single\-cell technologies reveal how cells interact within tissues\. A central question is*directional influence*: how does one cell type affect the molecular state of another? Answering this question is key to understanding tumor–immune interactions, stromal remodeling, and signaling in health and disease\. Graph neural networks \(GNNs\) model such data by representing cells as nodes and spatial proximity as edges\. While GNNs achieve strong predictive performance, they provide limited tools for quantifying directional effects between cell populations\. Existing approaches, such as attention weights, feature attribution, and correlation\-based cell–cell interaction analyses, capture associations but often fail to distinguish true directional influence from confounding factors such as shared microenvironmental signals or global expression trends\.

We introduce a framework for*structured counterfactual interventions*in graph\-based models to estimate directional influence\. Our principle is: if a sender population influences a receiver’s state, systematically removing that population should change the receiver’s predicted cellular state\. We operationalize this idea by training a*Neighbor Influence Model \(NIM\)*that predicts a cell’s transcriptomic state from its spatial neighborhood\. Then we apply counterfactual interventions that modify neighborhood composition while preserving key spatial and structural constraints, including neighborhood size, spatial weighting, and tissue context\. The resulting*Counterfactual Directionality Score \(CDS\)*quantifies the sensitivity of a receiver’s predicted state to targeted perturbations\. We further introduce within‑type interventions as an empirical baseline to distinguish the effects of cell‑type presence from individual cell‑state variability\. We provide a theoretical interpretation of CDS as a finite\-difference measure of local intervention sensitivity under structured neighborhood perturbations\. Our structured interventions isolate compositional effects while controlling for structural confounders\. To obtain valid uncertainty estimates, we develop a core‑level bootstrap procedure that accounts for dependencies within tissue samples\.

On synthetic spatial graphs with known directional structure, CDS recovers ground truth directional structure, remains well calibrated under null conditions, and is robust to confounding signals\. Compared to baseline methods, CDS uniquely combines sensitivity to true interactions with robustness to spurious correlations\. Preliminary results on spatial transcriptomics datasets show that CDS identifies biologically plausible sender–receiver relationships with consistent behavior across tissue cores\.

#### Contributions\.

Our main contributions are: \(1\) a framework for structured counterfactual interventions in spatial graphs that preserves key structural properties; \(2\) CDS, a model\-based measure of directional influence between cell types; \(3\) a theoretical characterization of CDS as a finite\-difference intervention measure under structured perturbations; \(4\) a core\-level bootstrap procedure for statistically principled inference; and \(5\) empirical validation showing that CDS recovers directional structure, remains well calibrated under null conditions, and is robust to confounding\.

## 2Related Work

Our work connects several research areas: interpretability for graph neural networks, attention\-based models, computational methods for cell–cell interaction inference, and causal machine learning with counterfactual reasoning\.

Several methods have been developed to explain predictions of GNNs\. GNNExplainer\[[11](https://arxiv.org/html/2607.13508#bib.bib1)\]identifies compact subgraphs and node features that are most influential for a GNN’s predictions\. PGExplainer\[[7](https://arxiv.org/html/2607.13508#bib.bib2)\]provides a parameterized approach that efficiently generates explanations for multiple instances\. While these methods reveal which graph components are important, they focus on attribution rather than estimating directional effects between node populations\. Our CDS explicitly perturbs neighbourhood composition in a structured way to estimate directional influence\.

Attention mechanisms in GNNs, such as Graph Attention Networks \(GATs\)\[[9](https://arxiv.org/html/2607.13508#bib.bib4)\], learn adaptive weights over neighbours and have been used to interpret which interactions are “important”\. However, attention weights reflect correlations within the observed data distribution and do not provide a counterfactual characterization of directionality\. Our work instead uses explicit counterfactual interventions to isolate the effect of a specific cell type’s presence\.

In computational biology, tools such as CellPhoneDB\[[4](https://arxiv.org/html/2607.13508#bib.bib5)\]and NicheNet\[[3](https://arxiv.org/html/2607.13508#bib.bib6)\]infer ligand–receptor interactions from single\-cell and spatial transcriptomics data using prior knowledge of molecular interactions\. More recent spatial CCI methods\[[2](https://arxiv.org/html/2607.13508#bib.bib7),[6](https://arxiv.org/html/2607.13508#bib.bib8)\]incorporate spatial proximity to score interactions\. Though these methods are powerful for hypothesis generation, they are primarily correlation\-based or rely on curated databases\. Our framework is model\-agnostic and learns a predictive model of cell state from neighbourhood composition\. This enables principled counterfactual estimation without requiring predefined interaction databases\.

Causal inference provides a formal language for estimating intervention effects\[[8](https://arxiv.org/html/2607.13508#bib.bib9)\]\. In machine learning, counterfactual explanations have been applied to image and tabular data \(e\.g\.,\[[10](https://arxiv.org/html/2607.13508#bib.bib10)\]\)\. For graphs, recent work has proposed counterfactual graph generation for model explanation \(e\.g\.,\[[1](https://arxiv.org/html/2607.13508#bib.bib11)\]\), but these approaches typically modify graph structure globally or focus on graph\-level classification\. To our knowledge, our work is the first to combine a learned graph\-based surrogate model with spatially constrained, cell\-type\-specific counterfactual interventions to estimate directional influence between node populations in spatial biological graphs\. The core\-level bootstrap further provides statistically principled inference that respects the dependence structure of tissue samples\.

## 3Problem Setup

We consider a spatially\-resolved single\-cell dataset represented as a graphG=\(V,E\)G=\(V,E\), where each nodei∈Vi\\in Vcorresponds to a cell, and edges\(i,j\)∈E\(i,j\)\\in Eencode spatial proximity between cells\. Each cell is associated with \(1\) a feature vectorxi∈ℝpx\_\{i\}\\in\\mathbb\{R\}^\{p\}representing its transcriptomic state, and \(2\) a discrete cell\-type labelτ​\(i\)∈𝒯\\tau\(i\)\\in\\mathcal\{T\}\. LetX∈ℝ\|V\|×pX\\in\\mathbb\{R\}^\{\|V\|\\times p\}denote the matrix of all cell features\.

#### Spatial neighborhood\.

For each cellii, we define a local neighborhoodN​\(i\)⊂VN\(i\)\\subset Vbased on spatial proximity\. In practice,N​\(i\)N\(i\)may be constructed viakk\-nearest neighbors in physical space, but our formulation is agnostic to the specific construction as long as it captures local cell–cell interactions\.

#### Sender and receiver populations\.

For a pair of cell types\(s,r\)∈𝒯2\(s,r\)\\in\\mathcal\{T\}^\{2\}, we define the sender and receiver populations:

S\\displaystyle S=\{i∈V:τ​\(i\)=s\},\\displaystyle=\\\{i\\in V:\\tau\(i\)=s\\\},\(1\)R\\displaystyle R=\{i∈V:τ​\(i\)=r\}\.\\displaystyle=\\\{i\\in V:\\tau\(i\)=r\\\}\.\(2\)

#### Interaction\-driven cell state\.

We view the observed transcriptomic state of a receiver celli∈Ri\\in Ras the result of both intrinsic factors and extrinsic influences from its local neighborhood:

xi=f​\(xiint,\{xj:j∈N​\(i\)\}\)\+ϵix\_\{i\}=f\\big\(x\_\{i\}^\{\\mathrm\{int\}\},\\\{x\_\{j\}:j\\in N\(i\)\\\}\\big\)\+\\epsilon\_\{i\}\(3\)wherexiintx\_\{i\}^\{\\mathrm\{int\}\}denotes intrinsic \(cell\-autonomous\) components and the neighborhood term captures microenvironmental influence\.ϵi\\epsilon\_\{i\}captures unmodeled variability\.

#### Problem objective\.

Our goal is to quantify the*directional influence*of sender\-type cells on receiver\-type cells\. Formally, for each receiver celli∈Ri\\in R, we define a neighborhood composition variableZi\(s\)Z\_\{i\}^\{\(s\)\}that summarizes the contribution of sender\-type cells inN​\(i\)N\(i\)\. We are interested in estimating the effect:

Δi\(s→r\)=f​\(xiint,Zi\(s\)\)−f​\(xiint,0\),\\Delta\_\{i\}^\{\(s\\to r\)\}=f\\big\(x\_\{i\}^\{\\mathrm\{int\}\},Z\_\{i\}^\{\(s\)\}\\big\)\-f\\big\(x\_\{i\}^\{\\mathrm\{int\}\},0\\big\),\(4\)which represents the change in the receiver cell state induced by removing sender\-type cells from the local neighborhood\.

#### Scientific question\.

This formulation allows us to address the central question:*Do neighboring sender\-type cells influence the transcriptomic program of receiver\-type cells, and if so, what is the magnitude and direction of this effect?*

## 4Neighbor Influence Model

### 4\.1Distance\-Weighted Neighborhood Aggregation

We model the influence of a cell’s local microenvironment through a distance\-weighted aggregation of its neighbors\. For each celliiand neighborj∈N​\(i\)j\\in N\(i\), letdi​jd\_\{ij\}denote their Euclidean distance in tissue space\. We define normalized weights:

wi​j=exp⁡\(−di​j/α\)∑j′∈N​\(i\)exp⁡\(−di​j′/α\),w\_\{ij\}=\\frac\{\\exp\\bigl\(\-d\_\{ij\}/\\alpha\\bigr\)\}\{\\sum\_\{j^\{\\prime\}\\in N\(i\)\}\\exp\\bigl\(\-d\_\{ij^\{\\prime\}\}/\\alpha\\bigr\)\},\(5\)whereα\>0\\alpha\>0controls the spatial decay of influence\.

Using these weights, we define the neighborhood representationai∈ℝpa\_\{i\}\\in\\mathbb\{R\}^\{p\}:

ai=∑j∈N​\(i\)wi​j​xj\.a\_\{i\}=\\sum\_\{j\\in N\(i\)\}w\_\{ij\}x\_\{j\}\.\(6\)

### 4\.2Modeling Neighborhood Influence

We approximate the unknown functionffintroduced in Section[3](https://arxiv.org/html/2607.13508#S3)using a graph\-based parametric modelfθf\_\{\\theta\}, referred to as NIM:

z^i=fθ​\(ai,τ​\(i\)\)∈ℝp\.\\hat\{z\}\_\{i\}=f\_\{\\theta\}\(a\_\{i\},\\tau\(i\)\)\\in\\mathbb\{R\}^\{p\}\.\(7\)
The model takes as input the neighborhood representationaia\_\{i\}and the receiver cell typeτ​\(i\)\\tau\(i\), and predicts the receiver’s transcriptomic state\. Learningfθf\_\{\\theta\}provides a data\-driven estimate of how local neighborhood composition maps to cellular gene expression\.

#### Architecture\.

fθf\_\{\\theta\}is implemented as a multi\-layer feedforward network with residual connections\. The inputaia\_\{i\}is projected into a latent space, processed byLLresidual blocks, and modulated by a cell\-type\-conditioned gating mechanism\. Specifically, each cell typeτ​\(i\)\\tau\(i\)is embedded into a vectoreτ​\(i\)e\_\{\\tau\(i\)\}, which is transformed into a feature\-wise gate:

gi=σ​\(Wtype​eτ​\(i\)\+btype\),g\_\{i\}=\\sigma\(W\_\{\\mathrm\{type\}\}e\_\{\\tau\(i\)\}\+b\_\{\\mathrm\{type\}\}\),\(8\)that multiplicatively modulates the hidden representation\.

This gating allows the model to capture heterogeneous response functions, i\.e\., different receiver cell types may respond differently to the same neighborhood signal\.

### 4\.3Training Objective

We trainfθf\_\{\\theta\}to minimize the discrepancy between predictionsz^i\\hat\{z\}\_\{i\}and observed statesxix\_\{i\}using the Huber loss:

ℒ​\(θ\)=1\|ℬ\|​∑i∈ℬℓδ​\(z^i,xi\),\\mathcal\{L\}\(\\theta\)=\\frac\{1\}\{\|\\mathcal\{B\}\|\}\\sum\_\{i\\in\\mathcal\{B\}\}\\ell\_\{\\delta\}\(\\hat\{z\}\_\{i\},x\_\{i\}\),\(9\)whereℓδ\\ell\_\{\\delta\}is defined as the Huber loss\.

#### Interpretation\.

A low prediction error implies that neighborhood composition contains predictive information about the receiver cell state\. Importantly, this learned mappingfθf\_\{\\theta\}serves as a surrogate model for evaluating counterfactual neighborhood perturbations \(Section[6](https://arxiv.org/html/2607.13508#S6)\), enabling estimation of sensitivity to sender\-to\-receiver perturbations\.

## 5CDS

### 5\.1Intuition

While NIM captures how neighborhood composition predicts cell state, it does not by itself identify which neighboring cell types drive this effect\. To address this, we introduceCDS, which measures the sensitivity of a receiver cell’s predicted state to targeted perturbations of its neighborhood\.

### 5\.2Counterfactual Predictions

Letfθf\_\{\\theta\}be the trained NIM\. For a receiver cellii, the original prediction is:

z^i=fθ​\(ai,τ​\(i\)\),\\hat\{z\}\_\{i\}=f\_\{\\theta\}\(a\_\{i\},\\tau\(i\)\),\(10\)whereaia\_\{i\}is the original neighborhood aggregation\.

We construct a counterfactual neighborhoodNcf​\(i\)N^\{\\mathrm\{cf\}\}\(i\)via an intervention operatorℐS→R\\mathcal\{I\}\_\{S\\rightarrow R\}that modifies the presence of sender\-type cells while preserving the spatial weighting structure\. The corresponding counterfactual prediction is:

z^icf=fθ​\(aicf,τ​\(i\)\),\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}=f\_\{\\theta\}\(a\_\{i\}^\{\\mathrm\{cf\}\},\\tau\(i\)\),\(11\)whereaicfa\_\{i\}^\{\\mathrm\{cf\}\}denotes the neighborhood aggregation computed from the counterfactual neighborhoodNcf​\(i\)N^\{\\mathrm\{cf\}\}\(i\)\.

### 5\.3Intervention Operators

We define two structured intervention operators that modify neighborhood composition while preserving key spatial and structural properties\. These interventions are designed to isolate the contribution of sender\-type cells to receiver predictions\.

#### Type\-Swap Intervention\.

This operator removes the presence of sender\-type cells from the receiver’s neighborhood\. For each receiveri∈Ri\\in Rand each neighborj∈N​\(i\)j\\in N\(i\)withτ​\(j\)=s\\tau\(j\)=s, we replacejjwith a cellj′j^\{\\prime\}sampled from the same tissue core such thatτ​\(j′\)≠s\\tau\(j^\{\\prime\}\)\\neq s\.

To ensure physically plausible perturbations, we perform*distance\-bin\-preserving replacement*\. Thekkneighbor slots are partitioned into spatial bins based on distance, and each sender slot is preferentially replaced by a non\-sender cell from the same bin; if no such candidate exists, the search expands to nearby bins\.

This construction preserves three structural invariants:

1. 1\.Degree preservation:\|Ncf​\(i\)\|=\|N​\(i\)\|=k\|N^\{\\mathrm\{cf\}\}\(i\)\|=\|N\(i\)\|=k\.
2. 2\.Microenvironment consistency: replacement cells are drawn from the same tissue core asii\.
3. 3\.Distance\-slot preservation: spatial weights\{wi​j\}\\\{w\_\{ij\}\\\}remain unchanged, so differences betweenaia\_\{i\}andaicfa\_\{i\}^\{\\mathrm\{cf\}\}arise solely from changes in neighborhood composition\.

As a result, changes in prediction are attributable primarily to changes in neighborhood composition rather than artifacts of graph structure or spatial weighting\.

#### Within\-Type Intervention\.

This operator help characterize the effect of*sender cell identity*while preserving sender cell type\. For each sender neighborj∈N​\(i\)j\\in N\(i\)withτ​\(j\)=s\\tau\(j\)=s, we replacejjwith a different cellj′j^\{\\prime\}of the same type, drawn from the donor pool:

𝒟i=\{j′∈V:τ​\(j′\)=s,core​\(j′\)=core​\(i\),j′≠i\}\.\\mathcal\{D\}\_\{i\}=\\bigl\\\{\\,j^\{\\prime\}\\in V:\\tau\(j^\{\\prime\}\)=s,\\;\\mathrm\{core\}\(j^\{\\prime\}\)=\\mathrm\{core\}\(i\),\\;j^\{\\prime\}\\neq i\\,\\bigr\\\}\.\(12\)
If𝒟i=∅\\mathcal\{D\}\_\{i\}=\\emptyset, the receiveriiis excluded from analysis\. This intervention preserves the presence of sender\-type cells while randomizing which specific cells occupy neighborhood slots\.

### 5\.4CDS Definition

For each receiver cellii, we define:

CDSi=1p​‖z^icf−z^i‖1\.\\mathrm\{CDS\}\_\{i\}=\\frac\{1\}\{p\}\\left\\\|\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i\}\\right\\\|\_\{1\}\.\(13\)This measures the average change in predicted state under the counterfactual intervention\.

A largerCDSi\\mathrm\{CDS\}\_\{i\}indicates that the perturbed cell type substantially influences the receiver’s transcriptomic state\.

We also define a signed variant:

CDSisgn=1p​∑m=1p\(z^i,mcf−z^i,m\),\\mathrm\{CDS\}\_\{i\}^\{\\mathrm\{sgn\}\}=\\frac\{1\}\{p\}\\sum\_\{m=1\}^\{p\}\\bigl\(\\hat\{z\}\_\{i,m\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i,m\}\\bigr\),\(14\)which captures the direction of influence \(upregulation vs\. suppression\)\.

Finally, we aggregate over receiver cells:

CDSS→R=1\|R∗\|​∑i∈R∗CDSi,\\mathrm\{CDS\}\_\{S\\rightarrow R\}=\\frac\{1\}\{\|R^\{\*\}\|\}\\sum\_\{i\\in R^\{\*\}\}\\mathrm\{CDS\}\_\{i\},\(15\)whereR∗R^\{\*\}includes receivers for which the intervention is valid\.

#### Interpretation\.

CDS provides a model\-based measure of directional influence: it quantifies how much the predicted state of receiver cells depends on the presence or identity of sender\-type cells in their local neighborhood\.

## 6Theoretical Properties of Counterfactual Directionality

CDS can be interpreted as a finite\-difference measure of the sensitivity of model predictions to structured neighborhood perturbations\. The differenceaicf−aia\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}captures changes in neighborhood composition induced by the intervention while preserving spatial structure and neighborhood size\.

Structured interventions preserve neighborhood size, spatial distances \(and therefore weightswi​jw\_\{ij\}\), and receiver type, so

aicf−ai=∑m∈N​\(i\)wi​m​\(xmcf−xm\),a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}=\\sum\_\{m\\in N\(i\)\}w\_\{im\}\\bigl\(x\_\{m\}^\{\\mathrm\{cf\}\}\-x\_\{m\}\\bigr\),thereby separating compositional changes from changes in graph topology\.

Under within\-type exchangeability,

𝔼​\[aicf−ai\]=0,\\mathbb\{E\}\[a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\]=0,providing an empirical null baseline\. Because CDS depends on the aggregated neighborhood representationaia\_\{i\}, it captures only effects mediated through this representation\. Iffθf\_\{\\theta\}isLL\-Lipschitz, then

CDSi≤Lp​‖aicf−ai‖2\.\\mathrm\{CDS\}\_\{i\}\\leq\\frac\{L\}\{p\}\\\|a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\\\|\_\{2\}\.
For statistical inference, we use a core\-level bootstrap procedure: we resampleMMtissue cores with replacement, compute bootstrap meansμ\(b\)\\mu^\{\(b\)\}, and construct percentile confidence intervals to account for dependence within tissue cores\.

#### Implementation\.

We train NIM using standard optimization procedures and compute CDS via structured counterfactual perturbations of local neighborhoods\. Full algorithmic details are provided in Appendix[A](https://arxiv.org/html/2607.13508#A1)\.

## 7Experiments

We evaluate whether CDS recovers directional influence under controlled synthetic settings and compare it against standard attribution and perturbation baselines\. Our evaluation focuses on directionality, robustness to confounding, and stability\. These experiments empirically validate the theoretical properties of CDS derived in Section[6](https://arxiv.org/html/2607.13508#S6), including its sensitivity to neighborhood perturbations and robustness under null exchangeability\.

### 7\.1Experimental Setup

#### Synthetic Spatial Generator

We construct spatial graphs with three cell types: sender \(SS\), receiver \(RR\), and background \(BB\)\. Cell coordinates are sampled uniformly and akk\-nearest neighbor graph is built with distance\-based softmax weights\.

Receiver features are generated under three regimes:

- •Positive \(true influence\): zi=μR\+ai​W\+ϵ,z\_\{i\}=\\mu\_\{R\}\+a\_\{i\}W\+\\epsilon,whereμR∈ℝp\\mu\_\{R\}\\in\\mathbb\{R\}^\{p\}is a receiver\-type baseline expression vector,aia\_\{i\}is the weighted aggregation of neighboring sender features,W∈ℝp×pW\\in\\mathbb\{R\}^\{p\\times p\}controls sender\-to\-receiver influence strength, andϵ∼𝒩​\(0,σ2​I\)\\epsilon\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}I\)is Gaussian noise\.
- •Null \(no influence\): zi∼𝒩​\(μR,σ2​I\),z\_\{i\}\\sim\\mathcal\{N\}\(\\mu\_\{R\},\\sigma^\{2\}I\),independent of neighboring cells\.
- •Spurious \(confounded\):A shared latent variable induces correlation across cell types without local directional influence\.

This design preserves feature marginals while controlling the underlying directional dependency structure\.

#### Predictive Model

We train NIM, a residual neural network that predicts cell features from neighborhood aggregates with cell\-type gating\. CDS is computed by comparing model predictions under original and counterfactual neighborhoods constructed via constrained neighbor replacement\.

### 7\.2Synthetic Validation

#### Directional Recovery

We first evaluate whether CDS correctly recovers directional structure under controlled synthetic settings\. Figure[1](https://arxiv.org/html/2607.13508#S7.F1)shows that in the positive regime, CDS consistently assigns higher values to the true direction \(S→RS\\rightarrow R\) than the reverse \(R→SR\\rightarrow S\) across all noise levels and random seeds\. This gap remains stable as noise increases, indicating that the learned model captures directional dependencies rather than symmetric correlations\. In contrast, in the null regime, both directions yield near\-zero CDS values, confirming that the method does not produce spurious directional signal \(Table[1](https://arxiv.org/html/2607.13508#S7.T1)\) and is well calibrated in the absence of true interaction\.

#### Null Control

Under the null setting, where receiver features are independent of neighboring sender cells, CDS values collapse to near zero for both directions \(Table[1](https://arxiv.org/html/2607.13508#S7.T1)\), confirming that the method does not introduce spurious directional signal under null conditions\.

#### Spurious Correlation Control

To test robustness to confounding, we introduce a shared latent signal affecting all cell types\. As shown in Figure[1](https://arxiv.org/html/2607.13508#S7.F1)and Table[1](https://arxiv.org/html/2607.13508#S7.T1), CDS values remain lower than in the positive regime, and the directional gap is reduced\.

This indicates that CDS preferentially captures structured local influence rather than global correlations, and does not incorrectly attribute strong directional effects under confounding\.

![Refer to caption](https://arxiv.org/html/2607.13508v1/1.png)Figure 1:Directional CDS across synthetic regimes\.Boxplots comparingCDSS→R\\mathrm\{CDS\}\_\{S\\rightarrow R\}andCDSR→S\\mathrm\{CDS\}\_\{R\\rightarrow S\}in positive, null, and spurious regimes\.

### 7\.3Quantitative Evaluation

We evaluate the ability of CDS to distinguish true interaction from null and confounded settings using ROC analysis\. As shown in Figure[2](https://arxiv.org/html/2607.13508#S7.F2), CDS achieves high AUC for distinguishing positive from null regimes \(AUC = 0\.97\) and maintains above\-chance discrimination under confounding \(AUC = 0\.79\)\. These results indicate that CDS captures structured directional signals beyond simple correlation\-based separation\. Table[1](https://arxiv.org/html/2607.13508#S7.T1)summarizes CDS magnitude and stability across regimes\. In the positive setting, CDS is substantially larger than in null and spurious regimes, with a low coefficient of variation across bootstrap samples\. Bootstrap analysis further shows that CDS estimates are tightly concentrated, with narrow confidence intervals, indicating stable estimation across tissue cores\. We compare CDS against random perturbation and gradient sensitivity baselines \(Table[2](https://arxiv.org/html/2607.13508#S7.T2)\)\. CDS consistently achieves higher scores in the positive regime while remaining near zero under null conditions\.

![Refer to caption](https://arxiv.org/html/2607.13508v1/15.png)Figure 2:Influence detection performance\.ROC curves for distinguishing positive influence from null and spurious controls using CDS scores\.Table 1:Summary of CDS performance across regimes\.Table 2:Baseline comparison\.Mean normalized scores across regimes\. Gradient sensitivity values are rescaled for comparability\.Table 3:Ablation study of structural constraints\.Removing structural constraints alters CDS magnitude and/or increases variability relative to the full intervention design\.
### 7\.4Ablation Study

We evaluate the importance of structural constraints in the counterfactual intervention\. Table[3](https://arxiv.org/html/2607.13508#S7.T3)shows that removing degree preservation, distance preservation, or within\-type baseline correction alters CDS magnitude and increases variability relative to the full method\. In particular, unmatched donor replacement introduces variability by altering neighborhood composition, while removing distance preservation disrupts spatial weighting, leading to noisier estimates\. These results are consistent with the theoretical motivation for the structured intervention design \(Section[6](https://arxiv.org/html/2607.13508#S6)\) and suggest that preserving neighborhood structure improves the stability and specificity of CDS estimates\.

### 7\.5Robustness Analysis

#### Noise and Interaction Robustness\.

As observation noise increases, CDS magnitude decreases but the directional gap betweenS→RS\\rightarrow RandR→SR\\rightarrow Spersists, showing graceful degradation\. Core\-level bootstrap yields low variance and tight confidence intervals, confirming stability \(Figure[3](https://arxiv.org/html/2607.13508#S7.F3)\)\. Increasing ground\-truth interaction strength produces monotonic increases in CDS, consistent with the finite\-difference sensitivity interpretation of CDS \(Section[6](https://arxiv.org/html/2607.13508#S6)\)\.

![Refer to caption](https://arxiv.org/html/2607.13508v1/5.png)

![Refer to caption](https://arxiv.org/html/2607.13508v1/6.png)

![Refer to caption](https://arxiv.org/html/2607.13508v1/7.png)

Figure 3:Robustness analysis of CDS\.\(a\) Noise robustness: meanCDSS→R\\mathrm\{CDS\}\_\{S\\rightarrow R\}andCDSR→S\\mathrm\{CDS\}\_\{R\\rightarrow S\}across increasing observation noise\. \(b\) Bootstrap stability: distribution of bootstrap CDS estimates across resampled tissue cores\. \(c\) Influence\-strength sensitivity: CDS as a function of ground\-truth sender\-to\-receiver interaction strength\.

### 7\.6Real\-World Data

We applied our CDS framework to two independent spatial transcriptomics datasets: breast cancer \(BRCA\) and lung cancer tissue microarrays\[[5](https://arxiv.org/html/2607.13508#bib.bib3)\]\. For each dataset, we considered three major cell populations: Tumor, Stromal, and Immune—resulting in four directional hypotheses, namely, Tumor→\\rightarrowStromal, Stromal→\\rightarrowTumor, Tumor→\\rightarrowImmune, and Stromal→\\rightarrowImmune\. We trained separate NIMs for each cancer type using the same architecture and hyperparameters, computed CDS via type\-swap interventions with core\-level bootstrap validation \(1000 iterations; 3 test cores per cancer type\), and report mean CDS values with 95% bootstrap confidence intervals\.

Figure[4](https://arxiv.org/html/2607.13508#S7.F4)shows CDS estimates for all four sender–receiver pairs in BRCA and lung cancer\. Tumor→\\rightarrowStromal influence \(CDS = 0\.0715 \[0\.0550, 0\.0758\]\) substantially exceeded Stromal→\\rightarrowTumor influence \(CDS = 0\.0247 \[0\.0191, 0\.0271\]\), indicating asymmetric directional signaling from malignant to stromal compartments\. We observed similar moderate scores for Tumor→\\rightarrowImmune \(0\.0657 \[0\.0598, 0\.0677\]\) and Stromal→\\rightarrowImmune \(0\.0687 \[0\.0576, 0\.0784\]\), suggesting that immune cell states are influenced by both tumor and stromal microenvironments\. The confidence intervals for Tumor→\\rightarrowStromal and Stromal→\\rightarrowTumor do not overlap, providing evidence for directional asymmetry\.

![Refer to caption](https://arxiv.org/html/2607.13508v1/11.png)Figure 4:Comparative CDS estimates for breast and lung cancer\.In the lung tissue, a clear directional signature can be seen \(Figure[4](https://arxiv.org/html/2607.13508#S7.F4)\)\. Tumor→\\rightarrowStromal influence \(0\.0917 \[0\.0837, 0\.0965\]\) was again the strongest observed effect, exceeding the BRCA magnitude and showing the largest directional gap relative to Stromal→\\rightarrowTumor \(0\.0315 \[0\.0277, 0\.0368\]\)\. However, immune\-related interactions diverged from BRCA: Tumor→\\rightarrowImmune \(0\.0358 \[0\.0328, 0\.0591\]\) and Stromal→\\rightarrowImmune \(0\.0523 \[0\.0495, 0\.0583\]\) were substantially lower in lung\. This suggests tissue\-specific immune interaction patterns\. The Stromal→\\rightarrowImmune signal remained moderately elevated relative to Tumor→\\rightarrowImmune, a pattern reversed from BRCA\. Higher immune\-directed CDS values in BRCA than in lung suggest that the framework captures biologically plausible differences in tumor–immune crosstalk across cancer contexts\.

## 8Conclusion

We introduced CDS, a framework for structured counterfactual interventions in spatial graphs to quantify directional influence between cell populations\. By combining neighborhood\-based predictive modeling with spatially constrained counterfactual perturbations, CDS provides a principled measure of sensitivity to sender\-to\-receiver perturbations while preserving key structural properties of the tissue microenvironment\. Across controlled synthetic experiments, CDS recovers directional structure, remains well calibrated under null conditions, and is robust to confounding signals\. Preliminary results on spatial transcriptomics datasets further demonstrate biologically plausible and consistent interaction patterns across tissue cores\. More broadly, our results suggest that structured counterfactual interventions provide a useful framework for studying directional dependencies in spatially organized biological systems\.

## References

- \[1\]\(2022\)Counterfactual generation for graph neural networks\.InProceedings of the 39th International Conference on Machine Learning \(ICML\),pp\. 42–57\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p5.1)\.
- \[2\]E\. Armingol, A\. Officer, O\. Harismendy, and N\. E\. Lewis\(2021\)Deciphering cell–cell interactions and communication from gene expression\.Nature Reviews Genetics22\(2\),pp\. 71–88\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p4.1)\.
- \[3\]R\. Browaeys, W\. Saelens, and Y\. Saeys\(2020\)NicheNet: modeling intercellular communication by linking ligands to target genes\.Nature Methods17\(2\),pp\. 159–162\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p4.1)\.
- \[4\]M\. Efremova, R\. Vento\-Tormo, S\. A\. Teichmann, and M\. Vento\-Tormo\(2020\)CellPhoneDB: inferring cell–cell communication from combined expression of multi\-subunit ligand–receptor complexes\.Nature Protocols15\(4\),pp\. 1484–1506\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p4.1)\.
- \[5\]Gene Expression Omnibus\(2025\)GSE308148\.Note:[https://www\.ncbi\.nlm\.nih\.gov/geo/query/acc\.cgi?acc=GSE308148](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE308148)Accessed: 2026\-02\-21Cited by:[§7\.6](https://arxiv.org/html/2607.13508#S7.SS6.p1.4)\.
- \[6\]J\. Liao, Y\. Lu, J\. Zhang, X\. Li, X\. Shen, and S\. Li\(2022\)SpaTalk: a graph attention network to infer cell\-type\-specific spatial signaling\.Nucleic Acids Research50\(12\),pp\. e70\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p4.1)\.
- \[7\]D\. Luo, W\. Cheng, D\. Xu, W\. Yu, B\. Zong, H\. Chen, and X\. Zhang\(2020\)Parameterized explainer for graph neural networks\.InAdvances in Neural Information Processing Systems 33 \(NeurIPS\),pp\. 19620–19631\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p2.1)\.
- \[8\]J\. Pearl\(2009\)Causality: models, reasoning, and inference\.2nd edition,Cambridge University Press\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p5.1)\.
- \[9\]P\. Veličković, G\. Cucurull, A\. Casanova, A\. Romero, P\. Liò, and Y\. Bengio\(2018\)Graph attention networks\.InInternational Conference on Learning Representations \(ICLR\),Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p3.1)\.
- \[10\]S\. Wachter, B\. Mittelstadt, and C\. Russell\(2017\)Counterfactual explanations without opening the black box: automated decisions and the gdpr\.Harvard Journal of Law & Technology31\(2\),pp\. 841–887\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p5.1)\.
- \[11\]R\. Ying, D\. Bourgeois, J\. You, M\. Zitnik, and J\. Leskovec\(2019\)GNNExplainer: generating explanations for graph neural networks\.InAdvances in Neural Information Processing Systems 32 \(NeurIPS\),pp\. 9244–9255\.Cited by:[§2](https://arxiv.org/html/2607.13508#S2.p2.1)\.

## Appendix AAlgorithms

We present the algorithms[1](https://arxiv.org/html/2607.13508#alg1),[2](https://arxiv.org/html/2607.13508#alg2), and[3](https://arxiv.org/html/2607.13508#alg3)that together constitute our pipeline:

Algorithm 1NIM Training1:Spatial graph

G=\(V,E\)G=\(V,E\); cell features

X∈ℝ\|V\|×pX\\in\\mathbb\{R\}^\{\|V\|\\times p\}; cell type labels

τ\\tau; training cores

𝒞train\\mathcal\{C\}\_\{\\text\{train\}\}; validation cores

𝒞val\\mathcal\{C\}\_\{\\text\{val\}\}; number of neighbors

k=20k=20; temperature

τagg=1\.0\\tau\_\{\\text\{agg\}\}=1\.0\.

2:Trained NIM parameters

θ∗\\theta^\{\*\}\.

3:// Step 1: Build spatial KNN graph

4:foreach cell

i∈Vi\\in Vdo

5:Compute

N​\(i\)←N\(i\)\\leftarrowkknearest neighbors of

iiby Euclidean distance over

\(xcentroid,ycentroid\)\(x\_\{\\text\{centroid\}\},y\_\{\\text\{centroid\}\}\)
6:Compute softmax weights:

wi​j=exp⁡\(−di​j/τagg\)∑j′exp⁡\(−di​j′/τagg\)\+εw\_\{ij\}=\\dfrac\{\\exp\(\-d\_\{ij\}/\\tau\_\{\\text\{agg\}\}\)\}\{\\sum\_\{j^\{\\prime\}\}\\exp\(\-d\_\{ij^\{\\prime\}\}/\\tau\_\{\\text\{agg\}\}\)\+\\varepsilon\}for each

j∈N​\(i\)j\\in N\(i\)
7:endfor

8:// Step 2: Compute aggregated neighborhood representations

9:foreach cell

i∈Vi\\in Vdo

10:

ai←∑j∈N​\(i\)wi​j​xja\_\{i\}\\leftarrow\\sum\_\{j\\in N\(i\)\}w\_\{ij\}\\,x\_\{j\}\(distance\-weighted average of neighbor states\)

11:endfor

12:// Step 3: Initialize NIMfθf\_\{\\theta\}

13:Initialize: BatchNorm

→\\toLinear

\(p,512\)\(p,512\)→\\to3

×\\timesResidualBlock

\(512\)\(512\)→\\toTypeGate

\(64\)\(64\)→\\toLinear

\(512,p\)\(512,p\)
14:// Step 4: Train with early stopping

15:

θ∗←θ\\theta^\{\*\}\\leftarrow\\theta,

best\_val←∞\\text\{best\\\_val\}\\leftarrow\\infty,

patience←0\\text\{patience\}\\leftarrow 0
16:foreach epoch

e=1,2,…,100e=1,2,\\ldots,100do

17:foreach mini\-batch

ℬ⊆𝒞train\\mathcal\{B\}\\subseteq\\mathcal\{C\}\_\{\\text\{train\}\}do

18:

z^i←fθ​\(ai,τ​\(i\)\)\\hat\{z\}\_\{i\}\\leftarrow f\_\{\\theta\}\(a\_\{i\},\\tau\(i\)\)for each

i∈ℬi\\in\\mathcal\{B\}
19:

ℒ←1\|ℬ\|​∑i∈ℬℓδ​\(z^i,xi\)\\mathcal\{L\}\\leftarrow\\frac\{1\}\{\|\\mathcal\{B\}\|\}\\sum\_\{i\\in\\mathcal\{B\}\}\\ell\_\{\\delta\}\(\\hat\{z\}\_\{i\},x\_\{i\}\)\(Huber loss,δ=1\.0\\delta=1\.0\)

20:Update

θ\\thetavia AdamW \(

η=3×10−4\\eta=3\\times 10^\{\-4\},

λ=10−4\\lambda=10^\{\-4\}\), gradient clipping at norm

1\.01\.0
21:endfor

22:Compute validation loss on

𝒞val\\mathcal\{C\}\_\{\\text\{val\}\}
23:ifval\_loss

<<best\_valthen

24:

θ∗←θ\\theta^\{\*\}\\leftarrow\\theta,

best\_val←val\_loss\\text\{best\\\_val\}\\leftarrow\\text\{val\\\_loss\},

patience←0\\text\{patience\}\\leftarrow 0
25:else

26:

patience←patience\+1\\text\{patience\}\\leftarrow\\text\{patience\}\+1
27:endif

28:ifpatience

≥15\\geq 15then

29:break\(early stopping\)

30:endif

31:endfor

32:return

θ∗\\theta^\{\*\}

Algorithm 2Type\-Swap CDS1:Trained NIM

fθ∗f\_\{\\theta^\{\*\}\}; spatial graph

GG; features

XX; cell type labels

τ\\tau; sender type

ss; receiver type

rr; tissue cores

𝒞test\\mathcal\{C\}\_\{\\text\{test\}\}; distance bins

ℬ=\{\[0,10\),\[10,20\),\[20,30\),\[30,40\),\[40,∞\)\}\\mathcal\{B\}=\\\{\[0,10\),\[10,20\),\[20,30\),\[30,40\),\[40,\\infty\)\\\}; bootstrap iterations

B=1000B=1000; min receivers per core = 20\.

2:Type\-Swap

CDSS→R\\mathrm\{CDS\}\_\{S\\rightarrow R\};

CI95%\\mathrm\{CI\}\_\{95\\%\}\.

3:// Step 1: Identify valid receivers

4:

R∗←\{i∈V:τ​\(i\)=r​and​∃j∈N​\(i\)​s\.t\.​τ​\(j\)=s\}R^\{\*\}\\leftarrow\\\{i\\in V:\\tau\(i\)=r\\text\{ and \}\\exists\\,j\\in N\(i\)\\text\{ s\.t\. \}\\tau\(j\)=s\\\}
5:// Step 2: Type\-Swap counterfactual \+ CDS per receiver

6:foreach receiver cell

i∈R∗i\\in R^\{\*\}do

7:Compute

ai=∑j∈N​\(i\)wi​j​xja\_\{i\}=\\sum\_\{j\\in N\(i\)\}w\_\{ij\}x\_\{j\}, predict

z^i=fθ∗​\(ai,τ​\(i\)\)\\hat\{z\}\_\{i\}=f\_\{\\theta^\{\*\}\}\(a\_\{i\},\\tau\(i\)\)
8:// Construct counterfactual neighborhood

9:

Ncf​\(i\)←N​\(i\)N^\{\\mathrm\{cf\}\}\(i\)\\leftarrow N\(i\),

nreplaced←0n\_\{\\text\{replaced\}\}\\leftarrow 0
10:foreach

j∈N​\(i\)j\\in N\(i\)with

τ​\(j\)=s\\tau\(j\)=sdo

11:Let

b←b\\leftarrowdistance bin of

di​jd\_\{ij\}in

ℬ\\mathcal\{B\}
12:forexpansion

=0,1,2=0,1,2do\(search binbb, then widen\)

13:Candidates

←\\leftarrowcells in same core as

ii,

τ≠s\\tau\\neq s, within expanded bin, not already used

14:ifCandidates

≠∅\\neq\\emptysetthen

15:Pick

j′←j^\{\\prime\}\\leftarrowrandom cell from Candidates

16:

Ncf​\(i\)​\[j\]←j′N^\{\\mathrm\{cf\}\}\(i\)\[j\]\\leftarrow j^\{\\prime\},

nreplaced\+=1n\_\{\\text\{replaced\}\}\\mathrel\{\+\}=1
17:break

18:endif

19:endfor

20:endfor

21:if

nreplaced=0n\_\{\\text\{replaced\}\}=0then

22:skip\(no valid replacement found\)

23:endif

24:Compute

aicf=∑j∈Ncf​\(i\)wi​j​xja\_\{i\}^\{\\mathrm\{cf\}\}=\\sum\_\{j\\in N^\{\\mathrm\{cf\}\}\(i\)\}w\_\{ij\}x\_\{j\}\(distances unchanged\)

25:Predict

z^icf=fθ∗​\(aicf,τ​\(i\)\)\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}=f\_\{\\theta^\{\*\}\}\(a\_\{i\}^\{\\mathrm\{cf\}\},\\tau\(i\)\)
26:

CDSi←1p​‖z^icf−z^i‖1\\mathrm\{CDS\}\_\{i\}\\leftarrow\\frac\{1\}\{p\}\\\|\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i\}\\\|\_\{1\},

CDSisgn←1p​∑p\(z^i,pcf−z^i,p\)\\mathrm\{CDS\}\_\{i\}^\{\\mathrm\{sgn\}\}\\leftarrow\\frac\{1\}\{p\}\\sum\_\{p\}\(\\hat\{z\}\_\{i,p\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i,p\}\)
27:endfor

28:// Step 3: Filter cores with too few receivers

29:Remove from

R∗R^\{\*\}all cells from cores with

<20<20valid receivers

30:// Step 4: Core\-level bootstrap

31:for

b=1b=1to

BBdo

32:Sample

𝒞\(b\)\\mathcal\{C\}^\{\(b\)\}from

𝒞test\\mathcal\{C\}\_\{\\text\{test\}\}with replacement

33:

μ\(b\)←mean​\(\{CDSi:core​\(i\)∈𝒞\(b\)\}\)\\mu^\{\(b\)\}\\leftarrow\\mathrm\{mean\}\\bigl\(\\\{\\mathrm\{CDS\}\_\{i\}:\\mathrm\{core\}\(i\)\\in\\mathcal\{C\}^\{\(b\)\}\\\}\\bigr\)
34:endfor

35:// Step 5: Report

36:

CDSS→R←1\|R∗\|​∑i∈R∗CDSi\\mathrm\{CDS\}\_\{S\\rightarrow R\}\\leftarrow\\frac\{1\}\{\|R^\{\*\}\|\}\\sum\_\{i\\in R^\{\*\}\}\\mathrm\{CDS\}\_\{i\}
37:

CI95%←\[μ^2\.5,μ^97\.5\]\\mathrm\{CI\}\_\{95\\%\}\\leftarrow\[\\hat\{\\mu\}\_\{2\.5\},\\,\\hat\{\\mu\}\_\{97\.5\}\]
38:return

CDSS→R\\mathrm\{CDS\}\_\{S\\rightarrow R\},

CI95%\\mathrm\{CI\}\_\{95\\%\}

Algorithm 3Within\-Type CDS1:Trained NIM

fθ∗f\_\{\\theta^\{\*\}\}; spatial graph

GG; features

XX; cell type labels

τ\\tau; sender type

ss; receiver type

rr; tissue cores

𝒞test\\mathcal\{C\}\_\{\\text\{test\}\}\.

2:Within\-Type

CDSS→Rwt\\mathrm\{CDS\}\_\{S\\rightarrow R\}^\{\\mathrm\{wt\}\}\.

3:// Step 1: Identify valid receivers

4:

R∗←\{i∈V:τ​\(i\)=r​and​∃j∈N​\(i\)​s\.t\.​τ​\(j\)=s\}R^\{\*\}\\leftarrow\\\{i\\in V:\\tau\(i\)=r\\text\{ and \}\\exists\\,j\\in N\(i\)\\text\{ s\.t\. \}\\tau\(j\)=s\\\}
5:// Step 2: Within\-type counterfactual \+ CDS per receiver

6:foreach receiver cell

i∈R∗i\\in R^\{\*\}do

7:Compute

ai=∑j∈N​\(i\)wi​j​xja\_\{i\}=\\sum\_\{j\\in N\(i\)\}w\_\{ij\}x\_\{j\}, predict

z^i=fθ∗​\(ai,τ​\(i\)\)\\hat\{z\}\_\{i\}=f\_\{\\theta^\{\*\}\}\(a\_\{i\},\\tau\(i\)\)
8:// Build donor pool: same type, same core

9:

𝒟i←\{j′∈V:τ​\(j′\)=s,core​\(j′\)=core​\(i\),j′≠i\}\\mathcal\{D\}\_\{i\}\\leftarrow\\\{j^\{\\prime\}\\in V:\\tau\(j^\{\\prime\}\)=s,\\;\\mathrm\{core\}\(j^\{\\prime\}\)=\\mathrm\{core\}\(i\),\\;j^\{\\prime\}\\neq i\\\}
10:if

𝒟i=∅\\mathcal\{D\}\_\{i\}=\\emptysetthen

11:skip\(no same\-type donor in this core\)

12:endif

13:// Replace each sender slot with a random donor

14:

Ncf​\(i\)←N​\(i\)N^\{\\mathrm\{cf\}\}\(i\)\\leftarrow N\(i\)
15:foreach

j∈N​\(i\)j\\in N\(i\)with

τ​\(j\)=s\\tau\(j\)=sdo

16:

j′←j^\{\\prime\}\\leftarrowrandom draw from

𝒟i\\mathcal\{D\}\_\{i\}
17:

Ncf​\(i\)​\[j\]←j′N^\{\\mathrm\{cf\}\}\(i\)\[j\]\\leftarrow j^\{\\prime\}
18:endfor

19:Compute

aicf=∑j∈Ncf​\(i\)wi​j​xja\_\{i\}^\{\\mathrm\{cf\}\}=\\sum\_\{j\\in N^\{\\mathrm\{cf\}\}\(i\)\}w\_\{ij\}x\_\{j\}\(distances unchanged\)

20:Predict

z^icf=fθ∗​\(aicf,τ​\(i\)\)\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}=f\_\{\\theta^\{\*\}\}\(a\_\{i\}^\{\\mathrm\{cf\}\},\\tau\(i\)\)
21:

CDSiwt←1p​‖z^icf−z^i‖1\\mathrm\{CDS\}\_\{i\}^\{\\mathrm\{wt\}\}\\leftarrow\\frac\{1\}\{p\}\\\|\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i\}\\\|\_\{1\}
22:endfor

23:// Step 3: Aggregate

24:

CDSS→Rwt←1\|R∗\|​∑i∈R∗CDSiwt\\mathrm\{CDS\}\_\{S\\rightarrow R\}^\{\\mathrm\{wt\}\}\\leftarrow\\frac\{1\}\{\|R^\{\*\}\|\}\\sum\_\{i\\in R^\{\*\}\}\\mathrm\{CDS\}\_\{i\}^\{\\mathrm\{wt\}\}
25:return

CDSS→Rwt\\mathrm\{CDS\}\_\{S\\rightarrow R\}^\{\\mathrm\{wt\}\}

## Appendix BTheoretical Properties of Counterfactual Directionality

We next clarify what the proposed CDS measures and why the structured intervention design helps isolate directional influence\. The following results show that CDS admits a first\-order interpretation as a local intervention effect and that the proposed interventions disentangle neighborhood composition from graph\-structural artifacts\.

### B\.1CDS as a Local Intervention Effect

Letai∈ℝpa\_\{i\}\\in\\mathbb\{R\}^\{p\}denote the aggregated neighborhood representation of receiver cellii, and letaicfa\_\{i\}^\{\\mathrm\{cf\}\}denote the counterfactual representation obtained via a structured intervention\. Letfθ:ℝp×𝒯→ℝpf\_\{\\theta\}:\\mathbb\{R\}^\{p\}\\times\\mathcal\{T\}\\to\\mathbb\{R\}^\{p\}be continuously differentiable in its first argument\. Define

z^i=fθ​\(ai,τ​\(i\)\),z^icf=fθ​\(aicf,τ​\(i\)\)\.\\hat\{z\}\_\{i\}=f\_\{\\theta\}\(a\_\{i\},\\tau\(i\)\),\\qquad\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}=f\_\{\\theta\}\(a\_\{i\}^\{\\mathrm\{cf\}\},\\tau\(i\)\)\.
###### Theorem B\.1\(First\-Order Interpretation of CDS\)\.

There exists a pointξi\\xi\_\{i\}on the line segment betweenaia\_\{i\}andaicfa\_\{i\}^\{\\mathrm\{cf\}\}such that

z^icf−z^i=Jf​\(ξi,τ​\(i\)\)​\(aicf−ai\),\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i\}=J\_\{f\}\(\\xi\_\{i\},\\tau\(i\)\)\\,\(a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\),whereJf​\(ξi,τ​\(i\)\)J\_\{f\}\(\\xi\_\{i\},\\tau\(i\)\)denotes the Jacobian offθ​\(⋅,τ​\(i\)\)f\_\{\\theta\}\(\\cdot,\\tau\(i\)\)evaluated atξi\\xi\_\{i\}\. Consequently,

CDSi=1p​‖Jf​\(ξi,τ​\(i\)\)​\(aicf−ai\)‖1\.\\mathrm\{CDS\}\_\{i\}=\\frac\{1\}\{p\}\\left\\\|J\_\{f\}\(\\xi\_\{i\},\\tau\(i\)\)\\,\(a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\)\\right\\\|\_\{1\}\.If‖aicf−ai‖2\\\|a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\\\|\_\{2\}is small, then

CDSi=1p​‖Jf​\(ai,τ​\(i\)\)​\(aicf−ai\)‖1\+O​\(‖aicf−ai‖22\)\.\\mathrm\{CDS\}\_\{i\}=\\frac\{1\}\{p\}\\left\\\|J\_\{f\}\(a\_\{i\},\\tau\(i\)\)\\,\(a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\)\\right\\\|\_\{1\}\+O\\\!\\left\(\\\|a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\\\|\_\{2\}^\{2\}\\right\)\.

This result shows that CDS measures the magnitude of a local, first\-order response of the learned model to perturbations in neighborhood composition\.

### B\.2Structural Invariance of Interventions

Recall that the aggregated neighborhood representation is defined as

ai=∑j∈N​\(i\)wi​j​xj,a\_\{i\}=\\sum\_\{j\\in N\(i\)\}w\_\{ij\}x\_\{j\},with weightswi​jw\_\{ij\}determined by spatial distances\.

###### Theorem B\.2\(Invariance under Structured Interventions\)\.

Assume that the counterfactual intervention preserves: \(i\) the number of neighborhood slotskk, \(ii\) the spatial distances associated with each slot \(and hence weightswi​jw\_\{ij\}\), and \(iii\) the receiver identityτ​\(i\)\\tau\(i\)\.

Then the perturbation satisfies

aicf−ai=∑m∈ℛiwi​m​\(xmcf−xm\),a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}=\\sum\_\{m\\in\\mathcal\{R\}\_\{i\}\}w\_\{im\}\\left\(x\_\{m\}^\{\\mathrm\{cf\}\}\-x\_\{m\}\\right\),whereℛi\\mathcal\{R\}\_\{i\}is the set of replaced neighborhood slots\. In particular, the perturbation depends only on the change in feature vectors occupying fixed spatial slots and is independent of graph degree, spatial weighting, or receiver identity\.

Consequently,z^icf−z^i\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i\}depends only on compositional changes induced by the structured intervention design underfθf\_\{\\theta\}\.

This property ensures that CDS captures compositional effects rather than artifacts arising from changes in graph topology or weighting\.

### B\.3Null Behavior under Within\-Type Exchangeability

###### Lemma B\.3\(Exchangeability Null\)\.

Assume that, conditional on tissue core and sender type, sender cells are exchangeable with respect to their contribution to the receiver prediction\. Then, for within\-type interventions,

𝔼​\[aicf−ai\]=0,\\mathbb\{E\}\\\!\\left\[a\_\{i\}^\{\\mathrm\{cf\}\}\-a\_\{i\}\\right\]=0,and, to first order,

𝔼​\[z^icf−z^i\]≈0\.\\mathbb\{E\}\\\!\\left\[\\hat\{z\}\_\{i\}^\{\\mathrm\{cf\}\}\-\\hat\{z\}\_\{i\}\\right\]\\approx 0\.Thus, the expected within\-type CDS reflects baseline variability due to finite sampling and sender\-state heterogeneity rather than systematic sender\-type effects\.

This result justifies the use of within\-type interventions as a baseline for interpreting type\-swap CDS values\.

### B\.4Approximation Bias from Neighborhood Aggregation

###### Lemma B\.4\(Aggregation Bias\)\.

Suppose the true receiver response depends on the full neighborhood configuration

zi⋆=g​\(\{\(wi​j,xj\)\}j∈N​\(i\),τ​\(i\)\),z\_\{i\}^\{\\star\}=g\\big\(\\\{\(w\_\{ij\},x\_\{j\}\)\\\}\_\{j\\in N\(i\)\},\\tau\(i\)\\big\),while the model uses only the aggregated representation

ai=∑j∈N​\(i\)wi​j​xj\.a\_\{i\}=\\sum\_\{j\\in N\(i\)\}w\_\{ij\}x\_\{j\}\.Ifggcannot be expressed solely as a function ofaia\_\{i\}, then no functionffexists such that

g=f​\(ai,τ​\(i\)\)g=f\(a\_\{i\},\\tau\(i\)\)for all neighborhoods\. Consequently, CDS may be biased relative to intervention effects defined on the full neighborhood\.

This highlights that CDS captures effects mediated through the chosen neighborhood summary and may not reflect higher\-order interactions among neighbors\.

### B\.5Statistical Inference via Core\-Level Bootstrap

Cells within the same tissue core are not statistically independent: they share a common microenvironment, patient\-level covariates, and technical batch effects\. Naively treating each cell as an independent observation would severely underestimate uncertainty\. We therefore perform statistical inference at the*tissue\-core level*using a nonparametric bootstrap\.

Let𝒞=\{c1,c2,…,cM\}\\mathcal\{C\}=\\\{c\_\{1\},c\_\{2\},\\ldots,c\_\{M\}\\\}denote the set ofMMtest tissue cores\. The bootstrap procedure is as follows: for each iterationb=1,…,Bb=1,\\ldots,B\(withB=1000B=1000\):

1. 1\.SampleMMcores*with replacement*from𝒞\\mathcal\{C\}to obtain a bootstrap core set𝒞\(b\)\\mathcal\{C\}^\{\(b\)\}\. Some cores may appear multiple times; some may not appear at all\.
2. 2\.Pool all cell\-level CDS values from cores in𝒞\(b\)\\mathcal\{C\}^\{\(b\)\}:𝒱\(b\)=\{CDSi:core​\(i\)∈𝒞\(b\)\}\\mathcal\{V\}^\{\(b\)\}=\\bigl\\\{\\mathrm\{CDS\}\_\{i\}:\\mathrm\{core\}\(i\)\\in\\mathcal\{C\}^\{\(b\)\}\\bigr\\\}\.
3. 3\.Compute the bootstrap mean:μ\(b\)=\|𝒱\(b\)\|−1​∑v∈𝒱\(b\)v\\mu^\{\(b\)\}=\|\\mathcal\{V\}^\{\(b\)\}\|^\{\-1\}\\sum\_\{v\\in\\mathcal\{V\}^\{\(b\)\}\}v\.

The resulting distribution\{μ\(b\)\}b=1B\\\{\\mu^\{\(b\)\}\\\}\_\{b=1\}^\{B\}approximates the sampling distribution ofCDSS→R\\mathrm\{CDS\}\_\{S\\rightarrow R\}under core\-level resampling\. Confidence intervals are computed via the*percentile method*:

CI95%=\[μ^2\.5,μ^97\.5\],\\mathrm\{CI\}\_\{95\\%\}=\\Bigl\[\\,\\hat\{\\mu\}\_\{2\.5\},\\;\\hat\{\\mu\}\_\{97\.5\}\\,\\Bigr\],\(16\)whereμ^α\\hat\{\\mu\}\_\{\\alpha\}denotes theα\\alpha\-th percentile of\{μ\(b\)\}\\\{\\mu^\{\(b\)\}\\\}\. A sender–receiver interaction is considered statistically meaningful if \(i\)CI95%\\mathrm\{CI\}\_\{95\\%\}for the signed CDS excludes zero, or \(ii\) the magnitude CDS is substantially elevated relative to the within\-type baseline, which quantifies the expected CDS under same\-type reshuffling\.