地磁场变化与人体心律参数之间的同步效应：植物神经系统的潜在作用

全文:

Background

One of the key interdisciplinary challenges in contemporary fundamental science is to understand how biological systems—from the molecular to the organismal level—respond to low-intensity environmental influences, including those related to space weather. Numerous studies have shown that strong solar flares, arrivals of plasma clouds to Earth, Forbush decreases, and planetary geomagnetic storms are accompanied by a sharp increase in the incidence of cardiovascular events such as myocardial infarction, stroke, and sudden cardiac death [1–5].

It has also been established that not only extreme space weather events but even moderate increases in geomagnetic activity (GMA) exert significant physiological effects on the human body. In such cases, external influences do not necessarily manifest as acute events or mortality but are instead reflected in significant changes in average values of physiological parameters related to various systems of the body, including the endocrine [6], nervous [7], and cardiovascular systems [8, 9]. Notably, such physiological responses have been observed not only in patients with functional impairments of the said systems but also in healthy individuals, including young adults [10].

These findings support the notion that space weather phenomena represent a significant environmental factor that warrants in-depth investigation, both for advancing fundamental understanding of living systems–environment interactions and for developing practical measures to protect human health from their adverse effects.

A major challenge in this area is the systemic nature of the organism’s response to external influences. For instance, during geomagnetic storms, significant alterations occur across a wide range of physiological indicators: blood pressure (BP) [9] and heart rate (HR) [3] increase, and parameters reflecting vascular tone—such as pulse wave velocity and endothelial function [8]—as well as microcirculation [11] also change.

The autonomic nervous system (ANS) has been repeatedly shown—at both the population and individual levels—to respond to changes in GMA. The clinical presentation of myocardial infarction associated with magnetic storms often includes a marked decrease in heart rate variability (HRV) [1, 6, 12].

Observational studies show that on geomagnetically disturbed days, the standard deviation of normal RR intervals (SDNN) index of HRV decreases by approximately 23% compared to “quiet” days. Total heart rate spectral power also declines, primarily due to a reduction in the power of the low-frequency (LF) and very-low-frequency (VLF) components, which reflect sympathetic modulation and baroreflex function [1, 12, 13]. At the same time, many researchers emphasize that HRV responses to GMA are highly individual [14, 15].

Experimental data also support these observations. In laboratory studies, a significant decrease in HRV was recorded in rabbits during simulated magnetic storms, which the authors attribute to baroreflex involvement [16]. The effects of artificial magnetic fields on human HRV were shown to depend on the characteristics of the field: under different conditions, exposure led to either increased or decreased stress levels [17]. However, due to the complexity of feedback and regulatory mechanisms, it is difficult in both observational and laboratory settings to determine which changes result directly from magnetic field variations and which are secondary.

About two decades ago, several studies reported frequency entrainment between various biological rhythms—primarily HR and electroencephalographic (EEG) activity—and GMF fluctuations with similar frequencies. This was initially observed in the hertz range, including frequencies of the fundamental Schumann resonances (8–14 Hz) and Pc1 geomagnetic pulsations (0.5–2.0 Hz) [18, 19], and later confirmed in laboratory experiments [15, 20–23].

Subsequently, our group reported a similar effect in the millihertz range (3- to 40-minute oscillation periods). We found that the dynamics of resting HR in healthy individuals showed statistically significant associations with variations in the GMF vector [24]. We later demonstrated that the dominant oscillation periods present in wavelet spectra of both HR and synchronous GMF variations largely overlapped during each 1- to 2-hour observation period [25, 26]. This synchronization effect was observed both in repeated intra-individual measurements [24, 25] and in recordings from groups of healthy volunteers and individuals with hypertension [24, 27]. An essential condition for detecting this effect was a state of resting wakefulness (but not sleep), which minimized interference from other HR-modulating factors. Given that the detected oscillation periods varied from one experiment to another, we ruled out the possibility of a random coincidence between the intrinsic frequencies of two oscillatory processes—biological and geophysical. We proposed a working hypothesis that the effect represents frequency entrainment of a specific biological process to concurrent GMF variations with matching frequencies [25]. This phenomenon appears to be analogous in nature to the effects described in previous works [1, 15, 18–23], but in the millihertz range, which had not been previously studied. We termed it the biogeophysical synchronization effect for the investigated frequency range of 0.5–5.0 mHz. However, it may be applicable to a much broader frequency range, extending from microwaves [28] to cosmic rhythms with periods spanning decades [29].

Aim: Given the extensive evidence implicating ANS involvement in the organism’s response to geomagnetic storms, we hypothesized that the mechanisms regulating autonomic balance may also contribute to the biogeophysical synchronization effect as one of the intermediate stages in the body’s response to such perturbations [24]. To test this hypothesis, we conducted a large-scale comparative analysis of the frequency of the biogeophysical synchronization effect in the dynamics of HR time series and in statistical parameters of HRV.

Methods

Experimental Data

