Mathematical modeling of endogenous and exogenously administered T cell recirculation in mouse and its application to pharmacokinetic studies of cell therapies

Introduction: In vivo T cell migration has been of interest to scientists for the past 60 years. T cell kinetics are important in the understanding of the immune response to infectious agents. More recently, adoptive T cell therapies have proven to be a most promising approach to treating a wide range of diseases, including autoimmune and cancer diseases, whereby the characterization of cellular kinetics represents an important step towards the prediction of therapeutic efficacy.

Methods: Here, we developed a physiologically-based pharmacokinetic (PBPK) model that describes endogenous T cell homeostasis and the kinetics of exogenously administered T cells in mouse. Parameter calibration was performed using a nonlinear fixed-effects modeling approach based on published data on T cell kinetics and steady-state levels in different tissues of mice. The Partial Rank Correlation Coefficient (PRCC) method was used to perform a global sensitivity assessment. To estimate the impact of kinetic parameters on exogenously administered T cell dynamics, a local sensitivity analysis was conducted.

Results: We simulated the model to analyze cellular kinetics following various T cell doses and frequencies of CCR7+ T cells in the population of infused lymphocytes. The model predicted the effects of T cell numbers and of population composition of infused T cells on the resultant concentration of T cells in various organs. For example, a higher percentage of CCR7+ T cells among exogenously administered T lymphocytes led to an augmented accumulation of T cells in the spleen. The model predicted a linear dependence of T cell dynamics on the dose of adoptively transferred T cells.

Discussion: The mathematical model of T cell migration presented here can be integrated into a multi-scale model of the immune system and be used in a preclinical setting for predicting the distribution of genetically modified T lymphocytes in various organs, following adoptive T cell therapies.

1 Introduction

Processes related to the migration and circulation of lymphocytes are critical for the functioning of the immune system and have been extensively studied over the past 60 years (1–3). Initially, T cells are produced in the thymus. Upon positive and negative selection, they reach the systemic circulation (4). A portion of naïve and effector T cells migrate from blood to tissue parenchyma (5); another portion of T cells migrate to lymph nodes (LNs) via the afferent lymphatic vessels, which collect lymph from non-lymphoid organs, although a majority of T cells enter LNs through high endothelial venules (HEVs) (6). This process is controlled by the expression of specific chemokine receptors on the surface of T cells, and only CCR7⁺ T cells may enter LNs through HEVs (7). In lymphoid tissue, T cells may undergo homeostatic proliferation in the absence of antigen or interact with antigen-presenting cells to get activated, proliferate (in the presence of antigen) and/or undergo apoptosis, and then leave LNs through efferent lymphatic vessels (8). Upon leaving the LNs via lymphatic vessels, all T cells flow to the thoracic duct, which returns them into the systemic circulation (9). The above scheme represents the most basic steps of T cell trafficking in the organism; these are further modulated via numerous and complex molecular and cellular interactions.

T cells are the main effector cells in the adaptive immune response and have been increasingly considered for use as therapeutic agents. Currently, such therapeutic interventions are based on the infusion of autologous or genetically modified T cells, such as tumor-infiltrating lymphocytes (TILs), gene-modified T cells expressing novel T cell receptors (TCR-based therapies), chimeric antigen receptor T cells (CAR-T therapies), CAR-NK therapies, and other treatment modalities. The multiple molecular modifications that occur in the immune engineering and manufacturing of T cell therapies significantly affect T cell migration characteristics, once these cells are infused back in vivo, thereby leading to specific T cell pharmacokinetic patterns and subsequent targeted and non-targeted pharmacological responses. A detailed, quantitative understanding of the migration processes of mixed population of exogenously delivered and endogenous T cells is, therefore, of fundamental importance in the development and optimization of translational and clinical strategies for such T cell-based therapies.

A multitude of in vitro and in vivo studies have been conducted, in the investigation of processes that regulate T lymphocyte migration. One foundational experimental analysis has been reported by Smith and Ford (10). In related studies, rats received radiolabeled thoracic duct lymphocytes intravenously, for an evaluation of radioactivity in various organs. Another series of foundational experiments has been performed by Sprent and Miller, who conducted studies in mice treated with labeled T or B cells (11, 12). Subsequently, the utilization of radiolabeled lymphocytes emerged as the predominant method for the in vivo examination of immune cell migration and distribution (13, 14). A number of recent studies have been reported, whereby T cell distribution in different tissues was quantitatively determined through Cr⁵¹-radiolabeled T cell kinetics (15), following adoptive therapy T cell administration (activated or gene-modified T cells) (16, 17).

Given the high degree of complexity of T cell migration processes - including numerous, intricate and nonlinear dynamic interactions - there have been multiple attempts to describe these mathematically. Mohler and Farooqi (18) and later Ganusov and Auerbach (19) described T cell migration processes based on data from experiments of Smith and Ford (10), using physiologically-based compartmental models. Specifically, Mohler and Farooqi (18) developed a family of related mathematical models. One model was formulated with delay differential equations (DDEs) and four other models made use of ordinary differential equations (ODEs) with variations in their structure, i.e., linear and non-linear models, time-variant and time-invariant models. In their model, the authors proposed a compartmental structure whereby blood plasma was considered as a central compartment, while other organs were secondary in terms of lymphocyte migration. Ganusov and Auerbach used their model to estimate lymphocyte residence time in various organs and proposed hypotheses about the mechanisms of lymphocyte arrest in LNs during activation (19). The limitation of the mathematical models calibrated using the data from the Smith and Ford studies relates to the lack of information on absolute numbers of lymphocytes, since the data and corresponding variables in the models were characterized as “percentages of the injected dose”, while the dose of infused lymphocytes was not explicitly provided. Later, Ganusov and Tomura used the same approach to describe in vivo T cell recirculation as a function of exogenously administered T cell kinetics (20). In that model, T cell migration in LNs was described with reference to cell fluorescence data measured in Kaede mice (21). Additional compartments were considered to describe the delay in T cell transition through LNs, and a two-compartment model was finally considered as optimal.

To investigate T cell migration through a specific organ such as the spleen or LNs, a set of experiments was conducted to further develop models of lymphocyte migration (22, 23). Later, Khot et al. (15) and Singh et al. (24) used a multi-compartmental, physiologically-based pharmacokinetic (PBPK) approach to develop mathematical models of T cell recirculation. In such models, each organ was typically divided into extravascular and intravascular spaces. This division allowed one to model the transition of T cells through a given organ, their exit from the blood circulation, and their migration through the lymphatic system. These models described the kinetics of exogenously administered T cells. However, they did not take into account endogenous T cell homeostasis and its impact on trafficking of the overall T cell population.

PBPK and population PK/PD models, combined with model-based simulations have become a quantitative companion tool in the decision-making process of novel therapies research and development, including small molecule and biologic therapies (25, 26). With the high number of cell-based therapies that are emerging, the development of a validated, flexible framework for PBPK modeling is critically needed (27). Consequently, the main objective of the present work was to develop a general and flexible PBPK framework, whereby, for the first time, quantitative aspects of homeostatic turnover and trafficking of resident T cells, under steady-state conditions, are combined with pharmacokinetic aspects of exogenously administered engineered T cells, within a single PBPK model. Additionally, we aimed at assessing the sensitivity of the model to key parameters controlling T cell migration processes.

2 Methods

2.1 Experimental data

We conducted a systematic review of published experimental data on T cell trafficking in mouse. Carefully curated datasets were used for model calibration and validation. In total, data from three preclinical studies that included 30 experimental measurements were used for the estimation of model parameters (i.e., 25 values of mean T-cell concentrations in time for different tissues, and 5 values for steady-state concentrations of endogenous T cells). Three additional studies, providing 24 experimental timepoint measurements were used for independent cross-validation. A multistep approach was implemented for model development, as shown in Figure 1A.

Figure 1

Figure 1 (A) Model development workflow; (B) Schematic of T cell migration model. The model describes the trafficking of two independent populations of T cells - exogenously administered T cells and endogenous T cells - across 6 compartments (blood, lungs, regional LN compartment, generalized LN compartment, spleen, and liver), with additional intra-tissue, time-related ‘delay compartments’ for exogenously delivered T cells located in the spleen, liver, and the regional and generalized LN compartments.

