We expanded our previously developed model describing the dynamics of murine heart transplant rejection (30) to examine the impact of adoptive transfer of Tregs on graft survival. The original model consists of 13 nonlinear first order ordinary differential equations (ODEs) tracking the following populations in a representative lymph node compartment (LN) and graft compartment (G): CD4 T cells ( T H L N and T H G ) , CD8 T cells ( T E L N and T E G ) , regulatory T cells ( T R L N and T R G ) , immature DCs ( A i m m ) , mature DCs ( A m a t L N and A m a t G ) , inflammatory macrophages ( A i n f ), pro- and anti-inflammatory cytokines ( C p and C a ) , and graft cells ( G ). An additional ODE is introduced to track naïve Tregs ( T R N L N ). All model equations and parameters appear in the Supplemental Digital Content, and all model variables and their initial values are provided in Supplementary Table S1. A complete list of model assumptions is provided in (30). While this model is parametrized for mouse heart transplantation, the predicted interactions between the immune response and transplanted organ would be similar for other transplants.

In this study, the model equations tracking the immature and mature DCs are adapted to include more realistic representations of experimental observations (36). Specifically, the decay rate of immature DCs is assumed to depend on the remaining graft mass (Supplementary Eq. S10, second term). In addition, DC maturation is assumed to occur in the presence of pro-inflammatory cytokines or CD4 T cells (Supplementary Eq. S10, third term; Supplementary Eq. S11, first term). For verification of the modified model and updated predictions for host immune dynamics without adoptive transfer treatment, see the Supplemental Digital Content (Supplementary Figures S2, S3).

