A dynamic time‐to‐event model for prediction of acute graft‐versus‐host disease in patients after allogeneic hematopoietic stem cell transplantation

Abstract Background Acute graft‐versus‐host disease (aGvHD) is a major cause of death for patients following allogeneic hematopoietic stem cell transplantation (HSCT). Effective management of moderate to severe aGvHD remains challenging despite recent advances in HSCT, emphasizing the importance of prophylaxis and risk factor identification. Methods In this study, we analyzed data from 1479 adults who underwent HSCT between 2005 and 2017 to investigate the effects of aGvHD prophylaxis and time‐dependent risk factors on the development of grades II–IV aGvHD within 100 days post‐HSCT. Results Using a dynamic longitudinal time‐to‐event model, we observed a non‐monotonic baseline hazard overtime with a low hazard during the first few days and a maximum hazard at day 17, described by Bateman function with a mean transit time of approximately 11 days. Multivariable analysis revealed significant time‐dependent effects of white blood cell counts and cyclosporine A exposure as well as static effects of female donors for male recipients, patients with matched related donors, conditioning regimen consisting of fludarabine plus total body irradiation, and patient age in recipients of grafts from related donors on the risk to develop grades II–IV aGvHD. Additionally, we found that higher cumulative hazard on day 7 after allo‐HSCT are associated with an increased incidence of grades II–IV aGvHD within 100 days indicating that an individual assessment of the cumulative hazard on day 7 could potentially serve as valuable predictor for later grades II–IV aGvHD development. Using the final model, stochastic simulations were performed to explore covariate effects on the cumulative incidence over time and to estimate risk ratios. Conclusion Overall, the presented model showed good descriptive and predictive performance and provides valuable insights into the interplay of multiple static and time‐dependent risk factors for the prediction of aGvHD.

Relevant covariates were gathered from previous publications on the development of aGvHD [1][2][3][4][5][6][7][8][9] or selected by exploratory analysis of covariate-stratified Kaplan-Meier curves using the existing dataset (Supplementary Table 1).We performed multivariable covariate analysis on the base model using a significance level of 5% and analyzed both longitudinal data of cell counts and CsA blood levels, as well as static data of patient characteristics, donor characteristics, and the HSCT procedure.Effects of covariates were either implemented by directly modulating the base hazard or the hazard function parametrization.Continuous covariates were centred around the population median and tested via linear or power functions.Missing values of patient characteristics were replaced by population mode for discrete variables and by population median for continuous variables.If more than one measurement per day and patient was available, the median was used.
Missing longitudinal laboratory data were imputed using Last Observation Carried Forward  2).Fludarabine was always combined with TBI (2-12 Gy, FluTBI), fludarabine plus treosulfan was always administered without TBI, and busulfan plus fludarabine was primarily administered without TBI (99.6%, n=255/256).The effects of drug combinations and the use of TBI were tested separately.

Cross validation
To assess the final model regarding overfitting and stability, it was additionally evaluated by fivefold cross validation.For this, the complete dataset (n=1479) was split randomly into five equitably sized subgroups (subgroup 1-4: 296 IDs, subgroup 5: 295 IDs) using the R base function sample().Subgroups were checked for comparable numbers of events (observed numbers of events: 113, 115, 118, 112 and 107 for subgroup 1-5, respectively).Subsequently, the parameters of the final model were re-estimated for each fold.
The parameter distributions of the five re-estimated parameter sets and the parameters of the final model were compared (see supplementary Figure 2).Here, parameter distribution were simulated based on the estimated parameters and the corresponding residual standard error (RSE) using the mvrnorm() function from the R package MASS.
For model assessment, performance metrics were obtained by simulating each of the five test datasets using the corresponding parameters sets and apply log-rank test for simulated vs.
observed events.Subsequently an assessment of the model's performance was done by counting the number of significant log-rank tests results (p<0.05).