Parameters of endogenous T cell recirculation were estimated based on experimental data describing T cell concentrations under non-inflammatory conditions and under steady-state physiological conditions in mouse tissues (28, 29). The estimated maximal rate of homeostatic T cell proliferation in the spleen was verified against data in mice which underwent thymectomy (30). Parameters of exogenously administered T cell recirculation were estimated based on murine experimental data from Khot et al. (15). In that study, autologous Cr⁵¹-labeled spleen T cells were injected intravenously, and radioactivity levels (per gram of tissue) were measured in various organs. To incorporate these data into the model, we parameterized T cell concentrations using Equation S4 (see Supplementary Material). Area-under-the-curve (AUC) estimates show that exogenously administered cells quickly transfer from blood and lung (AUC=(0.089 ± 0.057)*${10}^8$ and (2.26 ± 0.54)*${10}^8$ cells/ml), to then accumulate in liver and lymph nodes (AUC=(6.75 ± 2.21)*${10}^8$ and (1.62 ± 0.60)*${10}^8$ cells/ml), to ultimately home to the spleen (AUC= (30.65 ± 11.37)*${10}^8$ cells/ml).

Experimental data describing T cell kinetics following an intravenous administration of Zr⁸⁹-labeled T/CAR-T cells in mouse were used for an external cross-validation of the model. In those experiments, radioactivity detection was performed via positron-emission tomography (PET) for various types and doses of T cells as follows: (i) infusion of labeled T cells in a dose of 20.0* ${10}^6$ cells (31); (ii) infusion of CAR-T cells in doses of 1.90-, 7.10-, and 16.8*${10}^6$ cells (16); and (iii) infusion of CAR-T cells in a dose of 1.00*${10}^6$ cells (32).

2.2 Mathematical model

To capture the kinetics and inter-organ recirculation of endogenous and exogenously administered T cells, we developed a physiologically-based compartmental model of T cell migration, formulated with a system of sixteen ODEs. The model scheme is presented in Figure 1B. It includes 6 compartments corresponding to the following organs: blood, lungs, spleen, liver, ‘regional LNs’ and ‘generalized LNs’ (the latter representing secondary lymphoid structures). One of the main assumptions in our model relates to the independent circulation of endogenous vs. exogenous T cells in the various organs. Also, we assumed that exogenously administered and endogenous T cells displayed identical rates of organ-to-blood migration. For other model parameters, we considered quantitative differences in parameter values, given that experimental settings of exogenous cell product manufacturing resulted in variously modified T cells. Firstly, we assumed different rates of migration to the spleen, lung and LNs, for endogenous vs. exogenously administered T cells. Secondly, we assumed that the migratory behavior of exogenously administered T cells would not depend on disease status (i.e., identical T cell distribution patterns in healthy vs. tumor-bearing mice, in non-tumor organs). Thirdly, we assumed that exogenously administered T cells undergo degradation in the lungs [via phagocytosis by alveolar macrophages, based on a model by Zhu et al. (33)] and in blood – thereby accounting for all apoptotic processes in the whole organism.

Importantly, we did not include any proliferation potential of exogenous T cells in the model, since we did not consider the mechanisms required to stimulate antigen-specific proliferation (including exogenous T cell activation and interactions with antigen-presenting cells) in the present study.

Model equations were formulated using the following phenomenological view on T cell systemic and tissue dynamics. Endogenous T cells are produced in the thymus and migrate to blood (model parameter V_{inf_thymus}) and various organs. The rate of homeostatic T cell proliferation in the spleen was described using a Michaelis-Menten like equation. The corresponding EC₅₀ value for homeostatic proliferation in the spleen was estimated as the degree of T cell depletion in the mouse spleen following irradiation, upon which proliferating T cells were observed (34). T cells can migrate from blood to tissues and back. For lung and liver compartments, the model considered T cell trafficking to lymph nodes through afferent lymphatic vessels. The rate of T cell migration from LNs to blood was assumed to be the same as the rate of T cell migration through the thoracic duct. Small distal LNs were represented in the model by a regional LN compartment, with T cells being transported from there to a generalized LN compartment via efferent lymphatic vessels.

Overall, intravenously administered T cells were set to follow the same migration patterns as endogenous T cells (see Figure 1B). However, characteristic residence times were observed for exogenously administered T cells, in specific organs such as the spleen and the liver (10, 35). Various mechanistic explanations have been put forward to explain these observations; e.g., Ganusov and Tomura proposed that a tight network of sinusoids in liver and spleen may slow down cell migration at supra-high T cell numbers (20). An alternative biological explanation is that the rate of T cell migration may well be heterogeneous and dependent on morphological aspects of the spleen, e.g., the white pulp, red pulp, and marginal zones (22). To integrate these biophysical features into the model, we included additional variables representing the exogenous T cell abundance in the spleen, liver and LNs compartments $(T_{ex\_del\_spleen},T_{ex\_del\_liver},T_{ex\_del\_ln},T_{ex\_del\_regln},$ respectively). Exogenously administered T cells were set to migrate out of blood and the lungs according to the following clearance processes: systemic apoptosis in blood and degradation by alveolar macrophages in lungs.

The parameters of the model of the model are listed in Table 1. Overall, the system of equations was formulated as follows:

Table 1

Table 1 Model parameter values.

$\begin{array}{ll}\begin{array}{c}\frac{d{T}_{blood}}{dt}=V_{inf\_thymus}-T_{blood}\ast ((1-f_{ccr7})\ast {\mathit{\mu }}_{12}+{\mathit{\mu }}_{13}+(1-f_{ccr7})\\ \ast {\mathit{\mu }}_{14}+f_{ccr7}\ast {\mathit{\mu }}_{15}+f_{ccr7}\ast {\mathit{\mu }}_{16})+{\mathit{\mu }}_{21}\ast T_{lung}+{\mathit{\mu }}_{41}\\ \ast T_{liver}+{\mathit{\mu }}_{51}\ast T_{ln}+{\mathit{\mu }}_{31}\ast T_{spleen}-k_{apo}\ast T_{blood}\end{array}& (1)\end{array}$

$\begin{array}{ll}\begin{array}{c}\frac{d{T}_{ex\_blood}}{dt}=-T_{ex\_blood}\ast ((1-f_{ccr7\_ex})\ast {\mathit{\mu }}_{12ex}+(1-f_{ccr7\_ex})\ast {\mathit{\mu }}_{14}\\ +{\mathit{\mu }}_{13ex}+f_{ccr7\_ex}\ast {\mathit{\mu }}_{15ex}+f_{ccr7\_ex}\ast {\mathit{\mu }}_{16ex})+{\mathit{\mu }}_{21}\begin{array}{l}\ast T_{ex\_lung}+{\mathit{\mu }}_{41}\ast T_{ex\_liver}+{\mathit{\mu }}_{51}\ast T_{ex\_ln}+{\mathit{\mu }}_{31}\ast T_{ex\_spleen}\\ -k_{apo}\ast T_{ex\_blood}\end{array}\end{array}& (2)\end{array}$

$\begin{array}{ll}\frac{d{T}_{lung}}{dt}=(1-f_{ccr7})\ast {\mathit{\mu }}_{12}\ast T_{blood}-{\mathit{\mu }}_{21}\ast T_{lung}-{\mathit{\mu }}_{25}\ast T_{lung}& (3)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_lung}}{dt}=(1-f_{ccr7\_ex})\ast {\mathit{\mu }}_{12ex}\ast T_{ex\_blood}-{\mathit{\mu }}_{21}\ast T_{ex\_lung}-{\mathit{\mu }}_{25}\ast T_{ex\_lung}-k_{deg}\ast T_{ex\_lung}& (4)\end{array}$

$\begin{array}{ll}\frac{d{T}_{liver}}{dt}=(1-f_{ccr7})\ast {\mathit{\mu }}_{14}\ast T_{blood}-{\mathit{\mu }}_{41}\ast T_{liver}-{\mathit{\mu }}_{45}\ast T_{liver}& (5)\end{array}$

$\begin{array}{ll}\begin{array}{c}\frac{d{T}_{ex\_liver}}{dt}=(1-f_{ccr7\_ex})\ast {\mathit{\mu }}_{14}\ast T_{ex\_blood}-{\mathit{\mu }}_{41}\ast T_{ex\_liver}\hfill \\ \hspace{1em}\hspace{1em}-{\mathit{\mu }}_{45}\ast T_{ex\_liver}+k_{off4}\ast T_{ex\_del\_liver}\hfill \\ \hspace{1em}\hspace{1em}-k_{del4}\ast T_{ex\_liver}\hfill \end{array}& (6)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_del\_liver}}{dt}=k_{del4}\ast T_{ex\_liver} -k_{off4}\ast T_{ex\_del\_liver}& (7)\end{array}$

