武汉新启迪生物科技有限公司

新启迪-您的生物科研好伙伴！2008-2024

Wuhan Xinqidi Biotech Co.Ltd

本企业通过iso9001质量体系认证

定量分析揭示了小鼠胸腺细胞和胸腺环境恢复动力学的相互作用规律

二维码

发表时间：2019-12-01 13:25作者：武汉新启迪Xinqidi来源：www.qidibio.com

定量分析揭示了小鼠胸腺细胞和胸腺环境恢复动力学的相互作用规律

摘要

胸腺串扰是胸腺细胞与胸腺环境之间的一套相互调节机制，对组织合适的胸腺细胞发育以及从各种外源性侮辱中恢复胸腺具有重要意义。本文在对辐射诱导的胸腺细胞和胸腺上皮细胞的恢复动力学进行定量分析和建模的基础上，推导出胸腺细胞和胸腺上皮细胞相互作用对胸腺串扰的影响。该分析确定了与已知分子证据相一致的调节相互作用，并揭示了它们在恢复过程中的动态作用。此外，分析还预测，并随后的实验证实，以前未被承认的调节CD4+CD8+双阳性胸腺细胞，这暂时提高了他们的增殖速度，因为他们的群体规模缩小。我们的模型建立了一个关键步骤，以动态理解胸腺串扰作为一个调节网络系统。

导言

胸腺是一个器官，负责产生大量的T细胞，并有适当的储备。1...然而，它对压力、病毒感染、辐射和其他刺激的侮辱相对敏感。2,3...虽然正常情况下，健康动物的胸腺可以从这些损伤中恢复，但胸腺恢复过程相对较长的过程可能会损害T细胞介导的免疫，因为恢复期内天真的T细胞储备的补充减少。3,4.

摘要以小鼠亚致死剂量辐射作为小鼠毁损后胸腺再生的实验模型。5,6...电离辐射也被广泛用于造血移植和癌症治疗。7,8全身照射可引起急性胸腺损伤和胸腺功能恢复缓慢。一些研究表明，辐射不仅降低了胸腺细胞的细胞密度，而且降低了胸腺上皮细胞(TECs)的细胞密度，而TECs是胸腺环境的主要成分。5,9,10...因为胸腺细胞和TECs之间的相互作用支持了胸腺的形成。11了解胸腺恢复需要描述胸腺细胞与TECs之间的相互调节。

与此同时，各种追踪、扰乱和量化参与这些事件的细胞的技术使我们能够定量地描述它们的动态。12,13,14,15...通过将数学模型与这些定量数据相结合，胸腺发生的动态方面被提炼成详细的动力学信息，如增殖率、死亡率和分化率。12,16...Mehr等人17利用常微分方程建立了第一个胸腺细胞发育动力学模型。18...自这项开创性工作以来，通过考虑胸腺细胞的详细细胞结构和发育状态，以及结合不同的实验条件，逐步完善了胸腺发生的动力学模型。19,20,21,22,23.

然而，以往的工作只关注胸腺细胞。胸腺发育和胸腺恢复不是胸腺细胞自主的，而是由胸腺环境支持的。在过去的十年里，我们积累了大量的分子生物学证据，证明胸腺环境本身是由胸腺串扰、胸腺细胞和胸腺环境之间的双向相互作用来维持的。11,24,25...在构成胸腺环境的几个细胞中，皮质和髓质胸腺上皮细胞(cTECs和mTECs)在诱导和控制胸腺细胞的增殖、凋亡和谱系承诺中起着不可或缺的作用。11,26,27,28,29,30...胸腺细胞也通过调节其成熟和增殖来调节TECs。5,31,32,33...尽管胸腺串扰对胸腺发生和胸腺恢复有着明显的相关性和重要性，但胸腺细胞与TECs之间相互调节的动力学方面尚未得到澄清。

本研究结合数学模型，结合对照射后胸腺细胞和TECs数量的定量测量，研究了胸腺细胞和TECs的联合动力学。恢复动力学由我们的数学模型再现，在数学模型中，我们识别了胸腺细胞与TECs之间的相互作用，这些相互作用与恢复相关，并且与胸腺串扰相一致。此外，我们还证明了该模型为种群规模的动态变化机制提供了一个解释。特别是，我们的模型预测，并随后的实验证实，以前未被承认的调节CD4+CD8+双阳性(DP)胸腺细胞，这暂时提高了他们的增殖速度，因为他们的群体规模缩小。