A dosing function, D ( t ) , for the adoptive transfer of Tregs is defined in Eq. 1 (also in Supplementary Eq. S12). Parameters f G , f L N , and f N correspond to the fraction of the Treg dose that enters the graft as activated Tregs (Eq. 2), the lymph node as activated Tregs (Eq. 3), and the lymph node as naïve Tregs (Eq. 4), respectively. For simulations involving adoptive transfer, f G + f L N + f N = 1 . These three parameters are set to zero when adoptive transfer is not simulated. The dosing function is composed of one or multiple exponentially decaying functions of time, where D ₀ indicates the dosing rate, n indicates the total number of doses delivered, t i indicates the post-operative day (POD) of delivery of the ith Treg injection, and β indicates the decay rate. A decay rate of β = 2  day⁻¹ is used to simulate the relatively fast absorption of Tregs during adoptive transfer. We assume around 50% of injected cells are lost or distribute in non-draining lymphoid tissues (which are excluded from the model), and we define the dose magnitude C to be the total number of injected cells that localize to the modeled lymph node and graft compartments. Thus, for the majority of this study, we simulate adoptive transfer with C   =   5   × 10 5 Tregs. The dose magnitude is the area under the curve of our dosing function, which is given by C =   ∫ 0 ∞ ( D 0 ∑ i = 1 n d i ( t ) )   d t =   n D 0 β cells. A single dose with D 0 =   10 6 cells/day or n doses with D 0   =   10 6 n   cells/day both administer a total of C   =   5   × 10 5 Tregs. Figure 1 shows the dosing rate over time for a single dose administered on POD0 with D 0 = 10 6   cells/day (Figure 1A) and 5 doses administered on POD0, 1, 2, 3, and 4 with D 0 = 10 6 5 cells/day (Figure 1B). The dose magnitude C is varied during model analysis. D ( t ) =   D 0 Σ i = 1 n d ( t )   ,   w h e r e     d i ( t ) = { 0   ,     t < t i   e − β ( t − t i )   ,     t ≥ t i (1) d T R G d t = k e R T R L N − μ R T R G + r R G T R G ( T E G + T H G ) α 6 + T R G + f G D ( t ) (2) d T R L N d t = a R T R N L N A m a t L N γ 2 + A m a t L N − μ R T R L N + r R T R L N ( T E L N + T H L N ) α 2 + T R L N − e R T R L N + f L N D ( t ) (3) d T R N L N d t = μ R N   ( T 0   −   T R N ) + f N D ( t ) (4)

The activation state, accumulation site, magnitude, timing, and frequency of adoptive transfer are varied in the model to simulate the impact of each of these factors on the immune response to the graft. Parameters f G , f L N , and f N are varied to determine the impact of activation state (i.e., naïve or activated) and accumulation site of Tregs. The effects of dose magnitude are assessed by varying C . The timing and frequency of Treg doses are evaluated by varying the start day of the dose ( t 0 ) and the number of doses, n. If multiple doses are administered, it is assumed that the doses are given at 1-day intervals. Following our previous work, graft rejection is defined as a 75% reduction in the original number of graft cells (30). The time at which the model predicts graft rejection is used in this study to identify optimal dosing strategies for the adoptive transfer of Tregs.

To calibrate the model parameters in the updated model (see parameters shaded in Supplementary Table S2), we employed a non-dominated sorting-based multi-objective evolutionary algorithm (37), called NSGA-II (Non-dominated Sorting Genetic Algorithm II), which can greatly improve the efficiency in constrained multi-objective optimization tasks. NSGA-II is a popular non-domination based genetic algorithm for multi-objective optimization and parameter estimation. NSGA-II improves elitism and there is no need to choose sharing parameters a priori.

To analyze the effect of Treg adoptive transfer using our model, we first evaluated the impact of the activation status of the injected cells on graft survival. It is well recognized that pre-activated Treg (even after resting) are an order of magnitude more suppressive than naïve Tregs (38, 39). For this analysis, we simulated the equivalent of delivering C   =   5   × 10 5 Tregs on POD0. As expected, accumulation of naïve Tregs into the lymph node was not as effective as the adoptive transfer of pre-activated Tregs (Supplemental Digital Content, Supplementary Figure S4). This confirmed the benefit of using ex-vivo expanded Tregs, which is a necessary step in most clinical preparations of antigen-specific cells that renders homogeneous populations of pre-activated Treg. For the remainder of the study, we simulated the use of pre-activated Tregs.

Next, we focused on the question: does the tissue distribution of Tregs between the lymphoid compartment and the graft post-adoptive transfer impact transplant survival? Figure 2A depicts the change in the graft mass over time for the scenarios where the transferred Tregs exclusively accumulate in the lymphoid versus graft compartments, as well as for the case of no adoptive transfer (black curve). There is a substantial benefit with the accumulation of Tregs in the graft in comparison to lymphoid tissues, with estimated rejection on POD74 vs. POD25, respectively. The model can simulate a wide range of Tregs dose magnitudes; for almost all ranges, accumulation of Tregs directly to the graft (Figure 2B, blue curve) is more effective at extending survival than localization to the lymph node compartment (Figure 2B, red curve). For very small doses ( C < 2.9 × 10 4 cells), delivering activated Tregs to the lymph node is predicted to be more effective than delivering to the graft (Figure 2C); though, it improves graft survival by no more than 2.4 days.

Having shown the benefit of graft localization, we studied the impact of the timing of a single injection of varying amounts of Tregs. Figure 3 shows a non-monotonic relationship between rejection time and day of Treg administration, observed for some Treg dose magnitudes (Figure 3A). Specifically, when using C   =   5   × 10 5 cells (blue curve), we observed that delaying Treg administration to POD1.5 yields optimal graft survival time. The non-monotonic behavior between dose timing and graft survival is highlighted for physiological dose magnitudes ( C < 10 6 cells) in Figure 3B. Once the dose magnitude is increased above C   =   7.5   × 10 5 cells, a monotonic relationship is re-established in which graft survival time decreases with the delay of Treg administration.

Although a theoretical model can assume that all Tregs administered to an individual accumulate exclusively in the graft or lymphoid tissues, physiological in vivo constraints and Treg properties dictate a variable partitioning between the two compartments. Thus, in Figure 4, we show the simulated impact of varying the fraction of Tregs that enter the graft (f _g ) between 0 and 1; the remaining fraction (1–f _g ) is assumed to enter the draining lymphoid tissue compartment. As indicated by the red curve, the model predicts that if more than 10% of the injected Tregs locate to the graft, delaying injection until POD1.5 is the most beneficial strategy to prolong the time to graft rejection.

The model allows us to compare the protective effect of a single injection of   C   =   5   × 10 5 Tregs versus splitting that total number of cells among multiple daily injections (an approach that can also represent sustaining the presence of a smaller amount of Treg over time). As depicted in Figure 5, distributing a fixed dose prolongs graft survival. The relationship between the number of consecutive doses and graft survival was not monotonic, with graft survival peaking around 10–15 equivalent doses, depending on the first day of injection (Figure 5A). Significantly, when enough doses are administered, delaying the start of adoptive transfer treatment is no longer beneficial. The model predicts that administering C   =   5   × 10 5 cells divided evenly among 14 daily doses beginning on POD0 is the most effective treatment, extending graft survival until POD115. As shown in Figure 5B (graft mass over time), it is noteworthy that in addition to extending transplant survival, multiple doses of Treg provide better protection of the graft (i.e., slower rate of mass decline). Figure 5C shows that the impact of multiple doses compared with a single dose is apparent even if only a small fraction of Tregs (>5%) enters the graft.

Differently from bioinformatic approaches (e.g., machine learning), mathematical modeling enables the analysis of the dynamics of host immune cells that determine the simulated outcome. This valuable property can be used to gain insight into why certain adoptive transfer treatments are more effective than others. We investigated in detail the difference in effects between using a single injection of graft-infiltrating Tregs versus splitting the total number over multiple injections. Figure 6 depicts the predicted dynamics of host immune cells in the graft with single dose adoptive transfer (blue) and 14 dose adoptive transfer (red). For a similar investigation into the impact of site of accumulation and timing of single injections, please see the Supplemental Digital Content (Supplementary Figures S5, S6).

Delivering a large single dose of Tregs minimizes the maturation of graft infiltrating DCs (Figure 6B) and therefore limits the activation and graft accumulation of cytotoxic CD8 T cells (Figure 6C) and delays the activation of macrophages (Figure 6D). Therefore, delivering a substantial portion of Tregs soon after transplantation is effective at limiting the early inflammatory response. Nevertheless, the regulation of the destructive capacity of those cells is not sustained with a single administration.

When the injected cells are split among 14 doses, elevated Treg levels are maintained during the surge of cytotoxic T cells and inflammatory macrophages (Figures 6A,C,D, red). In contrast, for single dose adoptive transfer, Tregs in the graft decay before the surge of these graft destructive subsets (Figures 6A,C,D, blue). The lack of Treg-mediated control in the last scenario allows each cytotoxic T cell and inflammatory macrophage to cause greater graft damage despite a lower accumulation level. This phenomenon is highlighted in the boxed regions in Figures 6C–F. We note that Figures 6E,F show the value of each term on the right-hand side of the differential equation for graft cells that contributes to graft destruction (Supplementary Eq. S9). Therefore, higher values of the curve in panel E or F indicate a greater contribution to graft destruction. For example, from POD7.2 to POD11.9, although more cytotoxic T cells are predicted to be present with the 14-dose regimen (Figure 6C; the 14-dose curve is above the 1-dose curve in the dotted boxed region), their effector functions are regulated and less graft destruction occurs (Figure 6E; the 1-dose curve is above the 14-dose curve). Similarly, from POD5.3 to POD7.7, although more inflammatory macrophages are predicted with the 14-dose regimen (Figure 6D), less graft destruction is mediated by this subset (Figure 6F).

Overall, analyzing the population dynamics for various Treg dose frequencies indicates that, to optimize adoptive transfer, cell delivery must 1) maintain elevated Treg levels during the surge of cytotoxic T cells and inflammatory macrophages in the graft and 2) deliver enough Tregs early after transplantation to inhibit DC maturation and therefore limit T cell activation. These same principles hold true when examining dose timing and can explain why administering a single dose on POD1.5 is more effective than delivering a single dose on POD0 (see Supplemental Digital Content, Supplementary Figure S6). When examining site of Treg accumulation, minimizing DC maturation is again critical to understand the benefit of delivering Tregs to the graft rather than the lymph node (see Supplemental Digital Content, Supplementary Figure S5).

