The Epidemic Risk of Dengue Fever in Japan: Climate Change and Seasonality

Dengue fever is a leading cause of illness and death in the tropics and subtropics, and the disease has become a threat to many nonendemic countries where the competent vectors such as Aedes albopictus and Aedes aegypti are abundant. The dengue epidemic in Tokyo, 2014, poses the critical importance to accurately model and predict the outbreak risk of dengue fever in nonendemic regions. Using climatological datasets and traveler volumes in Japan, where dengue was not seen for 70 years by 2014, we investigated the outbreak risk of dengue in 47 prefectures, employing the temperature-dependent basic reproduction number and a branching process model. Our results show that the effective reproduction number varies largely by season and by prefecture, and, moreover, the probability of outbreak if an untraced case is imported varies greatly with the calendar time of importation and location of destination. Combining the seasonally varying outbreak risk with time-dependent traveler volume data, the unconditional outbreak risk was calculated, illustrating different outbreak risks between southern coastal areas and northern tourist cities. As the main finding, the large travel volume with nonnegligible risk of outbreak explains the reason why a summer outbreak in Tokyo, 2014, was observed. Prefectures at high risk of future outbreak would be Tokyo again, Kanagawa or Osaka, and highly populated prefectures with large number of travelers.


Introduction
Dengue fever is a mosquito-borne infectious disease, seen in most tropical and subtropical areas in the world. It is reported by the World Health Organization (WHO) that dengue cases have increased dramatically around the world, identifying Americas, South-East Asia, and Western Pacific regions as the most seriously affected regions [1]. In 2019, the Philippines reported 420,000 cases and Malaysia reported 131,000 suspected cases of dengue.
In Japan, dengue epidemic was continuously observed from 1942 to 1945, and, afterwards, the country remained dengue-free for about 70 years. In these days, Japan is not a dengue-endemic country, and approximately 50 to 200 imported and confirmed cases are annually reported, which is considered to be increasing steadily [2]. In 2013, the possibility of autochthonous dengue in Japan was suspected, because a German tourist returning from Japan and without any recent history to dengue-endemic countries was diagnosed as dengue fever when the patient was back in Germany. Dengue has become a threat to many nonendemic countries where the competent vectors such as Aedes albopictus and Aedes aegypti are abundant. WHO has reported that local transmission was seen for the first time in France and Croatia in 2010 [1]. Subsequently, the first local outbreak in Tokyo was also reported in 2014, involving 160 confirmed dengue cases from August to October [3]. In this way, imported cases can contribute to inducing local transmissions of dengue in Japan, indicating that any temperate zone countries that welcome travelers from endemic countries face the threat of local dengue outbreak.
Owing to an increase of the traveler volume and the importation of dengue in Japan, the investigation of the outbreak risk is increasingly recognized as more important than before. Nakamura et al. studied the risk of dengue among Japanese travelers and demonstrated that the risk is greater during epidemic season in the tropics [4]. Fukusumi et al. studied the monthly and yearly notification trends of dengue among Japanese travelers and found that the trend among Japanese travelers closely reflected the transmission dynamics in the travelers' destinations, varying seasonally and annually [2]. Yuan and Nishiura estimated that the size of infected travelers may be more than 20 times the notified number of confirmed imported cases [5]. ese published studies underscore the critical importance of travelers from dengue-endemic countries in characterizing the outbreak risk in Japan.
Climatological factors have been reported to be influential to many mosquito-borne diseases including dengue fever. It can not only affect the life cycle of mosquitoes but also affect the transmission probability and incubation period. e effective reproduction number varying with time, and thus with climatological variables, has been studied by various researchers, accounting for the relationship between climatological factors and epidemiological parameters, such as the biting rate and transmission probability per bite [5][6][7][8][9][10]. Assuming that the basic reproduction number is the mean of the offspring distribution of a Galton-Watson branching process, the probability of extinction or outbreak of a disease can also be assessed [11,12]. However, such assessment did not take place in the published studies of dengue fever, while seasonal variation of the reproduction number was dealt with. e purpose of the present study is to analyze the risk of local outbreak of dengue by prefecture in Japan, considering epidemiological impact of climatological factors and the number of imported cases on the risk of outbreak. We thereby provide fundamental insights into the time-dependent outbreak risk to which the future control strategy can be aligned.

