Prediction of intravenous immunoglobulin retreatment in children with Kawasaki disease using models combining lymphocyte subset and cytokine profile in an East Asian cohort

Abstract

Objectives

For children with Kawasaki disease (KD) at high risk of developing coronary artery lesions and requiring retreatment with intravenous immunoglobulin (IVIG), the availability of accurate prediction models remains limited because of inconsistent variables and unsatisfactory prediction results. We aimed to construct models to predict patient's probability of IVIG retreatment combining children's individual inflammatory characteristics.

Methods

Clinical manifestations and laboratory examinations of 266 children with KD were retrospectively analysed to build a development cohort data set (DC) and a validation cohort data set (VC). In the DC, binary logistic regression analyses were performed using R language. Nomograms and receiver operating curves were plotted. The concordance index (C index), net reclassification index, integrated discrimination improvement index and confusion matrix were applied to evaluate and validate the models.

Results

Models_5V and _9V were established. Both contained variables including the percentages of CD8⁺ T cells, CD4⁺ T cells, CD3⁺ T cells, levels of interleukin (IL)-2R and CRP. Model_9V additionally included variables for IL-6, TNF-α, NT-proBNP and sex, with a C index of 0.86 (95% CI 0.79–0.92). When model_9V was compared with model_5V, the NRI and IDI were 0.15 (95% CI 0.01–0.30, P < 0.01) and 0.07 (95% CI 0.02–0.12, P < 0.01). In the VC, the sensitivity, specificity and precision of model_9V were 1, 0.875 and 0.667, while those of model_5V were 0.833, 0.875 and 0.625.

Conclusion

Model_9V combined cytokine profiles and lymphocyte subsets with clinical characteristics and was superior to model_5V achieving satisfactory predictive power and providing a novel strategy early to identify patients who needed IVIG retreatment.

Introduction

Kawasaki disease (KD) is an acute, auto-inflammatory and systemic form of vasculitis that occurs in children. Cardiovascular sequelae are the primary causes of acquired paediatric heart disease.^1,2 Evidence of immunological imbalance during the acute phase of the disease has been previously reported, specifically increased cytokine production, such as tumor necrosis factor-α (TNF-α) and interleukin-1 (IL-1), and abnormal distribution of certain lymphocyte subsets, such as CD8⁺ T cells.¹ The standard therapy for KD has now been established as a single high dose of intravenous immunoglobulin (IVIG) together with aspirin during the initial phases of KD. Although the application of this regimen greatly reduces inflammation and arterial damage in children with KD, 15–20% of patients continue to have persistent or recurrent fever after primary treatment who are identified as resistant to IVIG and at high risk of developing coronary artery abnormalities (CAA) such as aneurysms or thrombosis.^1–4 Among various therapeutic approaches for rescue, retreatment with IVIG remains the predominant option, which achieves confirmed curative effect.^1,3,5,6 However, in clinical practice, the decision to administer IVIG again will not be made until the symptoms of KD have recurred several days after the initial therapy, usually 2–5 days later, as persistent or recurred fever and symptoms remain the primary criteria for retreatment with IVIG, which increases the risk of inflammatory aggression and delayed treatment. Therefore, patients who require IVIG retreatment should be recognised as early as possible to prevent the development of inflammation and coronary artery lesions.

Previous studies have established predictive models for children with KD at risk of IVIG resistance, identifying various risk factors including serum sodium, C-reactive protein (CRP) and total bilirubin (TBil) and others.^7–10 However, there is significant variation in the identified independent risk factors for KD among children from different cities within the same region.^7–11 In addition, these prediction models tended to focus more on clinical features and indicators rather than on items associated with specific aspects of inflammation.^11–15 The diverse clinical perspectives on KD among researchers result in different factors being included in predictive models, and the varying number of cases can also lead to the development of distinct models. In fact, increasing evidence shows the pathogenesis of KD is mainly promoted by abnormal and unbalanced innate and adaptive immune responses.¹⁶ Cytokines including IL-1, IL-2, IL-6, TNF-α, IFN-γ and so on have been studied in many researches of KD.¹⁷ Previous investigations, including our own work, revealed that during the acute phase of KD and before treatment, a child's count of CD19⁺ B cells increases markedly and CD4⁺ T cells decreases significantly suggesting the predominance of humoral immunity.^11,18,19 Hyperactive cytokines and unbalanced distribution of lymphocyte subgroups constitute the inflammatory feature of KD. Therefore, cytokines and subgroups of lymphocytes are routinely tested as part of a comprehensive laboratory evaluation of KD in most of the medical centres.^16,20 However, the reported investigation of KD lacks an overall understanding of the fluctuations of these inflammatory cytokines and lymphocyte subsets. This study, conducted at our paediatric medical centre in East Asia, analysed numerous variables such as clinical features, standard laboratory indicators and inflammatory markers of KD patients. The goal was to develop a predictive model for IVIG retreatment applicable to children with KD, facilitating the rapid development of a personalised IVIG treatment plan during the early stages of the disease.

Results Baseline characteristics of children with KD in the DC and in the VC were similar

The baseline profile, clinical characteristics, laboratory indices, percentage of lymphocyte subsets and cytokine levels in both DC and VC are summarised in Tables 1 and 2, respectively. In DC, out of 236 children with KD, 38 patients (16.1%) received IVIG retreatment. In the VC, 6 out of 30 children with KD (20%) required IVIG retreatment. In children with KD in the DC, the percentage of neutrophils and the levels of CRP and NT-proBNP were significantly higher in the IVIG retreatment group than those in the single-dose IVIG group. In contrast, PLT and serum sodium levels were significantly lower in those in the single-dose IVIG group. In terms of cytokines, IL-2R, IL-10, IL-6 and TNF-α levels were significantly higher in patients receiving IVIG retreatment than those in the single-dose IVIG group. In the VC (Table 2), although there was no significant difference in the changes in cytokine levels between the two groups, the cytokine levels in the IVIG retreatment group showed an increasing trend, similar to the cytokines in the DC, except for IL-2R. The IL-2R levels in the IVIG retreatment group were significantly higher than those in the single-dose IVIG group, which was similar to that in the DC group (P = 0.001).

Table 1 Baseline characteristics and differences between patients in the group of single-dose IVIG and in the IVIG retreatment group in the DC

	All children with KD	Single dose of IVIG	IVIG retreatment	P-value
