Monthly data of subjects from 15 countries, who had positive specimens for influenza viruses from January 2000 to December 2015 were obtained from the database of the WHO Global Influenza Programme (World Health Organization, 2015). Monthly data of multivariate El Niño-southern Oscillation Index (MEI) were downloaded from the website of National Oceanic and Atmospheric Administration (National Oceanic and Atmospheric Administration, 2015a). The MEI, which was computed from sea-level pressure, zonal and meridional components of the surface wind, sea surface temperature, surface air temperature, and total cloudiness fraction of the sky of the South Pacific Ocean (Wolter and Timlin, 1998), was used for analysis because it provides a single index of both atmospheric and oceanic components, rather than Southern Oscillation Index, which measures the atmospheric component, or Sea Surface Temperature, which measures the oceanic component.

Subjects were grouped into three regions in the north hemisphere, but into two regions in the southern hemisphere. Regional groupings, which were of contiguous countries as much as possible are North America, which included United States of America and Canada; Nordic/Europe, which included Spain, Portugal, France, Germany, Denmark, Norway and Sweden; Asia, which included China and Japan; South America, which included Argentina and Chile; and Oceania, which included Australia and New Zealand. Angular histograms of monthly distribution of influenza cases were drawn to determine the seasonality of influenza in each region of study. The Rayleigh's test for non-uniformity of circular distribution was performed to determine significance.

ENSO time series has been shown to be non-linear and (Elsner and Tsonis, 1993; Mukhin et al., 2015) and non-stationary (Timmermann et al., 2003). Linear and stationary methods are not appropriate to analyze time series with non-linearities. Presence of non-linearities in atmospheric and oceanic components was determined using the third order moment method (Barnett, 2007). The highest p-value for rejecting the null hypothesis of linearity and stationarity was < 0.001. The synchrosqueezing transform (Daubechies et al., 2011), a time-frequency reassignment algorithm, was applied to determine the spectra of time-varying oscillatory components of ENSO and influenza time series.

Non-linear, non-stationary signals which contain single time-varying amplitude and phases can be modeled as s(t) = A(t)cos[2πϕ(t)] (Boashash, 1992), where A(t) is time varying amplitudes, and ϕ(t) is time varying phases. Natural signals like the ENSO, which have multiple oscillatory components, can be modeled as sum of oscillatory components and noise, s(t)=∑k=1kA(k)cos[2πϕ(k)]+η(t). Such signals must, however, be decomposed to monocomponents before physically meaningful instantaneous phases (ϕ_i(t)), frequencies (A_i(t)), and amplitudes (ω_i(t)) can be generated. Empirical mode decomposition (Huang et al., 1998) separates fast and slow oscillatory components of non-linear and non-stationary signals into intrinsic mode functions (IMFs), unlike harmonic analysis, which decomposes signals into components that are integer multiples of the fundamental frequency. To align the oscillatory components of ENSO and influenza time series and reduce mode mixing, noise assisted multivariate empirical mode decomposition (Rehman and Mandic, 2010a) which applies multivariate empirical mode decomposition algorithm (Rehman and Mandic, 2010a), was used. The algorithm of multivariate EMD (Rehman and Mandic, 2010a) is as follows:

The synchrosqueezing transform was applied to determine the oscillatory components in each IMF. The IMFs with annual oscillatory component were selected for further analyses to determine time-dependent intrinsic correlations, phase space, recurrence and joint recurrence plots, synchronization, and direction of causality.

Time dependent intrinsic cross correlation measures the cross correlation (TDIC) of oscillatory components of non-linear, non-stationary time series (Chen and Huang, 2010). TDIC is a three stage method of empirical mode decomposition, Hilbert transform of IMFs to instantaneous amplitudes and frequencies, and TDIC of instantaneous amplitudes and frequencies (Huang and Schmitt, 2014). The Hilbert transform of IMF₂ was done to generate analytical signal from which instantaneous amplitude (ω_i(t)), phase (ϕ_i(t)), and frequencies (A_i(t)) were extracted (Huang et al., 1998). TDIC of instantaneous frequencies and amplitudes of ENSO and influenza for IMF2 were performed for each region.