Stochastic simulations
To evaluate covariate effects on the cumulative incidence of grades II-IV aGvHD over time for 100 days after HSCT, we simulated 1000 replicates of the training dataset with the final TTE model.To isolate each covariate effect, all other covariates of the dataset were held constant.To examine the influence of covariates on the development of grades II-IV aGvHD within 100 days after transplantation, we calculated risk ratios (RRs) from simulations of 1000 replicates of the training dataset for each covariate, with the parameters of the other covariates set to "no effect".
For  Simulation of grade II-IV aGvHD cumulative incidence were conducted for patients without risk factors (younger patients, patients with later increase in WBC, patients treated with FluTBI, patients with MRD and female patients; blue curve) and with the respective risk factor (red curve).CsA blood levels were assumed to be constant at 200 ng/ml for both groups.To match the incidences between the two groups, CsA blood levels were increased for the risk groups as indicated in the legend (green curve).CsA=cyclosporine A, FluTBI= fludarabine plus total body irradiation, MRD=matched related donor, MUD=matched unrelated donor, WBC=white blood cells

(
LOCF).If data was missing at the time of transplantation, we used Next Observation Carried Backward (NOCB).CsA whole blood concentrations were obtained during clinical routine and quantified via antibody-conjugated magnetic immunoassay (ACMIA), hence, measurement frequency varied between individual patients.To avoid confounding effects of other immunosuppressive drugs, patients that switched from CsA prophylaxis to another immunosuppressive drug were censored at the last time of CsA measurement.The conditioning regimens consisted of 32 different drug combinations (23% fludarabine, 22% fludarabine/treosulfan, 17% busulfan/fludarabine) with or without total body irradiation (TBI) resulting in 51 different conditioning regimens (see Supplementary Table

Supplementary Figure 3 :Supplementary Figure 4 :
Performance evaluation on the cumulative incidence of developing grades II-IV aGvHD by five-fold cross validation.Model performance of the five iterations of cross validation process (sets 1-5) on the corresponding test (upper panel) and training dataset (lower panel).The red lines show the observed cumulative incidence.The blue lines show the median of model simulations.The blue shaded areas show the 95% confidence interval calculated from stochastic simulations of 1000 replicates.Potential CsA adjustment for patient with risk factors to achieve comparable grade II-IV aGvHD cumulative incidence compared with patients without risk factors.

Table 1 :
Tested covariates for aGvHD time-to-event analysis.
each discrete covariate, we divided simulations into a group containing individuals with the covariate (covariate group) and a group without the covariate (reference group) to calculate RRs.For continuous covariates, we grouped simulations by covariate cut-off values.To calculate RRs, we divided the proportion of individuals with an event in the covariate group by the proportion of individuals with an event in the reference group.Subsequently, we calculated the median, 2.5%, and 97.5% quantiles of all RRs.Supplementary † Time-dependent covariates were measured repeatedly and static covariates were assessed once before allo-HSCT.‡Early stages: De-novo AML in 1st remission, ALL in 1st remission, MDS with single lineage

Table 2 :
Overview of conditioning regimens of the training and test datasets grouped by frequency of application and sorted alphabetically.In total, the XplOit dataset consists of 52 different conditioning regimens.Infrequently observed regimens with different cumulative doses of radiation were summarized.

Table 5 :
Used ATG cumulative dosisdepending on donor type

Table 6 :
Used ATG cumulative dosisdepending on stem cell source

Table 7 :
Parameter estimates with relative standard error (%) of the final model (training dataset) in comparison to the cross validation fold (sets 1-5) CsA: cyclosporine A, f: scale parameter, ktr: hazard transit rate, kel: shape parameter (hazard elimination), RSE: relative standard error MRD *λ Sex mismatch *λ Flu * † γ WBC = effect of WBC count on h(t), λ MRD = effect of matched related donor on h(t), λ sex mismatch = effect of sex mismatch on h(t), λ Flu = effect of fludarabine on h(t) with h(t) = h 0 (t)*λ

Supplementary Figure 2: Parameters distribution of five-fold cross validation.
Parameter distribution of estimated model parameters for each iteration of cross validation (sets 1-5) and for the final model were simulated and overlaid.Different colors represent different parameter sets.