Altered thymic niche synergistically drives the massive proliferation of malignant thymocytes

The details of the virtual thymus model implementation and parameter estimation have been extensively described previously (see Supplementary Material of Aghaallaei et al., 2021), and in the following, we will only delineate its main features. The model considers both spatial and temporal aspects of early thymus development at multiple scales: Cells are physically represented as one or multiple spheres that mechanically interact with adhesive and repulsive forces to capture cell crowding; each cell has an independent internal state that changes dynamically based on subcellular signaling modeled with ordinary differential equations or with phenomenological rules; and at the tissue level, a partial differential equation is used to solve the diffusion of extracellular molecules such as cytokines. To reduce computational cost, the computational representation focuses on a 5 µm deep slice of the lower half of the radially symmetric organ, approximately 1/10 of the total volume of a medaka embryonic thymus (Figure 1—figure supplement 1A and B). Consistent with observations in vertebrates, including in medaka (Bajoghli et al., 2009; Bajoghli et al., 2015; Aghaallaei et al., 2021; Aghaallaei et al., 2022), ETPs within our virtual thymus mechanically and biochemically interact with TECs, proliferate and differentiate in response to signals from the thymic niche, and subsequently diverge into two distinct T-cell sublineages before selection and exit from the thymus (Figure 1A; see Appendix 1 and Aghaallaei et al., 2021, for more technical details). In this work, we define cells of a lineage derived from a single founder ETP as a ‘clone’, and we do not consider differences in immunological clonotypes that emerge after recombination of the T-cell receptor.

Figure 1 with 3 supplements see all

Download asset Open asset

In the virtual thymus model, factors governing thymocyte motility, including cell speed (Figure 1B) and directionality (Figure 1C), were fine-tuned based on quantitative noninvasive imaging of the thymus using multiple medaka transgenic reporter lines (Bajoghli et al., 2009; Bajoghli et al., 2015; Aghaallaei et al., 2021). Based on these experimental observations, we modeled entry into the organ by creating new cells at the bottom of the simulated domain (ventral to the organ). These cells move into the thymus and, some time after differentiation, either leave the organ by exiting from the ventral side as a naïve T cell or undergo cell death by gradually shrinking and eventually being removed from the simulation (Figure 1—videos 1–2). The model accounts for the time taken for a single progenitor to enter and exit the thymus, which is estimated at approximately 3 days in zebrafish (Hess and Boehm, 2012) and medaka (Bajoghli et al., 2015). For time calibration, each simulation step was set to represent 15 s. With a total of 48,000 steps per simulation, the entire simulation spans 720,000 s, equivalent to about 8.3 days (Figure 1D). This time interval was chosen to recapitulate thymic development in medaka from 2.5 to 11 days post-fertilization (dpf), because at 2.5 dpf the first ETP enters the embryonic thymus (Bajoghli et al., 2019). Correspondingly, each simulation starts with entering the first ETP into the virtual thymus and ends with a homeostatic cell population of roughly 100 cells (mean ± standard deviation: 98±7; Figure 1—video 1). At 11 dpf, the medaka thymus contains 877±146 (N=7) thymocytes. Considering that the virtual thymus represents 1/10 of the medaka embryonic thymus, the cell population in the virtual thymus model reproduces the temporal dynamics in medaka within the margin of error.

Each cell’s internal state and decision to proliferate or differentiate depends on intrinsic and extrinsic factors that are integrated by signaling pathways. The virtual thymus includes Delta-like 4 (Dll4), the ligand of Notch1 receptor, and interleukin-7 (IL-7) cytokine as two cell-extrinsic factors provided by the TECs (Figure 1E). Consistent with medaka WT embryonic thymus (Aghaallaei et al., 2021; Aghaallaei et al., 2022), all virtual TECs in this model express the Dll4 (Figure 1E, top-right panel), with a subset spatially releasing IL-7 into the environment (Figure 1E, top-left panel), creating a short-ranged IL-7 cytokine gradient in the thymic cortex niche (Aghaallaei et al., 2021). In terms of cell-intrinsic factors, all ETPs uniformly express the Notch1 receptor (Figure 1E, middle-right panel), and each cell lineage expresses the IL-7 receptor (IL7R) at a constant level chosen randomly between 0 and 1 (Figure 1E, middle-left panel). In each simulation, new ETPs continually enter the thymic niche at regular intervals. Following engagement with TECs and receiving Notch1 and IL-7 signals (Figure 1E, bottom panels), they proliferate and differentiate (step 3 in Figure 1A). The cell cycle was modeled with an average duration of 7 hr and was subdivided into G1, S, G2, and M phases. In line with in vivo observations (Ruijtenberg and van den Heuvel, 2016), virtual cells require external pro-proliferative signals during the G1 phase to commit to S through M. Once differentiated, cells can no longer enter the cell cycle. The dynamics of pro-proliferative signals and the proliferation stop after differentiation allow clones in silico to undergo up to 4 rounds of cell division (Aghaallaei et al., 2021). Thymic selection is modeled as a probabilistic event and is not regulated by signaling dynamics.

In this work, we implemented a code that enables us to meticulously trace the unique identity of each cell upon entering the virtual thymus and subsequently track the identity of its descendants (Figure 1F, left panel; Figure 1—figure supplement 1C). This tool identified various waves of clones because of the constant influx of new ETPs and negative selection and efflux of fully differentiated T cells (Figure 1F, arrows in the right panel, Figure 1—video 1). We note a slight variability in the clone sizes, which could arise from three processes. First, random variations in the expression level of IL7R impact IL-7 signal transduction and, consequently, proliferation (Aghaallaei et al., 2021). Second, due to random cell motility and cell crowding effects, certain cells could coincidentally remain in close contact with IL-7-secreting TECs, maximizing their exposure to pro-proliferative signals such as IL-7 and Notch ligand, and thus increasing their chances of entering the cell cycle. Third, each cell determines the duration of its next cell cycle by drawing an Erlang-distributed random variable with a mean of 7 hr and a standard deviation of ≅0.99 hr.

After normalizing all clones based on their time of entry into the thymus, the pattern of developmental phases became clear across all clones (Figure 1G). In this visualization, the now near-synchronous rounds of cell division across clones can be seen as bumps in the graph indicating clonal expansion (Figure 1G, gray arrowheads). During the commitment phase, which has a minimum duration of 24 hr but can be extended due to insufficient Notch signaling, thymocytes were competent for proliferation, leading to an increase in the total cell number per clone. In the differentiation phase, which has a fixed duration of 24 hr, clones reached a maximum size. Note that differentiation prevents further entry into the cell cycle but permits cells that passed G1 and therefore committed to S through M phases to complete their division. Finally, clones underwent selection and exit from the virtual thymus, thus leaving the simulation. Together, the enhanced virtual thymus model facilitated the investigation of the heterogeneity and temporal dynamics of individual clones originating from a founder progenitor cell.

