A physical wiring diagram for the human immune system

Abstract

The human immune system is composed of a distributed network of cells circulating throughout the body, which must dynamically form physical associations and communicate using interactions between their cell-surface proteomes¹. Despite their therapeutic potential², our map of these surface interactions remains incomplete^3,4. Here, using a high-throughput surface receptor screening method, we systematically mapped the direct protein interactions across a recombinant library that encompasses most of the surface proteins that are detectable on human leukocytes. We independently validated and determined the biophysical parameters of each novel interaction, resulting in a high-confidence and quantitative view of the receptor wiring that connects human immune cells. By integrating our interactome with expression data, we identified trends in the dynamics of immune interactions and constructed a reductionist mathematical model that predicts cellular connectivity from basic principles. We also developed an interactive multi-tissue single-cell atlas that infers immune interactions throughout the body, revealing potential functional contexts for new interactions and hubs in multicellular networks. Finally, we combined targeted protein stimulation of human leukocytes with multiplex high-content microscopy to link our receptor interactions to functional roles, in terms of both modulating immune responses and maintaining normal patterns of intercellular associations. Together, our work provides a systematic perspective on the intercellular wiring of the human immune system that extends from systems-level principles of immune cell connectivity down to mechanistic characterization of individual receptors, which could offer opportunities for therapeutic intervention.

Fig. 1. A leukocyte receptor network by systematic protein interaction mapping.

a, SAVEXIS enables efficient and high-throughput screening for protein binding interactions between recombinant extracellular domains. b, Schematic showing the diverse structural architectures of leukocyte surface proteins within the pan-leukocyte library of 630 proteins. The number of proteins from each class is noted above, and the recombinant expression strategy is illustrated below. c, Summarized matrix of protein–protein pairs for immune receptors with interactions either identified by screening or previously reported in the literature. The average signal intensity for a given bait–prey measurement orientation across the primary and secondary screens is indicated by the shaded intensity, and the colour indicates which interactions are novel. d, Screening successfully finds most previously reported interactions with minimal false positives. Receiver operating characteristic (ROC) curve for average measurements of protein–protein pairs against reference sets of expected positive and randomized negative interactions. AUC, area under the curve. e, Organized interaction network of immune receptor interactions. The colour indicates which interactions are novel, and the line thickness is proportional to the magnitude of evidence from the screening measurements. Source data.

Fig. 2. Validating and assembling a quantitative immune interactome.

a, Novel receptor interactions are detectable on the surfaces of live human cells. For six examples that encompass different architectural classes, flow cytometry traces are shown for the binding of fluorescent-conjugated protein to HEK293 cells overexpressing its identified counter-receptor (blue) or control cells (red). b, SPR substantiates and quantifies the binding of novel leukocyte receptor–ligand pairs. For the same six example protein pairs, sensorgram data (black) are shown with Langmuir model fitting curves overlaid (red) for all interactions for which a robust fit could be calculated. Ig-SF, immunoglobulin superfamily; LRR-SF, leucine rich repeat superfamily. c, The quantitative interactome of immune cell-surface proteins. Proteins are shown as circular charts indicating the proportion of expression in each leukocyte population. Binding affinity between proteins is indicated by the size and intensity of red edges (expressed in terms of the binding dissociation constant (K_D), where smaller values reflect stronger binding). Abbreviations for cell type names are defined in Supplementary Table 5. d, Immune cell subsets use related but varying distributions of binding affinities when communicating with other cell types. For each pairing of two cells, a histogram of inferred interactions is shown alongside a colour shade that indicates the average affinity. e, Inflammatory activation broadly reconfigures receptors towards those with less-transient binding kinetics. After differential expression testing of surface proteins between activated and stimulated leukocytes (n = 4 samples per condition), the binding affinities of interactions involving downregulated (downreg) receptors are compared to the binding affinities of upregulated (upreg) receptors. Data are shown as Tukey box plots with Holm-corrected P values from a two-sided Welch’s t-test. f, Intercellular connectivity can be mathematically predicted by integrating protein expression, binding affinity and cell parameters using physics-based equations. A detailed description of the model can be found in the Supplementary Equations. g, Model predictions for baseline rates of immune cell association agree with published data measuring in vitro immune cell association. Each data point has two colours that correspond to the two physically interacting cell types. Shading depicts the 95% compatibility interval for the least-squares linear regression fit. Source data.

