Clocks in the Wild: Entrainment to Natural Light

Entrainment denotes a process of coordinating the internal circadian clock to external rhythmic time-cues (Zeitgeber), mainly light. It is facilitated by stronger Zeitgeber signals and smaller period differences between the internal clock and the external Zeitgeber. The phase of entrainment ψ is a result of this process on the side of the circadian clock. On Earth, the period of the day-night cycle is fixed to 24 h, while the periods of circadian clocks distribute widely due to natural variation within and between species. The strength and duration of light depend locally on season and geographic latitude. Therefore, entrainment characteristics of a circadian clock vary under a local light environment and distribute along geoecological settings. Using conceptual models of circadian clocks, we investigate how local conditions of natural light shape global patterning of entrainment through seasons. This clock-side entrainment paradigm enables us to predict systematic changes in the global distribution of chronotypes.


INTRODUCTION
Many organisms on Earth are subject to a precise 24 h light-dark rhythm. To predict the cycle, most organisms including humans maintain their own timekeeping system, the circadian clock running at a close-to-24 h period without external time cues. An organism's internal clock synchronizes to the daily cycles of the external environment through light as a time-giving cue, called Zeitgeber. The proper coordination of intrinsic rhythms with environmental cycles upon entrainment contributes to a better evolutionary fitness and health (Ouyang et al., 1998;Dodd et al., 2005;Chen et al., 2010). Elucidating entrainment under realistic conditions of intrinsic (clock properties) and extrinsic (light) factors can provide chronobiological insights into geoecological dynamics. Natural light conditions on Earth vary by latitude and seasons. Given the individual variability of circadian periods, either due to natural variation within the same species or due to differences between species, it is likely that entrainment produces a global patterning of circadian phases that are rich in diversity. We address this question by combining mathematical core clock models with geophysical data on seasonal and latitudinal light intensities.

Mathematical Clock Models
At the organismal level, circadian clocks exhibit properties of endogenous limit cycle oscillators, i.e., their rhythms have a well-defined period as well as amplitude under constant conditions and are robust against perturbations. In the mammalian core clock, tissue-level coupling leads to a precise pacemaker (Herzog et al., 2004), even though cellular rhythms of its constituent single neurons are relatively imprecise and partially arrhythmic (Hirata et al., 2019). Perturbations of these organismic rhythms, by light pulses, pharmacological treatments or jet-lag usually decay at relatively fast time scales of a few days (Spoelstra et al., 2004;Kiessling et al., 2010;Ono et al., 2017). Under particular circumstances, external perturbations under constant conditions can even stop the clock (singularity behavior) (Winfree, 1970;Engelmann et al., 1978).
Conceptual models make predictions on the dynamical behavior of circadian clocks based on heuristic principles (e.g., the notion of limit cycle behavior) without considering the detailed molecular machinery that leads to circadian rhythm generation at the cellular level. Such conceptual modeling has a long tradition in chronobiology and helped to understand the response of clock dynamics to external perturbations (Engelmann et al., 1978;Peterson, 1980), entrainment schedules (Wever, 1964), or synchronization processes (Kronauer et al., 1982;Liu et al., 1997), to name a few.
In contrast to this, contextual or detailed biochemical models aim to decipher the molecular regulatory machinery (i.e., the interaction of genes, mRNAs, proteins, post-translational modifications, etc.) that leads to rhythm generation within a biological context, specific to different tissues (Pett et al., 2018) or different organisms such as cyanobacteria, plants, fungi, flies, and mammals (Goldbeter, 1995;Hong et al., 2006;Locke et al., 2006;Kim and Forger, 2012;Woller et al., 2016).