This study was designed as an observational time series investigation and included 673 long-term recordings of cardiointervalograms (CIGs), each lasting from 60 to 120 minutes. Data were collected from two groups of volunteers who met the following inclusion criteria: age between 20 and 55 years; assignment to health groups I (excellent to good) or II (fair) based on preventive medical examination results; and willingness to participate in long and/or repeated cardiointervalographic recordings. Exclusion criteria comprised signs of hypertension or its complications, cardiac arrhythmias, and pulmonary diseases. Participants were not taking drugs affecting the cardiorespiratory system and refrained from intense physical activity (such as gym workouts). Group 1 consisted of 8 volunteers who underwent repeated recordings (at least 10 per person, total: 622 sessions). Group 2 included 39 volunteers who participated in 1 to 3 sessions each (total: 51 recordings). Comparing these two study designs—longitudinal (multiple recordings per subject in group 1) and cross-sectional (single measurements across individuals in group 2)—is essential, as heliobiological responses may be influenced both by inter-individual variability and by environmental conditions at the time of measurement.

Recordings were conducted from 2012 to 2024 in the Moscow, Leningrad, and Arkhangelsk regions of Russia. Summary information on group 1 participants is presented in Table 1. Group 2 included 39 individuals (14 men and 25 women) with a mean age of 38 ± 15 years. In post-hoc quality control, time segments associated with temporary health deviations—such as marked fatigue, acute respiratory infections, or psychoemotional stress—were excluded based on participants’ self-monitoring diaries. Sessions were also excluded if the participant had consumed coffee within 4 hours prior to CIG recording.

CIGs were derived from ECG recordings using standard lead I, with the subject in a supine position and a state of resting wakefulness following a 10-minute adaptation period.

 

Table 1. List of group 1 volunteers, anamnestic data, and median values of measured parameters Values are presented as Me (1^st quartile; 3rd quartile)

Volunteer number	Sex	Age, years	n	HR	RMSSD	SDNN	AMo	SI
V1	F	59	333	69.3 (65.2; 73.0)	21.9 (16.5; 27.9)	26.7 (20.9; 32.4)	59.8 (53.88; 66.0)	334.9 (239.6; 485.8)
V2	F	45	165	61.6 (59.6; 63.6)	36.9 (29.4; 45.9)	31.7 (27.7; 35.5)	52.6 (47.3; 57.2)	197.1 (157.3; 234.8)
V3	F	30	64	63.8 (60.2; 69.2)	45.3 (35.8; 52.9)	48.4 (41.3; 52.3)	38.7 (36.7; 44.0)	100.1 (83.4; 141.8)
V4	M	37	19	78.7 (75.6; 79.7)	18.8 (16.4; 21.6)	39.0 (36.0; 41.6)	46.6 (44.5; 48.3)	215.0 (185.5; 232.6)
V5	F	53	10	80.1 (74.3; 80.3)	19.5 (17.5; 23.4)	28.0 (26.7; 30.4)	55.1 (51.8; 57.6)	308.8 (260.4; 358.7)
V6	M	59	10	62.5 (60.5; 63.1)	18.2 (16.5; 21.6)	25.0 (23.0; 29.0)	60.9 (57.9; 63.6)	309.6 (249.0; 380.7)
V7	F	42	11	71.9 (66.5; 73.5)	42.1 (36.2; 46.4)	45.8 (42.9; 49.9)	42.3 (38.5; 43.1)	128.6 (103.2; 146.6)
V8	F	27	10	77.3 (73.5; 79.0)	40.6 (34.4; 47.4)	58.8 (53.8; 65.9)	33.8 (32.1; 35.8)	88.7 (73.9; 106.9)

Note: F. female, HR, heart rate; M, male; RMSSD, root mean square of successive differences; SDNN, standard deviation of normal NN intervals; AMo, mode amplitude; SI, stress index.

 

Ethical Approval

The study was conducted in compliance with all principles of ethics and humanity (WMA Declaration of Helsinki, 2013) and posed no risk to participants. It was approved by the Bioethics Committee of the Institute of Theoretical and Experimental Biophysics, Russian Academy of Sciences (Protocol No. 06/2012, dated June 1, 2012). A written informed consent was obtained from all participants.

Based on the CIG recordings, time series (each lasting 60–120 points, i.e., minutes) of minute-by-minute values for the following HRV parameters were generated [30]:

1) Heart rate (HR) in beats per minute;

2) Root mean square of successive differences between adjacent NN intervals (RMSSD) in ms; reflects vagal modulation of heart rhythm;

3) Standard deviation of normal RR intervals (SDNN) in ms; reflects total HRV and vagal tone in short-term recordings;

4) Amplitude of mode (AMo) in %, which is the percentage of RR intervals corresponding to the modal value and reflects sympathetic heart rhythm modulation;

5) Stress index (SI) calculated as: SI = AMo / (2 × Mo × MxDMn), where AMo is expressed as percentage, Mo and MxDMn in seconds; reflects sympathetic tone.

These parameters were compared with synchronous time series of minute-by-minute values of the X and Y components of the geomagnetic field (GMF) vector (in nanotesla) recorded at the geophysical station closest to the measurement site.