In vivo, individual TECs exhibit a distinctive star-shaped morphology, and protrusions of neighboring TECs contact each other forming a three-dimensional network (Figure 2A, left panel). Using confocal imaging, we estimated the relative area occupied by thymocytes (53 ± 3%; N=4) and TECs (47 ± 3%; N=4). A similar morphology and tissue density was replicated in our in silico model (Figure 2A, right panel), with thymocytes taking up on average 56±2% and TECs 44±2% of the volume in homeostasis. We wondered to what extent TEC morphology and density could impact thymocyte population size. This aspect is difficult to study experimentally and is, therefore, an ideal use case for the virtual thymus model. In a series of simulations, we combinatorially varied (1) the size of TECs, (2) the number of their protrusions, and (3) the TEC cell density (Figure 2B). The average homeostatic thymocyte population size in each simulation was then used as a readout. The results from 26 different tested conditions predicted that a higher density of larger TECs with more protrusions led to an almost twofold increase in the number of thymocytes (200±13, N=3; scenario 26 in Figure 2C, Figure 2—figure supplement 1A, bottom) compared to the reference condition (98±7, N=19; reference in Figure 2C). Conversely, scenarios where these parameters were reduced, particularly a reduced TEC size, also diminished the homeostatic thymocyte population size (scenarios 1–14 in Figure 2C). The condition with the least amount of cells, scenario 1, had a reduced TEC size, but an unchanged number of protrusions and density. The condition diametrically opposite to scenario 26 was scenario 6, with a reduced TEC size, reduced number of protrusions, and reduced TEC density (Figure 2—figure supplement 1A, top). In both scenario 1 (45±1, N=3) and scenario 6 (51±2, N=3), the thymic population size was only about half of that in the reference condition (98∓7, N=19). Statistical analysis confirmed that all three parameters and their interactions were significant predictors of total thymocyte population size (Figure 2—figure supplement 1B and C). The statistically significant interaction can be explained by the fact that higher TEC density increases the number of TECs and thereby amplifies the effect of an increased TEC radius or number of protrusions. Thus, these three alterations of the thymic niche morphology act synergistically to modulate thymocyte population size.

Figure 2 with 2 supplements see all

Download asset Open asset

There are two possibilities that could explain the variations in total thymocyte numbers in simulations: either over- or under-proliferation of a few thymocyte lineages or the cumulative impact of subtle changes across several thymocyte clones. To distinguish between these possibilities, we used our clonal analysis tool for scenarios 6 and 26 (Figure 2D). This analysis revealed that variations in total thymocyte numbers were a consequence of changes in proliferation across all clones to a similar degree. These results predict that the TEC architecture has a direct impact on thymocyte proliferation rate, which fits well with the fact that TECs act as the main source of ligands, growth factors, cytokines, and chemokines (Gameiro et al., 2010). Any changes in the TEC architecture could therefore influence the amount of cytokine production. In our virtual thymus model, this affects the extracellular spatial distribution of IL-7 cytokine (Figure 2E). To further explore this, we tested two additional scenarios. In scenario 27, TECs maintained their reference architecture and density but all expressed IL-7, leading to a more uniform cytokine distribution in the thymic environment (Figure 2E). The extracellular level of IL-7 (Figure 2F) and thymocyte population size (155±11, N=5; Figure 2G) in scenario 27 was similarly high as in scenario 26, where the size, protrusion number, and density of TECs were increased. Conversely, when TECs did not express IL-7 (scenario 28), the thymocyte population size (Figure 2G) was low (42±1, N=5), akin to scenario 6, where the size, protrusion, and density of TECs were reduced. Interestingly, although the average extracellular IL-7 level is at its highest in scenario 27 with ubiquitous expression, the spatially restricted elevated IL-7 concentration in scenario 26 is more effective at driving thymic population growth (Figure 2E–G). Thymocytes first enter the organ at its periphery from below (at the ventrolateral site), where IL-7 levels are highest, and tend to migrate to the IL-7-depleted center of the thymus as they become non-proliferative (Aghaallaei et al., 2021; see also Figure 1—video 1). Thus, in scenario 26, thymocytes in their proliferative phase immediately encounter elevated IL-7 at the organ periphery, stimulating cell proliferation. In contrast, in scenario 27, thymocytes that enter the organ periphery encounter comparable levels of IL-7 to the reference scenario, and only after entering deeper into the organ can the ubiquitous expression of IL-7 unfold its effects. This delay, contingent to random cell migration (which is further exacerbated by cell crowding posing an obstacle), reduces the impact of elevated IL-7 in scenario 27. These results highlight how the spatial structure of the thymic niche can impact thymocyte population dynamics.

Our simulations predict that the availability of IL-7 provided by TECs has a direct impact on thymocyte population size. To validate this prediction in vivo, we conducted two functional analyses using the medaka model organism (Figure 2H). In one experiment, employing the CRISPR-Cas9 technique, we knocked out the il7 gene (Figure 2—figure supplement 2A and B) in a medaka transgenic line, where thymocytes express green fluorescent protein (GFP). Consistent with our in silico outcomes, il7 crispant embryos displayed fewer thymocytes and a smaller thymus size (Figure 2I, Figure 2—figure supplement 2C and D). To estimate the extent of cell proliferation, we counted mitotic cells throughout the entire thymus using the M phase marker phospho-histone 3 (pH3) and normalized values to the mean of the control. This analysis further supported that a reduction in IL-7 availability in the thymus reduces thymocyte proliferation (Figure 2I and J). In a second experimental setup (Figure 2H), we artificially increased il7 levels in the thymic niche by injecting the ccl25a:il7 construct into embryos (hereafter called TEC:IL-7^HI), resulting in a strong upregulation of the levels of il7 produced by the TECs (Aghaallaei et al., 2021). Compared to WT, the thymus of TEC:IL-7^HI embryos appeared visibly larger, and the number of pH3 positive cells was significantly increased (Figure 2I and J). Therefore, the in vivo results confirm the critical role of IL-7 for thymocyte proliferation in medaka, a mechanism that is evolutionarily conserved among vertebrates (Iwanami et al., 2011). Additionally, the alterations in thymic size suggest a regulatory crosstalk between thymocyte proliferation and TEC architecture.

Given the important contribution of extracellular IL-7 to thymocyte population size, we next evaluated the impact of (i) extracellular IL-7 depletion (e.g. via ligand internalization), (2) the rate of signal transduction activation upon IL7R and IL-7 binding, and (3) the rate of decay of signaling activity (Figure 3A and B). In our virtual thymus model, the IL-7 signaling activity (σ_IL7) of a thymocyte over time (t) is modeled phenomenologically using a function involving the IL7R concentration ([IL7R]), the average extracellular IL-7 concentration (〈[IL-7_ex]〉), a parameter (a_IL-7) that scales the strength of signal transduction activation, and another parameter (d_IL-7) that scales the rate at which the signaling activity diminishes over time (Figure 3B). As before, we combinatorially perturbed these three factors, simulating 17 different scenarios, to assay their impact on thymocyte population size (Figure 3C). As expected, we observed a positive correlation between the level of signal transduction activation a_IL-7 and thymocyte numbers. Conversely, variations in signaling decay rate d_IL-7 showed the opposite effect. The rate of signal transduction deactivation models cellular short-term memory, thus a lower signaling deactivation rate indicates that cells retain their IL-7 stimulus for a longer duration. Because IL-7 signaling activity σ_IL7 promotes cell cycle entry (see Appendix 1 section X. Subcellular scale: cell proliferation model for an in-depth explanation), we expect that parameter changes that increase stimulus duration should promote cell division. Indeed, the scenario with the highest number of thymocytes occurred when signaling deactivation was low and activation was high (274±11, N=4; scenario 45 in Figure 3C). The diametrically opposite combination had the least amount of cells (42±1, N=4; scenario 29 in Figure 3C). Statistical analysis confirms that signal transduction activation a_IL-7 and the signal transduction deactivation rate d_IL-7 were both statistically significant predictors of thymocyte population size (Figure 3—figure supplement 1). There was no statistically significant effect of parameter interactions, likely because these parameters exert their effect independently. There was only a limited reduction in population size with IL-7 depletion (91±11, N=4; scenario 36 in Figure 3C) compared to no depletion in our reference setting (98±7, N=19). Overall, IL-7 depletion did not have a significant effect on thymocyte population (Figure 3—figure supplement 1). However, we noted that population size difference was more pronounced at higher cell numbers, e.g., comparing between scenario 44 (237±7, N=4) and scenario 45 (274±11, N=4). Similarly, the mean extracellular IL-7 content in the virtual thymus barely changed when comparing the reference to scenario 36 but was noticeably reduced in scenario 44 compared to scenario 45 (Figure 3D). Therefore, thymocyte proliferation could self-inhibit via depletion of extracellular IL-7, but we expect this effect to be small unless cell numbers are massively increased.

Figure 3 with 1 supplement see all

Download asset Open asset