Patients	236	198 (83.90%)	38 (16.10%)
Age, months	25.00 (13.00–44.75)	24.50 (13.00–44.25)	32.50 (12.00–46.00)	0.51
Sex
Male	148 (62.71%)	126 (63.64%)	22 (57.89%)	0.58
Female	88 (37.29%)	72 (36.36%)	16 (42.11%)
Incomplete Kawasaki disease	78 (33.05%)	68 (34.34%)	10 (26.32%)	0.45
Duration of fever before the first IVIG, days	6.00 (5.00–7.00)	6.00 (5.00–7.00)	6.00 (5.00–7.00)	0.55
Laboratory data at diagnosis
Neutrophils (50–70%)	67.85% (55.9%–77.10%)	66.70% (54.10%–75.40%)	75.85% (67.75%–83.48%)	0.00020
C-reactive protein (< 8 mg L−1)	70.00 (34.00–107.75)	66.00 (33.00–96.25)	106.00 (55.25–161.00)	0.00050
Erythrocyte sedimentation rate (0–15 mm h−1)	57.71 (28.22)	58.22 (28.56)	55.05 (26.58)	0.53
Platelets (the lowest level among tests before the first IVIG; 100–300 × 109 L−1)	273.50 (217.00–341.75)	278.50 (223.75–349.75)	238.00 (194.50–311.00)	0.0091
Serum sodium (135–145 mmol L−1)	136.30 (133.73–138.00)	136.50 (134.20–138.40)	134.45 (131.88–136.15)	0.00030
NT-proBNP (< 285 pg mL−1)	603.25 (262.93–1419.75)	554.45 (225.30–1409.25)	812.50 (435.63–1722.25)	0.034
Abnormal cases of cardiovascular ultrasound	31 (13.14%)	28 (14.14%)	3 (7.89%)	0.43
Inflammatory cytokine
IL-2R (223–710 U mL−1)	2128.00 (1310.00–3302.50)	1928.50 (1257.00–2985.25)	3649.00 (2558.50–5522.25)	< 0.0001
IL-10 (< 9.1 pg mL−1)	18.05 (7.11–54.73)	16.15 (6.53–39.60)	63.45 (21.08–202.25)	< 0.0001
IL-6 (< 5.9 pg mL−1)	80.55 (33.83–201.75)	65.80 (26.35–152.75)	198.00 (78.98–370.75)	0.00020
IL-1β (< 5 pg mL−1)	8.99 (4.99–21.70)	8.82 (4.99–21.75)	9.93 (4.99–19.55)	0.59
IL-8 (< 62 pg mL−1)	59.10 (17.73–384.75)	54.10 (16.40–415.75)	70.75 (28.08–261.00)	0.44
TNF-α (< 8.1 pg mL−1)	25.00 (17.60–34.88)	24.10 (17.03–33.75)	29.15 (22.38–43.23)	0.027
Lymphocyte subset
CD19+ B cell (6.8–15.8%)	27.94% (22.33–39.09%)	26.69% (21.52–36.25%)	41.79% (27.69–51.99%)	< 0.0001
CD16+CD56+ T cell (9.5–23.5%)	7.60% (4.41–11.16%)	7.60% (4.42–11.52%)	7.65% (3.35–9.58%)	0.35
CD8+ T cell (18.2–32.8%)	16.85% (13.80–21.93%)	17.18% (14.26–22.54%)	14.04% (10.90–18.69%)	0.00080
CD4+ T cell (29.4–45.8%)	36.45% (9.97%)	37.69% (9.50%)	29.99% (9.96%)	< 0.0001
CD3+ T cell (60.8–75.4%)	59.49% (51.31–67.32%)	60.90% (52.57–68.37%)	48.05% (38.05–58.94%)	< 0.0001

DC, development cohort data set; IL, interleukin; IQR, interquartile range; IVIG, intravenous immunoglobulin; KD, Kawasaki disease; NT-proBNP, N-terminal precursor brain natriuretic peptide; TNF, tumor necrosis factor.

Data are presented as median (IQR), n (%) or mean (SD), unless otherwise specified.

Table 2 Baseline characteristics and differences between patients in the group of single-dose IVIG and IVIG retreatment in the VC

	All children with KD	Single dose of IVIG	IVIG retreatment	P-value
Patients	30	24 (80.00%)	6 (20.00%)
Age (months)	41.70 (25.75)	40.42 (26.05)	46.83 (26.18)	0.59
Sex
Male	17 (56.67%)	12 (40.00%)	5 (16.67%)	0.20
Female	13 (43.33%)	12 (40.00%)	1 (3.33%)
Laboratory data at diagnosis
Percentage of neutrophil (50–70%)	69.83 (15.83)	69.38 (10.38)	71.67 (30.86)	0.76
C-reactive protein (< 8 mg L−1)	77.00 (28.75–114.50)	65.00 (28.25–104.50)	135.00 (62.00–201.00)	0.12
Erythrocyte sedimentation rate (0–15 mm h−1)	64.37 (33.18)	67.58 (32.04)	51.50 (37.64)	0.30
Platelet (the lowest level among several tests before the first IVIG) (100–300 × 109 L−1)	369.00 (283.75–545.00)	409.92 (161.32)	520.33 (220.03)	0.32
Serum sodium (134–143 mmol L−1)	135.70 (4.23)	136.62 (3.21)	132.05 (36.44)	0.015
NT-proBNP (< 285 pg mL−1)	377.05 (156.98–1265.41)	261.95 (146.13–544.85)	463.38 (127.00–912.87)	0.13
Inflammatory cytokine
IL-2R (223–710 U mL−1)	2372.50 (1246.75–3802.25)	1835.50 (976.00–3184.25)	4281.00 (3628.00–7501.00)	0.0010
IL-10 (< 9.1 pg mL−1)	19.10 (5.36–44.60)	10.13 (4.90–43.58)	30.75 (22.73–145.45)	0.073
IL-6 (< 5.9 pg mL−1)	48.15 (13.28–122.75)	34.00 (13.13–118.00)	94.55 (11.60–204.25)	0.43
IL-1β (< 5 pg mL−1)	14.25 (7.36–25.65)	14.25 (6.79–23.23)	18.65 (9.93–44.43)	0.60
IL-8 (< 62 pg mL−1)	39.60 (11.95–96.38)	32.00 (11.73–87.63)	57.10 (25.03–418.50)	0.25
TNF-α (< 8.1 pg mL−1)	23.80 (15.98–32.40)	22.55 (14.85–33.68)	26.10 (22.75–37.99)	0.25
Lymphocyte subset
CD19+ B cell (6.8–15.8%)	31.49% (12.06%)	30.49% (11.13%)	35.49% (15.82%)	0.37
CD16+CD56+ T cell (9.5–23.5%)	6.03% (3.86–12.40%)	6.62% (4.09–12.34%)	5.62% (2.64–12.66%)	0.67
CD8+ T cell (18.2–32.8%)	18.21% (5.84%)	19.20% (5.61%)	14.25% (5.45%)	0.062
CD4+ T (29.4–45.8%)	35.39% (8.68%)	34.88% (8.36%)	37.43% (10.44%)	0.53
CD3+ T cell (60.8–75.4%)	57.98% (11.35%)	58.72% (10.94%)	54.99% (13.53%)	0.48

IL, interleukin; IVIG, intravenous immunoglobulin; KD, Kawasaki disease; NT-proBNP, N-terminal precursor brain natriuretic peptide; TNF, tumor necrosis factor; VC, validation cohort data set.

Several studies have previously revealed significant alterations in lymphocyte subsets during the development of inflammation in KD.^19,21–24 In the present paper, as shown in Table 1, the percentage of CD19⁺ B cells in the IVIG retreatment group was significantly higher than that in the single-dose IVIG group, whilst the percentages of CD8⁺, CD4⁺ and CD3⁺ T cells were significantly lower compared with those in the single-dose IVIG group. The distribution of lymphocyte subgroups in the VC tended to be similar to that in the DC, except for the changing trend of CD4⁺ T cells but without significance between the single-dose IVIG group and retreatment group.

Variables analysis and prediction model construction

In the DC, we constructed ‘model_5V’ (abbreviated as mod_5V), which includes five variables, to estimate the likelihood of IVIG retreatment of children with KD. The percentage of CD8⁺ T cells, CD4⁺ T cells, CD3⁺ T cells, and serum levels of IL-2R and CRP were identified as significant independent predictors of IVIG retreatment requirement, producing the following function: g (X) = 0.2046–0.2708 × CD8⁺ T cells (%) – 0.2479 × CD4⁺ T cells (%) + 0.1759 × CD3⁺ T cells (%) + 0.0003 × IL2R + 0.0074 × CRP.

In this formula, g(X) is the logistic function of IVIG retreatment probability, and the odds ratios of the variables are shown in Table 3. The results revealed that lower percentages of CD8⁺ (OR = 0.7628, P = 0.0047) and CD4⁺ T cells (OR = 0.7804, P = 0.0012), higher percentages of CD3⁺ T cells (OR = 1.1923, P = 0.017), and increased levels of IL-2R (OR = 1.0003, P = 0.02) and CRP (OR = 1.0075, P = 0.067), are significantly associated with the requirement for IVIG retreatment. Mod_5V is shown as a nomogram in Figure 1a. Figure 1b shows an example of a 15-month-old girl using a five-item predictive nomogram. A girl with KD received 298 points, and the predicted probability of IVIG retreatment was 0.0387. The patient received only one dose of IVIG. Figure 1c shows another example of a boy with KD, aged 5 years and 2 months, whose total points were calculated to be 334, with a predicted probability of IVIG retreatment of 0.924. The patient received a second dose of IVIG 6 days after the first dose.

