In Type 1 diabetes (T1D), insulin-producing β-cells are killed by islet-infiltrating immune cells in a process called “insulitis”. T1D results in a critical requirement for exogenous insulin and affects over eight million individuals world-wide with an estimated 0.5 million new diagnoses each year (Gregory et al. 2022).

Recently, the first immunotherapy for delaying T1D onset, teplizumab (a human anti-CD3 monoclonal antibody) was approved by the US Food and Drug Administration (Hirsch 2023). However, this treatment and other immunotherapies help only a small fraction of patients and are significantly less effective after patients progress to overt T1D (Herold et al. 2013; Perdigoto et al. 2019; Herold et al. 2019; Pescovitz et al. 2009; Orban et al. 2011; 2014; Bluestone, Buckner, and Herold 2021). A better understanding of T1D pathogenesis is essential to building on this progress.

One of the challenges of studying human T1D pathology is the availability of suitable tissue samples. Obtaining pancreatic biopsies raises the risk of surgical complications and the progressive nature of T1D would necessitate serial, longitudinal studies over time, which is prohibitive (Krogvold et al. 2014). Fortunately, the Juvenile Diabetes Research Foundation (JDRF) Network for Pancreatic Organ Donors with Diabetes (nPOD) have provided human pancreatic tissues from cadaveric donors for this study (Campbell-Thompson et al. 2012; Pugliese et al. 2014). nPOD has enabled substantial progress towards characterizing the pathology of human T1D (Wilcox et al. 2016; Arif et al. 2014; Leete et al. 2016; Martino et al. 2015; Korpos et al. 2021).

Our understanding of key features of human T1D pathology remains limited. Although the cellular composition of insulitis, inflammation specifically of the islets, has been studied extensively, this has been done in separate studies looking at different tissue sections, prohibiting an understanding of how the numerous cellular programs in insulitis are coordinated throughout disease. This was recently addressed using Imaging Mass Cytometry (IMC), which uncovered alterations in β-cell phenotypes, immune composition, vascular density, and basement membrane (Damond et al. 2019; Wang et al. 2019). However, these studies did not deeply phenotype islet-infiltrating CD8⁺T cells, believed to be major driver of β-cell elimination.

Recently, intriguing differences in the extra-islet spaces of T1D and healthy controls have been reported. First, the abundance of multiple types of immune cells outside islets are increased in T1D patients compared to non-T1D controls (Rodriguez-Calvo et al. 2014; Campbell-Thompson, Rodriguez-Calvo, and Battaglia 2015; Bender et al. 2020). Second, HLA-DR expression is increased on ductal cells in T1D tissue donors, hinting at a functional link with CD4⁺T cells (Fasolino et al. 2022). Third, peri-insulitis, the accumulation of immune cells outside islets, is observed in tissues from patients with T1D (Korpos et al. 2013), indicating that not all T cells enter the pancreas directly via islet microvasculature (Savinov et al. 2003). Fourth, in human T1D, but less so in most animal models, islets in different regions of the pancreas are infiltrated at strikingly different rates for reasons that are unknown (In’t Veld 2014). This suggests that the extra-islet compartment could be responsible by governing the targeting of islets. Finally, tertiary lymphoid structures (TLS) – dense aggregates of lymphoid cells indicative of local immune activation – are observed outside islets in T1D patients (Korpos et al. 2021).

In summary, analyzing both compartments simultaneously could help identify how these extra-islet factors influence islet pathogenesis. However, to date, multiplexed imaging studies have only examined islets. A comprehensive, spatially resolved cellular analysis of both compartments in T1D is lacking.

Here, we investigated the islet and extra-islet pancreas together. We used CO-Detection by indexing (CODEX) with an antibody panel targeting 54 antigens to samples from a cohort of T1D patients with insulitis as well as non-T1D individuals with and without islet-specific autoantibodies (AA- and AA+ respectively) obtained through the JDRF nPOD program. We analyzed approximately 2000 islets and broad swaths of the extra-islet tissue to evaluate local and distal spatial architecture. We then used pseudotime analysis to characterize insulitis sub-states based on the activation states of islet-infiltrating CD8⁺T cells. We further investigated the cellular changes in niches and lobules beyond islets. Our results implicate both the islet microenvironment and inflammation at distal sites within the pancreas in the progression of insulitis.

Cohort curation, image acquisition, and cell annotation

The JDRF nPOD is a national registry of cadaveric pancreases donated by T1D patients that has transformed the ability of researchers to investigate the pathways underlying the progression of human T1D (Campbell-Thompson et al. 2012; Pugliese et al. 2014). Insulitis is present in only a small fraction of T1D cases, including nPOD’s (In’t Veld 2011; Atkinson and Mirmira 2023). Although nPOD has close to 200 T1D cases, at the time of our study, nPOD had only 17 with documented insulitis. For some of these 17, tissues had been heavily sectioned and the current tissue blocks did not still contain insulitis. AA+ and T1D cases were selected by performing triple-immunohistochemistry for Insulin, Glucagon, and CD3. T1D and AA+ cases that had CD3+ staining in islet or peri-islet spaces and tissue still available were selected for our study. The final cohort included two AA+ cases, eight T1D cases, and three non-T1D controls. Given that insulitis is not detected in non-T1D cases (Bruggeman et al. 2021), the blocks from controls were selected randomly. The T1D cases varied in the time between diagnosis and death from 0 years (diagnosed at death) to 6 years (Figure 1.A, left). The causes of death were mostly unrelated to T1D complications (Supplemental Table 1). Therefore, the time since diagnosis is not a reflection of the severity or aggressiveness of the individual’s disease.

Figure 1