Together, the outcomes of our simulations reveal that molecular and cellular changes, specifically those capable of increasing the extrinsic factor IL-7 in the niche – such as a higher density of TECs or an elevated expression level of IL-7 by individual TECs – directly contribute to an increase in thymocyte population size. Tuning thymocytes’ sensitivity to IL-7 by manipulating IL7R signaling can further amplify proliferation. Finally, we expect that cytokine internalization by thymocytes has only a mild inhibitory effect on excessive population growth by limiting the availability of the pro-proliferative IL-7 in the extracellular space.

Our 45 tested scenarios thus far reached a maximum threefold induction of thymocyte population size when parameter changes promoted increased IL-7 signaling. While the activation of IL-7 signaling is linked to T-ALL development, this pathology is associated with proliferation of clones carrying somatic mutations (Oliveira et al., 2019; Aghaallaei et al., 2022; Silva et al., 2021). We therefore modified the code of our virtual thymus model to enable the introduction of parameter changes exclusively in a single clone and all of its progeny (hereafter referred to as the lesioned clone). Inspired by observations that 10% of T-ALL patients exhibit dominant active mutations in the IL7R gene (Liu et al., 2017; Shochat et al., 2011; Zenatti et al., 2011), and a significant subset of T-ALL patients with active IL7R signaling display elevated IL7R expression levels (Silva et al., 2021), two types of lesions were introduced: (1) a dominant active IL7R lesion (hereafter called IL7R^DA lesion), (2) an IL7R overexpression lesion (hereafter called IL7R^HI lesion). For the IL7R^DA lesion, we set the IL7R level of the lesioned clone to the maximum WT level of 1 and allowed it to activate IL-7 signaling regardless of extracellular IL-7 cytokine, mimicking the dominant activation mutations in the IL7R gene (Figure 4A, middle panel). For the IL7R^HI lesion, we mimicked overexpression of the IL7R gene by assigning a value of 10 for the IL7R level in the lesioned clone, i.e., 10-fold higher than the maximum attainable in the WT population (Figure 4A, right panel). In both scenarios, non-lesioned thymocytes retained reference IL7R levels randomly ranging between 0 and 1 (Aghaallaei et al., 2021). In each simulation, a single lesioned clone entered the thymic niche after the 60th hour of a simulation, i.e., shortly after the establishment of a homeostatic population size (Figure 4B, Figure 4—video 1).

Figure 4 with 8 supplements see all

Download asset Open asset

Figure 4—source data 1

Source data file for visualization in Figure 4E and Figure 4—figure supplement 1B.

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.

https://cdn.elifesciences.org/articles/101137/elife-101137-fig4-data1-v1.csv

Download elife-101137-fig4-data1-v1.csv

First, we evaluated the impact of IL7R^HI and IL7R^DA lesions in the reference condition, meaning no additional modifications were made to our model. The clonal analysis of simulations revealed that introducing a lesioned clone did not markedly alter the division behavior of non-lesioned clones in the same simulation (scenario with no lesion: 1.04±1.26 rounds of cell division; scenario with IL7R^DA: 0.98±1.22 rounds of cell division; scenario with IL7R^HI: 1.04±1.23 rounds of cell division; Figure 4B and C). The large variance in non-lesioned clones is consistent with our previous work (Aghaallaei et al., 2021), which showed that clones expressing endogenously high levels of IL7R undergo more rounds of cell division in the virtual thymus, whereas clones with lower levels of IL7R proliferate very little or not at all. In contrast, lesioned clones averaged just above 4 rounds of cell division (4.20±0.34 for IL7R^DA and 4.10±0.13 for IL7R^HI). Consequently, the size of IL7R^DA and IL7R^HI clones was nearly eightfold larger than an average non-lesioned clone (Figure 4C). This result suggests that lesioned clones acquired a distinct proliferative advantage over non-lesioned clones in our virtual thymus model. Nevertheless, since the expansion of lesioned clones amounted to, at best, only one additional round of cell division compared to the upper range of the non-lesioned clone distribution (Figure 4C), our virtual thymus model predicts that a single lesion in the IL7R gene will not be clinically impactful.

Therefore, we next explored whether additional modifications in both cell-autonomous and non-autonomous factors might substantially enhance the clonal expansion. To identify these factors, we decided to undertake a systematic approach and tested combinations of scenarios affecting the TEC architecture and IL-7 signaling parameters, as tested in Figures 2 and 3, together with changing other parameters affecting proliferation and differentiation of all clones (Figure 4D). Furthermore, we considered scenarios in which only lesioned clones exhibited additional modifications, such as reduced cell motility, slower differentiation time, or autocrine IL-7 production (Table 1). The inclusion of the latter was inspired by a study showing that certain malignant T cells derived from T-ALL patients possess the ability to ectopically express IL-7 cytokine (Buffière et al., 2019). In theory, there are 104,976 possible permutations of these scenarios. To reduce the space of possibilities and thus the computational cost, we prioritized alterations expected to increase proliferative potential based on outcomes shown in Figures 2 and 3. We also decided to focus on simulating the extremes; e.g., low and high levels shown in Figure 4D. In addition, we simulated a small sample of scenarios expected to decrease proliferative potential, such as a reduction in TEC density. In the end, 1580 scenarios were simulated. The clonal analysis tool was used to compare clone size and the number of cell division rounds between lesioned and non-lesioned clones under the same conditions. Overall, we observed a similar proliferation pattern among IL7R^HI and IL7R^DA clones in most of the tested scenarios. Likewise, non-lesioned clones displayed similar trends as clones in control simulations that lacked a lesioned clone (Figure 4—figure supplement 1). The lesioned clones consistently showed a higher proliferation rate than their non-lesioned counterpart (Figure 4E; Figure 4—source data 1). This trend persisted even when parameters known to generally increase cell division for all cells, such as mean cycle duration and proliferation duration, were modified (Figure 4—figure supplement 2, Figure 4—video 2). In fact, lesioned clones disproportionally gained up to 6 rounds of cell division by modulating proliferative parameters, while normal clones only gained 2 rounds of cell division in the same scenarios. This difference results from cell crowding creating a disproportional disadvantage to non-lesioned clones: In the scenario with increased proliferation, the first clones to colonize the organ can massively proliferate and completely surround the TECs, preventing later arriving clones from accessing TEC-derived pro-proliferative signals (Figure 4—video 2, right column). Thus, these few ‘lucky’ early colonizers prevent late cells from attaining their full proliferative potential. Lesioned clones carry an intrinsic pro-proliferative advantage due to their lesions in the IL7R, enabling proliferation despite a lack of direct contact with TECs. In contrast, modifying parameters that affected IL-7 activity only enhanced the proliferation rate of non-lesioned clones but showed no impact on either IL7R^HI or IL7R^DA clones (Figure 4—figure supplement 3). While this outcome is expected for IL7R^DA clones, which are insensitive to IL-7, for IL7R^HI clones, it indicates that 10-fold receptor overexpression is enough to activate downstream pro-proliferative effects of the pathway regardless of the parameter permutations that we tested. In contrast, proliferation of normal clones is strongly affected in these conditions, gaining up to 3 rounds of cell division between least and most proliferative conditions.

In our systematic approach, we did not detect an added effect on clonal expansion from combining reduced cell motility or autocrine IL-7 production with IL7R^HI or IL7R^DA lesions (Figure 4—figure supplement 4). Added IL-7 secretion by lesioned clones with the autocrine lesion had only a small effect on proliferation of non-lesioned clones in the same simulation. The most critical lesion in conferring a proliferative advantage to lesioned clones was a delay in differentiation. This delay essentially doubled the duration of the proliferative phase, and indeed we observed up to 8 rounds of cell division, double the level of clones having only a single lesion in the IL7R. This result fits well with in vivo data, suggesting that perturbations in differentiation have been implicated in contributing to the initiation and progression of cancer (Ruijtenberg and van den Heuvel, 2016).

