Seasonal transmission dynamics of varicella in Japan: The role of temperature and school holidays

1.   Introduction

Varicella (chickenpox) is a common childhood disease caused by varicella-zoster virus (VZV), which is a member of the alphaherpesvirus subfamily ^[1]. Healthy children mostly experience relatively mild symptoms, characterized by an itchy vesicular rash. However, VZV infection can be more serious, especially in infants, adolescents, adults and immunocompromised individuals ^[2]. Varicella vaccine is highly effective in preventing varicella: A meta-analysis showed that one-dose and two-dose varicella vaccination reduced varicella incidence by 81 and 92%, respectively ^[3]. In Japan, a two-dose national varicella immunization program was implemented in 2014, targeting children aged 12 to 36 months. The varicella incidence declined by nearly 80% in 2019 compared with that before implementation of the universal vaccination program ^[4]; epidemic patterns after 2019 under the COVID-19 pandemic require careful investigation ^[5]. One concern with initiation of a universal varicella vaccination program is a shift in the average age of primary infection to older age groups ^[6,7], who are at risk of severe varicella; therefore, many countries choose voluntary or risk-based varicella vaccination ^[8]. Prevention of varicella remains a public health challenge, and understanding its transmission mechanisms is crucial for controlling epidemic spread.

Many childhoods infectious diseases, including measles, mumps and varicella, show seasonal epidemic patterns ^[9,10]. In the case of varicella in Japan, major and minor bimodal seasonal patterns of varicella have been observed ^[5]. The underlying mechanisms of seasonality may be influenced by both social patterns and environmental factors. Several transmission modeling studies have indicated that transmission seasonality is driven by the school calendar, with the incidence declining during school holidays and increasing with school reopening ^[11,12,13]. A published modeling study evaluated the effect of the school term on varicella incidence using long-term data from England and Wales during 1967–2008 and concluded that the reduction in contact among children during the summer holiday resulted in a decline in varicella incidence ^[14]. Climate may also affect transmission patterns of varicella. For example, most children contract varicella by age 10 years in temperate climates, whereas children in tropical regions acquire varicella at a later age ^[15]. Another modeling study evaluating the climate effects on varicella indicated that a future decrease in relative humidity driven by climate change would result in an increase in the incidence of varicella ^[16].

Despite efforts to model the seasonal patterns of childhood infectious diseases, the underlying mechanisms that give rise to bimodal and moderately bimodal epidemic patterns over different geographic regions require further investigation. In this study, we aimed to evaluate the effects of the school term and temperature on the incidence of varicella in Japan. Using the reported number of varicella cases from 2000 to 2009, which encompasses the period of the voluntary varicella vaccination policy, we quantified the transmission rate as a function of underlying mechanisms by reconstructing the transmission dynamics.

2.   Materials and methods

2.1.   Epidemiological data

Japan is an island country, located between approximately 20 and 45 degrees north latitude and between approximately 123 and 154 east longitude ^[17,18]. The northern part of Japan (e.g., Hokkaido) has a subarctic climate with heavy snow in winter, whereas the southern part (e.g., Okinawa) has a subtropical climate with hot and humid summers. Japan mainly comprises the four large islands of Hokkaido, Honshu, Shikoku and Kyushu and is divided into 49 prefectures. To evaluate the effect of climate on the varicella dynamics in Japan, we selected seven specific prefectures: (ⅰ) Hokkaido, the northernmost main island; (ⅱ) Miyagi, in the northern part of Honshu; (ⅲ) Tokyo, on the eastern side of Honshu; (ⅳ) Kyoto, on the western side of Honshu; (ⅴ) Fukuoka in northern Kyushu and Kagoshima as the southernmost prefecture of Kyusyu; and (ⅵ) Okinawa, the southernmost main island (Figure A1).

Varicella is designated a category V disease, monitored in sentinel surveillance and notified to the National Epidemiology Surveillance for Infectious Diseases (NESID) system in compliance with the Communicable Disease Prevention Law through March 1999 and, subsequently, the Infectious Diseases Control Law. Approximately 3000 pediatric sentinel sites (representing around 10% of pediatric clinics and hospitals in Japan) report cases to the NESID on a weekly basis ^[19].