Large regions averaging 55 mm² were imaged with CODEX as previously described (Schürch et al. 2020; Phillips, Matusiak, et al. 2021; Hickey et al. 2021). Regions were selected to capture islets and the surrounding region simultaneously (Figure 1.A, center and right). CellSeg was used to segment cell nuclei and quantify marker expression from CODEX images as previously described (Lee et al. 2022). In total, our dataset consisted of 7.0×10⁶ cells across all donors (ranging from 3.0×10⁵ to 9.8×10⁵ cells per donor). Twenty-one cell types were identified with Leiden clustering and manual merging and visualized using Uniform Manifold Approximation and Projection (UMAP) (Figure 1.B, Supplemental Table 2). Endocrine cells were manually gated from UMAP embeddings derived from Proinsulin, Glucagon, and Somatostatin to identify b-cells, a-cells, and d-cells respectively. Immune cells were sub-clustered with the Leiden algorithm using immune-specific markers (Figure 1.C, Supplemental Table 2). To verify the accuracy of our annotations, we overlaid labelled data onto the original images (Supplemental Figure 1.A). Of note, we could not accurately identify macrophage subsets or distinguish dendritic cells from macrophages due to the panel design, complex combinations of co-expression, and the difficulty in segmenting and quantifying markers on myeloid populations due to their morphology. Therefore, we refer to this cluster as “macrophage/DCs” In addition, we identified a cell population that could not be definitively annotated, expressing high levels of CD45, CD69, Granzyme-B, and CD44, intermediate levels of CD16, S100A6, Galectin-3, and Hyaluronan, but not expressing CD3, CD20, CD56, CD57, CD15, or MPO. We confirmed from the raw images that CD3, CD4, and CD8 were not internalized, indicating activation, nor did these cells express other T cell activation markers CD45RA, CD45RO, PD-1, or LAG-3 (Figure 1.C, Supplemental Figure 1A, bottom right). This population could represent a type of innate lymphoid cell (Dalmas et al. 2017) and was labeled “Granzyme-B⁺/CD3-”.

Islet- and extra-islet regions are altered in T1D

We first sought to identify cellular changes in T1D within islets specifically. Previous reports observed that insulin-containing islets are significantly more common in recent-onset T1D cases than cases with diabetes durations of greater than one year (In’t Veld 2011; Campbell-Thompson et al. 2016; Richardson and Pugliese 2022). Similarly, we found that samples from patients who had been diagnosed with T1D for 0-2 years had significantly reduced β-cell frequencies compared to non-T1D controls. Furthermore, samples from subjects with disease durations of 5-6 years had minimal remaining β-cell mass (Figure 1.D). Whereas one AA+ case had β-cell mass comparable to those of cases with disease duration of 0-2 years, the other AA+ case was comparable to non-T1D controls (Figure 1.D). The total islet area imaged was comparable across all donors (Supplemental Figure 1.B).

Next, we investigated how the abundances of non-endocrine cell types inside islets differed across donors. We performed Principal Component Analysis (PCA) on the donors using the frequencies of non-endocrine cell types located in islets. Donors were clearly separated into two groups by the first two principal components; one group included all T1D cases and one AA+ case and the second included all non-T1D cases and the other AA+ case (Figure 1.E). In this analysis, we did not consider β-cells, ⍺-cells, and δ-cells. Thus, donors were stratified by disease duration strictly according to the abundances of immune and other pancreatic, non-endocrine cell types in the islets.

We next considered only cells located outside islets. Again, donors were clearly separable by the first two principal components (Figure 1.F). The first principal component separated cases with times since diagnosis between 0-2 years from non-T1D, AA+, and cases with diabetes durations of 5-6 years (Figure 1.F). The second principal component separated cases with diabetes durations of 5-6 years from the rest (Figure 1.F). Therefore, both the islet and extra-islet spaces of T1D and non-T1D cases were distinct.

Many cell types were increased in T1D cases with times since diagnosis of 0-2 years relative to non-T1D controls (Supplemental Figure 1C). In T1D cases with times since diagnosis of 5-6 years, the abundance of different cell types either remained higher than non-T1D controls or returned to baseline (Supplemental Figure 1.C). This trend was present in both islet and extra-islet regions. These data demonstrate that the immune activity between the islet and extra-islet compartments are coordinated.

Pseudotemporal reconstruction of islet pathogenesis identifies a conserved trajectory of insulitis

In human T1D, β-cell destruction does not occur simultaneously across all islets and even neighboring islets can be at different stages of destruction (In’t Veld 2011; 2014; Damond et al. 2019). We therefore used pseudotime analysis to infer the most likely progression of a single islet through disease space (Damond et al. 2019). To develop a pseudotemporal map, we quantified the cellular composition of each islet, including cells within 20 µm of the islet’s boundary, and applied the pseudotime algorithm PArtition-based Graph Abstraction (PAGA) (Figure 2.A; Figure 2.B), (Wolf et al. 2019). PAGA was selected because it is a high-performing algorithm able to identify multiple trajectories, if they exist, while making minimal assumptions about the true structure (Saelens et al. 2019).

Figure 2

Displayed in Figure 2.B is the PAGA force-directed layout where each point represents an islet. Each islet’s color reflects the pseudotemporal distance from the centroid of non-T1D islets. As expected, the islets from different donor groups (no T1D, AA+, T1D) had different distributions across the PAGA map (Figure 2.C). In the PAGA map, a continuum is apparent from islets abundant in insulin-expressing β-cells on the left of the map to islets depleted in β-cells on the right (Figure 2.D, Figure 2.E, Figure 2.F, top row). PAGA uses Leiden clustering internally, enabling the following regions of the pseudotime map to be labelled objectively: 1) Islets with low pseudotime values on the left of the map (PAGA-internal Leiden clusters 0 and 5 in Supplemental Figure 2.A) were labelled “Normal” even if they originated from T1D donors. 2) Islets in the middle of the map (PAGA-internal Leiden clusters 6, 2, and 8 in Supplemental Figure 2.A) were elevated in HLA-ABC (MHC Class I) expression, CD8⁺T cells, and macrophage/DCs (Figure 2.D, Figure 2.E, Figure 2.F, rows 2-4) and were labelled “Inflamed”. 3) Islets with late pseudotime values on the right of the map (PAGA-internal Leiden clusters 1,3,7, and 4 in Supplemental Figure 2.A) were devoid of β-cells and were labelled “Insulin-Depleted” (Figure 2.D, Figure 2.E, Figure 2.F, top row).

In addition, islets lacking β-cells occasionally contained CD8⁺T cells and were labelled “Insulin-Depleted + Immune Islets” (Figure 2.D, Figure 2.E, Figure 2.F, rows 2-4). The presence of these islets suggests that signals retaining CD8⁺T cells in islets linger after β-cells die. The distribution of all cell types across pseudotime is reported in Supplemental Figure 2.B.

Islets from non-T1D controls and one of the AA+ donors (6314) were primarily in the Normal group to the left of the map (Figure 2.C, Supplemental Figure 2.C). Islets from subjects who had had T1D for of 5-6 years (cases 6195 and 6323) were primarily in the Insulin-Depleted group to the right of the map (Figure 2.C, Supplemental Figure 2.C). All the remaining T1D donors and the other AA+ donor were distributed broadly throughout the map (Figure 2.C, Supplemental Figure 2.C).