$\begin{array}{ll}\frac{d{T}_{spleen}}{dt}={\mathit{\mu }}_{13}\ast T_{blood}-{\mathit{\mu }}_{31}\ast T_{spleen}+\frac{Vma{x}_{pro}\ast T_{spleen}}{EC{50}_{pro}+T_{spleen}}& (8)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_spleen}}{dt}={\mathit{\mu }}_{13ex}\ast T_{ex\_blood}-{\mathit{\mu }}_{31}\ast T_{e\text{x}\_spleen}-k_{del3}\ast T_{e\text{x}\_spleen}+k_{off3}\ast T_{ex\_del\_spleen}& (9)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_del\_spleen}}{dt}=k_{del3}\ast T_{ex\_spleen}-k_{off3}\ast T_{ex\_del\_spleen}& (10)\end{array}$

$\begin{array}{ll}\begin{array}{c}\frac{d{T}_{ex\_ln}}{dt}=f_{ccr7\_ex}\ast {\mathit{\mu }}_{15ex}\ast T_{ex\_blood}+{\mathit{\mu }}_{45}\ast T_{ex\_liver}\\ \hspace{1em}\hspace{1em}+{\mathit{\mu }}_{25}\ast T_{ex\_lung}-{\mathit{\mu }}_{51}\ast T_{ex\_ln}\\ \hspace{1em}\hspace{1em}-k_{del5}\ast T_{ex\_ln}+k_{off5}\ast T_{ex\_del\_ln}+{\mathit{\mu }}_{65}\ast T_{ex\_regln}\end{array}& (11)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ln}}{dt}=f_{ccr7}\ast {\mathit{\mu }}_{15}\ast T_{blood}+{\mathit{\mu }}_{45}\ast T_{liver}+{\mathit{\mu }}_{25}\ast T_{lung}-{\mathit{\mu }}_{51}\ast T_{ln}+{\mathit{\mu }}_{65}\ast T_{regln}& (12)\end{array}$

$\begin{array}{ll}\frac{d{T}_{regln}}{dt}=f_{ccr7}\ast {\mathit{\mu }}_{16}\ast T_{blood}-{\mathit{\mu }}_{65}\ast T_{regln}& (13)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_regln}}{dt}=f_{ccr7\_ex}\ast {\mathit{\mu }}_{16ex}\ast T_{ex\_blood}-{\mathit{\mu }}_{65}\ast T_{ex\_regln}-k_{del6}\ast T_{ex\_regln}+k_{off6}\ast T_{ex\_del\_regln}& (14)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_del\_ln}}{dt}=k_{del5}\ast T_{ex\_ln}-k_{off5}\ast T_{ex\_del\_ln}& (15)\end{array}$

$\begin{array}{ll}\frac{d{T}_{ex\_del\_regln}}{dt}=k_{del6}\ast T_{ex\_regln}-k_{off6}\ast T_{ex\_del\_regln}& (16)\end{array}$

The first equation describes the influx of endogenous T cells from the thymus to the blood compartment (first term: $V_{inf\_thymus}$). T cell migration from blood ($T_{blood})$ to all other compartments (lungs, spleen, liver, generalized LNs and regional LNs with rate constants of, respectively, µ₁₂, µ₁₃, µ₁₄, µ₁₅, µ₁₆) is described by the second term. The selective transition of CCR7⁺ T cells to LNs is taken into account by multiplying the rate of T cell migration from blood to LN compartments (regional and generalized) by f_ccr7, and the transfer rate to non-lymphoid organs, where non-naive T cells migrate to, by (1-f_ccr7). The next five terms in Equation 1 describe T cell migration from, respectively, lungs, liver, spleen, generalized LNs and regional LNs, to blood. The last term of Equation 1 represents the rate of T cell death via apoptosis.

Equation 2 describes the rate of change of exogenously administered T cell abundance in blood, T_{ex_blood}, following the infusion. This rate is determined by the balance of the rates of cell transitions from blood to lungs, liver, spleen, generalized LNs and regional LNs, and back to blood. The last term in Equation 2 depicts the apoptotic rate of exogenously administered T cells.

Equation 3 describes the rate of change of T cell levels in the lungs. The three terms in Equation 3 represent the kinetics of T cell transition from blood to lungs, from lungs to blood, and from lungs to the generalized LN compartment.

According to Equation 4, exogenously administered T cells migrate to the lungs, and back, via the same route as endogenous T cells. They degrade in the lungs, according to the last term in Equation 4. This T cell degradation feature was taken from the model by Zhu et al. (33).

Equation 5 is similar to Equation 3; it describes the rates of T cell migration from blood to liver and back, and T cell flow within afferent lymphatics to LNs.

Equation 6 represents the migration of exogenously administered T cells from blood to liver and the migration from liver to blood and LNs. The last two terms describe the transitions of exogenously administered T cells through a ‘delay-causing’ compartment within the liver. The rates of T cell inflow and outflow via the delay compartment in Equation 7 are assumed to be linear.

The spleen is a lymphoid organ; we considered T cell homeostatic proliferation, as described by the last term of Equation 8.

Equation 9 depicts processes which affect exogenously administered T cell numbers in the spleen: migration from blood to spleen, and back. To reproduce T cell kinetics in the spleen, another delay compartment was introduced. The last two terms of Equation 9 represent rates of exogenously administered T cell transitions to the delay compartment and back.

Equation 10 describes the transition of exogenously administered T cells within the spleen to the delay compartment and back.

The inflow of CCR7⁺ T cells from blood to LNs through HEVs was taken into account in the first terms of Equations 11, 14 (exogenously administered T cells) and of Equations 12, 13 (endogenous T cells). The first three terms in Equations 11, 12 are similar; they represent the inflow of T cells from blood and the transition of T cells from non-lymphoid organs (liver and lung, respectively). The last terms in these two equations describe the rate of T cell inflow via efferent lymphatics to the generalized LNs. Equation 11 also includes delay-related terms, to capture cell accumulation in the LNs, as shown in previous research (19, 20). The first two terms in Equations 13, 14 represent rates of T cell inflow from blood - through HEVs - to regional LNs and migration via efferent lymphatics to the generalized LNs compartment. For exogenously administered T cells in the regional LN compartment, we also considered cell accumulation via a delay compartment, as described by the last two terms of Equation 14.

Equations 15, 16 represent the dynamics of exogenously administered T cells in the delay compartments for, respectively, the generalized and regional LNs.

Overall, the model includes 30 parameters which are summarized in Table 1. Some parameters were assumed to differ between endogenous vs. exogenously administered T cells (e.g., rate constants for migration from blood to spleen and generalized and regional LNs). Parameters which are in common were sub-divided into three groups: a) those verified against data describing steady-state levels of CD3⁺ T cells in mouse tissues (28, 29); b) those estimated based on half-lives of T cells in various organs (19); and c) those estimated using additional experimental data. The last group included parameters for endogenous T cell homeostatic proliferation, which were estimated using data in irradiated mice, and parameters for exogenously administered T cells, which were adjusted against experimental data on exogenously administered T cell kinetics in different tissues (15).

2.3 Model calibration

A nonlinear fixed-effect modeling approach was used for the calibration of model parameters with a log-normal distribution assumed for the experimental data. The following log-likelihood cost (LLC) function was maximized to find an optimal set for the fitted parameter values:

$$L{L}_y(\widehat{\theta })=log(L_y(\widehat{\theta }))\widehat{=}\log (p(y;\widehat{\theta }))$$

where $\widehat{\theta }$ is the vector of parameter estimates for the model being considered, y represents the observed data, and $p(y;\widehat{\theta })$ is the probability density function of the observed data given the parameter estimates. The model with the highest value of $L{L}_y(\widehat{\theta })$ was selected as the model of choice.

In addition, relative standard errors (RSE) for model parameters $\widehat{\theta }$ and the residual error model were estimated. For all parameters, RSE values lower than 50% were a prerequisite for model selection.

Weighted residuals (WRES) were calculated as the normalized difference between observations and their expected values:

$$WRES=V^{-1/2}(y-E(\text{y}))$$

where $V$ is the variance-covariance matrix of $y$, and E(y) is the mathematical expectation of the model solutions obtained using parameters $\widehat{\theta }$ perturbed by a random value sampled from the residual error model.

2.4 Software