Photoperiodic Entrainment of Circadian Clocks
Unidirectional synchronization of an endogenously oscillating system like the circadian clock to an external periodic forcing signal (Zeitgeber) is termed entrainment. Upon entrainment, a circadian clock with an internal free-running period τ is forced onto the external Zeitgeber period T (period locking). By this means, a stable phase of entrainment ψ emerges as the internal rhythm aligns with the external rhythm (phase locking). The range of entrainment denotes the set of periods that are able to entrain to a Zeitgbeber at a given strength. It usually gets wider (or narrower) with increasing (or decreasing) Zeitgeber strength. The entrainment region within the parameter plane spanned by the period detuning τ − T and Zeitgeber strength Z adopts a tongue like geometry, called Arnold tongue (Arnold, 1987).
Entrainment ranges, phases of entrainment ψ and entrained amplitudes vary systematically not only with respect to Zeitgeber but also with intrinsic clock properties (Aschoff, 1960). Entrainment protocols have therefore been extensively used to compare circadian clocks in different tissues (Abraham et al., 2010), in different genetic background (clock mutants) within the same species (Rémi et al., 2010;Erzberger et al., 2013), or to compare circadian systems of different species (Aschoff and Pohl, 1978).
We recently extended the concept of Arnold tongues to account for seasonality (Schmal et al., 2015). If we consider a circadian system entrained by Zeitgbeber signals of varying period T and photoperiod ̹ (Figure 1A), the region of entrainment adopts an onion shaped geometry (Figure 1B), called Arnold onion. Here, the photoperiod ̹ has been defined as the fraction (in percent) of illumination time (daylength) over the period T of the Zeitgeber cycle. While the entrainment range is largest during equinoctial Zeitgeber cycles (̹ = 50%) it tapers toward the free-running periods under constant darkness (̹ = 0%) and constant light (̹ = 100%) for increasingly extreme photoperiods. Since the free running periods are generally not the same under constant light and constant darkness as predicted by Aschoff 's rule (Aschoff, 1960), the Arnold onion appears to be tilted (compare Figure 1B).
A similar but conceptually different situation occurs if we consider entrainment of organisms with various circadian freerunning periods τ to Zeitgeber signals of fixed period T but varying ̹ (Figure 1C). Diverse circadian periods exist due to natural variation within a population of the same species or traceable clock gene mutations (Konopka and Benzer, 1971). The entrainment range in the free-running period τ and photoperiod ̹ parameter plane adopts an onion shaped geometry similar to the above described case while being tilted in the opposite direction, compare Figures 1B,D.
Here, we extend our analysis to the Earth's natural light conditions from photoperiodic entrainment as originally described in Schmal et al. (2015), where we investigated entrainment to square wave Zeitgeber signals, corresponding to typical laboratory conditions. Using a conceptual model of the circadian clock in combination with a simple celestial mechanicsbased approximation of local light signals, we illustrate circadian entrainment under ecologically relevant conditions. This reveals a global map of how organisms entrain under different seasons and at different latitudes. It also proposes a possibility that the distribution of chronotypes, a quantitative proxy of the phase of entrainment, can change throughout the year and at different positions on the Earth.

The Poincaré Oscillator: A Conceptual Model for Limit Cycle Behavior
Circadian clocks share several properties of limit cycle oscillators at the tissue or organismal level. They cause rhythms of gene expression and physiological processes that are stable in period and amplitude, robust against (small) external perturbations. In order not to restrict ourself to a specific biological context (i.e., a particular organism or cell type) we omit detailed biochemical reactions and exploit a generic amplitude-phase model, whose dynamics is determined by in radial coordinates. Equations (1, 2) describe a Poincaré oscillator (Glass and Mackey, 1988), which exhibits limit cycle oscillations of explicitly defined steady state amplitude A and free-running period τ . Here, r(t) describes the dynamics of the radial component while ϕ(t) describes the angular component.  (4) and (5), applying a square-wave Zeitgeber signal of peak-trough amplitude 0.1 to the x-variable and by using oscillator parameters A = 1 and λ = 0.5 h −1 as described in Schmal et al. (2015).
Radial-relaxation rate λ determines the time-scale of transient dynamics, i.e., the time a Poincaré oscillator needs to return to its steady-state orbit after an amplitude perturbation (e.g., through light pulses) (Figure 2A). While parameter dependencies of oscillator properties such as the amplitude or free-running period are usually intertwined in complex molecular models, the Poincaré oscillator (1, 2) conveniently treats these features as independent parameters. Thus, the impact of internal clock parameters on entrainment and synchronization properties can be studied in a straightforward manner as demonstrated in several studies (Abraham et al., 2010;Granada et al., 2013;Gu et al., 2014;Tokuda et al., 2015;Myung et al., 2018;Schmal et al., 2019). In the following we will represent the circadian clock of an organism by a Poincaré oscillator (1, 2) and study its entrainment under natural light environments, using a simple celestial mechanics derivation.