Surprisingly, the expansion of lesioned clones was amplified when the TEC niche was at its sparsest – a result opposite to the non-lesioned clones which instead profit from a dense TEC niche (Figure 4—figure supplement 5). While non-lesioned clones gained 2 rounds of division with an increase in TEC density, lesioned clones lost 2 rounds of division. Indeed, parameter combinations with reduced density and size of TECs were overrepresented in the scenarios that most increased proliferation (Figure 4—figure supplement 6A–C). For example, in the most extreme scenario observed in our systematic approach, IL7R^HI or IL7R^DA clones underwent up to 18 rounds of cell division – far above the 12 rounds predicted from modifying each group of parameters in isolation (+6 from proliferation, +4 from added lesions, +2 from reduced TEC density, assuming an additive effect; Figure 4—figure supplements 2–5). This result suggests that a defective niche amplifies the effect of other pro-proliferative modulations such as slower differentiation and shorter cell cycle, conferring a synergistic advantage to IL7R^HI and IL7R^DA clones compared to their normal counterparts. Note that at 18 rounds of division, lesioned clones were so massive that the total lesioned cell volume was over a 1000-fold larger than the volume of the simulated organ slice (Figure 4—figure supplement 6D). In comparison, in the homeostatic condition, the total volume occupied by both thymocytes and TECs amounted to only 0.7-fold of the available tissue volume. Though our simulations by default include physical volume exclusion effects that prevent dense packing, we did not implement density-dependent feedback on proliferation, enabling cells to duplicate their volume at cell division and proliferate uncontrolled despite the lack of available space. Presumably, similar conditions in vivo would lead to thymus hyperplasia. In support of this view, our experimental manipulations of IL-7 levels in the thymus had a direct effect on organ size (Figure 2I, Figure 2—figure supplement 2C and D).

Identifying a sparse TEC network as a new factor that could influence the clonal expansion of lesioned clones was an unexpected outcome. This is because, in the WT situation, a sparser TEC network results in decreased extracellular IL-7, leading to reduced clonal size (Figure 2). While being disadvantageous for non-lesioned clones, this compromised thymic niche proved to be optimal for a massive expansion of lesioned clones. Conversely, we observed that a denser TEC network – characterized by a higher density of TECs with more protrusions and larger size – positively influenced the population size of non-lesioned cells but had a negative impact on the IL7R^HI and IL7R^DA lesioned clones (Figure 4—figure supplement 5). Similarly, in conditions of high cell crowding, lack of contact with TECs put non-lesioned clones at a disadvantage, while lesioned clones appeared to benefit (Figure 4—figure supplement 2, Figure 4—video 2). This outcome might be explained by the signaling dynamics within our model implementation, where the engagement of the NOTCH1 receptor on thymocytes with the DLL4 ligand on TECs, on the one hand, promotes proliferation and, on the other hand, is a prerequisite for the differentiation process until commitment to a nondividing fate. In the model, if the NOTCH1 receptor fails to engage with DLL4 on TECs due to their low density or cell crowding effects, cells will remain in an undifferentiated state for a longer time. A closer inspection of time-normalized clones in scenario 6 (very sparse TECs) and in scenario 26 (very dense TECs) indeed confirms that the thymocyte population has a longer turnover time when TECs are sparse (almost 80 hr, Figure 2D, middle panel) and a shorter turnover when TECs are dense (less than 60 hr, Figure 2D, right panel), indicative of the effect of Notch signaling. Together, cells in an environment with fewer TECs and thus lower Notch signaling will experience a prolonged proliferative phase but still require pro-proliferative signals to commit to the cell cycle. Owing to the higher levels of IL7R signal intrinsic to the lesion, lesioned clones will be more easily competent to proliferate even without Notch signaling and, therefore, benefit from a lower TEC density.

In our virtual thymus model, ETPs simultaneously assess whether the combined sum of IL7R and NOTCH1 signals surpasses a threshold required for entry into the cell cycle (Aghaallaei et al., 2021; Figure 5A). Besides its effect on proliferation, NOTCH1 signaling also promotes ETP differentiation (Aghaallaei et al., 2021; Aghaallaei et al., 2022), and we used this fact and the observations in notch1b medaka mutant phenotypes to include an accelerating effect of NOTCH1 on differentiation in our model (Aghaallaei et al., 2021; see also Appendix 1). Thus, NOTCH1 has a dual effect in thymocytes, which we implemented in our model as an incoherent feed-forward loop (Figure 5A): NOTCH1 signaling promotes proliferation and accelerates differentiation, while differentiation inhibits proliferation by preventing further entry into the cell cycle. IL7R signaling promotes proliferation independently of NOTCH1 and has no direct or indirect effect on the differentiation process. To further explore the interplay between these two factors, we compared scenarios where lesioned clones exhibited only constitutive activation of the NOTCH1 receptor (hereafter called NOTCH1^DA) or in combination with either IL7R^DA or IL7R^HI. The simulations predicted that lesioned clones with only NOTCH1^DA modification displayed a very small advantage compared to non-lesioned clones; this increase in proliferation was only marginally higher than expected for a clone expressing the highest endogenous levels of IL7R (hereafter IL7R^WT; Figure 5B). In contrast, lesions that specifically delay differentiation led to higher proliferation. Indeed, the addition of NOTCH1^DA modification to IL7R^DA or IL7R^HI lesioned clones did not show any effect on clonal expansion, but combining delayed differentiation with IL7R^DA or IL7R^HI doubled proliferation. This result indicates that modulating the effect of NOTCH1 globally does not affect proliferation in our model, which is likely due to feedback inhibition, but specifically targeting the differentiation process leads to clonal expansion.

Figure 5 with 1 supplement see all

Download asset Open asset

Figure 5—source data 1

Detailed statistics for phenotype comparison using Fisher’s exact test and pairwise comparisons; adjusted p-values were calculated with the Benjamini and Hochberg method.

Values reported in figures are the adjusted p-values. Significance threshold was set to 0.05. Data rounded to 4 decimal places.

https://cdn.elifesciences.org/articles/101137/elife-101137-fig5-data1-v1.docx

Download elife-101137-fig5-data1-v1.docx

To test the outcome of the virtual model, we performed a series of in vivo experiments to induce NOTCH1 and IL7R signaling in thymocytes. To constitutively activate NOTCH1 signaling in thymocytes, we took advantage of the previously cloned medaka notch1b intracellular domain (NICD) construct (Aghaallaei et al., 2021; Aghaallaei et al., 2022), which was driven by a thymocyte-specific promoter (Bajoghli et al., 2015). Furthermore, the destabilizing PEST domain from the NICD was removed in this construct (hereafter referred to as NICD^∆PEST; Figure 5C). This genetic modification was made due to the known impact of the PEST domain on protein stability, as nonsense mutations lacking the PEST domain of the human NOTCH1 gene have been frequently observed in T-ALL patients (Breit et al., 2006; Ferrando, 2009; Weng et al., 2004; Figure 5—figure supplement 1A). To mimic the in silico IL7R^HI clones, we employed a previously developed construct wherein a thymocyte-specific promoter drives the medaka full-length il7r cDNA (Aghaallaei et al., 2021). We then performed a mutation in the extracellular juxtamembrane-transmembrane region of the medaka il7r cDNA (Figure 5—figure supplement 1D) to develop a dominant active IL7R form, akin to the NPC mutation found in the IL7R gene in some T-ALL patients (Zenatti et al., 2011; Oliveira et al., 2022), hereafter called il7r^NPC. In our experimental setup, the promoter also co-expressed GFP, which allowed us to (i) identify thymocytes expressing the oncogenes, (2) determine the clonal expansion of cells expressing the oncogene, and (3) assess thymus hyperplasia and infiltration into other organs, a characteristic feature for T-ALL. DNA constructs were then injected into blastomeres of embryos at one-cell stage, and they were observed during their development using live imaging, with a focus on two specific phenotypes: First, we assessed whether the thymus exhibited enlargement beyond the normal size for its developmental stage. In particular, thymus lymphoma was defined as hyperplasia with an increase of more than twice the organ’s typical size. Second, we identified T-ALL by observing infiltration of GFP-co-expressing malignant cells in other organs, including the brain, intestine, or heart (Figure 5D). To align our in vivo experiments with the in silico conditions, we limited our observations to 11 dpf, concluding the experiments at the freshly hatched yolk-sac larval stage. Since thymopoiesis begins at 3 dpf in medaka embryos (Bajoghli et al., 2015; Bajoghli et al., 2009), this means we monitored thymus growth for a period of 8 days.