Table 3 Regression coefficients, odds ratio and the 95% confidence interval of the prediction model_5V (mod_5V) and model_9V (mod_9V)

		Regression coefficients	Odds ratio (95% confidence interval)	P-value
Mod_5V	Intercept	0.2046 (1.2065)	1.2270 (0.1116–13.4611)	0.87
	CD8	−0.2708 (0.0958)	0.7628 (0.6266–0.9194)	0.0047
	CD4	−0.2479 (0.0767)	0.7804 (0.6677–0.9095)	0.0012
	CD3	0.1759 (0.0735)	1.1923 (1.0281–1.3836)	0.017
	IL2R	0.0003 (0.0001)	1.0003 (1.0000–1.0005)	0.020
	CRP	0.0074 (0.0041)	1.0075 (1.0000–1.0156)	0.067
Mod_9V	Intercept	−0.5431 (1.5385)	0.5809 (0.0283–12.2770)	0.72
	Sex	1.0739 (4.8042)	2.9268 (1.1639–7.7837)	0.025
	CD8	−0.2727 (1.0096)	0.7613 (0.6196–0.9277)	0.0069
	CD4	−0.2697 (0.0830)	0.7636 (0.6456–0.9016)	0.0012
	CD3	0.1951 (0.0788)	1.2155 (1.0359–1.4237)	0.013
	IL2R	0.0005 (0.0001)	1.0005 (1.0002–1.0007)	0.0018
	TNF-α	−0.0249 (0.0145)	0.9754 (0.9434–0.9973)	0.086
	IL6	0.0020 (0.0012)	1.0020 (1.0002–1.0044)	0.086
	CRP	0.0090 (0.0042)	1.0090 (1.0008–1.0175)	0.032
	NT-proBNP	−0.0001 (0.0001)	0.9999 (0.9996–1.0000)	0.14

Figure 1. Nomogram predicting intravenous immunoglobulin (IVIG) retreatment probability using model_5V (mod_5V) and examples. For each predictor, we read the points assigned on a 0–100 scale at the top before adding these points. Subsequently, find the number on the ‘Total Points’ scale and read the corresponding predictions of ‘Prediction probability f+98-**61or IVIG retreatment’. *P [less than] 0.05; **P [less than] 0.01. (a) Nomogram for mod_5V. (b) Nomogram using mod_5V for one negative example (child treated with IVIG one dose). A 15-month-old girl with KD had a C-reactive protein (CRP) level of 93 mg L−1, interleukin-2 receptor (IL-2R) level of 1675 U mL−1, CD8+ T-cell percentage of 18.51%, CD3+ T-cell percentage of 63.12% and CD4+ T-cell percentage of 43.08%. Using the mod_5V nomogram, 298 points of 298 were ts of 334 was ac obtained, which was equal to a prediction probability of 3.87% for IVIG retreatment. (c) Nomogram using mod_5V for one positive example (child treated with IVIG retreatment). An example of a boy with KD aged 5 years and 2 months: CRP level, 201 mg L−1; IL-2R level, 7501 U mL−1, CD8+ T-cell percentage, 10.58%; CD3+ T-cell percentage, 26.57% and CD4+ T-cell percentage, 12.78%. A total point of 334 was achieved by mod_5V, and the prediction probability was 92.4% for IVIG retreatment, indicating that it was highly likely that this child would need IVIG retreatment.

Subsequently, we developed another predictive model incorporating nine variables using stepwise backward or backward and forward elimination in multivariate logistic regression, selected from the initial pool of 23 variables. The results showed that sex (OR = 2.9268, P = 0.025), percentage of CD8⁺ T cells (OR = 0.7613, P = 0.0069), CD4⁺ T cells (OR = 0.7636, P = 0.0012), CD3⁺ T cells (OR = 1.2155, P = 0.013), IL-2R (OR = 1.0005, P = 0.0018), TNF-α (OR = 0.9754, P = 0.086), IL-6 (OR = 1.0020, P = 0.086), CRP (OR = 1.0090, P = 0.032) and NT-proBNP (OR = 0.9999, P = 0.14) were independent influencing factors for IVIG retreatment (Table 3). The prediction model, termed ‘model_9V’ (mod_9V), utilises the following formula: g (X) = −0.5432 + 1.0739 × sex – 0.2727 × CD8⁺ T cells (%) – 0.2697 × CD4⁺ T cells (%) + 0.1951 × CD3⁺ T cells (%) + 0.0005 × IL-2R + 0.009 × CRP – 0.0249 × TNF-α + 0.002× IL-6 – 0.0001 × NT-proBNP. The regression coefficient of each factor is expressed as the change in the probability of IVIG retreatment when the predictor increased by one unit at the same time the value of other predictors in the model is fixed after computer fitting using a logistic regression model.

In addition to the five variables in mod_5V, sex (OR = 2.9268, P = 0.025), IL-6 level (OR = 1.002, P = 0.086), NT-proBNP level (OR = 0.9999, P = 0.14) and TNF-α level (OR = 0.9754, P = 0.086) were significant predictors of IVIG retreatment in mod_9V. Sex is a binary variable with 0 representing boy and 1 representing girl. Mod_9V showed that when a girl patient's percentages of CD8⁺ T cells and CD4⁺ T cells were low, the percentage of CD3⁺ T cells was high, the levels of IL-2R, IL-6 and CRP were increased, and the levels of TNF-α and NT-proBNP were decreased, the possibility of IVIG retreatment was greater. Figure 2a shows the nomogram for Mod _9V. Figure 2b represents the prediction probability from the same patient as Figure 1b, in which the prediction probability of IVIG retreatment was 0.0723. Figure 2c shows an example of a boy with KD aged 6 months who yielded a prediction probability of IVIG retreatment of 1. The patient's clinical situation was consistent with his prediction probability, as he received IVIG retreatment 6 days after the first IVIG therapy.

Figure 2. Nomogram for predicting IVIG retreatment probability for model_9V (mod_9V) and examples. For each predictor, we read the points assigned on a 0–100 scale at the top before adding these points. Find the number on the ‘Total Points’ scale and read the corresponding predictions of ‘Prediction probability for IVIG retreatment’. *P [less than] 0.05; **P [less than] 0.01. (a) Nomogram for mod_9V. (b) Nomogram using mod_9V for one negative example (child treated with IVIG one dose). Example of a 15-month-old girl with KD: N-terminal pro-B-type natriuretic peptide (NT-proBNP) level of 557.5 pg mL−1, CRP level, 93 mg L−1; IL-2R level, 1675 U mL−1; tumor necrosis factor-α (TNF-α) level of 32.6 pg mL−1, IL-6 level, 84 pg mL−1; CD8+ T-cell percentage, 18.51%; CD3+ T-cell percentage, 63.12%; and CD4+ T-cell percentage, 43.08%. The patient achieved a total point of 361, with a probability of 7.23% for IVIG retreatment. (c) Nomogram using mod_9V for one positive example (child treated with IVIG retreatment). Example of a boy with KD aged 6 months: NT-proBNP level of 672.4 pg mL−1, CRP 10 mg L−1, IL-2R level of 452 U mL−1, TNF-α level of 38.9 pg mL−1, IL-6 level of 7501 pg mL−1; CD8+ T-cell percentage, 2.65%; CD3+ T-cell percentage, 16.18% and CD4+ T-cell percentage, 12.4%. A total of 423 points was calculated, indicating that the probability of IVIG retreatment was 100%.

Comparison and evaluation of prediction models

The maximum likelihood ratio test was conducted for both models. The χ² values were calculated as 48.34 for mod_5V and 61.69 for mod_9V. The P-values of mod_5V and mod_9V were 3.02 × 10⁻⁹ and 6.32 × 10⁻¹⁰, respectively, both of which were lower than 0.0001, indicating that the two models were significantly different from the empty model, which contained none of the variables (Table 4).