Geophysical Data

For each recording site, the X and Y components of the GMF vector were obtained at 1-minute resolution from geophysical stations located as close as possible to the corresponding experimental location: Borok station (BOXX; 58.070°N, 38.230°E) for the Moscow region (55°45´N/ 37°36´E); Nurmijärvi station (NUR; 60.500°N, 24.600°E) for the Leningrad and Arkhangelsk regions (59°57´N/ 30°19´E). All geomagnetic data were obtained via the INTERMAGNET network (https://imag-data.bgs.ac.uk/GIN_V1/GINForms2).

We focused on horizontal GMF components because their spatial variability with distance is relatively low, a fact we verified separately [26]. In contrast, minute-by-minute variations in the vertical (Z) component are highly dependent on local ground conditions at the measurement site. Therefore, in cases where the distance between the biological recording site and the geophysical station was considerable, we considered the use of Z-component data—and, consequently, the full GMF vector—to be inappropriate. However, in earlier papers where biological measurements were conducted in close proximity to geophysical stations, variations in the vertical component and full vector were included in the analysis [24, 27, 31].

While earlier investigations of the biogeophysical synchronization effect used HRV data recorded exclusively under geomagnetically undisturbed conditions [24, 26, 27, 31], this study did not differentiate by levels of geomagnetic disturbance. This decision was based on prior findings [25] indicating that the frequency of synchronization does not depend on the GMF disturbance level, as assessed by daily Kp-index values. Furthermore, the physiological parameters compared in this study (HR and HRV indices) were measured under identical space weather conditions.

Data Analysis Algorithm

All calculations were performed in MATLAB R2018 using built-in functions and custom-developed applications.

The analysis combined cross-correlation and wavelet-based approaches. The full algorithm is described in detail in our previous work [25].

Before analysis, both physiological and geophysical time series were preprocessed using a bandpass filter to remove trends and ultra-low-frequency oscillations.

Correlation Analysis: Since both biological and geophysical time series often failed to meet the normality criterion, we used Spearman's rank correlation coefficient to assess the strength of correlation. This metric is robust in deviations from normal distribution.

As previously noted, one manifestation of the biogeophysical synchronization effect is the simultaneous presence of quasi-periodic oscillations with similar frequencies in both time series, with an a priori unknown phase shift. This formed the alternative hypothesis (H1), while the null hypothesis (H0) stated that no association exists between the series. To detect such relationships, we calculated correlation coefficients between the biological and geophysical series at time lags ranging from −5 to +5 minutes (11 total lags). The highest absolute correlation coefficient was selected and its p-value was calculated.

To address the increased probability of false positives due to multiple testing of 11 time lags instead of one, we applied the Bonferroni correction. The Bonferroni correction method states that to reduce the likelihood of false-positive results, hypotheses should be rejected if the p <α/m, where m is the number of hypotheses tested (in this case, m=11). This correction ensures that the family-wise error rate (FWER) remains below α, as derived from Boole's inequality, which holds that the probability of at least one event occurring in a finite or countable set of events does not exceed the sum of the individual event probabilities. Accordingly, if each individual test is evaluated at a significance level of α/m, the overall significance level for the family of hypotheses is maintained at α. Therefore, correlation coefficients are considered statistically significant if p <α/m=0.05/11=0.0045.

Because the length of time series varied between 60 and 120 values across experiments, direct comparison of correlation coefficients would have been inappropriate. Instead, we compared p-values. To facilitate analysis and graphical representation, a logarithmic transformation of the p-value accounting for the sign of the correlation coefficient was used: Ks=–sign(r_s)×lg(p). This format offers several advantages over the traditional reporting of paired r_s and p-values, particularly when analyzing large datasets. First, using Ks allows for a single composite metric instead of two. Second, it enables comparisons across time series of different lengths. Third, Ks increases (rather than decreases, as p does) with stronger correlation, which is more intuitive. This transformation simplifies result interpretation without loss of data, as there is a one-to-one correlation between Ks and the original r_s and p-values. In this context, Ks values greater than 1.3 or less than –1.3 (where 1.3=–log(0.05)) indicate statistically significant positive or negative correlations, respectively, at the p <0.05 level. Values of |Ks| >2 correspond to p <0.01, while |Ks| <1.3 indicates no statistically significant correlation. In this study, applying the Bonferroni correction yielded a critical threshold of |Ks|=–log(0.0045)=2.35, corresponding to α=0.0045.

Wavelet Spectrum Similarity Analysis: For each of the 673 analyzed experiments, the time series of HR and HRV parameters, as well as the X and Y components of the GMF, were processed according to the following algorithm:

1) Wavelet coefficient matrices W(h)_i, W(x)_i, W(y)_i were computed for each experiment i=1...673. These matrices represented spectral power density values and had dimensions of 50×D_i, where 50 is the number of tested periods ranging from 1 to 50 minutes, and D_i is the duration of the i^th experiment in minutes. A standard complex Morlet wavelet function was used for the transformation.