结果

胸腺细胞和TECs恢复动力学的定量研究

为了定量研究胸腺细胞与TECs之间的动力学关系以及胸腺恢复机制的建立，我们采用亚致死剂量4.5Gy对胸腺细胞和TECs群体进行了人工扰动，并测定了照射后3周内胸腺细胞和TECs群体大小的动态变化。1A和表1)。图形1B总结细胞数量的变化，这些变化是根据胸腺细胞的常规标记进行分类的(如图所示)。1C和补充图。1A)和TECs(图1.1D和补充图。1B)。图形1B结果表明，照射后各类型胸腺细胞和TECs均以不同速率呈指数下降。随后，胸腺细胞和胸腺细胞都在10天内开始恢复，CD4−CD8−双阴性(DN)胸腺细胞和cTECs在5天内开始恢复，而CD4+CD8−单阳性(SP)胸腺细胞和mTECs需要更长的时间间隔，反映了胸腺细胞从DN到CD4+SP(SP4)细胞从cTECs到mTECs的时间顺序。

在恢复后，除SP细胞外，其他细胞的种群规模在15天左右达到高峰，并最终恢复到稳定的数量，几乎相当于或至少相当于照射前原始种群规模的一半。这种过度行为表明胸腺细胞和TECs的数量是通过相互作用动态和相互调节的(表)。1).

表1辐照后各时间点的样品数。

全尺寸表

数学模型可以再现恢复动力学

为了推断动力学背后的调节相互作用，我们用常微分方程建立了胸腺细胞和TECs种群动力学的数学模型，其中明确包括五种细胞类型：

$i \in C = {D N, D P, S P 4, c T E C, m T E C}$ ...考虑辐射对细胞的急性影响，细胞类型的总数。i时时t(日)， $n_{i}^{t o t} (t)$ 分解为两部分； $n_{i}^{x} (t)$ 以指数的方式表示受照射的死亡细胞，并且ni(t)表示辐照后存活或新生的细胞。 $n_{i}^{x} (t)$ 假设以一个恒定的速率呈指数下降，ωi(日)−1)，如 $n_{i}^{x} (t) = n_{i}^{t o t} (0) (1 - p_{i}) e^{- ω_{i} t}$ ，在哪里pi是辐照后存活细胞的比例；我们模拟了ni(t)常微分方程。因此，单元格类型的总数i, $n_{i}^{t o t}$ ，我们通过实验观察到，它被描述为 $n_{i}^{t o t} (t) = n_{i}^{x} (t) + n_{i} (t)$ .

时间变化ni(t)是由这种类型的内流、增殖、死亡和外流之间的不平衡所驱动的。i单元，每个单元取决于其他单元的数量。 $n (t) := {[n_{D N} (t), n_{D P} (t), n_{S P 4} (t), n_{c T E C} (t), n_{m T E C} (t)]}^{T}$ ，在哪里T表示转置。而流入量可能与类型的数量无关。i在本质上，单元格应该取决于现有类型的数量。i细胞，ni(t)。这使得我们可以一般地表示ni(t)

\frac{d n_{i} (t)}{d t} = ϕ_{i} (n (t)) + f_{i} (n (t)) n_{i} (t),

在非负数的情况下， $ϕ_{i} (n (t)) \geq 0$ 而被边缘化的扩散、死亡和外流率， $f_{i} (n (t))$ ，可以是积极的，也可以是消极的。实际价值 $f_{i} (n (t))$ 取决于增殖、细胞死亡和类型外通量之间的平衡。i细胞。为了得到一个复杂度最小的最小模型，我们假设两者都是 $ϕ_{i} (n (t))$ 和 $f_{i} (n (t))$ 最多是线性的，相对于n(t)可能存在恒定的时间延迟。因此，我们的常微分方程模型作为一个整体，最多具有二次非线性。考虑到辐照后回收动力学的重现性和与已知分子证据的一致性，我们得到了描述如下的整个模型：