The phase space of time series can be reconstructed from scalar data using the method of delayed embedding vector coordinates, which derives from the Taken's theorem (Takens, 1981). Scalar time series {x_i, i = 1, 2, …} is converted to its vector form {X_i = x_i, x_i+L, x_i+2L, …, x_i+(m−1)L}, where m is embedding dimension, and L is the delay or lag (Packard et al., 1980). Theoretically the embedding dimension should be m ≥ 2⌈d⌉+1, where d is the fractal dimension, and ⌈d⌉ is the lowest integer greater than d to preserve the dynamical properties of the attractor. To characterize the phase spaces of ENSO and of influenza in five geographical regions, embedding dimensions m were calculated for IMF2s using Cao's algorithm (Cao, 1997), and lags τ were estimated using mutual information algorithm (Martinerie et al., 1992). The Lyapunov exponent, which measures the sensitivity of dynamic systems to small changes in initial conditions (Wolf et al., 1985), was used to detect the presence of chaotic attractors.

Differentiable graphs have topological dimension d, while non-differentiable graphs have fractal dimension, which have values between the topological dimension d and d+1 (Gneiting et al., 2012). Fractal time series, which are not differentiable, can be characterized using the Hurst exponent (H), which measures long range dependence or memory. Hurst exponent, which is expressed as x_(t) = a^Hx(at), is related to fractal dimension D as D = 2−H (Mandelbroit and van Ness, 1968). While the Hurst exponent of monofractal time series is independent of time and space, multifractal time series have time varying Hurst exponents (Ihlen, 2012). Detrended fluctuation analysis (DFA) algorithm estimates Hurst exponent of long-range temporal correlations, while multifractal detrended fluctuation analysis (Kantelhardt et al., 2002) detects multiple scaling of fractal time series. To characterize the geometry of ENSO and influenza attractors, multifractal detrended fluctuation analysis (Ihlen, 2012) of annual oscillatory components of ENSO and influenza in five geographical regions were performed.

The recurrence plot of dynamic systems is generated from recurrence matrix, which is defined as Ri,j(ϵ)=Θ(ϵ-∥x→i-x→j∥),i,j=1,…,N, where N is the number of measured points x→i, ϵ is the threshold distance, Θ(.) is the Heaviside function, and ∥·∥ is the norm (Marwan et al., 2007). When a trajectory in phase space, xi→i=1N, visits a state xi→≈xj→ the recurrence R_{i, j} is 1, but when the trajectory visits a state xi→≉xj→ the recurrence R_{i, j} is 0. Since trajectories do not return to exact regions of phase space xi→xj→, the threshold ϵ defines the region of phase that is the neighborhood (Marwan et al., 2007). Recurrence plots were produced to characterize the recurrence of annual oscillatory components of ENSO and influenza.

Chaos-periodic and periodic-chaos transitions of non-linear dynamics systems have been detected using metrics of recurrence quantification analysis, which are computed for windows or sub-matrices along the line of identity of recurrence plots (Marwan, 2011). These metrics, which complement visual interpretation of recurrence plots, are measures of complexities that are derived from diagonal and vertical structures of recurrence plots (Marwan et al., 2002, 2007). Determinism (DET), which is the ratio of recurrence points that form diagonal lines of at least l_min and all recurrence plots, and divergence (DIV), which is the inverse of the longest diagonal line L_max (Marwan, 2011) are metrics of diagonal structures, which show maximal values at periodic-chaos and chaos-periodic transitions (Marwan et al., 2007; Marwan, 2011). To characterize the regimes in ENSO and influenza dynamics, divergence and determinism were calculated for window sizes of 36 months, which is approximately the median period of recurrence of El Niño events. To determine the statistical significance of the metrics, the distributions of determinism and divergence for 5,000 surrogate time series of influenza and ENSO were determined and compared with the distribution of 5,000 bootstrap samples of determinism and divergence of annual oscillatory components of influenza and ENSO. Each geographical region was tested separately.