2) From the resulting matrices W(h)_i, W(x)_i, W(y)_i, the mean spectral power across periods was computed by averaging the values in each row (1 to 50). This yielded vectors [h]_i, [x]_i, [y]_i, sized 1×50, reflecting the intensity of each period in the HR, X, and Y series, respectively, for the i^th experiment.

3) To quantify the similarity or difference between the sets of periods represented in the wavelet spectra for a given pair of time series (e.g., HR–Y), the scalar product of the normalized vectors [h]_i and [y]_i was calculated: Qy_i=(h_i,y_i)/|h_i|∙|y_i|.

Mathematically, the value of the Q_y parameter is equivalent to the cosine of the angle between the vectors [h] and [y], or the correlation coefficient between them. However, because adjacent values in these vectors are not independent, standard methods for assessing statistical significance are not applicable. Therefore, the threshold for considering the Q_x and Q_y parameters as indicating directional similarity—and thus spectral similarity—was empirically set at Q ≥0.4.

Results

Fig. 1 and Fig. 2 illustrate the key steps of the two applied analysis algorithms—correlation analysis and wavelet spectrum similarity evaluation—using a single experiment as an example (volunteer V2, recording started on June 11, 2013, at 07:00 UT).

 

Fig. 1. Illustration of the correlation-based method for assessing synchronization of physiological parameters—heart rate (HR), RMSSD, and SI—with variations in the X component of the geomagnetic field (GMF): (a), superimposed raw time series of physiological parameters (red) and the horizontal GMF component from the Borok geophysical station (BOXX, blue); (b), superimposed filtered time series; (c), cross-correlation functions between values of each physiological parameter and the GMF component. Ks=–log₁₀(p)×sign(r), where r is the Spearman rank correlation coefficient and p is its statistical significance level. The red dashed line indicates the threshold of statistical significance at p=0.0045 (|Ks| >2.35).

  

Fig. 2. Illustration of the wavelet spectrum comparison method. Left: wavelet spectra of BOXX geomagnetic field, heart rate (HR), RMSSD, and SI time series. Right: mean spectra of corresponding series along the ordinate axis.

 

Table 2 presents the numerical results of time series comparison using both methods. Time series were considered synchronous if the correlation analysis yielded |Ks| >2.35 or if the wavelet similarity criterion produced Q_x >0.4.

 

Table 2. Example of results assessing the similarity between time series of physiological parameters and the geomagnetic field vector in the experiment shown in Figs. 1 and 2

Physiological parameters	|Ks|	Qx
HR	6.53	0.522
RMSSD	3.65	0.574
SI	2.43	0.472

Note: HR, heart rate; RMSSD, root mean square of successive differences; SI, stress index.

 

As shown in Fig. 2, all four time series demonstrate peak spectral power at a period of 18–19 minutes, with an additional smaller peak at approximately 9–10 minutes. Table 2 shows a monotonic decrease in correlation strength with GMF from HR to SI. In the case of SI, the correlation coefficient is close to the statistical significance threshold. Conversely, based on wavelet spectrum similarity, the degree of alignment appears roughly equal across the three analyzed physiological parameters.

Fig. 3 presents the results of all 673 experiments analyzed using cross-correlation (Fig. 3a) and wavelet spectrum comparison (Fig. 3b). The y-axis indicates the frequency of synchronization detection, N (i.e., the relative number of experiments in which a given physiological parameter was synchronized with a specific GMF component: N=N_k/n, where N_k is the number of experiments showing synchronization per the respective criterion, and n = 673 is the total number of analyzed experiments).

 

Fig. 3. Cumulative distribution of the frequency of biogeophysical synchronization between heart rate (HR) and heart rate variability (HRV) parameters with each horizontal component of the geomagnetic field (GMF) across all experiments: (a), cross-correlation analysis; (b), wavelet spectral similarity analysis. *p <0.05; **p <0.01; ***p <0.001. Asterisks next to the HRV parameter bars indicate the level of statistical significance for differences in synchronization frequency between HR and the respective HRV parameter with each GMF component.

 

As shown in Fig. 3a, the frequency of synchronization events (N) between HR and each selected component of the GMF, as determined using the correlation analysis method, was approximately 32%, while for HRV parameters it ranged from 9% to 17%, representing a difference of 2-fold or more. According to the wavelet spectrum similarity criterion, HR synchrony with GMF components was observed in 40% of cases, whereas HRV parameters showed synchrony in 24%–28% of experiments.

According to the χ² test, the synchronization frequency N for HR in both analytical methods differed significantly (***, p <0.001) from that of each of the four HRV parameters, and this result held for both GMF components.

Fig. 4 displays partial sample distributions of synchronization frequency values (N), derived from correlation analysis for each of the eight volunteers in group 1, similar to the full-sample distribution in Fig. 3a. Fig. 5 presents the corresponding distributions obtained via wavelet spectrum comparison (similar to Fig. 3b). Fig. 6 presents the distributions of N values obtained from both analytical methods for the 39 volunteers in group 2. Collectively, the distributions in Fig. 4–6 constitute the corresponding distributions presented in Fig. 3.

 

