Quantitative analysis reveals reciprocal regulations underlying recovery dynamics of thymocytes and thymic environment in mice

Thymic crosstalk, a set of reciprocal regulations between thymocytes and the thymic environment, is relevant for orchestrating appropriate thymocyte development as well as thymic recovery from various exogenous insults. In this work, interactions shaping thymic crosstalk and the resultant dynamics of thymocytes and thymic epithelial cells are inferred based on quantitative analysis and modeling of the recovery dynamics induced by irradiation. The analysis identifies regulatory interactions consistent with known molecular evidence and reveals their dynamic roles in the recovery process. Moreover, the analysis also predicts, and a subsequent experiment verifies, a previously unrecognized regulation of CD4+CD8+ double positive thymocytes which temporarily increases their proliferation rate upon the decrease in their population size. Our model establishes a pivotal step towards the dynamic understanding of thymic crosstalk as a regulatory network system.

Thymic crosstalk, a set of reciprocal regulations between thymocytes and the thymic environment, is relevant for orchestrating appropriate thymocyte development as well as thymic recovery from various exogenous insults. In this work, interactions shaping thymic crosstalk and the resultant dynamics of thymocytes and thymic epithelial cells are inferred based on quantitative analysis and modeling of the recovery dynamics induced by irradiation. The analysis identifies regulatory interactions consistent with known molecular evidence and reveals their dynamic roles in the recovery process. Moreover, the analysis also predicts, and a subsequent experiment verifies, a previously unrecognized regulation of CD4+CD8+ double positive thymocytes which temporarily increases their proliferation rate upon the decrease in their population size. Our model establishes a pivotal step towards the dynamic understanding of thymic crosstalk as a regulatory network system. https://doi.org/10.1038/s42003-019-0688-8 OPEN 1 Graduate School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-kuTokyo 113-8656, Japan. 2 RIKEN Center for Integrative Medical Sciences, 1-7-22 Suehiro-cho, Tsurumi-ku, Yokohama, Kanagawa 230-0045, Japan. 3 Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguroku, Tokyo 153-8505, Japan. 4 PREST, Japan Science and Technology Agency (JST), 4-1-8 Honcho Kawaguchi, Saitama 332-0012, Japan. 5 These authors contributed equally: Kazumasa B. Kaneko, Ryosuke Tateishi. *email: taishin.akiyama@riken.jp; tetsuya@mail.crmind.net T he thymus is an organ responsible for producing a large portion of T cells with appropriate repertoires 1 . However, it is relatively sensitive to insults from stress, viral infection, radiation, and other stimuli 2,3 . While a thymus in a healthy animal can be normally recovered from these damages, a relatively prolonged process of thymic recovery may impair T-cellmediated immunity due to a reduced replenishment of naïve Tcell repertoire during the recovery period 3,4 .
Sub-lethal dose radiation on mice has been utilized as an experimental model of the thymic regeneration after insults 5,6 . Ionizing irradiation is also broadly used for hematopoietic transplantation and cancer therapy 7,8 , and total body irradiation causes acute thymic injury and slow recovery of thymopoiesis. Several studies have shown that irradiation reduces cellularity, not only of thymocytes but also of thymic epithelial cells (TECs), which are major constituents of the thymic environment 5,9,10 . Because thymopoiesis is supported by interactions between thymocytes and TECs 11 , understanding thymic recovery requires characterization of the reciprocal regulations between thymocytes and TECs.
Concomitantly, various techniques to trace, perturb, and quantify cells involved in these events have enabled us to quantitatively characterize their dynamics [12][13][14][15] . By combining mathematical models with such quantitative data, dynamic aspects of thymopoiesis have been distilled into the form of detailed kinetic information, e.g., rates of proliferation, death, and differentiation 12,16 . Mehr et al. 17 developed the first kinetic model of thymocyte development using ordinary differential equations 18 . Since this seminal work, kinetic models of the thymopoiesis have been progressively refined by considering detailed cellularity and developmental states of the thymocytes, as well as by incorporating different experimental conditions [19][20][21][22][23] .
However, previous works have focused only on thymocytes. Thymic development and thymic recovery are not thymocyteautonomous but rather are supported by the thymic environment. In the last decade, we have accumulated molecular-biological evidence that the thymic environment itself is homeostatically maintained by thymic crosstalk, bidirectional interactions between the thymocytes and the thymic environment 11,24,25 . Among several cells comprising the thymic environment, cortical and medullary thymic epithelial cells (cTECs and mTECs) play integral roles in inducing and controlling proliferation, apoptosis, and lineage commitments of thymocytes 11,[26][27][28][29][30] . Thymocytes also regulate TECs by modulating their maturation and proliferations 5,[31][32][33] . Despite the evident relevance and importance of thymic crosstalk for the thymopoiesis and the thymic recovery, kinetic aspects of the reciprocal regulations between the thymocytes and the TECs have not yet been clarified.
In this work, we investigate the joint dynamics of thymocytes and TECs by combining a mathematical model with a quantitative measurement of the number of thymocytes and TECs during recovery after irradiation. Recovery dynamics are reproduced by our mathematical model, in which we identified reciprocal interactions between thymocytes and TECs that are relevant for recovery and consistent with thymic crosstalk. Furthermore, we demonstrate that the model provides an explanation for the mechanism of the dynamical change in population size. Particularly, our model predicts, and a subsequent experiment verifies, a previously unrecognized regulation of CD4+CD8+ doublepositive (DP) thymocytes, which temporarily increases their proliferation rate upon the decrease in their population size.