We quantified the fraction of each cell type in swaths at varying distances from the islet edge. We found that for B Cells, CD4⁺T cells, CD8⁺T cells, macrophage/DCs, neutrophils, and plasma cells, the fraction of the given cell type in the islet relative to outside the islets increased (Supplemental Figure 2.D), demonstrating that the inflammation in islets was distinct from the inflammation in the extra-islet tissue. Thus, immune cells were targeting islets specifically.

Together, these results illustrate a single, non-branching progression from Normal Islets to Insulin-Depleted Islets via Inflamed Islets, consistent with previous pseudotime analyses (Damond et al. 2019).

IDO expression on islet vasculature is linked to T cell infiltration

We observed islets with vasculature staining positive for indoleamine 2, 3-dioxygenase 1 (IDO). In the tumor microenvironment, IDO is commonly expressed by myeloid cells and suppresses CD8⁺T cell activity through multiple mechanisms including the induction of FOXP3⁺ regulatory T cells and the inhibition of CD8⁺T cell function (Munn and Mellor 2016). In islets, IDO co-stained with CD31⁺ vasculature but not CD45⁺ immune cells adjacent to vasculature (Figure 2.G). IDO was not expressed by endocrine or other cell types in islets or by vasculature or any cell type outside islets (Supplemental Figure 2.E). We manually quantified vascular expression of IDO on islets throughout pseudotime and found that all but two IDO⁺ islets were in the Inflamed group (Figure 2.H). Therefore, IDO expression by islet vasculature was tightly associated with insulitis.

A potent inducer of IDO expression is Interferon-γ (IFN-γ), a cytokine highly expressed by activated T cells and macrophages (Munn and Mellor 2016). Therefore, we hypothesized that IDO expression was induced by infiltrating immune cells during insulitis. We compared the frequency of CD8⁺T cells and macrophage/DCs in islets from the Inflamed group with and without IDO⁺ vasculature and found that CD8⁺T cells were significantly more abundant in islets with IDO⁺ vasculature than islets without IDO+ vasculature (Figure 2.I). However, the abundance of γδ-T cells and CD4⁺T cells was less strongly associated with IDO+ vasculature and the abundance of macrophage/DCs was not significantly associated with IDO⁺ vasculature (Supplemental Figure 2.F).

In summary, IDO expression by islet vasculature is positively associated with T cell infiltration and may be an immunoregulatory checkpoint.

Insulitis has sub-states, defined by functional states of CD8⁺T cells

CD8⁺T cells are a major component of insulitis (row 4 of Figure 2.E and Figure 2.F) and are capable of directly and indirectly killing β-cells. A comprehensive description of the activation profiles of CD8⁺T cells could provide insight into their roles in T1D pathogenesis. To obtain extremely high-quality marker quantification, we trained a neural network on manually labelled images of single T cells (Figure 3.A, Supplemental Figure 3.A). Using our neural network, we quantified the expression of T cell markers on islet CD8⁺T cells (Figure 3.B).

Figure 3

PD-1, TOX, CD45RO, CD69, and CD44—markers of antigen experience—were the most commonly expressed by islet-infiltrating CD8⁺T cells (Figure 3.B, Supplemental Figure 3.B). CD8⁺T cells expressing CD45RA (which are either naïve or terminally differentiated effector memory cells (TEMRA)) were detectable in islets, as previously reported (Damond et al. 2019) (Figure 3.B). In addition, we observed a rare population of CD45RO⁺ CD8⁺T cells co-expressing Lag-3, Granzyme-B, and ICOS (Figure 3.B bottom clade). Lastly, a rare population of CD57⁺ CD8⁺T cells was present but these cells rarely co-expressed LAG-3, Granzyme-B, or ICOS (Figure 3.B top clade). These populations resemble the two exhausted T cell populations identified in the peripheral blood of T1D patients that were associated with responsiveness to alefacept (Diggins et al. 2021). Therefore, the activation profiles of islet-infiltrating T cells are heterogeneous.

We reasoned that the islet microenvironment may dictate the activation state of CD8⁺T cells by specifically recruiting T cells of a particular state and/or inducing changes after they enter the islet. If so, islets would contain specific combinations of CD8⁺T cell states. To interrogate this, we performed UMAP only on Inflamed islets, using the frequencies of CD8⁺T cells expressing each functional marker. We identified four inflamed sub-clusters, I-IV, (Figure 3.C top). Here, the term “sub-cluster” is used to highlight that these groups were all contained within the previously defined “Inflamed” cluster and the roman numerals do not imply a temporal ordering. Inflamed-I contained only CD8⁺T cells that did not express any of the functional markers analyzed (Figure 3.C bottom, top row). Inflamed-II was characterized by a high frequency of CD45RA⁺CD8⁺T cells (Figure 3.C bottom, second row from top and Figure 3.D top row). Inflamed-III was characterized by a low frequency of CD45RA⁺ cells and high frequency of CD45RO⁺ and PD-1⁺ cells (Figure 3.C bottom, third row from top and Figure 3.D middle row). Inflamed-IV was characterized by an enrichment of CD8⁺T cells expressing CD57, LAG-3, ICOS, Granzyme-B, PD-1 or CD45RO (Figure 3.C bottom, bottom row and Figure 3.D bottom row). In summary, the phenotypes of CD8⁺T cells are coordinated across islets.

Regulation of insulitis sub-states by the islet microenvironment

To identify cellular or molecular factors that regulate the state of CD8⁺T cells in islets, we first inspected the distribution of inflamed sub-clusters in each patient. Each donor possessed islets that belonged to multiple inflamed islet sub-clusters (Figure 3.E). Therefore donor-level factors such as genetics, the location within the pancreas (i.e. head, body, or tail), or time since T1D onset, are not associated. Instead, these insulitis sub-states are conserved among T1D patients.