None of the embryos injected with the NICD^∆PEST construct (N=32) displayed thymus hyperplasia at 11 dpf (Figure 5E), supporting our simulation outcome showing that global activation of NOTCH1 signaling in thymocytes alone is insufficient to result in clonal expansion within a short time frame of 8 days. The MYC oncogene is an endogenous downstream target of NOTCH1 and plays a major role in NOTCH1-induced transformation (Sanchez-Martin and Ferrando, 2017). Therefore, to attempt to shift the balance of NOTCH1 action toward proliferation, we decided to introduce the medaka mycn cDNA into our construct. Consequently, we observed that 25% of injected embryos with NICD^∆PEST and mycn (N=36) displayed thymus hyperplasia at 11 dpf (Figure 5E). Among them, 8% also exhibited a massive infiltration of GFP-expressing cells in other organs such as the brain and gut. Further analysis of sorted GFP-expressing thymocytes revealed ectopic il7 expression, while the expression level of endogenous il7r remained unchanged (Figure 5—figure supplement 1B). Whole-mount in situ hybridization (WISH) analysis further confirmed a robust upregulation of il7 expression in the thymus of these embryos (Figure 5—figure supplement 1C). MYCN is a basic helix-loop-helix transcription factor that is downstream of several pro-proliferative signaling pathways (Ruiz-Pérez et al., 2017), including NOTCH1 in the context of thymocytes (Sanchez-Martin and Ferrando, 2017). Together with our experimental results, it is therefore likely that MYCN could enact the proliferative effect of NOTCH1, shifting the balance in the incoherent feed-forward loop toward proliferation. Interestingly, our results also indicate that transformed thymocytes acquire the ability to release IL-7 (Figure 5—figure supplement 1B and C), which may further promote proliferation in an autocrine fashion.

To test whether constitutive IL7R activation could enhance the development of thymus hyperplasia and T-ALL, we injected various constructs to overexpress either il7r or the dominant active il7r^NPC alone, or in combination with NICD^∆PEST and mycn. Thymus hyperplasia was found in a slightly higher frequency of 26% (N=49) or 59% (N=29) of embryos, respectively, when only il7r or il7r^NPC was overexpressed in thymocytes. Additionally, we FACS-sorted thymocytes and performed a qPCR for lck, since this gene was shown to be upregulated in zebrafish after IL7R^DA T-ALL development using RNASeq data (Oliveira et al., 2022). Indeed, we detected a significant upregulation of lck upon il7r^NPC overexpression when compared to the WT counterparts. However, upregulation was not significantly different when the WT was compared to the il7r^HI condition (Figure 5—figure supplement 1D). Combining il7r^NPC and NICD^∆PEST and mycn yielded a notably elevated frequency of ~73% (N=33) for thymus hyperplasia (Figure 5E). Of this group, 33% also exhibited a T-ALL phenotype, suggesting an additive effect of these factors in the T-ALL development (Figure 5—source data 1). GFP and pH3 double staining of embryos exhibiting the T-ALL phenotype further confirmed that the combination of NICD^∆PEST, mycn, and il7r overexpression in thymocytes resulted in a higher number of mitotically active cells within the thymus (Figure 5F). Notably, we observed many pH3-stained cells in the inner zone of the medaka thymus (Figure 5G, arrows), an area where WT thymocytes are mitotically quiescent (Bajoghli et al., 2015; Aghaallaei et al., 2021).

The prediction that higher NOTCH1 signaling does not lead to increased proliferation is likely due to the incoherent feed-forward loop downstream of NOTCH1, which shuts down excessive proliferation by triggering differentiation to a nondividing cell fate (Figure 5B). We observed a similar effect in our in vivo experiments, where activation of NOTCH1 alone did not lead to thymus hyperplasia within the observed time window of 8 days (Figure 5E). In contrast, combining NOTCH1 with its downstream pro-proliferative effector mycn led to thymus hyperplasia. Intriguingly, in other in vivo systems, long-term activation of NOTCH1 signaling in thymocytes is one of the main drivers of T-cell leukemogenesis and T-ALL development (Weng et al., 2004; Lin et al., 2006; Chen et al., 2007; Liu et al., 2017; Neumann et al., 2015). Several in vivo studies have demonstrated that constitutive activation of the NOTCH1 receptor leads to leukemia development after several weeks in zebrafish (Chen et al., 2007; Blackburn et al., 2012) or months in mice (Chiang et al., 2008; Hu et al., 2009; Sharma et al., 2006; Wendorff and Ferrando, 2020). We did not observe thymus hyperplasia when upregulating NOTCH1 signaling alone, which may stem from the shorter time window of our in vivo experiments (i.e. 8 days). However, in the virtual thymus model, a sole upregulation of NOTCH1 cannot produce thymus hyperplasia regardless of the duration of pathway activation. Different possibilities could reconcile these discrepancies. First, the balance between pro-proliferative and pro-differentiation effects of NOTCH1 in vivo could be different from what we modeled. Second, there could be additional negative feedbacks between proliferation and differentiation, which our model does not consider. Third, proliferation in vivo could continue even past differentiation if NOTCH1 signaling is sufficiently stimulated. Finally, long-term pathway misregulation could act as a stressor that promotes additional compounding lesions.

Overall, our in vivo results support the outcomes of our simulations, showing that overexpression of IL7R alone, but not NOTCH1 alone, is sufficient to induce thymus hyperplasia in a very short time period. However, dual activation of IL7R and the NOTCH1-MYCN axis can promote clonal expansion in vivo, leading to the rapid development of thymus hyperplasia and T-ALL in the medaka model system within 8 days.

One of the most unexpected outcomes of our simulations was that a defective TEC network provides lesioned clones with a substantial advantage in proliferation. Given the crucial role of IL-7 in proliferation of thymocytes, and considering that TECs, and not thymocytes, serve as the primary source of IL-7 cytokine in the normal thymic milieu, we next asked whether elevated il7 expression in TECs alone is sufficient to induce thymus hyperplasia and to enhance T-ALL development in a thymus where thymocytes express high levels of il7r. This question was also driven by the results of scenario 27, which showed a twofold increase in thymocyte population size when all TECs express IL-7 (Figure 2E–G). In this thymic niche enriched with IL-7, lesioned clones either with IL-7^HI or extended differentiation times showed a significant advantage in clonal expansion compared to their WT counterparts (Figure 6A). Further, the impact of alterations in the thymic niche on T-ALL development has not been previously addressed. Therefore, we monitored the thymus development of the TEC:IL-7^HI embryos until 11 dpf to mimic scenario 27 in vivo. We found that 12% of these embryos (N=24) displayed signs of thymus hyperplasia (Figure 6B); however, we did not observe massive infiltration of GFP-expressing cells into other organs. The frequency of thymus hyperplasia was increased to 41% (N=32) when DNA constructs designed to overexpress il7 in TECs were co-injected with the construct overexpressing il7r in thymocytes (Figure 6B). Comparably, thymus hyperplasia was only observed in 26% of embryos injected with the IL7R overexpression construct in thymocytes alone (N=49). A synergistic effect was observed when il7 was overexpressed in TEC along with NICD and mycn in thymocytes (Figure 6B). Interestingly, the frequency of embryos with T-ALL development was increased to 25% (N=48) in this group. Taken together, our in vivo results reveal that an excess of IL-7 in the thymic environment, combined with alterations that affect the differentiation status of thymocytes, triggers the T-ALL development in a short period of time.