Fig. 3. An interactive atlas of immune cell connections across the human body.

a, Systematic integration of single-cell datasets to map cellular connectivity across tissues with substantial immune populations. Cell types are positioned around each circle, with each position along it marking a cell-surface protein. Linkages formed by physically interacting surface proteins between cell types are marked by curved lines, coloured by interaction identity. Full abbreviations are listed in Supplementary Table 5. b, Functionalities available through our interactive atlas of physical immune interactions (https://www.sanger.ac.uk/tool/immune-interaction/immune-interaction). c, Myeloid cells act as interaction hubs in immune tissues. Eigenvector centrality metrics of myeloid cells compared to all other populations after converting the total interaction count for all cell–cell pairs into a weighted undirected graph. Data are shown as Tukey box plots with Benjamini–Hochberg P values calculated from a two-sided Welch’s t-test. d, Spatial transcriptomics of a human lymph node confirms that our identified interaction partners are physically colocalized in situ. An example data point of a tissue section analysed for the JAG1 + VASN interaction is shown. The percentage of measured spots in which the expression of one protein of an interacting pair is spatially connected to the other protein of the pair is compared for previously reported interactions (green), novel interactions (blue) and a negative control of the same proteins with interaction links randomly permuted (yellow) (n = 100). Data are shown as Tukey box plots with P values calculated from a Tukey's honest significance test. e, Single-molecule RNA hybridization on human lymph nodes defines regions in which newly identified interaction partners are expressed in spatially bordering cells. A single lymphoid follicle enriched in CD45⁺ leukocytes is magnified (left), showing the zonation of JAG1- and VASN-expressing cells into the corona and the germinal centre, respectively (middle). An inset (right) highlights a region of bordering cells expressing each marker. Scale bars, 200 μm (left); 100 μm (middle); 50 μm (right). Source data.

Fig. 4. Multiplex leukocyte assays identify functional pathways for receptor proteins.

a, High-content microscopy set-up for perturbing human peripheral blood mononuclear cells (PBMCs) with recombinant proteins and measuring changes in cellular activation and connectivity. Scale bar, 30 µm. b, Proteins with identified receptor interactions elicit responses on lymphocyte action. Polarization of lymphocyte populations in resting and weakly activating background conditions (y axis) after addition of soluble protein extracellular domains (x axis). Stimulation (red) or inhibition (blue) relative to control is shown of a cell polarization marker of lymphocyte activation. Phenotypes that have P values below the adjusted significance threshold are outlined in bold. n = 10 samples. c, Interacting cellular communities can be extracted from high-content imaging data. Representative microscopy fields (left) and computed physical cell contacts (right, white lines) are depicted for leukocytes perturbed with recombinant SEMA4D and SIRPA as examples. Scale bar, 30 µm. d, Rewiring of cellular interactions by perturbing receptor pathways. Measured changes in cell–cell interactions (x axis) induced by recombinant proteins (panels) across measurement time points and background conditions (y axis). The same colour scale as in c is used to identify cell pairs along the x axis. n = 10 samples. e, Observed interaction changes conform to mathematical model predictions. Average magnitudes of cell–cell interaction changes (y axis) after the addition of recombinant proteins are compared for cell pairs predicted by the model to be likely to change after perturbing that surface protein (‘true’) and those predicted not to change (‘false’). Each panel considers a different recombinant protein added in the experiment and the corresponding model predictions for that same protein. The two colours for each data point depict the identity of the cell pair according to the colour scale in c. n = 10 samples for the experimental data. Source data.

Extended Data Fig. 1. Developing SAVEXIS.

a. Design of expression vectors and recombinant protein constructs for SAVEXIS screening of divergent architectural classes of cell-surface proteins. For heterodimers, the exact formulation of each chain will depend on the receptor subunit’s topology (e.g. using the Type I vector for integrins, and Type II for CD94/NKG2). b. Empirically gauging streptavidin multimerization stoichiometry by ELISA. Schematic of the procedure for measuring tetramerization around streptavidin by titrating soluble streptavidin (SAV) against a fixed concentration of biotinylated protein before transferring to a streptavidin-coated plate for an ELISA. The measured dilution at which signal ceases represents the optimal tetramerization stoichiometry. c. Calculated stoichiometric equivalence points of 6 example biotinylated proteins incubated with streptavidin. The 4:1 stoichiometric equivalence point inferred mathematically based on molecular mass calculations is indicated on the x axis as “1”. A dashed line indicated the empirically measured median equivalence point. A set of 6 Danio rerio Jam proteins were measured, with the average indicated by a dashed line. d. Titrating prey concentrations identifies a common prey activity with optimal sensitivity. For 4 known zebrafish receptor–ligand interactions of varying affinity (colour shades, ordered by known binding affinity expressed as dissociation half-life), the ratio of raw absorbance signal for the specific interaction against non-specific interactions (y axis) is measured across prey concentrations (y axis). Error bars represent the standard error of the mean. n = 6 technical replicates. e. Applying soluble desthiobiotin greatly enhances the assay signal to noise ratio by sealing unoccupied biotin-binding sites. Measurements of a set of 3 example interactions compared without applying any sealing step in between bait and prey incubations (“Original”), applying biotinylated rCD4 tag at either IC50 or plate-saturation concentrations, or applying a molar equivalent of saturation with either biotin or desthiobiotin. The ratio between mean signal to noise is indicated above each condition. Error bars represent the standard error of the mean. n ≥ 5 independent wells. f. Immobilizing lower quantities of bait protein can reduce off-target signals. Example interaction assays with baits saturating each well (top) or baits at half-maximal dilutions (bottom). A dilution series of prey is applied across each column. The top two rows are specific interactions (binding half-lives listed in green), and the bottom two are known non-interactors. g. Assay miniaturization into 384-well plates can retain and enhance performance. Different assay miniaturization strategies to adapt from 96-well format to 384-well format are indicated along the x axis, including reducing by a half or a third all volumes and protein amounts proportionally (left two), or only reducing volumes while concentrating the samples so the total protein quantities applied are conserved (right two). The ratio between mean signal to noise is indicated above each condition. Error bars represent the standard error of the mean. n ≥ 12 independent wells. h. Representative screening example for human receptors. The appearance of the raw screening plate (left) is shown alongside absorbance values following median polish normalization (right). Bold borders indicate interactions that are expected based on literature publications. One protein, corresponding to a construct later found to be incorrectly labelled, is omitted. i. SAVEXIS consumes small enough quantities of protein that thousands of assays are possible from small input sizes. For a sample set corresponding to all measured proteins from the leukocyte surface receptor library, the amount of interaction tests that could be performed for each protein based on its expression yield was calculated. As a typical case illustration in which protein is purified from a 30 mL cell culture, the median number of interaction tests possible is indicated by a dashed line.

Extended Data Fig. 2. Interaction screen quality controls testify to robust measurements.

a. Recombinant proteins produced at scale match their expected molecular masses. Observed molecular masses from denaturing protein gel electrophoresis (x axis) are compared against computationally predicted molecular mass. Predictions were made by taking the known masses of each amino acid in the protein after processing of its signal peptide, with 2.5 kDa added per predicted N-linked glycosylation site. Shading indicates the 95% compatibility interval for the least-squares linear regression fit. Full images for all gel electrophoresis samples are provided in Supplementary Fig. 1. b. Quantitative protein concentration measurements by Bradford assay agree with qualitative estimates of protein concentration based on densitometry of Coomassie-stained protein gels. Measured concentration percentiles (y axis) are compared against discretized expression categories based on staining intensity (y axis and colour shade). c. Control wells included on every screening plate indicate high consistency across the primary interaction screen. Boxplots of plate measurements for negative control wells (blank baits and tag-only rCD4 baits), positive control wells (the known interaction between P. falciparum P12 and P41 at either a 1x or 3x dilution), and loading control wells (OX68 antibody) that capture prey proteins by recognizing their rCD4 tag. n = 1,262 wells per condition. d. Positive interactions from the primary screen were reproducible in a secondary screen with independently produced and measured proteins. For each protein–protein pair measured, the processed signal in the primary screen (x axis) is correlated against the signal in the secondary screen (y axis). Pairs previously described as being interactors are denoted by red colouration.

Extended Data Fig. 3. Interactome validation summaries.

a. Precision-recall curve corresponding to Fig. 1d. Colour shading indicates the cut-off for the summed screen signal across both bait–prey orientations. The performance of a random classifier is shown by the dotted line, and grey shading indicates the valid range between perfect performance and a random classifier. This curve only considers proteins for which expression was detectable, and defines a positive set based on previously published interactions and a negative set based on randomized interaction pairs. b. ROC curves of screen performance are consistent across possible definitions of positive and negative sets. These curves consider all proteins regardless of whether they were detectably expressed. Columns provide different positive reference sets, and rows delineate possible negative reference sets. High-quality support refers to interactions with experimental support by SPR, isothermal titration calorimetry, analytical ultracentrifugation, or a co-crystal structure. c. Overview of evidence for newly identified interactions. Data from both measured bait–prey orientations in the primary and secondary screen are indicated in pink. Results from cell-binding experiments and SPR are categorized for each interaction along a simple qualitative scale of green to red for ease of comparison. This includes whether a binding response was detectable in SPR equilibrium experiments, if the binding response in SPR experiments was sufficiently quantifiable that 1:1 binding models could be fit, and whether gains in cell-surface binding were observed in cell-based assays when the counter-receptor was overexpressed. The full experimental results that are summarized here can be found in Extended Data Fig. 4.

Extended Data Fig. 4. Orthogonal binding assays to confirm each interaction.

a. Full set of cell-binding traces, extended from Fig. 2a. For each interaction pair, the protein named on the right of the title was transfected into human cells and the protein named on the left was tested for binding as a fluorescently linked recombinant protein tetramer. Tetramer binding to cells (x axis) was measured by flow cytometry at different tetramer concentrations (y axis). The traces in blue are cells overexpressing the indicated receptor, and red shows binding to the mock-transfected control cells. Because cDNAs may express to widely varying levels or not at all, and some proteins may bind to endogenously expressed HEK293 surface proteins, some experiments give inconclusive binding data. For example, HLA-F is known to be predominantly sequestered intracellularly^,, whereas soluble APP is known to already have strong baseline binding activity to cell lines. b. Full set of SPR sensorgrams, extended from Fig. 2b. For each interacting pair, the sensorgram on the left side shows kinetic binding measurements, and the sensorgram on the right shows equilibrium binding measurements. The protein used as the analyte is named on the left and the immobilized ligand is on the right. Association and dissociation constants from a 1:1 binding model fit (red line) are displayed where applicable on the kinetic traces. All analytes were resolved by gel filtration immediately prior to use in binding experiments to reduce the influence of protein aggregates, which otherwise can dominate binding kinetics (Supplementary Fig. 2). Some analytes show clear evidence of two-phase binding kinetics such as PTPRF and MCAM.

Extended Data Fig. 5. Distribution of immune receptor binding affinities.

a. Antigen-presenting cells vary in their receptor contacts with circulating T lymphocytes. Each dot is a cell–cell pair comprising T cells (including all CD4 and CD8 subsets) with either a B cell (red distribution) or dendritic cell (blue distribution). Statistics are overlaid for a Welch’s t-test. b. The overall distribution of immune receptor binding affinities is centred in the range of micromolar dissociation constants. All quantified immune interactions in our network are plotted as a log-scale distribution of binding dissociation constants (K_D). c. Cell surface receptor–ligand binding affinities weakly correlate with protein expression level. For all immune cell types measured by proteomics and all protein interaction pairs measured in our study, the equilibrium binding constant (y axis) is correlated to the summed expression of the two proteins. Weaker interactions (corresponding to higher K_D values) are associated with higher expression. Owing to the number of points, a hex plot is shown with the number of unique combinations of cell type and protein pairs represented by colour shading. A least-squares linear regression line is overlaid in blue. d. Within individual cell types, a surface protein’s expression level weakly and variably correlates with its binding affinity. Instead of all pair combinations, each cell type measured by proteomics is individually presented. One point is placed for each detectable protein with an interaction in that cell type. In cases in which a single protein has multiple interactions of varying affinities, one point is drawn per binding interaction that protein participates in. Least-squares linear regression fit lines and 95% compatibility intervals are overlaid for each cell type. p-values for the Pearson correlation fit are accentuated for conditions where p < 0.10.

Extended Data Fig. 6. Average receptor affinity increases after immune activation.

a. Extended version of Fig. 2e, showing all cell subtypes identified in the original proteomics study instead of higher-level cell-type categories. Each point represents the strength of a protein’s interaction for proteins that are either differentially downregulated upon activation (“Down”) or differentially upregulated upon activation (“Up”). To be classified as differentially expressed, the protein must have more than a 2-fold change upon activation. Data are shown as Tukey boxplots with Holm-corrected p-values calculated from a two-sided Welch’s t-test. b. Combined analysis of all cell types for which paired activated and resting expression data is available. As in panel a, differential expression was defined by more than 2-fold changes upon activation (n = 4 blood donors). Statistics are overlaid for a two-sided Welch’s t-test. c. Same as panel b except setting as the criterion for differential expression that the protein must have a corrected p-value below 0.05 across the 4 proteomics replicates available per condition. Statistics are overlaid for a two-sided Welch’s t-test. d. Same analysis as panel b performed on an independent dataset of RNA-seq measurements on sorted and stimulated immune cell populations. Statistics are overlaid for a two-sided Welch’s t-test (n = 4 blood donors). e. Same analysis as panel a, performed on an independent dataset of RNA-seq measurements (Calderon et al. 2019). Every leukocyte subpopulation measured in activated and resting states is show as a separate box. Data are shown as Tukey boxplots with Holm-corrected p-values calculated from a two-sided Welch’s t-test.

Extended Data Fig. 7. Internal representations of the quantitative cellular interaction model and additional validations.

a. Interaction spectra showing how individual receptor interactions (x axis) contribute to the overall connection strength between cell pairs. The model output per pair, which corresponds to the calculated density of proteins in a bound configuration at a theoretical equilibrium, is shown on the y axis in log₁₀ scale. b. Relative contributions of different protein–protein interactions to a cell pair’s overall connection scores. The colour shading indicates the percent of the total calculated interactions between a given cell pair (x axis) that are attributable to each specific protein–protein interaction (y axis). c. Null distribution of correlations to published data when the model’s cell-surface protein interaction pairs are randomized. Histogram bins representing correlations equal to or greater than the fit of the true model are shaded in light blue after 1,000 random permutations were performed. d. Distribution of model correlation coefficients to previously published leukocyte binding data following a complete leave-one-out analysis of each surface protein. The histogram bin that includes the observed correlation in Fig. 2g is shaded in lighter blue. The two proteins leading to the greatest change in correlation (ICAM1 and PECAM1) are labelled. e. Model fits remain robust on independently measured datasets of leukocyte cellular contacts. Following an analogous approach to Fig. 2g, the kinetic model’s predictions for baseline leukocyte interaction rates were compared to empirical data generated during our pharmacoscopy experiments. Only negative control conditions treated with PBS instead of recombinant protein were included, with interaction scores calculated following the same methodology as the previously published study described in Fig. 2g. Shading indicates the 95% compatibility interval of the least-squares linear regression fit. f. Differential equation model output simulating cell pairs reaching a binding equilibrium. Each colour is either a cell type in an unbound state or a cell pair. The absolute proportions in blood (y axis) were tracked over an arbitrary time scale (x axis) until equilibrium is reached.

Extended Data Fig. 8. Integrating single-cell expression atlases with cellular interactions.

a. Interactively searching immune interactions in human tissues through our web tool. This screenshot depicts the core features of the website, including a drop-down menu to select one of eight different tissue datasets, adjustable sliders for expression cut-off thresholds when determining interactions, and multiple tabs offering different kinds of views for both cellular and molecular interaction types. b. Myeloid-lineage cells do not express any greater quantity of interaction-capable surface proteins than non-myeloid cells. Each dot represents a specific measured cell type in the indicated dataset, which are grouped into categories as myeloid or non-myeloid (x axis). The absolute number of cell-surface proteins for which at least one interaction is annotated with another cell-surface protein (y axis) is compared between categories. p-values from a two-sided Welch’s t-test are shown alongside their corresponding false discovery rate-corrected q values. Tissues match those shown in Fig. 3c. c. Changes in cellular interaction frequencies among immune cells isolated from cancerous versus healthy kidney implicate cellular contacts with potential relevance to pathology. The total number of interactions detected in a single-cell sequencing dataset of paired healthy kidney and kidney tumours are compared across immune populations. d. Comparisons of interactions detected in paired samples of healthy and diseased tissue can suggest functional targets. Human lung tissue from healthy donors and patients with asthma were processed in identical ways through our web tool’s functions. At the indicated expression cut-off which requires the mRNA encoding an interacting surface protein be detected in at least 5% of cells from a given cell type, interactions between cells were categorized based on whether they were present in both healthy and asthmatic samples (purple), only healthy (green), or only asthma (red). Although more qualitative than differential expression based tests, this approach may have utility in conducting more sensitive exploratory analyses of interaction sets. e. Intercellular signalling pathways in tumour-infiltrating immune cells inferred by NicheNet analysis. Two analyses are shown on a single-cell RNA-seq dataset of immune cells isolated from kidney for cell-surface signals that differentially regulate gene expression in plasmacytoid dendritic cells (pDCs) and helper T cells within kidney tumours when compared to adjacent healthy kidney tissue. In the larger box, intercellular signals being received by pDCs are shown (left), matched to genes inferred to be regulated by those signals. In the smaller inset box, the finding that tumour pDCs upregulate JAML is expanded by analysing pDC communication specifically with helper T cells. Targets in our gene regulatory analysis were filtered to exclude those which recurred non-specifically in more than half of all cases.

Extended Data Fig. 9. Microscopy readouts identify that receptor binding partner transcripts are colocalized in human lymph nodes, and that polarization of B and NK cells serves as a marker of classical activation pathways.

a. Single-molecule fluorescent in situ hybridization in human lymph node demonstrates that cells transcribing genes encoding surface proteins found to physically interact in biochemical assays are also physically colocalized in the lymph node. A lymphoid follicle is shown for 4 different transcript pairs encoding proteins that directly interact. Each experiment was repeated on tissue from two donors. The scale bar is 100 μm and applies to all images. b. Polarization of B cells relative to PBS-treated controls after treatment with cytokines and other immunomodulatory molecules. As expected IL-4 and IL-6 are strong B-cell specific activators^,, whereas IL-15 and LPS activate both NK and B cells. Data are shown as Tukey boxplots with Holm-corrected p-values calculated from a one-sided t-test. n = 7 blood donors. c. Polarization of NK cells relative to PBS-treated controls after treatment with cytokines and immunomodulatory molecules. As has been reported for NK cell activation, IL-15 invokes the strongest activation of NK cells, whereas there are no effects from IL-4 and IL-6. The NK cells are also inhibited by the steroid dexamethasone (Dex), consistent with known pharmacology. Data are shown as Tukey boxplots with Holm-corrected p-values calculated from a one-sided t-test. n = 7 blood donors.

Extended Data Fig. 10. High-content microscopy datasets show human leukocyte phenotypes after infusion of purified proteins.

a. Full dataset of changes to cellular state proportions, extended from Fig. 4b. The proportions of measured cell types and cell states (x axis) are compared across different protein doses, timepoints, and with or without background LPS activation (y axis). Points are sized by their adjusted p-values, and shaded to show relative change compared to controls. n = 5 wells. b. Observed immunomodulatory phenotypes match the expression profiles of the cell types expressing the applied protein or its identified binding receptor. For every statistically significant phenotype upon protein application that is described in Fig. 4b (grey top bars), the protein itself and all its identified interaction partners (columns) had their expression in lymphocyte populations (x axis) compared according to quantitative proteomics measurements (y axis). The lymphocyte populations are coloured blue if they match the cell type for which the significant phenotype was found. c. Clustering on the full set of cell-to-cell connectivity changes identifies recurring modules of cellular shifts. All protein conditions and timepoints (y axis) were hierarchically clustered by the complete linkage method. As above, changes in each cell pair’s interactions (y axis) are indicated by the provided colour scale and sized based on their adjusted p-values. n = 10 wells.