Four pieces of information were used in the analysis: (ⅰ) prefecture-dependent weekly number of varicella notifications, (ⅱ) average weekly temperature in the capital city of each prefecture, (ⅲ) prefecture-dependent number of newborns and (ⅳ) the national level vaccination coverage. The average temperature of the observatory station in each capital city was retrieved from the Meteorological Agency of Japan ^[20]. Here, we specifically examined temperature as the climatological input variable because variations in the bimodal patterns of varicella in Japan have been captured well by temperature using time-series analysis ^[21]. The yearly number of live births was used to account for an inflow of susceptibles to the population and was retrieved from each prefectural report, which is based on the national estimate and is available from the Statistics Bureau of Japan ^[22]. We assumed that birth events were evenly distributed throughout the year, and the weekly number of newborns in each prefecture was estimated as the annual number of newborns by prefecture, divided by 52 weeks. Vaccination coverage was estimated as the annual sales of varicella vaccine divided by the number of newborns, per our earlier study ^[4,23], because varicella vaccination coverage was not routinely monitored before the routine immunization program began in 2014.

2.2.   Mathematical model

In the following model, we considered a susceptible–infected–recovered (SIR)-type model with discrete time step as a method to reconstruct the transmission dynamics of varicella using prefecture-dependent notifications and demographic data. The following assumptions were adopted: (ⅰ) All susceptible individuals eventually contract varicella at some point in their life (or become immune by vaccination), (ⅱ) lifelong immunity is attained following natural infection or vaccination (for simplicity and owing to data limitations, we did not differentiate between one- and two-dose vaccination), and (ⅲ) only birth events can contribute to the recruitment of susceptibles. Considering these assumptions, the balance equation of susceptibles by prefecture is written as

where t = 1, 2, 3… denotes the calendar week starting from week 1 of 2000, and p identifies each prefecture. ${S}_{t, p}$, ${B}_{t, p}$ and ${I}_{t, p}$ represent the estimated number of susceptible individuals, the number of newborns and the estimated varicella incidence in calendar week t in prefecture p, respectively. ${v}_{t}$ expresses the vaccination coverage in calendar week t, assuming uniform vaccination coverage across the country. The incidence of infection is assumed as proportional to the number of susceptibles. Let β represent the transmission coefficient, and the standard SIR model describes the incidence as

where ${\gamma }_{t, p}$ is referred to as the force of infection in calendar week t and prefecture p, and N_{t, p} is the population size of the prefecture in week t. The standard SIR model assumes that the force of infection is proportional to the incidence I_{t, p}. in the corresponding time step. As it is difficult to quantify the system using the mass action term, we intend to quantify ${\gamma }_{t, p}$ for the rest of this paper. That does not assume that the model is static, but rather ${\gamma }_{t, p}$ is considered to capture the time-dependent transmission dynamics reflecting time-varying β as well as time-dependent value of I_{t, p}.

In addition to the transmission dynamics, we account for the observation process. Prefecture dependent weekly reported number of cases can be expressed as follows:

where ${\delta }_{p}$ is the prefecture-dependent parameter that represents biases associated with ascertainment as well as the reporting rate.

In our previous study, evaluating the effect of public health and social measures targeting the coronavirus disease 2019 (COVID-19) pandemic on varicella dynamics in Japan ^[5], we implemented trigonometric functions to capture the bimodal seasonal pattern of varicella in Japan. The force of infection at week t in entire Japan, ${\gamma }_{t}$, was expressed as follows:

where n is the order of the trigonometric function, and ${a}_{0}$, ${a}_{l}{\mathrm{a}\mathrm{n}\mathrm{d}b}_{l}$ are parameters that govern the function. As an alternative way to capture the bimodal seasonal pattern by prefecture, we modeled the force of infection as the interaction between two factors: temperature variation and the school calendar.

Instead of using the actual temperature value, we chose to use the deviation from the threshold temperature, assuming the existence of the preferable temperature for varicella transmission. We presumed that temperature directly affected the infection event, and the force of infection changes with temperature variation. Moreover, as the present study intends to fit the above-mentioned model to the observed incidence data using Eq (3) over time, we use the power-law approximation in the transmission term rather than adhering to the mass action principle in Eq (2). Employing the power law approximation with the power exponent α to the incidence, the prefecture-dependent force of infection in calendar week t., ${\gamma }_{t, p}$, is defined as

where ${\beta }_{i, p}$ is the same as Eq (2), i.e., the prefecture-dependent transmission coefficient with respect to the change in contact patterns during the school term (i = 1) and the summer holiday (i = 2) that starts in late July and ends at the end of August. We assumed these transmission rates took constant values (i.e., a seasonal two-step function) from 2000 to 2009. The effect of temperature on the transmission rate is expressed as the exponential function of the temperature deviation from the threshold temperature, $\stackrel{-}{t}.$ In other words, $\stackrel{-}{t}$ is the most preferable temperature for varicella transmission. $\mu$ is a scaling factor of the temperature deviation that is common to all prefectures. ${temp}_{t, p}$ is the temperature in calendar week t in prefecture p. We introduced the parameter $\alpha$ to capture nonlinearities of the model that reflected the heterogeneity in contact patterns. This method has been conventionally used for the time series SIR models ^[24].