Result
Quantification of recovery dynamics of thymocytes and TECs. To quantitatively investigate the kinetic relationship between thymocytes and TECs as well as the establishment of thymic recovery, we artificially perturbed populations of thymocytes and TECs in thymi by using sub-lethal 4.5 Gy irradiation, and measured the dynamic changes in their population sizes over 3 weeks following irradiation ( Fig. 1a and Table 1). Figure 1b summarizes the changes in cell numbers, which were sorted based on conventional markers of thymocytes ( Fig. 1c and Supplementary  Fig. 1a) and TECs ( Fig. 1d and Supplementary Fig. 1b). Figure 1b shows that all types of investigated thymocytes and TECs decreased exponentially at different rates immediately after the irradiation. Then, both thymocytes and TECs started recovering within 10 days at the longest; the CD4−CD8− double-negative (DN) thymocytes and the cTECs began recovery within 5 days, whereas the CD4+CD8− single-positive (SP) thymocytes and the mTECs required longer intervals, reflecting the temporal order of the thymocyte development from DN to CD4+ SP (SP4) cells through interactions from cTECs to mTECs.
Upon recovery, the population sizes of all but the SP cells peaked around 15 days, and eventually returned to stationary numbers, which were almost equivalent to or at least half of the original population sizes before irradiation. Such overshooting behaviors suggest that the numbers of thymocytes and TECs are dynamically and mutually regulated via reciprocal interactions (Table 1).
Mathematical model can reproduce recovery dynamics. To infer regulatory interactions behind the dynamics, we constructed a mathematical model for the population dynamics of the thymocytes and the TECs using ordinary differential equations, which explicitly include five cell types: i 2 C ¼ fDN; DP; SP4; cTEC; mTECg. To account for the acute influence of irradiation on the cells, the total number of the cell type i at time t (day), n tot i ðtÞ is decomposed into two parts; n x i ðtÞ represents exponentially dying cells by the irradiation and n i (t) represents cells that survived or were newly generated after irradiation. n x i ðtÞ is assumed to decrease exponentially at a constant rate, ω i (day −1 ), as n x i t ð Þ ¼ n tot i 0 ð Þ 1 À p i ð Þe Àω i t , where p i is the proportion of survived cells after irradiation; we modeled the dynamics of n i (t) with ordinary differential equations. Therefore, the total number of cell type i, n tot i , which we observed experimentally, is described as n tot i t ð Þ ¼ n x i t ð Þ þ n i ðtÞ. The temporal change in n i (t) is driven by the imbalance among influx, proliferation, death, and outflux of the type i cells, each of which depends on the numbers of other cells n t ð Þ :¼ n DN t ð Þ; n DP t ð Þ; n SP4 t ð Þ; n cTEC t ð Þ; n mTEC t ð Þ ½ T , where T denotes transpose. While the influx may be independent of the number of type i cells, the other should, in nature, depend on the number of existing type i cells, n i (t). This allows us to generally represent the ordinary differential equations for n i (t) as where the influx should be non-negative, ϕ i n t ð Þ ð Þ! 0, whereas the marginalized rate of proliferation, death, and outflux, f i ðnðtÞÞ, can be either positive or negative. The actual value of f i ðn t ð ÞÞ is determined by the balance among proliferation, cell death, and outflux of type i cells. To obtain a minimal model with minimal complexity, we assume that both ϕ i n t ð Þ ð Þand f i ðnðtÞÞ are at most linear with respect to n(t) with possible constant time delays. Therefore, our ordinary differential equation model as a whole has, at most, quadratic nonlinearity. Considering reproducibility of the recovery dynamics after the irradiation and consistency with previously known molecular evidence, we obtained the whole model described as: a diagrammatic representation of which is shown in Fig. 2a. Based on this model with candidate parameter values as the initial condition, we conducted a nonlinear least square estimation of the whole parameter values in Eq. (1), fn tot i 0 ð Þg i2C , fω i g i2C , and p i f g i2C so that the whole model can reproduce all the experimental data at once ( Fig. 2b and Table 2). As shown in Fig. 2b, our model, Eq. (1), nicely reproduced the experimentally observed recovery, demonstrating that the interactions depicted in Fig. 2a sufficiently account for the dynamics. Moreover, to reevaluate the importance and statistical confidence of several parameters, we statistically estimated the potential variability of the estimated values by conducting a bootstrap parameter estimation (Figs. 2c and 3, and Supplementary Table 1). As shown in Fig. 3, most parameter values statistically fluctuate around single peak, whereas a few parameters, e.g., the influx rate of DN, ϕ 1 , have multiple peaks in their estimates.
DN thymocytes and cTECs form a negative feedback. Our estimated model indicates that DN thymocytes and cTECs form a negative feedback. DN cells marginally work to increase the number of cTECs because μ c in Eq. (1) is positive, whereas cTECs effectively inhibit the increase in DN cells because −μ 1 in Eq. (1) is negative (Fig. 2a). This negative feedback is the source of the overshooting behaviors in the recovery dynamics, and can account for slower onset of cTECs recovery, which lagged a few days behind DN cells.  Table 1. The right panel shows violin plots of the numbers of thymocytes and TECs without perturbation (n = 15 for thymocytes, n = 16 for TECs). c Typical flow cytometric profiles of the thymocytes after the sub-lethal dose radiation. Thymocytes were analyzed by staining with anti-CD4 and anti-CD8α. Percentage of each fraction is shown in the panels. d Typical flow cytometric profiles of TECs after the sub-lethal dose radiation. TECs (EpCAM + CD45 − TER119 − ) were analyzed by staining with a combination of UEA-1 lectin and anti-CD80. Percentages of UEA-1 + cells (mTECs) and UEA-1cells (cTECs) are shown in the panels.
These interactions inferred from the quantitative recovery data are also consistent with previously identified molecular evidence. On one hand, the positive interaction from DN thymocytes to cTECs may be interpreted as induction of cTEC proliferation by DN cells, evidenced by the fact that the number of mature cTECs decreases if DN differentiation is blocked at early stages 31,34 . On the other hand, our model suggests that cTECs work to decrease the number of DN cells. This negative interaction is a marginal effect of induced cell death, induced differentiation from the DN to the DP stages, and inhibition of DN proliferation by cTECs. This negative regulation of DN cells by cTECs is consistent with the lineage commitment of DN cells to the DP stage mediated by cTECs in Notch1-Delta-like4-dependent manner 27,28 . It should be noted, however, that our model does not exclude other possibilities of additional molecular interactions, as long as their marginal influences are consistent with the diagram in Fig. 2a.
To further analyze the consistency of our model with the underlying dynamics of DN subpopulations (DN1, DN2, DN3, and DN4), we additionally quantified the dynamics of these populations after irradiation ( Supplementary Fig. 4). We also modified Eq. (1) (denoted here as a coarse-grained model) to include DN subpopulations (denoted as a detailed model and shown in Methods), the parameter values of which were similarly estimated. As demonstrated in Fig. 4a, b, the detailed model reproduces the dynamics of the DN subpopulations ( Fig. 4a) with only small deviation from the coarse-grained model in which DN subpopulations are lumped together (Fig. 4b). We should mention that our estimates of DN1 and DN2 subpopulations can be overestimates because an additional cell surface marker, CD117, is required to exclude non-T-lineage fractions 35 . This overestimate, however, has little effect on the inferred dynamics of the total DN population (Fig. 4b), because the major fraction of DN cells consists of DN3 and DN4 cells. The estimated parameter values were also consistent with those of the coarse-grained model except for the DN1 influx rate ϕ 1 estimate, which was much smaller than that of coarse-grained model. Because the peak other than that around the optimal value in the bootstrap estimate of ϕ 1 was at the lower bound of its estimation range (Fig. 3), a value smaller than the lower bound may also reproduce the same recovery dynamics. To verify whether the estimate of ϕ 1 in the detailed model is also consistent with the coarse-grained model, we simulated coarse-grained model by replacing the value of ϕ 1 in the model with the estimate from the detailed model. As shown in Fig. 4b, the trajectories were almost unaffected by this replacement. Moreover, from the biological viewpoint that the number of the influx DN progenitors is quite small 36 , this value of ϕ 1 is also reasonable. Altogether, analysis of the detailed model revealed that the smaller value of ϕ 1 , which cannot be selected only from the analysis of the coarse-grained model, is more relevant.   Days after irradiation  0  1  4  7  9  11  12  13  14  15  17  19  49  Number of samples (thymocytes)  4  6  3  3  3  3  2  3  3  3  3  6  3  Number of samples (TECs)  4  6  3  6  3  3  2  3  6  3  3  6  3 DP recovery by temporal increase in proliferation rate. The kinetic component characteristic to the DP dynamics is its much faster recovery compared to DN cells (Fig. 1b), which strongly suggests that the DP recovery is achieved by self-proliferation rather than by the influx from the DN population. However, evidence about the self-proliferation ability and speed of DP cells is inconsistent and may depend on strains 36,37 ; some studies showed that DP cells proliferate little 22,38 while others have suggested that DP cells can proliferate faster than other types of thymocytes 37,39 . Our model coordinates these properties with auto-inhibitory regulation of the DP proliferation, represented by the logistic term θ 2 ð1 À n DP ðtÞ=K 2 Þ in Eq. (1), which can realize fast proliferation during the recovery period and slowdown at the steady state. Nevertheless, such auto-inhibitory regulation in DP proliferation has not yet been reported.
To experimentally verify this prediction by our model, we estimated the fraction of proliferating DP cells under the same condition as in Fig. 1a by staining the DP population with proliferation marker Ki67 (Fig. 4c). We observed that the fraction of the proliferating DP cells transiently increased and peaked at day 7 after irradiation, coinciding perfectly with the timing of exponential increase in DP cells during recovery. Selfproliferation ceased when the number of the DP cells recovered to the normal population size before the irradiation. This result strongly supports that the proliferation rate of DP cells is inhibited by total population size to maintain homeostasis. Further, this autoregulatory mechanism is consistent with the previous observations that DP cells proliferate little when their numbers are at the steady state 38 .
While the autoregulatory proliferation of DP cells is necessary for reproducing fast recovery, it cannot solely account for the overshooting behavior of DP cells, which suggests that other cells regulate DP cells. Supported by well-established evidence that cTECs engage in positive selection of DP cells, our model includes a negative influence of cTECs to DP cells with a time delay, which can nicely reproduce the overshoot of DP cell count. This negative interaction with a time delay can be interpreted as the marginal effect of an induced apoptosis of DP cells with nonfunctional T cell receptors (TCRs) and the differentiation of DP cells into SP cells upon apoptosis rescue. The existence of the time delay may be interpreted by the sequential and multiple interactions of DP cells with cTECs that are required for positive selection.
Our model estimates that the stable rate of DP cells to differentiate into CD4 SP cells, r 24 μ 2 n Ã cTEC , ranges from 6.0 × 10 −2 to 2.5 × 10 −1 (day −1 ), overlapping the range of the previous estimates of 1.2 × 10 −2 to 9.9 × 10 −2 (day −1 ) ( Table 2). The estimated value of r 24 varied from 4.9% to 30%, within the range of the previous estimates that 0.02-65% of DP cells survive and differentiate into CD4 SP via positive and negative selections ( Table 2). This result supports the interpretation that r 24 is the fraction of rescued DP cells that differentiated into CD4 SP, and the remaining fraction 1 − r 24 of DP cells undergoes apoptosis. However, we should note that an apoptosis rate cannot be directly estimated solely by population size dynamics. This may be the major reason why the estimated fraction of the rescued DP cells varies in our and previous studies. DP and CD4 SP thymocytes incoherently regulate mTEC recovery. Compared with other thymocytes and TECs, CD4 SP cells recovered much slower, with less pronounced overshooting (Fig. 1b). This slow recovery of CD4 SP cells is consistent with their lack of proliferation capacity 13,37 , which leads to prolonged recovery. The CD4 SP dynamics can be reproduced by assuming no proliferation and mTEC-dependent death and outflux Àμ 4 n mTEC ðtÞ, which may represent the negative selection of SP cells by mTECs (Fig. 2a).
In contrast, the mTEC recovery was initiated almost concurrently with that of cTEC (Fig. 1b). While interactions with CD4 SP cells have been proven essential for the maturation of mTECs 40 , the prolonged CD4 SP recovery is insufficient for We note that our point estimate of the DP residence time 230 h may be an overestimate, although the previous estimates overlap the statistically confident range of the values obtained by our bootstrap analysis. This is because the residence time was estimated only from the output flux rate, due to the fact that the apoptosis rate cannot be estimated in our model reproducing earlier onset of mTEC recovery. Our model incorporates an auto-inhibitory regulation of mTEC proliferation r m ð1 À n mTEC t ð Þ=K m Þ and its negative regulation by DP cells with a time delay Àγ mp n DP ðt À τ m Þ as in Eq. (1). The auto-inhibitory regulation is necessary because without it, we obtained biologically inconsistent parameter values in mTEC dynamics (Fig. 5a, e). The negative regulation by DP cells is also responsible for mTEC overshooting. Preceding experimental investigations 5,41 support these mechanisms. Metzger et al. reported that the percentage of Ki67hi mTECs increases only after the depletion of mTECs 41 , suggesting auto-inhibitory regulation. Based on a depletion experiment of DP cells, Dudakov et al. suggested that DP cells negatively regulated TEC proliferation in an IL22dependent manner 5 . However, the DP-dependent regulation was not the sole interaction that could explain the early onset of mTEC recovery. We also found that a DN-dependent regulation could reproduce it (Fig. 5b, f). However, this possibility was excluded in our model because we lack molecular evidence supporting the long-range interaction from DN cells to mTECs, which reside in spatially segregated areas of a thymus.
Along with regulated proliferation, our model assumes reciprocal regulations between mTECs and CD4 SP cells to account for evidence that mTEC maturation is also related to CD4 SP cells. According to Williams et al. 33 , mTECs express ligands CD80 and CD86 and a receptor, CD40; the corresponding ligand and receptor of CD4 SP cells are mainly CD28 and CD40L. A knockout of CD80, CD86, and CD40 was shown to decrease the number of mTECs and double the number of CD4 SP cells. We substituted smaller values of μ 4 and ϕ m4 than the estimated values into our model to reproduce the experiment in ref. 33 by assuming that the knockout of CD80, CD86, and CD40 corresponds to this substitution. The result qualitatively reproduced the knockout mutant result in ref. 33 ; the stable number of CD4 SP cells doubled whereas the number of mTECs was decreased, as shown in Fig. 4d.

