Risk Factors and Predictive Models for Intravenous Immunoglobulin Resistance in Children with Recurrent Kawasaki Disease

Introduction

Kawasaki disease (KD) is an acute and systemic vasculitis, predominantly affecting children under 5 years.¹ Systemic inflammatory response plays a crucial role in KD development and progression.² Coronary artery lesions (CALs) are the most significant complications of this disease, leading to KD being the main cause of acquired heart diseases in developed countries.³ Previous studies suggested the KD incidence, with an increasing trend and a marked geographic variation, is expected to be higher in Asian countries compared to western countries.⁴ Similarly, the rate of KD recurrence is higher in certain Asian countries, such as Japan and Korea (3–4%)^5–7 while lower in Caucasian population in US and Canada (1.5–1.7%).^8,9 Sudo et al have analyzed the data from a nationwide survey in Japan to compare the incidence rates of recurrent KD between different years and identify that it has remained largely unchanged over the past 30 years.⁵ The etiology of KD as well as its recurrence remains elusive. It seems to result from the interplay of genetic and environmental susceptibility factors with infectious triggers.¹⁰

High-dose intravenous immunoglobulin (IVIG) combined with aspirin is the current most effective therapy.¹¹ Timely initiating IVIG treatment remarkably reduces the rate of CALs from 25% to 4%.¹¹ However, approximately 10% to 20% of KD patients are resistant to IVIG treatment and are at increased risk for CALs.^11,12 The most common treatment regimens for IVIG-resistant KD are a repeat infusion of IVIG, retreatment with IVIG plus prednisone, or infliximab.¹³ Predicting the KD patients with high risk for IVIG resistance may help those patients get benefits from timely adjunctive therapy.¹¹ Several scoring systems including Kobayashi score, Egami score, and Sano score have been proposed to predict IVIG resistance in patients with initial KD onset.^14–19 Kobayashi score is the earliest well-known model and comprises seven variables including serum sodium, days of illness at initial treatment, aspartate aminotransferase (AST), neutrophil percentage (N%), c-reactive protein (CRP), age, and platelet count.¹⁴ Several Japanese clinical studies have demonstrated that in the patients with predicted high risk for IVIG resistance based on Sano score, Egami score, and Kobayashi score, adding corticosteroids to initial IVIG treatment is effective to prevent CAL development.^20–22 Based on the results, the guidelines from America and Japan recommend a combination of corticosteroids with IVIG for the initial treatment of KD patients with high risk for IVIG resistance.¹¹ Kobayashi score is also used in our practice.²³

Currently, however, no study has examined the risk factors for IVIG resistance in recurrent KD patients. In the present study, we retrospectively enrolled the recurrent KD patients over a 12-years observation period to explore the associations of clinical features/laboratory tests with IVIG resistance and further establish predictive models for IVIG resistance in recurrent KD patients.

Materials and Methods

Study Subjects

The clinical records of the consecutive KD children treated at our hospital from January 2010 to December 2021 were retrospectively reviewed. KD was diagnosed based on the guidelines proposed by the American Heart Association in 2004.²⁴ Recurrent KD was defined as intervals between the first and second episodes of more than 2 months.²⁵ The patients were treated with initial standard IVIG (2 g/kg) therapy. IVIG resistance was defined as persistent fever (an axillary temperature > 38°C) after 36 h post-initial IVIG treatment or recurrent fever within two weeks of IVIG treatment accompanied by one or more KD clinical manifestations.²⁴ Inclusion criteria: (1) <18 years; (2) diagnosed as complete or incomplete KD. Exclusion criteria: (1) with another rheumatic or infectious vasculitis; (2) recurrent within 2 months after the first episode; (3) suffer more than two episodes; (4) not receive initial IVIG therapy; (4) with incomplete data. This study was approved by the Ethics Committee of Beijing Children’s Hospital (2022-E-015-R) and complied with the Declaration of Helsinki. All data were fully anonymized and the requirement for informed consent was waived.

Data Collection