In this way, the system is parameterized, and parameters that govern the transmission kinetics in Eq (5) are estimated. Using Eqs (3) and (5), the weekly reported number of cases can be described as follows.

which is log-linearized as

It should be noted that Eq (7) can now be used to be fitted to the observed weekly data to estimate unknown parameters $\left({\beta }_{i, p}, \mu , \stackrel{-}{t}, \alpha \right)$. That is, we fitted Eq (7) to the weekly number of varicella notifications using a generalized linear model and assuming that the variation is sufficiently captured by Gaussian distribution (verifying the assumption by checking residuals).

Before estimating parameter values $\left({\beta }_{i, p}, \mu , \stackrel{-}{t}, \alpha \right)$, we calculated the parameter value ${\delta }_{p}$, using linear regression analysis between the cumulative yearly number of newborns and cumulative number of varicella notifications during 2000–2009. This calculation involves an assumption that all children are born as susceptible and eventually become immune once vaccinated or infection, and the slope of the regression model can be used as the parameter value ${\delta }_{p}$(4, 5). Subsequently, ${\delta }_{p}$ was dealt with as a known parameter when using Eq (7) to fit to the observed notification data. An epidemiological study in Japan indicated that VZV was activated at a temperature ranging from 5 to 20 ℃. We were unable to jointly estimate $\stackrel{-}{t}$ explicitly with other parameters; therefore, we performed regression with a fixed value of $\stackrel{-}{t}$ in the range 5 to 20 ℃ and specifically chose 11.0 ℃, which commonly yielded the smallest Akaike information criterion.${S}_{t, p}$ was fixed using the estimated proportion of susceptibles estimated in our earlier study ^[5] (i.e., ${S}_{p, t}$ was calculated as the product of the proportion of susceptibles in calendar week t and ${N}_{p, t}$). Using the estimated parameters, we reconstructed the epidemiological dynamics and calculated the transmission rate and force of infection during 2000–2009. All analyses were performed with R software version 4.0.2 (The R Project for Statistical Computing, Vienna, Austria)

2.2.1.   Ethics approval of research

This study used publicly available data that contained no private information. Therefore, ethical approval was not required.

3.   Results

Figure 1 shows the epidemic curve for the weekly number of varicella notifications, the average weekly temperature and average weekly temperature deviation from the threshold (11.0 ℃) by prefecture from 2000 to 2009. In Hokkaido, where the lowest average weekly temperature was below 0 ℃ in the winter season, a clear bimodal pattern in the epidemic curve was visually identifiable. A similar bimodal pattern was seen in the average weekly temperature deviation from the threshold value. These bimodal patterns flattened toward the southern part of Japan. Eventually, a unimodal epidemic pattern was seen in Okinawa, in southernmost Japan, with low annual climate variability. Table 1 shows the estimated values for the transmission parameters during the school term $\left({\beta }_{1, p}\right)$ and summer holiday $\left({\beta }_{2, p}\right)$, and their ratio. An increasing trend was observed in the transmission rate ratio of the summer holiday and other estimated parameter values $\left({\delta }_{p}, \mu , \alpha \right)$, shown in Table A1.

Figure 1.  Prefecture-dependent weekly number of varicella notifications (a)–(g), weekly average temperature (h)–(n) and weekly average temperature deviation from the 11.0 ℃ threshold (o)–(u), from 2000 to 2009.

DownLoad: Full-Size Img PowerPoint

Table 1.  Estimated values for the transmission parameters during the school term and summer holiday, and their ratio for the seven prefectures investigated. Numbers in parentheses represent 95% confidence intervals derived from profile likelihood.

	${\beta }_{1, p}$ (School term)	${\beta }_{2, p}$ (Summer holiday)	${\beta }_{2, p}/{\beta }_{1, p}$ (Ratio)
Hokkaido	60.18 (58.10, 62.25)	46.82 (44.72, 48.91)	0.78
Miyagi	45.23 (43.17, 47.28)	33.09(31.01, 35.16)	0.73
Tokyo	73.02 (70.93, 75.11)	64.92 (62.82, 67.03)	0.89
Kyoto	51.83 (49.76, 53.89)	47.72 (45.63, 49.81)	0.92
Fukuoka	59.67 (57.59, 61.74)	56.19 (54.10, 58.29)	0.94
Kagoshima	36.09 (34.03, 38.15)	35.52(33.44, 37.60)	0.98
Okinawa	52.05 (49.98, 54.11)	44.80 (42.71, 46.87)	0.86

