Baseline scenarios of heat-related ambulance transportations under climate change in Tokyo, Japan

Epidemiological data

To construct a forecast model, we used heat-related ambulance transportation data for the period from 2015 to 2019. The latest two years, 2020 and 2021, were omitted to disregard the impact of coronavirus disease (COVID-19) on associated behaviors (e.g., a stay-home policy was in place during the pandemic). Specifically, we used a dataset on the daily number of people transported by ambulance for heatstroke, which records every single transportation event and is openly published by the Fire & Disaster Management Agency (FDMA) (2022a). The ambulance transportation case data were structured by age group; in the following analysis, we analyze the data by dividing the population into three age groups: 0–17 years, 18–64 years, and 65 years or older. In addition to case data, data on the daily maximum WBGT were obtained from the Ministry of the Environment (2022). In addition, the population size for Tokyo for the period from 2015 to 2019 in the abovementioned three age groups was obtained (Tokyo Metropolitan Government, 2022) to calculate the risk of heat-related ambulance transportations. For instance, in 2019, the population sizes of 0–17 years, 18–64 years, and 65 years or older were 1.86 million, 8.25 million and 3.08 million, respectively, and of these, 553, 2,322 and 3,171 heatstroke cases called for ambulance transportation.

The Intergovernmental Panel on Climate Change (IPCC) released its AR6 on August 9, 2021. They found that, by at least the late 2030s, global average temperatures are likely to have increased by 1.5 °C above pre-industrial levels (IPCC, 2021). In the present study, we obtained climate change scenarios based on CMIP6, published by the National Institute for Environmental Studies (Ishizaki, 2021). Of the five different scenario models available, we selected two: (i) Model for Interdisciplinary Research on Climate version 6 (MIROC6), which was cooperatively developed by a Japanese modeling community to precisely reflect the meteorological conditions in Japan, and (ii) Meteorological Research Institute Earth System Model version 2.0 (MRI-ESM-2.0), which was developed by the Japan Meteorological Agency and offers accurate forecasts of rainfall and humidity. MRI-ESM-2.0 was adopted because humidity plays a critical role in calculating the WBGT. For each model, three different time-series scenarios—Representative Concentration Pathway (RCP) 2.6, RCP4.5, and RCP8.5—with daily climatological data from 1900 to 2100 were obtained for each 1-km square of geographic mesh. Using the two models and three time-series scenarios, the daily maximum temperature, relative humidity, total solar radiation, and average wind speed in a geographic space that contains the Imperial Palace of Japan, located in central Tokyo, were extracted to calculate the WBGT for each date (Ono & Tonouchi, 2014). The reason why we chose the area around the Imperial Palace as the location of our analysis is that data on the daily maximum WBGT were consistently available for this particular geographic space over the whole time period.

Mathematical model

