Causal Integrated Gradients: Scalable Feature Attribution Along Causal Paths

Joran Michiels; Anis Ismail; Maarten De Vos; Alejandro Sifrim

cs.LGOct 15, 2026

Causal Integrated Gradients: Scalable Feature Attribution Along Causal Paths

Authors: Joran Michiels, Anis Ismail, Maarten De Vos, Alejandro Sifrim

Organizations: Laboratory of Multi-omics Integrative Bioinformatics, Department of Human Genetics, KU Leuven, Leuven, Belgium · KU Leuven Institute for Single Cell Omics (LISCO), KU Leuven, Leuven, Belgium · Department of Electrical Engineering, KU Leuven, Leuven, Belgium · Department of Development and Regeneration, KU Leuven, Leuven, Belgium

Abstract

Existing feature attribution methods that account for feature dependencies do not scale to high dimensions. We propose causal integrated gradients (CIG) and causal expected gradients (CEG), which extend integrated gradients to incorporate prior causal structure among input features while retaining its computational efficiency and, under conditions on the causal graph, its axiomatic guarantees. In many settings, some causal structure between features is already known, for example, that a transcription factor causally regulates a target gene's expression. Yet existing gradient-based methods cannot exploit this structure, instead crediting a feature's downstream effects over its true upstream cause, such as ranking the target gene above the transcription factor that regulates it. Our key idea is to embed this causal structure by routing the integration path along the causal graph. The method thereby decomposes each feature's attribution into a direct effect and effects mediated through its causal descendants. We evaluate CIG in three settings of increasing realism. First, on a synthetic benchmark with a known causal model, CIG's attributions best match the ground-truth causal effect from Pearl's do-calculus, matching the accuracy of causal Shapley values at a fraction of the computational cost. Second, on simulated single-cell gene expression data with a known gene regulatory network (Dyngen), CIG most accurately recovers the transcription factors driving the simulation. Third, on real interferon- $β$ stimulation data, where only an approximate causal graph is available, CIG better prioritizes the transcription factors driving the measured response than DeepSHAP and integrated gradients. Together, these results show that incorporating causal structure into gradient-based attributions recovers upstream drivers that existing methods miss, without sacrificing scalability.

Figures & tables

Panel of Figure 1 — Figure 1: Simulated example inspired by attribution of gene regulation. Existing graph-agnostic methods assign the highest attribution scores to the directly correlated target genes, missing the true upstream regulators; our graph-aware method correctly pushes the transcription factors to the top.

Panel of Figure 2 — Figure 2: Chain rule diagrams illustrating the reparametrization of path functions to incorporate causal dependencies. Left: general $n$ -feature case. Right: two-feature example. Black edges are standard path method dependencies; blue edges are the causal dependencies introduced by our reparametrization.

Panel of Figure 3 — Figure 3: Left: Spearman correlation with the CATE on root nodes only. Right: wall-clock runtime (log–log scale) to explain all 100 samples; the dip for the Shapley-based methods reflects their reduced sampling budget at larger feature counts. Both as a function of the number of features. Shaded areas denote the standard deviation over 5 runs. Results on all nodes are shown in Figure F.1 .

Figure 4: Evaluation of feature attribution methods on Dyngen. (A) Area under the recall curve for backbone transcription factors, normalized so that a random ranking scores 0 and the optimal ranking 100. (B) DCG over the top-

K

ranked genes, with the DCG of the optimal ranking shown for reference; relevance scores follow Table G.6 . Error bars denote the standard deviation over 5 seeds.

Figure 5: Recovery of stimulation-responsive regulators on the IFN-

\beta

dataset. (A) Distribution of attribution ranks for genes grouped by hop distance from the roots of the SCENIC validation graph, median indicated by a horizontal line within each box. (B)

z

-score of the Mann–Whitney

U

test statistic per method and cell type, averaged across 5 seeds; higher values indicate responsive TFs are ranked higher than other genes; error bars denote standard deviation across seeds;

\ast

: FDR

<0.05

for every seed.

Appendix figures & tables42 assets

Supplementary material from the paper’s appendix.

Appendix