Approximation of the Perceived Light Intensity on Earth
For a sessile observer at a given position on Earth, the Sun obeys an apparent movement across the celestial sphere. This movement leads to changes in light intensity, predictable from celestial mechanics. While the rotation of the Earth around its own axis leads to daily changes of light intensity, the tilt ε (obliquity of the ecliptic) of its rotation axis with respect to the plane that is spanned by the Earth's orbit around the Sun (the ecliptic) leads to seasonal variations ( Figure 2B).
A simple celestial mechanics based model for the light intensity I, perceived by an organism at latitude φ and on a given day N of the solar year, assumes that the declination δ of the Sun (see Figure 2B) varies sinusoidally across the year, thus neglecting the elliptical motion of the Earth around the Sun (Khavrus and Shelevytsky, 2012). Here, 365.2422 denotes the mean number of days in a year, i.e., the mean orbital period of the Earth. The ecliptic ε is subject to fluctuations at timescales that are extremely long (Milanković cycles) compared to any organism's life expectancy such that we consider it constant. It currently takes a value of about ε = 23.44 • . In addition to celestial constellations between the Earth and the Sun, scattering and absorption of light within the atmosphere further affects the solar radiation that reaches an observer at the Earth's surface. The extent of absorption, reflection and scattering depends on the path length within the Earth's atmosphere, called the optical air mass (AM), which in turn depends on latitude as well as season (Khavrus and Shelevytsky, 2012).
The orientation of the light-perceiving surface (e.g., the retina or plant leaves) with respect to the direction of the Sun affects the solar irradiance as well. Combining all of these considerations, we can give a closed form expression for the solar irradiance perceived at any position on Earth in dependence on latitude φ, time of year N, time of day and the orientation toward the Sun, as given by Equation (11) in section 4. Figure 3 illustrates dependencies of the solar irradiance (A), daylength (B), total daily insolation (C) and standard deviation of the daily solar irradiance (D) on geographical latitudes and seasons for sun rays reaching a horizontally oriented surface. While the daylength (see Equation 9) is constantly 12 h throughout the year at the equator as well as at vernal and autumnal equinoxes for all latitudes φ, it adopts the highest variation between different latitudes at summer and winter solstices ( Figure 3B). The total daily insolation (i.e., the summed irradiance per day) is determined jointly by daylength and irradiance. At summer solstice, for example, the daily insolation adopts similar values at the equator and the arctic pole since a high irradiance is perceived over a shorter daylength near the equator, while a low irradiance is perceived over a long 24 h photoperiod at the arctic pole ( Figure 3C).
Critical for circadian entrainment is the variation of a Zeitgeber signal, which can be quantified by the standard deviation of the solar irradiance over one day ( Figure 3D). Its highest values are reached during the equinoxes at the equator and during summer solstice at temperate latitudes, while the lowest variation is found during polar night. These local geographical Zeitgeber variations hint at global variations of the entrainment range. (B) Illustration of the celestial constellation between the Earth and Sun during summer in the northern hemisphere. Here, ε denotes the obliquity of the ecliptic while δ is the declination of Sun. (C) Global distribution of daylengths on a summer day in the northern hemisphere, with geographical locations of representative cities used in this study: Kampala (equator), Taipei (sub-tropics), Berlin (temperate), and Ny-Alesund (above arctic circle). Note that the irradiance is lower in the northern extreme despite the longer daylength. LD, long daylength; 12:12, equinox; SD, short daylength. (D) Combining a conceptual oscillator model with a celestial mechanics based approximation of the light intensity, we compute entrainment ranges and phases of entrainment ψ under different seasons and latitudes. By this we can study impact of season and latitude on the distribution of chronotypes.  (9) in section 4. The daylength corresponds to the distance between sunrise and sunset as exemplified in (A) for three different latitudes. (C) Daily insolation for different latitudes and times of year. The daily insolation is defined as the daily sum over the solar irradiance as depicted in (A) for three exemplary latitudes. Please note that the solar irradiance in (A) and solar insolation in (C) have been calculated for sun rays reaching a horizontal plane. The direction of the surface (e.g., of plant leaves, the retina, etc.) has an impact on the solar irradiance perceived. (D) Standard deviation of the solar irradiance, calculated over a complete day (24 h), for different latitudes and times of year.

