Coupling the circadian rhythms of population movement and the immune system in infectious disease modeling

The dynamics of infectious diseases propagating in populations depends both on human interaction patterns, the contagion process and the pathogenesis within hosts. The immune system follows a circadian rhythm and, consequently, the chance of getting infected varies with the time of day an individual is exposed to the pathogen. The movement and interaction of people also follow 24-hour cycles, which couples these two phenomena. We use a stochastic metapopulation model informed by hourly mobility data for two medium-sized Chinese cities. By this setup, we investigate how the epidemic risk depends on the difference of the clocks governing the population movement and the immune systems. In most of the scenarios we test, we observe circadian rhythms would constrain the pace and extent of disease emergence. The three measures (strength, outward transmission and introduction speeds) are highly correlated with each other. For example of the Yushu City, outward transmission speed and introduction speed are correlated with a Pearson’s correlation coefficient of 0.83, and the speeds correlate to strength with coefficients of −0.85 and −0.75, respectively (all have p < 0.05), in simulations with no circadian effect and R0 = 1.5. The relation between the circadian rhythms of the immune system and daily routines in human mobility can affect the pace and extent of infectious disease spreading. Shifting commuting times could mitigate the emergence of outbreaks.


Introduction
Infectious diseases are a major burden to global health and will continue to be that for the foreseeable future. They hinder us in our everyday lives, and thus stops society from operating at well as it could. There is also a looming threat of pandemic of novel and hard-to-cure pathogens. There are both positive and negative factors affecting the long-term threats of infectious diseases-the continuous advances of medicine are a factor that keeps reducing the threat. Factors working in the opposite directions include: the increase of drug resistant pathogens, and urbanization. Since a few years back, most people in the world live in cities [1], and the fraction of city dwellers keep increasing. Bringing people together like a city does, increases the rate people meet close enough to transmit infectious diseases [2]. This makes epidemic thresholds lower, outbreaks more frequent and faster. For these reasons, stopping the epidemic spreading is almost synonymous with stopping them in cities. Disease spreading is a phenomenon that is affected by a multitude of factors. Some of these occurs at different time scales-like evolution of pathogens is usually slower than the microscopic pathogenesis in individuals, that typically is slower than the mixing and interaction of people. In such a case, in epidemic modeling, one can separate the timescales and treat the slower changing feature as static. In other cases, the timescales are similar and then one cannot simplify a model like that. The life in cities follow a 24-hour, circadian rhythm. Many people commute at the same times every day [3][4][5]. People using public transport could be very close together at the rush hours of the day. The immune systems of people also follow a circadian rhythm. Quoting Ref. [6]: "Nearly every arm of the immune response (innate and adaptive) has been reported to oscillate in a circadian manner." The best evidence comes from studies of herpes and influenza virus infection associated with the circadian clock [7]. Typically, the immune system is most effective in the midnight [8,9]. This means that the synchronization of the travel and the immune systems could affect the predictions of disease models and should ideally not be ignored. Still, to the best of our knowledge, no study takes this coupling into account but assume epidemiological characteristics (i.e., the transmission rate) are static over the course of hours and days [10,11]. If the immune system and the travel patterns were the only phenomena following circadian patterns, then taking their coupling into account would be straightforward. One could simulate the disease propagation at a time resolution of days using standard models, only the parameter values would need to be adjusted. However, this is not the case. Many types of infections have a faster development than could be captured by such a coarse time resolution [12]. Thus, for accurate modeling one need to simultaneously take the travel pattern, the strength of the immune systems, and the disease propagation into account. This is what we set out to do in this paper.
As a basis for our study, we use real data of daily travel patterns as our substrate for interpersonal contacts. The data comes from mobility surveys from two counties (roughly the second level administrative division of China, after provinces)-Yushu City and Nong'an County -in China's northeastern Jilin province (Fig 1). From the data, we can see how many people that travelled between subdivisions of the counties at different times. This type of data lends itself well to a meta population type modeling, where one simulates assumes each location to be well-mixed and use the observed population flow to model the disease spreading between the locations [13,14]. We use a susceptible-exposed-infectious-removed (SEIR) model with realistic parameters for influenza [15] and a susceptible-infectious-removed (SIR) model with realistic parameters for varicella [16]. To incorporate the circadian patterns of the immune system, we use a sinusodial correction factor to the transmission rate parameters to model the circadian rhythms of the immune system. With this setup we investigate different scenarios like: the effect of including the circadian effects of the immune system (compared to modeling it as constant through the day) and the effect of shifting the mobility patterns relative to the immune system.