Figure 6

Download asset Open asset

Overall, several of our experimental observations indicate that there might be additional negative feedback from pro-proliferative signaling to differentiation (Figure 6C). Most strikingly, manipulations in IL7R or in thymic IL-7 levels are sufficient to induce hyperplasia and T-ALL (Figures 5E and 6B) – an outcome that is much more severe than the model predicts (Figures 5B and 6A). In simulations, we could only obtain such massive overproliferation by modifying multiple parameters, including delaying differentiation, accelerating cell divisions, and reducing TEC density (Figure 4—figure supplement 5), as differentiation in our model is a hard stop to proliferation regardless of pro-proliferative stimuli. Similarly, the experimental observation that NOTCH1-MYCN upregulation leads to thymus hyperplasia and T-ALL development also suggests that proliferation can proceed unbridled by differentiation in vivo when sufficiently stimulated. A negative feedback loop from pathways regulating proliferation (or cell cycle entry) on differentiation pathways could explain these discrepancies between the virtual thymus model and in vivo observations. We therefore propose that continued thymocyte proliferation can downregulate differentiation pathways and maintain cells in a constant proliferative state.

This study demonstrates the powerful synergy between computational modeling and in vivo experimental validation in dissecting the complex cellular and molecular interactions within an organ that drive disease initiation and progression. Further, we highlight a previously underexplored area – the role of TECs and the thymic niche in the development of T-ALL – thereby bridging a major gap in our understanding of the earliest stages of disease development.

By simulating over 1500 scenarios, we identified potential drivers such as alterations in the shape, density, and spatial organization of the TEC network – factors that are inherently difficult to manipulate in experimental systems. Consistent with these computational predictions, our in vivo results reveal that an enriched thymic IL-7 milieu can accelerate disease progression, particularly when thymocytes display enhanced IL7R levels or express constitutively active NOTCH1 and MYCN overexpression. These results challenge the traditional gene-centric paradigm of T-ALL and emphasize the critical influence of the spatial structure of the thymic microenvironment.

However, it is worth noting several limitations of our current virtual thymus model. The cell-based computational framework was calibrated to mimic the medaka embryonic thymus (Bajoghli et al., 2009; Bajoghli et al., 2015; Aghaallaei et al., 2021), and as such, it may not fully capture the cellular dynamics or thymic architecture of other species, particularly human. On the technical side, the simulation lacks the ability to represent changes in TEC shape over time and TEC proliferation, as all particles are treated as loose objects and hence movement of a TEC would lead to detachment of its protrusion particles. This limitation could be overcome with a more advanced biomechanical framework that includes bonded interactions between particles such as Tissue Forge (Sego et al., 2023), which would ensure that particles belonging to the same cell stay together as one part of the cell moves. Due to this limitation, the dynamic remodeling of the thymic niche is not fully accounted for in our current model. Furthermore, the model simplifies certain T-cell developmental stages as it does not consider cross-regulatory feedback loops between IL-7 and Notch pathways and does not include thymic growth over time. These simplifications are due to the lack of species-specific knowledge and the tools to investigate these questions in medaka. Nonetheless, we anticipate that adapting the framework for use in murine or human thymic contexts, where stage-specific markers and niche components are well characterized and where detailed morphometric studies of TEC architecture are emerging (Lagou et al., 2024) would enable us to explore these questions.

Despite its current simplifications, our integrative approach offers valuable insights into how microenvironmental factors influence leukemia expansion and provides a framework for exploring their role in chemoresistance and relapse. Ultimately, these findings open new avenues for developing therapeutic strategies aimed at disrupting the leukemic-supportive niche, with the potential to improve treatment outcomes in T-ALL.

To develop, implement, and simulate the virtual thymus model, we used the modeling and simulation software EPISIM (Sütterlin et al., 2013; Sütterlin et al., 2017). This multiscale simulation software uses an agent-based paradigm to represent cells as spheres or ellipsoids in three-dimensional space, enabling each of the cells to individually perform internal processes based on flow diagrams and logic rules (e.g. if/else statements) or differential equations, and also implements a partial differential equation solver to simulate diffusion of chemicals in the extracellular space. The virtual thymus model has been comprehensively described in the supplementary material to our previous study (Aghaallaei et al., 2021). In this work, we used this model as a baseline to introduce a small number of additions, as explained in the following. A summary of the full model implementation and parameter values is provided in Appendix 1.

Clone size was defined as the number of terminal leaves in a cell lineage tree (see scheme in Figure 4C); terminal leaves were counted regardless of the ultimate fate of those cells (i.e. positive or negative selection). All clones of the same genotype (WT or lesioned) belonging to the same simulated scenario were averaged across all replicate simulations of that scenario. This calculation gives the average number of cells per clone in a given scenario. The number of cells was log2-transformed to represent rounds of cell division.

Medaka (O. latipes) husbandry was performed in accordance with the German animal welfare standards (Tierschutzgesetz §11, Abs. 1, Nr. 1, husbandry permit no. 35/9185.46/Uni TÜ). The transgenic line (tg) lck:gfp was described previously (Bajoghli et al., 2015). All experiments conducted in medaka embryos were performed prior to the legal onset of animal life stages under protection, utilizing both males and females of the WT line Cab (Loosli et al., 2000).

To overexpress il7 in TECs, we used the construct ccl25a:il7, which was described previously (Aghaallaei et al., 2021). To overexpress il7r in thymocytes, the full-length medaka il7r cDNA (Aghaallaei et al., 2021) was cloned into vectors containing a medaka thymocyte-specific promoter (Bajoghli et al., 2015) that drives a fluorescent protein (sfGFP or TagRFP). To generate a dominant active form of il7r, we introduced three amino acids, asparagine (N), proline (P), and cysteine (C) (known as NPC mutation, as shown in Zenatti et al., 2011, and Oliveira et al., 2022), in the medaka il7r extracellular juxtamembrane-transmembrane interface region after position 266 using site-directed mutagenesis (Figure 5—figure supplement 1D), this corresponds to the human NPC mutation in the position 242 of IL7R gene. To activate Notch signaling in thymocytes, we utilized the medaka NICD, as described previously (Aghaallaei et al., 2021), and removed its PEST domain to enhance protein stability, mimicking the situation reported in most T-ALL patients harboring gain-of-function mutations in the NOTCH1 gene (Weng et al., 2004). For overexpressing the mycn oncogene, we first isolated the full-length medaka mycn cDNA (accession number: ENSORLG00000022362) and introduced a mutation at position 44 from proline to leucine (P44L) using site-directed mutagenesis, as identified in T-ALL patients (Liu et al., 2017). The DNA fragments were then cloned into vectors containing a medaka thymocyte-specific promoter, either alone or in combination (Bajoghli et al., 2015).

Plasmids at concentration 10–25 ng/µl together with I-SceI meganuclease and NEB buffer (New England BioLabs) were co-injected into the blastomere at one-cell stage embryos. Fluorescent signals based on the constructs sfGFP, tagRFP, or mTurquoise were used to select positive embryos.

The CRISPR-Cas9 approach was employed to generate medaka il7 crispant. CRISPR RNA (crRNA AGTAGACTGATGCAAAGAAG) was designed using the CCTop website (Stemmer et al., 2015; Stemmer et al., 2017) and ordered from IDT. RNP complexes were prepared as described before (Hoshijima et al., 2019) using Alt-R tracrRNA (IDT) and Alt-R S.p. Cas9 nuclease (IDT). The injection mixture contained a total of 25 µM crRNA:tracrRNA duplex and 25 µM Cas9. Injection was performed into the blastomere at the one-cell stage transgenic embryos carrying the lck:gfp reporter. Injected embryos were raised until 8 dpf and then fixed in 4% PFA/2x PBS+0.1% Tween20 for whole-mount immunostaining. To correlate the phenotype with the genotype, each embryo was genotyped using PCR and Sanger sequencing.