Entrainment Ranges Vary Systematically With Season and Latitude
Circadian clock properties vary among different species and even among individuals of the same species as a result of natural (genetic and epigenetic) variation. This eventually translates into variations of entrainment properties such as the phase of entrainment ψ, even within the same species that receive the same entrainment cues. Such emergence of different "temporal phenotypes" can be represented as a spread of entrainment phases ψ in a population of individuals of same species. One such representation is commonly known as chronotypes (Pfeffer et al., 2015). In the following we investigate how organisms with different internal clock properties as described by the conceptual oscillator model (1, 2) entrain to natural light signals at different seasons and geographical locations. Figure 4 depicts ranges of internal free-running periods τ , for which entrainment upon natural light signals given by Equation (11) can be observed, at different times of the year N and at four different representative latitudes φ, namely in Kampala close to the equator (A), Taipei in the sub-tropics (B), Berlin at a temperate latitude (C), and Ny-Alesund (D) above the arctic circle. Here, we only consider the effect of free-running period τ . Intrinsic clock properties other than τ are held constant at A = 6 and λ = 0.1 h −1 throughout all simulations. Phases of entrainment ψ have been color-coded within the regions of entrainment. Solar noon, the highest position of the Sun at the celestial sphere, is used as a reference point such that ψ = 0 h is adopted if the acrophase of the circadian clock oscillations (given by x(t)) coincides with solar noon.
Close to the equator, there is little per-day change of daylength, insolation and variation (standard deviation; SD) throughout the year in comparison with higher latitudes (Figures 3B,D). Thus, entrainment ranges remain relatively constant throughout the year (Figure 4A). Within these small variations, the largest range of entrainment occurs twice a year at vernal and autumnal equinoxes while it becomes the smallest around summer and winter solstices. This is in line with the profiles of solar insolation and SD of daily irradiance, peaking twice per year close to the equator (Figures 3C,D).
At higher latitudes, this half-year period in insolation and SD of daily irradiance is absent, leading to entrainment ranges that are maximal during summer solstice and minimal at winter solstice (Figures 4B,C). The seasonal effect is stronger as the latitude φ is higher. At high latitudes, organisms with a circadian clock whose free running period τ substantially deviates from 24 h might be able to well entrain only during summer but could fail to entrain during winter season, as described for European hamsters (Cricetus cricetus) kept under natural lighting conditions (Monecke and Wollnik, 2005).
At latitudes above the arctic circles, the Sun never rises during polar night and never sets during midnight sun ( Figure 3B). Since our approximation of light intensity ignores diffraction effects by the atmosphere, organisms are unable to entrain during the polar night ( Figure 4D). Between the vernal and autumnal equinoxes, entrainment is generally possible even during a 24 h midsummer day, that is because the solar altitude above the horizon still exhibits diurnal variations. In line with these results, reindeers (Rangifer tarandus) in Spitzbergen do not show diurnal rhythmicity during winter but some residual diurnal activity throughout the summer season (Arnold et al., 2018). Other species such as Svalbard ptarmigans (Lagopus mutus hyperboreus) show behavioral arrhyhtmicity during polar night and midnight sun (Williams et al., 2015).
So far, we studied the range of entrainment throughout the year at four different exemplary latitudes (Figure 4). In Figure 5, we investigate the effect of varying latitudes at three different times of the year. Since the daylength is approximately 12 h at all latitudes ( Figure 3B) and the dependency of the Zeitgeber intensity (SD of irradiance) on latitude is symmetric for the northern and southern hemispheres, the entrainment region during vernal equinox (and analogously at autumnal equinox) is perfectly symmetric (Figure 5A). The largest entrainment range is adopted at the equator (φ = 0 • ) and tapers toward the polar regions. At summer solstice in the northern hemisphere, an asymmetric entrainment region can be observed ( Figure 5B). Due to polar night no entrainment is possible above the antarctic circle during this season, but it is possible above the arctic circle. The largest entrainment range is adopted around φ ≈ 30 • as the SD of daily irradiance becomes maximal (compare Figure 3D). Since summer solstice in the northern hemisphere corresponds to winter solstice on the southern hemisphere and vice versa, the entrainment region at winter solstice in the northern hemisphere Frontiers in Physiology | www.frontiersin.org is a mirror image of the summer solstice entrainment region ( Figure 5B) with respect to the equator ( Figure 5C).