Setting	CIG	IG
Synthetic, 100 features	$1.90\times 10^{-3}\pm 1.16\times 10^{-3}$	$5.77\times 10^{-8}\pm 6.45\times 10^{-9}$
Synthetic, 129 features	$1.40\times 10^{-3}\pm 3.78\times 10^{-4}$	$5.53\times 10^{-8}\pm 2.34\times 10^{-9}$
Synthetic, 166 features	$1.38\times 10^{-3}\pm 2.99\times 10^{-4}$	$5.96\times 10^{-8}\pm 4.62\times 10^{-9}$
Synthetic, 215 features	$1.60\times 10^{-3}\pm 4.11\times 10^{-4}$	$6.56\times 10^{-8}\pm 4.21\times 10^{-9}$
Synthetic, 278 features	$1.75\times 10^{-3}\pm 5.08\times 10^{-4}$	$6.59\times 10^{-8}\pm 3.51\times 10^{-9}$
Synthetic, 359 features	$1.96\times 10^{-3}\pm 3.64\times 10^{-4}$	$7.09\times 10^{-8}\pm 5.30\times 10^{-9}$

Appendix

Table C.1 : Relative completeness delta (WAPE), mean

\pm

std over seeds.

Figure C.1 : Per-sample absolute completeness delta $|\delta^{(n)}|$ at 774 features ( $n=500$ , pooled across 5 seeds). Median: CIG $1.48\times 10^{-2}$ , IG $5.74\times 10^{-7}$ . Proportional spread (max/min): CIG $1.71\times 10^{3}\times$ , IG $2.77\times 10^{3}\times$ — the same order of magnitude.

Figure C.2 : Per-sample absolute completeness delta

|\delta^{(n)}|

for the Dyngen dataset (

n=14{,}076

, pooled across 5 seeds).

Figure C.3 : Per-sample absolute completeness delta

|\delta^{(n)}|

for the IFN-

\beta

dataset (6,941 features;

n=25{,}450

, pooled across 5 seeds). Median: CIG

6.75\times 10^{-2}

, IG

3.34\times 10^{-3}

. Proportional spread (max/min): CIG

9.05\times 10^{5}\times

, IG

2.11\times 10^{5}\times

— the same order of magnitude.

Figure C.4 : Relative completeness delta of CIG as a function of the number of integration steps $K$ , at 100 synthetic features (log-log axes). The shaded band shows the standard deviation over seeds. The linear trend on log-log axes is the $O(1/K)$ signature of first-order discretization error.

Method	Cost
IG (vanilla)	$O(K\cdot d\cdot C)$
CIG (ours)	$O(K\cdot d\cdot C)$ for a sparse causal graph ( $E{=}O(d)$ ); $O(K\cdot d^{2})$ for a dense graph (worst case)
Causal Shapley values (sampled, as implemented)	$O(M\cdot d^{2})$
Conditional SHAP (sampled, as implemented)	$O(M\cdot d^{3})$ (one-time $O(d^{3})$ Cholesky)
Exact Shapley value (causal or conditional)	$O(2^{d})$

Appendix

Table D.1 : Complexity per explained instance.

Dataset	Method	Stage	Time (s)	Peak CUDA alloc. (MB)	Peak CUDA reserved (MB)
Dyngen	CIG (linear basis, 1st order)	Surrogate training	129.7 $\pm$ 4.2 ( $\sim$ 130 s)	224 $\pm$ 7	980 $\pm$ 958
		Attribution computation	175.4 $\pm$ 19.9 ( $\sim$ 175 s)	17818 $\pm$ 927	18680 $\pm$ 1333
		Total	305.1 $\pm$ 23.6 ( $\sim$ 305 s)	17818 $\pm$ 927	18680 $\pm$ 1333
Dyngen	CIG (linear basis, 2nd order)	Surrogate training	130.4 $\pm$ 5.0 ( $\sim$ 130 s)	221 $\pm$ 12	966 $\pm$ 970
		Attribution computation	349.6 $\pm$ 39.6 ( $\sim$ 350 s)	26653 $\pm$ 1416	27500 $\pm$ 2004
		Total	480.0 $\pm$ 44.6 ( $\sim$ 480 s)	26653 $\pm$ 1416	27500 $\pm$ 2004

Appendix

Table E.1 : End-to-end pipeline profiling: wall-clock time and peak CUDA memory per stage, for the full test set at

K{=}100

integration steps. Dyngen benchmark (

2{,}514

–

3{,}000

test cells per seed); IFN-

\beta

