The virtual thymus model enables spatiotemporally resolved clonal analysis.

(A) Schematic overview highlighting 6 key steps during T-cell development, which are integrated into the virtual thymus model: (1) ETPs enter the thymus niche; (2) they receive Dll4 and IL-7 signals from the TECs. Triggered by these signals, (3) thymocytes undergo proliferation and (4) differentiation into two distinct T-cell sublineages (5). Finally, thymocytes undergo selection, either leaving the organ or undergoing cell death (6). (B-C) Experimentally measured (“Medaka”) and computationally modeled (“Virtual”) cell speed (B) and cell directionality (C) for ETPs and thymocytes. Bars are mean values, and error bars standard deviation. (D) Timeline of medaka embryonic development (top) compared to simulated time (bottom). As indicated by the cartoon at the bottom, the first ETP enters the thymus roughly at 2.5 dpf (t = 0 h in the simulation). At roughly 60 h (roughly 5 dpf), the thymus reaches its homeostatic cell population, with both cell entry and exit. (E) Representative snapshots of the simulation at homeostasis showing spatial expression patterns and signaling activity of components of IL7R and Notch1 signaling pathways in TECs (top two panels) and thymocytes (bottom 4 panels). Data are shown in arbitrary units normalized to a scale from 0 to 1. (F) Left: Cartoon explaining the asynchronous arrival (homing) of new cells into the organ, their clonal expansion, and asynchronous cell death or exit. Right: Visualization of clonal diversity in a typical simulation; each color shade uniquely indicates a clonal lineage. The height of each colored patch indicates the number of cells in that clone, the length indicates simulation time. For illustrative purposes, a black arrow indicates the homing of a sample of lineages. (G) Data from the right panel in (F) was normalized to time of homing. This visualization highlights developmental phases of thymocyte lineages. Gray arrowheads indicate rounds of cell division. Scale bar 50 cells. Abbreviations: ETP = early thymic progenitor, T = Thymocyte, N = number of replicates (biological or computational).

TEC density influences thymocyte proliferation via IL7R signaling.

(A) Left: Detail of confocal section highlighting the tight spatial interaction between TECs and thymocytes (T). Scale bar 5 µm. Right: Detail from a three-dimensional render of a simulation snapshot; TECs cyan, thymocytes green. (B) Parameter permutations tested. (C) Impact of parameter permutations on average homeostatic cell population. Bars indicate means, error bars standard deviation. (D) Clonal lineages normalized to time of thymic entry for reference (left), a scenario with low average population (middle) and a scenario with high average population (right). Scale bar: 100 cells. (E) Extracellular IL-7 gradient in simulations; sum projection of all z-planes of a simulated confocal stack. Units are arbitrary. (F) Mean extracellular IL-7 concentration in the entire simulated volume. Note that these values are lower than in (E) because of averaging over the volume, including empty space (black areas in (E)). (G) Thymocyte population size over time averaged over several simulation runs. Shaded area indicates standard deviation. (H) Illustration of experimental setup for IL-7KO and TEC:IL-7HI. At least 3 technical replicates were done for the injection of each construct. (I) Representative images of embryos stained with GFP (green) and pH3 (red) of WT, IL-7KO and TEC:IL-7HI. Scale bar 20µm. (J) Numbers of pH3 positive cells per thymus normalized to mean of WT. N represents the number of individual fish. Data show means ± standard deviation.

High IL7R signaling promotes proliferation.

(A) Illustration of scenarios without extracellular IL-7 depletion (left) and with extracellular IL-7 depletion (right). (B) Ordinary differential equation used to model the IL7R signaling activity. Parameters that were modified are highlighted. (C) Effect of IL7R signaling-related parameter permutation. Bars indicate means and error bars standard deviation. (D) Left: Mean extracellular IL-7 concentration in the entire simulated volume. Right: Extracellular IL-7 gradient in simulations; sum projection of all z-planes of a simulated confocal stack. Units are arbitrary. The introduction of IL-7 depletion by thymocytes leads to a slight reduction in extracellular IL-7 availability.

Systematic screen shows that lower TEC density promotes malignant expansion.