The patients’ basic demographic information was collected from the medical records including gender, age at the first episode, age at the second episode, interval, IVIG resistance at the first episode, physical examination results, and clinical manifestations such as fever duration, conjunctival injection, erythema of oral mucosa, rash, cervical lymphadenopathy, swelling or redness of the extremities. Laboratory data before IVIG administration at recurrence were also collected including white blood cell count (WBC), neutrophils, neutrophil percentage (N%), lymphocyte, lymphocyte percentage (L%), serum sodium (Na), total bilirubin (TB), direct bilirubin (DB), C-reactive protein (CRP), platelet (PLT), hemoglobin (Hb), erythrocyte sedimentation rate (ESR), serum albumin (ALB), serum alanine aminotransferase (ALT), serum aspartate transaminase (AST), prothrombin time (PT), creatine kinase MB (CK-MB), and international standardized ratio (INR).

Model Development

The patients in the entire cohort were randomly divided into the training and validation cohorts with a ratio of 7:3. Potential predictors were first selected by standard least absolute shrinkage and selection operator (Lasso) regression.²⁶ The coefficients of the variables were compressed along with the λ changes and the optimal λ was selected by cross-validation (lLasso method). Given the continuous variables might exhibit a non-linear relationship to the outcome. The continuous variables were transformed by restricted cubic splines (RCS) with 3 or 4 knots²⁷ and were subjected to group Lasso for feature selection (gLasso method).²⁸ The Lasso scores were calculated by multiplying the selected variables with corresponding coefficients and summing them and were then input into a standard logistic regression to obtain the predicted probability and linear predictors. The two models were subsequently compared in both the training and validation cohorts. A detailed description of the model development is presented in Supplementary Data.

Model Performance

Model discrimination performance was evaluated by the area under the receiver operating characteristics curves (AUC). AUCs of 0.5, 0.5–0.7, 0.7–0.8, 0.8–0.9, and >0.9 reflected no, poor, acceptable, excellent, and outstanding discriminations, respectively.²⁹ The sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and diagnostic accuracy were also assessed. Calibration was evaluated by the calibration slope, the calibration curve, and the Hosmer-Lemeshow test with ten groups. Validation was further performed by assessing the model performance in the validation cohort. Furthermore, the models were recalibrated by the intercepts and slopes in the validation cohort to obtain the final models.

Statistical Analysis

R software (Version 4.1.2) was used to perform the statistical analyses. Continuous variables were expressed as mean ± standard deviation (M±SD) or median with an interquartile range (25–75%). The categorical variables were expressed as cases number with a percentage. Continuous data between two groups were compared by Student’s t-test, Kruskal–Wallis Rank Sum Test, paired t-test, or Wilcoxon matched-pairs signed-rank test based on the normality of variable distribution. Categorical data between two groups were compared by Chi-square test or Fisher exact test. The associations of the clinical/laboratory variables with IVIG resistance at recurrence were also evaluated by univariate logistic regression analyses. Nomograms were also constructed by the RMS package to facilitate the easy calculation of risk scores and the risk of IVIG resistance at recurrence.

Results

Demographic Data

During the study period of 12 years, a total of 6125 children with a diagnosis of KD were admitted to our hospital. Of them, 92 cases were recurrent KD. According to the inclusion and exclusion criteria, two cases were excluded due to incomplete data, and 90 cases were finally included in this study. Sixty-five were boys and 25 were girls. The median ages at the first episode and second episode were 22.5 and 34.0 months, respectively. The median interval between the two episodes was 8.5 months. Fifty-four (60.0%) cases recurred within 12 months after the first episode.

Risk Factor for IVIG Resistance in Recurrent KD

Of the 90 recurrent KD diseases, 74 cases were IVIG responsive while 16 cases were IVIG resistant (Table 1). Comparison of the clinical features and pre-treated laboratory data at the second episode according to IVIG response revealed the patients with IVIG resistance had higher age and IVIG resistance probability at the first episode, increased CRP levels, neutrophils count, N%, DB level, prothrombin time, and INR, and decreased lymphocyte count, L%, and serum sodium levels (P < 0.05, Table 1), indicating these variables might be risk factors for IVIG resistance in recurrent KD, which were also verified by the univariate Logistic regression analyses (Table 2).