Table 4 Evaluations of mod_5V and mod_9V

	mod_5V	mod_9V
Maximum likelihood ratio test (χ2)	48.34	61.69
P-value	< 0.0001	< 0.0001
Hosmer-Lemeshow test (χ2) for GOF	5.13	9.23
P-value	0.74	0.32
P-value of ANOVAa	0.0097
ROC auc (95% CI)	0.82 (0.74–0.90)	0.86 (0.79–0.92)
Specificity	81.82%	77.78%
Sensitivity	73.68%	78.95%
PRauc	0.52	0.61
Original C index (95% CI)	0.82 (0.74–0.90)	0.86 (0.79–0.92)
Bias-corrected C index	0.80	0.82
Brier score	0.096	0.093
AIC	171.97	166.63
BIC	192.76	201.27

AIC, Akaike information criterion; BIC, Bayesian information criterion; CI, confidence interval, GOF, goodness of fit; mod_5V, model_5V; mod_9V, model_9V; PR curve, precision-recall curve; ROC auc, receiver operating characteristic curve and the area under the curve.

^aANOVA was performed to compare the value of χ² for GOF between the mod_5V and the mod_9V.

The Hosmer-Lemeshow test was also applied to both predictive models to assess the goodness of fit (GOF). The χ² values of mod_5V and mod_9V were 5.13 and 9.23, respectively, and no significant difference was found between the observed values and the predicted values in either model. The P-values were 0.7432 (mod_5V) and 0.3236 (mod_9V; Table 4). Subsequently, ANOVA was performed to compare the value of χ² for the GOF between mod_5V and mod_9V, yielding a P-value of 0.0097, which was lower than 0.01, indicating that there was a significant difference between the two models (Table 4).

The ROC curves for the two models are depicted in Figure 3. The area under the curve for mod_5V was 0.819 (95% confidence interval (CI), 0.742–0.897) as presented in Figure 3a, and for mod_9V, it was 0.857 (95% CI, 0.789–0.925) as shown in Figure 3b, with details in Table 4. The specificities of mod_5V and mod_9V were 81.82 and 77.78%, respectively (Table 4). The sensitivity of mod_5V was 73.68%, while that of mod_9V was 78.95%. Precision-recall (PR) curves are plotted in Figure 4a and b. The PR AUC values for mod_5V and mod_9V were 0.521 and 0.61, respectively, indicating that mod_9V had a higher ability to predict real positive cases among the predicted positive cases than mod_5V. The original and bias-corrected C indices of mod_5V were 0.82 and 0.80, respectively. In contrast, those of mod_9V were 0.86 and 0.82 indicating the coherence between the predicted results and the actual results achieved by mod_9V was better than those by mod_5V. The lower and upper values of the C index for mod_5V and mod_9V were 0.74 to 0.90 and 0.79 to 0.92, respectively (Table 4). The calibration curve of mod_9V was better fitted with the standard curve than that of mod_5V, and the Brier scores of mod_5V and mod_9V were 0.096 and 0.093, respectively, indicating that, as a measure of prediction accuracy, the Brier scores in mod_9V were lower than those of mod_5V, and mod_9V displayed higher prediction accuracy (Figure 5a and b, Table 4). AIC (Akaike information criterion), which is a standard to measure the goodness of statistical model fit, can weigh the complexity of the estimated model against the goodness of the data that the model fits. As well as AIC, BIC (Bayesian information criteria) is used for model selection and introduces a penalty term to avoid the model complexity with to large dimension. Generally, the smaller the value of the ACI and BIC, the better the performance of the model. The AIC and BIC values for mod_5V were 171.97 and 192.76, respectively. However, the AIC and BIC values for mod_9V were 166.63 and 201.27, respectively. This indicates that although the BIC value of mod_9V is higher than that of mod_5V showing that mod_9V contains more variables and is therefore more complex, the AIC value of mod_9V is lower than that of mod_5V revealing that mod_9V has a better GOF.

Figure 3. Receiver operator characteristic curves (ROC) were used to evaluate mod_5V and mod_9V. (a) ROC analysis for mod_5V was conducted using CRP, IL-2R, and percentages of CD3+, CD4+ and CD8+ T-cell lymphocyte subsets, to predict the probability of IVIG retreatment. The area under the curve (AUC) was 0.819 with a 95% confidence interval ranging from 0.742 to 0.897. The cut-off value was 0.184, with associated specificity and sensitivity of 0.818 and 0.737, respectively. (b) ROC analysis for mod_9V included NT-proBNP, CRP, sex, IL-2R, TNF-α, IL-6, and percentages of CD3+, CD4+ and CD8+ T-cell lymphocyte subsets to predict the probability of IVIG retreatment. The AUC was 0.857, with a 95% confidence interval from 0.789 to 0.925. The cut-off value was 0.154, with related specificity and sensitivity of 0.778 and 0.789, respectively.

Figure 4. Precision-recall curve (PR curve) analysis for the prediction model of mod_5V and mod_9V. (a) PR curve for mod_5V. The area under the PR curves for mod_5V was 0.521, demonstrating model accuracy across the five variables, including the percentage of CD3+ T cells, CD4+ T cells, CD8+ T cells, serum levels of IL-2R and CRP. The plot shows the percentages of cases that yielded positive predictions. ‘Recall’ (or sensitivity) represented the number of cases that were accurately predicted to require IVIG retreatment in children with KD who did receive IVIG retreatment. (b) PR curve for mod_9V. The area under the PR curve for mod_9V was 0.611, demonstrating the model accuracy across the nine variables, including the percentage of CD3+ T cells, CD4+ T cells, CD8+ T cells, sex and serum levels of IL-2R, IL-6, TNF-α, CRP and NT-proBNP. The plot shows the percentages of cases that yielded positive predictions. ‘Recall’ (or sensitivity) represented the number of cases that were accurately predicted to necessitate IVIG retreatment in children with KD who did receive IVIG retreatment.

Figure 5. Calibration curve analysis and Brier Score value for the prediction IVIG retreatment for mod_5V and mod_9V. (a) Calibration curve and Brier Score for mod_5V. The plot shows the mean-squared difference between the actual value of IVIG retreatment and the corresponding predicted probability, which is 0.096. (b) Calibration curve and Brier Score for Mod _9V. The plot shows the mean-squared difference between the actual value of IVIG retreatment and the corresponding predicted probability, which is 0.093.

Decision curve analysis was performed to evaluate the clinical usefulness of the two prediction models. As shown in Figure 6a, based on a continuum of potential thresholds for IVIG retreatment (x-axis) and the net benefit of using the two models to assess the risk of patients needing IVIG retreatment (y-axis), assuming that no patient will require IVIG retreatment, the decision curve analysis (DCA) of mod_9V provided a larger net benefit across the range of IVIG retreatment compared with that of mod_5V.

Figure 6. Decision curve analysis (DCA) and clinical impact curve (CIC) for mod_5V and mod_9V. (a) DCA for mod_5V and mod_9V. The solid red line represents DCA for mod_5V, whereas the solid green line represents DCA for mod_9V. The area under the green solid curve was greater than the area under the red solid curve over the entire range of IVIG retreatment probabilities, as assessed by DCA, where the net benefit was higher for mod_9V, indicating that mod_9V provides superior clinical usefulness compared with mod_5V. (b) CIC for mod_5V. Of 1000 children with Kawasaki disease (KD), the red solid curve represents the total number of patients deemed high risk at each risk threshold according to mod_5V. The blue dashed line shows the number of patients who received IVIG retreatment. (c) CIC for mod_9V. For 1000 children with KD, the green solid curve shows the total number of patients considered high risk at each risk threshold according to mod_9V. The blue dashed line shows the number of patients who received IVIG retreatment.