The development of injected embryos was monitored using a NIKON SMZ18 stereo-fluorescent microscope. Only embryos with fluorescent signals in the thymus region were selected for further analysis. Thymuses of embryos before hatching or freshly hatched yolk sac larvae were examined. Thymuses with at least a twofold increase in size compared to transgenic embryos carrying the lck:gfp reporter construct (Bajoghli et al., 2015) were considered as having thymic hyperplasia. The identification of cells expressing fluorescent proteins and oncogenes outside the thymus, such as in the brain, gut, and heart, was considered indicative of the T-ALL phenotype.

Immunostaining procedures followed established protocols (Inoue and Wittbrodt, 2011; Aghaallaei et al., 2021). Briefly, mitotically active cells were identified using a rabbit anti-phosphohistone-3 antibody (Ser10, Millipore 06-570, 1:500 dilution), with a Cy3-donkey anti-rabbit immunoglobulin G secondary antibody (the Jackson Laboratory, 711-165-152; 1:500 dilution). GFP expression was detected using a goat Anti-GFP antibody (Abcam, ab5450, 1:500 dilution), with Alexa 488-donkey anti-goat IgG (Abcam, ab15129, 1:600 dilution) as the secondary antibody. To quantify the number of pH3-expressing cells, the entire thymus region was imaged using an LSM710 (Zeiss) confocal microscope with z-stacks (z=1 µm).

WISH in medaka embryos was carried out using digoxigenin-labeled RNA probes, as previously described (Aghaallaei et al., 2005). Probes used in this study for il7r and il7 were previously described (Aghaallaei et al., 2021).

Injected embryos were first smashed through a 40 µM strainer (Greiner) and cells were collected in 0.9x PBS with 500 UI Heparinum Natricum (Liquemin). Sorting was then performed using the Cell Sorter MA900 (Sony Biotechnology). WT embryos were used as negative sorting control. Only GFP positive cells were collected and used for RNA isolation. RNA of sorted cells was isolated by using NucleoSpin RNA XS (Macherey-Nagel) and 20 ng carrier RNA, following the manufacturer’s protocol. RNA was treated with rDNase and eluted in 14 µl dAqua. The first-strand cDNA synthesis was carried out with random hexamer primers and SuperScript III Reverse Transcriptase (Thermo Fisher Scientific) by following the manufacturer’s protocol. SYBR Green Kit (Applied Biosystems) was used for Quantitative PCR on the LightCycler 480 (Roche). The data was evaluated in Microsoft Excel using the ΔCt method and normalized to the housekeeping gene ef1a. The primers used for RT-PCR have been described previously (Aghaallaei et al., 2021).

Statistical analysis was performed using RStudio version 2024.4.2.764 (Posit team, 2024) using R version 4.4.0 (R Development Core Team, 2024), the R library jtools (Long, 2022) and GraphPad Prism version 8.0.2. The two-tailed Fisher’s exact t-test was used to compare the thymus hyperplasia phenotype with normal thymus. Welch’s t-test was used for thymus volume size comparison in WT and il7 crispants, and a multiple t-test was used for comparison of normal thymus with thymus hyperplasia and T-ALL phenotype. Mann-Whitney test was used to evaluate the qPCRs. p-Values<0.05 were considered as statistically significant. The number of biological samples for each experiment (N) is indicated in the figures, and bar graphs present the absolute numbers with mean ± SD, ± SEM, or the percentage as indicated in the figure legends.

The cell-based model of the virtual thymus was implemented in the platform EPISIM, which simulates center-based biomechanical interactions for spherical or ellipsoid cells (Sütterlin et al., 2013, Sütterlin et al., 2017).

A detailed description of the mathematical underpinning and biological inspiration of the model and its parametrization was published previously (Aghaallaei et al., 2021). Parameters were adapted from this previous work; a few additional scaling parameters were added to enable differential behavior in lesioned thymocytes (Appendix 1—Tables 1 and 2). The code and model files of the most recent version associated with this publication are freely available online on a Zenodo repository (Tsingos, 2024). In the following, we recapitulate the main features of the model implementation for the reader’s convenience.

The virtual thymus model represents processes in the thymus at the tissue, cellular, and subcellular scales (Appendix 1—figure 1). At the tissue scale, we model a thin three-dimensional (3D) slice of the thymus of medaka hatchlings.

At the cellular scale, we model TECs and thymocytes, and both cell types are represented as one or more spherical particles. Biomechanical interactions modeled at the cellular scale are cell motility and cell-cell contact mechanics. This scale also includes processes that change the number of cells such as cell influx into and efflux out of the tissue volume, cell proliferation, and cell death.

At the subcellular scale, we model biochemical interactions such as extracellular diffusion of cytokines, signaling between cells, intracellular signal transduction, and include phenomenological models of the thymocyte cell differentiation and proliferation programs. In the following sections, we provide an overview of each of these processes and refer to original publications for further details where appropriate.

The details of the mathematical implementation of contact forces have been published elsewhere (Sütterlin et al., 2017). In essence, contact forces are made up of repulsive pressure forces that push particles apart and attractive adhesion forces that pull neighboring particles toward each other. Two particles i and j are considered neighbors if

where d_ij is the Euclidean distance of the particle centroids, d_opt(i,j) is an optimal distance obtained from the respective dimensions of particle i and its neighbor particle j (for the spherical particle interaction radii that we use in this work the calculation of d_opt(i,j) simplifies to the sum of radii of i and j), and δ_adh > 1 is a parameter that defines the neighborhood interaction distance.

The contact forces for each particle-neighbor pair are given by a piecewise function with multiple distance thresholds based on d_ij and d_opt(i,j), which results in concentric interaction zones (illustrated in Appendix 1—figure 2):

Here, δ_{ol_max}, δ_ol, d_{ol_min} are parameters that scale the interaction distances. The parameter δ_{ol_max} ∈ (0, 1) defines the inner repulsive interaction zone where the repulsion force F_pr,exp exponentially increases as d_ij decreases; δ_ol ∈ (0, 1) defines the outer repulsive interaction zone where the repulsion force F_pr,lin linearly increases as d_ij decreases; d_{ol_min} ∈ (0, 1) defines a neutral zone with no net force; δ_adh, which is used to determine neighborhood in Equation A3 additionally defines the interaction zone where adhesive forces F_adh dominate (Appendix 1—figure 2). For the functional form chosen for the force terms, we refer to Sütterlin et al., 2017. Based on previous work, the values δ_{ol_max} = 0.5; δ_ol = 0.85; d_{ol_min} = 0.1 μm; δ_adh = 1.3 were chosen as they result in appropriate spacing of particles without mechanical instability (Sütterlin et al., 2017, Tsingos et al., 2019; Aghaallaei et al., 2021). Note that our approach is consistent with state-of-the-art models of force-based interactions between particles to simulate cells (Pathmanathan et al., 2009; Osborne et al., 2017).

During the simulation, each cell determines an active instantaneous displacement d_act from its instantaneous speed s and a unit vector v representing cell orientation

where t is time. Both s and v are functions of the cell’s internal state and its position in the tissue and are dynamically evaluated during the simulation. For an in-depth explanation, see the supplementary material of Aghaallaei et al., 2021. In brief, the functional form for the speed s was chosen as

where μ_S and σ_S are model parameters, X ∈ (–1, 1] is a random uniform number, and τ is either set to τ=1 for thymocytes entering or exiting the thymus, or set to τ = τ_diff for thymocytes currently residing in the thymus (for an explanation on τ_diff, see section IX. Subcellular scale: differentiation model). Essentially, setting τ = τ_diff enables speed to be linearly increasing with developmental stage, which is an assumption we made based on experimental measurements of thymocyte motility in vivo (Bajoghli et al., 2015). For thymocytes of lesioned clones with reduced cell speed, we set τ=0.5. Coefficients μ_S and σ_S were fit to the speed measured in confocal imaging data of thymocytes in medaka; different coefficients for μ_S were used for thymic immigrant/emigrants and thymic resident cells (Aghaallaei et al., 2021). Cells from a lesioned clone with slower cell speed had the value of their μ_S coefficient halved compared to WT cells.