Table 1 Comparison of the Clinical Features and Pre-Treated Laboratory Data According to IVIG Resistance in Children with Recurrent KD

Table 2 Associations of the Clinical/Laboratory Variables with IVIG Resistance in Patients with Recurrent Kawasaki Disease Analyzed by Univariate Logistic Regression Model

Construction of Predictive Models for IVIG Resistance in Recurrent KD

The patients were randomly divided into the training cohort (n = 63) for model development and the validation cohort (n = 27) for model validation. The data distribution in the training and validation cohorts were similar (Table S1). The continuous variables were first assumed to have a linear relationship to the IVIG resistance at recurrence and subjected to standard Lasso logistic regression analysis in the training cohort combined with other categorical variables (lLasso model). Five variables including age and IVIG resistance at the first episode, lymphocyte count, and serum sodium levels were finally selected at the optimal λ with the smallest cross-validation deviance (Figure 1). Next, all continuous variables were assumed to have a potential non-linear relationship with the outcome and expanded to restricted cubic splines. Then all the variables were selected by group Lasso (one original variable per group). Four variables including age and IVIG resistance at the first episode, N%, and CRP levels were selected at the optimal λ with the smallest cross-validation error based on deviance (Figure 2).

Figure 1 Variables selection and coefficient estimation by standard Lasso logistic regression analysis to construct the lLasso model. (A) A total of 33 variables were assumed to have a linear relationship to the IVIG resistance at recurrence and input to standard Lasso logistic regression analysis. The coefficient compress path along with penalty parameter λ changes was presented. The dashed vertical line indicated the optimal λ obtained from 10-fold cross-validation. (B) Ten-fold cross-validation was performed to select the optimal λ with the minimized binomial deviance. The dashed vertical line indicated the optimal λ.

Figure 2 Variables selection and coefficient estimation by group Lasso logistic regression analysis to construct the gLasso model. (A) The continuous variables were assumed to have a potential non-linear relationship to the IVIG resistance at recurrence and input to group Lasso logistic regression analysis combined with other categorical variables. The coefficient compress path along with penalty parameter λ changes was presented. The dashed vertical line indicated the optimal λ obtained from 10-fold cross-validation. (B) Ten-fold cross-validation was performed to select the optimal λ with the minimized cross-validation error, which was also based on binomial deviance. The dashed vertical line indicated the optimal λ.

Model Performance

ROC curves suggested the AUCs of lLasso and gLasso models were similar in the training cohort (0.895 vs 0.906, P = 0.790; Figure 3A and Table 3) and both models exerted excellent discrimination. Brier score and P-value examined by Hosmer-Lemeshow test of lLasso model were 0.094 and 0.480, respectively, and that of gLasso model were 0.085 and 0.530, respectively, indicating good performance in both models.

Table 3 Model Performances Comparison in the Training Cohort and the Validation Cohort

Figure 3 Comparison of the ROC curves of the lLasso and gLasso models in the training and validation cohorts. (A) Comparison of the ROCs of the lLasso and gLasso models in the training cohort. (B) Comparison of the ROCs of the lLasso and gLasso models in the validation cohort.

Validation