Our model predicts that Tregs accumulating to the transplant are more therapeutically effective than Tregs distributing to the lymph node (Figure 2). Importantly, analysis of immune dynamics provides the rationale behind such a different outcome: inhibition of DC maturation in the graft is more effective than direct inhibition of T cell activation in the draining lymphoid tissue at minimizing the number of CD8 and CD4 T cells that reach the graft (Supplemental Digital Content, Supplementary Figure S5). Interestingly, experimental reports in rodent models have suggested that an immediate accumulation of Treg in the graft promotes longer survival (44, 45). Overall, our analysis suggests that identifying methods to increase the fraction of Tregs that translocate early on to the graft will maximize the impact of adoptive transfer. It is noteworthy that the conditions of ex-vivo Treg expansion impart specific migratory capacity to the cells and these differences impact the therapeutic result obtained (45–49). Moreover, the breadth of genetic engineering that has become possible for clinical products (e.g., CRISPR/Cas9 based modifications) could provide useful tools to impart the ideal behavior in adoptively transferred Tregs.

When administering a single dose of Tregs, our simulations suggest that delaying cell delivery is more effective than administering cells on the day of transplantation (Figure 3B). The justification for this unexpected result is similar to that given for splitting Treg into multiple doses. Slightly delaying administration allows Tregs to directly inhibit the destructive functions of cytotoxic T cells and inflammatory macrophages while ensuring a timely reduction of DC activation and, therefore, T cell activation (Supplemental Digital Content, Supplementary Figure S5). Obviously, this scenario is specific to the unique conditions of administering Tregs without any additional manipulation of the recipient (like immunosuppression, see below), but it highlights the impact of examining cellular interactions at the system level.