Tracking of Dusk and Dawn Depends on Free-Running Period τ
From an evolutionary point of view, two aspects of circadian entrainment need to be considered. On the one hand, it has to be assured that the circadian program entrains to the environmental Zeitgeber signals at all, which is related to the range of entrainment. Of equal importance, on the other hand, is the existence of a proper phase of entrainment ψ, ensuring that physiological processes are executed at most beneficial times around the solar day.
An internal oscillator's activity onset phase can be locked to different phases around the Zeitgeber reference. This is because ψ depends on Zeitgeber properties as well as internal oscillator properties such as the free-running period and amplitude. Figure 6A depicts the evolution of ψ throughout the year for three different free running periods τ (23, 24, and 25 h) at the latitude of Taipei (φ = 25.03 • ), corresponding to vertical cross-sections through the entrainment regions as depicted in Figure 4B. While circadian oscillators with an internal period of 24 h (bold line, Figure 6A) appear to be phase locked in parallel to solar noon throughout the year, ψ for clocks with an intrinsic period τ longer than 24 h (dotted line, Figure 6A) occurs earlier during summer days compared to winter and has its earliest value at summer solstice. The opposite is true for circadian clocks with τ < 24 h, where ψ occurs later during summer in comparison to winter days (dashed line, Figure 6A). These results quantitatively confirm a similar qualitative description by Pittendrigh (Pittendrigh, 1988). The biological importance can be interpreted within the framework of Aschoff 's rule which states that day-active animals typically have a period longer than 24 h while night-active animals have a period shorter than 24 h (Aschoff, 1958). Thus, ψ of day active animals (τ > 24 h) follows the earlier sunrise as photoperiods increase in spring and summer, while nightactive animals (τ < 24 h) rather follow a later dusk (onset of night) as photoperiods increase (Pittendrigh, 1988) (compare Figure 6A).

Phases of Entrainment ψ Vary Systematically With Season and Latitude
The range of entrainment depends on Zeitgeber strength and photoperiod. It has been shown that, within the range of period detunings τ − T that allow for entrainment at a given fixed Zeitgeber strength and photoperiod, the phase of entrainment ψ varies by about 180 • (Wever, 1964;Granada et al., 2013;Schmal et al., 2015). This translates to a higher sensitivity of ψ with respect to period detunings τ − T for small ranges of entrainment. Figure 6B shows the dependency of entrainment ranges on latitude, exemplarily for March equinox as well as for northern hemisphere's summer and winter solstices. It can be seen that the entrainment range gradually decreases as the latitude gets higher. From this, we can infer that on March equinox, ψ has an increased sensitivity to τ variations as the latitude increases ( Figure 6C). This corresponds to horizontal cross-sections at March equinox (N = 79) in Figures 4A-D or to horizontal crosssections at the latitudes of Kampala (equator), Taipei, Berlin and Ny-Alesund (above arctic circle) in Figure 5A, respectively. Likewise, the phase of entrainment ψ has the highest sensitivity to τ variations (T = 24 h) at winter solstice, the lowest sensitivity at summer solstice, and an intermediate sensitivity at March equinox. This is illustrated for the latitude of Berlin (ψ = 52.51 • ) in Figure 6D, which corresponds to horizontal cross-sections at the respective times of the year in Figure 4C. A similarly strong dependency of the slope of ψ on the Zeitgeber strength has been described for the entrainment of lizards to temperature cycles (Hoffmann, 1969).
The sensitivity of ψ with respect to τ variations can be used to explain the behavior of chronotypes across different seasons and latitudes as investigated in the following section.