Joint recurrence plots were computed to compared recurrences of ENSO and influenza in phase space. Metrics of network analysis, which maps time series to networks with distinct topological characteristics (Campanharo et al., 2011), have been used to detect the presence of regimes in dynamic systems. These metrics are computed from adjacency matrix, which is defined as A_pq(ϵ) = R_pq(ϵ) − δ_pq, and has correspondence to the recurrence matrix. Global clustering coefficient (C), which is the mean of the local clustering coefficient of all vertices, and average path length (L), which is the mean geodesic distance between all pairs of vertices where paths exist (Donner et al., 2011; Zou et al., 2012), are network metrics that have been used to detect the presence of regimes and transitions from periodic to chaotic regimes and vice versa. To characterize the regimes in joint recurrence of ENSO and influenza using geometrical properties, global clustering and average path lengths were calculated for 36 months windows of joint recurrence networks of annual oscillatory components of ENSO and influenza. To determine statistical significance of the network metrics, the distribution of average path length and global clustering of joint recurrence network of 5,000 surrogate time series of ENSO and influenza were compared with 5,000 bootstrap samples of joint recurrence network of annual oscillatory components of ENSO and influenza. To determine the effect of chaos on severity of influenza epidemics, the five most severe peaks of influenza epidemics between 2000 and 2015 were mapped to joint recurrence network of annual oscillatory components of ENSO and influenza, Figure 4.

The phase space of periodic trajectories increases by 2π for each rotation, or ∥x(t+T)−x(t)∥ = 0, where T is the period (Romano et al., 2005). Although chaotic trajectories shift between unstable periodic orbits, recurrence occurs when each unstable periodic orbit is completed, which is an increment of 2π (Boccaletti et al., 2002). The unstable periodic orbits of chaotic trajectories are locked when they are in phase synchrony, which implies n_i : m_i relationship of their recurrence plot phases, where _i is time index (Romano et al., 2005). The statistical measure is how often ϕ₁ and ϕ₂ increase by 2π or its multiples within time interval τ (Romano et al., 2005). Phases of recurrence plot of chaotic trajectories, therefore, increase by 2πk within time interval τ. Phase synchrony is measured by comparing the probability P^(ϵ)(τ) that the trajectory returns to ϵ-neighborhood of a previous point on the trajectory. The cross correlation between P₁(τ) and P₂(τ) is defined as CPR=〈P-1(τ)P-2(τ)〉/(σ1σ2). When the trajectories are in phase synchrony the probability of recurrence are maximal at the same time, and CPR ≈ 1. To determine if the trajectories of ENSO and influenza are synchronized, the CPR index was calculated for the five regions.

Two dynamic systems are topologically equivalent or homeomorphic, if there is a function that maps every point on one trajectory with the other. If time series x(t) and y(t) are from the same n-dimensional manifold, there will be correspondence between their reconstructed manifold M_x and M_y (Deyle and Sugihara, 2011). Causal relationship can be unidirectional x(t) → y(t) or bidirectional x(t) ↔ y(t). To determine if there is causal relationship between the dynamics of ENSO and influenza, annual oscillatory components of influenza in five geographical areas were cross mapped with annual oscillatory components of ENSO using convergent cross mapping algorithm (Sugihara et al., 2012).

Matlab scripts, which contained relevant algorithms, were used to test for non-linearities in time series (Barnett, 2007), for noise assisted multivariate empirical mode decomposition (Rehman and Mandic, 2010b), for time-dependent intrinsic correlation (Chen and Huang, 2010), to generate surrogate time series (Leontitsis, 2004), and for multifractal detrended fluctuation analysis (Ihlen, 2012). All Matlab scripts were implemented in Matlab-Octave programming language.

Lyapunov spectrum, lag, embedding dimension, Takens, and phase space plots were implemented in the non-linearTseries package of R Statistical Programming and Environment, Austria, version 3.2.2, 2015 (R Core Team, 2016). Embedding dimensions m were calculated for IMF2s using Cao's algorithm (Cao, 1997), and lags τ were estimated using mutual information algorithm (Martinerie et al., 1992), both implemented in non-linearTseries package of R Statistical Programming and Environment, Austria, version 3.2.2, 2015 (R Core Team, 2016). Python programming language pyunicorn package was used to generate matrices of recurrence plots and networks, while its astropy package was used for bootstrap statistics. Python seaborn and matplotlib packages were used to produce distributions, time series, and angular histogram plots. Calculations of CPR index to determine synchronization of dynamics were performed as described in publication (Romano et al., 2005).