The model predicts that maintaining a prolonged influx of Tregs is the most effective method of lengthening graft survival. Delivering multiple, smaller doses of Tregs preserves the graft longer than delivering the same total number of Tregs in a single bolus. This regimen allows for elevated levels of Tregs to remain in the graft during the surge of graft-destructive cells (and control them), while also compromising on an “early enough” limitation of DC maturation. Although the clinical implementation of multiple daily dosing of Treg is unrealistic, this result highlights the important point of sustaining the survival and function of transferred Tregs. This is a contentious issue with the reported negative effect that many immunosuppressive drugs have on the homeostasis and function of Treg (50, 51). There is growing interest in devising approaches to sustain the persistence and function of Tregs. The use of IL-2/anti-IL2 complexes, IL-2 muteins, as well as the genetic engineering of Tregs to respond to “orthogonal” IL-2 are all examples of active investigations to promote Tregs persistence (10, 20, 52). In parallel, the utilization of biomaterials to promote the sustained accumulation of Tregs in proximity to the transplant represent a very promising strategy (53). Our model results provide the rationale to strongly support these ongoing efforts.

While the processes of sensitivity analysis and parameter estimation performed in this study (based on experimental data) have improved the accuracy of several model parameters, some parameters remain uncertain. For example, the persistence and proliferation of Tregs in the graft are unquantified variables that have a profound impact on the dynamics of graft infiltrating immune cells. Similarly, the relationship between number of Tregs injected and the number of cells that reach the graft is not quantified; the factors that influence such a relationship are poorly defined, and thus, further experimentation is needed (54).

While ours is the only ODE transplantation model to date that tracks Tregs independently of other T cells, several immunological mathematical models have recently been developed to elucidate the role of Tregs on self-tolerance and to identify key Treg interactions with other immune populations (55–58). The assumptions from these models may be used to improve our existing model to better replicate Treg behavior. In particular, given the importance of IL-2 to Treg survival, proliferation, and function, explicitly tracking IL-2 concentrations as in (57) would strengthen our current model and allow further investigation into IL-2 therapies as a method of extending Treg survival.

Although there is very limited quantitative information, other pre-clinical transplant models show important concordance with some of our model simulations. In a mouse model of pancreatic islet transplantation (45), Zhang et al. compared the i.v. infusion of Treg with the co-transplantation at the site of grafting. The significantly longer survival obtained in the latter case agrees with our model-suggested principle of a higher therapeutic impact when Treg can rapidly and directly modulate graft immune populations. Unfortunately, their report did not present different doses or timing of Treg administration. A similar scenario emerges from reports using Treg in mouse skin transplant models. The use of CAR-Treg (59), cell lines of Treg (60), or polyclonal Treg (61) injected on the day of the transplant promotes only a modest increase in transplant survival, a result that correlates with simulation of a fixed number of Treg that probably mostly home to the lymphoid compartment. In all cases, the combination with so called “adjunct therapies” (ranging from thymectomy and T cell depletion to irradiation and bone marrow co-transplantation) is demonstrated as necessary to achieve lasting therapeutic effects. Currently, we can only draw qualitative comparisons to these experimental models since theoretical model parameters would need to be adapted to each different transplant model.

Overall, our modified model is beneficial in identifying methods to maximize the benefits of Treg adoptive transfer and in generating hypotheses on the key immune dynamics that govern its outcome and that can be tested experimentally. However, as demonstrated in this study and by experimental evidence to date, Treg adoptive transfer alone is insufficient to prevent transplant rejection. These results highlight the need to understand what additional perturbations to the system would better support or enhance the protective function of Tregs. Our mathematical model provides the framework into which treatments like immunosuppression (existing or hypothetical) can be included and used to dissect their very complex effects. Combined with the promising technological advances in both the investigation and manipulation of cells, there is tangible optimism toward the ultimate goal of optimizing therapeutic strategies for transplantation.