Next, we asked if each T cell marker is enriched in CD8⁺T cells in islets compared to CD8⁺T cells in the peri-islet and exocrine space. We computed the frequencies of each CD8⁺T cell state inside islets of each inflamed sub-cluster and in separate swaths 0-25µm, 25-50µm and 50-100µm away from the islets (Supplemental Figure 3.C). We found that for islets of Inflamed-II, -III, and -IV, functional markers characterizing their CD8⁺T cells were expressed more frequently inside than in the surrounding tissue areas. Although the different functional markers were all enriched on islet-infiltrating T cells compared to T cells outside islets, the degree of this enrichment varied across the markers measured. The markers of Inflamed-IV, LAG-3, ICOS, and Granzyme-B, were highly enriched inside islets (∼20% of islet-infiltrating CD8⁺T cells vs <5% of extra-islet CD8⁺T cells). CD45RA and CD69 in Inflamed-II were slightly less enriched in islets (∼28% of islet-infiltrating and ∼20% of extra-islet CD8⁺T cells). CD45RO and PD-1 in Inflamed-III islets were the most abundant but had the least enrichment in islets (∼45% of islet-infiltrating CD8⁺T cells vs ∼35% of extra-islet CD8⁺T cells) (Supplemental Figure 3.C). This demonstrated that the specific differences in the compositions of CD8⁺T cell states in different islets were attributable to the islet microenvironment and not the surrounding extra-islet spaces.

Although macrophage/DCs are abundant in islets from the Inflamed group (Figure 2.E, Figure 2.F) and are capable of interacting with T cells through antigen presentation and cytokine secretion, neither the expression of markers of macrophage/DC activity nor macrophage/DC abundance was significantly associated with any of the inflamed sub-clusters (Supplemental Figure 3.D). Similarly, no other cell-type nor the vascular expression of IDO was linked to CD8⁺T cell programs in islets (Supplemental Figure 3.D). Accordingly, the four inflamed sub-clusters had identical distributions throughout the original PAGA force-directed layout (Figure 3.F). Therefore, the activation states of islet-infiltrating T cells are independent of the abundance of any of the other cell types we could identify.

Lastly, we compared CD8⁺T cells in islets from the Insulin-Depleted + Immune group to CD8⁺T cells in islets from the Inflamed group. Insulin-Depleted + Immune islets contained a higher frequency of CD45RA+ CD8⁺T cells and a lower frequency of CD45RO+ CD8⁺T cells than Inflamed islets (Figure 3.G). TOX was expressed by a higher frequency of CD8⁺T cells in Insulin-Depleted + Immune islets than CD8⁺T cells in Inflamed islets (Figure 3.G). Importantly, CD45RA+ and TOX were never co-expressed on the same CD8⁺T cell (Figure 3.B). These data indicate that the activation state or persistence of these two populations in islets depends on insulin expression.

Vasculature, nerves, and Granzyme-B⁺/CD3-cells outside islets are associated with lobular patterning

The destruction of islets in T1D is known to exhibit lobular patterning: that is, islets in the same lobule are likely to be in the same stage of disease (Gepts 1965). The reason why islets within the same lobule are synchronized is unknown. One explanation could be an islet-intrinsic mechanism where the expression of programs sensitizing beta cells to immune-killing (i.e. stress) are correlated across islets in the same lobule. Alternatively, it could be mediated by cells outside islets if they facilitate extravasation into the lobule or trafficking from one afflicted islet to the next.

To systematically investigate lobular patterning in T1D, we used a neural network to segment lobules and assign each single cell and islet, to its lobule. We first quantified the degree of lobular patterning within each donor using the intra-class correlation coefficient (ICC). The ICC ranges from 0 to 1 where cases with values closer to 1 have islets whose states are more synchronized within lobules (Figure 4.A). Islets of non-T1D cases and 6314, 6195, and 6323 did not have appreciable variability in their pseudotimes, but in the remaining cases, the ICCs ranged from 0.17 to 0.74 (Figure 4.B). This highlights that the magnitude of lobular patterning ranges widely across T1D cases with insulitis.

Figure 4

We employed hierarchical linear modeling (HLM), a statistical framework designed to identify relationships between levels of multi-level data. HLMs are standard in fields where multi-level data are common such as Education, in which students are grouped into classrooms, which are grouped into schools (Gelman et al. 2014) and have been applied in biomedical settings (Jerby-Arnon and Regev 2022; Yi et al. 2019). We were interested in cell-types if their abundance in a lobule correlated with the lobule’s average islet pseudotime. Importantly, we omitted the islet region itself from the calculation of a cell type’s lobular abundance because we were interested in identifying features in the extra-islet space that were associated with islet destruction. For each cell type, we estimated the effect of its total abundance in a lobule (the number of cells divided by the number of acinar cells to normalize for lobule area) on the pseudotimes of islets in that lobule. We performed this analysis in two-level HLMs for each donor and a three-level HLM considering all donors together (Figure 4.C).

We identified three cell types, vasculature, Granzyme-B⁺/CD3-cells, that were significantly associated with lobules across multiple T1D tissue donors. All three were more abundant in lobules with islets late in pseudotime (Figure 4.C boxed rows, Figure 4.D). Samples from cases 6323 and 6195 which had very few insulin-containing islets had increased abundances of vasculature, Granzyme-B/CD3-cells, and nerves in their extra-islet spaces compared to non-T1D controls (Supplemental Figure 1C), indicating these changes persist through the point when the entire tissue is afflicted. In addition, vasculature, Granzyme-B⁺/CD3-cells, and nerves were increased in Inflamed islets compared to Normal islets indicating that they may serve a role in islets in addition to their role in the extra-islet compartment (Supplemental Figure 4A). Note that Supplemental Figure 4A differs from Supplemental Figure 1C because islets are broken up according to pseudotime, not donor, and not all islets in T1D donors are undergoing insulitis.

Counterintuitively, although CD8⁺T cells and macrophage/DCs were higher in the extra-islet compartments of T1D cases vs non-T1D cases (Supplemental Figure 1C), they were not associated with lobular patterning (they were not more abundant in the extra-islet space of lobules with more advanced insulitis) (Figure 4.C). These data raise the possibility that vasculature, Granzyme-B/CD3-cells, and nerves outside islets help predispose lobules to insulitis or are affected by extensive insulitis.

Immature tertiary lymphoid structures are enriched in subjects with T1D

We hypothesized that T cells’ interactions with certain cell types in specific areas of the pancreas may be important for their functionality. Therefore, we identified Cellular Neighborhoods (CNs) (Schürch et al. 2020; Bhate et al. 2021), tissue regions that are homogeneous and have defined cell-type compositions. To identify CNs, single cells were clustered according to the cell-type composition of their twenty nearest spatial neighbors and automatically annotated with the names of enriched cell types using (Figure 5.A, See Methods). This resulted in 75 CNs. Throughout the manuscript, CNs are referred to with the nomenclature (Cell Type A|Cell Type B|…) to indicate all the cell types that are enriched in them (See Methods).