Regions of northern and southern hemispheres in the study are shown in Figure 1. Monthly time series of influenza in the five regions show consistent seasonal peaks, but multiple peaks occurred during the 2009–2010 pandemic in northern hemisphere regions (Figures 2A–E). Seasonal epidemics occurred from October to March in North America, October to April in Nordic/Europe, December to March in Asia, which range from middle of fall to middle of spring of northern hemisphere, but from May to November in South America, and June to October in Oceania, which range from late fall to middle of spring of southern hemisphere (Figures 3A–D). Seasonal epidemics peaked in January in North America and Asia, February in Nordic/Europe, July in South America, and August in Oceania, during the winter of both hemispheres (Figures 3A–D). The Rayleigh's test for non-uniformity of circular distribution was statistically significant at p < 0.0001 for all the regions. Of five El Niños which occurred in 2002–2003, 2004–2005, 2006–2007, 2009–2010, and 2015–2016, onsets of three were in summer, one in spring, and one in fall. The peaks of three were in winter (December), and two in fall (October and November).

The time series of ENSO and influenza of all regions show statistically significant non-linearities, p < 0.001. Synchrosqueezing transform of ENSO shows that time-varying oscillatory components with periodicities from less than a year to multi-decadal are present in the spectra, (Figures 4A–F). Influenza time series in North America, Nordic/Europe, Asia, South America, and Oceania also shows spectra of time-varying oscillations with periodicities from less than a year to multidecadal. Time-varying spectral of influenza time series showed concentration of power predominantly at 1.0 year periodicity in all regions except Asia, which has power concentrated at multiple periodicities.

Figures 5A–I shows intrinsic mode functions 2–4 for Asia, Oceania, and South America regions. Synchrosqueezing transform of intrinsic mode functions showed that the second intrinsic mode function (IMF2) is the annual oscillatory component. TDIC of instantaneous frequencies of ENSO and influenza in five regions show correlations > 0.9 for all windows, (Figures 6A–D). Similarly, time dependent intrinsic correlation of instantaneous amplitudes of ENSO and influenza show correlations > 0.9 for all windows, (Figures 7A–D).

Phase spaces of all oscillatory components of ENSO and influenza time series show dissipative dynamics, (Figures 8A–I). The trajectories of oscillatory components of ENSO and influenza graphically appear to run in the same manifold. The maximum Lyapunov exponents of annual oscillatory components of ENSO, and of influenza in all five regions are positive. Hurst spectra of multifractal detrended fluctuation analysis show that annual oscillatory components of ENSO and influenza in five regions of the study are multifractals.

Short diagonal lines in the recurrence plots of ENSO is typical of chaotic dynamics (Figure 9E). Recurrence plots of influenza in the five regions also show line structures typical of chaotic dynamics, (Figures 9A–D). The line structure of both dynamics, however, suggest presence of quasi-periodic and chaotic regimes. Peaks and troughs of divergence and determinism confirm the presence of regimes and transitions in both dynamics (Figures 9F–J). The bootstrap distributions of determinism and divergence for surrogate time series, influenza, and ENSO show that determinism and divergence of ENSO and influenza are highly statistically significantly different from surrogate data, p < 0.0001.

Joint recurrence plot of ENSO and influenza show line structures similar to the recurrence plots of ENSO and influenza, (Figures 10A–E). Global clustering coefficient, and average path length confirm the presence of quasi-periodic and chaotic regimes (Figures 11F–J). Bootstrap statistics shows statistically significant difference in global clustering coefficient and average path length of quasi-periodic and chaotic regimes, p < 0.0001.

CPR index peaked during the quasi-periodic regime in all five regions, but reached lowest values before 2012 in northern and southern America, but after 2012 in Europe, Asia, and Oceania, (Figures 11A–E). There is highly statistical significance difference of CPR index of ENSO and influenza compared with 5,000 samples of their surrogate time series (p < 0.0001), (Figures 11F–J).

The signature of annual oscillatory component of ENSO on influenza in North America is statistically significant, but not vice versa (Figure 12A). Convergence, which is improvement of cross-mapped estimates as the length of time series increases is present. Statistically significant signature of annual oscillatory component of ENSO in influenza is also present in South America, Europe, Asia, and Oceania (Figures 12B–E).

The five highest peaks of influenza epidemics occurred during chaotic regimes in all regions of the study. Of the five highest peaks two occurred during the 2009–2010 influenza pandemic in North America, South America, and Oceania, but one in Europe and four in Asia, (Figures 10A–E). There is highly statistical significant difference between the bootstrap distributions of 5,000 bootstrap samples of average path length and global clustering of windows, which contained the five highest influenza peaks compared with the remaining windows, (p < 0.0001).