(A) Schematic representation of the major lesions to IL7R that were implemented in the model. (B) Representative clonal plot highlighting in black the lesioned clone. (C) We defined clone size as the maximum number of terminal leaves of a lineage, regardless of eventual fate of each cell. Clone size distribution in the reference simulation (REF), compared to IL7RDA and IL7RHI. REF: N=20 simulations, with 2088 clones total: IL7RDA: N=13 simulations with 1330 non-lesioned clones and 13 lesioned clones. IL7RHI: N=13 simulations with 1336 non-lesioned clones and 13 lesioned clones. (D) Schematic of tested scenarios shown in (E). (E) Simulated permutations ordered along increasing mean lesioned clone size for the IL7RHI condition. Expected effect on proliferation is a qualitative estimate based on preliminary simulations. A total of 1580 permutations were tested with a varying amount of replicates per permutation. Error bars indicate standard deviation. Note that in (C) and (E) the axis is in base-2 logarithmic scale to better illustrate rounds of cell division.

Interplay of IL7R and NOTCH1 potentiates clonal expansion in vivo.

(A) Schematic of Notch1 and IL7R signaling and their regulation of thymocyte cell cycle entry, proliferation, and differentiation as implemented in the virtual thymus model. (B) Simulated data; rounds of cell division of normal and lesioned clones with different types of lesions in NOTCH1 or IL7R signaling. Note that the y-axis is base-2 logarithmic. Green boxes: Condition/lesion present. Gray boxes: Condition/lesion absent. (C) Illustration of plasmids for DNA micro-injection for transient overexpression of the different genes. Note that the expression of GFP was used to select successfully injected embryos. (D) Representative image: (left) of the WT control at 9 dpf where thymocytes are labeled in green with normal thymus size; (middle) of a thymus hyperplasia phenotype at 7 dpf; (right) same embryo developed T-ALL phenotype at 9 dpf, where infiltration in other organs was detected, namely in central nervous system (CNS) and gut. (E) Percentage of phenotypes: normal thymus, thymus hyperplasia and T-ALL at 11dpf. For each construct injection, we did at least three technical replicates. For detailed statistics see Figure 5 - Supplement Table 1. Statistical data rounded to 4 decimal places. (F) Number of pH3 positive cells per thymus counted during confocal imaging. (G) Representative confocal images of immunostaining against pH3 (red) and GFP (green). Scale bar is 20µm. N represents the number of biological samples.

IL-7 supplied by the niche can lead to thymus hyperplasia.

(A) Simulated data; rounds of cell division of normal and lesioned clones for different combinations of lesions and TECs expressing IL-7 ubiquitously. Note that some conditions are reproduced from Figure 5B for easier comparison. (B) Left panel: Frequency of thymus hyperplasia and T-ALL phenotype observed in the injected embryos at 11 dpf. For detailed statistics see Figure 5 - Supplement Table 1. Statistical data rounded to 4 decimal places. Note that some data from Figure 5E are repeated here to facilitate a better comparison. N represents the number of biological samples. Right panel: representative images at 10 dpf after co-injection of lck:gfp-NICDΔPEST-mycn and ccl25a:tagRFP-il-7 constructs. Arrows indicate GFP-expressing malignant thymocytes in the thymus, brain, and gut. Scale bar is 400 µm (C) We propose a hypothetical negative feedback loop from proliferation to differentiation, which could explain discrepancy between in silico predictions and in vivo observations.

The virtual thymus represents a slice of the organ.

Three-dimensional renders of the simulation. (A) Showing only TECs in cyan. (B) Showing both thymocytes (bright green) and TECs (cyan). (C) Detail showing clonal thymocyte lineages. Each lineage was assigned a randomly selected color, TECs in dark gray. In homeostasis, most clones consist of 1-8 cells.

Simulation time-lapse.

Three-dimensional renders of a typical simulation where each thymocyte lineage was assigned a randomly selected color. TECs are shown in dark gray. Simulation runs for 48000 steps, where 1 step is 15 seconds. Note how thymocytes can crawl over and under TECs.

Simulation time-lapse with fewer cells to highlight clonal dynamics.

In this simulation, a total of eight ETPs were introduced into the simulation at the same time to highlight clonal expansion and motility within the organ (00:10-02:00s), and later cell exit and cell death (starting at 02:08s). Three-dimensional renders of a simulation where each thymocyte lineage was assigned a randomly selected color. TECs are shown in dark gray. Each simulation step is 15 seconds.

Parameters varied in this work

For a full list of parameters and their values, please refer to Aghaallaei et al. 2021.

Statistical analysis of simulated data.