Fig. 4. Sample distributions of synchronization frequency between heart rate (HR) and heart rate variability (HRV) parameters with geomagnetic field components for group 1 volunteers using the correlation method. Legend is identical to that in Fig. 3.

 

Fig. 5. Sample distributions of synchronization frequency between heart rate (HR) and heart rate variability (HRV) parameters with geomagnetic field components for group 1 volunteers using the wavelet spectrum comparison method. Legend is identical to that in Fig. 3.

  

Fig. 6. Sample distributions of the frequency of synchronization events between heart rate (HR) and heart rate variability (HRV) parameters with components of the geomagnetic field for group 2 volunteers: (a), cross-correlation analysis; (b), wavelet spectrum comparison method. Legend is identical to that in Fig. 3.

 

As shown in Table 1, there is considerable heterogeneity in the number of experiments per individual in group 1: volunteer V1 contributed n = 333 experiments, whereas volunteers V5, V6, and V8 each contributed n=10. Our experience suggests that this is the minimum sample size per individual sufficient to reveal certain trends, if not consistent patterns—especially when similar results are observed across multiple individuals. This imbalance results in substantially greater variability in the distributions shown in Fig. 4 for volunteers V4–V8 compared to V1–V3. Nevertheless, individual-level distributions allow us to assess the extent to which the conclusions derived from the full dataset (Fig. 3) are reproduced when analyzing its independent, non-overlapping subsets.

In Fig. 4, the χ² test indicates that, for six volunteers (V1, V2, V3, V5, V6, V8), the frequency of synchronization events (N) for HR with at least one GMF component is statistically significantly higher (p < 0.05) than the corresponding N values for any HRV parameter. For the remaining two volunteers, the same trend is observed; however, the small sample size prevents the results from reaching statistical significance.

A downward trend in N values from RMSSD to SI was also observed in volunteers V1–V6, consistent with the overall distribution shown in Fig. 3a.

Analysis of the distributions in Fig. 5 indicates that statistically significant differences in N between HR and the HRV parameters were observed in volunteers V1–V4, all of whom had large or relatively large experimental sample sizes. For volunteer V7, the p-value did not reach the 0.05 threshold; however, the N value for HR was still higher than for the other physiological parameters. In contrast to the distributions in Fig. 4, no excess in HR synchronization frequency (N) relative to other physiological parameters was observed in volunteers V5, V6, and V8, most likely due to the limited number of experimental observations.

A similar pattern is observed in Fig. 6, which presents the analysis results of CIG recordings from volunteers in group 2. According to the correlation analysis (Fig. 6a), the frequency of HR synchronization events with both GMF vector components was significantly higher than that of the HRV parameters (p <0.05). At the same time, the N value for HR was higher than for the other parameters in the frequency distributions of N based on wavelet spectrum similarity (Fig. 6b); however, this difference did not reach statistical significance.

When comparing the various N values in Fig. 3a and 3b, it can be seen that in the former, the synchronization frequency for HR with each GMF component was 1.9–3.5 times higher than for each of the four HRV parameters, while in the latter, the ratio ranged from 1.5 to 1.8. Thus, the correlation criterion revealed more pronounced differences between HR and HRV parameters than did the spectral similarity criterion.

Nevertheless, across the overall cumulative distribution and each of the analyzed individual and group-level subsets, the same conclusion was consistently obtained: the frequency of HR synchronization with GMF variations was significantly higher than that observed for any of the four analyzed HRV parameters. No statistically significant differences were found among the HRV parameters (RMSSD, SDNN, AMo, SI).

Discussion

In this study, correlation analysis revealed that the frequency of HR synchronization with each GMF component was 32% for the entire dataset (673 recordings), while for HRV parameters it ranged from 9% to 17%, indicating a ≥2-fold difference in corresponding frequencies. Based on the wavelet spectrum similarity criterion, HR synchronization with GMF component variations was observed in 40% of cases, and HRV parameter synchronization in 24%–28%. Statistically significant differences in synchronization frequency were also identified in separate experimental subsamples, both in the longitudinal study design (repeated measurements in each of the eight volunteers) and in the cross-sectional design (single measurements in a group of 39 volunteers). In some cases, a gradual decline in synchronization frequency was observed across the sequence of HRV parameters: RMSSD–SDNN–AMo–SI.

The consistent detection of synchronization with each GMF component in approximately 35%–40% of experiments across various subsamples suggests that, on the one hand, the observed association between HR time series and the GMF vector is not random and that the biogeophysical synchronization effect is indeed real. On the other hand, the currently built phenomenological model of the effect may be incomplete and may lack important factors or include extraneous elements that obscure the signal, or both. The next step is to progressively refine the current model to achieve a more accurate characterization of the effect, with the goal of reproducing it under laboratory conditions for controlled investigation. Three main avenues for further investigation can be identified: studying the dynamics of heart rate regulation processes, analyzing the spectral characteristics of GMF variations, and refining the frequency–time parameters of the analytical algorithm.