The power spectra of ENSO of the past 21,000 years, which was derived from paleoclimate data, showed presence of oscillatory components from sub-annual to 64 years periodicities (Liu et al., 2014). The synchrosqeezing spectra show that ENSO is composed of time-varying sub-annual, annual, decadal, and multidecadal oscillations (Figures 4A–F). The periodicities of the oscillatory components are consistent with the multiple timescales that ENSO modulates the climate. Influenza time series of northern and southern hemisphere regions, which also show sub-annual, annual, decadal, and multi-decadal oscillations, indicate that influenza time series have oscillatory components of similar periodicities to ENSO (Figures 4A–F). Although the power of ENSO spectra is strong in multiple bandwidths from sub-annual to multi-decadal periodicities, that of influenza is strong in the bandwidth about 1 year periodicity except for Asia region, which has strong sub-annual to multi-decadal oscillatory components (Figures 4A–F). Thus, both ENSO and influenza are multicomponent systems with oscillations of similar periodicities.

Annual warming of equatorial east Pacific ocean during Christmas was documented in a lecture in 1895 (Cane, 1986), but references to this phenomenon, which lasts about 2 months, were made by Peruvian fishermen in the 18th century (Julian and Chervin, 1978). The annual oscillatory cycle of ENSO, which is tightly coupled to season cycles, has strong impact on global precipitation (Kessler et al., 1998; Stuecker et al., 2015). Very strong time dependent intrinsic correlations of instantaneous frequencies and amplitudes of annual oscillatory cycles of ENSO and influenza indicate that frequencies and amplitudes of both oscillations are coupled (Figures 6A–D, 7A–D). The ENSO, which is amplitude and frequency modulated (An and Wang, 2000; An and Jin, 2011), had periodicity of 3–4 years for El Niño events between 1872 and 1910, 5–7 years between 1911 and 1960, and 5 years between 1970 and 1972 (Torrence and Compo, 1998; Torrence and Webster, 1999). Thus, strong coupling of amplitude and frequency modulations of irregular ENSO with influenza indicate functional relationship.

Models of ENSO dynamics which predicted chaos include discharge-recharge oscillator (Jin, 1997), delayed oscillator (Battisti and Hirst, 1989), and linear annual subharmonic steps to chaos (Jin et al., 1994). In this study, however, ENSO dynamics was reconstructed from phase space embedding, which is derived from Taken's theorem (Takens, 1981), (Figures 8A–I). The phase space, state or vector space, of deterministic dynamic systems are usually defined by a set of first-order differential equations, which ensures uniqueness of the trajectories (Kantz and Schreiber, 2003). Every point in the phase space, which recovers the topology and preserves the properties or behavior of the dynamic system (Packard et al., 1980), specifies the state of the system (Kantz and Schreiber, 2003). Differential equations are, however, not available for ENSO and influenza, which are time series of scalar values. While the phase space volume of conservative dynamics system is preserved as time evolves, that of dissipative dynamic system contracts as time evolves t → ∞ (Grebogi et al., 1987). Figures 8A–I, which show that the phase space volumes of all oscillatory components of ENSO and influenza contract, indicate dissipative dynamics of both systems.

Trajectories of dissipative systems return to a region of phase space called attractor as time t → ∞. The geometry of phase space attractors characterizes dynamic systems. Phase space attractors with Euclidean geometry like fixed point, limit cycle, or torus are simple, while phase space attractors with fractal geometry are chaotic (Ruelle and Takens, 1971). Fractals are geometrical objects that have non-integer dimensions, which exceed their topological dimensions. This implies they are not locally linear when dilated or contracted, nor differentiable at all scales (Mandelbroit, 1989). The Lyapunov exponent, which measures the sensitivity of dynamic systems to small changes in initial conditions or the average rate of divergence of close trajectories (Wolf et al., 1985), is an indicator of the presence of chaos. A dynamic system is chaotic if the maximum Lyapunov exponent (λ) is > 0. Since the maximum Lyapunov exponents of annual oscillatory components of ENSO and influenza are positive, the attractors of both dynamics have fractal dimensions, which indicate chaos. Thus, the dynamics of annual oscillatory components of ENSO and influenza are chaotic as predicted by theoretical models (Vallis, 1986; Thornley and France, 2012).