Contact network from mobility data
We use mobility data from a public survey by two Chinese counties (Yushu City and Nong'an County) [5]. These data sets give information about the spatiotemporal co-location of individuals and their movement between locations in the survey. Specifically, our data consist of one-week location records of anonymized cellular phone users across two counties in the Changchun area of China's northeastern Jilin province. These data sets were sampled between August 7, 2017 and August 13, 2017. The two counties consist of a mix of country side and urban areas are subdivided into 49 geographical divisions (Yushu City, with 1.18 million inhabitants, have 27 subdivisions, and Nong'an County, with 0.97 million inhabitants, in 22 divisions), corresponding to an administrative division [5,17].
We divide each county into metapopulations (henceforth locations) corresponding to administrative divisions within a city as reported in the data (Fig 1). In this contact network, an edge denotes the flow of people between two locations and is weighted by the total number of users' movements in a specific hour of the data. In the raw data, the trips are reported as quadruplets (a, i, j, t) meaning that individual a travelled between locations i to j and arrived during the t'th hour of the data. From this data, we derive a series of 168 hourly mobility matrices W t , one for each hour t of the week. The ij-th entry of an hourly mobility matrix at hour t, w ij,t , denotes the number of users traveling from location i to location j at time t (at an hourly resolution).
For our analysis, we need a time-independent measure corresponding to strength [18] (the sum of a node's link weights) for our sequence of mobility matrices. To this end, we simply sum the weights of all the 168 matrices to get the strength

Epidemic model
Basing disease simulations of large data sets of human mobility is challenging. Even though one has access to the trajectories of individuals, running individual based simulations is difficult [19]. How to compensate for missing data and how to scale up the simulations to actual city sizes makes this a questionable approach. Instead researchers rather use metapopulation models where the individuals are assumed to be well mixed within a metapopulation [20] at one segment of time. With metapopulation models, one still takes the flow of people and structures in geography and time into account [21]. Many cases, including ours, one can argue that metapopulation models is modeling at the same level of detail as the original data set, since the sampled trajectories does not cover everyone moving in the area and this incompleteness hard to assess and compensate by adding simulated people. Just like individual-based epidemic models, metapopulation models divide the population into classes (compartments) with respect to the disease and models the dynamics of the individuals within these compartments. We base or simulation on the work of Wesolowski et al. [22], and thus considers four disease compartments: susceptible (individuals don not have, but could get, the disease), exposed (individuals who have got the disease, but cannot yet infect others), infective (individuals who have the disease and can spread it to susceptibles) and recovered (individuals who are immune or deceased and cannot get or spread the disease). For each location i, we denote the number of individuals in the three compartments at time t by S t,i , E t,i , I t,i and R t,i , respectively. In our simulations, time is effectively a discrete variable as the state of the system is reported at an hourly basis.
As epidemic outbreaks emerge, their transmission dynamics can vary much depending on the geographical locations of the initial infection [14,23]. This maybe the most important insight from the recent, data-driven research on epidemics in urban populations [13,14,24]. We chose a metapopulation model, that can capture such effects. In an entirely susceptible population at location i at time t, an outbreak happens with probability: where β t represents the transmission rate at hour t; x k,t and y j,t denote the fraction of the infected and susceptible populations at hour t at location k and location j, respectively, given by x k;t ¼ I k;t N k and y j;t ¼ S j;t N j where N k and N j are the population sizes of location k and location j according to the 2017 Changchun statistical yearbook [17]. Then we simulate a stochastic process introducing infections into completely susceptible metapopulations, in which I j,t+1 is a Bernoulli random variable with probability h(t, j).
After disease introduction, the epidemic model is deterministic: or with an additional exposed class: where σ is the incubation rate at which individuals move from the exposed to the infectious classes. γ denotes the recovery rate at which people move from the infectious to the recovery classes. An important parameter for the discussion is the basic reproduction number, R 0 . It is defined as R 0 = β t /γ and can be interpreted as the expected number of secondary infections if an infectious individual joins into a susceptible metapopulation. This reflects that the population in one location is well mixed. Everyone in a location at a given time has an equal probability to meet any other. The epidemics were simulated over all introduction locations and timings. In our setting, no matter when and where the epidemic begins, the outbreak will always reach the entire system, but the dynamics will vary with the parameter values.
To model the circadian rhythm of the immune system, we modify the transmission rate by a sinusodial function [7,9]: where t is the time of the day, t max is the time of the immune system maximum (both rescaled to the unit interval), β 0 is the average transmission rate over the day, C 2 [0, 1] is the a factor controlling the strength of the circadian effects (C = 0 meaning no effect, C = 1 meaning the at the maximum of the immune system the transmission rate is zero (so the person is perfectly immune to the disease). See the illustration in Fig 2. We include the parameter Δt to model the effect of hypothetically shifting the travel patterns (for example by changing office hours or the time zone). If Δt = 1 hour then the travel patterns would be delayed one hour compared to the observed patterns.

Simulation setup and analysis
In our model, outbreaks are characterized by five parameters-the transmission rate β 0 , the incubation rate σ, the recovery γ, the time of the immune system maximum t max , and the strength of circadian rhythm C. The temporal resolution is one hour. We set epidemiological parameters (β 0 , σ and γ) with respect to two well-studied infectious diseases: First, the 2009 H1N1 influenza pandemic (which we model by an SEIR model)-β 0 = 1/(2.25 × 24), γ = 1/ (3.38 × 24), σ = 1/(2.26 × 24), and R 0 = 1.5 [25]. Second, the 2010 Taiwan varicella by SIR model-β 0 = 1.50/24, γ = 1/(5 × 24) and R 0 = 7.5 [16]. During the time of the survey (mid-August), the time of the sunset (setting the immune system maximum) was 7:00 PM at both the data sets. How to parameterize the circadian effect on the immune system, is still rather unknown in the medical literature. There are some results about herpes and influenza virus infections, respectively, with respect to viral replication volumes in the 24-hour circadian clock [7] pointing at the range 0.5 < C < 1. We use three levels of the strength of the circadian impact on the immune system, spanning the entire range: high C = 1, medium C = 0.5 and low (or, rather, no) circadian effect C = 0.
There are in total 16, 464 possible combinations of parameter values (which can be seen by multiplying the two R 0 s (1.5 and 7.5), two cities, three circadian effect strengths, the 28 introduction timings (at the 8-th hour, 10-th hour, 16-th hour, and 21-th hour of each day of a week [12]), and 49 introduction sites. We make 20 runs for each parameter set, i.e. 329, 280 simulation runs with an initial infection seed and analyze these simulated time-series of the disease prevalence in two ways [12]: 1. We estimate the outward transmission speed Γ 10% as the time (in hours) following an introduction into a location until 10% of the locations have cases.
2. We assess the introduction speed χ as time (in hours) following introduction in a location until the focal location receives its first infection. This time is averaged over all introduction locations and all introduction times.

Ethics statement
This study was approved by the Ethics Review Board of the College of Computer Science and Technology, Jilin University, China.

Effects of including circadian rhythms in epidemic models
Our first investigation concerns the epidemiological importance of introduction timing with and without explicit modeling of the circadian rhythm of the immune system. We measure the time following an introduction of the disease at a location until 10% of the locations have infected people, i.e. the outward transmission speed Γ 10% . As shown in Fig 3, the epidemic spreads slowest in simulations without the circadian effects, fastest in the case of highest circadian dependence of the immune system. The effect is significant but rather weak. The observed pattern holds for both high and low R 0 but is somewhat more pronounced at low R 0 . Furthermore, for each introduction location, we measure the introduction speed χ-the number of hours following introduction in another location until the focal location receives its first infection over diverse introduction timings (see Fig 4). The speed of introduction location to spread

PLOS ONE
Coupling the circadian rhythms of population movement and the immune system in infectious disease modeling pathogens, as captured by Γ 10% , seems to be anti-correlated with the strength of the locations (the total amount of in-and outflow of people), confirming observations in Ref. [12].
In Fig 4, we measure the network structural effect directly. We plot the correlations of three measures: strength (mobility flow into a location), Γ 10% , and χ. The three measures are highly correlated with each other. For example of simulations with no circadian effect and R 0 = 1.5 in the Yushu City, χ and Γ 10% are correlated with a Pearson's correlation coefficient of 0.83, and the speeds correlate to strength with coefficients of −0.85 and −0.75, respectively (all have p < 0.05). Although not related to, or much affected by, the circadian effects, the centrality of locations (quantified by strength) is thus an important determinant of their role in the epidemics.

Effects of shifting time of the travel patterns
In the previous section, we established that there is a non-negligible effect of the circadian rhythms on the results of epidemic modeling. In this section, we look at the effects of the synchronicity between the immune system and the mobility patterns of people. Phrasing this in terms of parameter values, while in the previous section we were varying C and kept Δt constant, in this section we vary Δt and keep C constant. In Figs 5 and 6, we use a box-and-whiskers plots to show the effects of shifting the travel patterns forward or backwards in time. This effect is much bigger than the effect of tuning C. In particular, for the case of strong circadian rhythms Γ 10% gets delayed up to 10%, or seven hours, with one-hour delays of the travel rhythm of people Δt = 1 hour. (Relative changes are summarized in Table 1.) There is an asymmetry in this effect in that setting Δt = −1 gives the reverse effect, but of a lesser magnitude. This is a manifestation of what we mentioned above-that the effect of the coupling between travel patterns and the circadian rhythms of people is complex and trying to capture it by a general fudge factor might not be enough for an accurate modeling. Furthermore, we observe that the effects of different diseases are different depending on whether Δt is positive or negative. The SEIR simulation (with a lower R 0 ) of the influenza simulation makes the negative effects of advancing the travel patterns (negative Δt) much larger. For positive Δt the effect is more similar.

Discussion and conclusion
We have investigated the effects of the coupling of the circadian rhythms of the immune systems and the daily travel patterns of city dwellers. In line with recent studies of urban epidemiology [10,26], we used a metapopulation model based on real survey data that we parameterized to reflect two recent urban disease outbreaks of influenza and varicella [12,22]. We modeled the circadian response of the immune system by multiplying the transmission rate by a sinusodial dependence of the time of day.
With this set-up, we found a small effect (3% difference or less) in including the circadian dependence of the immune systems. On one hand, possibly this effect is smaller than other sources of error, and could thus be neglected. On the other hand, it is possible that for other situations-other diseases, in other populations, at other locations and times of the year-this effect would be larger and less negligible. If possible, we thus recommend modelers to include the effect of the immune system's circadian dependence. Here we study epidemics by shifting mobility with one hour earlier (panels a and b) or later (e and f), in contrast with that without shifting (c and d). For each city, the time until 10% of locations are infected (Γ 10% ) is averaged over epidemics for each of the different introduction times and sites. We compare simulations in settings with the time dependent transmission rate characterized by three levels of circadian effects (high C = 1, medium C = 0.5 and C = 0 with no circadian effect). Epidemic growth varies in different levels of circadian rhythms. Larger Γ 10% is observed with two patterns over R 0 s: First, the time dependent transmission rate than without, increasingly with higher C, when mobility is earlier one hour or not. Second, the time dependent transmission rate than without, increasingly with higher C, when mobility is delayed one hour. These two patterns can be distinguish by the relative transmission speed (Table 1). https://doi.org/10.1371/journal.pone.0234619.g005 We also studied effects of shifting the clock time one hour relative to the sun. These effects could be larger (up to 10%) in the scenarios we simulated. Once again, this is specific to the travel data we studied (two urban regions in northeastern China, in August 2017). Probably, in even larger cities, or cities with yet more pronounced rush hours, these effects could be much larger. By shifting the time zone of a location, one could thus, in theory, slow down the disease propagation considerably, which would give more time for public health authorities to raise awareness of the disease and employ countermeasures. That much said, shifting the times of our daily lives would of course be a political process of great inertia. There is also the caveat that some of our routines are guided by the sun rather than the clock, so shifting the time zone is not the same as shifting the travel patterns. Clearly there is room for more research in the interface of chronobiology, epidemiology and data science (cf. Ref. [27]).  Table 1. Relative transmission speed when shifting the travel patterns. We compare simulations in settings with the time dependent transmission rate characterized by three levels of circadian effects (high C = 1, medium C = 0.5 and C = 0 with no circadian effect) in three scenarios by advancing, delaying mobility one hour or not. Taking mobility without shifting as baseline, we estimate the relative measure of (Λ Δt ) for each scenario to illustrate the potential transmission speed, from which we can distinguish trends with respect to circadian effects (see Figs 5 and 6).

Earlier (Λ −1 )
Delayed (Λ 1 ) This work connects to a larger context of urban planning and social engineering. Earlier studies have found certain locations in cities to be extra important for epidemics [12]. This adds to the urban catalyst theory which explains the dynamics of cities as driven by certain points of interests [28,29]-urban locations that undergoes a rapid development and redirects the flows of people. Such localities do not only have a positive effect on socioeconomic factors [30][31][32], but since they per definition have high strength (many people passing through them), they will also be key locations for the transmission of infections. By taking this into account in city planning, public facilities built to lower disease transmission would be efficient and perhaps cost effective [33,34].