---

 | Show Table

DownLoad: CSV

Figure 2 shows the time series of the estimated transmission rate, including the effect of annual climate change and the force of infection, by prefecture from 2000 to 2009. The transmission rate and force of infection were lowest during the summer season in all seven prefectures, and both showed a bimodal pattern with a small trough during early spring in Hokkaido and Miyagi. The bimodal pattern became flat in Tokyo, Kyoto, Fukuoka and unimodal in Kagoshima and Okinawa. This trend was also observed in the average weekly temperature deviation from the threshold.

Figure 2.  Prefecture-dependent estimated transmission rate (a)–(f) and estimated force of infection (g)–(l) from 2000 to 2009. Black shaded areas represent 95% confidence intervals.

DownLoad: Full-Size Img PowerPoint

A comparison of the observed and estimated number of varicella notifications from 2000 to 2009 is shown in Figure 3. Our model considered both the school term and temperature effect overall, and it captured the observed bimodal patterns in northern prefectures, which gradually diminished and eventually become unimodal in the south. The estimated number of susceptibles by prefecture is shown in Figure A2.

Figure 3.  Comparison of observed and estimated number of notifications from 2000 to 2009. Black dots indicate the prefecture-dependent weekly number of notifications and blue lines show the estimated number of notifications. Blue shaded area represents 95% confidence interval for the fitted line.

DownLoad: Full-Size Img PowerPoint

Figure 4 depicts the decomposed mechanism of the bimodal and partly unimodal seasonal epidemic patterns of varicella in Japan. Because the extent of seasonal temperature deviation from the threshold is high in Hokkaido, a bimodal pattern is generated. However, the stable temperature in Okinawa results in a unimodal pattern. A duration of summer school holidays of more than 5 weeks acts as a suppression factor of transmission, temporarily reducing the varicella incidence in summer.

Figure 4.  Decomposed mechanism of bimodal and unimodal seasonal epidemic patterns in Japan using data in 2000. (a) Hokkaido, (b) Okinawa. Red lines and black dotted lines show the weekly number of notifications and weekly average temperature deviation from the threshold (11.0 ℃). Black horizontal steps show the estimated transmission rate. Dotted vertical lines indicate the summer holiday.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

In the present study, we investigated the impact of both school term and temperature on the transmission dynamics of varicella in Japan, using demographic, meteorological and surveillance data. First, we estimated the transmission parameters that reflected the summer holiday season by fitting a generalized linear model to the weekly number of varicella notifications. Using the estimated parameter values, we reconstructed the transmission dynamics of varicella from 2000 to 2009 in Japan. We set the two transmission parameters, allocating one each for the summer holiday and the remaining period, and modeled the transmission rate as the product of the transmission parameters and the exponential function of the temperature deviation from the threshold value. Our modeling method was simple, but it captured the bimodal pattern in the northern part of Japan, where the annual temperature variability from the threshold was large. The model also successfully reproduced the unimodal pattern in the south, where annual variability from the threshold is small. To the best of our knowledge, the present study was the first modeling study to mechanistically investigate the impact of temperature and school term on the seasonal dynamics of varicella transmission.

Our study findings highlight two important take-home messages. First, our study indicated that preferable temperatures exist for varicella transmission. In Hokkaido, where the temperature during winter season normally drops below 0 ℃, two main peaks of the temperature deviation from the threshold were observed in summer and winter. These peaks corresponded to the peak and trough (i.e., bottom) of the epidemic curve. However, in Okinawa, which has hot summers and warm winters, the temperature deviation from the threshold remained moderate. Consequently, the deviation pattern was unimodal, with relatively small amplitude, compared with other prefectures, and the epidemic curve also showed a unimodal pattern. A published time-series study evaluating the association between seasonal pattern of varicella and meteorological factors indicated that a second peak of the bimodal pattern in Japan was observed at a temperature of approximately 8.5 ℃ ^[21], which was not far from our threshold value (11 ℃), which accounted for the transmission mechanisms in an explicit manner.

Second, both the summer holiday and temperature affect the decrease in varicella incidence during the summer season. An interactive effect may exist between behavioral changes and climate conditions during the summer season, when the temperature deviation from the threshold is high. This finding has important implications for future studies of climate change. Once the temperature is further elevated, the earlier peak in the epidemic curve may be enhanced. Moreover, the impact of the summer school holidays in suppressing the epidemic during late summer may be further intensified in the future.