(A) Three-dimensional rendering of TECs in the simulation with the indicated parameter alterations. Orange: IL-7-secreting TECs. Cyan: Non-secreting TECs. Dotted outline indicates thymus outline in the simulation. (B) Partial residual effect plot showing linear model fit on log2-transformed thymic population size. Shaded area is 95% confidence interval. (C) Estimated coefficients and p-values of the linear model. Significance codes: <0.001 “***”; <0.01 “**”; <0.05 “*”. All parameters have a statistically significant effect on cell number. Parameter interactions are also statistically significant, which indicates that changing combinations has a greater effect than the sum of individual changes.

Reduced thymus volume upon IL-7KO.

(A) Experimental design: gRNA targeting the signal peptide of il-7 was co-injected with Cas9 into the 1-cell-stage blastomere of Tg(lck:gfp) medaka. The embryo was raised until 8 dpf. Whole-mount immunostaining was performed using anti-GFP and M-phase marker phospho-histone 3 (pH3). The thymus was imaged using confocal microscopy, and embryos were genotyped by sequencing. This panel was created using BioRender.com. (B) Schematic of the il-7 gene structure and binding site for the crRNA for CRISPR-Cas9 gene editing in the signal peptide region of exon 2. Zoom in: Representative sequencing data analyzed with decodr.org (Bloh et al. 2021) software. (C) Representative image of immunostaining of Tg(lck:gfp) in WT and IL-7KO and evaluation of the thymus size using the area measurement tool in Imaris software (yellow). Scale bar 20 µm. (D) Quantification of thymus size: Thymus area in µm2 in WT and IL-7KO. N represents the number of biological samples. For statistics, unpaired t-test with Welch’s correction was used; Data show means ± standard deviation. Significance threshold was set to p= 0.05. ** means p < 0.01.

Statistical analysis of simulated data.

(A) Partial residual effect plot showing linear model fit on log2-transformed thymic population size. Shaded area is 95% confidence interval. (B) Estimated coefficients and p-values of the linear model. Significance codes: <0.001 “***”; <0.01 “**”; <0.05 “*”. Although the population size is slightly reduced for IL-7 depletion, the effect is not statistically significant. In contrast, IL7R signaling activation and deactivation have a statistically significant effect. None of the interactions have statistically significant effects.

Clonal lesions and their effects

Source data file for visualization in Figure 4 E.

In the data file, each row corresponds to the log2-transformed average clone size for one cell genotype (normal or lesioned) from one single independent replicate simulation. Additionally, the parameter values used in that run are listed in separate columns. The “grouping” column corresponds to the simulated scenarios (parameter permutation). The “simulation_datafile” column is the unique identifier of a given simulation run.

Simulation time-lapse with lesioned clone.

Two-dimensional projections of the virtual thymus volume showing time-lapse of simulations representing conditions in Figure 4B. Each thymocyte lineage was assigned a unique color; the lesioned clone was colored in black. Representative position of TECs shown in light blue (non-secreting) and light orange (IL-7-secreting). Movies were synchronized to the time of entry of the first thymocyte.

Modulation of proliferation-related parameters affects lesioned thymocytes more strongly.

(A) Clones with IL7RHI lesion. (B) Clones with IL7RDA lesion. (C) Control simulations without lesioned clones. Note that the y-axis of all plots is base-2 logarithmic to better illustrate rounds of cell division. Error bars indicate standard deviation. By modulating proliferation parameters, lesioned clones gain up to 6 rounds of division, while non-lesioned clones gain only 2 rounds of division. Expected effect on proliferation is a qualitative estimate based on preliminary simulations.

Simulation time-lapse.

Two-dimensional projections of the virtual thymus volume showing time-lapse of simulations representing conditions in Figure 4 – Supplement 1A-C. Specifically, left column represents the condition with reduced proliferation, and the right column the condition with increased proliferation. Note how, with increased proliferation, newcomer cells struggle to enter deeper into the organ due to cell crowding. Each thymocyte lineage was assigned a unique color; the lesioned clone was colored in black. Representative position of TECs shown in light blue (non-secreting) and light orange (IL-7-secreting). Movies were synchronized to the time of entry of the first thymocyte.

Modulation of parameters affecting IL7-signaling have almost no effect on lesioned clones.

