Joint models quantify associations between immune cell kinetics and allo-immunological events after allogeneic stem cell transplantation and subsequent donor lymphocyte infusion

Alloreactive donor-derived T-cells play a pivotal role in alloimmune responses after allogeneic hematopoietic stem cell transplantation (alloSCT); both in the relapse-preventing Graft-versus-Leukemia (GvL) effect and the potentially lethal complication Graft-versus-Host-Disease (GvHD). The balance between GvL and GvHD can be shifted by removing T-cells via T-cell depletion (TCD) to reduce the risk of GvHD, and by introducing additional donor T-cells (donor lymphocyte infusions [DLI]) to boost the GvL effect. However, the association between T-cell kinetics and the occurrence of allo-immunological events has not been clearly demonstrated yet. Therefore, we investigated the complex associations between the T-cell kinetics and alloimmune responses in a cohort of 166 acute leukemia patients receiving alemtuzumab-based TCD alloSCT. Of these patients, 62 with an anticipated high risk of relapse were scheduled to receive a prophylactic DLI at 3 months after transplant. In this setting, we applied joint modelling which allowed us to better capture the complex interplay between DLI, T-cell kinetics, GvHD and relapse than traditional statistical methods. We demonstrate that DLI can induce detectable T-cell expansion, leading to an increase in total, CD4+ and CD8+ T-cell counts starting at 3 months after alloSCT. CD4+ T-cells showed the strongest association with the development of alloimmune responses: higher CD4 counts increased the risk of GvHD (hazard ratio 2.44, 95% confidence interval 1.45-4.12) and decreased the risk of relapse (hazard ratio 0.65, 95% confidence interval 0.45-0.92). Similar models showed that natural killer cells recovered rapidly after alloSCT and were associated with a lower risk of relapse (HR 0.62, 95%-CI 0.41-0.93). The results of this study advocate the use of joint models to further study immune cell kinetics in different settings.


JOINT MODEL I
Joint models only consider measurements taken prior to the occurrence of the clinical events of interest. Occasionally, the measurement time and event time coincide: for example, T-cell counts may be recorded on the same day as the start of therapeutic systemic immunosuppression for Graft-versus-Host-Disease (GvHD). In order to retain the information of the measurements taken at event times, we set the time of these measurements to one day earlier, which assumes that the measurement at the event time was representative of the T-cell counts the day before the event. However, we excluded measurements at time of relapse, since the presence of blasts in the peripheral blood could lead to incorrect counts of the normal T-cells. We also excluded measurements at time of autologous recovery, as donor-derived T-cells were no longer present, and therefore also no potentially alloreactive T-cells capable of inducing GvHD or Graft-versus-leukemia (GvL) effect.

Model formulation
The longitudinal submodel assumes that the true underlying (log) immune cell counts (either CD3, CD4, CD8, or NK) for the i th patient are given by with random effects vector b i ∼ N (0, D). The observations for the i th patient at timepoints t ij (j = 1, . . . , n i ) are given by Risk i , Donor i and CMV i respectively represent the dummy variables for baseline disease risk (the intention-to-treat variable, high-risk compared to non-high risk), donor type (unrelated compared to related donor) and patient/donor Cytomegalovirus (CMV) serostatus at baseline (any one of patient or donor positive, compared to patient and donor both negative).
Time since allogeneic stem cell transplantation (alloSCT) was modelled flexibly assuming restricted (natural) cubic splines with two internal knots placed at the 33.3% and 66.7% percentiles of the measurement times. This is represented above by B q (t), corresponding to the q th basis function of the spline. The fixed effects part of the model posits a three-way interaction between time, donor type and baseline disease risk, as well as a main effect of patient/donor CMV status. The three-way interaction was constructed to a) capture the slower expected average trajectory of patients with an unrelated donor, due to the use of anti-thymocyte globulin (ATG) in this group; and b) to test for a difference in average trajectories between baseline disease risk groups.
In terms of random effects, this models assumes random slopes b iq (one for each basis function), and a fixed intercept. This fixed intercept was justified given that this cohort underwent T-cell depleted (TCD) alloSCT, and all patients were therefore expected to start follow-up with immune cell counts close to zero. The random slopes were assumed to be normally distributed with mean zero, with unstructured covariance matrix D.
The time-to-event submodel was composed of multiple cause-specific proportional hazards models as where the h ki (t) for k ∈ {1, 2, 3} respectively represent the cause-specific hazards of GvHD, relapse, and other failures. The cause-specific baseline hazards h k0 (t) were approximated on the log scale using cubic B-splines with three internal knots. The above corresponds to the 'current value' parametrization of the joint model, where the exp(α k ) would represent the hazard ratio (for cause k) when comparing two patients (with same covariates) whose 'true' (model-based) underlying log immune cell values at a particular timepoint m i (t) differ by one. The γ kp coefficients are interpreted analogously to main effects in standard cause-specific Cox proportional hazards models.
In addition to the current value parametrization, we also ran the models assuming a time-dependent slopes association structure as α k1 m i (t) + α k2 {dmi(t)/dt}.

Standardized Residuals
Above we present standardized residuals plots, which summarize how well the model fits the data overall (i.e. across all observations) -both for the average and subject-specific trajectories. The fitted (i.e. log immune cell counts predicted by the model) values are plotted against the standardized distance between the observed measurement and the predicted value. The blue line is a smoothed average of the standardized residuals as a function of the fitted values, and should ideally be horizontal at 0.

Model formulation
For model II, the time scale was no longer from alloSCT, but instead from time of early low-dose donor lymphocyte infusion (DLI). Therefore, this model was only run among the subset that did in fact receive an early low-dose DLI before the occurrence of other competing events. Furthermore, some patients did not have a T-cell measurement on the day of DLI but only a few days prior. For these patients, we used the measurement closest to DLI taken within the last week before DLI as the measurement at time of DLI (time 0). The longitudinal submodel was again a linear mixed-effects model, where the true underlying log T-cell counts are given by N (0, D). Observations for i th patient are again given by where ϵ ij ∼ N (0, σ 2 ) are independent random error terms.
Time was again modelled with restricted cubic splines, but in contrast to model I, we used a single internal knot. The focus on a shorter timespan resulted in a reduced sample size, and fewer measurements per person. For consistency with model I, this average trajectory was allowed to differ across donor types (two-way interaction). In this model, disease risk at baseline was redundant as we ran the model among those having actually received an early low-dose DLI. A fixed effect for patient/donor CMV serostatus was also added to the model. This model comprised both random intercepts b i0 and random slopes b iq , assumed to follow normal distributions with mean zero and unstructured covariance matrix.
Due to a limited number of events, relapse and other failures were merged into a composite endpoint. The time-to-event submodel was therefore specified as where the h ki (t) for k ∈ {1, 2} respectively represent the cause-specific hazards of GvHD and the composite of relapse and other failures for subject i. The cause-specific baseline hazards h k0 (t) were approximated on the log scale using cubic B-splines with two internal knots. In this joint model, only the current value parametrization was explored.