Several limitations of this study should be noted. First, we did not consider the age-dependent transmission rate, even though we investigated the effect of the school calendar. Because most varicella notifications were among children under age 15 years, the population-level transmission rate was sufficient to capture the effect of the school term on seasonal trends. It should be remembered that age-dependent heterogeneity was not addressed in the present study. Second, we only focused on temperature in the meteorological data. Climate conditions, including humidity, may also affect the seasonal pattern of infectious diseases. Among directly transmitted diseases, the impact of humidity and temperature on the transmission of influenza has been well investigated, with drier and colder climates presenting more suitable conditions for transmission of the influenza virus ^{[25,26,27,28,29,30,31]}. Regarding varicella, one modeling study conducted in Mexico on the association between relative humidity and the incidence of varicella concluded that a decrease in relative humidity is favorable for varicella ^[16]. In our study, we tried to include the effect of humidity in the model, but it did not improve the model fitting, and we chose to include the effect of temperature only in the model. Interestingly, the model captured the bimodal pattern in the northern part of Japan (i.e., Hokkaido and Miyagi), while it did not capture the peak of the unimodal pattern in the southern part of Japan (i.e., Kagoshima and Okinawa). The kinetics of temperature may vary depending on climate conditions, or unknown mechanism might exist. Additional evidence, including the results of experimental studies, is required regarding the effect of climate conditions on the transmission of varicella. Third, we only investigated the effect of behavioral change during the summer holiday. In Japan, there are also school holidays in winter and spring, although the durations are shorter. The summer holiday lasts for approximately 6 weeks, whereas the winter and spring holidays are only 2 weeks long each, which is too short to evaluate the effect of behavioral changes. Lastly, we did not discuss the effect of climate on the seasonality of herpes zoster, which is caused by reactivation of latent VZV infection, especially among the elderly. The incidence of herpes zoster in Japan shows an opposite seasonal pattern to varicella: It increases in summer and decreases in winter ^[32]. Evaluating the seasonality of herpes zoster using a similar method to what we used in the present study would be the subject for future consideration.

5.   Conclusions

In the present study, we quantitatively assessed the impacts of school term and annual change in temperature on the transmission dynamics of varicella in Japan. Our findings suggest the existence of preferable temperatures for VZV transmission and an interactive effect of school contact patterns and temperature. Further studies are required to investigate the underlying mechanism of transmission as well as the effect of future climate change. Assessing how meteorological conditions, including humidity and the amount of rainfall, affect varicella transmission is an important topic for future investigation.

Acknowledgments

HN received funding from a Health and Labor Sciences Research Grant (20CA2024, 20HA2007, 21HB1002, and 21HA2016), the Japan Agency for Medical Research and Development (AMED; JP20fk0108140, JP20fk0108535, and JP21fk0108612), the Japan Society for the Promotion of Science (JSPS) KAKENHI (21H03198 and 22K19670), Environment Research and Technology Development Fund (JPMEERF20S11804) of the Environmental Restoration and Conservation Agency of Japan, Kao Health Science Research, Daikin GAP fund program of Kyoto University and the Japan Science and Technology Agency (JST) SICORP program (JPMJSC20U3 and JPMJSC2105) and RISTEX program for Science of Science, Technology and Innovation Policy (JPMJRS22B4). We thank local governments, public health centers and institutes for surveillance, laboratory testing, epidemiological investigations and data collection. The funders had no role in the study design, data collection and analysis, decision to publish or preparation of the manuscript. We thank Analisa Avila, MPH, ELS, of Edanz (https://jp.edanz.com/ac) for editing a draft of this manuscript.

Conflict of interest

The authors declare no conflicts of interest.

Appendix

Table A1.  Estimated parameter values for the model.

$\delta$ (Hokkaido)	$\delta$ (Miyagi)	$\delta$ (Tokyo)	$\delta$ (Kyoto)	$\delta$ (Fukuoka)
0.425 (0.423, 0.427)	0.415 (0.413, 0.417)	0.110 (0.109, 0.112)	0.286 (0.285, 0.287)	0.400 (0.397, 0.404)
$\delta$(Kagoshima)	$\delta$ (Okinawa)	$\mu$	$\alpha$
0.526 (0.525, 0.530)	0.334 (0.333, 0.335)	-0.021 (-0.018, -0.023)	0.791 (0.775, 0.807)

---

 | Show Table

DownLoad: CSV

Figure A1.  Location of the seven study prefectures: Hokkaido, Miyagi, Tokyo, Kyoto, Fukuoka, Kagoshima and Okinawa.

DownLoad: Full-Size Img PowerPoint

Figure A2.  Estimated number of susceptibles by prefecture from 2000 to 2009.

DownLoad: Full-Size Img PowerPoint