Figure 5

Next, we identified CNs that were more abundant in T1D than non-T1D tissues (Figure 5.B). The top three CNs (fold change of abundance in T1D relative to abundance in non-T1D) were (CD8⁺T cells|B Cells), (Macrophage|Stromal Cells|B Cells), and (Vasculature|B Cells) (Figure 5.B, Figure 5.C). We asked whether these three CNs were commonly adjacent to each other as this could indicate that they act as components of a larger structure (Bhate et al. 2021). Measuring the frequency with which the three B cell CNs were adjacent to each other throughout the tissues demonstrated that the (CD8⁺T cells|B Cells) CN is predominantly found adjacent to both the other CNs but that (Macrophage|Stroma|B Cells) and (Vasculature|B Cells) are less commonly adjacent to each other (Figure 5.D).

We next asked whether these CN assemblies corresponded to either peri-vascular cuffs (Agrawal et al. 2013; Wekerle 2017) or tertiary lymphoid structures (TLSs) (Korpos et al. 2021; Rovituso et al. 2016; Agrawal et al. 2013), as these are two B Cell-rich structures commonly present in autoimmune conditions. Although the (CD8⁺T cells|B Cells) CN was adjacent to vessels (Figure 5.D, Figure 5.E), it was not in the fluid-filled perivascular space, as is the case with perivascular cuffs (Figure 5.E). In our samples, the (CD8⁺T cells|B Cells) CN did not have segregated T cell and B cell zones as seen in mature TLSs, consistent with (Korpos et al. 2021).

In summary, the (CD8⁺T cells|B Cells) CN is more abundant in T1D tissues from patients with diabetes durations of 0-2 years compared to non-T1D tissues and T1D tissue from patients who had T1D for more than 4 years.

Immature tertiary lymphoid structures are in the extra-islet pancreas, and are enriched in CD45RA⁺ /CD8⁺ T cells

We next asked whether the (CD8⁺T cells|B Cells) CN had high endothelial venules (HEV), specialized blood vessels that are commonly found in TLSs that enable naïve lymphocytes to extravasate into peripheral tissues. We observed expression of peripheral lymph node addressin (PNAd), an HEV marker, in the vessels associated with the (CD8⁺T cells|B Cells) CN (Figure 5.E left image) but not in other vessels (data not shown). Although we could not assess the presence of other TLS traits such as follicular dendritic cells, fibroblastic reticular cells, or follicular helper T cells, the aggregation of B cells and presence of HEVs, but the lack of compartmentalized B and T cell zones indicate that instances of the (CD8⁺T cells|B Cells) CN are immature TLSs.

Next, we asked if immature TLSs could support the entry of naïve CD8⁺T cells into the pancreas. We observed CD8⁺T cells co-expressing CD45RA and CD62L (the ligand for PNAd) near PNAd⁺ vasculature (Figure 5.E, middle and right image respectively). Thus, CD45RA⁺/CD8⁺T cells in the pancreas can adhere to HEVs. Furthermore, CD45RA⁺ was enriched three-fold on CD8⁺T cells in the (CD8⁺T cells|B Cells) CN relative to CD8⁺T cells in the tissue as a whole (Figure 5.F), providing additional evidence that CD45RA+ T cells may enter the pancreas through HEVs.

We found immature TLS both adjacent (Figure 5.G.1) or not adjacent (Figure 5.G.2) to islets. Quantifying the frequency of this adjacency revealed that fewer than half of the immature TLS were adjacent to islets (Figure 5.H). We did not identify any differences in the cellular composition of the immature TLS that were or were not adjacent to islets (Supplemental Figure 5). We reasoned that even if immature TLSs were far from islets, extravasating cells may migrate to islets. Accordingly, islet-adjacent CD45RA⁺ CD8⁺T cells (that were not in islet-adjacent TLSs) co-expressed CD62L, suggesting that they originated from the (CD8⁺T cells|B Cell) CN (Figure 5.I). Consistent with this, in one notable tissue donor, regions of the pancreas with Insulin-Depleted islets were enriched in the (CD8⁺T cells|B Cell) CN relative to regions of the pancreas with β-cell containing islets (Figure 5.J). This spatial correlation between the (CD8⁺T cells|B Cells) CN and the destruction of islets implicates immature TLS with islet pathology even if they are not adjacent to islets (Figure 5.J).

Discussion

We have performed CODEX imaging and comprehensive computational analysis of whole cadaveric pancreata from T1D subjects. Our data support several conclusions.

First, our results are consistent with the model previously proposed by Damond et al, who proposed a single trajectory for insulitis, characterized by an enrichment in HLA-ABC expression, CD8⁺T cells, and macrophage/DCs (Damond et al. 2019).

Second, we are the first to report that IDO+ vasculature is present in inflamed islets but not in normal islets or islets that have lost insulin-expression (Figure 2.G,Figure 2.H). Furthermore, islets with IDO+ vasculature contained higher frequencies of CD8⁺T cells, but not macrophage/DCs compared to inflamed islets that did not contain IDO+ vasculature, suggesting that IDO is induced by a cytokine produced by infiltrating CD8⁺T cells such as IFN-γ (Figure 2.I). Given IDO’s established tolerogenic role, these data suggest that the loss of IDO on vasculature could be a prerequisite for β-cell death. Leveraging this checkpoint to protect transplanted β-cells from rejection has shown promise (Alexander et al. 2002) and could be combined with similar approaches using programmed death-ligand 1 (Yoshihara et al. 2020; Castro-Gutierrez et al. 2021)