Season and Latitude Affect the Spread of Chronotypes
In our analysis of generic clocks, ψ is referenced to the point of midday. This referencing is arbitrary and is easily translatable to chronotypes. Human chronotype reflects the phase of entrainment, and is commonly referenced to the mid-sleep phase. Variations in the free-running period τ have been shown to translate into changes of the phase of entrainment ψ (Winfree, 2001;Bordyugov et al., 2015). Genetically modified strains of Neurospora crassa with relatively shorter and longer τ , respectively, exhibit early and late values of ψ under light entrainment (Rémi et al., 2010). Likewise, the free-running period of fibroblasts is positively correlated with phase ψ under temperature entrainment (Brown et al., 2008). This has also been seen in weakly coupled Bmal1 oscillators in the suprachiasmatic nucleus explants cultures (Myung et al., 2012).
Along these lines, it has been hypothesized that small variations of free-running periods τ within a population of the same species can lead to the emergence of different chronotypes as variations in the phase of entrainment (Phillips et al., 2010;Granada et al., 2013). Chronotypes based on sleepwake time questionnaires can be a proxy to the phase of entrainment ψ, although in higher organisms other factors such as a homeostatic sleep drive interferes with pure circadian effects. "Strong" clocks with a small entrainment range and a large sensitivity of ψ will map a given spread of freerunning periods τ (described by standard deviation σ τ ) into a relatively large spread of chronotypes (σ ψ ) when compared to weak clocks with a large entrainment range and a small sensitivity of ψ.
Under a forced desynchronization protocol, it has been shown that human females adopt a distribution of internal free-running periods of mean µ τ = 24.09 h and standard deviation σ τ = 0.2 h (Duffy et al., 2011). Taking this data as an illustrative example, we sampled 1, 000 periods from a normal distribution with this mean and standard deviation (see Figure 7A) and computed the corresponding entrainment phases ψ at different seasons and latitudes. This task corresponds to a convolution of the period distribution with the phase of entrainment profiles as determined in Figures 4, 5, 6C,D. The effect of season on the distribution of chronotypes in our model has been exemplarily calculated for the latitude of Berlin ( Figure 7B). In summer (Jun 20th), a strong Zeitgeber signal (large SD of irradiance) leads to a wide entrainment range and a small sensitivity of ψ with respect to τ variations. In other words, the distribution of chronotypes is relatively narrow (σ ψ ). In winter (Jan 1st), light becomes a relatively weak Zeitgeber, leading to a narrow entrainment range and a high sensitivity of ψ. This ultimately leads to a larger spread in chronotypes (σ ψ ) in winter compared to summer.
If the same population would travel during winter season (Jan 1st) from the latitude of Berlin to latitudes substantially closer to the equator, for example Taipei, the population would find the distribution of chronotypes to become narrower ( Figure 7C). That is due to the fact that the Zeitgeber is stronger toward the equator in winter (Figure 3D), eventually leading to a wider entrainment range and a smaller sensitivity of ψ with respect to τ variations (see Figure 5C in the northern hemisphere).
These theoretical observations imply that the spread in chronotypes (σ ψ ) varies throughout the year. Examples in Figure 7D have been calculated, based on the τ distribution in Figure 7A, the temperate latitude of Berlin and the latitude of Taipei closer to the equator. It can be seen that the seasonal effect is relatively smaller at the latitude of Taipei than in Berlin. At the higher latitude, a standard deviation of σ τ = 0.2 h in free running periods leads to a standard deviation of up to almost σ ψ = 3 h in the distribution of ψ during winter solstice.
The combined geographical distributions of entrainment ranges throughout the seasons can be mapped on the globe, using the scheme outlined in Figure 2. This visualizes which regions on the globe are easier to entrain to local day-night conditions (Figure 8).