Both models were assessed in the validation cohort. ROC curves suggested the AUC of gLasso was slightly bigger than that of lLasso but without significant difference (0.909 vs 0.855, P = 0.439), indicating both models still exerted excellent discrimination in the validation cohort. Brier score, calibration slope, and P-value examined by Hosmer-Lemeshow test of gLasso model were 0.104, 0.834, and 0.380, respectively, indicating an acceptable calibration, which was also supported by the calibration curves (Figure 4), and was slightly prior to lLasso model whose Brier score, calibration slope, and P-value examined by Hosmer-Lemeshow test were 0.109, 0.788, and 0.168, respectively (Figure 4A and B). Then the models were recalibrated by the slope and intercept from the validation cohort to obtain the final models. After recalibration, Brier score and P-value examined by Hosmer-Lemeshow test of gLasso model were 0.088 and 0.491, respectively, and that of lLasso model were 0.084 and 0.173, respectively (Figure 4C and D). The linear predictor (LP) of the lLasso model can be calculated by the formula (Figure S1A): LP_lLasso=0.0253*age at the first episode - 0.3149*lymphocytes counts + 0.0210*CRP - 0.0248*serum sodium levels + 1.4797* IVIG resistance at the first episode - 0.8411. LP of the gLasso model can be calculated by the formula (Figure S1B): LP_gLasso=0.1081*age at the first episode - 0.4660*[(age at the first episode)’] + 1.1593*[(age at the first episode)”] - 0.0299*N% + 0.0652*[(N%)’] + 0.0091*CRP + 0.0129*[CRP’] + 1.8384*IVIG resistance at the first episode - 4.7177. The continuous variables, age at the first episode, N%, and CRP, were modeled using RCS with 4, 3, and 3 knots, respectively, resulting in extra parameters, which were dependent on the original variables and could be calculated by the formula in Figure S1B.

Figure 4 Calibration curves of lLasso and gLasso models in the validation cohort before and after recalibration. (A) The lLasso model before recalibration. (B) The gLasso model before recalibration. (C) The lLasso model after recalibration. (D) The gLasso model after recalibration.

Application of the Risk Models

Given lLasso and gLasso models exhibited similar performance; thus, two nomograms were established based on the two models to facilitate the potential application of the predictive models (Figure 5). With the nomograms, the probability of IVIG resistance at recurrence can be estimated as follows: firstly, draw a straight line from the predictor up to the “Point” line to obtain the points of each predictor. Secondly, sum all the rewarded points to get total points. Thirdly, draw a straight line from the “Total Points” line down to the “Risk of IVIG resistance” line to obtain the chance of the patient developing IVIG resistance at recurrence.

Figure 5 Nomograms constructions based on lLasso model and gLasso model. (A). The lLasso model consists of five variables, including age at the first episode, lymphocyte, CRP, serum sodium, and IVIG resistance at the first episode. (B). The gLasso model consists of four variables, including age at the first episode, neutrophils percentage, CRP, and IVIG resistance at the first episode. The probability of IVIG resistance at recurrence can be estimated as follows: firstly, draw a straight line from the predictor up to the “point” line to obtain the points of each predictor. Secondly, sum all the rewarded points to get total points. Thirdly, draw a straight line from the “total points” line down to the “risk of IVIG resistance” line to obtain the chance of the patient developing IVIG resistance at recurrence.

Discussion

Previous studies suggested the KD recurrence rate varies among different regions and ethnicities, with a range from 1.5% in Canada,⁹ 1.7% in the USA,⁸ to 3–4.6% in Japan.^5,6,30,31 In the USA, the recurrence rates are higher in children of Asian descent compared to the Caucasian population, similar to that in Asian regions. In China, large epidemiologic surveys of recurrent KD are rare. Individual cohort studies indicated that the recurrence rate of KD ranges from 1.34% to 3.10% in China.^32–35 The KD recurrence rate in our study was 1.5%, consistent with the previous reports. The risk factors of KD recurrence are poorly understood. Yang et al found longer fever duration, lower Hb levels, and higher AST levels are independent risk factors for KD recurrence.³³ In addition, Sudo et al identified male sex, young age, and initial IVIG resistance as risk factors for KD recurrence.⁵