\frac{d n_{D N} (t)}{d t} = ϕ_{1} + (δ_{1} - μ_{1} n_{c T E C} (t)) n_{D N} (t), \frac{d n_{D P} (t)}{d t} = μ_{1} r_{1} n_{c T E C} (t) n_{D N} (t) + {θ_{2} (1 - \frac{n_{D P} (t)}{K_{2}}) - μ_{2} n_{c T E C} (t - τ_{2})} n_{D P} (t), \frac{d n_{c T E C} (t)}{d t} = ϕ_{c} + (- δ_{c} + μ_{c} n_{D N} (t)) n_{c T E C} (t), \frac{d n_{S P 4} (t)}{d t} = μ_{2} r_{24} n_{c T E C} (t - τ_{2}) n_{D P} (t) - μ_{4} n_{m T E C} (t) n_{S P 4} (t), \frac{d n_{m T E C} (t)}{d t} = ϕ_{m} + ϕ_{m 4} n_{S P 4} (t), + {r_{m} (1 - \frac{n_{m T E C} (t)}{K_{m}}) - γ_{m p} n_{D P} (t - τ_{m})} n_{m T E C} (t),

(1)

图示如图所示。2A...在此模型的基础上，以候选参数值为初始条件，对方程中的整个参数值进行了非线性最小二乘估计。(1), ${n_{i}^{t o t} (0)}_{i \in C}$ , ${ω_{i}}_{i \in C}$ ，和 ${p_{i}}_{i \in C}$ 这样整个模型就可以一次复制所有的实验数据(图)。2B和表2)。如图所示。2B，我们的模型，Eq。(1)，很好地再现了实验观察到的恢复，证明了图中描述的相互作用。2A充分解释了动态。此外，为了重新评估几个参数的重要性和统计置信度，我们通过自举参数估计(图1)对估计值的潜在变异性进行了统计估计。2C和3，及补充表1)。如图所示。3，大多数参数值在单峰附近有统计上的波动，而一些参数，例如DN的流入率，ϕ1，在它们的估计中有多个峰值。

表2.估计动力学率与以往研究的比较。

全尺寸表

DN胸腺细胞和cTECs形成负反馈

我们估计的模型表明DN胸腺细胞和cTECs形成负反馈。DN单元可以稍微增加cTECs的数量，因为μc在等式里。(1)是阳性的，而cTECs有效地抑制了DN细胞的增加，因为−μ1在等式里。(1)是阴性的。2A)。这种负反馈是恢复动力学中超调行为的来源，可能是cTECs恢复缓慢的原因之一，cTECs的恢复滞后于DN细胞数天。

从定量恢复数据推断的这些相互作用也与先前确定的分子证据相一致。一方面，DN胸腺细胞与cTECs之间的正向相互作用可能被解释为DN细胞对ctc增殖的诱导作用。31,34...另一方面，我们的模型表明cTECs可以减少DN细胞的数量。这种负性相互作用是诱导细胞死亡、诱导DN向DP期分化和cTECs抑制DN增殖的边缘效应。cTECs对DN细胞的这种负性调节，与cTECs介导的dp期的遗传承诺是一致的，呈Notch 1-Delta样4依赖关系。27,28...但是，应该指出的是，我们的模型并不排除其他额外分子相互作用的可能性，只要它们的边际影响与图中的图表一致。2A.

为了进一步分析我们的模型与DN亚群体(DN1、DN2、DN3和DN4)潜在动力学的一致性，我们还量化了照射后这些种群的动态(补充图)。4)。我们还修改了方程。(1)(此处表示为粗粒度模型)，包括DN子种群(表示为详细模型，并在方法中显示)，其参数值也是类似估计的。如图所示。4A，b，详细的模型再现了DN亚群体的动态。4A)与粗粒度模型只有很小的偏差，在粗粒度模型中DN子种群集中在一起(图1)。4B)。我们应该指出，我们对dn1和dn2亚居群的估计可能过高，因为需要额外的细胞表面标记cd 117来排除非T系分数。35...然而，这一过高估计对总DN种群的推断动态几乎没有影响(如图所示)。4B)，因为DN细胞的主要部分由DN3和DN4细胞组成。除dn1流入率外，估计的参数值与粗粒度模型的参数值也是一致的。ϕ1估计值比粗粒度模型小得多.的最优值以外的峰值。ϕ1在其估计范围的下界(图1)。3)，小于下限的值也可以再现相同的恢复动态。以核实ϕ1在详细的模型中，与粗粒度模型也是一致的，我们通过替换粗粒度模型的值来模拟粗粒度模型。ϕ1在模型中用详细的模型进行估计。如图所示。4B，轨迹几乎不受这一替代物的影响。此外，从生物学的角度来看，流入的DN祖细胞数量相当少。36的价值ϕ1也是合理的。总之，对详细模型的分析表明，ϕ1它不能仅从粗粒度模型的分析中选择，更具有相关性。