The Monolix software, version 2020R1, was used for fixed-effects modeling, while data visualisation and forward simulations were performed in the R software, version 4.0.2 (R-project, www.r-project.org). Model sensitivity analyses were performed using the “sensitivity” R package; model simulations were performed using the “RsSimulx” R package; visualization of model solutions was performed using “ggplot2” and “tidyverse” R packages.

3 Results

3.1 Model calibration

We used a sequential approach to curate and integrate all available quantitative information into the model (see the flowchart in Figure 1A).

In a first stage, we calibrated parameters pertaining to endogenous T cell migration (V_{inf_thymus}, µ₁₂, µ₁₃, µ₁₄, µ₁₅, µ₂₅, µ₄₅, µ₁₆), based on data describing the steady-state levels of T cells in various compartments. In addition, we determined a ratio between rates of T cell migration through HEVs (V_{blood_ln}) vs. afferent lymphatics (V_{lung_ln}, V_{liver_ln}), based on the knowledge that only ∼10% of T cells migrate to lymph nodes via afferent lymphatics (6). The rate constants for T cell migration from various tissues to blood (µ₂₁, µ₃₁, µ₄₁, µ₆₁) were calculated using estimates of cellular half-lives provided by Ganusov and Tomura (20). The parameter value for µ₅₁was calculated as the rate of T cell flux to the thoracic duct. As shown in Figure 2F, the model consistently described the experimental data, i.e., endogenous T cells would accumulate in lymphoid organs at concentrations of 1.57*${10}^8$ cells/ml in the spleen and 3.36*${10}^8$ cells/ml in LNs, whereas T cell levels in blood and lungs at steady-state were significantly lower (1.70*10⁶ cells/ml and 4.22*${10}^6$ cells/ml, respectively); and liver exhibited the lowest level of T cell infiltration (∼4.20*${10}^5$ cells/ml).

Figure 2

Figure 2 (A–E) Model calibration against data on exogenously administered T cells in: (A) blood, (B) lungs, (C) spleen, (D) regional LN compartment, and (E) liver. Points: Experimental data (mean ± error), calculated from (15). Curves: Model predictions; yellow uncertainty prediction bands: 95% confidence intervals (CI). (F) Model calibration against data on endogenous T cells. Dark blue columns: Model predictions; yellow columns: Experimental data, calculated from (28, 29, 36).

In a second stage of model calibration, we fixed all parameters for endogenous T cells, all rate constants for exogenous T cell migration from organs to blood (assumed to be the same as for the endogenous cell migration), as well as the parameter that characterizes the rate of T cell apoptosis. We then estimated parameters µ_13ex, µ_12ex, µ_15ex, µ_16ex and parameters describing the delay in T cell migration into liver, spleen and LNs (k_del3, k_del4, k_del5, k_del6) using experimental data on exogenously administered T cell kinetics (15). Since parameters k_del3, k_del5, k_del6 represent the rate constants for the accumulation of cells in lymphoid organs, we assumed the same value for these three parameters. The resultant calibrated model consistently described available experimental data (Figures 2A–E). In particular, the model adequately reproduced the fast T cell elimination from blood (Figure 2A) and the comparably rapid transitory kinetics in the highly vascularized lungs (Figure 2B), as well as the slower kinetics in the spleen (with the peak value C_max of ∼12.6*10⁶ cells/ml at 9 hrs) and the liver (C_max of ∼3.4*10⁶ cells/ml at 8 hrs) (Figure 2E). The model also satisfactorily described the cellular kinetics in the regional LN compartment, with a C_max of 0.5*10⁶ cells/ml at 45.5 hrs post administration (Figure 2D). Goodness-of-fit plots (Figure 3) provide further support for an adequate description of the experimental data by the model. In particular, the moderate over-estimation of C_max in blood holds for only some of the reported blood measurements (squares on Figure 3A), with weighted residuals (WRES) lying within a [−2, +2] range (Figure 3B).

Figure 3

Figure 3 (A) Goodness-of-fit for exogenously administered T cell concentrations in spleen (square with a cross inscribed), liver (circles), lungs (asterisks), blood (squares) and regional LN compartment (triangles). (B) Weighted residuals for exogenously administered T cells in all compartments.

It is worth noting that, initially, we had considered a model structure featuring only one lumped LN compartment, i.e., without any regional LN compartment (Supplementary Figure 1, Supplementary Material). The corresponding model could not describe the available experimental LN data for exogenously administered T cells, yet it adequately captured the cellular dynamics in tissues other than LNs (Supplementary Figure 2, Supplementary Material).

3.2 Model validation

To validate the model and assess its predictive power, we performed simulations of external data which had not been used in the model calibration process. In particular, the model was validated against experimental data in mice which received infusions of radiolabeled T cells or CAR-T cells, in doses of 1.0*10⁶, 1.9*10⁶, 7.1*10⁶, or 16.8*10⁶ cells per animal (16, 31, 32). Due to the different methods used for T cell labeling, with corresponding differing detections and units of radioactivity in a region of interest (ROI) in the experimental studies, we adjusted the data using a scaling coefficient, ω, in Equation 17 as follows:

$\begin{array}{ll}C(oncentration)ofTcellsinorgan(scaled)\left\{cells/ml\right\}=\omega *CofTcellsinorgan\left\{cells/ml\right\}& (17)\end{array}$

An averaged calibration coefficient value of 2.7 was determined for ω, although there are several reasons for significant variability across studies due to the characteristics of the various experiments (See Section 5, Supplementary Material).

As shown in Figure 4A, the model adequately described experimental data for doses of 1.9*10⁶ and 7.1*10⁶ CAR-T cells. However, it slightly over-predicted peak levels following doses of 1.0*10⁶ and 1.9*10⁶ (C_max values of 0.6*10⁶ and 1.1*10⁶ cells/ml predicted by the model vs. observed C_max values of 0.3*10⁶ and 0.5*10⁶ cells/ml, respectively). At the same time, T cell levels at later timepoints, following a dose of 16.8*10⁶ cells were under-estimated, with a predicted value C_max of 1.02*10⁷ cells/ml vs. a higher C_max of 1.50*10⁷ cells/ml observed in the experimental data. These differences between observed and model-predicted values may have been caused by uncertainties related to the calibration of the radioactivity labeling data. Since the actual calibration curves for the various experimental settings were not available, we applied an empirical approach to map the data onto a uniform scale, which may have introduced an additional bias. To further validate the model, we also used external data on exogenous T cell levels in the spleen, under various doses (13). The linear dependence between T cell counts in the spleen and the administered doses was successfully reproduced by the model, for this external set of experimental data (Figure 4B).

Figure 4

Figure 4 (A) Model validation against data on exogenously administered, radiolabeled T/CAR-T cells in mouse spleen. Solid curves: model predictions; bands: 95% confidence intervals; points: experimental data with confidence intervals taken from publications of Sta Maria et al. and Wang et al. (16, 31, 32). Brown shadows and datapoint circles correspond to the administration of CAR-T cells (10⁶ cell dose); violet shadows and datapoint squares: administration of CAR-T cells (1.9*10⁶ cell dose); green shadows and datapoint triangles: administration of CAR-T cells (7.1*10⁶ cell dose); dark blue shadows and datapoint rhombi: administration of T cells (12.0*10⁶ cell dose); yellow shadows and datapoint asterisks: administration of CAR-T cells (16.8*10⁶ cell dose). (B) Dependence of exogenously administered T cell numbers in the spleen vs. given dose of labeled T cells. Black line: model simulations; green line: experimental data taken from (13).

3.3 Sensitivity analysis

Two types of sensitivity analyses were conducted: a global sensitivity analysis for parameters and variables describing endogenous T cell homeostasis, and a local sensitivity analysis for those describing exogenously administered T cell numbers and accumulation.

The Partial Rank Correlation Coefficient (PRCC) method was used to perform the global sensitivity assessment of the model (47). Model parameters were sampled using a Latin hypercube sampling (LHS) procedure. 5000 sets of parameters were sampled for single estimation of PRCC. The procedure was repeated 200 times to estimate standard errors. The lower and upper bounds of the parameters were assumed to be determined by values corresponding to a two-fold change from steady-state T cell concentrations. Model solutions were computed over a large time window (up to 80,000 hrs), to ensure steady-state conditions. Under such a setting, the PRCC may characterize the relative importance or effect of a parameter, along with its positive/negative correlation mode, vs. model output variability.