(

5{,}090

cells

\times

6{,}941

genes). Values are mean

\pm

std over 5 seeds. CEG averages over

100

random baselines; all CIG rows use a zero baseline.

Figure F.1 : Spearman correlation with the CATE averaged over all nodes, as a function of the number of features. Shaded areas denote the standard deviation over 5 runs.

Panel of Figure 14 — Figure F.2 : Spearman correlation with ground-truth CATE as a function of the number of features: CIG using the SCM’s true mechanisms directly (oracle SCM surrogate) vs. CIG using our learned, expanded-basis surrogate (labeled “GAM” in the legend) vs. graph-agnostic integrated gradients. Shaded areas denote the standard deviation over 5 seeds.

	Cells		Genes		Filtered cells per split
Seed	Raw	Filtered	Raw	Filtered	Train	Val.	Test
42	10,000	9,790	2,548	2,378	4,797	2,056	2,937
420	10,000	8,380	2,548	2,445	4,106	1,760	2,514
100	10,000	8,871	2,548	2,320	4,346	1,863	2,662
1000	10,000	10,000	2,548	2,388	4,900	2,100	3,000
2000	10,000	9,876	2,548	2,494	4,839	2,074	2,963

Appendix

Table G.2 : Dyngen regression datasets across the five simulation seeds. Cells and genes are reported before and after quality-control filtering; the train/validation/test split is a nested 70/30 split of the filtered cells (49%/21%/30%).

Table 17

Method	Used	Baseline	Notes
DeepSHAP [ Lundberg and Lee, 2017 ]	$\checkmark$	Random set	captum implementation
GradientSHAP [ Lundberg and Lee, 2017 ]	$\checkmark$	Random set	captum implementation
Integrated gradients [ Sundararajan et al., 2017 ]	$\checkmark$	Zero	captum implementation
DeepLIFT [ Shrikumar et al., 2017 ]	$\checkmark$	Zero	captum implementation
Saliency [ Simonyan et al., 2014 ]	$\checkmark$	$\times$	captum implementation
Guided backpropagation [ Springenberg et al., 2015 ]	$\checkmark$	$\times$	captum implementation

Appendix

Table G.5 : Feature attribution methods evaluated in the simulated single-cell experiment.

Gene Category	Score	Description
Backbone transcription factors (TFs)	+1.00	TFs defining main differentiation trajectory in GRN
Direct targets of backbone TFs	+0.25	Genes regulated by one backbone TF
Joint targets with opposing interaction signs	+0.25	Genes jointly regulated by two TFs with opposite regulatory signs
Burn-in phase TFs	-0.25	TFs active only during burn-in phase
Targets of burn-in TFs	-0.25	Genes regulated by burn-in phase TFs
Housekeeping genes	-0.25	Non-informative constitutively expressed genes

Appendix

Table G.6 : Gene relevance scoring scheme based on the simulated gene regulatory network (GRN). The DCG at cutoff

K

\mathrm{DCG}@K=\sum_{i=1}^{K}\mathrm{rel}_{i}/\log_{2}(i+1)

, where

\mathrm{rel}_{i}

is the score of the gene at rank

i

Table 20

Min. seeds	1	2	4	6	8	10
Edges	74,426	58,082	42,448	33,172	26,379	19,055
TFs	252	229	193	174	156	131

Appendix

Table H.3 : Edges kept at each recurrence threshold (full data).

Figure I.1 : Distribution of attribution ranks for genes grouped by their hop distance from the roots of the OmniPath prior graph. Each boxplot represents the distribution of ranks for genes within a specific hop-distance group, with the median rank indicated by a horizontal line within each box.

	Roots	1-hop	2-hop	3-hop	4-hop	Isolates
CIG	0.1220	0.1085	0.1126	0.1197	0.1445	0.0959
DeepSHAP	0.1312	0.1194	0.1292	0.1335	0.1598	0.1428
Integrated gradients	0.0980	0.0891	0.0956	0.0987	0.1065	0.1017

Appendix

Table I.1 : Median standard deviations of percentile rank across seeds (OmniPath); chance

=0.278

, perfectly stable

=0.0

	Roots	1-hop	2-hop	3-hop	4-hop	5-hop	Isolates