增殖率随时间增长而恢复的DP

与DN细胞相比，DP动力学的动力学成分特征是其恢复速度快得多(如图所示)。1B)，这有力地表明，DP的恢复是通过自我增殖而不是通过DN人群的流入来实现的。然而，关于dp细胞自我增殖能力和速度的证据并不一致，可能取决于菌株。36,37一些研究表明dp细胞增殖缓慢。22,38而另一些人则认为dp细胞的增殖速度比其他类型的胸腺细胞更快。37,39...我们的模型将这些特性与dp增殖的自抑制调节相协调，用Logistic项表示。 $θ_{2} (1 - n_{D P} (t) / K_{2})$ 在等式里。(1)，能实现恢复期的快速增殖和稳态的减速。然而，对DP增殖的这种自抑制性调控尚未见报道.

为了用我们的模型验证这一预测，我们估算了在相同条件下增殖的DP细胞的比例，如图所示。1A用增殖标记Ki 67对DP群体进行染色。4C)。我们观察到，在照射后第7天，增殖的DP细胞的比例瞬间增加，并在第7天达到高峰，完全符合恢复过程中DP细胞指数增长的时间。当DP细胞数量恢复到照射前的正常群体时，细胞的自我增殖就停止了。这一结果有力地支持了DP细胞的增殖速度被总种群大小所抑制，从而维持稳态。此外，这种自我调节机制与先前观察到的dp细胞在处于稳定状态时几乎不增殖是一致的。38.

虽然DP细胞的自我调节增殖是复制快速恢复所必需的，但它不能完全解释DP细胞的过度行为，这表明其他细胞调控DP细胞。在cTECs积极选择DP细胞的充分证据的支持下，我们的模型包含了cTECs对具有时滞的DP细胞的负面影响，可以很好地再现DP细胞计数的超调。这种与时间延迟的负性相互作用可以解释为非功能T细胞受体(TCRs)诱导DP细胞凋亡和诱导DP细胞向SP细胞分化对细胞凋亡的边缘效应。时间延迟的存在可以用DP细胞与cTECs的连续和多重相互作用来解释，而cTECs是正选择所必需的。

我们的模型估计DP细胞分化为CD4 SP细胞的稳定率， $r_{24} μ_{2} n_{c T E C}^{*}$ ，范围为6.0×10−2降至2.5×10−1(日)−1)，与以前估计的1.2×10的范围重叠−2至9.9×10−2(日)−1)(表)2)。估计值r24从4.9%到30%不等，在先前估计的范围内，0.02%-65%的dp细胞通过阳性和阴性选择存活并分化为CD4 sp(见表)。2)。这一结果支持了以下解释：r24是已获救的dp细胞分化为CD4 sp的比例，而其余的dp细胞则是1−。r24DP细胞发生凋亡。然而，我们要注意的是，凋亡率不能仅仅用种群规模动态来直接估算。这可能是被拯救的DP细胞的估计部分在我们和以前的研究中不同的主要原因。

DP和CD4 SP胸腺细胞不一致调节mTEC恢复

与其他胸腺细胞和TECs相比，CD4 SP细胞的恢复要慢得多，并且没有明显的过度(图1)。1B)。CD4 SP细胞的这种缓慢恢复与其缺乏增殖能力是一致的。13,37导致长期的恢复。假设没有增殖和mTEC依赖的死亡和外通量，CD4 SP的动力学可以再现。 $- μ_{4} n_{m T E C} (t)$ ，可能代表了mTECs对SP细胞的阴性选择(图1)。2A).