DISCUSSION
The major finding of the present study is the organization of circadian phases in dependence of the geographical and seasonal context and individual clock properties (personal chronotypes). This organization can be assessed by a construct we call Arnold onion (Schmal et al., 2015). The clock's ability to entrain to natural light-dark cycles, i.e., the entrainment range, reorganizes across the year, which can be most intuitively visualized for different latitudes on the globe (Figure 8). By combining a generic clock model with a geophysical model of local light intensities, we could examine season-and latitude-dependent entrainment properties for circadian systems with different clock properties.
The range of internal periods τ , capable of entraining to daily rhythms of the outside world, remains relatively stable close to the equator but shows strong seasonal variations at high latitudes, having a narrower entrainment range during winter (Figure 4). This theoretical prediction is consistent with experimental studies, showing stronger entrainment during summer compared to winter where entrainment might even fail at higher latitudes (Monecke and Wollnik, 2005;Arnold et al., 2018).
For a given circadian phenotype with a well-defined internal period τ and given amplitude, our model predicts systematic variations of the phase of entrainment ψ across the seasons. Similar to earlier studies (Pittendrigh and Daan, 1976b;Pittendrigh, 1988), dusk tracking is predicted for clocks shorter than a day (τ < 24 h) while dawn tracking occurs for clocks with τ > 24 h ( Figure 6A). Such seasonal variations in the timing of daily activity have been described for entrainment under natural conditions in several species including birds and rodents (Aschoff and Wever, 1962;Daan and Aschoff, 1975). In some species, other entrainment cues such as social synchronization and temperature can overtake photic entrainment (Robbers et al., 2015;Fuchikawa et al., 2016).
Our theoretical examination predicts that there can be latitudinal and seasonal changes in the global distribution of entrainment phases ψ, which can be reflected in chronotype measures. Season and latitude can affect the mean (Figure 6) and standard deviation (Figure 7) of the chronotype distribution as the average internal free-running periods (µ τ ) deviate from 24 h. In humans, distributions of the midpoint of sleep on free days (MSF), a common measure of chronotypes, are different between male and female populations and vary also with age (Foster and Roenneberg, 2008). Humans adopt their latest chronotype roughly at the age of 20, and become earlier again with increasing age (Roenneberg et al., 2004). Reliable estimates of population  Figure 2, where the entrainment ranges correspond to those depicted in Figure 6B.
variance through seasonal and latitudinal variations are difficult to obtain amidst these internal heterogeneities, in addition to confounding cultural influences amongst studies including societies of different lifestyles (Pilz et al., 2018). We expect that finer details of chronotype diversity can be predicted in defined subpopulations of age and gender when more precise physiological mechanisms of chronotypes are known. It should be noted that our model has been developed to exploit general principles of circadian entrainment and has not been finetuned to account for human entrainment characteristics. Human chronotypes are commonly associated with the "preferred sleepwake schedule, " and the sleep-wake cycle is not linearly related to the underlying circadian cycle. Along with the circadian regulation, human sleep patterns are simultaneously affected by homeostatic processes (Phillips et al., 2010;Skeldon et al., 2016). Previous conceptual models of human circadian entrainment more specifically elaborated on light-sensing mechanisms and mimicked human phase-response behavior Jewett et al., 1999). It would be interesting to investigate in future studies how these human-specific constraints affect the season and latitude dependent responses of chronotype distributions. However, the occurrence of chronotypes is not exclusive to humans and has been described for other species as well (Pfeffer et al., 2015). Our theoretical predictions could be useful for further studies on the effect of detailed Zeitgeber characteristics upon the distribution of chronotypes or entrainment ranges within a controlled setting, using populations of well-studied model organisms at different geographical locations and seasons.
Parameters other than free-running periods τ such as the amplitude have been proposed to further affect entrainment characteristics (Lakin-Thomas et al., 1991). Throughout this study we held amplitudes (A) and radial relaxation rates (λ) constant and focused on the impact of τ on entrainment properties. Analogous to previous studies (Abraham et al., 2010;Bordyugov et al., 2011) we find that decreasing amplitudes and relaxation rates broaden the entrainment range with respect to τ variations throughout the year (Figures S1, S2). Therefore, decreasing amplitudes and relaxation rates will lead to a reduced spread of entrainment phases. A careful interpretation of entrainment characteristics among different circadian phenotypes (organisms, tissues, mutants, etc.) will require inclusion of a thorough analysis of all relevant internal oscillator parameters.
Properties of the circadian clock have been shown to be plastic and diverse within the same species and across an organism's domicile environment and life span. Chronic conditions of external entrainment can alter the intrinsic period of a clock (Pittendrigh, 1960;Pittendrigh and Daan, 1976a). Across several organisms, including mice and human, it has been shown that entrainment at different photoperiods (Ciarleglio et al., 2011;Myung et al., 2015) and different Zeitgeber periods T (Scheer et al., 2007;Azzi et al., 2014Azzi et al., , 2017 reversibly or irreversibly changes the free-running period τ , a phenomenon termed aftereffect or imprinting, respectively. Additionally, it has been found that properties of the clock change with increasing age (Pittendrigh and Daan, 1974). Although our conceptual oscillator model does not include light-or age-dependent parameter plasticities, Arnold onions can be used to predict changes of entrainment properties upon plasticities in the circadian clock behavior by highlighting subpopulations of τ . Assuming that plasticities occur either under relatively short (after-effects) or long (aging) time scales in comparison to seasonal changes of light intensity, one can predict plasticity dependent changes of ψ in a latitude and season dependent manner by tracing relevant τ subpopulations expected from plasticities within an Arnold onion as depicted in Figures 4, 5.
Latitudinal clines of oscillator properties have been observed in various organisms. Plant leaf movements and oviposition or eclosion rhythms in Drosophila show shorter free-running periods in northern compared to southern strains (Mayer, 1966;Hut et al., 2013). The free-running period in tomato decelerates during domestication, most probably to adapt to longer summer days after its transferal from the equator to northern latitudes (Müller et al., 2015). Our conceptual clock model, extended to account for photoperiodic entrainment under natural conditions, can help predict how properties of the internal clock affect the phase of entrainment and thus can help to untangle mechanisms underlying latitudinal clines.
Several approximations have been done with respect to the calculated light intensity that entrains the conceptual clock model. We assumed that a given organism perceives light similar to a horizontal plane. Young sunflowers follow the path of the sun across the celestial sphere during the day, driven by antiphasic patterns of stem elongation (Vandenbrink et al., 2014;Atamian et al., 2016). By this means, a larger amount of daily insolation can be perceived as compared to a horizontal leaf orientation (Shell et al., 1974; Figure S3A). The variation of the light intensity per day, on the other hand, increases with solar tracking (Figure S3B) which facilitates the entrainment of the circadian system especially during winter ( Figure S3C). Active orientation toward the Sun can therefore lead to broader entrainment ranges and a narrowed distribution of chronotypes. Future studies might investigate the impact of orientation toward the light source on entrainment characteristics in more detail.
Light intensity as perceived by an organism on Earth has been calculated as the sum over the whole spectrum of sunlight throughout this study. It is known that some species harbor different wavelength of natural light differentially (Roenneberg and Foster, 1997;Fankhauser and Staiger, 2002). Changes of spectral composition throughout the course of the day have been proposed to serve as an entrainment cue (Walmsley et al., 2015). In modern societies, anthropogenic artificial light interferes with the natural light-dark profiles. Self-selected light during the night has been described to affect the sleep-wake behavior in humans (Skeldon et al., 2017;Pilz et al., 2018). In addition, changing weather and movement between differentially shaded habitats affect the amount of light that a given organism perceives. An extension of our current study could investigate how above described effects on the effective Zeitgeber strength alters our computationally estimated entrainment properties under natural conditions.
The entrainment phase depends on the detuning (τ −T) of the clock and Zeitgeber period, on the ratio of the Zeitgeber strength to clock amplitudes, and on relaxation rates (Wever, 1964;Granada et al., 2013;Bordyugov et al., 2015;Schmal et al., 2015). These generic dependencies theoretically allow us to predict entrainment properties as a function of season and latitude. It is therefore expected that our modeling approach will provide a framework to analyze empirical data on chronotype distributions under different geographical and ecological contexts.