CIG	0.1057	0.1026	0.1007	0.1041	0.1145	0.0872	0.1095
DeepSHAP	0.1940	0.1326	0.1242	0.1295	0.1502	0.1323	0.1535
Integrated gradients	0.1309	0.0974	0.0929	0.0969	0.1085	0.0535	0.1026

Appendix

Table I.2 : Median standard deviations of percentile rank across seeds (SCENIC); chance

=0.278

, perfectly stable

=0.0

Figure I.2 : Top-

k

enrichment of stimulation-responsive SCENIC regulons. The

z

-score of the hypergeometric test statistic is shown for each method, cell type, and top-

k

cutoff, averaged across five random seeds. Higher values indicate that responsive TFs are over-represented among the top-

k

genes. Error bars denote the standard deviation across seeds.

P

-values were BH-corrected across cell types and

k

. Only significant enrichment bars are shaded.

Cell type	Method	Mean Spearman $\rho$	S.d. across seeds
B cells	CIG	0.3098	0.0072
B cells	DeepSHAP	0.3724	0.0015
B cells	Integrated gradients	0.3465	0.0039
CD14 + monocytes	CIG	0.3530	0.0343
CD14 + monocytes	DeepSHAP	0.3945	0.0215
CD14 + monocytes	Integrated gradients	0.4119	0.0395

Appendix

Table I.3 : Whole-panel agreement between attribution and differential expression. For each cell type and method, the Spearman correlation between the absolute attribution difference

\delta_{g}

(stimulated vs. control) and the absolute DESeq2 Wald statistic is computed over the full measured gene panel. Values are the mean and standard deviation across five random seeds.

F1-score at $k=n_{\text{DEGs}}$
Cell type	CIG	DeepSHAP	Integrated gradients
B cells	$47.12\pm 0.48$	$57.56\pm 0.44$	$54.78\pm 0.40$
CD14 + monocytes	$43.70\pm 1.56$	$46.41\pm 0.83$	$46.78\pm 1.72$
CD4 T cells	$35.56\pm 1.07$	$38.10\pm 1.06$	$39.39\pm 1.06$
CD8 T cells	$49.20\pm 1.90$	$67.15\pm 1.55$	$65.55\pm 0.80$
Dendritic cells	$35.55\pm 0.83$	$36.19\pm 0.93$	$35.52\pm 0.78$

Appendix

Table I.4 : Recovery of pseudobulk DEGs at a cutoff equal to the number of DEGs (

k=n_{\text{DEGs}}

). For each cell type, genes are ranked by the absolute attribution difference

\delta_{g}

and scored against that cell type’s DESeq2 DEG set. The recovery AUC is normalized so that a random ranking scores

0

and the optimal ranking scores

1

; NDCG is graded by the absolute Wald statistic. At this cutoff, precision, recall, and F1 coincide by construction, so only F1 is reported. All values are percentages, given as mean

\pm

s.d. across five random seeds. Megakaryocytes are omitted as they contributed no significant DEGs.

Precision
Cell type	CIG	DeepSHAP	Integrated gradients
B cells	$52.00\pm 13.04$	$100.00\pm 0.00$	$100.00\pm 0.00$
CD14 + monocytes	$48.00\pm 8.37$	$100.00\pm 0.00$	$100.00\pm 0.00$
CD4 T cells	$34.00\pm 13.42$	$96.00\pm 5.48$	$96.00\pm 5.48$
CD8 T cells	$40.00\pm 12.25$	$100.00\pm 0.00$	$100.00\pm 0.00$
Dendritic cells	$48.00\pm 13.04$	$100.00\pm 0.00$	$96.00\pm 5.48$

Appendix

Table I.5 : Recovery of pseudobulk DEGs per cell type at a top-

k

cutoff of

k=10

. Genes are ranked by the absolute attribution difference

\delta_{g}

and scored against that cell type’s DESeq2 DEG set. NDCG is graded by the absolute Wald statistic. All values are percentages, given as mean

\pm

s.d. across five random seeds. Megakaryocytes are omitted as they contributed no significant DEGs.

Precision
Cell type	CIG	DeepSHAP	Integrated gradients
B cells	$56.80\pm 12.46$	$96.00\pm 0.00$	$96.00\pm 0.00$
CD14 + monocytes	$60.80\pm 11.45$	$100.00\pm 0.00$	$97.60\pm 3.58$
CD4 T cells	$49.60\pm 8.29$	$96.00\pm 0.00$	$96.00\pm 0.00$
CD8 T cells	$48.80\pm 5.93$	$97.60\pm 2.19$	$94.40\pm 3.58$
Dendritic cells	$50.40\pm 8.29$	$98.40\pm 2.19$	$93.60\pm 6.07$