相比较而言，mTEC的恢复几乎与CTEC的恢复同步进行(如图所示)。1B)。然而，与CD4 SP细胞的相互作用已被证明是mTECs成熟所必需的。40CD4 SP恢复时间较长，不足以早期恢复mTEC。我们的模型包含了一种对mtc增殖的自我抑制调节。 $r_{m} (1 - n_{m T E C} (t) / K_{m})$ 具有时滞的dp细胞对其负调控作用 $- γ_{m p} n_{D P} (t - τ_{m})$ 就像在Eq一样。(1)。自抑制调节是必要的，因为没有它，我们在mTEC动力学中得到了生物不一致的参数值。5A，e)。DP细胞的负调控也是mTEC过度调节的原因之一。前期实验研究5,41支持这些机制。Metzger等人报告说，Ki67hi mtecs的百分比只有在mtecs耗尽后才会增加。41，提示自抑制性调节。基于DP细胞耗尽实验，Dudakov等人。提示dp细胞以IL 22依赖方式负调控tc增殖。5...然而，依赖dp的调节并不是解释mTEC早期恢复的唯一相互作用。我们还发现依赖DN的调节可以复制它(图1)。5B，f)。然而，这种可能性在我们的模型中被排除了，因为我们缺乏支持从DN细胞到mTECs的长期相互作用的分子证据，而mTECs位于胸腺的空间分离区域。

随着受调控的增殖，我们的模型假设mTEC和CD4 SP细胞之间相互调节，从而证明mTEC的成熟也与CD4 SP细胞有关。根据Williams等人的说法。33mTECs表达CD 80、CD 86配体和CD 40受体，CD4 SP细胞的相应配体和受体主要为CD 28和CD40L。CD 80、CD 86和CD 40被敲除后，mTECs数量减少，CD4 SP细胞数增加一倍。我们用较小的μ4和ϕm4比我们的模型中的估计值更好地再现了实验结果。33假设CD 80、CD 86和CD 40的敲除对应于这一替代。该结果定性地再现了基因敲除突变体在参考文献中的结果。33如图所示，CD4 SP细胞的稳定数量增加了一倍，而mTECs的数量却减少了一倍。4D.

讨论

根据亚致死性X射线照射后胸腺细胞和TECs恢复的定量时间序列数据，建立了胸腺细胞和TECs恢复动力学的数学模型。该模型较好地再现了细胞种群大小的瞬态动力学，模型所确定的大多数相互作用与已知的分子证据是一致的。

由于以往关于胸腺细胞发育的定量表征的建模工作仅侧重于胸腺细胞动力学，我们的工作还包括胸腺细胞的动力学和与tecs的相互作用，因此我们的工作可以看作是这些工作的延伸。17,19,20,21,42,43...我们验证了模型中胸腺细胞参数的估计值与先前工作中的估计值基本一致(见表)。2)。参数值不匹配也可归因于实验环境和条件的差异。但是，值得注意的是，我们比较了我们的模型的表观增殖率和先前的估计，即由于增殖和凋亡之间的不平衡而导致的种群规模变化的边际率，因为我们的模型无法估计纯的增殖或凋亡率。为了揭示细胞凋亡率的动态变化，我们必须开发一种新的方法，将我们的倒数调控模型与直接量化胸腺细胞增殖、凋亡和分化的实验方法相结合。14,44...另外，TECs的参数是首先由模型估计出来的，需要独立的研究来验证。特别是辐射所造成的损伤会影响胸腺组织的结构，这可能会导致当计算tecs时会产生系统的偏差。45,46...虽然我们的模型有效地吸收了这一系统偏差，但当我们将TEC参数估计值与其他估计值进行比较时，必须考虑到这种潜在的标度。此外，为了更仔细地评估这一问题，我们必须开发一种新的图像分析方法，能够准确地检测和计数通过先进成像和组织清除获得的三维组织图像中的细胞。47,48.

胸腺串扰包括多种信号通路，表明胸腺细胞和TECs种群大小控制背后存在复杂的调控机制。由于这种复杂性，我们的模型可能包含缺失的交互，或者可能包含不同的规则，其中一些规则在我们的模型识别过程中进行了测试。仅凭有限的数据不能排除这种可能性；因此，我们利用先前获得的分子生物学证据和定量估计来评估可能的模型。

例如，cTECs通过正向选择从凋亡中拯救DP胸腺细胞，从而导致DP群体的数量增加。同时，正向选择还可诱导胸腺细胞从DP期向SP期分化，使DP群体规模减小。这些矛盾的相互作用引入了cTECs增加DP胸腺细胞数量的可能性，而不是像我们的模型所假设的那样减少它。我们通过引入cTECs对DP种群大小的增加效应来检验这一可能性，并得出结论：我们的模型所假定的递减效应更有效，因为在增加效应的情况下，基于先前的工作，模型的参数值比预期的要高得多(图1)。5C，g).