Figure 6b and c shows the clinical impact curves (CIC) for mod_5V and mod_9V, respectively, obtained based on Figure 6a. The two images show the estimated number of patients who would be deemed high risk for each risk threshold compared with that of patients who were actually re-treated with IVIG (true positives). As shown in Figure 6b and c, the red solid curve presents the number of patients requiring IVIG retreatment, as predicted by mod_5V, whereas the green solid curve represents the number of patients predicted by mod_9V. The blue curve of the dashed line shows the actual number of children with KD who received IVIG retreatment. As shown in Figure 6b and c, as the threshold high risk increases, the number of children with KD needed IVIG retreatment predicted by mod_9V was more consistent with the actual number of patients than that predicted by mod_5V.

To evaluate the improvement between mod_5V and mod_9V, NRI and IDI were calculated.²⁵ When mod_9V was compared with mod_5V, the NRI was 0.15, with a 95% CI of 0.01–0.30 and a P-value of 0.0039, which was < 0.01. The IDI was 0.07, with a 95% CI of 0.02–0.12 and a P-value of 0.0076, indicating a good improvement in mod_9V over mod_5V.

Validation of the prediction models in the DC and the VC

The confusion matrix is a 2 × 2 table that summarises model prediction results in machine learning, categorising samples from the data set based on their actual and predicted categories. The confusion matrices of the two models from the DC are shown in Figure 7a and b. Both prediction models were validated in the VC, and the results are shown in Figures 7c and d. When the two prediction models were applied to the VC, the recall, specificity, PPV, NPV and Kappa values were 0.8333, 0.8750, 0.6250, 0.9545 and 0.6296 for mod_5V, respectively (Table 5). While for mod_9V, the recall, specificity, PPV, NPV and Kappa values were 1, 0.8750, 0.6667, 1 and 0.737, respectively (Table 5). The Kappa values of both mod_5V and mod_9V on the DC were higher than 0.4, indicating medium stability of the models, whereas the value of both on VC was higher than 0.6, indicating very good stability.²⁶ The F₁ score is the harmonic average of the sensitivity and precision. In the DC, the F₁ scores of both mod_5V and mod_9V were higher than 0.53. The F₁ score of mod_5V reached 0.714 and that of mod_9V was 0.8, indicating that mod_9V was more generalizable to different databases than mod_5V.

Figure 7. Confusion matrix for mod_5V and mod_9V based on the DC and the VC. (a) Based on the DC, confusion matrix for mod_5V. (b) Based on the DC, confusion matrix for mod_9V. (c) Based on the VC, confusion matrix for mod_5V. (d) Based on the VC, confusion matrix for mod_9V.

Table 5 The performance of mod_5V and mod_9V for prediction of IVIG retreatment probability on the DC and the VC

	DC		VC
	mod_5V	mod_9V	mod_5V	mod_9V
TP	28	30	5	6
TN	162	154	21	21
FP	36	44	3	3
FN	10	8	1	0
Sensitivity/recall	73.68%	78.95%	83.33%	100%
Specificity	81.82%	77.78%	87.50%	87.50%
Accuracy (95% CI)	80.51% (74.87–85.36%)	77.97% (72.13–83.08%)	86.67% (69.28–96.24%)	90% (73.47–97.89%)
Precision/PPV	43.75%	40.54%	62.50%	66.67%
NPV	94.19%	95.06%	95.45%	100%
Kappa	0.43	0.41	0.63	0.74
F1	0.55	0.54	0.71	0.80

DC, development cohort data set; FN, false negative; FP, false positive; NPV, negative predictive value; PPV, positive predictive value; TN, true negative; TP, true positive; VC, validation cohort data set.

Discussion

Numerous diagnostic risk scores have been developed to predict IVIG non-responders in KD. However, Tremoulet et al. found that the sensitivity of the Egami score for predicting IVIG resistance in Asian children was only 33.3%.²⁷ Some scholars in China used the Egami and Sano scores to evaluate children with KD in Shanghai and predicted IVIG non-responders. The area under the ROC curve of the Egami score was 0.656 and its sensitivity and specificity were 8.2% and 93.5%, respectively. The area under the ROC curve of the Kobayashi score was 0.650 and the sensitivity and specificity were 61.6% and 66.5%, respectively, which were lower than those of the newly established models. It can be seen that the evaluation and application of the model still vary greatly in different regions, even in different hospitals in the same region. Combined with the testing system of the hospital where the patients live in, the prediction of non-responders to IVIG is more valuable. Therefore, the present study retrospectively screened children with KD and constructed models for the prediction of IVIG retreatment, combined with routine diagnosis and treatment and integrated variables representing demographic information, clinical indicators, inflammatory responses and related immunology indices.

In the DC and the VC, the percentages of KD patients with IVIG retreatment were respectively 16.1% and 20%, which is in accordance with reported data.^1–4 The laboratory indicators, the level of cytokines, and the tests of lymphocyte subsets except for the CD4⁺ T cells, showed a similar trend between the IVIG retreatment group and the IVIG single-dose group in both the DC and the VC, although there was no significant differences in the VC, indicating that the DC and VC could well present the characteristics of children with KD. Although the percentage of CD4⁺ T cells revealed a different trend in the VC than in the DC, there was no obvious difference.

Several studies have reported that KD tends to occur more frequently in boys than in girls, although there was no obvious difference in the number observed between boys and girls with KD.^1,3 However, in previous studies, the number of enrolled patients was more in boys.^11,28 Therefore, it is important to consider the role of sex. The present study revealed that the female sex had a significant impact on the probability of IVIG retreatment. It has been reported that higher IgG levels in patients with KD are associated with more effective immunomodulatory effects and better clinical outcomes.²⁹ We speculated that IgG levels may be lower in girls than in boys; therefore, the effect of IVIG on girls with KD was lower than that on boys. This may be the reason that the female sex had an increased risk of IVIG retreatment. However, the mechanism requires further investigation.

The notable elevation of cytokines levels in KD and the novel therapeutic prospects they present have garnered increasing attention, including the study of IL-1. Two case reports described successful treatment of KD in children using IL-1 antibodies.^30,31 The conditions of the two patients were rare and presented with extreme severity. One was a neonate and the other was a relapsed KD occurring giant CAA. In our data, we did not encounter any KD patients with a similar profile. No difference in IL-1β levels was observed between the patients with a single dose of IVIG and IVIG retreatment, so IL-1β was automatically removed by R language when the model was built.

IL-6 is typically synthesised during the initial stages of inflammation³² and can either promote inflammatory responses through IL-6-trans-signalling or mediate anti-inflammatory reactions by binding to the IL-6 receptor.^33,34 The degree of alterations in IL-6 levels in an animal model of KD vasculitis was similar to that observed in children with KD, which increased during the acute phase of KD, followed by a rapid reduction during the course of convalescence. IL-6 fluctuates sensitively along with the inflammation development or inhibition by therapy in the acute phase of KD.^13,35,36 Therefore, changes in IL-6 levels may be considered sensitive indicators of the extent of inflammation. In the initial phase of KD, the significantly higher levels of IL-6 in the group of children who received IVIG retreatment than in the group single dose of IVIG indicated a heightened KD activity and a more intensive inflammatory response in these patients suggesting enhanced IL-6 levels is significantly associated with IVIG retreatment. A small, prospective pilot study of tocilizumab has been conducted for IVIG-resistant KD including totally four children.³⁷ All patients' clinical and laboratory measures rapidly improved after the start of tocilizumab therapy. However, two patients still developed progressive giant coronary artery aneurysms suggesting the clinical application of tocilizumab in KD needs more experience.