As shown in Figure 5, the PRCC coefficients for the rate constants of T cell migration from blood to spleen and back (µ₁₃ and µ₃₁) were about 4-times higher in terms of their impact on steady-state levels of T cells in spleen vs. blood. A similar relationship was found for the lungs (see Supplementary Figure 3 in the Supplementary Material): the rate constants of migration from blood to the lungs and back (µ₁₂ and µ₂₁) were more important for the lungs than for blood, with corresponding PRCC coefficients about 10-times higher. According to this sensitivity analysis, the key systemic parameters with comparable impact on T cell homeostasis in all compartments were V_{inf_thymus} and k_{_apo}. Parameters of homeostatic proliferation Vmax_pro and EC50_pro were also important for T cell levels in all compartments.

Figure 5

Figure 5 Global sensitivity analysis for steady-state concentrations of endogenous T cells in: (A) blood; and (B) spleen. Bars depict PRCCs (partial rank correlation coefficients) for each parameter (mean ± standard error).

Although the sensitivity analysis predicted rather low PRCC coefficients for T cell inflow rate parameters from the lungs and liver to the LNs via afferent lymphatic vessels (µ₂₅ and µ₄₅), in relation to T cell levels in the blood and the spleen, it revealed that these processes were essential in driving T cell levels in the organs themselves. This is demonstrated by the 5x higher values of the corresponding PRCC coefficients for the lungs and liver vs. other compartments (excluding LNs for µ₂₅), which highlights the importance of considering the model’s sensitivity to the lymphatic system parameters. The T cell inflow constant from lungs to LNs, µ₂₅, carried significance for homeostasis in the lung and LNs, as shown in Supplementary Figures 3A, C of the Supplementary Material. Also, it should be noted that the parameter $f_{ccr7}$was predicted to play a minor influence on T cell homeostasis in blood and spleen, while the levels of T cells in liver, the lungs, and both regional and generalized LN compartments appeared to be more sensitive to $f_{ccr7}$.

To estimate the impact of kinetic parameters on exogenously administered T cell dynamics, a local sensitivity analysis was conducted. Parameter values were varied in the following way: decreased by 50% or increased by 100% from their best-fit estimates (i.e., two-fold changes). Figure 6 shows the sensitivity profiles for exogenous T cell in blood (A), spleen (B), liver (C) and regional LN (D). Supplementary Figure 6 in the Supplementary Material presents the dependence of the total number of exogenous T cells vs. time, which was calculated as the sum of T cell numbers in all compartments. For all compartments excluding blood, the most important parameters were the rate constants of T cell inflow from blood to organ: changes in µ_13ex, µ₁₄ and µ_16ex resulted in proportional changes of the AUC in spleen, liver, and regional LNs, respectively. The same behavior was shown to hold for the lungs (see Supplementary Figure 4A in the Supplementary Material).

Figure 6

Figure 6 Local sensitivity analysis for exogenous T cell steady-state concentrations in: (A) blood; (B) spleen; (C) liver; and (D) regional LN compartment.

For the generalized LN compartment (see Figure S4B), the most important parameter was the T cell inflow rate from organ to blood (µ₅₁). Interestingly, the rate constant of T cell transition to the lungs, µ_12ex, was significant with respect to all variables, including the level of T cells in blood. This rate relates to the change in the fraction of total exogenous T cells trafficking to the lungs, where a relatively high level of T cell degradation is set in the model (Supplementary Figure 6B), red and grey curves, Supplementary Material). Also, an increase in the rate constant, k_del3, which describes the T cell transition rate to the delay compartment in the spleen, resulted in a stronger accumulation of exogenous T lymphocytes in the delay compartment and a higher AUC in the spleen (Figure 6B). The same sensitivity pattern was observed for other delay compartments: increases in k_del4, k_del5 and, to a lesser degree, in k_del6 resulted in higher T-cell AUC values in the liver as well as in the generalized LN and regional LN compartments. On the contrary, increases in k_off parameters resulted in a faster migration to, and exit from the corresponding compartments and, as a consequence, in lower AUC values (as shown for k_off3, k_off4, k_off6 in Figures 6B–D), respectively). The rate constant capturing the inflow of cells to the regional LN compartment and the delay-related parameter for that compartment (µ_16ex, k_del6) were important only within this compartment; inflow rate constants and delay parameters in all other compartments had a low impact on T cell AUC values. We also explored the importance of these parameters with respect to their effects on C_max values (Supplementary Figure 5, Supplementary Material) and obtained similar results. Except for parameters µ_12ex, µ_13ex, µ_14, µ_16ex, which were the most important parameters determining C_max in the corresponding compartments, the inflow rates from blood to other compartments were more significant than the rate constants reflecting transitions to delay compartments. For example, for the value of C_max in the spleen, the rate constants of T cell inflow to the lungs and liver (µ_12ex, µ₁₄) were more significant than k_del3, k_off3.

Overall, these kinetic behaviors reflect the inter-dependence of T cell levels in various organ and tissue compartments, which are inherently linked in a continuous flow-like system such as the whole organism.

3.4 Model simulations

Model simulations were conducted to investigate the effects of specific factors such as exogenous T cell dose and CCR7 expression levels, on T cell kinetics and distribution. To explore more specifically the dose dependency on cellular kinetics, we performed the following numerical experiments representing intravenously administered T cells at doses of 1.0, 3.0, 5.0, 10, 20, 50, or 100 million cells per mouse (Figure 7). Low doses of 1.0, 3.0, 5.0 or 10 million T cells resulted in C_max values in blood of, respectively: {0.46, 1.40, 2.30, 4.65}*10⁶ cells/ml. High doses of 20, 50, or 100 million T cells resulted in C_max values of, respectively, {9.30, 23.3, 46.5}*10⁶ cells/ml. Consequently, there is a linear relationship between the dose of adoptively transferred cells and T cell concentrations in the spleen (see Figure 4B, 72 hr post-infusion timepoint), which is observed experimentally and predicted by the calibrated model. A similar relationship was revealed in other compartments (see Figures 7B–F). We also simulated T cell concentration in the generalized LN compartment (Figure 7F) and found a linear dependence between the dose and the accumulation of lymphocytes in this compartment.

Figure 7

Figure 7 Model simulations of exogenously administered T cell recirculation kinetics upon given cell doses of 1.00*10⁶ to 100*10⁶ cells in: (A) blood; (B) lungs; (C) spleen; (D) liver; (E) generalized LN compartment; and (F) regional LN compartment.

In order to investigate distribution effects of CCR7 expression levels on exogenously administered T cells, the parameter f_{ccr7_ex} was varied over a range of 0.1 to 0.9 (Figure 8). As summarized in Table 2, a change in f_{ccr7_ex} from 0.1 to 0.3 did not significantly affect the exogenous cell AUC in spleen and blood (respectively, {9.91,11.09,12.59}*10⁸; and {0.11, 0.12, 0.14}*10⁸ cells/ml). A further increase in CCR7 expression levels (with increasing f_{ccr7_ex} values of 0.7, 0.8, 0.9) led to a higher accumulation of T cells in spleen and blood, with AUC levels of {27.52, 39.1, 67.38}*10⁸ in the spleen and {0.26, 0.36,0.61}*10⁸ in blood (Figure 8C). Interestingly, this non-linear relationship is not observed in the liver and lung compartments, where the effects of varying CCR7 expression levels were less pronounced. Indeed, there were only minor differences in kinetics, whether low (0.1 to 0.5) or high (0.7 to 0.9) values of f_{ccr7_ex} were tested (Figures 8B, D). The C_max value in non-lymphoid organs, reflective of the distribution phase, was inversely related to f_{ccr7_ex}: it increased from 2.63*10⁶ up to 4.01*10⁶ cells/ml for the liver, and from 7.67*10⁶ to 31.3*10⁶ cells/ml for the lungs, with decreasing values of f_{ccr7_ex} (from 0.9 to 0.1). However, during the elimination phase, higher C_max values led to a slightly lower trough concentrations in the lungs at 360 hrs, with trough cell levels decreasing from 7.7*10⁴ to 7.5*10⁴ cells/ml. Therefore, the lowest value of C_max (f_{ccr7_ex} = 0.9) led to the highest trough level of T cells in the lungs, although the difference was rather modest.

Figure 8

Figure 8 Model simulations of exogenous T cell recirculation kinetics following an administration of T cells in doses of 20 million cells with different CCR7⁺ T cell proportions (0.1 - 0.9), (A) blood; (B) lungs; (C) spleen; (D) liver; (E) generalized LN compartment; and regional LN compartment (F).