我们还研究了dp胸腺细胞促进mtecs和ctecs恢复的模型。5...我们发现DP细胞与cTECs相互作用的估计参数几乎为0。5D，h)，这不支持dp细胞在我们实验条件下对cTECs恢复的重要贡献。

我们的模型可以解释恢复动力学中出现特异性动力学的机制及其潜在的生物学功能，DN胸腺细胞和cTECs的超调可能源于它们之间的负反馈，可能有助于从影响胸腺细胞和TEC数量的各种扰动中迅速恢复。同样，由于DP细胞的数量比DP细胞小得多，减少DP增殖对DP细胞的抑制作用有利于DP细胞的迅速恢复，而Dp细胞的迅速恢复不能仅仅通过DN细胞的流入来实现。我们的模型提供了胸腺串扰作为调节网络的综合观点，并为全面理解胸腺发育中的稳态提供了一个起点。

然而，我们的模型仍有改进的空间，通过提供更详细的关于胸腺常驻细胞(如B细胞、树突状细胞和胸腺内皮细胞)细胞的信息。与胸腺细胞和TECs相比，这些细胞在胸腺内稳态的动态调节中可能具有不同的作用，尽管我们没有明确地将它们包括在内，因为我们假定它们对胸腺细胞或TECs数量的影响相对较小或恒定，这是由我们的模型中的常数参数隐式地模拟出来的。事实上，内皮细胞在辐射后产生BMP 4，这有助于细胞的恢复，据报道，当以胸腺的大小进行归一化时，BMP 4的生成是恒定的。49...明确地整合这些细胞对于将我们的模型扩展到其他实验环境以及建立一个更综合和全面的胸腺发育和稳态模型可能是至关重要的。除其他外，胸腺细胞的储存特别相关。TECS负责控制胸腺细胞的数量以及选择具有适当储备的胸腺细胞。即将面临的挑战可能是通过定量测量和高通量测序相结合，对胸腺内稳态进行整合建模和分析，无论是在细胞数量上还是在数量上都是如此。50.

方法

道德声明

本研究中使用的动物按照日本生理学会制定的“生理学领域动物护理和使用指导原则”进行维护。所有的动物实验都得到了RIKEN动物研究委员会的批准。

小鼠X线照射与流式细胞术

Balb/CA小鼠购自日本Clea。雌性小鼠(7周龄)接受X射线照射(4.5Gy)。在照射后的每个采样点，处死小鼠，用胸腺进行流式细胞术分析。每只胸腺被切下，在2ml RPMI-1640(Sigma-Aldrich，St.Louis，MO，USA)中轻轻搅拌，释放胸腺细胞进行流式细胞分析。测量的天数和抽样老鼠的数量如图所示。1A...其余胸腺组织在37°C下用RPMI 1640(Wako)中的自由酶消化30 min。用流式细胞仪分析胸腺基质富营养化组分，检测TEC群体。流式细胞术用抗CD 16和CD 32(Biolegend)预处理20 min，然后用含3%胎牛血清的磷酸盐缓冲液中荧光标记抗体染色。染色细胞用CantoⅡ(BD)分析。胸腺基质丰富组和胸腺细胞组的细胞之和决定胸腺细胞总数。TECS定义为CD 45-TER119-EpCAM+细胞.UEA-1染色分离mTECs和cTECs.DN胸腺细胞染色用CD 25和CD 44抗体染色分离阴性细胞。由于DN群体含有树突状细胞等其他小细胞群体，在稳态条件下，根据这些细胞在DN中所占的平均比例(16.6%)，从DN细胞的数量中减去这些细胞的数量(补充图)。1A)。PECy7-anti-CD4(clone RAM4-4，used as×200 dilution)，FITC-anti-CD4(clone RAM4-4，×200 dilution)，APCCy7-anti-CD8(clone 53-6.7，×200)，APCCy7-anti-CD45(clone 30 F-11，×200)，APCCy7-anti-TER119(clone TER-119，×200)，FITC-anti-EpCAM(BioLegend，clone G8.8，×400)，PE-anti-CD80(clone 16-10A1，×400)，Biotin-anti-mouse Ly-6G/Ly-6C(Gr-1)(×400)，从Biolegnd购买生物素抗小鼠/人CD45R/B 220(×400)、生物素抗小鼠(ter-119、×400)、生物素结合抗小鼠CD11b(×400)、PE抗小鼠/人CD 44(克隆IM 7，×400)、APC抗小鼠CD 25(克隆PC 61，×400)、链霉亲和素PE-Cyanine7(×400)和链霉亲和素-PECy 7(×400)。UEA-生物素(×400)来自载体实验室(伯林加梅，CA)。