As a component of the acute-phase immune response, TNF-α is trace or undetectable in healthy individuals, but its levels increase during the inflammatory state.³⁸ Serum TNF-α levels have been reported to be notably elevated in children with KD.¹² However, the efficacy of the anti-TNF-α antibody infliximab in preventing the development of CAA in the pathogenesis of KD is controversial.^4,14,39 The presented study revealed that serum TNF-α levels were elevated in both the single-dose IVIG group and the retreatment group, with a more pronounced increase observed in the latter. However, mod_9V showed that along with the increase in serum TNF-α, the IVIG retreatment probability decreases, indicating that the ability of body to produce TNF-α early, rather than the value of TNF-α, is the significant effector of IVIG retreatment, which reflects the body's ability to initiate inflammatory reaction of the immune system. As reports have identified that TNF-α is crucial for effective and controlled inflammatory responses when it is produced in the right place, time and context.⁴⁰ Recently, it was proposed that TNF also has the function of preserving body homeostasis during disease.⁴¹ Therefore, a high TNF-α level in the initial phase of KD is a protective factor for patients to avoid receiving another IVIG dose.

During the first 10 days of KD, myocarditis has been identified in the entire cardiac tissue,⁴² which has been deemed to be one of the earliest features of KD.⁴³ A previous study revealed that BNP release under certain cardiac stress situations precedes the clinical manifestation of myocardial ischemia or infarction.⁴⁴ The present study showed that NT-proBNP levels in patients of IVIG retreatment increased significantly higher than those in patients who received only one dose of IVIG in the early phase of KD. Mod_9V showed that increased serum NT-proBNP was associated with the avoidance of IVIG retreatment. The OR of NT-proBNP was < 1, and the upper limit of 95% CI reached 1, suggesting that the level of serum NT-proBNP had remained a significant impact on IVIG retreatment. We presume that as well as TNF-α, the early increase in NT-proBNP reflects a kind of stress capacity from the cardiac tissue in children with KD other than the intensity of inflammation for IVIG retreatment prediction. Therefore, the body's ability to produce high level of NT-proBNP in the initial phase of KD is a protective factor against IVIG retreatment.

Decreased T-cell counts, CD3-induced T-cell proliferation and increased levels of pro-inflammatory cytokines were observed during the acute stages of KD.⁴⁵ The present study revealed that during the initial phase of KD and before the first administration of IVIG, a significant decrease in the percentage of CD3⁺, CD4⁺ and CD8⁺ T cells and an increase in the percentage of CD19⁺ B cells were observed in patients from IVIG retreatment group compared with those in the single-dose IVIG group. These results are consistent with those of recent reports from other research groups in Shanghai.^11,28,46 Changes in the distribution of lymphocyte subsets during the acute phase of KD before treatment indicated that patients who underwent IVIG retreatment exhibited higher humoral immunity activity, stronger suppression of cellular immunity and more severe inflammatory responses than patients with single-dose IVIG. Mod_9V showed the patient with lower percentages of CD4⁺ and CD8⁺ T cells and the higher the percentage of CD3⁺ T cells were more likely to receive IVIG retreatment suggesting the normal cellular immunity activity is necessary for the KD inflammation suppression. Logically, the model should include B cell subsets that mainly mediate humoral immunity, while the fact is that comparison between the two groups of patients showed, in the patients with IVIG retreatment, in addition to the obvious increase in the percentage of CD19⁺ B cells, the percentage of CD3⁺, CD4⁺ and CD8⁺ T cells were all significantly reduced, although CD4⁺ T cells contain Th₂ cells that can mediate humoral immunity, the changing trend of CD4⁺ and CD3⁺ T cells still significantly decreased, which suggests that in patients with IVIG retreatment, the degree of suppression of cellular immunity may be far greater than the degree of hyperactivity of humoral immunity. Therefore, when different lymphocyte subsets are analysed together, The R program automatically prioritized CD3⁺, CD4⁺ and CD8⁺ T-cell subsets, which exhibited more pronounced changes, while excluding the CD19⁺ B cell subsets. A number of previous studies, including ours, have observed that the levels of IL-2R as a cytokine receptor change more sensitively, rising or falling rapidly along with disease changes, compared with the counts of lymphocyte subsets.^11,15,47 The level of IL-2R is not only associated closely with alterations in the counts of CD8⁺ T cells and CD19⁺ B cells but also fluctuates rapidly along with the intensity of inflammation before and after IVIG treatment. A decrease in IL-2R levels in KD children following effective IVIG treatment was inversely associated with the pre-treatment percentage of CD8⁺T cells, indicating a significant suppression of cellular immunity prior to treatment.¹¹ In the presented study, children with IVIG retreatment often in a higher level of the IL-2R by more than 2000 units compared with children with IVIG one dose (Table 2). The OR value of IL-2R was 1.005 in mod_9V, that is, the risk of IVIG retreatment will increase by 1.005 times for every 1 unit increase in IL-2R (Table 3) suggesting the increase of IL-2R can be a sign of increased cellular immune suppression in children with KD and children with a high level of IL-2R are more likely to undergo IVIG retreatment. Therefore, relatively high levels of CD4⁺ T cells and CD8⁺ T cells reduced the risk for KD patients to receive IVIG retreatment. While the elevated percentage of CD3⁺ T cells and high level of IL-2R increases the risk of IVIG retreatment.

Of its sensitivity, CRP has received most research attention as a clinical marker of inflammation in cardiovascular disease and even in asymptomatic individuals.⁴⁸ In the present study, CRP levels showed a clear difference between patients re-treated with IVIG and those who received only one dose. In both prediction models, higher CRP was found to be significantly associated with IVIG retreatment.

The suite of evaluation and validation metrics for both predictive models demonstrated a feasible overall accuracy in predicting the need for IVIG retreatment in children with KD. However, there were significant differences between mod_5V and mod_9V by comparison, and the datasets from both DC and VC presented mod_9V as the superior model over mod_5V. Mod_9V may be limited by the number of IVIG retreatment cases in the DC data such that the nine variables may seem to be redundant. However, higher values of the ROC AUC, sensitivity (TPR, true predictive rate), C indices (both the original C index and the bias-corrected version), PPV, NPV and lower AIC scores, were achieved by mod_9V than those by mod_5V. These results suggest that mod_9V has greater predictive power in addition to a higher concordance between the actual and predicted probabilities. What is more, the PR curve was used as a complementary evaluation to assess the difference between the predicted cases of IVIG retreatment and actual cases by recall; the larger the area under the PR curve, the better the model performance.^49,50 The results revealed that compared with mod_5V, mod_9V had a larger AUC for PR in the DC, where the recall value of mod_9V even reached 1 in the VC, outperforming mod_5V. Brier score curves were plotted to evaluate the overall prediction accuracies of the models. Lower values of the mean-squared difference were found in mod_9V than in mod_5V, suggesting a superior Brier score for the effect calibration for mod_9V. DCA and clinical impact curves also showed that mod_9V predicted IVIG retreatment cases with superior clinical benefits and utility compared with mod_5V. The positive NRI and IDI values indicate an improvement in mod_9V over mod_5V. In the confusion matrix, the values of PPV, kappa, TPR (recall) and NPV of mod_9V were all higher than those of mod_5V in the VC, verifying the superiority of mod_9V over mod_5V. In particular, the kappa values of mod_9V in DC and VC were 0.41 and 0.737, achieving medium consistency and good consistency, respectively and reflecting the stability of the model. At the same time, the consistency index C index of mod_9V reached 0.8567, whereas 0.8192 of mod_5V indicated better prediction results mod_9V would obtain than mod_5V under the premise of maintaining good stability.

In a clinical setting, the top priority when treating patients with KD is to determine the optimal therapeutic regimen for IVIG. It is of higher importance to predict whether IVIG needs to be used once more in the future according to individual clinical characteristics, in addition to confirming whether there is sufficient pharmacological effect of IVIG on children with KD. The rationale for this is to promptly suppress inflammation and alleviate cardiovascular lesions as early as possible. In fact, most KD children experienced complete remission of symptoms after the second dose of IVIG. Therefore, the present study set the IVIG retreatment as the prediction endpoint based on the actual condition, rather than ‘IVIG resistance’. The nine variables contained in mod_9V study contain a combination of specific inflammatory indicators signing the inflammatory reaction process, demographic information and reflecting accurately the patient's individual clinical situation. Although there are unequal tightness correlations among the nine factors, relevance based on sign signal paths is originally existed among various biological indices. The distribution of independent variables indeed revealed that the IVIG retreatment is the result mainly of inflammation and immunoreaction.