Discussion
From quantitative time-series data of thymocytes and TECs recoveries after sub-lethal X-ray irradiation, we constructed a mathematical model for the recovery dynamics of thymocytes and TECs. The model reproduces the transient dynamics of the cell population sizes fairly well, and most of the interactions identified by the modeling are consistent with known molecular evidence.
Since previous modeling works on quantitative characterizations of thymocyte development focused only on thymocyte dynamics, our work, which additionally includes both the dynamics of and the interactions with TECs, can be viewed as an extension of those works 17,[19][20][21]42,43 . We validated that the estimated parameter values of thymocytes in our model are mostly consistent with those estimated in previous works ( Table 2). Few parameter value mismatches may also be attributed to differences in the experimental setting and conditions. However, it should be noted that we compared the apparent proliferation rates of our model to previous estimates, which are the marginal rates of population size change owing to the imbalance between proliferation and apoptosis because pure rates of proliferation or apoptosis cannot be estimated from our model. To reveal dynamical change in apoptosis rates, we must develop a new approach that combines our reciprocal regulation model with experimental methods that can directly quantify thymocyte proliferation, apoptosis, and differentiation 14,44 . Further, the parameters of TECs are the first to be estimated by modeling, and should be verified by independent research. Particularly, damage caused by irradiation can affect the thymic tissue structure, which may result in a systematic bias when counting TECs 45,46 . While this systematic bias is effectively absorbed in our model by parameters scaling, this potential scaling must be considered when we compare our TEC parameter value estimates with others. Moreover, to assess this problem more carefully, we must develop a new image analysis method that can accurately detect and count cells in 3D tissue images obtained by advanced imaging and tissue clearing 47,48 .
Thymic crosstalk includes various signaling pathways, indicating complex regulations behind the population size control of thymocytes and TECs. Because of this complexity, our model may contain missing interactions or possibly different regulations, some of which were tested during our model identification process. Such possibilities cannot be excluded by the limited amount of data alone; therefore, we employed previously obtained molecular-biological evidence and quantitative estimates to evaluate the possible models.
For example, cTECs rescue DP thymocytes from apoptosis via positive selection, which leads to increase in the DP population size. Concurrently, positive selection also induces differentiation of thymocytes from the DP stage to the SP stage, causing DP population size to decrease. These contradicting interactions introduce the possibility that cTECs increase the DP thymocyte population size, rather than decreasing it as assumed in our model. We examined this possibility by introducing the increasing effect of the DP population size by cTECs and concluded that the decreasing effect assumed by our model is more valid because the model with the increasing effect resulted in much higher parameter values than expected based on the previous works (Fig. 5c, g).
We also investigated a model in which DP thymocytes contribute to the recovery of both mTECs and cTECs 5 . We found that the estimated parameter for the interaction from DP cells to cTECs was almost 0 (Fig. 5d, h), which does not support a major contribution of DP cells to cTECs recovery under our experimental condition. d In silico evaluation of the impact from the disturbed crosstalk between SP4 thymocytes (light green) and mTECs (brown). Thick solid curves are simulated trajectories of SP4 thymocytes and mTECs with parameter values mimicking the experimental condition in ref. 33 , γ 4 = 5.0 × 10 −6 and ϕ m4 = 0. The thin broken curves are those obtained with the optimal parameter values used in Fig. 2b for comparison.
Our model can explain the mechanisms by which specific dynamics appear in recovery dynamics and their potential biological functions; overshoots of DN thymocytes and cTECs may originate from negative feedback between them and may contribute to prompt recovery from various perturbations affecting thymocyte and TEC numbers. Similarly, the disinhibition of DP proliferation upon DP population size decrease facilitates the swift recovery of DP cells, which could not be achieved solely by the influx of DN cells, as they have a much smaller population size than DP cells. Our model provides an integrative view of Fig. 5 Possible regulatory mechanisms that are capable of reproducing the recovery dynamics of the data, but that are biologically less relevant than the proposed model shown in Fig. 2a. Differences between each model and the one shown in Fig. 2a are designated by red solid lines if additionally included or red broken lines if excluded. The equations corresponding to the models are shown in Methods. The model in a excludes the autoinhibitory regulation of mTECs. The model in b includes an interaction from DN cells to mTECs influx instead of the inhibition from DP cells. The model in c assumes that cTECs promote DP cell proliferation, rather than inducing DP cell differentiation or cell death. The model in d includes an inhibitory regulation of cTECs from DP cells similar to that of mTECs. e, f, g, and h show corresponding trajectories of the models in a, b, c, and d. thymic crosstalk as a regulatory network and serves as a starting point for comprehensively understanding homeostasis in thymic development. However, our model still has room for future improvement by accommodating more detailed information on the cellularity of the thymic resident cells, such as B cells, dendritic cells, and thymic endothelial cells. These cells may have different roles in the dynamic regulation of thymic homeostasis than thymocytes and TECs, although we did not explicitly include them by presuming that their effects to the number of thymocytes or TECs are relatively small or constant, which was implicitly modeled by the constant parameters in our model. Actually, BMP4 production by endothelial cells after irradiation, which can contribute to TECs recovery, was reported constant when normalized by the size of thymus 49 . Explicitly incorporating these cells may be crucial to extending our model to other experimental settings as well as for deriving a more integrative and comprehensive model of thymic development and homeostasis. Among others, the repertoires of thymocytes are particularly relevant. TECs are responsible for controlling the number of thymocytes as well as for selecting thymocytes with appropriate repertoires. An upcoming challenge may be integratively modeling and analyzing thymic homeostasis, both in cell numbers and repertoires by combining quantitative measurement and high-throughput sequencing 50 .