增殖性DP细胞比例的估算

用抗CD 16和CD 32(Biolegend)预处理胸腺细胞，然后在含3%FBS的磷酸盐缓冲液中用抗CD4和抗CD8抗体染色。用Foxp 3/转录因子染色缓冲液集(EBioscience)固定和渗透细胞。固定和渗透后，用PE标记的抗Ki 67抗体(Biolegend)染色，然后用Canto II(BD)进行分析。

统计与重现性

在每个实验条件下，我们分别测量了至少两只小鼠胸腺细胞和TECs的数量。采用“置信区间自举”一节中的Bootstrap方法分析了数学模型参数值的变化。

胸腺细胞的数学模型及TEC动力学

我们假设类型的总数i细胞， $n_{i}^{t o t}$ ，是因辐射而死亡的细胞之和。 $n_{i}^{x}$ 存活或新生细胞ni:

n_{i}^{t o t} (t) = n_{i}^{x} (t) + n_{i} (t), i \in C := {D N, D P, S P 4, c T E C, m T E C},

哪里C单元格类型的集合。

我们描述了辐照细胞以指数衰减的方式减少，假定细胞以恒定的速度死亡。ωi辐照后：

n_{i}^{x} (t) = n_{i}^{x} (0) e^{- ω_{i} t} .

在模型中， $n_{i}^{t o t} (0)$ 表示类型的初始总体大小。i细胞和pi被认为是在t≤0 AS