Table 2

Table 2 Area-under-the-curve (AUC) values and biodistribution coefficients (BC) as a function of CCR7⁺ cell proportions (f_{ccr7_ex}) in exogenously administered T cells.

This observed non-linearity in cell kinetics also led to significant differences in the biodistribution coefficient, BC (Table 2). An increase in f_{ccr7_ex} resulted in an increase in BC in the regional LN compartment, yet in a reduced T cell concentration relatively to blood in the generalized LN compartment (with a lower BC value in LNs). BC in non-lymphoid organs was also diminished; an increase in $f_{ccr7\_ex}$ resulted in a higher BC in the spleen compartment. Interestingly, the AUC value appeared to change faster in blood vs. LNs.

To illustrate key differences in the kinetic behavior of exogenous vs. endogenous T cell populations, we performed simulations which simultaneously tracked the kinetics of these two T cell populations. Concentrations of endogenous T cells were at steady-state in each compartment (see Figure 2F and Supplementary Figure 7, dark red lines, Supplementary Material). Since exogenous T cells were administered in high numbers, their kinetics appeared as temporary ‘perturbations’ in T cell numbers in various compartments, with higher or lower exposures depending on the presence of delay compartments (see Figures 2A–E) and Supplementary Figure 7, dark red lines, Supplementary Material). However, even the largest dose of exogenously administered T cells would not exceed ∼1% of the total lymphocyte count in mouse; exogenous T cells would therefore not be expected to significantly perturb the dynamics of endogenous T cells.

4 Discussion

Various approaches have been used to study T cell migration and distribution processes in the body: experimental methods include the use of exogenously administered T cells (labeled, pre-activated, with subsequent antigen injection, etc.), whereas theoretical approaches include the development and analysis of mathematical models describing T cell migration throughout the whole organism and within individual organs. The works by Mohler and Farooqi (18), Ganusov and Auerbach (19) and others have yielded significant insights into the kinetics of T lymphocyte migration and distribution. However, given the growing interest in both fundamental immunology and practical applications of preclinical and clinical research, there is a pressing need for an analytical tool that can be used to quantitatively predict proportions of T cells which may enter a target organ in response to, e.g., adoptive T cell therapy or immune activation. Therefore, the main objective of our work was to characterize, via modeling, the migratory behavior and inter-compartmental trafficking of T cells, under homeostatic conditions and upon intravenous infusion.

We thus developed a physiologically-based compartmental model describing the dynamics of exogenously administered T cells and endogenous T cells. The model was used to quantitatively capture homeostatic levels of T cells in mouse tissues and the recirculation kinetics of exogenously administered T cells following intravenous infusion. The model includes time-dependent variables representing CD4⁺ and CD8⁺ T lymphocytes as a single population of CD3⁺ T cells, even though it is known that CD4⁺ and CD8⁺ T cells differ in their rates of migration, apoptosis, and proliferation (48). It is established that the expression of CCR7 is a key determinant for the migration of T lymphocytes to the LNs via HEVs (6). We considered this feature by including parameters f_ccr7 and f_{ccr7_ex} into the model, which represent the fractions of CCR7⁺ T cells in the populations of endogenous and exogenously administered T lymphocytes, respectively. In the model, the entry of T cells into the regional and generalized LN compartments was CCR7-dependent. We neglected any potential effects of CCR7 expression on the inflow to the spleen, in agreement with experimental data showing that mainly CCR7^- T cells were detectable in the spleen (49). However, this receptor is significant for further migration through the spleen to zones (e.g., white pulp) where active immune responses take place and may, therefore, be important when modeling T cell activation processes in the spleen.

Model parameters were verified based on data from exogenously administered T cell studies.

In an initial version of the model, only one generalized LN compartment was included. The model adequately reproduced T cell kinetics in all compartments with the exception of the LNs, where the model solution predicted higher levels of T cells. The main reason for this discrepancy lied in the structural difference between our model-based description of LNs, as a lumped compartment, vs. distinct LNs with specific features underlying the experimental data. In that initial model, the lumped LN compartment would collect T cells in the lymph from a number of organs (spleen and lungs), with a high rate of efflux into blood (as calculated from the rate of T cell flux in the thoracic duct). In contrast, the experimental data described T cell kinetics in distinct, specific lymph nodes, which collect T cells from the mouse’s hind leg. With the addition of a regional LN compartment representative of these smaller, distal LNs to the system of model compartments, we were able to consistently describe the experimental data. At the same time, the generalized LN compartment provided an opportunity to estimate the inflow of T cells from the lymphatic system to blood more adequately.

In order to assess the predictive power of the model, we used experimental data on Cr⁵¹-labeled T cell distribution (15, 16, 31), since Cr⁵¹ labeling is a broadly used measurement technique in T cell distribution studies. More recent studies have also been conducted, with an estimation of T cell radioactivity using PET measurements in the live mouse (16, 31, 32). In the absence of calibration curves across the various radioactivity measurement experiments, we had to introduce an empirical scaling factor, ω, into the model, to map the data from these heterogeneous studies onto a common scale, although this may have introduced an inaccuracy in the dose vs. T cell number relationship. Our mathematical model adequately described the T cell recirculation kinetics following low and intermediate doses of administered CAR-T/T cells. However, for the highest administered dose, the model performed well only during an early time window following cell infusion (first 48 hours). There may be two reasons to explain such discrepancy: (i) a non-linear relationship between the dose (absolute number of T cells) and radioactivity may be in effect, as has been shown using another radioactive label (17); and (ii) altered intrinsic properties of CAR-T cells engineered ex vivo. In the Sta Maria et al. study, mice were infused with CAR-T cells, whereby the CAR included a 41-BB co-receptor (16). There is indirect evidence indicating that persistence of CAR-T with 41-BB may differ vs. other T/CAR-T cells (50).

We also conducted a global sensitivity analysis to rank the mathematical model parameters with respect to their impact on changes in endogenous T cell dynamics in tissues. We determined that the rates of influx and efflux from tissues were more important for T cell homeostasis in organs than in blood. Systemic parameters affected homeostatic T cell levels in blood, lymphoid and non-lymphoid tissues. It is known that only ∼10% of T cells reach LNs via afferent lymph vessels. The mathematical model shows that this flow is a more important determinant of homeostatic levels of T cells in organs from which lymphocytes were released vs. the steady-state level of T cells in LNs (6). The parameter f_ccr7, which captures the proportion of CCR7⁺ endogenous T cells, had a higher influence on T cell levels in non-lymphoid organs and in regional LNs vs. in the spleen, generalized LNs, or blood. The local sensitivity analysis we performed for exogenous T cells also showed the importance of the rates of influx and efflux from tissues. Importantly, the rate of inflow to the lungs also had a strong effect on exogenous T cell accumulation. In addition, the maximal level of T cells (C_max value) depended on rates of influx into other organs, reflective of a continuous network-type and flow-like system.

Using the model, we simulated T cell kinetics and distribution following an intravenous administration of various doses of T cells engineered ex vivo. The model predicted a linear dose-dependence in blood and spleen, in full accordance with results reported by Zatz et al. (13). We conducted further simulations of engineered T cell infusions with differing fractions of CCR7⁺ T cells. Increasing the fraction of CCR7⁺ T cells in the infused dose contributed to an increase in T cell numbers in the spleen. This is supported by indirect evidence from the literature: when CAR-T cells were preincubated with IL-7/IL-15, CCR7⁺ (TCM, TSCM) cell counts in the infusate subsequently increased (51), and IL-7/IL-15-preincubated CAR-T cells were present in higher numbers in the spleen, 3 days post-infusion (52). While CCR7 expression may be associated with cytotoxic activity, cytokine production, and proliferation potential, its effect on migratory features of T cells has already been exploited in CAR-T development.

Our modeling study also presents some limitations. Firstly, we assumed that trafficking of endogenous vs. exogenously administered T cell populations would occur independently of one another. These two subsets of T cells, however, share similar cytokine/chemokine receptors and could, for example, compete for corresponding ligand molecules or interact with each other while migrating through the same vessels and tissues. Also, these two T cell populations interact during in vivo expansion and functional cytotoxic responses. However, including such exogenous vs. endogenous T cell interaction processes would require direct experimental data. Additional experimental data would be needed for the quantitative introduction and explicit inclusion of various T cell sub-populations (e.g., CD4⁺ and CD8⁺ T cells) into the model. Such an extension would be warranted to take into account processes related to antigen stimulation, whereby distinct T cell sub-populations would need to be described separately (proliferating/non-proliferating, naïve, effector, memory, resident T cells). Also, we did not consider the proliferation of exogenous T cells in this model, since mechanisms of exogenous T cell activation and interactions with antigen-presenting cells were not included in the present study; experimental data (Cr⁵¹-radiolabeled T cells) used for model calibration did not include quantitative information on newly developing T cells. Despite such limitations, the quantitative multi-compartmental PBPK model presented here holds the potential for integration into a broader multiscale immune system model. It may be further used to address questions underlying the rational development of immune cell-based therapies and anti-infection therapies of an immunological nature.