Methods
Ethics statement. Animals used in the present study were maintained in accordance with the "Guiding Principles for Care and Use of Animals in the Field of Physiological Science" set by the Physiological Society of Japan. All animal experiments were approved by the Animal Research Committees of RIKEN.
Estimation of the fraction of proliferating DP cells. Thymocytes were pretreated with anti-CD16 and CD32 (Biolegend) and subsequently stained with anti-CD4 and anti-CD8 antibodies in phosphate buffered saline containing 3% FBS. The cells were fixed and permeabilized with Foxp3/Transcription Factor Staining Buffer Set (eBioscience) according to the manufacturer's protocol. After fixation and permeabilization, the cells were stained with a PE-labeled anti-Ki67 antibody (Biolegend) and subsequently analyzed by Canto II (BD).
Statistics and reproducibility. For each experimental condition, we measured the numbers of thymocytes and TECs from at least two mice, independently. The variations of parameter values for the mathematical model were analyzed by the bootstrap method described in the section "Confidence interval by bootstrap".
Mathematical modeling of thymocyte and TEC dynamics. We assume that the total number of the type i cells, n tot i , is the sum of cells dying by irradiation n x i and survived or newly generated cells n i : where C is the set of the cell types.
We describe the decrease in the irradiated cells by an exponential decay, which assumes that cells die at a constant rate ω i after irradiation: In the model, n tot i 0 ð Þ represents the initial population size of type i cells and p i is assumed to be the fraction of survived cells at t ≤ 0 as Given these initial conditions, the model of Eq. (1) was implemented on MATLAB (R2018a; The MathWorks, Natick, MA) and was numerically simulated by 'dde23' function or on Mathematica (version 11.2; Wolfram research, Champaign, Illinois) and simulated by 'NDSolve' function.
Parameter estimation. In the parameter estimation, ω i , n tot i 0 ð Þ, p i , and all parameters appearing in Eq. (1) were simultaneously estimated. Parameters were estimated by minimizing the sum of the squares of difference between the logarithms of the observed data and simulated values of the model. Because the orders of the parameters are different, and this caused difficulty in the minimization, we decomposed the parameters as θ ¼ θ c θ p , where θ c is a coefficient vector to estimate, θ p is a constant vector of a power of 10, and denotes elementwise multiplication. For the observed time points t Ã ¼ ½t 1 ; Á Á Á ; t m and the corresponding data points N i (t j ) for all i 2 C, the estimated parameter setθ was obtained by solving To solve this minimization, we used the 'lsqnolin' function in MATLAB Optimization Toolbox in which parameters were estimated by Trust Region Reflective method. The initial parameter values in the estimation were given, so that the result converges to moderate values considering the results of related previous works. The searching range of each parameter, except for p i , r 1 , and r 24 , was set between 10 and 0.1 times the initial value. Since p i , r 1 , and r 24 represent fractions, their searching ranges were set between 0 and 1. The symbols, descriptions, and estimated values of the parameters are listed in Supplementary  Table 1.
Confidence interval by bootstrap. We calculated the confidence intervals of the estimated parameter values by a bootstrap method 51 .
First, for type i cells, we modeled the statistical variation of the data points using a Gaussian random variable ε i $ N ð0; σ 2 i Þ with mean 0 and variance σ 2 i as ln N i t ð Þ ð Þ¼lnðn tot i ðt; b θÞÞ þ ε i : We estimated σ 2 i by the sample variance aŝ ½ln N i ðt j Þ À ln n tot i ðt j ;θÞ 2 : We obtained the kth bootstrapped sample of the time point t j , N b k i ðt j Þ by using a random number ε k i;j $ N ð0;σ 2 i Þ as ln N b k i t j ¼ ln n tot i t j ;θ þ ε k i;j : The kth bootstrapped parameter setθ b k was obtained by solving the same optimization problem of the previous section by replacing the data with the kth The total number of the bootstrapped samples generated was B = 1000. The two-sided α × 100% confidence interval of the lth parameter was calculated as ½θ   Supplementary Fig. 3. The trajectories of cells obtained from 100 samples of the bootstrap parameter sets are shown in Fig. 2c.
Detailed model of DN thymocytes. We additionally measured dynamic changes in the population sizes of DN1, DN2, DN3, and DN4 cells after irradiation.
To estimate the DN subpopulation dynamics in the original data (Fig. 1b), we utilized the DN subpopulations data as follows. First, we calculated the average proportions of the DN subpopulations at each time point. Subsequently, assuming that the dynamics of the DN subpopulation proportions were the same as the original data, we multiplied the number of DN cells from the original data with the calculated DN subpopulation proportions at each time point. At time points where we did not have corresponding DN subpopulation data (days 12 and 14), we used the average proportions of neighboring time points (days 11 and 13 for day 12, and days 13 and 15 for day 14).
Employing the obtained estimates of the DN subpopulation dynamics, we estimated the parameter values of the following detailed model of DN1, DN2, DN3, DN4, DP, and cTEC: The parameter estimation procedure was the same as for the coarse-grained model. Because the detailed model has parameters common to the coarse-grained model, ϕ 1 and the model parameters of DP and cTECs, except for r DN4 , μ DN4 , and μ cTEC;i , we first fixed those parameter values to the estimates from the coarsegrained model and estimated the remaining parameter values. However, the detailed model with the estimated parameter values did not reproduce the DN1 dynamics ( Supplementary Fig. 5). To obtain the parameter values capable of reproducing the dynamics of all cell types, we estimated parameter values including ϕ 1 while other common parameter values were fixed (Fig. 4a and Supplementary  Table. 2).