(A) Clones with IL7RHI lesion. (B) Clones with IL7RDA lesion. (C) Control simulations without lesioned clones. Note that the y-axis of all plots is base-2 logarithmic to better illustrate rounds of cell division. Error bars indicate standard deviation. By modulating IL-7 signaling parameters, lesioned clones do not gain extra rounds of division, while non-lesioned clones gain up to 3 rounds of division. Expected effect on proliferation is a qualitative estimate based on preliminary simulations.

A secondary lesion delaying differentiation doubles proliferative potential.

(A) Clones with IL7RHI lesion. (B) Clones with IL7RDA lesion. (C) Control simulations without lesioned clones. Note that the y-axis of all plots is base-2 logarithmic to better illustrate rounds of cell division. Error bars indicate standard deviation. Secondary clonal lesions lead to a gain of 4 rounds of division in the lesioned clones, with no effect on non-lesioned clones. Expected effect on proliferation is a qualitative estimate based on preliminary simulations. (D) Mean extracellular IL-7 concentration over time for selected simulations shown in (A)-(C). The increase in extracellular IL-7 by autocrine activity of lesioned clones is not enough to promote substantial proliferation of non-lesioned cells.

Parameters affecting TEC distribution have opposing effects on lesioned and non-lesioned clones.

(A) Clones with IL7RHI lesion. (B) Clones with IL7RDA lesion. (C) Control simulations without lesioned clones. Note that the y-axis of all plots is base-2 logarithmic to better illustrate rounds of cell division. As TEC architecture is made denser, lesioned clones lose up to 2 rounds of cell division, while non-lesioned clones gain 2 rounds of cell division. Expected effect on proliferation is a qualitative estimate based on preliminary simulations.

Sparse TEC distribution is among permutations that most enabled lesioned clone expansion.

Shown are top 8 conditions from Figure 4 (E). (A) Clones with IL7RHI lesion. (B) Clones with IL7RDA lesion. (C) Control simulations without lesioned clones. Note how lesioned clones with autocrine IL-7 production favor expansion of non-lesioned clones as well. The y-axis of plots in (A)-(C) is base-2 logarithmic to better illustrate rounds of cell division. Points are mean values, error bars indicate standard deviation. Expected effect on proliferation is a qualitative estimate based on preliminary simulations. (D) Volume occupied by thymocytes (1 thymocyte ≅ 33.5 μm3) of a hypothetical lesioned clone compared to the volume of the virtual thymus (≅ 7854 μm3). The y-axis is scaled in millions of cubic micrometers. The volumes strongly diverge after 12 rounds of cell division indicating extremely dense cell packing that would result in hyperplasia or metastasis in vivo.

Increased il7 expression upon NICDΔPEST-mycn overexpression.

(A) Schematic depiction of NOTCH1 signaling (left) and overexpression of NICD (right). The engagement of NOTCH1 receptor with its ligand (DLL4) leads to its cleavage, and release and translocation of the NICD in the nucleus. NICD through interaction with other co-factors and triggers the expression of the target genes. Overexpression of NICD leads to activation of NOTCH signaling independent of the ligand. (B) Relative expression of il7r and il7 of isolated thymocytes from WT compared to those expressing both NICDΔPEST and mycn. Each dot represents data from an individual sorted embryo. Data are log10 means ± standard error of the mean. Statistics: unpaired Mann-Whitney test, significance threshold was set to 0.05. (C) WISH with the il7 probe of freshly hatched yolk sac WT and NICDΔPEST and mycn larvae. Scale bare is 100µm. (D) Pairwise sequence alignment was done using EMBOSS Needle and the Needleman-Wunsch algorithm (Madeira et al. 2022), highlighting the partial medaka IL7R protein sequence and IL7RNPC insertion mutation of amino acids NPC at position 266 in the extracellular juxtamembrane-transmembrane interface. (E) Relative lck expression quantified by qPCR of sorted WT thymocytes compared to those in IL7RHI and IL7RNPC. The lck expression was normalized to the housekeeping gene ef1a. Each dot represents an individual sorted embryo. Data are log10 means ± standard error of the mean. Statistics: unpaired Mann-Whitney test. Significance threshold was set to p=0.05. * means p < 0.05.

Detailed statistics for phenotype comparison using Fisher’s exact test and pairwise comparisons; adjusted p values were calculated with the Benjamini & Hochberg method. Values reported in figures are the adjusted p values. Significance threshold was set to 0.05. Data rounded to 4 decimal places.