\begin{array}{ll} n_{i} (t) = {\begin{matrix} n_{i}^{t o t} (0), t < 0 \\ p_{i} n_{i}^{t o t} (0), t = 0 \end{matrix}, \\ n_{i}^{x} (t) = {\begin{array}{l} 0, t < 0 \\ (1 - p_{i}) n_{i}^{t o t} (0), t = 0 \end{array} . \end{array}

给出了这些初始条件，建立了方程的模型。(1)在MATLAB(R2018a；MathWorks，Natick，MA)上实现，并用‘dde 23’函数或Mathematica(Version 11.2；Wolfram Research，Champaign，伊利诺伊州)进行数值模拟，并用NDSolve函数进行仿真。

参数估计

在参数估计中，ωi, $n_{i}^{t o t} (0)$ , pi，所有参数出现在方程中。(1)同时估计。通过最小化观测数据对数与模型模拟值之间差值平方之和来估计参数。由于参数的阶次不同，给最小化带来了困难，因此我们将参数分解为 $θ = θ_{c} \circ θ_{p}$ ，在哪里θc是要估计的系数向量，θp是幂为10的常数向量，并且 $\circ$ 表示元素乘法。对于观察到的时间点 $t^{*} = [t_{1}, \dots, t_{m}]$ 以及相应的数据点Ni(tj)所有人 $i \in C$ ，估计参数集 $\hat{θ}$ 通过求解

{\hat{θ}}_{c} = a r g min_{θ_{c}} \sum_{j = 1}^{m} \sum_{i \in C} {[\ln (n_{i}^{t o t} (t_{j}, θ_{c} \circ θ_{p})) - \ln (N_{i} (t_{j}))]}^{2}, \hat{θ} = {\hat{θ}}_{c} \circ θ_{p} .

为了解决这一极小化问题，我们在MATLAB优化工具箱中使用了“lsqnlin”函数，其中参数的估计采用了信赖域反射法。给出了估计中的初始参数值，使估计结果收敛到适中值，并考虑到相关工作的结果。每个参数的搜索范围，除pi, r1，和r24，设置为初始值的10至0.1倍。自pi, r1，和r24表示分数，它们的搜索范围设置在0到1之间。参数的符号、描述和估计值列在补充表中。1.

置信区间

我们用引导法计算了估计参数值的置信区间。51.

首先，类型i单元，我们用高斯随机变量模拟数据点的统计变化。 $ε_{i} \sim N (0, σ_{i}^{2})$ 平均0和方差 $σ_{i}^{2}$ 如

\ln (N_{i} (t)) = \ln (n_{i}^{t o t} (t, \hat{θ})) + ε_{i} .

我们估计 $σ_{i}^{2}$ 由样本方差

{\hat{σ}}_{i}^{2} = \frac{1}{m - 1} \sum_{j = 1}^{m} [\ln N_{i} (t_{j}) - \ln n_{i}^{t o t} (t_{j}, \hat{θ})]^{2} .

我们得到了k时间点的引导样本tj, $N_{i}^{b_{k}} (t_{j})$ 通过使用随机数 $ε_{i, j}^{k} \sim N (0, {\hat{σ}}_{i}^{2})$ 如

\ln N_{i}^{b_{k}} (t_{j}) = \ln n_{i}^{t o t} (t_{j}, \hat{θ}) + ε_{i, j}^{k} .

这个kTH自举参数集 ${\hat{θ}}^{b_{k}}$ 通过将数据替换为kTH自举样本 $N_{i}^{b_{k}} (t_{j})$ 如

\hat{θ}^{b_{k}} = \arg min_{θ} \sum_{j = 1}^{m} \sum_{i \in C} [\ln (n_{i}^{t o t} (t_{j}, θ)) - \ln (N_{i}^{b_{k}} (t_{j}))]^{2} .

生成的引导样本总数为B=1000。双面α的×100%置信区间lTH参数计算为 $[{\hat{θ}}_{l}^{(B (1 - α) / 2)}, {\hat{θ}}_{l}^{(B α / 2)}]$ ，在哪里 ${\hat{θ}}_{l}^{(x)}$ 是x的最小值l从引导样本中获得TH参数。我们用α=0.95。每个参数的置信区间在补充表中显示。1...自举估计的成对散点图如附图所示。3...从自举参数集的100个样本中获得的细胞轨迹如图所示。2C.

DN胸腺细胞的详细模型

此外，我们还测量了照射后DN1、DN2、DN3和DN4细胞群体大小的动态变化。

估计原始数据中的DN亚群体动态(图1)。1B)，我们如下所示地使用了DN亚群体数据。首先，我们计算了DN亚群体在每个时间点的平均比例。随后，假设DN亚群体比例的动态与原始数据相同，我们在每个时间点将原始数据中的DN单元数与计算的DN亚群体比例相乘。在没有相应DN亚群体数据的时间点(第12和14天)，我们使用相邻时间点的平均比例(第12天和第13天，第14天第13和15天)。

利用所获得的DN亚群体动力学估计，我们估计了DN1、DN2、DN3、DN4、DP和CTEC的详细模型的参数值：

\frac{d n_{D N 1}}{d t} = ϕ_{1} + (δ_{D N 1} - μ_{D N 1} n_{c T E C} (t)) n_{D N 1} (t), \frac{d n_{D N i}}{d t} = μ_{D N i - 1} n_{c T E C} (t) n_{D N i - 1} (t) + (δ_{D N i} + μ_{D N i} n_{c T E C} (t)) n_{D N i} (t), \frac{d n_{D P}}{d t} = r_{D N 4} μ_{D N 4} n_{c T E C} (t) n_{D N 4} (t) + {θ_{2} (1 - \frac{n_{D P} (t)}{K_{2}}) - μ_{2} n_{c T E C} (t - τ_{2})} n_{D P} (t), \frac{d n_{c T E C}}{d t} = ϕ_{c} + (- δ_{c} + \sum_{j = 1}^{4} μ_{c T E C, j} n_{D N j} (t)) n_{c T E C} (t), f o r i = 2, 3, 4.

参数估计过程与粗粒度模型相同.由于详细模型的参数与粗粒度模型相同，ϕ1和dp和ctecs的模型参数，除 $r_{D N 4}$ , $μ_{D N 4}$ ，和 $μ_{c T E C, i}$ 首先，我们将这些参数值固定到粗粒度模型的估计中，并估计出剩余的参数值。然而，带有估计参数值的详细模型没有再现DN1动力学(补充图)。5)。为了获得能够再现所有单元格类型的动态的参数值，我们估计了参数值，包括ϕ1而其他共同的参数值是固定的(如图所示)。4A和补充表。2).