The orientation vector v is a unit vector obtained by combining a directional vector u pointing from the current cell position toward a target location and a random vector w that introduces a random bias. The directional u components and random w components are weighted differently depending on cell location, developmental stage, and each cell’s intrinsic expression levels of the cytokine receptor Ccr9b (schematically summarized in Appendix 1—figure 3). For a detailed explanation of the implementation, choice of parameter values, and a sensitivity analysis of tissue homeostasis with respect to parameters, we refer to our previous publication (Aghaallaei et al., 2021).

The active displacement d_act is used to obtain the cell motility force F_act

where γ is the friction or damping constant of the environment. In the implementation, the active force F_act obtained in the previous simulation step with (Equation A7) is used in the equation of motion (Equation A2) at the beginning of each new simulation step. Because contact forces F_contact to neighboring cells may hinder cell motility, the effective displacement of cells that results from Equation A2 may differ from the displacement calculated in Equation A5.

The model includes two mechanisms for adding new cells to the simulated volume: (1) influx of early undifferentiated thymocytes from areas elsewhere in the body and (2) increase in cell number due to thymocyte proliferation.

In medaka, thymocytes enter the organ from the ventral side at an estimated influx rate of 9 thymocytes per hour (Bajoghli et al., 2015). The model implements cell influx by instantiating new cells at two positions at the bottom of the simulation box with rate k_spawn. Since the model represents about one-tenth of the full organ’s volume, we scaled down the influx rate accordingly to k_spawn = 0.45 thymocytes per hour, which adds up to 0.9 thymocytes per hour.

Cell proliferation follows a model of the cell division cycle, where progression through G1 phase and commitment to cell division depends on intracellular signaling at the subcellular level (see section X. Subcellular scale: cell proliferation model for an in-depth explanation). Once the subcellular program of the cell division cycle is complete, the implementation instantiates a new cell adjacent to the cell that completed the cycle at a random orientation; this simulates the close positioning of daughter cells just after cytokinesis with a random spindle orientation.

Two processes account for cell removal from the simulation: (1) efflux of fully differentiated thymocytes that are positively selected, (2) death of fully differentiated thymocytes that are negatively selected. To our knowledge, there is no evidence of substantial cell death at earlier stages of differentiation in medaka (Bajoghli et al., 2015); hence, we do not include other mechanisms of cell removal.

Once fully differentiated, each thymocyte undergoes thymic selection (see subsection IX. Subcellular scale: differentiation model for details). Positively selected thymocytes migrate out of the thymus and eventually leave the simulation box, while negatively selected thymocytes stay in position and slowly shrink until their radius reaches 25% of the original value, at which point they are removed from the simulation.

X. Subcellular scale: cell proliferation model

The cell proliferation model is summarized in Appendix 1—figure 6. Thymocytes in the model are only competent for proliferation if they are undifferentiated (τ_diff < 1). Each cell draws a random Erlang-distributed random number E with mean μ=7 hr and shape value k=50. This random number is then set as the minimum duration of the next cell cycle – i.e., the fastest possible cell cycle if sufficient pro-proliferative signals are present. We subdivide the randomly chosen cell cycle interval E into four phases using the assumption that a cell spends κ_S = 50% of the cell cycle duration in the S phase and t_M = 30 min in the M phase, as follows:

(A17) ${\displaystyle {\displaystyle \begin{array}{ll}{t}_{\mathrm{G}1}& =\frac{E}{2}\left(1-\frac{{\kappa }_S}{2}\right)\\ t_{\mathrm{S}}& =E{\kappa }_S\\ t_{\mathrm{G}2}& =\frac{E}{2}\left(1-\frac{{\kappa }_S}{2}\right)-t_{\mathrm{M}}\\ t_{\mathrm{M}}& =t_{\mathrm{M}}\end{array}}}$

where t_G1 is the time spent in G1 phase, t_S the time spent in S phase, t_G2 the time spent in G2 phase. The model also includes a G0-like quiescence due to the absence of sufficient pro-proliferative signals. At the beginning of a new cycle, each cell sets the value of its internal variable τ_cycle = 0, which tracks the progression through the cell cycle. This variable is incremented as

(A18) ${\displaystyle {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}{\tau }_{\text{cycle }}=\{\begin{array}{ll}0& \text{if }{\tau }_{\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}}\le t_{G1}\text{ and }{\sigma }_{\text{IL-7}}+{\sigma }_{\mathrm{N}\mathrm{o}\mathrm{t}\mathrm{c}\mathrm{h}}>{\theta }_{prol}\\ \frac{1}{E}& \text{otherwise }\end{array}}}$

where σ_IL-7 is the IL-7 pathway signaling given by Equation A9 (or Equation A11 for clones with dominant-active receptor lesion), σ_Notch is the Notch pathway signaling given by Equation A13, and θ_prol = 1.4 is a threshold parameter.

When τ_cycle ≥E, the cell divides into two daughter cells, and for both daughters, the value is reset to τ_cycle = 0. Essentially, Equation A18 models how the IL-7 and Notch pathways act as independent permissive pro-proliferative signals on G1 phase progression; if these signals are absent or too weak, then the cell cycle is delayed as τ_cycle does not increase and the cell remains in a G0-like quiescent state.

During the initial creation of the model, the sensitivity of the model with respect to the parameter values in Equation A17 and Equation A18 was studied by parameter scan (Aghaallaei et al., 2021). To choose appropriate values for the mean cell cycle duration μ and the threshold proliferation-permissive level θ_prol, the cell population size and the percentage of simulated mitotic events at homeostasis were compared to experimental measurements of pH3 staining from confocal slices (Aghaallaei et al., 2021).

Appendix 1—figure 1

Download asset Open asset

Appendix 1—figure 2

Download asset Open asset

Appendix 1—figure 3

Download asset Open asset

Appendix 1—figure 4

Download asset Open asset

Appendix 1—figure 5

Download asset Open asset

Appendix 1—figure 6

Download asset Open asset

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

The authors thank Anna-Sophia Hellmuth for assistance with il7 gRNA injection; the Institute of Medical Virology and Microbiology for continuous support in confocal microscopy; Larissa Doll, Julia Skokowa, and Karl Welte for support and encouragement. Initial computational work presented here was performed using the ALICE computer resources provided by Leiden University; later work was performed using the computer lab in the Theoretical Biology and Bioinformatics Department of Utrecht University. This work was supported by the Madeleine Schickedanz Kinderkrebsstiftung (grant number D.30.28666) and Deutsche Forschungsgemeinschaft (BA 5766/5-1). ET is funded by the Dutch Research Council (NWO) in the NWO Talent Programme with project number VI.Veni.222.323. We thank members of the Theoretical Biology and Bioinformatics division for their constructive feedback on the manuscript. We acknowledge support from the Open Access Publishing Funds of the University of Tübingen and the Faculty of Science of Utrecht University for supporting this research.

Medaka (Oryzias latipes) husbandry was performed in accordance with the German animal welfare standards (Tierschutzgesetz §11, Abs. 1, Nr. 1, husbandry permit no. 35/9185.46/Uni TÜ.). All experiments conducted on medaka embryos took place before the legal start of the stages of animal life that are protected by law.

1. Preprint posted: June 29, 2024
2. Sent for peer review: July 10, 2024
3. Reviewed Preprint version 1: October 9, 2024
4. Reviewed Preprint version 2: July 16, 2025
5. Version of Record published: September 16, 2025

You can cite all versions using the DOI https://doi.org/10.7554/eLife.101137. This DOI represents all versions, and will always resolve to the latest one.

© 2024, Tsingos, Dick et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.