Conceptual Amplitude-Phase Model
Using the identities x = r cos(ϕ) and y = r sin(ϕ), we can rewrite Equations (1, 2) of the main text as in Cartesian coordinates. Similar to previously published studies Schmal et al., 2015), we choose that light affects the circadian clock by adding the light intensity I(t) at a given time t to the x direction on the phase plane such that Equation (4) reads as dx dt = λx(A − r) − 2π τ y + I(t).
Using a common approximation for the optical air mass (AM) (Kasten and Young, 1989) AM = cos(z) + 0.50572(96.07995 − z) −1.6364 −1 (10) where z is the zenith angle in degrees, we can write the total irradiance I that reaches a plane horizontally-aligned to the Earth's surface as I = 1.353 kW m 2 × 1.1 × 0.7 AM 0.678 × cos(z) if cos(z) > 0 0 otherwise (11) where 1.353 kW m 2 is the average solar constant and factor 1.1 is for diffuse irradiance (Khavrus and Shelevytsky, 2012).

Numerics
Equations (4, 5) in combination with Equation (11) have been numerically solved using the SCIentificPYthon function odeint at a step size of t = 0.01 h. Entrainment ranges as well as phases of entrainment have been determined as described earlier (Schmal et al., 2015). The highest position of the Sun at the celestial sphere served as a reference point of the Zeitgeber cycle such that the phase of entrainment ψ ∈ [−12h, 12h] is defined as the distance between the acrophase of x(t) oscillations and solar noon.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
CS designed the study, performed the analysis, visualized the data and drafted the manuscript. JM designed the study, made a figure, and drafted the manuscript. All authors interpreted the data and wrote the manuscript. All authors read and approved the submitted version.