Tracing Tokyo’s air sources to identify Kawasaki Disease’s etiological triggers

This document intends to expand, provide context and reproducible analysis code on the study presented in a poster for the 20th International Vasculitis and ANCA Workshop, which takes place in Dublin, Ireland, from the 3rd to the 6th of April 2022.


Kawasaki Disease (KD) is a systemic vasculitis that mainly affects children younger than 5 years old. Although KD cases have been registered in over 60 countries across several continents, its incidence is highest in East Easia, particularly in Japan, where the highest annual incidence rate was recorded in 2018: 359 per 100 000 children aged 0-4 years [Ae et al., 2020]. After more than five decades since its discovery and active research, the etiology of KD is yet to be elucidated. Recent studies have analyzed the association between KD and diverse environmental factors, with some advances pointing towards a relevant role of the atmospheric transport of a wind-borne agent triggering the disease [Rodó et al., 2011, Rodó et al., 2014]. The specific nature of this agent(s) is still unknown, with biological (bacterial, fungal or viral) and chemical (pollution) elements being the strongest candidates to initiate the autoimmune cascade which leads to the disease.

One of the most striking features of the epidemiological dynamics of Kawasaki Disease is the defined seasonal structure across multiple countries, with broad coherence in fluctuations of cases across the Northern Hemisphere extra-tropical latitudes [Burns et al., 2013].

In the case of Japan, the number of admissions has been recorded in the epidemiological records since 1970, and show the following pattern:

KD Records Japan

Three main features are notable:

  • The epidemic peaks of 1979, 1982 and 1986 (and to a lesser degree, 1984).

  • The increasing trend starting around 1995 until today.

  • The marked seasonality from 2000 onwards.

All of this prompts several questions towards the nature of the drivers of these changes in incidence for a disease that is non-communicable in nature.

In this work, we focus on the period from 2011 to 2018 for the prefecture of Tokyo, and try to associate the seasonal maxima and minima to atmospheric transport patterns.


Kawasaki Disease


We collected data on date of admission and reported date of symptoms onset for a total of 13970 KD cases in hospitals belonging to the Tokyo prefecture from 2011 to 2018. This data was originally collected by the (22nd to 25th) nationwide epidemiological surveys of Kawasaki disease in Japan [Ae et al., 2020, Makino et al., 2015, Makino et al., 2018, Makino et al., 2019].

Computation of KD maxima and minima

The data was aggregated to generate daily counts of hospital admissions and registered onsets.

Daily counts

The daily counts were then aggregated by week to obtain a weekly estimate of case onsets per day (that is, the daily average for each consecutive week in a year).

To generate the yearly maxima and minima, we selected the 5 weeks with the most daily cases and the 5 weeks with the least daily cases for each year, respectively, and are shown in the following figure:

Weekly counts

Trajectory Generation

To model air mass back-trajectories, we used the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model version 5.2 [Stein et al., 2015], which we operated programmatically via the Python package PySPLIT [Warner, 2018] to generate a high amount of trajectories.

For each day of the period from the 1st of January 2011 to the 31st of December 2018, we generated 4 backtrajectories up to the previous 96h starting at 00:00, 6:00, 12:00 and 18:00, with the initial point 10 meters over sea surface in central Tokyo (139.65E, 35.68N).

Below, an extract showing the individual trajectories generated at different points in time:

Weekly counts

Meteorology data

To model the air trajectories, HYSPLIT requires a set of gridded meteorological data at all pressure levels. In this case, we downloaded the GDAS 1x1°, 3 hour resolution dataset for every week from December 2010 (to be able to generate backtrajectories starting on January 2011) to December 2018. These data can directly be accessed through NOAA’S ARL FTP Server.

Differential Trajectory Analysis

As a first step, we generated a grid of frequency of intersection between lat-lon grid cells and the trajectories associated to KD maxima and KD minima dates, which are a total of 1120 trajectories associated to each group.

We first generate a grid of 200 by 200 cellsvbased on the boundaries of the generated trajectories, which ends up looking like a 0.8x0.8° resolution grid like the one displayed in the following image:

Map Grid

Then, we count the number of intersections of the trajectories of each KD group for each of the grid cells, which allows us to visualize longer term or aggregated patterns, like in this figure that displays the different sources associated to KD maxima and minima per calendar year:

Trajectories per year

By comparing, cell-wise, the ratio between KD maxima and KD minima intersections, we can obtain an overview of the source areas overrepresented in association with the timings of either phenomena.

To generate a summarized figure, we use the Log2 transformation of the ratios to generate the image, and colour each cell grid according to this. The transformation defined as:

\[R = \text{Log}_2\cfrac{\text{Intersections}_{max}}{\text{Intersections}_{min}}\]

This allows for a symmetric scale on both sides, generating negative values for areas overrepresented for KD minima, and positive values for areas overrepresented in KD maxima.


The main findings support the hypothesis that winds from Northeastern Asia are associated and synchronized with periods of high KD Maxima in Tokyo, as the main figure portrays:

Main Figure

Red values represent those areas of more common occurence during weeks of KD maxima with respect to those of KD minima, with blue values representing the opposite.

All in all, there seem to be two main different routes associated to periods of


Temporal association study

There are several limitations when trying to relate studies to connect temporal changes between different processes, as finding associations doesn’t necessarily imply a causal link between both phenomena, and external drivers connecting both might be the reason of the link.

Choice of temporal variability