Appendix

Table I.6 : Recovery of pseudobulk DEGs per cell type at a top-

k

cutoff of

k=25

. Metrics and conventions are as in Table I.5 .

Precision
Cell type	CIG	DeepSHAP	Integrated gradients
B cells	$62.80\pm 4.82$	$98.00\pm 0.00$	$95.60\pm 2.19$
CD14 + monocytes	$70.00\pm 8.60$	$98.80\pm 1.10$	$95.20\pm 3.35$
CD4 T cells	$59.60\pm 5.55$	$96.00\pm 0.00$	$96.00\pm 0.00$
CD8 T cells	$54.40\pm 9.84$	$94.80\pm 1.79$	$93.20\pm 2.28$
Dendritic cells	$61.20\pm 4.60$	$97.20\pm 2.28$	$94.40\pm 5.37$

Appendix

Table I.7 : Recovery of pseudobulk DEGs per cell type at a top-

k

cutoff of

k=50

. Metrics and conventions are as in Table I.5 .

Precision
Cell type	CIG	DeepSHAP	Integrated gradients
B cells	$64.00\pm 1.87$	$93.80\pm 2.68$	$90.40\pm 2.97$
CD14 + monocytes	$74.60\pm 6.23$	$95.80\pm 1.30$	$92.60\pm 4.39$
CD4 T cells	$61.20\pm 4.60$	$91.00\pm 0.71$	$89.20\pm 1.79$
CD8 T cells	$57.00\pm 3.61$	$78.40\pm 2.30$	$75.40\pm 1.82$
Dendritic cells	$67.40\pm 2.19$	$92.40\pm 3.65$	$89.40\pm 4.16$

Appendix

Table I.8 : Recovery of pseudobulk DEGs per cell type at a top-

k

cutoff of

k=100

. Metrics and conventions are as in Table I.5 .

Precision
Cell type	CIG	DeepSHAP	Integrated gradients
B cells	$56.48\pm 1.43$	$76.64\pm 1.15$	$73.28\pm 1.37$
CD14 + monocytes	$71.28\pm 4.21$	$84.88\pm 1.80$	$81.60\pm 3.97$
CD4 T cells	$52.64\pm 1.00$	$66.72\pm 1.58$	$66.00\pm 1.52$
CD8 T cells	$37.60\pm 1.96$	$45.84\pm 0.83$	$45.28\pm 0.87$
Dendritic cells	$64.08\pm 2.05$	$79.28\pm 1.18$	$75.76\pm 2.31$

Appendix

Table I.9 : Recovery of pseudobulk DEGs per cell type at a top-

k

cutoff of

k=250

. Metrics and conventions are as in Table I.5 .

Term	Adjusted P-value	Odds Ratio	Combined Score	Genes
TNF-alpha Signaling via NF-kB	0.0003	14.1062	156.8329	JUN, KLF6, IRF1, MYC, FOS, NFE2L2
Hypoxia	0.0003	17.0927	178.5641	JUN, KLF6, FOS, NR3C1, ETS1
Interferon Gamma Response	0.0003	12.1222	124.9307	IRF4, IRF1, STAT2, IRF2, IRF8, HIF1A
G2-M Checkpoint	0.0041	11.7922	86.3737	MYC, CTCF, E2F4, HIF1A
IL-2/STAT5 Signaling	0.0041	10.6567	74.2016	KLF6, IRF4, MYC, IRF8
Inflammatory Response	0.0041	10.5641	73.2283	KLF6, IRF1, MYC, HIF1A

Appendix

Table I.10: Over-representation analysis of the top-100 genes unique to CIG in dendritic cells. Terms are gene sets from the MSigDB Hallmark 2020 collection, and the measured gene panel was used as the background set. The ten most significant terms are listed.

P

-values come from a hypergeometric test and are Benjamini–Hochberg adjusted across the collection; the combined score is

\ln(p)

multiplied by the

z

-score of the deviation from the expected rank. The last column gives the genes shared between the top-100 set and the term. Adjusted

P