Approximately 10–20% of cases are resistant to IVIG therapy in children with initial KD onset.^36,37 Rare studies have reported the IVIG resistance rates in recurrent KD and its first episode.^33,38 Yang et al³³ and Kang et al³⁸ included 22 and 25 recurrent KD cases in their publications, respectively, and found the IVIG resistance rates in the recurrent episode and initial episode are similar (18.1% vs 22.7% and 20% vs 16%, respectively). In our study, we found 17.8% of cases were IVIG resistant at the second episode and not significantly different from that at the first episode (11.1%). In children with initial KD onset, male, swelling of extremities, rash, decreased Hb, PLT, ALB, and serum sodium, and increased neutrophils, ESR, CRP, total bilirubin, ALT, and AST have been found to be potential risk factors for IVIG resistance.³⁹ Based on the risk factors, several scoring models have been established to predict IVIG resistance in children with initial KD onset.^{14–16,18,40} Up to now, no study has been performed to investigate the risk factors for IVIG resistance in children with recurrent KD. In the present study, the patients with IVIG resistance at recurrence had higher age at the first episode, CRP levels, neutrophils count, neutrophils percentage, DB level, prothrombin time, international normalized ratio, and IVIG resistance proportion at first episode while having lower lymphocyte counts, lymphocyte percentage, and serum sodium levels.

Due to the rare KD recurrence (1.5% in our study), only 90 children with recurrent KD were included during a long observation period of 12 years in one of the largest pediatric hospitals in China. Given the limited sample size and up to 33 clinical/laboratory variables as candidate predictors were included in our study, model development based on conventional methods such as stepwise regression based on the Akaike information criterion (AIC) might be not appropriate. Thus, we performed variable selection and coefficient estimation using standard Lasso and its variant, group Lasso, and those methods could also be used for biomarker selection from high-dimensional data. The continuous variables were assumed to have a linear relationship to outcome in standard Lasso analysis while were assumed to have a potential non-linear relationship with the outcome and expanded to restricted cubic splines in group Lasso analysis. Five variables including age at the first episode, lymphocytes count, serum sodium levels, CRP levels, and IVIG resistance at the first episode were finally selected in the lLasso model whereas four variables including age at the first episode, N%, CRP levels, and IVIG resistance at first episode were selected in the gLasso model. The two models shared three variables. Through the gLasso model, we found age at the first episode and N% were associated with the outcome in a non-linear way and an approximately linear relationship was observed for CRP. The risk of IVIG resistance at recurrence increased quickly along with age at the first episode before about 20 months, remained similar till to 60 months, and increased again. The patients with N% from 40% to 60% at a normal reference range had a lower risk of IVIG resistance while high or low N% levels especially high N% were associated with increased IVIG resistance risk. Notably, N%, serum sodium, and CRP have been identified as predictors in one or more previous models developed to predict IVIG resistance in patients with initial KD onset. Both lLasso and gLasso models exerted good discrimination in both the training cohort (AUC, 0.895 vs 0.906) and the validation cohort (AUC, 0.855 vs 0.909). The gLasso model performed slightly better than the lLasso model in the validation cohort; however, the difference was not significant (DeLong test, P = 0.439). As the patients were collected from a long observation period (12 years), the two models were compared in subgroups of the early period (2010–2018) and the late period (2019–2021) and we found both lLasso and gLasso models performed well in either early period subgroup (AUC, 0.857 vs 0.860) or late period subgroup (AUC, 0.889 vs 0.960). Notably, all the variables included in lLasso and gLasso models were conventional and could be easily obtained. Thus, we recommended validating both models in external cohorts in the future. And then, the models can be used to identify the recurrent KD patients with high risk of IVIG resistance, who may benefit from adding corticosteroids to the initial IVIG treatment after KD recurrence.

To our knowledge, this was the first study to develop and validate predictive models for IVIG resistance in children with recurrent KD; however, there were some shortages in the study that needed to be considered. Firstly, although we collected the study subjects from one of the largest children’s hospitals in China over a 12-years observation period, the sample size was small due to the extremely low recurrence rate of KD. Secondly, this was a retrospective study performed in a single center and only internal validation was performed. Further external validation in prospective and multi-centric studies with a large sample size was needed.

Conclusion

Age and IVIG resistance at the first episode and some laboratory variables at recurrence are identified as potential risk factors for IVIG resistance in patients with recurrent KD. Furthermore, two predictive models for IVIG resistance in children with recurrent KD are successfully constructed and internally validated. However, external validation should be performed before clinical use.