Taking this into account, the temporal changes in KD incidence in Tokyo has three main sources of variability:

  • The increasing trend.

  • The seasonal variation (mainly, of yearly frequence).

  • The anomalies once removed of trend and seasonal variability.

The methods employed here, given the selection of the yearly maxima and minima, only allow us to associate or identify the phenomena related to yearly variation, and would fail to identify the processes that drive the general increase in cases for the last 20 years and the anomalies unrelated to the seasonal variability.

Backtrajectory modelling

The use of HYSPLIT for modelling the trajectories implies that there is a certain degree of error in the estimation of the trajectories, and the magnitude of these is known to increase as the trajectories surpass the planet boundary layer. By not using single trajectories to assess any of the effects but aggregates of many of them, we attempted to mitigate these issues. The use of higher resolution meteorology data might also allow us to improve the accuracy of the modelling.

In this study, we’ve also just considered the 2D intersections of the trajectories with the latitudinal-longitudinal plane, without regard for the actual altitude of the trajectories. To be more thorough, we should include the altitude of the trajectories in each point in time, paired with calculations of moisture uptake rates in order to able to actually pinpoint with better accuracy the areas where the air masses that arrive to the destination are collecting matter.

Conclusions and outlook

This study further confirms, with updated data, the hypothesis connecting the tropospheric winds from NE Asia with higher Kawasaki Disesase incidence periods which was discussed in [Rodó et al., 2014], allowing us now to pinpoint with higher accuracy the specific locations that might source the (or one of) the etiologic agents of the disease.

Further studies should (and are being done) on the characterization of these air masses to understand the potential immune-triggering compounds carried in them.

Stay posted for more updated on the issue and feel free to comment below with any doubts or comments.


Click to expand!


Ryusuke Ae, Nobuko Makino, Koki Kosami, Masanari Kuwabara, Yuri Matsubara, and Yosikazu Nakamura. Epidemiology, treatments, and cardiac complications in patients with kawasaki disease: the nationwide survey in japan, 2017-2018. The Journal of Pediatrics, 225:23–29.e2, 2020. URL:, doi:


Xavier Rodó, Joan Ballester, Dan Cayan, Marian E. Melish, Yoshikazu Nakamura, Ritei Uehara, and Jane C. Burns. Association of Kawasaki disease with tropospheric wind patterns. Sci Rep, 1(1):152, November 2011. Number: 1 Publisher: Nature Publishing Group. URL:, doi:10.1038/srep00152.


Xavier Rodó, Roger Curcoll, Marguerite Robinson, Joan Ballester, Jane C. Burns, Daniel R. Cayan, W. Ian Lipkin, Brent L. Williams, Mara Couto-Rodriguez, Yosikazu Nakamura, Ritei Uehara, Hiroshi Tanimoto, and Josep-Anton Morguí. Tropospheric winds from northeastern China carry the etiologic agent of Kawasaki disease from its source to Japan. Proceedings of the National Academy of Sciences, 111(22):7952–7957, June 2014. Publisher: Proceedings of the National Academy of Sciences. URL:, doi:10.1073/pnas.1400380111.


Jane C. Burns, Lauren Herzog, Olivia Fabri, Adriana H. Tremoulet, Xavier Rodó, Ritei Uehara, David Burgner, Emelia Bainto, David Pierce, Mary Tyree, Daniel Cayan, and for the Kawasaki Disease Global Climate Consortium. Seasonality of Kawasaki Disease: A Global Perspective. PLOS ONE, 8(9):e74529, September 2013. Publisher: Public Library of Science. URL:, doi:10.1371/journal.pone.0074529.


Nobuko Makino, Yosikazu Nakamura, Mayumi Yashiro, Ryusuke Ae, Satoshi Tsuboi, Yasuko Aoyama, Takao Kojo, Ritei Uehara, Kazuhiko Kotani, and Hiroshi Yanagawa. Descriptive epidemiology of kawasaki disease in japan, 2011–2012: from the results of the 22nd nationwide survey. Journal of epidemiology, pages JE20140089, 2015.


Nobuko Makino, Yosikazu Nakamura, Mayumi Yashiro, Takashi Sano, Ryusuke Ae, Koki Kosami, Takao Kojo, Yasuko Aoyama, Kazuhiko Kotani, and Hiroshi Yanagawa. Epidemiological observations of kawasaki disease in japan, 2013–2014. Pediatrics International, 60(6):581–587, 2018.


Nobuko Makino, Yosikazu Nakamura, Mayumi Yashiro, Koki Kosami, Yuri Matsubara, Ryusuke Ae, Yasuko Aoyama, and Hiroshi Yanagawa. Nationwide epidemiologic survey of kawasaki disease in japan, 2015–2016. Pediatrics International, 61(4):397–403, 2019.


A. F. Stein, R. R. Draxler, G. D. Rolph, B. J. B. Stunder, M. D. Cohen, and F. Ngan. NOAA’s HYSPLIT Atmospheric Transport and Dispersion Modeling System. Bulletin of the American Meteorological Society, 96(12):2059–2077, December 2015. Publisher: American Meteorological Society Section: Bulletin of the American Meteorological Society. URL: (visited on 2022-03-23), doi:10.1175/BAMS-D-14-00110.1.


Mellissa S.C. Warner. Introduction to PySPLIT: A Python Toolkit for NOAA ARL’s HYSPLIT Model. Computing in Science Engineering, 20(5):47–62, September 2018. Conference Name: Computing in Science Engineering. doi:10.1109/MCSE.2017.3301549.