Correlation dimension of 3.5 for sea surface temperature (SST) of east Pacific ocean (Tziperman et al., 1994) showed that ENSO has fractal geometry. Scale invariance of multivariate El Niño Index (MEI) between 1950 and 2008 (Mazzarella and Giuliacci, 2009), and detrended fluctuation analysis of southern oscillation index (SOI) between 1866 and 2000 also showed that ENSO is fractal (Ausloos and Ivanova, 2001). Similarity of ENSO waveforms during five influenza pandemics from 1876 to 2015, which suggests scale invariance (Oluwole, 2016), is indicative of fractal time series (van Rooij et al., 2013). While chaotic trajectories have fractal geometry, not all fractals are self-similar under magnification (Mandelbroit, 1989). Fractal time series, which do not have the same reduction ratio of shape at all scales (Mandelbroit, 1989), are multifractals. Multifractality of the annual oscillatory components of ENSO and influenza of five regions of the study indicate statistical self-similarity, which is typical of natural fractals.

Recurrence, which describes the return of trajectories to regions of phase space that have been previously visited, is characteristic of dynamic systems (Kantz and Schreiber, 2003). Poincaré recurrence theorem states that trajectories of dynamic systems return to states that are close to the initial state as time evolves sufficiently long (Kantz and Schreiber, 2003). The recurrence plot is a graphical tool to display the pattern of recurrence of dynamic systems (Eckmann et al., 1987). Its line structures reveal the dynamics of trajectories in phase space (Eckmann et al., 1987). Periodic dynamics show long uninterrupted diagonal lines in recurrence plots, while chaotic dynamics like the Lorenz system show short diagonal lines (Gao and Cai, 2000). There are no uninterrupted diagonal lines in the recurrence plot of annual oscillatory component of ENSO, but periods of relatively long and short diagonal lines which resemble recurrence plots of theoretical chaotic dynamics like Lorenz and Rössler systems (Figure 9E). The recurrence plots of annual oscillatory components of influenza in five regions of the study have line structures similar to that of ENSO (Figures 9A–D). Thus, recurrence of annual oscillatory components of ENSO and influenza show chaotic dynamics, but regimes are not homogeneous.

Chaotic regimes of dynamic system can be preceded by period doubling, intermittency, crisis, and Ruelle-Takens, which describe changes of attractors from simple to chaotic (Argyris et al., 1993). The dynamics of climate data has been shown to have periodic, quasi-periodic, and chaotic regimes (Donges et al., 2011). Periodic or quasi-periodic regimes that have long diagonal lines, which are described as unstable periodic orbits, are present in chaotic systems like the Lorenz and Rössler system (Marwan et al., 2007). Quasi-periodic regimes preceded chaos in theoretical models of ENSO which included season cycles (Tziperman et al., 1994). Determinism (DET), which is the ratio of recurrence points that form diagonal lines of at least l_min and all recurrence plots, and divergence (DIV), which is the inverse of the longest diagonal line L_max (Marwan, 2011) are metrics of diagonal structures, which show maximal values at periodic-chaos and chaos-periodic transitions (Marwan et al., 2007; Marwan, 2011). Peaks of divergence and determinism confirm the presence of regimes and transitions in the dynamics of ENSO and influenza (Figures 9F–J). The highly statistically significant difference between the distributions of determinism and divergence of ENSO and influenza and their surrogate time series indicate these observations are not random. Thus, the annual oscillatory components of ENSO and influenza follow quasi-periodic route to chaos.

Joint recurrence plot assesses the probability that similar points in phase space are visited by two chaotic dynamics, which are from physically different time series (Marwan et al., 2007). Figures 10A–E, which shows that joint recurrence plots of ENSO and influenza in five regions have similar line structures as recurrence plots of ENSO and influenza, indicate that the phase spaces of both systems are similar. Global clustering and average path length, which are metrics of network analysis, have been used to detect transitions in chaotic dynamics (Donner et al., 2011; Zou et al., 2012). Global clustering coefficient (C) is the mean of the local clustering coefficient of all vertices, while average path length (L) is the mean geodesic distance between all pairs of vertices where paths exist (Donner et al., 2011; Zou et al., 2012). Global clustering coefficient is high during unstable periodic orbits, while the average path length is high during chaotic regimes. These metrics indicate that quasi-periodic and chaotic transitions are also present in joint recurrence of ENSO and influenza (Figures 10F–J). The highly statistically significant differences between the distributions of global clustering and average path length metrics of ENSO and influenza compared with surrogate time series indicate that the observations are not random. Thus, the trajectories of annual oscillatory components of ENSO and influenza run in similar geometrical phase space.