Data Source.
Fukusumi et al. [2] indicated that 70-90% of imported cases in Japan are from India, Indonesia, ailand, and Philippines. Adding Malaysia to this list, the analysis of the cumulative risk from 2006 to 16 by Yuan and Nishiura [5] shows that these five countries are the five leading countries of origin of dengue for Japan. us, here we investigate the dataset of travelers from these five countries. Monthly datasets of travelers arriving from India, Indonesia, Philippines, ailand, and Malaysia from 2013 to 2016 were obtained from the Japan National Tourist Organization (http://www.jnto.go.jp/jpn/statistics/visitor_ trends/indexhtml). Annual dataset of travelers' destination in Japan, 2013, was also extracted from the Japan National Tourist Organization. e distribution of prefectures that were visited by travelers was obtained. Combining two pieces of data, the distribution of prefectures by country of origin and month was obtained. e annual total numbers of imported cases from 2013 to 2016 were derived from nationwide surveillance data (National Epidemiological Surveillance of Infectious Diseases (NESID)), as shown in Table 1. e dataset enabled us to capture yearly variations of imported cases which may be caused by different epidemic dynamics in countries of origin, associated with dominant serotypes and genotypes. Subsequently, the monthly number of imported cases by prefecture was calculated by multiplying the total number of imported cases to the distribution of travelers by month and prefecture. e monthly temperature by prefecture from 2013 to 2016 was also retrieved from Japan Meteorological Agency (http://www.jma.go.jp/jma/index.html).

e Effective Reproduction Number.
Let μ be the mortality rate of mosquito vector, r be the recovery rate of human host, m be the vector-to-host ratio, a be the biting rate of vector, b be the transmission coefficient from human to vector, and let c be the transmission coefficient from vector to human. en, the effective reproduction number can be written as follows [13][14][15]: where T is the monthly mean temperature and EIP is the extrinsic incubation period of the virus. Aedes albopictus is the major vector species of dengue transmission in Japan. Of the six temperature-dependent parameters in the expression of the effective reproduction number, to our knowledge, only the mortality rate and the biting rate are available for Aedes albopictus as the published evidence [16,17], and they are quantified as where μ A � 0.02 and 1− μ A stands for the maximum adult mortality rate. Due to the shortage of studies on Aedes albopictus, known temperature-dependent relationship over the transmission probability [6,8] and the extrinsic incubation period were derived from those for Aedes aegypti. Besides, a reduction factor of 0.7 was multiplied to the probability of transmission per bite to human (c) for Aedes albopictus relative to Ae. Aegypti, based on experimental evidence [6,15,18]. e relationships between temperature and the transmission probability and extrinsic incubation period of Aedes aegypti are described as Expected length of the extrinsic incubation period duration is [15] e infectious period of dengue is estimated to range from 4 to 12 days [1], and it is usually assumed to be 5 days [19]. e vector host ratio m was set to be 0.37, by fitting the default value of the reproduction number at 3 in August in Tokyo.

A Model for Dengue Extinction.
We consider a statistical model to assess the probability of extinction of dengue in Japan. As shown in published studies [11,12,20], the branching process model has been employed to estimate the probability of extinction of various diseases by assuming that the basic reproduction number to be the mean of the offspring distribution. To account for the effect of climatological factors on the transmission of dengue, we assume that the reproduction number varies with temperature. We thus employ the time-varying branching process to approximate the probability of extinction of dengue over the course of calendar time.
Denote Z n is the number of cases of generation n; Y n,i is the number of secondary cases generated by the i-th case in generation n.
Suppose that the offspring distribution of generation n is a distribution with mean R n where R n is the effective reproduction number of generationn. Denote the probability generating function (p.g.f.) for Y n,j by G n (s) � ∞ l�0 p n,l s l , with probability P(Y n,i � l) � p n,l , and let G Z n (s) � ∞ k�0 q n,k s k be the p.g.f. for Z n . e following composition is obtained: One finds that this is a p.g.m. for Y n,1 + Y n,2 + · · · + Y n,Z n and just Z n+1 . So, we have G Z n (s) � G Z n− 1 (G n− 1 (s)). By descending the recursion, with taking Z 0 � 1 being encoded by G Z n (s) � G 0 (G 1 (. . . G n− 1 (s)) . . .), the mean of Z n is Let p e n be the probability that the process is extinct by generation n; that is, because by the definition of the p.g.f., we have If Z 0 � n 0 , then the extinction probability by generation n will be (p e n ) n 0 . So, the probability that the disease is not extinct by generation n can be shown as 1 − (p e n ) n 0 , which can be also denoted by the outbreak probability. If there are I imported cases, each imported case can be treated as the initial of a branching process. Denote p e ni , i � 1, 2, . . . , I to be the extinction probability of the ith imported case in generation n. en, the total extinction probability will be p e n1 p e n2 . . . p e nI .
us, the outbreak probability will be 1 − p e n1 p e n2 . . . p e nI .
Here, the offspring distribution of generation n is denoted to be a negative binomial distribution with mean R n and dispersion parameter k, which includes the Poisson distribution (1/k ⟶ ∞) and geometric distribution (k � 1) as special cases. en, the probability generating function (1), the probability of extinction if one dengue infectious individual is introduced is subsequently obtained. Figure 1 shows the distribution of travelers from India, Indonesia, Philippines, ailand, and Malaysia from 2013 to 2016. As shown in the figure, distributions of travelers in these four years share common qualitative patterns. Travelers were dominated by visitors from ailand and Malaysia. ere are two peaks of traveler volumes in one year. One is from March to May and the other is from October to December. Comparing four subfigures, the number of travelers from 2013 to 2016 has increased with time.

Traveler Volumes.
To investigate the geographic distribution of the travelers by season, the map of traveler volume in 2015 was illustrated ( Figure 2). According to seasonal climate conditions in Japan, one year is divided into four seasons: spring (from March to May), summer (from June to August), autumn (from September to November), and winter (from December to February). Figure 2 reveals that the traveler volume in spring is the largest and the traveler volume in summer is the lowest. Moreover, travelers' favorite place is the North and Midland of Japan, and Hokkaido, Tokyo, and Osaka are the most common destinations. e average monthly volumes of these three prefectures are higher than 5000 persons regardless of the season. Besides, Chiba, Kanagawa, Shizuoka, Kyoto, and Hyogo are also identified as prefectures with high number of travelers. e monthly number of travelers of the most prefectures in the southern Japan are smaller than 2000 persons in summer, autumn, and winter.

e Effective Reproduction in Tokyo.
Substituting the monthly mean temperature in 2013-2016 into the expression of parameters in (2)-(5), we can get the monthly variation of the biting rate, probability of transmission, extrinsic incubation period and the mortality rate. Of these, monthly biting rate and mortality rate are theoretical values and computationally obtained from (2). e monthly variation of the effective reproduction number in Tokyo is shown in Figure 3(a). It follows from this figure that the effective reproduction number is almost zero in the winter season and early spring (from November to April) and reaches its peak in July and August. e generation interval consists of four components, that is, (i) the mean extrinsic (mosquito) incubation period (5-15 days), (ii) the mean intrinsic incubation period (4-7 days), (iii) mean host infectious period (3-7 days), and (iv) mean adult mosquito lifespan (6-15 days) [21]. We assume the generation interval to be a constant 30 days, and then there would be exactly one generation per month. Let the effective reproduction number to be the same within one month. Using the expression of the extinction probability for the branching process as shown in (1), we can calculate the extinction probability if one infectious individual is imported into Tokyo. Figure 3(b) shows the probability of extinction from July to December with different values of the dispersion parameter k, when an infectious individual is imported in June. is figure indicates that the probability of extinction from July to November is similar, and it finally reaches 1 in December due to the low temperature; such dependence is preserved when the parameters under multiple generations and varying temperature are being taken into consideration. Besides, large dispersion parameter leads to a big extinction probability. So, in the following, we mainly study the extinction probability in November considering different importation time.
To understand the month in which the outbreak risk is the highest, the probability of extinction in November versus different importation time, when one infectious individual is imported into Tokyo is shown in Figure 3(c). Effects of different values of dispersion parameter on the extinction probability are also shown in the figure. It indicates that the probability of extinction in November is 1 if an infectious individual is imported when the temperature is very low. If the infectious individual is imported in summer, the extinction probability would be far smaller. e lowest extinction probability reaches when the infectious individual is imported in July. Besides, the value of the dispersion parameter can affect the extinction probability greatly. When k � 0.5, the extinction probability by November if the infectious individual is imported in July is only 0.44. However, when k � 10, the extinction probability by November if the infectious individual is imported in July is 0.87. e figure also shows that the difference of the extinction probability among these different values of k is smaller when the infectious individual is imported in the month with relative low temperature.

e Geographical Distribution of the Effective Reproduction Number.
To compare the effective reproduction number by prefecture and season, we drew the heat map of the mean reproduction number in 2013-2016 in different seasons for the 47 prefectures as shown in Figure 4. It follows from this figure that the reproduction number in summer is much higher than that in other seasons. In summer, the effective reproduction numbers in all prefectures except for    Hokkaido and Aomori are higher than 1. Besides, in autumn, only in Kagoshima and Okinawa, the effective reproduction number is higher than 1. e reproduction numbers in spring and in winter in all prefectures are lower than 1. Obviously, the reproduction number in winter is the lowest, which is below 0.1.

3.4.
e Geographical Distribution of the Outbreak Risk. e heat map of the mean outbreak probability for the 47 prefectures using data from 2013 to 2016 is shown in Figures 5-7, in which the dispersion parameter k � 0.5, 1, and 5, respectively. According to the results of Figure 3, the probability of extinction from July to November is similar, and it finally reaches 1 in December because of the low temperature. e outbreak risk in this study is defined as the probability that the dengue fever is not extinct until November. From Figures 5-7, we can see that importation of infectious individual in July results in the highest risk of outbreak. Besides, comparing these three figures, we can see that the outbreak risk when k � 5 is the lowest. In other words, small dispersion parameter leads to a high probability of outbreak. When k � 5, the outbreak risks in all prefectures are lower than 0.3 even in summer. However, when k � 0.5, the outbreak risks in the south of Japan are much higher. We can see clearly that there is a high risk of outbreak in prefectures in the south of Japan, especially in the southern coastal areas. Miyazaki, Kagoshima, Nagasaki, and Okinawa are always the prefectures with highest risk (above 0.4). However, in northern coastal areas such as Tottori, Shimane, and Fukui, although they are also in the south of Japan, there is a lower probability of outbreak (below 0.3). Besides, it is interesting that Shiga and Nara are also with relatively low probabilities of outbreak although prefectures around them are at high risk. Shizuoka, Gifu, Aichi, and Kanagawa form the boundary of the high risk and low risk areas. All prefectures located in the north of these four prefectures have a very low probability of dengue outbreak.
In Figures 2 and 5, we have shown the traveler volume and the outbreak risk when one individual is imported. To combine the traveler volume into the outbreak risk, the actual outbreak risk is assessed in the following. According to the geographical distribution and the monthly distribution of travelers and incorporating the number of imported cases every year as shown in Table 1, we can calculate the average imported cases in each month and each prefecture. As shown in Figures 5-7, the outbreak risk is higher when k is smaller. We consider the case of k � 0.5 in the following analysis. e map of the actual probability that dengue fever sustains until November can be shown in Figure 8. Comparing four figures, we identify that the outbreak risk in 2015 was the lowest.

Discussion
In this study, we investigated the outbreak risk of dengue fever in Japan, accounting for the traveler volume and temperature of every prefecture. First of all, according to the travel data, the distribution of traveler was obtained over month, prefecture and country of origin. en, using the relationship between temperature and transmission parameters, the effective reproduction number was estimated by month and prefecture. To assess the outbreak risk of dengue fever, the theory of branching process was used to calculate the probability of extinction, thereby yielding the outbreak risk of dengue across all 47 prefectures in Japan. e time distribution of travelers in Japan showed that the numbers of travelers in March, April, May, October, November, and December in 2015 are larger than 150000  persons, while the number of travelers in summer is far smaller. However, the map of the effective reproduction number indicates that the reproduction number reaches its peak in summer, which is bigger than 1 in most prefectures in Japan. Also, it is higher in autumn than in spring, which is similar to the results of the previous paper studying the relative vectorial capacity of dengue in Japan 2004-2013 [3].
In published studies, the Galton-Watson branching process has been used to estimate the basic reproduction number and the probability of extinction of some diseases, such as influenza and pneumonic plague [12,20,22]. e subcritical, critical, and supercritical cases can be explicitly examined. In subcritical and critical cases, the extinction probability is 1, while in supercritical cases, there is a positive probability of extinction. However, in fact, the reproduction number is time-varying due to many factors, such as climatological factors and control measures. In the present study, we used time-varying branching process to model the transmission of dengue, assuming that the offspring distribution was negative binomially distributed with mean R n and dispersion parameter k. As shown in Figure 4, the effective reproduction number varies largely by season and by prefecture. So, the probability of extinction if one infectious individual is imported varies greatly by variations in the time (season) of importation and geographic position.
According to our results, July is at the highest risk of outbreak, and the south and central areas of Japan are high risk areas (as shown in Figure 5). However, incorporating effect of traveler volumes, the actual outbreak risks in many prefectures as shown in Figure 8 are much higher than those shown in Figure 5, especially in Chiba, Tokyo, Kanagawa, Shizuoka, Kyoto, Osaka, and Hyogo. Moreover, there are also some prefectures not at high risk as expected, such as Nagano, Toyama, Fukui, Shiga, Tottori, Shimane, and Tokushima. Although the reproduction numbers in these prefectures are also very large, the traveler volumes are not  so high, leading to relatively small outbreak risks compared with urban locations. Hokkaido is a tourist attraction in Japan with more than 5000 persons per month throughout the year; it has very low outbreak risk due to low temperature. Tokyo was always the prefecture at the highest risk, not only because of high temperature but also due to large travel volumes. e finding is consistent with observing the local outbreak in Tokyo in 2014 [23][24][25].
Furthermore, according to the results of Figure 4.1, the traveler volume has increased over time from 2013 to 2016. Larger traveler volumes would lead to larger outbreak risk. However, as shown in Figure 4.8, the outbreak risk in 2015 is much lower than that in 2013 and 2014, due to the low temperature in 2015. Besides, it is shown in Figure 5 that tourists in June and July need special attention, and according to Figure 1, there are more tourists  is indicates that it is not advisable to judge the probability of outbreaks solely based on the number of imported cases or climatological factors. Multiscale model with multiple factors, especially mobility, would act as the key element of dengue prediction [26].
e variance-to-mean ratio for the negative binomial distribution is 1 + (R 0 /k), so the smaller values of k, the greater heterogeneity of infectious individual. To study effect of heterogeneous transmission, sensitivity analysis for the dispersion parameter k was carried out by investigating the probability of extinction under different values of k. Our results indicate that greater heterogeneity of infectious individual leads to a larger probability of extinction. Also, the differences in the probability of extinction under different values of k is larger, if the reproduction number becomes greater. ese results are consistent with the published study of individual level variations in disease transmission [20,27,28]. Besides, according to the actual situation of Japan, there is no local outbreak of dengue fever before 2014, so the outbreak probability should not be too large, and to give a lower limit of k to interpret risk maps, it is appropriate to assume k � 0.5.
Our study provides a method to assess the outbreak risk of vector borne disease, considering both climatical factors and traveler volumes. It can predict where would be at the highest risk of outbreak and when would be the most hazardous time of importation. As the number of travelers increases, increasing the strength of surveillance targeting travelers would be required, and mosquito control may have to be intensified. Geographical distributions of outbreak risk can provide fundamental insights for public health departments into developing control plans. e proposed method could also be applicable to many other seasonally varying infectious diseases. It can allow us to consider effects of control measures on the probability of extinction. However, there are also some limitations. We assumed that the distribution of imported cases mirrored the distribution of traveler volumes. In fact, dengue fever is seasonal disease, and the number of imported cases may not only depend on the traveler volume but also depend on the prevalence in the country of origin [2,5]. More precise estimation can be attained if the prevalence is taken into account. As an important assumption that is conventionally adopted in other published studies, the ratio of human to mosquito and the generation interval were assumed as constant. More realistic and quantitative models are called for, and the presented framework in the present study would act as the basis for such future improvement.

Data Availability
e present study rests on mathematical modelling. Simulation results data are available upon request. Disclosure e funders played no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflicts of Interest
e authors have declared that no conflicts of interest exist.

Authors' Contributions
H. N. conceptualized the study, developed the methodology, performed project administration, performed supervision, and reviewed and edited the article. H. N. and X. W. performed formal analysis, were responsible for funding acquisition, performed investigation, provided the resources, and performed validation. X. W. wrote the original draft.