The present study addresses the first of these directions and is based on the working hypothesis that a certain rhythmic process exists in the human body—a mediator process (possible examples are discussed below)—that is sensitive to GMF variations and involved in heart rate regulation [25]. Naturally, multiple such processes may exist and may be integrated into the regulatory system either sequentially or in parallel. Within this working model of synchronization, the instability in detecting the effect may be attributed to internal regulatory processes of the body: the greater the current contribution of this magnetosensitive mediator process to HR regulation, the stronger the observed synchronization between HR and the GMF vector. The intervals of synchronization and desynchronization may alternate in a quasiperiodic or nearly random fashion and may last from minutes to hours or days. These intervals may include a circadian component or depend on the presence of a third factor. Identifying these specific features of the effect is the objective of future research.

When comparing the present findings with earlier results reported by other authors, it is crucial to note that in the vast majority of studies on the ANS sensitivity to geomagnetic variations [1, 8, 12, 13], researchers used conventional 5-minute HRV recordings, typically performed once per day. From these recordings, a single value of HR and HRV indices was calculated for each experiment and then compared with global geomagnetic disturbance levels (e.g., Kp and Ap indices, integrated intensity of the first Schumann resonance, etc.). Because measurements were conducted once daily, the effective temporal resolution of HRV assessments in those studies corresponded to the daily timescale, and the data were discrete. In that context, the parameter of interest was the change (shift) in the mean value of HRV indices in response to varying levels of GMA. In contrast, in our study, HRV parameters and GMF component values were computed every minute, with 60–120 observations per experiment. Moreover, the effect under investigation was not a shift in the mean physiological parameter but rather a frequency adjustment of oscillatory dynamics.

Therefore, our experiments addressed much higher-frequency (ultradian) and lower-amplitude manifestations of the ANS response to GMF variations than those explored in earlier studies. Based on the characteristic timing of responses in each case, it can be hypothesized that the minute-resolution effects observed in our data represent one of the early phases in the development of this physiological response. Meanwhile, the larger-scale and longer-term responses observed during geomagnetic storms—such as significant shifts in HRV indices at the daily level—are indicative of systemic physiological adaptations associated with specific and nonspecific stress responses. Hence, direct comparison of findings should only be made with studies that used the same or comparable temporal resolution of data acquisition.

Vasin et al. [15] conducted experiments on the effects of millihertz-range magnetic fields (f₁ = 1.67 mHz and f₂ = 1.11 mHz) in healthy volunteers at rest. Participants were exposed to the magnetic field for 1 hour, and HRV parameters were calculated over successive 5-minute intervals. Changes in various HRV metrics were assessed using two approaches: first, by evaluating shifts in the mean values resulting from magnetic field exposure; second, by analyzing changes in the spectral power density of each HRV parameter within the 0.833–3.333 mHz frequency band (corresponding to oscillation periods of 5–20 minutes), which is close to Pc5–Pc6 geomagnetic pulsations.

Analysis of changes in mean HRV values revealed that exposure most notably affected pNN50, SDNN, LF/HF, and VLF. Thus, the artificial magnetic field, with frequency and amplitude characteristics resembling natural GMF variations, induced statistically significant shifts in mean values similar to those observed during geomagnetic storms.

However, for our purposes, the more relevant findings from this experiment pertain to changes in the spectral power density of various HRV parameters, as these directly relate to our results. First, the spectral range assessed (5–20 minutes) overlaps with the frequency band analyzed in our study for spectral coincidence between HR and GMF variations (3–40 minutes). Second, one of the applied frequencies (f₁ = 1.67 mHz) corresponds to the 10-minute period we previously identified [25, 32], around which the synchronization effect between HR and minute-scale GMF variations was most pronounced. Finally, the authors of paper [30] reported an increase in spectral power for only two HRV parameters: meanNN (equivalent to HR) and LF/HF, whereas RMSSD and SDNN showed slight decreases, and AMo and SI were not included in their analysis. Because spectral HRV parameters were not examined in our study, the conclusions align fully with the findings of Vasin et al. [15]: among time-domain HRV parameters in the 5- to 20-minute oscillatory range, increased spectral power was observed for meanNN (HR) relative to RMSSD and SDNN.

In the present study, spectral HRV parameters and their synchronization with GMF rhythmicity were not analyzed; this remains a subject for future investigation. However, previous experiments involved 30-minute HRV recordings in groups of healthy volunteers and individuals with impaired vascular tone (i.e., elevated or decreased blood pressure), followed by assessment of synchronization frequency between GMF components and various HRV metrics in each group. It was found that the frequency of statistically significant correlations between GMF components and HRV parameters was higher in the group with blood pressure dysregulation compared with the healthy group. This difference was most pronounced for HRV indices reflecting vagal activity (RMSSD and HF) [31]. A similar conclusion was drawn when analyzing only the subgroup of individuals with arterial hypertension within the blood pressure dysregulation group: synchronization between GMF component variations and HRV fluctuations was more frequent in hypertensive participants than in healthy individuals, particularly for HF (60% vs 8.7%, respectively; p < 0.05) and RMSSD (50% vs 13%, respectively; p < 0.05) [33]. Conversely, synchronization frequency with GMF variations for the LF parameter, which reflects baroreflex activity, was significantly higher in participants with normal blood pressure compared with those with hypertension [27]. Thus, the synchronization frequency of HRV parameters with GMF variations observed in healthy individuals in our earlier papers—up to 20%, with a slight predominance for RMSSD—is consistent with the values obtained in the present research. In this context, the observed synchronization of the LF parameter, which reflects baroreflex activity, with GMF component variations was interpreted as adaptive. In this context, other HRV parameters retained greater independence from GMF fluctuations in individuals with normal vascular tone than in those with dysregulation (i.e., hypertension).