For the purpose of prediction, we exploited the relationship between the number of heat-related ambulance transportations per million population and the daily maximum WBGT. These heatstroke data were used to predict the per capita risk because we subsequently used the estimate multiplied by the predicted age-specific population size during the forecasting process. The model was fitted by age group (i.e., 0–17 years, 18–64 years, and 65 years or older). Several different statistical models were employed to describe the dose–response relationship. The first is a classical model that captures a monotonic increase in the risk of heatstroke as a linear function of the WBGT after the WBGT exceeds a threshold value; this is referred to as the hockey stick regression method (Yanagimoto & Yamamoto, 1979). This classical model was specifically revisited because there is an existing notion of a threshold around a WBGT of 28 °C. Using the daily maximum WBGT (T), the number of heat-related ambulance transportations per million population n (T) was modeled as (1)${\displaystyle E\left(n\left(T\right)\right)=\left\{\begin{array}{cc}{\beta }_1,\hfill & \mathrm{for}T<T_w\hfill \\ {\beta }_1+{\beta }_2\left(T-T_w\right),\hfill & \mathrm{for}{T}_w\le T\hfill \end{array}\right.}$

where T_w is the WBGT threshold above which the risk of heatstroke abruptly increases. β₁ is the constant risk at a WBGT below T_w; β₂ is the gradient of risk when the WGBT exceeds T_w. The model was independently fitted to each age group, and T_w was treated as a parameter because different age groups were expected to have different threshold levels (i.e., different physical endurance levels).

The second model is an extension of the first. Frequently, the WBGT is expressed as one of five discrete levels, with the highest level being a WBGT of 31 °C or more; this implies another threshold that boosts the risk of heatstroke. With two threshold temperatures, T_w1 and T_w2, the expected number of transportations per million population n (T) was modeled as (2)${\displaystyle E\left(n\left(T\right)\right)=\left\{\begin{array}{cc}{\beta }_1,\hfill & \mathrm{for}T<T_{w1}\hfill \\ {\beta }_1+{\beta }_2\left(T-T_{w1}\right),\hfill & \mathrm{for}{T}_{w1}\le T<T_{w2}\hfill \\ {\beta }_1+{\beta }_2\left(T_{w2}-T_{w1}\right)+{\beta }_3\left(T-T_{w2}\right),\hfill & \mathrm{for}{T}_{w2}\le T\hfill \end{array}\right.}$

where β₁ is the constant risk at a WBGT below T_w1 and β₂ and β₃ are the gradients of risk when the WBGT exceeds 28 °C and 31 °C, respectively.

The third model is a phenomenological model that assumes an exponential increase in the risk of heatstroke. We allowed a threshold for this model as well. Let r be the exponential rate of increase in risk and let β₁ be the constant risk below the WBGT threshold. We have (3)${\displaystyle E\left(n\left(T\right)\right)=\left\{\begin{array}{cc}{\beta }_1,\hfill & \mathrm{for}T<T_w\hfill \\ {\beta }_{1}exp\left(r\left(T-T_w\right)\right),\hfill & \mathrm{for}{T}_w\le T\hfill \end{array}.\right.}$

Assuming that the variations in the observed counts of heat-related ambulance transportations are sufficiently captured by a Poisson distribution with the expectation given by (1), (2), or (3), maximum likelihood estimation was performed to obtain optimal values of the parameters. The Akaike information criterion (AIC) was calculated to compare the model fit. The mean absolute error (MAE) of the prediction from the observed values was also calculated to evaluate the accuracy of prediction.

Future prediction scenarios of heat-related ambulance transportations

We consistently used the following formula to calculate the WBGT from readily available climatological data (Ono & Tonouchi, 2014): (4)${\displaystyle T_{WBGT}=0.735\times T_{asmax}+0.0374\times RH+0.00292\times T_{asmax}\times RH+7.619\times SR-4.557\times S{R}^2-0.0572\times WS-4.064}$

which is most frequently used in Japan to estimate the WBGT. T_asmax is the maximum daily temperature (^∘C), RH is the relative humidity (%), SR is global solar radiation (kW/m²) measured by horizontally installed solar radiation meter , and WS is the average wind speed (m/s). A validation study of this approximate equation has excellently shown that 98.3% to 99.8% of WBGT estimate involved a bias less than 1.0 °C. In Japan, WBGT is physically measured in 11 different locations (weather observatories), and an approximate WBGT using Eq. (4) has been officially adopted in other 829 observatory stations. The MIROC6 and MRI-ESM-2.0 meteorological data were substituted into (4) to calculate the daily WBGT values from 2000 to 2100. Because the bifurcation of temperature change begins in 2015 in the two climate change models, data from 2000 to 2014 were calculated using common data.

In each future year, the proportion of days with a WBGT value falling into each of five discrete ranges was calculated. This was calculated as the number of days with a WBGT value in each range divided by the total number of days (365). Because climatological variables fluctuate from year to year, the abovementioned proportions were calculated for every five years.

We examined the projected trends in the number of heat-related ambulance transportations, explicitly accounting for future demography. To do so, we first calculated the WBGT from May 1 to September 30 in each year, and used models (1), (2), or (3) to calculate the risk per million population. Even in the long-term future, we assumed that the heatstroke risk would occur only between May and September. The predicted risk was embedded onto the projected age-specific populations aged 0–14 years, 15–64 years, and 65 years or older in the period 2015–2100 (The Climate Change Adaptation Information Platform, 2022). To calculate the 95% confidence intervals of our predictions of heat-related ambulance transportations, the parametric bootstrapping method was employed, because parameters were inferred by maximum likelihood method and the variance–covariance matrix was reasonably computed. To account for parameter uncertainties of models (1), (2), or (3), 5,000 bootstrap resampling experiments were conducted and the 2.5th and 97.5th percentile points were taken.

Data sharing statement

All datasets, including modeled climatological data for the period 2015–2100 and empirical data on heat-related ambulance transportations for the period 2015–2019 in Tokyo are openly shared as the online supporting material (Tables S1–S4).

Ethical considerations

The present study used publicly available information. Because no private information was used, ethical approval was not required.