The model may also be applied, preclinically, to perform predictive simulations of genetically modified T lymphocyte levels in various organs, following the administration of adoptive T/CAR-T cells. In recent years, various gene modifications of autologous T cells (‘armored’ T cells) have been developed and tested in preclinical and clinical studies (53). The effects of such modifications, e.g., on T cell persistence, may be captured using the model structure presented here, based on changes in the parameter $k_{apo}$. Gene modifications resulting in changes in cell migration (e.g., CCL19 expression by CAR-T cells) could be investigated by adding a tumor compartment in the model and assuming an ‘LN-like’ influx rate into such a compartment, with CCL19-producing CAR-T cells and all other CCR7⁺ T cells (exogenous and endogenous) migrating, to some degree, to the same locations.

One further application of this model relates to its translation to human studies, e.g., to address the kinetics of immune cells distribution throughout the human body to complement a recent systematic analysis presented in (36). Using human data on T cell homeostasis and dynamics of exogenously administered T cells, we may scale parameters from the present model, calibrated for mouse, to derive a corresponding model for human. Such a model may then allow for predictive simulations of cell kinetics following adoptive T/CAR-T cell therapies, an effective tool for dose optimization purposes with such therapeutic approaches.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

Author contributions

AN: Investigation, Writing – original draft, Visualization. GH: Validation, Writing – review & editing. KP: Conceptualization, Supervision, Validation, Writing – review & editing. GB: Conceptualization, Supervision, Validation, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Academic leadership program Priority 2030 proposed by the Federal State Autonomous Educational Institution of Higher Education, the I.M. Sechenov First Moscow State Medical University of the Ministry of Health of the Russian Federation (Sechenov University), the Ministry of Science and Higher Education of the Russian Federation (Agreement 075-10-2021-093, Project MMD-RND-2266). GB was supported by the Moscow Center of Fundamental and Applied Mathematics at INM RAS (agreement with the Ministry of Education and Science of the Russian Federation No. 075-15-2022-286). AN and KP were supported by the Russian Science Foundation (Grant Number 23-71-10051). The authors declare that this study received funding from Modeling & Simulation Decisions FZ - LLC. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Acknowledgments

We would like to acknowledge the contribution of all our colleagues at M&S Decisions and the MID3 Research Center at Sechenov University for their continuous support.

Conflict of interest

Author KP is an employee and shareholder of Modeling & Simulation Decisions, FZ - LLC. Author GH is an employee of Biorchestra US Inc., Cambridge, Massachusetts, USA.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1357706/full#supplementary-material

References

1. Gowans JL. The recirculation of lymphocytes from blood to lymph in the rat. J Physiol. (1959) 146:54–69. doi: 10.1113/jphysiol.1959.sp006177

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bhan AK, Reinisch CL, Levey RH, McCluskey RT, Schlossman SF. T-cell migration into allografts. J Exp Med. (1975) 141:1210–5. doi: 10.1084/jem.141.5.1210

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Chauveau A, Pirgova G, Cheng HW, De Martin A, Zhou FY, Wideman S, et al. Visualization of T cell migration in the spleen reveals a network of perivascular pathways that guide entry into T zones. Immunity. (2020) 52:794–807. doi: 10.1016/j.immuni.2020.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Janeway CA Jr, Travers P, Walport M, Shlomchik MJ, Travers P. Immunobiology. 5th ed. New York: Garland Science (2001).

Google Scholar

5. Cose S. T-cell migration: a naive paradigm? Immunology. (2007) 120:1–7. doi: 10.1111/j.1365-2567.2006.02511.x

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Hunter MC, Teijeira A, Halin C. T cell trafficking through lymphatic vessels. Front Immunol. (2016) 7:613. doi: 10.3389/fimmu.2016.00613

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Vella G, Guelfi S, Bergers G. High endothelial venules: A vascular perspective on tertiary lymphoid structures in cancer. Front Immunol. (2021) 12:736670. doi: 10.3389/fimmu.2021.736670

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Delves P. J., Martin S. J., Burton D. R., Roitt I. M. Roitt’s Essential Immunology. 13th Edition. Chichester, West Sussex; Hoboken, [NJ]: John Wiley & Sons, Inc. (2017).

Google Scholar

9. Hampton HR, Chtanova T. Lymphatic migration of immune cells. Front Immunol. (2019) 10:1168. doi: 10.3389/fimmu.2019.01168

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Smith ME, Ford WL. The recirculating lymphocyte pool of the rat: a systematic description of the migratory behaviour of recirculating lymphocytes. Immunology. (1983) 49:83–94.

PubMed Abstract | Google Scholar

11. Sprent J. Circulating T. and B lymphocytes of the mouse. I. Migratory properties. Cell Immunol. (1973) 7:10–39. doi: 10.1016/0008-8749(73)90180-9

CrossRef Full Text | Google Scholar

12. Sprent J. Fate of H2-activated T lymphocytes in syngeneic hosts. I. Fate in lymphoid tissues and intestines traced with 3H-thymidine, 125I-deoxyuridine and 51chromium. Cell Immunol. (1976) 21:278–302. doi: 10.1016/0008-8749(76)90057-5

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zatz MM, Lance EM. The distribution of 51Cr-labeled lymphocytes into antigen-stimulated mice. J Exp Med. (1971) 134:224–41. doi: 10.1084/jem.134.1.224

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Austyn JM, Kupiec-Weglinski JW, Morris PJ. Migration patterns of dendritic cells in the mouse. Adv Exp Med Biol. (1988) 237:583–9. doi: 10.1007/978-1-4684-5535-9_89

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Khot A, Matsueda S, Thomas VA, Koya RC, Shah DK. Measurement and quantitative characterization of whole-body pharmacokinetics of exogenously administered T cells in mice. J Pharmacol Exp Ther. (2019) 368:503–13. doi: 10.1124/jpet.118.252858

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Sta Maria NS, Khawli LA, Pachipulusu V, Lin SW, Zheng L, Cohrs D, et al. Spatio-temporal biodistribution of 89Zr-oxine labeled huLym-1-A-BB3z-CAR T-cells by PET imaging in a preclinical tumor model. Sci Rep. (2021) 11:15077. doi: 10.1038/s41598-021-94490-0

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Volpe A, Lang C, Lim L, Man F, Kurtys E, Ashmore-Harris C, et al. Spatiotemporal PET imaging reveals differences in CAR-T tumor retention in triple-negative breast cancer models. Mol Ther. (2020) 28:2271–85. doi: 10.1016/j.ymthe.2020.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Farooqi ZH, Mohler RR. Distribution models of recirculating lymphocytes. IEEE Trans BioMed Eng. (1989) 36:355–62. doi: 10.1109/10.19856

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ganusov VV, Auerbach J. Mathematical modeling reveals kinetics of lymphocyte recirculation in the whole organism. PLoS Comput Biol. (2014) 10:e1003586. doi: 10.1371/journal.pcbi.1003586

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Ganusov VV, Tomura M. Experimental and mathematical approaches to quantify recirculation kinetics of lymphocytes. In: Molina-París C, Lythe G, editors. Mathematical, Computational and Experimental T Cell Immunology. Cham: Springer International Publishing (2021). p. 151–69. doi: 10.1007/978-3-030-57204-4_10

CrossRef Full Text | Google Scholar