Severity of seasonal influenza epidemics vary interannually from mild to severe (Doshi, 2008; Frieden et al., 2010) even when the circulating influenza virus strains are unchanged. Seasonal influenza epidemics in the USA was moderately severe in 2003–2004 and 2007–2008, moderate in 2002–2003, but mild in 2011–2012 (Bresee and Hayden, 2013). The virulent H3N2 influenza virus, which circulated during the moderately severe seasons of 2007–2008 and the mild season of 2011–2012 (Bresee and Hayden, 2013), show that virulence is not the only determinant of severity. Five most severe peaks of influenza epidemics which map to chaotic regime of ENSO and influenza joint recurrence indicate that the most severe influenza epidemics occur during chaos (Figures 10A–E). Thus, chaotic regime has greater impact on seasonal influenza epidemics than quasi-periodic regimes.

Synchronization occurs when oscillating systems adjust to a common rhythm due to coupling (Pikovsky et al., 2001). Although two chaotic oscillators with similar initial conditions diverge exponentially and become uncorrelated, chaotic dynamics like Lorenz and Rössler systems have been shown to adjust to common rhythm. Synchronization of chaotic oscillators occurs when a property of their motion, like phase or amplitude, adjusts to common behavior due to coupling or forcing such that their trajectories evolve in a common phase (Boccaletti et al., 2002). The cross correlation index (CPR), which compares the probability P^(ϵ)(τ) that the trajectory returns to ϵ-neighborhood of a previous point on the trajectory, shows that the annual oscillatory components of ENSO and influenza are phase coherent during quasi-periodic and chaotic regimes in five regions in northern and southern hemispheres (Figures 11F–J). The highly statistically significant differences between the distributions CPR index compared with null values indicate that the observations are not random. Thus, ENSO and influenza dynamics are phase coherent in northern and southern hemispheres.

Synchronization is unidirectional, when one oscillator drives the other, but bidirectional when both oscillators adjust to a common manifold (Pikovsky et al., 2001). The concept of temporal precedence of time series was developed by Wiener, who proposed in 1956 that if prediction of a time series X by its past values is improved by the past values of another time series Y, then time series Y has causal influence on time series X (Wiener, 1956). Granger's formulation of this concept as autoregression models led to the Wiener-Granger causality (Granger, 1969). The trajectories of causally related dynamic systems share the same attractor manifold (Deyle and Sugihara, 2011). If time series x(t) causally influences time series y(t), the signature of x(t) exists in y(t) (Sugihara et al., 2012). Historical values of y(t) can, therefore, be used to estimate x(t). Convergent cross mapping compares the coordinates of reconstructed manifolds M_x and M_y for correspondence to determine if signature of x(t) exists in y(t) (Sugihara et al., 2012). If the average correlation between x(t) and its estimated values x(t) is high, it is concluded that there is sufficient information of x(t) in y(t) (Sugihara et al., 2012). Wiener-Granger causality, which is applicable to dynamic systems that are separable and have non-zero entropy rate is, therefore, not applicable to coupled non-linear dynamic systems, which contain information about each other and are not separable. Comparison of the reconstructed manifolds of annual oscillatory components of ENSO and influenza in North America region, which show statistically significant signature of ENSO in influenza but not vice versa, indicate unidirectional causal effect of ENSO on influenza (Figure 12A). Convergence, a necessary condition for causation (Sugihara et al., 2012), is fulfilled by the improved cross-mapped estimates as the length of time series increased (Figure 12A). Statistically significant signature of annual oscillatory component of ENSO in influenza, which is also present in South America, Europe, Asia, and Oceania regions, indicate that ENSO has unidirectional dynamic causal relationship with influenza in disparate regions of northern and southern hemispheres (Figures 12B–E). Thus, seasonality of influenza epidemics is driven by the dynamics of ENSO, although biological causation is due to viruses, while impaired immunity and overcrowding, which are contributory causes, are also induced by climate.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.