We did not observe IDO expression on β-cells, in contrast to Anquetil et al. that report this in healthy patients (Anquetil et al. 2018). Our observations are consistent with data demonstrating that IDO needs to be induced. First, it has been shown with western blot and RT-PCR that human islets require cytokine stimulation to express IDO (Sarkar et al. 2007). Second, the human protein atlas has tested multiple IDO antibodies and demonstrated that IDO is negative in human islets via immunohistochemistry: (https://www.proteinatlas.org). In addition, the human protein atlas’s single-cell RNAseq atlas only report IDO transcripts in immune cells in healthy pancreas.

Third, we performed the first high-dimensional spatial phenotyping of CD8⁺T cells in T1D islets. We found that most T cells were antigen experienced. A small population expressed CD45RA and CD69, which could be naïve or TEMRA cells (Figure 3.B). Another population expressed Lag-3, Granzyme-B, and ICOS. It is notable that only a small population of islets had Granzyme-B-expressing T cells. This could indicate that alternative mechanisms are contributing to the elimination of β-cells.

Fourth, the insulitis trajectory is comprised of four sub-clusters, each characterized by the activation profile of the islet-infiltrating CD8+T cells (Figure 3.C). Multiple of these inflamed sub-clusters were present in all T1D donors, indicating that the sub-clusters are capable of inter-converting (Figure 3.E). The factors that regulate the conversion of a given islet between sub-clusters could correspond to immunoregulatory checkpoints that are critical to the progression of T1D. Unfortunately, our search for such features failed to generate any candidates (Supplemental Figure 3.D). By phenotyping T cells at different distances from the islet edge, we were able to determine that the T cell activation profiles characterizing each sub-cluster were only present in the islet, not in the surrounding tissue (Supplemental Figure 3.C). Unfortunately, from our data, we cannot speculate whether the insulitis sub-clusters arose due to differential stimulation of T cells that had already entered islets, differential recruitment of pre-activated T cells, or both.

Fifth, pancreatic lobules affected by insulitis are characterized by distinct tissue markers. We discovered that lobules enriched in β-cell-depleted islets were also enriched in nerves, vasculature, and Granzyme-B⁺/CD3-cells, suggesting these factors may make lobules permissive to disease (Figure 4.C). The role of islet enervation in T1D has been studied but such work has focused on nerves in the islet rather than on nerves in the exocrine tissue (Christoffersson, Ratliff, and von Herrath 2020). The Granzyme-B+/CD3-cells could be natural killer cells; if so, they are most likely of the CD56dim subset as CD56 was not detected on these cells. It is noteworthy that the cell types linked with direct islet invasion were distinct from those linked to lobule patterning even though both sets of cell types were found across islet and extra-islet regions. Therefore, for insulitis to consume every islet, crosstalk may be required between the cell types in the islet and extra-islet compartments. Conversely, inhibiting this interaction might contain pathology to isolated lesions.

Finally, we identify immature TLS away from islets where CD45RA⁺ CD8⁺T cells aggregate. We also observed an inflamed islet-subcluster characterized by an abundance of CD45RA⁺/CD8⁺ T cells. It will be important to determine whether the CD45RA⁺ T cells localized around islets may have originated from immature TLS. In mice, blocking immune egress from lymph nodes reduced the size of TLSs and halted diabetes (Penaranda et al. 2010). Thus, therapeutic targeting of immune cell trafficking to TLS could help mitigate autoimmunity in human T1D.

Together, these data illuminate relationships between insulitis, the local islet microenvironment and inflammation at distal sites.

A major limitation for the study is the cohort size. Cases with documented insulitis are very rare, significantly limiting the feasibility of curating large cohorts (Campbell-Thompson et al. 2016). Due to this limitation, factors such as the donors’ histories of drug use, durations of stay in the intensive care unit, and BMIs could not be balanced or statistically adjusted for but should be considered when relevant and are associated with exocrine inflammation but not the prevalence of insulitis (Bruggeman et al. 2021). In addition, only one of the AA+ cases, Case 6267, had detectable insulitis (Figure 2.C). The other case, 6314 tested positive for only one autoantibody, GADA+, and therefore had a significantly lower probability of progressing to overt T1D (Ziegler et al. 2013).

Another limitation is our limited perspective on myeloid cell populations. Although antibodies in our panel detect numerous myeloid markers, we failed to identify any heterogeneity in myeloid populations during insulitis. This was likely due in part to the difficulty of segmenting myeloid cells and quantifying marker expression due to their morphology. Spatial transcriptomics could be used in future studies to better define the myeloid populations and inform additions to future CODEX panels.

Lastly, our samples are 2-dimensional sections which could affect some of the adjacency relations.

In conclusion, using a data-driven approach, we mapped conserved sub-states of insulitis and integrated the spatial pathology of islet and extra-islet regions into a single model of T1D pathogenesis. The tools and computational pipelines developed here will enable further investigation of immune pathology at the tissue scale that may lead to development of therapies for T1D.

Methods

Human tissues

Cadaveric pancreatic FFPE tissue sections were obtained through the nPOD program, sponsored by the Juvenile Diabetes Research Fund. Case numbers cited herein are assigned by nPOD and comparable across nPOD-supported projects. 17 cases in the nPOD biorepository had been previously documented to contain insulitis. For each of these 17 cases, we examined the triple stained immunohistochemistry images (CD3, Insulin, and Glucagon) using nPOD’s online pathology database to select blocks in which insulitis was present. To ensure that the tissue regions still contained insulitis (and had not been sectioned extensively after their images were uploaded to the nPOD pathology database), we re-sectioned and visualized CD3, Insulin, and Glucagon. The use of cadaveric human tissue samples is approved by Stanford University’s Institutional Review Board.

CODEX data collection

CODEX Antibody Generation and Validation

Oligonucleotides were conjugated to purified, carrier-free, commercially available antibodies as previously described (Schürch et al. 2020; Kennedy-Darling et al. 2021). For validation experiments, human tonsils and non-diabetic pancreata were co-embedded in FFPE blocks so both tissues could be stained and imaged simultaneously. Each antibody in the CODEX panel was validated by co-staining with previously established antibodies targeting positive and negative control cell types. Once validated, the concentration and imaging exposure time of each antibody were optimized. The tissue staining patterning was compared to the online database, The Human Protein Atlas, and the published literature. The specificity, sensitivity, and reproducibility of CODEX staining has been previously validated (Schürch et al. 2020; Kennedy-Darling et al. 2021; Black et al. 2021; Phillips, Schürch, et al. 2021; Phillips, Matusiak, et al. 2021)

CODEX Staining

Staining and imaging was conducted as previously described (Schürch et al. 2020; Kennedy-Darling et al. 2021; Phillips, Schürch, et al. 2021; Black et al. 2021). Briefly, FFPE tissues were deparaffinized and rehydrated. Heat-induced epitope retrieval (HIER) antigen retrieval was conducted in Tris/EDTA buffer at pH9 (Dako) at 97°C for 10 minutes. Tissues were blocked for 1 hour with rat and mouse Ig, salmon-sperm DNA, and a mixture of the non-fluorescent DNA oligo sequences used as CODEX barcodes. Tissues were stained with the antibody cocktail in a sealed humidity chamber overnight at 4°C with shaking. The next day, tissues were washed, fixed with 1.6% paraformaldehyde, 100% methanol, and BS3 (Thermo Fisher Scientific), and mounted to a custom-made acrylic plate attached to the microscope.

CODEX Imaging

Imaging was conducted using the Keyence BZ-X710 fluorescence microscope with a CFI Plan Apo λ 20x/0.75 objective (Nikon). “High resolution” mode was selected in Keyence Navigator software, resulting in a final resolution of .37744 µm/pixel. The exposure times are listed in Supplemental Table 3. Regions for imaging were selected by rendering HLA-ABC, Proinsulin, and CD8 and selecting large regions (averaging 55mm²). The full antibody panel and cycle-ordering is detailed in (Supplemental Tables 3 and 4). Biotinylated hyaluronan-binding protein was rendered by adding streptavidin-PE at 1:500 concentration to the 96 well plate containing fluorescent oligos in the last cycle and running the CODEX program normally. DRAQ5 was added to the last cycle because we found it stained nuclei more evenly than HOECHST which slightly improved segmentation. Each tissue took between 3 and 7 days to image depending on the tissue area.

Characterization of Tertiary Lymphoid Structures

A serial section from case 6228 was imaged in a separate CODEX experiment using an antibody panel tailored for characterizing tertiary lymphoid structures, as described in Supplemental Table 3. These data were acquired and analyzed identically to the main dataset.

Image Pre-processing

Drift compensation, deconvolution, z-plane selection was performed using the CODEX Toolkit uploader (github.com/nolanlab/CODEX, Goltsev et al. 2018). Cell segmentation using the DRAQ5 nuclear channel and lateral bleed compensation was performed with CellSeg (Lee et al. 2022).

Cell Type Clustering and Annotation

Marker expression was z-normalized within each donor and subsequently clustered in two steps. First, cells were projected into 2 dimensions using the markers indicated in Supplemental Table 2 and Parametric Uniform Manifold Approximation and Projection (pUMAP)(Sainburg, McInnes, and Gentner 2021) was applied on a downsampled dataset. The fit model was used to transform the remaining cells. Cell types were gated using Leiden clustering and manual merging. The cluster containing immune cells was sub-clustered using the markers detailed in Supplemental Table 2. Acinar cells contaminating the Immune cluster were gated out and merged with the Acinar cluster from the previous step. The Endocrine class was sub-clustered into ⍺-, β-, and δ-Cells using Glucagon, Proinsulin, and Somatostatin respectively. Clusters were annotated according the heatmap marker expression, and overlaying annotations onto raw images (Supplemental Figure 1.A) using custom scripts in Fiji.

Islet Segmentation and Pseudotime

Preprocessing

Windows consisting of the twenty nearest spatial neighbors surrounding each single cell were clustered according to their cell-type composition using Mini Batch K Means with k=200. For this analysis, ⍺-, β-, and δ-Cells were combined into one ‘Endocrine’ cell type. One cluster was highly enriched in Endocrine cells and accurately defined the islet area. Individual islets were identified using the connected components algorithm and filtering out islets that had fewer than ten total cells.

PAGA Analysis

For each islet, the number of each cell type inside the islet and between the islet edge and 20µm beyond were extracted. To adjust for variation due to the islet size, the cell type counts were divided by the number of endocrine cells inside the islet. Data were then log-transformed. The PAGA embedding was computed using the default parameters except for the following: The neighborhood search was performed using cosine distance and 15 nearest neighbors; Leiden clustering used a resolution of 1. For computing the pseudotime values (used in the colormap in Figure 2B, the x-axis in Figure 2F, and Figure 4), the path through the inflamed islet was isolated by temporarily omitting 25 islets positioned in the middle of the map between Normal and Insulin-Depleted islets. Only 9 of these were from T1D or AA+ donors.

Quantification and Validation of Functional Marker Gating

Annotation of Ground-Truth Dataset

4000 CD8⁺T cells were labelled for fifteen markers by an immunologist familiar with the staining patterns of each marker using VGG Image annotator (Abhishek Dutta and Andrew Zisserman 2019).

Automated Thresholding

For each functional marker of interest, the lateral-bleed-compensated mean fluorescence (Lee et al. 2022) of cell types known to not express the marker in question were used to calculate a background distribution. Marker-positive cells were defined as those whose expression was greater than the 99^th percentile of the background distribution.

Gating with Neural Network

22µm x 22µm cropped images of each single cell were used as training data. The marker that the image corresponded to was not included as an input in the neural network and one classifier was trained for all markers. Cells were split into training, validation, and test splits (60/15/25 respectively). ResNet50 architecture and initial weights were imported from the Keras library pre trained on ImageNet. Image augmentation consisted of random flips, rotations, zooms, contrast, and translation (+/- ten pixels only). All weights were unfrozen, and the model was trained for 100 epochs (see accompanying source code for training details).

Sub-clustering of Inflamed Islets with Cell-Type specific Functional Markers

For each Inflamed Islet (n=351), the frequency of each marker expressed by CD8⁺T cells was computed. Single cells inside the islet and within 20µm from the islet’s edge were combined before the frequency was measured. The subsequent matrix underwent z-normalization followed by UMAP gating using Bokeh. Insulin-Depleted + Immune Islets were defined as islets without β-cells with greater than two CD8⁺T cells and greater than seven macrophage/DCs. These thresholds correspond to the 95^th percentiles of CD8⁺T cells and macrophage/DCs in Normal islets.

Identification of Cellular Neighborhoods

Previously, CNs (Schürch et al. 2020) were identified by, for each single cell, defining its ‘window’ as the 20 spatial nearest neighbors. Cells were clustered according to the number of each cell type in their windows using Mini Batch K-Means. The output clusters corresponded to CNs. To ensure our method was sensitive to rare neighborhoods, we adapted this algorithm by intentionally over-clustering, using k=200 in the K-Means step rather than using a k ranging from 10-20 as used elsewhere (Bhate et al. 2021; Phillips, Matusiak, et al. 2021; Shekarian et al. 2022). Next, to determine which cell types were characteristic of each cluster, we identified, for each cluster, the set of cell-types that were present in more than 80% of the windows allocated to that cluster. We named the clusters according to this set of cell-types and merged all clusters with the same name, resulting in seventy-five CNs. Acinar cells and epithelial cells were used in the kNN graph and in the clustering but were not considered when merging clusters. Note that this method does not differentiate neighborhoods that have the same combination of cell types but different stoichiometries.

Lobule Segmentation

A training dataset was generated by manually tracing the edges of lobules in ImageJ using the ROI function. The ROI were then floodfilled in Python and used as masks for training. For each tile, the blank cycle was selected to distinguish tissue from background coverslip. A U-Net model was trained for 10 epochs (see attached source code for training details). After stitching together all masks, the resulting images required slight refinement where lobules were not completely separated, and this was done manually in ImageJ. The connected components in the stitched image defined the lobule instances. Cells were assigned to a lobule by indexing the lobule mask with their X and Y coordinates. Cells in the inter-lobular space were assigned to one “edge” lobule. This resulted in 464 lobules.

Formulation of Hierarchical Linear Models

For each lobule, the number of each cell type in the extra-islet space was divided by the number of acinar cells in the extra-islet space. For all HLMs, the lme4 package for R was used (Bates et al. 2015) and statistical significance was computed using the lmerTest package for R (Kuznetsova, Brockhoff, and Christensen 2017). Lobular cell-type abundance was z-normalized within each donor and the pseudotime was z-normalized across the entire dataset prior to fitting.

The ICC was computed using the model: pseudotime_islet ∼1|lobuleID with the performance package in R. A value of 0 indicates that the variation in pseudotimes of islets within the same lobule is equal to the variation across all islets in the donor and a value of 1 indicates that the variation in pseudotimes of islets within the same lobule is much smaller than that of all islets in the donor.

For each cell type, a two-level, random intercept HLM within each donor was constructed with the following formulation (in R formula syntax): pseudotime_islet ∼ celltype_lobule + (1|lobuleID) and a three-level random intercept, random slope HLM including islets from all donors was formulated: pseudotime_islet ∼celltype_lobule +(1+celltype_lobule|donorID)+(1|lobuleID). Here, pseudotime_islet equals the pseudotime of each islet, celltype_lobule equals the number of the given cell type in a particular lobule divided by the number of acinar cells in that lobule, z-normalized within each donor, and lobuleID and donorID are categorical variables specifying the lobule and donor that the given islet belongs to.

Neighborhood Adjacency

The adjacency between neighborhoods was computed as in (Bhate et al. 2021). The only modification was that neighborhood instances were identified using connected components of the k-NN graph with k=5 rather than from the thresholded images.

Data Availability

All primary code is attached with the submission files and will be posted on the Nolan Lab github upon acceptance. Processed single cell dataframes are accessible in the dropbox link provided in the Key Resources Table. We are in communication with the Human Cell Atlas Data Portal for hosting the raw imaging data which is on the order of 20TB. Currently, we can provide access to raw data but cannot do so anonymously.

Data Availability

Processed single cell dataframes are accessible at: https://www.dropbox.com/sh/8sqef3iwjb6f4wp/AAAQzFgDOgMBiWh9Jqcb81U8a?dl=0.

The raw and stitched, processed data is hosted by the BioImage Archive (https://www.ebi.ac.uk/bioimage-archive/) with the accession number S-BIAD859.

Code Availability

All primary code is attached with the supplementary information and will be posted on the Nolan Lab’s github upon acceptance. The code is completely open source. The code is not intended as a software tool and so no small example data set is applicable.

Acknowledgements

This research was performed with the support of the Network for Pancreatic Organ donors with Diabetes (nPOD; RRID:SCR_014641), a collaborative type 1 diabetes research project supported by JDRF (nPOD: 5-SRA-2018-557-Q-R) and The Leona M. & Harry B. Helmsley Charitable Trust (Grant#2018PG-T1D053, G-2108-04793). The content and views expressed are the responsibility of the authors and do not necessarily reflect the official view of nPOD. Organ Procurement Organizations (OPO) partnering with nPOD to provide research resources are listed at http://www.jdrfnpod.org/for-partners/npod-partners/. Research reported in this publication was supported by the National Cancer Institute and National Institute of Allergy and Infectious Diseases of the National Institutes of Health under Award Numbers K99CA246061, 5U54CA209971-05, 5U2CCA233195-02, 1U2CCA233238-01, 5U2CCA233195-02, 5U01AI101984-09. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. C.M.S. was supported by the Swiss National Science Foundation (P300PB_171189, P400PM_183915). We like to thank Yury Goltsev, Pauline Chu, Sarah Black, Gustavo Vazquez, Aviv Hargil (Stanford University), and Irina Kusmartseva (nPOD) for excellent assistance. We like to thank Dr. Xavier Rovira-Clavé (Stanford University) for critical comments on the manuscript.

Author information

G.L.B curated data, wrote software, and analyzed data. D.P. and C.M.S contributed significantly to validating the antibody panel and assisted with analysis. S.D., S.S.B, A.Y., H.A.M., G.K.K., N.N., S.R., and J.M. assisted with analysis. G.L.B., J.A.B., G.P.N., and P.L.B., conceptualized the study and wrote the manuscript. All authors have read and approved the final version of the manuscript.

Ethics declarations

P.L.B.: Founder, Halo Biosciences.

N.N.: Founder, Halo Biosciences.

PLB, NN and GK have filed intellectual property around 4-MU. PLB, NN and GK hold a financial interest in Halo Biosciences, a company that is developing 4-MU for various indications.

G.P.N. has received research grants from Vaxart and Celgene during the course of this work and has equity in and is a scientific advisory board member of Akoya Biosciences. Akoya Biosciences makes reagents and instruments that are dependent on licenses from Stanford University. Stanford University has been granted US patent 9909167, which covers some aspects of the technology described in this paper.

J.A.B.: Board of director for Gilead and CEO and President of Sonoma Biotherapeutics; scientific advisory boards of Arcus Biotherapeutics and Cimeio Therapeutics; consultant for Rheos Medicines, Provention Bio; stockholder in Rheos Medicines, Vir Therapeutics, Arcus Biotherapeutics, Solid Biosciences, Celsius Therapeutics; Gilead Sciences, Provention Bio, Sonoma Biotherapeutics.

C.M.S.: Scientific advisory board of, stock options in, research funding from Enable Medicine, Inc.

Supplemental Information

Supplemental Figure 1 Related to Figure 1.

Supplemental Figure 2 Related to Figure 2.

Supplemental Figure 3 Related to Figure 3

Supplemental Figure 4 Related to Figure 4

Supplemental Figure 5 Related to Figure 5

Supplemental Table 1:

Supplemental Table 2:

Supplemental Table 3:

Supplemental Table 4: