Long-Term Study of Heart Rate Variability Responses to Changes in the Solar and Geomagnetic Environment

This long-term study examined relationships between solar and magnetic factors and the time course and lags of autonomic nervous system (ANS) responses to changes in solar and geomagnetic activity. Heart rate variability (HRV) was recorded for 72 consecutive hours each week over a five-month period in 16 participants in order to examine ANS responses during normal background environmental periods. HRV measures were correlated with solar and geomagnetic variables using multivariate linear regression analysis with Bonferroni corrections for multiple comparisons after removing circadian influences from both datasets. Overall, the study confirms that daily ANS activity responds to changes in geomagnetic and solar activity during periods of normal undisturbed activity and it is initiated at different times after the changes in the various environmental factors and persist over varying time periods. Increase in solar wind intensity was correlated with increases in heart rate, which we interpret as a biological stress response. Increase in cosmic rays, solar radio flux, and Schumann resonance power was all associated with increased HRV and parasympathetic activity. The findings support the hypothesis that energetic environmental phenomena affect psychophysical processes that can affect people in different ways depending on their sensitivity, health status and capacity for self-regulation.

Standing wave field-line resonances are typically classified as Pc3 to Pc5 waves corresponding to frequencies between 1 mHz and 100 mHz. Pc1 and 2 oscillations are classified as traveling waves, which can have frequencies up to 5 hertz, and are typically excited by geomagnetic substorms 76 . Studies have found that increases in field-line resonances can affect the cardiovascular system, which may be because the Pc waves are in a comparable range with those of the autonomic nervous and cardiovascular systems 77 . Khabarova and Dimitrova demonstrated that the magnetic waves in the 2-10 mHz region had the higher correlations with increased blood pressure than geomagnetic indices 69 . In addition, Zenchenko et.al found there was a significant degree of synchronization between HRV rhythms and the osculation's in geomagnetic field in the frequency range between 0.5 to 3.0 mHz in two-thirds of their experiments of over periods of 4 to 30 minutes 78 . It has also been demonstrated that ANS activity not only reacts to shifts in geomagnetic and solar activity, it can also synchronize with rhythms in the time-varying magnetic fields related to the Schumann resonances (SR) and geomagnetic field-line resonances 15,79 .
Pobachenko et al. 82 examined the SR and the EEGs in a group of participants over six-weeks and found that during the daily cycle, changes in the EEG were similar to variations in the SR. The highest correlations between the SRs and brain rhythms occurred when the magnetic activity was increased. Persinger et al. have also studied SR and EEG activity the in real-time and have shown that many of the SR frequencies can be observed in the power spectrums of most human brain activity 83,84 . They have also shown that the spectral profiles within the EEG activity displayed recurrent transient segments of real-time coherence (synchronization) with the first three resonant frequencies of the SRs (7-8 Hz, [13][14][19][20]. This suggests that under certain conditions, variables affecting the Schumann parameters (such as solar wind) may affect brain activity, such as modifications of perception and dream-related memory consolidation 84 . Altered EEG rhythms in response to changing magnetic fields have also been observed by Belov et al., with low frequency magnetic oscillations (around 3 Hz) having a sedative effect 85 .
In the study reported here, the potential correlations between solar and magnetic factors and the time course or lags in autonomic nervous system responses to changes in solar and geomagnetic activity were examined by collecting HRV data for seventy-two consecutive hours each week over a five-month period using ambulatory HRV recorders. The intent was to examine the group's ANS responses, as reflected in HRV, to changes in solar and geomagnetic background activity over an extended period. The maximum Ap Index level for each day during the study period is summarized in Table 1 and provides an overview of the general level of geomagnetic activity during the study period (April 1 -August 31, 2012).

Methods and Procedures
Participants. Eighteen females, all healthy employees of the Prince Sultan Cardiac Center in Hofuf, Saudi Arabia (8 nursing staff, 6 housekeeping and, 4 from the research department) volunteered to participate in this study. The average age was 32 ± 8 years, ranging from 24-49 years. Inclusion criteria were no known physical or mental health disorders and were fulltime employees. Exclusion criterion was a known health disorder, or taking any medications known to affect autonomic function. All participants signed informed consent and were free to withdraw from the study at any time. Two participants experienced uncomfortable irritation at the ECG electrode sites and dropped out of the study. The research was performed in accordance with all relevant guidelines and met all applicable regulations for the ethics of experimentation and was approved by the Ministry of Health Research Centers ethical decision party (R.C. 100/25). Data Collection. All participants underwent weekly 24-72 hour ambulatory HRV recordings with Bodyguard HRV recorders (Firstbeat Technologies Ltd, Finland). The recorder samples the ECG at 1000 Hz and calculates the inter-beat-interval (IBI), which is the time in milliseconds between consecutive heartbeats. IBI data was stored locally in the device memory and uploaded loaded to an FTP site at the end of each week. Participant recordings were generally 72-hours in length and were scheduled once a week over a 5-month period between April and the end of August 2012. A total of 960 twenty-four hour long HRV recordings were obtained.  Table 2 summarizes frequency band ranges for the HRV and magnetic field measures used in this study. The HeartMath Institute maintains a network of highly sensitive induction coil magnetometers (Zonge ANT-4; sensitivity 10 −12 T) as part of a special project called the Global Coherence Initiative 87 . Each site includes two magnetometers positioned in the north-south and east-west axis to detect local time varying magnetic field strengths over a relatively wide frequency range (0.001-50 Hz) while maintaining a flat frequency response. The data acquisition infrastructure collects and timestamps all data using GPS time signals before uploading to a common server. Each magnetometer is continuously sampled at a rate of 130 Hz. Figure 1 shows the time domain data for the environmental data activity across the study period.
Statistical Analysis. Regression Analysis. These data were analyzed by means of multivariate linear regression. Each HRV variable is considered separately as a single dependent variable. In the first analysis stage, each participant's dataset is handled separately. For each HRV variable and each participant dataset, regression coefficients were computed to minimize the total squared error of a linear model. The assumed model is that where HRV is the dependent HRV variable, E1, E2, etc., are the independent environment variables, and R1, R2, etc., are the regression coefficients. Before the model results were calculated, hourly observations, which have missing data in any variable, were discarded so that only observations with a complete set of environmental measures were used. This so-called "complete cases" approach is a standard method for dealing with missing data in regressions. In addition, the data were normalized so as to eliminate circadian confounds. A side effect of the circadian normalization is that it forces both dependent and independent variables to have an overall mean of zero over the observation period, so that the regression model can safely ignore "intercept" terms and compute only "slopes", as implicit in the model described above. In addition to the considerations above, the analysis also evaluated possible time-offset delayed effects, since it is known that HRV measures may show a reaction to a stimulus several hours after the stimulus is applied. The hourly timestamps in both HRV and environment data allow each set of HRV measurements to be linked unambiguously to a set of environment measurements. These can also be used to examine behavior with a lag in response time simply by offsetting the timestamps that are to be matched. Temporal offsets ranging from zero to 40-hours were used, meaning that the calculation described above was actually repeated 41 times (counting the 0-offset model for matched times). Figure 2 illustrates how time-offsets are used to examine lags in response times.
The result of all of these calculations is a multidimensional array of regression coefficients with their associated standard errors, computed by the R functions "lm" and "summary.lm". There is one regression coefficient (and computed statistical standard error) for each of six HRV variables, as influenced by each of nine environment variables, for each of 16 participants at each of 41 time lags. The next analysis step after computing the models was to compute a cross-participants' average for the regression coefficient, thereby eliminating the "participants" dimension of the multidimensional array. This average is the mean of the participants' regression coefficients, individually weighted by their associated standard errors in the usual least-squares formula: where M is the weighted mean, x i are the individual observations, and σ i are their associated measurement uncertainties. (This weighting is "standard" because it minimizes the total squared error between the mean and the individual observations). The result of this calculation is that the array of individual participant observations can be replaced by a single value, the weighted mean, with an associated measurement uncertainty. Environmental data activity across the study period. There was a large increase in the Kp and Ap indexes that occurred on July 14 th , which resulted from a coronal mass ejection that hit the earth's magnetic field at approximately 1800 UT that day.
In addition to these two values, the statistical evaluation of the participant-averaged data also computed the Z-score M the p-value of each Z-score, a chi-squared value expressing the amount of excess variation between participants beyond that expected from the regression uncertainties, and the p-value of this chi-squared. These last few values relate only to inter-operator variability, which is beyond the scope of the current analysis. The Z-score, however, was used extensively in both analysis and visualization, since it is directly related to the statistical significance of a given result and provides a scale-independent measure for comparing the contributions of different environment variables. Z-scores are sometimes referred to as standard normal deviates. Strictly speaking, this assumes that the errors are normally distributed, but since the calculation discussed here is the mean of sixteen independent models computed for the sixteen participants, any departures from normality will be very small regardless of the underlying distribution of observational errors for an individual participant. Multiple Analysis Corrections: After averaging across the contributions of the participants, one is still left with a multidimensional array of results for nine environment variables, by six HRV measures, by 41 time lags, for a total of 2214 analysis results. As insurance against introducing analysis errors by selecting significant results that appear purely by chance, we applied a Bonferroni correction for multiple analyses to this array of analysis results, or to any subset that we may be examining in more detail. The usual Bonferroni correction is made by multiplying the p-value of the most significant analysis by the number of analyses. Strictly speaking, the proper correction for N analyses, with the strongest result having a p-value of P, is This exact formula expands, however, into The simple formula corr can thus be seen to be conservative, always larger than the actual p-value, with an error that is negligible for P ≪ 1. For the record, the strongest individual result in this full analysis array has a p-value of 1.4 × 10 −32 , which becomes 3.1 × 10 −29 after Bonferroni correction for 2214 analyses.
Lag Analysis. As noted above, the multiple-analyses corrections must include the consideration that 41 different time lags were examined by the same regression methods. When allowing for the fact that the actual relationship between geomagnetic environmental variables and an HRV measure might involve an unknown time lag, the most straightforward way of dealing with the multiple analysis problem is to apply a 41-test Bonferroni correction to the lag results for a single HRV-to-environment regression coefficient. This allows one to draw a valid conclusion for whether that coefficient is displaying a statistically significant value at any of the time lags examined. In graphical visualizations of the lag results, this is often done implicitly by plotting the unmodified series of Z-scores and placing a 5% significance envelope at ±3.24, the two-tailed 5% threshold for the largest of 41 standard normal deviates (1% significance at ±3.67). The analysis has 54 such coefficients, analyzing six HRV measures by nine independent variables. Some evaluations, such as an overall significance for the entire model, require that a factor of 54 Bonferroni correction be applied to the most significant individual results. Other types of analysis, such as evaluating how many of the regression coefficients are individually significant at some time lag, require no further Bonferroni correction beyond the one applied for the multiplicity of lags.

Results
Environmental and HRV measure correlation results are provided in Tables 3 and 4. In agreement with previous studies, solar wind speed was highly correlated with Kp (r 0.50, p < 0.01) and Ap (r 0.35, p < 0.01) indexes, ULF power and negatively correlated with cosmic ray counts (r −0.15, p < 0.01) [88][89][90] . As expected, the solar radio flux (F10.7) was also highly correlated with the number of sunspots (r 0.81, p < 0.01). The Schumann Resonance Power was negatively and highly correlated with cosmic ray counts (r −0.58, p < 0.01). ULF power was positively correlated with solar wind speed (r 0.44, p < 0.01), Kp (r 0.58, p < 0.01), Ap (r 0.61, p < 0.01) and PC(N) (r 0.43, p < 0.01) indexes. The correlations among the HRV variables were as expected and in agreement with other studies with the exception of HF power, where correlations to IBI, Total Power, and VLF power were all higher than seen in individual 24-hour recordings (Table 4) 91 .
Results of the multivariate linear regression between the environmental and HRV variables revealed a number of significant findings with different HRV metrics responding at different time lags with changes in environmental factors. Z-scores for regression coefficients are shown in Figs 3, 4 and 5. Detailed Z-score tables were too large for inclusion here and can be found in the supplemental information. There were significant autonomic nervous system responses reflected in the HRV to changes in cosmic ray counts (Fig. 3a), which were strong and consistent. The TP (Z 7.30 to Z 10.49, p < 0.01), VLF (Z 5.20 to Z 8.19, p < 0.01), LF (Z 8.49 to Z 11.88, p < 0.01) and HF (Z 6.63 to Z 10.16, p < 0.01), all quickly and strongly responded, and were positively correlated across the entire 40-hour period. IBIs responded from hour 4 through 12, and again at hours 36, 37 and 40 (Z 3.32, p < 0.01 to Z 4.11, p < 0.01). The negatively correlated LF/HF ratio was significant only at hour 38 (Z −3.33, p < 0.05).
For changes in the number of sunspots (Fig. 3b), IBIs were positively correlated in hour 12 through 15 (Z 3.42, p < 0.05 to Z 3.84, p < 0.01). TP was significant at hour 10, and again during hours 12 to 15 (Z 3.25, p < 0.05 to Z 3.73, p < 0.01). HF was also positively correlated from hour 10 to hour 15 (Z 3.28, p < 0.05 to Z 4.13, p < 0.01). The VLF was nearly the same, but was not significant until the 13 th hour and remained significant through hour 15 (Z 3.43 to Z 3.52, p < 0.05). The LF and LF/HF ratio were non-significant.
For solar radio flux (F10.7) (Fig. 3c), there was an immediate, robust positively correlated response in the IBIs, TP, VLF and HF power. The IBI response was significant for the first two hours (Z 3.73 to Z 3.96, p < 0.01). The TP power response remained significant over the first 8-hours, and again during the 19 th and 20 th hour and from hour 34 to the end of the analysis period (Z 3.25, p < 0.05 to Z 5.26, p < 0.01). The VLF response was significant     Figure 4a shows the correlations for Schumann Resonance Power (SRP) and Fig. 4b ULF Power. There were a number of significant correlations. For the SRP, the IBIs were the first to respond, and became significant in the 4 th and remained significant through hour 31 and again at the 36 th to 40 th hour (Z 3.37, p < 0.05 to Z 6.76, p < 0.01). The Total Power response was significant during hours 8 to 18 and from hours 21 to 37 (Z 3.30, p < 0.05 to Z 6.23, p < 0.01). VLF power was significant from hours 10 to 17 and again at hours 22 through 33 (Z 3.46, p < 0.05 to Z 5.14, p < 0.01). LF power became significant starting at hour 9 and continuing through hour 16 and again between hours 21 and 38 (Z 3.45, p < 0.05 to Z 5.74, p < 0.01). In addition, HF power was positively correlated between hours 9 to 15, hours 23 to 30, and hours 37 to 39 (Z 3.28, p < 0.05 to Z 5.76, p < 0.01). The LF/HF ratio was correlated only during the first hour (Z 3.86, p < 0.01).
In response to changes in ULF power, there was an immediate positively correlated response lasting 3 hours in the IBIs (Z 4.42 to Z 5.26, p < 0.01), and HF power (Z 3.51, p < 0.05 to Z 4.54, p < 0.01). HF power became negatively correlated during the last two hours (Z −3.49, p < 0.05 to Z −3.69, p < 0.01). The LF/HF ratio was negatively correlated during the first three hours (Z −3.72 to -Z 5.02, p < 0.01).
Keeping in mind the relatively few magnetic disturbances that occurred during the study period, there were very few correlations with the indices of geomagnetic disturbances. The Polar Cap activity (PC.N) (Fig. 5b) had a significant negative correlation with TP in the 7 th hour (Z-3.42, p < 0.05), and IBI in the 8 th hour (Z −339, p < 0.05). There were no significant correlations to the Kp Index. The Ap index (Fig. 5a), was positively correlated with TP, VLF, and LF in the 9 th hour (Z 3.42, Z 3.37 and Z 3.30, p < 0.05) respectively, and LF was also significant in the 6 th hour (Z 3.37, p < 0.05). The LF/HF ratio was negatively correlated during the 38 th hour (Z −3.29, p < 0.05).

Discussion
Overall, this study confirms that daily ANS activity, as reflected by HRV measures, reacts to changes in geomagnetic and solar and activity during periods of normal undisturbed activity. Furthermore, these ANS responses  are initiated at different times after the change in the various environmental factors and persist over different lengths of time. It should be noted that a common finding is that different individuals respond differently to changes in the same environmental variable 69 . In a separate analysis of the data presented here, collaborators at the Lithuanian University of Health Sciences developed a new measure for assessing sensitivity to magnetic field variations 92 . Although different participants had different responses at the individual level, the analysis presented in this paper concerns responses at the group level.
It is clear that a major driver of changes and disturbances in Earth's magnetic field environment are the sun and solar wind 88,93 . Consistent with these findings, the solar wind speed was highly correlated with Kp, Ap and the PC(N) all of which reflect magnetic field disturbances. We also found that ULF power, which is related to magnetic field-line resonances, was positively correlated with solar wind speed, and indices of field disturbance was negatively correlated with cosmic ray counts which is consistent with the well-known inverse action of solar and geomagnetic activity and cosmic ray counts at the Earth's surface 90 .
Regarding HRV responses, IBIs have an inverted relationship to heart rate where larger IBIs equated to a lower heart rate. Heart rate and IBIs are an ideal indicator of changes in the relative balance between parasympathetic and sympathetic activity and how the autonomic system responds and adapts to various types of stressors or challenges 40 . If an environmental variable is negatively correlated with IBIs, it indicates that heart rate increases with increases in that variable, which suggests a physiological stress reaction occurred. On the other hand, a positive correlation with IBIs indicates a lower heart rate. There were robust positive correlations between IBIs and Schumann Resonances and to a lesser degree with cosmic rays.
The positive correlation found between HF power and solar radio flux indicates an enhancement of parasympathetic nervous system activity during periods increased solar radio flux. This was of particular interest because a previous study with 1,643 participants in 51 countries found that the solar radio flux index was positively correlated with reduced fatigue, improved positive affect, and mental clarity while increases in solar wind speed had the opposite effects 94 . The potential beneficial effects of the solar radio flux was also observed in several studies that looked at death rates from various causes which found a strong and inverse relationship between the F10.7 and death rates 73,95 . The solar radio flux may be an important mediator of the anticipatory reactions observed by Tchijevsky which can occur several days before increases in the solar wind reach Earth and create magnetic disturbances. Of course, other sources of radiation such as X-ray, cosmic rays and UV emissions from the sun during coronal mass ejections are also likely aspects of the anticipatory reaction.
The strong and positive correlations between HF power and cosmic rays, as well as TP and VLF power, suggest a favorable response to increased cosmic rays, at least in a healthy population. The ANS response to increases in cosmic rays was immediate and continued throughout the forty-hour analysis window. It is clear from Fig. 3 that the largest ANS relationships were to cosmic rays, which supports the perspective of Stoupel, that cosmic rays are emerging as a principal factor of the environmental forces that affect human physiology 73 . The observations in this study, which suggests a positive reaction to cosmic rays, may seem at odds with the finding that increased cosmic rays, in conjunction with low geomagnetic activity, is positively associated with the timing of sudden cardiac deaths 72,95 and strokes 73 . However, it may be that in only looking at death rates, which primarily reflect a sick and elderly population, that important aspects of how these environmental factors affect healthy populations could be overlooked. For example, in a large study of a younger population with potential inflammatory related problems who had serum C-reactive protein (CRP) tests taken, found a robust, inverse correlation between C-reactive protein levels and cosmic rays 96 . Interestingly, inflammation and higher levels of C-reactive protein have also been shown to be associated with lower levels of HRV, especially with a reduced VLF band 97,98 .
The other environmental variable that was strongly associated with increased the HF, LF and VLF power and TP of the HRV measures was Schumann resonance power (SRP). This was accompanied by the positive correlation to IBIs (lower heart rate), which was also significant during most of the analysis period. Also supporting a beneficial effect of enhanced SRP is a study that found reduced systolic, diastolic and mean arterial blood pressure during times of higher SRP 99 . Persinger and colleagues have conducted a number of studies showing that not only are the base rhythms of the brain similar, but real-time coherence between Schumann resonances and brainwaves can occur in participants globally, and that the intensity of the Schumann resonance is linearly related to the amount of coherence 83,100 . They also have proposed that information transfer can occur between human brains and the Schumann resonances 84 .
We have suggested that although the specific mechanisms are not yet clear, the energetic environmental factors discussed above either can directly or indirectly affect human psychophysiology and behaviors in different ways depending on the health status and maturity of the individuals 87 . This perspective is supported by the finding that increases in solar radio flux, cosmic rays and Schumann resonance power are all associated with increased HRV and parasympathetic activity. The ANS also appears to respond quickly to changes in cosmic rays, Schumann resonance power and the solar radio flux. These may well be some of the key drivers of Tchijevsky's Index of Mass Human Excitability that clearly tracks the solar cycle 6 . On a larger societal scale, an increase in these energetic factors are associated with increased social unrest 3 , motivation 101 and human flourishing 8 . The findings that both extremely high as well as extremely low values of geomagnetic activity are associated with the timing of increased death rates 72 also suggests that our energetic environment affects people's energy levels and that low activity or disturbances can act as triggers in sensitive and unhealthy populations and serves to motivate and facilitate human activity.
Limitations. Although this type of study is correlational due to the independent measures of interest, it is a limitation. Additional merging of empirical results from a larger number of studies using different populations, designs, geographic locations and sampling increments, adds support to our findings. Although it is a truism that correlation does not imply causation, relationships generally must exist to be consistent with theories about causality.
SCIeNTIfIC REpoRTS | (2018) 8:2663 | DOI:10.1038/s41598-018-20932-x Consideration of the logically possible causal relationships between correlated variables A and B allows only four alternatives. Either A and B are causally unrelated and their correlated variation is happenstance, or A causes B, or B causes A, or the correlated variations in A and B are jointly caused by some external factor. In the current analysis, the "happenstance" hypothesis is the one that is ruled out by the high significance level of the results. It is at least plausible that external electromagnetic effects contribute causally to changes in HRV, while it seems utterly absurd to posit the reverse causal relationship that HRV changes in a small population of individuals in one geographical region causes changes in solar and geomagnetic activity. The fourth case, that both phenomena are caused by yet another factor external to either, suffers from the lack of any theoretical candidate for such an external factor. As a result, the hypothesis that environmental electromagnetic effects cause changes in ANS dynamics appears to be the most plausible causal interpretation of the observations. Because no one study or form of evidence can be considered as definitive, support for causality can best be argued when various classes of evidence all converge on the same conclusion. We have therefore reviewed and documented multiple lines of evidence that have rigorously tested the hypothesis that daily ANS activity, as reflected by HRV, responds to changes in geomagnetic and solar activity at different times and persists over differencing periods of time. Moreover, the results presented are consistent with and extend previously published results. Unique to our findings is the observation that changes in the various ambient magnetic environment affect the human autonomic nervous system differently with specific temporal response patterns.

Conclusions
Overall, this study strongly confirms that daily ANS activity, as reflected by HRV measures, responds to changes in geomagnetic and solar activity primarily during periods undisturbed by solar activity. Furthermore, these ANS responses are initiated at different times after the change in the various environmental factors and continue over different lengths of time. Solar wind was negatively correlated with IBIs indicating that heart rate increases with increases in solar wind that suggests a physiological stress reaction occurred. It appears that increased cosmic rays, solar radio flux, and Schumann resonance power are all associated with increased HRV and increased parasympathetic activity, and the ANS responds quickly to changes in these environmental factors. These may well be some of the key drivers of Tchijevsky's Index of Mass Human Excitability that clearly tracks the solar cycle. These findings support the hypothesis that these energetic environmental factors act as energy sources that outplay in different ways depending on an individual's health status and maturity level and capacity of self-regulation.