It is important to note that LF power reflects baroreflex activity mediated by baroreceptors, which feature ion channels with piezoelectric properties (particularly Piezo2), as well as vagal influences [34, 35]. Some researchers refer to aortic baroreceptors as “low-pressure” baroreceptors. Mechanical stretch impulses are transmitted to the right atrium, where they initiate the mechanism of cardiac contraction [36]. Theoretically, oscillations in baroreceptor activity may be modulated by GMF variations through the modulation of subthreshold membrane potential oscillations, which in turn influence sinoatrial node function and thereby HR.

Subthreshold membrane potential oscillations exhibit intrinsic rhythmicity. Under certain conditions (e.g., inflammation, metabolic disturbances), the rhythmic bursting activity of these oscillations changes, triggering action potentials. Such intrinsic rhythmicity has been identified in brain cells, including circadian neurons of the suprachiasmatic nucleus and the retrotrapezoid nucleus of the brainstem [37]. This oscillatory activity is mediated by the transient receptor potential cation channel subfamily M member 4 (TRPM4). TRPM4 is involved in subthreshold oscillations that support pacemaker activation of neurons in the retrotrapezoid nucleus of the brainstem, which is essential for basal respiratory activity, CO₂-stimulated breathing, and state-dependent respiratory control. This receptor is also present in cardiomyocytes and plays a critical role in the regulation of bioelectrogenesis in the myocardium [38]. It is therefore hypothesized that a resonance effect may occur between GMF variations and subthreshold membrane potential oscillations in excitable structures, mediated through altered ion channel activity via TRPM4 receptors directly in the myocardium. In this scenario, autonomic nervous system activity as an intermediary in sinoatrial node excitation during GMF variations may be minimal. An alternative mechanism may involve the transmission of excitation to heart rhythm through TRPM4 receptor activity changes in the brain. Arterial and cardiopulmonary baroreceptor afferents also converge within the vagus and glossopharyngeal nerves, transmitting signals to autonomic control centers in the brainstem. This could result in longer-period oscillations in both HR and HRV parameters.

Another promising mechanism for direct GMF influence on myocardial bioelectrogenesis involves ephaptic (non-synaptic) signal transmission between excitable structures in the myocardium—a phenomenon demonstrated in the mammalian neocortex. It has been hypothesized that biomagnetic fields of astroglia, associated with transient changes in Ca²⁺ concentrations, may participate in ephaptic neuronal communication through direct magnetic modulation of intercellular local field potentials [39]. Ephaptic impulse transmission in myocardial cells in vitro was demonstrated as early as the 1980s [40] and remains an area of active research. It is believed that electrical impulses can propagate to neighboring cardiomyocytes not only via gap junctions (nexuses between sarcolemmas) but also through specialized adjacent spaces known as perinexi [41]. These bioelectrical processes occurring within perinexal spaces are considered the basis for ephaptic impulse transmission between cardiomyocytes in the myocardium. Given that perinexal structures in the myocardium are relatively distant from sympathetic and vagal regulatory inputs, HR responses to external electromagnetic influences may be mediated by ephaptic mechanisms that bypass autonomic nervous system modulation. This may occur within specific frequency ranges.

Thus, the present findings may be interpreted as follows. The HRV parameters calculated over 1-minute intervals may offer only limited insight into ANS function. Autonomic regulation of heart rhythm may be more sensitively detected through alternative approaches, such as rhythmic assessments of catecholamine and acetylcholine secretion into synaptic clefts or systemic circulation. These secretory rhythms may influence HR oscillations, which in turn synchronize with GMF variations. However, such measurements are technically challenging to implement in human experimental settings.

A promising avenue for future research lies in extending synchronization analysis to spectral HRV parameters (HF, LF, VLF, LF/HF), including assessments using 2- to 3-minute HRV segments. These would allow more appropriate estimation of moving averages for spectral HRV indices [27]. This analysis will require substantial refinement of both the algorithm used to transform RR interval series into spectral HRV time series and, potentially, revalidation of the similarity assessment algorithm parameters to accommodate the characteristics of these new time series.

A key limitation of this study was the exclusive analysis of time-domain HRV parameters. Spectral analysis of low-frequency (LF) and very low-frequency (VLF) HRV components requires specialized mathematical preprocessing of cardiointervalograms.

Therefore, future research directions may be outlined as follows:

1. Analyze the potential manifestation of the biogeophysical synchronization effect in spectral HRV parameters by comparing the likelihood of spectral overlap between HRV indices and GMF components.
2. Investigate potential spectral power redistribution across different HRV frequency bands under varying geomagnetic conditions.
3. Develop protocols and conduct experiments for time-series recording of biochemical markers that reflect autonomic balance.
4. Analyze how the occurrence of the biogeophysical synchronization effect depends on specific geomagnetic conditions.

Based on these investigations, further refinement is warranted in methodologies for assessing individual cardiovascular magnetic sensitivity, considering health status (presence or absence of hypertension). Moreover, it is necessary to improve models of biotropic oscillatory periods in HR and related physiological parameters to enhance the simulation of cardiac regulatory mechanisms under external electromagnetic field variations.

Conclusion

The method we developed for studying the synchronization of HR with GMF variations in the millihertz frequency range represents an effective tool for testing various hypotheses within the broader fundamental investigation of physiological pathways involved in the organism’s response to low-intensity external factors.

The analyzed time-domain HRV parameters, as indicators of autonomic regulation of heart rhythm, demonstrated lower sensitivity in detecting statistical associations with GMF parameters compared with HR. Nevertheless, the autonomic nervous system may function as a mediating link in the influence of GMF on HR fluctuations, potentially through the rhythmic activity of arterial and cardiopulmonary baroreceptors, as well as the rhythmic secretion of catecholamines and acetylcholine into synaptic clefts and systemic circulation—an assumption that warrants further experimental validation.

Additional information

Authors’ contribution. T.A. Zenchenko — development of the research concept, data analysis, preparation and writing of the article; L.V. Poskotinova — development of the research concept, data collection, editing of the article; N.I. Khorseva — data collection, editing of the article; T.K. Breus — literature review, collection and analysis of literary sources, writing of the text and editing of the article. All authors confirm that their authorship meets the international ICMJE criteria (all authors have made a significant contribution to the development of the concept, research and preparation of the article, read and approved the final version before publication).

Ethical expertise. The study was approved by the Bioethics Committee of the ITEB RAS (Protocol No. 06/2012 dated 06/01/2012).

Consent for publication. All study participants voluntarily signed an informed consent form before being included in the study.

Funding sources. The work was carried out within the framework of the State Assignment of ITEB RAS No. 075-00223-25-02, State Assignment of IKI RAS, topic "Plasma", State Assignment of IBCP RAS (44.1. state topic number 0084-2019-004) and State Assignment of FGBUN FITC RAS Ural Branch No. FUUW-2025-0032.

Disclosure of interests. The authors have no relationships, activities or interests for the last three years related with for-profit or not-for-profit third parties whose interests may be affected by the content of the article.

Statement of originality. In creating this work, the authors did not use previously published information (text, illustrations, data).

Data availability statement. The editorial policy regarding data sharing does not apply to this work, and no new data was collected or created.

Generative AI. Generative AI technologies were not used for this article creation.

Provenance and peer-review. This paper was submitted to the journal on an unsolicited basis and reviewed according to the usual procedure. Two external reviewers, a member of the editorial board, and the scientific editor of the publication participated in the review.

Acknowledgments. The results presented in this paper were obtained using geophysical data collected by Nurmijarvi and Borok observatories. The authors thank the Finnish Meteorological Institute and Borok Geophysical Observatory for providing the data and for their work within the INTERMAGNET project to disseminate high standards of geophysical observations. The authors also thank M.E. Diatroptov, A.A. Stankevich and A.E. Elfimova for their assistance in data collection.

作者简介

Tatiana A. Zenchenko

Institute of Theoretical and Experimental Biophysics of the Russian Academy of Sciences; Space Research Institute of Russian Academy of Sciences

编辑信件的主要联系方式.
Email: zench@mail.ru
ORCID iD: 0000-0002-0520-2029
SPIN 代码: 8974-6685

Dr. Sci. (Biology), Cand. Sci. (Physics and Mathematics)

俄罗斯联邦, 3 Institutskaya st, Pushchino, Moscow region, 142290; Moscow

Liliya V. Poskotinova

N. Laverov Federal Center for Integrated Arctic Research of the Ural Branch of the Russian Academy of Sciences

Email: liliya200572@mail.ru
ORCID iD: 0000-0002-7537-0837
SPIN 代码: 3148-6180

Dr. Sci. (Biology) MD, Cand. Sci. (Medicine), Assistant Professor

俄罗斯联邦, Arkhangelsk

Nataliya I. Khorseva

Institute of Biochemical Physics of the Russian Academy of Sciences

Email: sheridan1957@mail.ru
ORCID iD: 0000-0002-3444-0050

Cand. Sci. (Biology)

俄罗斯联邦, Moscow

Tamara K. Breus

Space Research Institute of Russian Academy of Sciences

Email: breus36@mail.ru
ORCID iD: 0000-0003-4057-0844
SPIN 代码: 1267-8561

Dr. Sci. (Physics and Mathematics)

俄罗斯联邦, Moscow

参考

补充文件