Therefore, between the 5-variable model and the 9-variable model, we choose and recommend the latter between the five- and nine-variable models.

To the best of our knowledge, the present study was the first study to use logistic regression to construct a predictive model for IVIG retreatment in combination with lymphocyte subsets, inflammatory cytokine levels, laboratory indicators and demographic characteristics, which is a strength of this study. Cytokines and subgroups of lymphocytes delineate the substantive feature of individual child with KD and nowadays have been routinely detected and obtained in most children's cardiovascular departments in medical centres once the patient was diagnosed with KD. Mod_9V was mainly constituted of these variables and could achieve good results according to the individual patient with KD.

In the present study, two different and completely uncorrelated datasets DC and VC were constructed. DC was used to construct the prediction models, and VC was used to verify the models. The validation indexes guarantee the good performance of the model, which is another strength of this study.

The present study had some limitations. All children with KD in the present study were treated at the Xinhua Hospital affiliated to Shanghai Jiao Tong University School of Medicine, which is a single-center study. Although Xinhua Hospital is a medical center located in East China, the patients who come for treatment are from all over the country, but primary from the eastern region, which may present some geographical limitations in terms of patients' origins. The number of cases in the present study was also limited, especially cases of IVIG retreatment. One possible reason for this is that Chinese parents tend to pay more attention to their children's health conditions. Upon noticing symptoms such as fever, rash or a poor spirit, parents or guardians may seek medical attention for the child promptly, providing an opportunity for early treatment during the initial stages of KD. Timely therapy objectively suppresses the inflammation development and improves the prognosis of KD while also reducing the incidence of retreatment with IVIG, which could be a contributing factor to the lower number of children requiring IVIG retreatment in this study. In addition, the newly established electronic medical record system in our hospital was not fully functional during the early stages of use, resulting in incomplete medical data for some patients. Therefore, while 796 patients were screened, the actual number of patients included in our study, with complete and compliant medical data, was 266.

However, we conducted a series of model evaluations and a validation dataset to reduce the possible bias because of the limited number of cases. We compared the predictive models with an empty model and calculated the maximum likelihood value, which showed that the models had significant improvements in the null hypothesis.⁵¹ The GOF of the model was tested and found that there was good agreement between the actual and predicted values of mod_9V, including the χ² value of the GOF, C index and AIC value. In addition, we evaluated the stability of the model by the value of kappa. The agreements of mod_5V and mod_9V in the DC and VC were moderate and good as the kappa coefficient was higher than the value of 0.6 and 0.8, respectively. The AIC value of mod_9V is smaller than that of mod_5V, and the other values of the evaluation indices of mod_9V are better than those of mod_5V, including sensitivity (recall), C index and NPV. Thus, we believe that the results of mod_9V are more robust and reliable than those of mod_5V.

Conclusion

During the initial phase of KD in children, lymphocyte subsets and cytokines provide specific inflammatory profiles and immune features and play a dominant role in characterising individual KD. Based on these variables, mod_9V was constructed for IVIG retreatment prediction in real-life cohort and proved superior to mod_5V. This allowed the achievement of satisfactory predictive power and provided a novel strategy for identifying children with KD most in need of IVIG retreatment and for determining an early therapeutic regimen.

Methods Data acquisition

Data were obtained from Xinhua Hospital Affiliated with Shanghai Jiao Tong University School of Medicine from 1 April 2017 to 31 May 2022. In total, the medical records of 796 children with KD were screened, and patients whose medical records met the requirements of the present study were enrolled consecutively. Children with KD were divided into two data sets in chronological order: the development cohort data set (DC) and the validation cohort data set (VC). From 1 April 2017 to 28 December 2021, data sets of 236 patients were included to develop the prediction model by setting up the DC. In contrast, from 29 December 2021 to 31 May 2022, VC was established to validate the prediction model using data from the medical records of 30 children with KD. As a result of the relatively small number of positive cases (patients who received IVIG retreatment), cases were included in this study along with the analysis, in chronological order and the later included cases were used for validation (Figure 8).

Figure 8. The flowchart of inclusion and study design. A total of 796 cases were screened from 1 April 2017 to 31 May 2022, enrolled continually in chronological order when meeting the requirement of this study, and 266 cases used to construct the development cohort data set (DC). From 29 December 2021 to 31 May 2022, the data enrolled were used to construct the validation cohort dataset (VC), 30 cases in total. KD, Kawasaki disease.

Infants and children aged between 4 months and 17 years who met the criteria for KD as stated by the Guidelines from the Scientific Statement for Health Professionals from the American Heart Association (AHA) were included in the study.¹ The detailed criteria for the diagnosis of KD were as follows: (1) The presence of fever for > 5 days and the presence of greater than or equal to four of the five principal clinical features of KD, including extremity changes such as erythema of the palms and soles, rash, conjunctivitis, oral changes and cervical lymphadenopathy. (2) In practice, using the abundance of experience and other typical symptoms exhibited by patients, clinicians can occasionally make an earlier diagnosis of KD even when the duration of the child's fever is only 3–4 days. The prognostic outcomes of all patients included in this study were categorized as either improved or recovered.

The exclusion criteria for cases were as follows: (1) children with incomplete clinical medical or laboratory indicator data; (2) those with serious cardiovascular, hepatic or renal diseases; primary disease associated with tumors; hematological diseases; congenital malformations; genetic metabolic diseases; children with primary myocarditis or other primary diseases of major organs and (3) children who had previously been diagnosed with KD and were readmitted due to relapse, requiring retreatment.

Once the diagnosis of KD was confirmed, all children were monitored and treated according to the standard routine regimen established in the American Heart AHA.¹ Child with KD received the initial therapy, which included one dose of IVIG (2 g kg⁻¹; Shanghai RAAS Blood Products Co., Ltd, Shanghai, China) and a moderate oral dose of aspirin (30–50 mg kg⁻¹ per day; Bayer HealthCare Manufacturing, Garbagnate Milanese, Italy), until the child's fever resolved. The subsequent regimen was a low oral dose of aspirin (3–5 mg kg⁻¹ per day for ≥ 2 months). Any patient whose fever persisted for ≥ 36 h or who had recrudescent fever for ≥ 48 h after the first dose of IVIG would receive a second IVIG infusion at 2 g kg⁻¹ as a rescue treatment method until the fever resolved. The temperature of the children was measured at all the time points (three times a day and once every 8 h), mainly by the axillary route, including the oral and rectal routes. Temperatures of > 37.5 °C in the axillary and 38 °C according to oral measurements were defined as fever. Venous blood samples were collected before primary treatment according to the diagnosis and treatment rules established by the Department of Paediatric Cardiology in our hospital. Routine blood examinations, C-reactive protein (CRP), erythrocyte sedimentation rate (ESR), platelet count (PLT), N-terminal precursor brain natriuretic peptide (NT-proBNP), serum electrolytes including serum sodium, inflammatory cytokines, subsets of lymphocytes, laboratory indicators including indicators of the liver and kidney and cardiovascular ultrasounds were performed. A total of 48 h after the primary IVIG treatment, the necessary indices, including cytokine levels and echocardiography, were detected again. Relevant records of the aforementioned items were extracted from the medical records of children with KD, which were used to establish DC and VC.

In the present retrospective study, all research protocols were approved by the Ethics Committee of Xinhua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, in accordance with the Declaration of Helsinki. This study conformed to the exemption from informed consent standards because of the use of general patient information, with the medical history and laboratory examination data being obtained from the medical records and database. We made an application for exemption from informed consent and obtained approval (Approval No. XHEC-D-2022-180). Written consent was not required for this study because of its retrospective nature.

Laboratory tests

