Clustering and climate associations of Kawasaki Disease in San Diego County suggest environmental triggers

Abstract

Kawasaki Disease (KD) is the most common cause of pediatric acquired heart disease, but its etiology remains unknown. We examined 1164 cases of KD treated at a regional children's hospital in San Diego over a period of 15 years and uncovered novel structure to disease incidence. KD cases showed a well-defined seasonal variability, but also clustered temporally at much shorter time scales (days to weeks), and spatiotemporally on time scales of up to 10 days and spatial scales of 10-100 km. Temporal clusters of KD cases were associated with strongly significant regional-scale air temperature anomalies and consistent larger-scale atmospheric circulation patterns. Gene expression analysis further revealed a natural partitioning of KD patients into distinct groups based on their gene expression pattern, and that the different groups were associated with certain clinical characteristics that also exhibit temporal autocorrelation. Our data suggest that one or more environmental triggers exist, and that episodic exposures are modulated at least in part by regional weather conditions. We propose that characterization of the environmental factors that trigger KD in genetically susceptible children should focus on aerosols inhaled by patients who share common disease characteristics.

Figure 1

(A) Location of primary residence for KD cases treated at Rady Children’s Hospital San Diego, a regional pediatric medical center, between 1 January 2002 and 15 April 2017. Background imagery from the Wolfram Mathematica Knowledgebase. (B) Average number of cases per day for each month of the year in the San Diego region. Error bars show the standard deviations over the sample period. (C) Number of cases for each day (black dots) over the study period; no points are plotted for zero cases. Blue vertical bars show temporal clusters, defined as 7-day windows including onsets of 4 or more KD cases (see Materials and Methods, Table 1).

Figure 2

Robustness of temporal and spatiotemporal clustering of KD Cases in San Diego County. (A) The cluster size distribution of the construction shown in Fig. 1C. The number of clusters of with more than S cases are plotted against S. The result is compared with the corresponding curves obtained from a Monte Carlo simulation of synthetic time series that have the same average incidence and average seasonality as the San Diego KD time series, but for which there is no temporal structure apart from a regular seasonal pattern. The error bars are ±1 standard deviation. (B) Excess K statistic for spatiotemporal distribution of KD cases. Red indicates significant spatiotemporal clustering with p < 0.05, and yellow indicates 0.05 <p< 0.15.

Figure 3

Composites of (A) 700 hPa geopotential height (700 mb height) and (B) daily maximum temperature for a sequence of days before and after each KD case within a given group of cases, from 4 and 5 days before (DOO - 5) to 4 and 5 days after (DOO + 5) date of onset (DOO). The groups of cases shown are: (2nd row from top of each panel) all KD cases; KD cases not falling into a temporal cluster (3rd row from top); and KD cases falling into a temporal cluster (bottom row). For comparison, composites associated with a random sampling of days having no KD cases is shown in the upper row in each panel. In both sets of figures, stippling indicates statistical significance of an anomaly composite in a given location. On the bottom (local) panel, the blue dots in the date of onset figure (center column) show where the KD cases occurred. Here we show only cases (and random, non-KD DOO equivalent days) from the main KD season, between December and April (DJFMA). Equivalent plots for other seasons and other variables (U-wind, V-wind, Minimum Daily Temperature, and Precipitation) are shown in the SI. The high pressure anomalies in the global figures for Cluster KD Case Days are consistent with reduced wind speeds across the Pacific and the high temperatures in the San Diego region.

Figure 4

Unique grouping of gene expression patterns in KD patients. (Left) Heat map of the top 100 genes that were significantly differentially expressed between acute and convalescent phases in KD patients who were not in clusters, adjusted for illness day of acute blood sample collected on days 2–10 after fever onset. Rows represent 100 genes and columns represent 118 acute phase patient samples. Each cell value represents the intensity of mRNA expression in log-scale. Two patient phenotypes, age and GGT values, are displayed at the top. Annotation at top: Horizontal line is median age (2.8 yrs.), High GGT = > 3 SD above upper limit of normal for age reference range (Right). Demographic and clinical characteristics associated with 118 KD patients grouped based on unsupervised clustering of gene expression patterns. Clinical laboratory data were obtained before treatment with IVIG and at the same phlebotomy as the acute samples for gene expression. Values are median and IQR except those marked with *, which are ratios. Comparisons by Mann-Whitney test for continuous variables, and a χ² test for proportions. (A larger version of the heat map is included in the Supplementary Information).

Figure 5

(A) Shows the median age in a 10-day window conditioned on the age of a reference patient with onset immediately prior to the beginning of the time interval. Fine black dots show all cases, blue circles show median values for a given bin. The p-value indicates the significance of a non-zero slope in the linear regression. (B) As in (A) but for erythrocyte sedimentation rate (ESR). (C) As in (A,B), but for absolute neutrophil count (ANC). R² values for the three panels are 0.55, 0.29, and 0.47, respectively.