21. Tomura M, Yoshida N, Tanaka J, Karasawa S, Miwa Y, Miyawaki A, et al. Monitoring cellular movement in vivo with photoconvertible fluorescence protein “Kaede” transgenic mice. Proc Natl Acad Sci U S A. (2008) 105:10871–6. doi: 10.1073/pnas.0802278105

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Hammond BJ. A compartmental analysis of circulatory lymphocytes in the spleen. Cell Tissue Kinet. (1975) 8:153–69. doi: 10.1111/j.1365-2184.1975.tb01217.x

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Regoes RR, Barber DL, Ahmed R, Antia R. Estimation of the rate of killing by cytotoxic T lymphocytes in vivo. Proc Natl Acad Sci U S A. (2007) ;104:1599–603. doi: 10.1073/pnas.0508830104

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Singh AP, Zheng X, Lin-Schmidt X, Chen W, Carpenter TJ, Zong A, et al. Development of a quantitative relationship between CAR-affinity, antigen abundance, tumor cell depletion and CAR-T cell expansion using a multiscale systems PK-PD model. MAbs. (2020) 12:1688616. doi: 10.1080/19420862.2019.1688616

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Green D. Physiologically Based Pharmacokinetic Analyses — Format and Content Guidance for Industry. Available online at: https://www.fda.gov/media/101469/download.

Google Scholar

26. Madabushi R, Seo P, Zhao L, Tegenge M, Zhu H. Review: role of model-informed drug development approaches in the lifecycle of drug development and regulatory decision-making. Pharm Res. (2022) 39:1669–80. doi: 10.1007/s11095-022-03288-w

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chaudhury A, Zhu X, Chu L, Goliaei A, June CH, Kearns JD, et al. Chimeric antigen receptor T cell therapies: A review of cellular kinetic-pharmacodynamic modeling approaches. J Clin Pharmacol. (2020) 60 Suppl 1:S147–59. doi: 10.1002/jcph.1691

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Boyer SW, Rajendiran S, Beaudin AE, Smith-Berdan S, Muthuswamy PK, Perez-Cunningham J, et al. Clonal and quantitative in vivo assessment of hematopoietic stem cell differentiation reveals strong erythroid potential of multipotent cells. Stem Cell Rep. (2019) 12:801–15. doi: 10.1016/j.stemcr.2019.02.007

CrossRef Full Text | Google Scholar

29. Misumi I, Mitchell JE, Lund MM, Cullen JM, Lemon SM, Whitmire JK. T cells protect against hepatitis A virus infection and limit infection-induced liver injury. J Hepatol. (2021) 75:1323–34. doi: 10.1016/j.jhep.2021.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

30. den Braber I, Mugwagwa T, Vrisekoop N, Westera L, Mögling R, de Boer AB, et al. Maintenance of peripheral naive T cells is sustained by thymus output in mice but not humans. Immunity. (2012) 36:288–97. doi: 10.1016/j.immuni.2012.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Maria NS, Khawli L, Lin S, Pachipulusu V, Zheng L, Cohrs D, et al. Simultaneous PET/MRI measurement of primary T cells biodistribution in naïve mice following adoptive cell transfer. J Nucl Med. (2019) 60:277–7.

Google Scholar

32. Wang XY, Wang Y, Wu Q, Liu JJ, Liu Y, Pan DH, et al. Feasibility study of 68Ga-labeled CAR T cells for in vivo tracking using micro-positron emission tomography imaging. Acta Pharmacol Sin. (2021) 42:824–31. doi: 10.1038/s41401-020-00511-5

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Zhu H, Melder RJ, Baxter LT, Jain RK. Physiologically based kinetic model of effector cell biodistribution in mammals: implications for adoptive immunotherapy. Cancer Res. (1996) 56:3771–81.

PubMed Abstract | Google Scholar

34. Dummer W, Niethammer AG, Baccala R, Lawson BR, Wagner N, Reisfeld RA, et al. T cell homeostatic proliferation elicits effective antitumor autoimmunity. J Clin Invest. (2002) 110:185–92. doi: 10.1172/JCI0215175

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Sprent J, Miller JF. Fate of H2-activated T lymphocytes in syngeneic hosts. II. Residence in recirculating lymphocyte pool and capacity to migrate to allografts. Cell Immunol. (1976) 21:303–13. doi: 10.1016/0008-8749(76)90058-7

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Cosgrove J, Hustin LSP, de Boer RJ, Perié L. Hematopoiesis in numbers. Trends Immunol. (2021) 42:1100–12. doi: 10.1016/j.it.2021.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Wehrmann F, Lavelle JC, Collins CB, Tinega AN, Thurman JM, Burnham EL, et al. γδ T cells protect against LPS-induced lung injury. J Leukoc Biol. (2016) 99:373–86. doi: 10.1189/jlb.4A0115-017RR

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ma Q. Severe pneumonia induces immunosenescence of T cells in the lung of mice. Aging (Albany NY). (2023) 15:7084–97. doi: 10.18632/aging.v15i14

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Shrewsbury MM. Rate of flow and cell count of thoracic duct lymph in the mouse. Proc Soc Exp Biol Med. (1958) 99:53–4. doi: 10.3181/00379727-99-24244

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Deaton JG. Thoracic duct lymph drainage in the mouse: a technique for producing lymphocyte-depleted animals. Lymphology. (1972) 5:115–20.

PubMed Abstract | Google Scholar

41. Ionac M, Laskay T, Labahn D, Geisslinger G, Solbach W. Improved technique for cannulation of the murine thoracic duct: a valuable tool for the dissection of immune responses. J Immunol Methods. (1997) 202:35–40. doi: 10.1016/S0022-1759(96)00226-8

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Shrewsbury MM. Thoracic duct lymph in unanesthetized mouse: method of collection, rate of flow and cell content. Proc Soc Exp Biol Med. (1959) 101:492–4. doi: 10.3181/00379727-101-24992

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Vianna PHO, Canto FB, Nogueira JS, Nunes CFCG, Bonomo AC, Fucs R. Critical influence of the thymus on peripheral T cell homeostasis. Immun Inflammation Dis. (2016) 4:474–86. doi: 10.1002/iid3.132

CrossRef Full Text | Google Scholar

44. Bjorkdahl O, Barber KA, Brett SJ, Daly MG, Plumpton C, Elshourbagy NA, et al. Characterization of CC-chemokine receptor 7 expression on murine T cells in lymphoid tissues. Immunology. (2003) 110:170–9. doi: 10.1046/j.1365-2567.2003.01727.x

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Shah DK, Betts AM. Towards a platform PBPK model to characterize the plasma and tissue disposition of monoclonal antibodies in preclinical species and human. J Pharmacokinet Pharmacodyn. (2012) 39:67–86. doi: 10.1007/s10928-011-9232-2

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Jiao Y, Zhang J, Yan J, Stuart J, Gibson G, Lu L, et al. Differential gene expression between wild-type and Gulo-deficient mice supplied with vitamin C. Genet Mol Biol. (2011) 34:386–95. doi: 10.1590/S1415-47572011005000031

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol. (2008) 254:178–96. doi: 10.1016/j.jtbi.2008.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

48. De Boer RJ, Homann D, Perelson AS. Different dynamics of CD4+ and CD8+ T cell responses during and after acute lymphocytic choriomeningitis virus infection. J Immunol. (2003) 171:3928–35. doi: 10.4049/jimmunol.171.8.3928

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Sharma N, Benechet AP, Lefrançois L, Khanna KM. CD8 T Cells Enter the Splenic T Cell Zones Independently of CCR7, but the Subsequent Expansion and Trafficking Patterns of Effector T Cells after Infection Are Dysregulated in the Absence of CCR7 Migratory Cues. J Immunol. (2015) 195:5227–36. doi: 10.4049/jimmunol.1500993

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Guedan S, Posey AD, Shaw C, Wing A, Da T, Patel PR, et al. Enhancing CAR T cell persistence through ICOS and 4-1BB costimulation. JCI Insight. (2018) 3:e96976, 96976. doi: 10.1172/jci.insight.96976

CrossRef Full Text | Google Scholar

51. Xu XJ, Song DG, Poussin M, Ye Q, Sharma P, Rodríguez-García A, et al. Multiparameter comparative analysis reveals differential impacts of various cytokines on CART cell phenotype and function ex vivo and in vivo. Oncotarget. (2016) 7:82354–68. doi: 10.18632/oncotarget.v7i50

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Xu Y, Zhang M, Ramos CA, Durett A, Liu E, Dakhova O, et al. Closely related T-memory stem cells correlate with in vivo expansion of CAR.CD19-T cells and are preserved by IL-7 and IL-15. Blood. (2014) 123:3750–9. doi: 10.1182/blood-2014-01-552174

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Foeng J, Comerford I, McColl SR. Harnessing the chemokine system to home CAR-T cells into solid tumors. Cell Rep Med. (2022) 3:100543. doi: 10.1016/j.xcrm.2022.100543

PubMed Abstract | CrossRef Full Text | Google Scholar