The differential count of peripheral blood leucocytes was determined using a fully automatic haematology analyser. Inflammatory cytokines were detected using a solid-phase enzyme-labelled chemiluminescent immunoassay (Siemens Healthcare Diagnostics Products Limited, Gwynedd, UK). Interleukin-2 receptor (IL-2R) is a cytokine receptor; however, it is classified as a cytokine by the assay method of Siemens AG and in the present paper the IL-2R was also classified as a cytokine for conciseness of expression. NT-proBNP was detected using an NT-proBNP kit with electrochemiluminescence (Roche Diagnostics GmbH, Mannheim, Germany). Flow cytometry was performed to detect the lymphocyte subgroup counts (Becton, Dickinson and Company, BD Biosciences, California, USA). For individual patients, the absolute values of flow cytometry compared with the normal values will be more meaningful. However, given that we studied children of varying ages, we chose to use percentages to avoid this important difference because of the inconsistencies in lymphocyte maturity in children with KD of different ages. What is more, using a percentage of lymphocytes compared to an absolute count of lymphocytes makes the value less variable and more reliable, avoiding fluctuations and biases.

Establishment and verification of IVIG retreatment prediction models

In the DC, a total of 23 variables (demographic information, clinical manifestations and laboratory indicators) were assessed, and univariable logistic regression analysis was performed on each variable. Variables with P < 0.05, were then subjected to multivariate logistic regression by forward selection and backward elimination. In the final model, risk factors were created with P < 0.1. Cut-off points were selected based on the receiver operating characteristic (ROC) curve and outcome, which is located in the highest Youden index value calculated according to sensitivity and specificity. Another prediction model was constructed by multivariable logistic regression analysis including a total of 23 variables using either backward stepwise selection or forward and backward elimination. Variables with P ≥ 0.2 were removed from this analysis. Two prediction models were obtained and presented as nomograms, which are based on multivariable regression analysis. Multiple prediction indicators are integrated and then scaled line segments are drawn on the same plane in a certain proportion. It can be used to express the relationship between the variables in the prediction model. Significant variables were expressed as odds ratios (OR) and 95% confidence intervals (CI).⁸ The predictive ability of the models was assessed by the ROC curve and area under the curve (AUC). Precision-recall (PR) curves were plotted before calculating the associated AUC. The discriminative ability of the nomograms was measured using a concordance index (C index). Discrimination between the actual and predicted values was evaluated using a calibration curve followed by the Hosmer-Lemeshow test presented with goodness of fit. The Brier score was also used to assess the overall predictive performance using the calculated value, which was defined as the mean-squared difference between the predicted and actual values.⁵² In addition, the AIC and BIC were applied for model comparison. Clinical usefulness and net benefit were estimated using a DCA and a CIC.⁵³ To assess the difference between the two models in terms of usefulness, a net reclassification index and an integrated discrimination improvement were performed.^54,55 The data were also classified through a confusion matrix to derive core metrics, including the positive prediction value (PPV) and negative prediction value (NPV).⁵⁶ All P-values were two-tailed, and statistical significance was set at P < 0.05, unless stated otherwise.

Both prediction models were assessed and validated in the DC and VC by a confusion matrix using the ‘caret’ package in R to identify predictive power.

All data were analysed using SPSS software (version 25.0, IBM Corp., USA). Modelling and plotting were performed in R language (version 4.1.2, USA), and the packages used in R are listed in Table 6.

Table 6 List of R packages

Purposes	Functions	Packages
Data input	‘reade_excel’	readxl
Modelling, nomogram, Calibration curve Goodness of fit	‘glm’, ‘lm’, ‘lrm’, ‘nomogram’, ‘Calibrate’	‘rms’, ‘ResourceSelection’
Precision-recall curve plot	‘aupr’	‘modEvA’
Stepwise regression analysis	‘stepAIC’	‘MASS’
Hosmer and Lemeshow test	‘hoslem.test’, ‘Mann–Whitney non-parametric test’	‘rcompanion’
Multivariable logistic regression	‘vglm’	‘VGAM’
Receiver operating curve plots, AUC values and test	‘roc’	‘pROC’
Decision curve analysis plots, clinical impact curve plot	‘decision curve’	‘rmda’
Diagnosis values	‘confusionMatrix’	‘caret’
Models comparison, improvement estimation	‘NRI’	‘nricens’
Calibration, calculation of IDI	‘plotCalibration’, ‘IDI’, ‘reclassification’	‘PredictABEL’

AUC, area under the curve; IDI, integrated discrimination improvement; NRI, net reclassification index.

Statistical analysis

Children with KD were divided into the IVIG single-dose group and the IVIG retreatment group. Continuous variables are expressed as medians with interquartile ranges (M [Q1, Q3]) or as mean ± standard deviation (mean ± SD). Categorical variables are summarised and presented as frequencies (N) and percentages (%). The baseline characteristics of the patients, including the clinical, laboratory and echocardiographic results, were recorded, described and compared between the single-dose IVIG and IVIG retreatment groups. For continuous variables, if the data followed a Gaussian distribution, the unpaired t-test was used; otherwise, the Mann–Whitney U-test was applied to assess the statistical significance between the two groups. For categorical data, the χ² test or Fisher exact test was used.

In the DC, the following 23 variables were analysed: demographic factors of age, sex and weight; clinical characteristics of the type of KD (complete or incomplete KD); days of fever at diagnosis; laboratory indices of the neutrophil count, platelets, CRP, ESR, NT-proBNP and result of the cardiovascular ultrasound; the levels of a range of cytokines, including IL-2R, IL-10, IL-6, IL-1β, IL-8 and TNF-α and the percentages of lymphocyte subgroups, including cluster of differentiation 19 (CD19)-positive CD19⁺ B cells, CD16⁺CD56⁺ NK cells, CD8⁺ T cells, CD4⁺ T cells and CD3⁺ T cells. In the VC, according to the results from the DC, all variables included in the DC were also compared between the two groups, except for weight, type of KD, number of days with fever at diagnosis and result of the cardiovascular ultrasound, which have been identified as not having a significant impact on IVIG retreatment in the DC. All variables were analysed to test whether they are significant independent risk factors, with IVIG retreatment designated as the dependent variable.

Acknowledgments

We thank Professor Liang Zhu for overall discussions. This study was supported by the Xinhua Hospital Investigator-initiated Clinical Research Project (No. XHEC-C-2019-115-2), International Cooperation Project of the National Natural Science Foundation of China (No. 81720108003), National Natural Science Foundation of China (No. 82071936), National Key Research and Development Program of China (No. 2018YFC1002403), Collaborative Innovation Program of the Shanghai Municipal Health Commission (No. 2020CXJQ01) and Translational Medicine Cross-Research Fund of Shanghai Jiao Tong University (No. YG2021ZD11).

Author contributions

Chun Zhang: Conceptualization; data curation; formal analysis; funding acquisition; investigation; methodology; project administration; resources; software; supervision; validation; visualization; writing – original draft; writing – review and editing. Sun Chen: Conceptualization; data curation; funding acquisition; investigation; project administration; resources; supervision; validation; visualization. Yan Bian: Data curation; investigation; project administration; validation. Xiaohua Qian: Formal analysis; methodology; software. Yurui Liu: Data curation; formal analysis; investigation; software; validation. Liqing Zhao: Data curation; investigation. Jia Shen: Data curation; investigation. Jiani Song: Investigation; validation. Peng Zhang: Methodology. Lun Chen: Validation. Limin Jiang: Methodology.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could influence the work reported in this study.

Data availability statement

The data that support the findings of this study are available from the Xinhua Hospital, but restrictions apply to the availability of these data, which are used under licence for the current study, and so are not publicly available. However, data are available from the authors upon reasonable request and with the permission of the Xinhua Hospital.

Ethics approval and consent to participate

This study was approved by the Ethics Committee of Xinhua Hospital (Approval No. XHEC-D-2022-180). In particular, the detection of cytokines and subgroups of lymphocytes are routinely performed in the child cardiovascular department in Xinhua Hospital according to the KD Diagnosis and treatment routine, which was permitted and signed informed consent by children's guardians when they are hospitalised.