Skip to main content

Systemic lupus Erythematosus and geomagnetic disturbances: a time series analysis

Abstract

Background

To examine the influence of solar cycle and geomagnetic effects on SLE disease activity.

Methods

The data used for the analysis consisted of 327 observations of 27-day Physician Global Assessment (PGA) averages from January 1996 to February 2020. The considered geomagnetic indices were the AP index (a daily average level for geomagnetic activity), sunspot number index R (measure of the area of solar surface covered by spots), the F10.7 index (measure of the noise level generated by the sun at a wavelength of 10.7 cm at the earth’s orbit), the AU index (upper auroral electrojet index), and high energy (> 60 Mev) proton flux events. Geomagnetic data were obtained from the Goddard Space Flight Center Space Physics Data Facility. A time series decomposition of the PGA averages was performed as the first step. The linear relationships between the PGA and the geomagnetic indices were examined using parametric statistical methods such as Pearson correlation and linear regression, while the nonlinear relationships were examined using nonparametric statistical methods such as Spearman’s rho and Kernel regression.

Results

After time series deconstruction of PGA averages, the seasonality explained a significant fraction of the variance of the time series (R2 = 38.7%) with one cycle completed every 16 years. The analysis of the short-term (27-day) relationships indicated that increases in geomagnetic activity Ap index (p < 0.1) and high energy proton fluxes (> 60 Mev) (p < 0.05) were associated with decreases in SLE disease activity, while increases in the sunspot number index R anticipated decreases in the SLE disease activity expressed as PGA (p < 0.05). The short-term correlations became statistically insignificant after adjusting for multiple comparisons using Bonferroni correction. The analysis of the long-term (297 day) relationships indicated stronger negative association between changes in the PGA and changes in the sunspot number index R (p < 0.01), AP index (p < 0.01), and the F10.7 index (p < 0.01). The long-term correlations remained statistically significant after adjusting for multiple comparisons using Bonferroni correction.

Conclusion

The seasonality of the PGA averages (one cycle every 16 years) explains a significant fraction of the variance of the time series. Geomagnetic disturbances, including the level of geomagnetic activity, sunspot numbers, and high proton flux events may influence SLE disease activity. Studies of other geographic locales are needed to validate these findings.

Peer Review reports

Introduction

The Earth is immersed in a geomagnetic field (GMF) that deflects the flowing ionized particles of solar winds. In doing so, the GMF is impressed and altered by the solar winds [1]. These alterations in GMF are called geomagnetic disturbances. The GMF has different strength and direction in different areas of the Earth. Its magnitudes vary from ~ 35,000 nano Tesla (nT) at the equator to ~ 70,000 nT at the magnetic poles [2]. Several periodic variations have also been shown in GMF, mostly related to various solar cycles (especially the 11 year cycle), and others related to Earth’s rotation and situation in its orbit [2].

The idea that variations in solar activity can influence biological phenomena originated with the work of Alexander Chizhevsky in the early twentieth century. He carried out a systematic study of the influence of electric, magnetic and electromagnetic perturbations in the external physicochemical environment on biological, social, and historical events [3]. The advent of space-based observatories such as the joint USA/European Solar and Heliospheric Observatory during the late twentieth century and early twenty-first century facilitated further research into the effect of variations in the sun and earth’s magnetic field on human health. High solar activity associated with excess death from myocardial infarction [4]. Geomagnetic disturbances were directly correlated with erythrocyte aggregation studies and blood velocity of patients with acute myocardial infarction [5] and with blood pressure variability [6]. A strong non-linear correlation between sudden infant death syndrome and geomagnetic activity was described in Ontario, Canada [7]. Geomagnetic disturbances have been postulated to explain the relapsing-remitting nature, chronobiology, and latitudinal prevalence of multiple sclerosis [8]. Giant cell arteritis and rheumatoid arthritis incidence rates were shown to be highly correlated with the AL index, a proxy of the westward auroral electrojet and a measure of geomagnetic activity [9].

We previously described significant seasonal variation in SLE disease activity with more arthritis activity in the spring and summer months, and an increase in renal activity in winter months, significantly higher anti-dsDNA antibody titers in the fall, and a significant variation of global disease activity as measured by SELENA-SLEDAI through the year [10]. We also performed a spatial-time cluster analysis of the Hopkins Lupus Cohort and detected space-time organ specific lupus flare clusters which had multi-year cluster patterns that did not conform to any known pattern of infectious disease or environmental exposure [11].

The only previous study assessing the effect of geomagnetic disturbances in SLE analyzed the fluctuation of antiphospholipid antibodies over time and found that the anti-beta 2 glycoprotein (p < 0.0001) and lupus anticoagulant (p < 0.05) were more commonly detected on days of major geomagnetic storms [12].

We thus hypothesized that geomagnetic disturbances influence SLE disease activity.

Methods

Patients and activity indices

As previously described [13], the Hopkins Lupus Cohort is a prospective cohort study of predictors of lupus flare, atherosclerosis, and health status in SLE. The study cohort includes all patients at the Hopkins Lupus Center who have a clinical diagnosis of SLE and have given informed consent to participate in the study. Subjects enrolled in the cohort are followed quarterly or more frequently if clinically necessary. The clinical features, laboratory testing, and damage accrual data are recorded at the time of entry into the cohort and are updated at each subsequent visit. The Hopkins Lupus Cohort has been approved by the Johns Hopkins University School of Medicine Institutional Review Board and complies with the Health Insurance Portability and Accountability Act. All patients gave written informed consent.

Ninety-five percent of patients fulfilled 4 or more of the American College of Rheumatology (ACR) 1982 revised classification criteria for SLE [14, 15] and/or the SLICC classification criteria for SLE [16]. Disease activity was measured with the Physician Global Assessment (PGA) [17]. The data used for the analysis consisted of 327 observations of 27-day PGA averages from 24 January 1996 to 29 February 2020. A 27-day PGA average was used for analysis due to the specific periodicity of the OMNI dataset [18].

The OMNI data were obtained from the Goddard Space Flight Center/Spatial Physics Data Facility OMNIWeb interface at https://omniweb.gsfc.nasa.gov. OMNI is an open-source 20-spacecraft compilation of near-Earth, solar wind magnetic field and plasma (proton and, since 1971, alpha particle) data, energetic proton fluxes, and geomagnetic activity indices (Kp, Dst, AE) and sunspot number. We used the Ap index, sunspot number index R, the F10.7 index, the AU index, and high energy (> 60 Mev) proton flux events for our analyses. Detailed descriptions of these geomagnetic indices can be located in supplementary materials.

Statistical methodologies

As a first step we studied the components of the 27-day PGA average time series, specifically its trend, seasonality and noise. The trend was assumed to be a linear function of time, and its parameters were estimated by Ordinary Least Squares (OLS). The seasonality was assumed to be a sinusoidal function of time, and its parameters were also estimated by OLS. The behavior of the noise was consistent with an autoregressive moving average (ARMA) process, and its parameters were obtained by Maximum Likelihood Estimation (MLE).

As a second step we investigated the relationship between changes in the 27-day average PGA and changes in the 27-day averages of the considered geomagnetic indices. We examined both the contemporaneous relationships (specifically, the relationship between changes in the PGA and changes in the geomagnetic indices over the same 27-day period) and the lead-lag relationships (specifically, the relationship between changes in the PGA over a given 27-day period and changes in the geomagnetic indices over the previous 27-day period). We used Pearson’s linear correlation coefficient and univariate linear regression to investigate the presence of linear relationships, while we use Spearman’s rank correlation coefficient and univariate Kernel regression to investigate the presence of nonlinear relationships.

Results

Physician global assessment (PGA) time series decomposition

The plot of the 27-day PGA average time series is shown in panel (a) of Fig. 1.

Fig. 1
figure1

27-day PGA average time series decomposition. a 27-day PGA average. b 27-day PGA average and fitted trend component. c 27-day PGA average and fitted trend and seasonality components. d 27-day PGA average and fitted trend, seasonality, and noise components

Trend. In what follows yt denoted the value of the 27-day PGA average on period t, with t = 1, 2, …, 327. For estimating the trend, yt was assumed to be described by a linear function of time as in Eq. (1).

$$ {y}_t={\beta}_0+{\beta}_1t+{\eta}_t\kern1.5em {\eta}_t\sim \mathcal{N}\left(0,{\sigma}^2\right) $$
(1)

The parameters in Eq. (1) were estimated by Ordinary Least Squares (OLS). The regression results were reported in panel (a) of Table 1. The estimated intercept (\( {\hat{\beta}}_0 \)) is in line with the sample mean of the time series and was statistically significantly different from zero. The estimated slope (\( {\hat{\beta}}_1 \)) was instead not statistically significantly different from zero. The R2 of the regression was 0.1%. Overall, the regression results indicated that the trend explained only a negligible fraction of the variance of the time series, and therefore the trend was not a significant component of the time series.

Table 1 27-day PGA average time series parameter estimation

In what follows the estimated trend was denoted by

$$ {\hat{T}}_t={\hat{\beta}}_0+{\hat{\beta}}_1t $$
(2)

The plot of the estimated trend is shown in panel (b) of Fig. 1. The Root Mean Square Error (RMSE) obtained when predicting the time series (i.e. yt) using only the estimated trend component (i.e. \( {\hat{T}}_t \)) was 0.084. This indicates that, when using only the estimated trend to predict the lupus disease activity, the prediction error was 0.084 on average.

Seasonality. In what follows st denoted the value of the seasonality function on period t, with t = 1, 2, …, 327. The seasonality function st was assumed to be a sinusoidal function of time as in Eq. (3).

$$ {s}_t={\theta}_0+{\theta}_1\cos \left(2\pi \left(\omega t+\phi \right)\right) $$
(3)

The interpretation of the parameters in Eq. (3) was the following:

  • θ0 was the mean, which represented the level around which the seasonality function fluctuates.

  • θ1 was the amplitude, which represented the width of the fluctuations of the seasonality function.

  • ω was the seasonal frequency, which represented the number of cycles of the seasonality function over a given seasonal period.

  • ϕ was the phase, which defined the time when the seasonality function attained its peak.

The seasonality function in Eq. (3) can be reparametrized as follows

$$ {s}_t={\theta}_0+{\theta}_1\cos \left(2\pi \omega t+2\pi \phi \right) $$
$$ ={\theta}_0+{\theta}_1\left[\cos \left(2\pi \omega t\right)\cos \left(2\pi \phi \right)-\sin \left(2\pi \omega t\right)\sin \left(2\pi \phi \right)\right] $$
$$ ={\theta}_0+\left[{\theta}_1\cos \left(2\pi \phi \right)\right]\cos \left(2\pi \omega t\right)+\left[-{\theta}_1\sin \left(2\pi \phi \right)\right]\sin \left(2\pi \omega t\right) $$
$$ =C+A{x}_{1,t}+B{x}_{2,t} $$

where C = θ0, A = θ1 cos(2πϕ), B =  − θ1 sin(2πϕ), x1, t = cos(2πωt) and x2, t = sin(2πωt). The parameters A, B and C could be estimated by OLS by fitting a linear regression model to the detrended time series (i.e. \( {y}_t-{\hat{T}}_t \)) as in Eq. (4).

$$ \left({y}_t-{\hat{T}}_t\right)=C+A{x}_{1,t}+B{x}_{2,t}+{\epsilon}_t\kern1.5em {\epsilon}_t\sim \mathcal{N}\left(0,{\gamma}^2\right) $$
(4)

The regression results are reported in panel (b) of Table 1. The estimated intercept \( \hat{C} \) was not statistically significantly different from zero, which was expected since the detrended time series has zero mean. The estimated coefficients of the sinusoids \( \hat{A} \) and \( \hat{B} \) were instead both statistically significantly different from zero. The R2 of the regression is 38.7%. Overall, the regression results indicated that, differently from the trend, the seasonality explained a significant fraction of the variance of the time series.

The estimates of θ0, θ1 and ϕ were derived from the estimates of A, B and C as follows [19]

$$ {\theta}_0=\hat{C} $$
$$ {\theta}_1=\sqrt{{\hat{A}}^2+{\hat{B}}^2} $$
$$ \hat{\phi}=\left\{\begin{array}{c}\begin{array}{cc}\frac{1}{2\pi }{\tan}^{-1}\left(-\frac{\hat{B}}{\hat{A}}\right)&\ \hat{A}>0\\ {}\frac{1}{2\pi}\left({\tan}^{-1}\left(-\frac{\hat{B}}{\hat{A}}\right)-\pi \right)& \kern0.5em \hat{A}<0,\hat{B}>0\end{array}\\ {}\begin{array}{cc}\frac{1}{2\pi}\left({\tan}^{-1}\left(-\frac{\hat{B}}{\hat{A}}\right)+\pi \right)& \hat{A}<0,\hat{B}\le 0\\ {}-\frac{1}{4}& \hat{A}=0,\hat{B}>0\end{array}\\ {}\begin{array}{cc}\kern0.5em +\frac{1}{4}&\ \hat{A}=0,\hat{B}<0\\ {}\kern10.5em & \kern0.5em \end{array}\end{array}\right. $$

The estimated mean θ0 was − 0.003, the estimated amplitude θ1 was 0.074, and the estimated phase ϕ was − 0.475.

The seasonal frequency ω was not estimated directly, but was instead inferred by re-estimating the regression model in Eq. (4) using as input different values of ω and by then choosing the value of ω which provided the highest R2. The value of ω which maximized the R2 is 0.062, indicating that the seasonality function completed one cycle every 16.05 years. This is in line with the behavior of the autocorrelation function (ACF) of the detrended time series reported in panel (a) of Fig. 2 (supplementary material), which showed a periodic pattern with a frequency of approximately 220 periods, where the length of each period was 27 days. The most recent seasonal peak was in 2019, and therefore the next seasonal peak is expected in 2035.

Fig. 2
figure2

27-day PGA average time series analysis. a ACF of detrended 27-day PGA average. b ACF of detrended and deseasonalized 27-day PGA average. c ACF of ARMA(1,1) residuals. d Q-Q plot of ARMA(1,1) residuals

In what follows the estimated seasonality function was denoted by

$$ {\hat{S}}_t={\hat{\theta}}_0+{\hat{\theta}}_1\cos \left(2\pi \left(\hat{\omega}t+\hat{\phi}\right)\right) $$
(5)

The plot of the estimated trend and seasonal components (i.e. \( {\hat{T}}_t+{\hat{S}}_t \)) is shown in panel (c) of Fig. 1. The RMSE obtained when predicting the time series (i.e. yt) using both the estimated trend component and the estimated seasonal component (i.e. \( {\hat{T}}_t+{\hat{S}}_t \)) was 0.065. This indicated that, when using both the estimated trend and the estimated seasonality to predict the lupus disease activity, the prediction error was 0.065 on average.

Noise. The detrended and deseasonalized time series followed an ARMA (1, 1) process. The ARMA (1, 1) specification was consistent with the plot of the ACF of the detrended and deseasonalized time series shown in panel (b) of Fig. 2 (supplementary material). The appropriateness of the ARMA (1, 1) model choice was confirmed by the results of a grid search, which showed that the ARMA (1, 1) model had the lowest Akaike Information Criterion (AIC) compared to alternative ARMA(p,q) models. The ARMA process was:

$$ \left({y}_t-{\hat{T}}_t-{\hat{S}}_{\mathrm{t}}\right)={\psi}_1\left({y}_{t-1}-{\hat{T}}_{t-1}-{\hat{S}}_{t-1}\right)+{\lambda}_1{\kappa}_{t-1}+{\kappa}_t\kern1.5em {\kappa}_t\sim \mathcal{N}\left(0,{\mathcal{X}}^2\right) $$
(6)

The parameters of the ARMA (1, 1) model were obtained by Maximum Likelihood Estimation (MLE). The estimation results are reported in panel (c) of Table 1. The estimated autoregressive parameter (\( {\hat{\psi}}_1 \)) and moving average parameter (\( {\hat{\lambda}}_1 \)) were both statistically significantly different from zero and the behavior of the residuals indicated that the model was correctly specified. In particular, the Ljung-Box test indicated that the residuals were serially uncorrelated, while the Jarque-Bera test indicated that the residuals were normally distributed. The ACF plot and Q-Q plot of the residuals are reported in panels (c) and (d) of Fig. 2 (supplementary material). The plot of the estimated trend, seasonal and noise components is shown in panel (d) of Fig. 1.

Relationship between physician global assessment (PGA) and geomagnetic indices

The time series of the AU index had 27 missing values, the time series of Proton Flux (> 60 Mev) had 12 missing values, while all other time series had no missing values.

The results of the KPSS test in (supplementary material) indicated that the levels time series (i.e. the time series of the levels of the variables on each period) were non-stationary, while the differences time series (i.e. the time series of the absolute changes of the variables between one period and the next) were stationary, with the exception of the lupus disease activity time series in which case both the levels and the differences were stationary. In view of the results of the KPSS test, the analysis presented in the following sections was conducted using the differences time series instead of the levels time series. The descriptive statistics of the time series are reported in (supplementary material), while the plots of the time series are shown in Figs. 7, 8 and 9 (supplementary material).

Short-term linear relationships

Table 2 reports the Pearson’s correlation coefficients between changes in PGA and changes in the geomagnetic indices over 1 period (27 days). The Pearson’s correlation coefficient measures the strength of the linear relationship between two variables. The first column reports the correlations of the changes in PGA with the changes in the geomagnetic indices in the same period (i.e. between the change in PGA over period t and the change in each geomagnetic index over the same period t), while the second column reports the correlations of the changes in PGA in the current period with the changes in the geomagnetic indices in the previous period (i.e. between the change in PGA over period t and the change in each geomagnetic index over the previous period t − 1).

Table 2 Short-term (27-day) correlations of the changes in PGA with the changes in the geomagnetic indices in the same period (first column) and of the changes in PGA in the current period with the changes in the geomagnetic indices in the previous period (second column)

The results in Table 2 indicate that, over short time periods:

  • Increases in R (Sunspot No.) could anticipate decreases in the PGA, i.e. positive ΔRt − 1 tend to be associated with negative ΔPGAt (p < 0.05).

  • Increases in the AP index were associated with decreases in the PGA, i.e. positive ΔAPt tended to be associated with negative ΔPGAt (p < 0.1).

  • Increases in the Proton Flux (> 60 Mev) were associated with decreases in the PGA, i.e. positive ΔPFt tended to be associated with negative ΔPGAt (p < 0.05).

Overall the values of the correlation coefficients reported in Table 2 were relatively small, indicating that the linear relationship between changes in PGA and changes in the geomagnetic indices was in general weak. All the correlations in Table 2 are statistically insignificant after adjusting for multiple comparisons using Bonferroni correction (p > 0.1).

Figure 3 (supplementary material) shows the fitted linear regressions. Note that, when only one independent variable was included in the regression, the R2 was equal to the square of the correlation coefficient between the dependent variable and the independent variable, i.e. the R2 reported in Fig. 3 (supplementary material) are the squares of the correlation coefficients reported in Table 2.

Fig. 3
figure3

Short-Term Linear Regressions (27-day differences). a Scatter plot of ΔPGAt vs ΔRt. b Scatter plot of ΔPGAt vs ΔRt-1. c Scatter plot of ΔPGAt vs ΔAPt. d Scatter plot of ΔPGAt vs ΔAPt-1. e Scatter plot of ΔPGAt vs ΔF10.7t. f Scatter plot of ΔPGAtvs ΔF10.7t-1. g Scatter plot of ΔPGAt vs ΔAUt. h Scatter plot of ΔPGAt vs ΔAUt-1. i Scatter plot of ΔPGAt vs ΔPFt. j Scatter plot of ΔPGAt vs ΔPFt-1

Long-term linear relationships

Table 3 reports the Pearson’s correlation coefficients between changes in PGA and changes in the geomagnetic indices over 11 periods (297 days).

Table 3 Long-term (297 days) correlations of the changes in PGA with the changes in the geomagnetic indices in the same period (first column) and of the changes in PGA in the current period with the changes in the geomagnetic indices in the previous period (second column)

The results in Table 3 indicate that, over long time periods:

  • Increases in R (Sunspot No.) are associated to decreases in the PGA, i.e. positive

    ΔRt are associated with negative ΔPGAt (p < 0.01). Furthermore, increases in R (Sunspot No.) anticipate decreases in the PGA, i.e. positive ΔRt − 1 are associated with negative ΔPGAt (p < 0.01). Both correlations are still statistically significant after adjusting for multiple comparisons using Bonferroni correction (p < 0.05).

  • Increases in the AP index are associated with decreases in the PGA, i.e. positive ΔAPt are associated with negative ΔPGAt (p < 0.01). The correlation is still statistically significant after adjusting for multiple comparisons using Bonferroni correction (p < 0.05).

  • Increases in the F10.7 index are associated with decreases in the PGA, i.e. positive ΔF10.7t are associated with negative ΔPGAt (p < 0.01). Furthermore, increases in the F10.7 index anticipate decreases in the PGA, i.e. positive ΔF10.7t − 1 are associated with negative ΔPGAt (p < 0.05). Both correlations are still statistically significant after adjusting for multiple comparisons using Bonferroni correction (p < 0.05 and p < 0.1 respectively).

The fitted linear regression lines and the corresponding R-squared are reported in Fig. 4 (supplementary material). The interpretation of the R-squared is the same as in Section 2.1.

Fig. 4
figure4

Long-Term Linear Regressions (297 day differences). a Scatter plot of ΔPGAt vs ΔRt. b Scatter plot of ΔPGAt vs ΔRt-1. c Scatter plot of ΔPGAt vs ΔAPt. d Scatter plot of ΔPGAt vs ΔAPt-1. e Scatter plot of ΔPGAt vs ΔF10.7t. f Scatter plot of ΔPGAt vs ΔF10.7t-1. g Scatter plot of ΔPGAt vs ΔAUt. h Scatter plot of ΔPGAt vs ΔAUt-1. i Scatter plot of ΔPGAt vs ΔPFt. j Scatter plot of ΔPGAt vs ΔPFt-1

Short-term nonlinear relationships

Table 4 reports the Spearman’s rank correlation coefficients between changes in PGA and changes in the geomagnetic indices over short time periods (27 days). The Spearman’s rank correlation coefficient measures the strength of the monotonic relationship between two variables. The Spearman’s rank correlation coefficient is more general than Pearson’s correlation coefficient, as it also captures nonlinear relationships.

Table 4 Short-term (27 days) Spearman’s rank correlation coefficients of changes in PGA with the changes in the geomagnetic indices in the same period (first column) and of the changes in PGA in the current period with the changes in the geomagnetic indices in the previous period (second column)

As in the previous section, the first column of Table 4 reports the correlations of the changes in PGA with the changes in the geomagnetic indices in the same period (i.e. between the change in PGA over period t and the change in each geomagnetic index over the same period t), while the second column reports the correlations of the changes in PGA in the current period with the changes in the geomagnetic indices in the previous period (i.e. between the change in PGA over period t and the change in each geomagnetic index over the previous period t − 1).

The results in Table 4 indicate that, even after allowing for nonlinear monotonic relationships, over short time periods the association between changes in PGA and changes in the geomagnetic indices was still weak. None of the correlations in Table 4 is statistically significant after adjusting for multiple comparisons using Bonferroni correction (p > 0.1).

Figure 5 (supplementary section) shows the fitted Kernel regressions. Kernel regression is more general than linear regression as it allows for nonlinear relationships between the dependent variable and the independent variables. Furthermore, in Kernel regression the functional form of the nonlinear relationship is estimated nonparametrically, i.e. it does not need to be specified in advance.

Fig. 5
figure5

Short-Term Kernel Regressions (27-day differences). a Scatter plot of ΔPGAt vs ΔRt. b Scatter plot of ΔPGAt vs ΔRt-1. c Scatter plot of ΔPGAt vs ΔAPt. d Scatter plot of ΔPGAt vs ΔAPt-1. e Scatter plot of ΔPGAt vs ΔF10.7t. f Scatter plot of ΔPGAt vs ΔF10.7t-1. g Scatter plot of ΔPGAt vs ΔAUt. h Scatter plot of ΔPGAt vs ΔAUt-1. i Scatter plot of ΔPGAtvs ΔPFt. j Scatter plot of ΔPGAt vs ΔPFt-1.

Note that in this case there was no relationship between the R2 reported in Fig. 5 (supplementary section) and the correlation coefficients reported in Table 4. Note also that, as explained above, Spearman’s rank correlation coefficient measures only monotonic relationships, while as shown in Fig. 5 (supplementary section), Kernel regression also captures non-monotonic relationships.

Long-term nonlinear relationships

Table 5 reports the Spearman’s rank correlation coefficients between changes in PGA and changes in the geomagnetic indices over long time periods (297 days).

Table 5 Long-term (297 days) Spearman’s rank correlation coefficients of changes in PGA with the changes in the geomagnetic indices in the same period (first column) and of the changes in PGA in the current period with the changes in the geomagnetic indices in the previous period (second column)

The results are similar to the ones reported in Table 3:

  • Increases in R (Sunspot No.) are associated with decreases in the PGA, i.e. positive

    ΔRt are associated with negative (p < 0.05). However, this correlation is not statistically significant after adjusting for multiple comparisons using Bonferroni correction (p > 0.1). Furthermore, increases in R (Sunspot No.) anticipate decreases in the PGA, i.e. positive ΔRt − 1 tend to be associated with negative ΔPGAt (p < 0.01). This correlation is still statistically significant after adjusting for multiple comparisons using Bonferroni correction (p < 0.05).

  • Increases in the AP index were associated with decreases in the PGA, i.e. positive ΔAPtare associated with negative ΔPGAt(p < 0.01). This correlation is still statistically significant after adjusting for multiple comparisons using Bonferroni correction (p < 0.05).

  • Increases in the F10.7 index are associated with decreases in the PGA, i.e. positive ΔF10.7t are associated with negative (p < 0.01). Furthermore, increases in the F10.7 index anticipate decreases in the PGA, i.e. positive ΔF10.7t − 1 are associated with negative (p < 0.01). Both correlations are still statistically significant after adjusting for multiple comparisons using Bonferroni correction (p < 0.05).

Figure 6 (supplementary section) shows the fitted Kernel regressions. The interpretation of the R-squared is the same as in Section 2.3.

Fig. 6
figure6

Long-Term Kernel Regressions (297-day differences) a Scatter plot of ΔPGAt vs ΔRt. b Scatter plot of ΔPGAt vs ΔRt-1. c Scatter plot of ΔPGAt vs ΔAPt. d Scatter plot of ΔPGAt vs ΔAPt-1. e Scatter plot of ΔPGAt vs ΔF10.7t. f Scatter plot of ΔPGAt vs ΔF10.7t-1. g Scatter plot of ΔPGAt vs ΔAUt. h Scatter plot of ΔPGAtvs ΔAUt-1. i Scatter plot of ΔPGAt vs ΔPFt. j Scatter plot of ΔPGAt vs ΔPFt-1.

Fig. 7
figure7

Levels time series. a PGAt. b Rt. c APt. d F10.7t e AUt. f PFt

Fig. 8
figure8

27-day differences time series. a PGAt. b Rt. c APt. d F10.7t e AUt. f PFt

Fig. 9
figure9

297-day differences time series. a PGAt. b Rt. c APt. d F10.7t e AUt. f PFt

Discussion

Overall, this study demonstrates the presence of a seasonality cycle of 16 years as an important factor underlying the variance of the SLE global disease activity time series, and suggests an association of SLE global disease activity with geomagnetic activity, sunspot numbers, and high energy proton fluxes.

The time series decomposition of the 27-day PGA averages over a period of 24 years revealed that seasonality (one cycle every 16 years) explains a significant fraction of the variance of the time series. A previous study from the Hopkins Lupus Cohort showed no evidence of short-term (within 1 year) seasonality of global SLE activity expressed by PGA [10]. While many short-, mid-, and long-term periodicities in global climate, solar and geomagnetic activity have been described [20, 21], only the N-S (North-South) asymmetry periodicity which is 16.5 years for both small and large sunspots is similar to the periodicity described in our PGA time series. There are a growing number of examples where the immune system undergoes complex temporal oscillations, some well-known such as periodic fevers during malarial infection, familial Mediterranean fever, cyclic neutropenia, and other less known examples including oscillations in the levels of NFκB [22], HeS1 [23], and the tumor suppressor protein p53 [24]. Apart from these, the immunological literature includes a variety of isolated reports of oscillations [25]. Overall, however, there is still relatively little understanding of the causes and implications of oscillatory behavior of the immune system and its clinical implications.

Oscillatory behavior is an inherent characteristic of solar output and thus of geomagnetic activity. Yearly sunspot number analysis reveals periodicities, (in order of decreasing power) as 11-years, 38.6-years, 22.7-years [20], but analysis of their double-peaked behavior demonstrated dominant periods of 9, 12, and 16.5 years [26]. In geomagnetic data, the periodicities (in order of decreasing power) are 11-years, 5.3-years, 30-years and 46.5-years and some others [20]. The majority of great geomagnetic storms take place in the descending phase of the solar cycles [27]. When a solar cycle is strong, the great geomagnetic storms concentrate during the period from the 2 years before the solar cycle peak time to the 3 years after the solar cycle peak time [27]. In a database [4] of 129,205 deaths from myocardial infarction recorded in Minnesota from 1968 to 1996, the yearly incidences of myocardial infarction were categorized by solar cycle stage: maximum, descending stage, minimum, and ascending stage. An excess of 220 deaths from myocardial infarction per year was seen during years of maximal vs. minimal solar activity [4]. On the other hand, during periods of solar minima the penetration of cosmic rays in Earth’s atmosphere is higher which results in atmospheric air ionization, which has been associated with atmospheric electricity, cloudiness, and climate, all affecting human health [28]. For example, in a populational retrospective study from Sao Paulo, strong correlations between cosmic ray exposure and total mortality, infectious diseases, maternal mortality, and perinatal mortality rates were observed [29].

We show that short-term (27-day) increases in geomagnetic activity Ap index (p < 0.1) and high energy proton fluxes (> 60 Mev) (p < 0.05) were associated with decreases in SLE disease activity, while short-term increases in the sunspot number index R anticipated decreases in the SLE disease activity expressed as PGA (p < 0.05). The long-term (297 day) relationships showed a stronger negative correlations between changes in the PGA and changes in the sunspot number index R (p < 0.01), AP index (p < 0.01), and the F10.7 index (p < 0.01), which remained statistically significant after adjusting for multiple comparisons using Bonferroni correction.

We can only speculate on potential underlying mechanisms. The combined effects of sunspot-induced changes in solar irradiance and increases in atmospheric greenhouse gases have been shown to offer the best explanation for the rise in average global temperature over the last century [30], which is also a tempting hypotheses for its effect on SLE disease activity, as temperature changes have been associated with SLE organ specific flares [11].

Increased geomagnetic disturbances are likely to activate the photo/magneto-reception system, including the ferromagnetic receptors and the cryptochrome protein [31,32,33], affect cell membrane excitability by the enhancement of ionic motion, signaling, and accumulation in ion channels, with primary focus on intracellular calcium [34], and ultimately regulate melatonin secretion and 24-h circadian rhythm disruption [35]. These effects have been postulated to underlie the association of geomagnetic disturbances with influenza pandemics [36]. Another hypothesis includes the possibility that diurnal geomagnetic variation is a zeitgeber for biological circadian rhythms, and animals could perceive geomagnetic storms as a disruption of the usual circadian rhythm [37].Whether these mechanisms could affect SLE disease incidence or disease activity remains unknown, but warrants further study.

The shortcoming of our study is the inclusion of data from only one geographical area. The short time span of the data, covering only 2–3 solar cycles, may have contributed to the loss of statistical power. Future studies should analyze SLE data covering a longer interval from many locations and organize them in geomagnetic coordinates [38].

Availability of data and materials

The datasets during and/or analyzed during the current study available from the corresponding author on reasonable request.

References

  1. 1.

    Hanslmeier A. The sun and space weather; 1970. p. 233–49. https://doi.org/10.1007/978-3-642-11341-3_14.

    Book  Google Scholar 

  2. 2.

    Zhadin MN. Review of russian literature on biological action of DC and low-frequency AC magnetic fields. Bioelectromagnetics. 2001;22(1):27–45. https://doi.org/10.1002/1521-186x(200101)22:1<27::aid-bem4>3.0.co;2-2.

    CAS  Article  Google Scholar 

  3. 3.

    Lovetsky G, Kosushkin V, Samylov P, Aleksandrov M. Natural scientific foundations of Historiometry of Chizhevsky: Atlantis Press; 2019. p. 1116–9. https://doi.org/10.2991/cesses-19.2019.250.

  4. 4.

    Cornelissen-Guillaume GG, Halberg F, Breus T, Syutkina EV, Baevsky R, Weydahl A, Watanabe Y, Otsuka K, Siegelova J, Fiser B, Bakken EE. Non-photic solar associations of heart rate variability and myocardial infarction. J Atmospheric Sol-Terr Phys. 2002;64(5–6):707–20. https://doi.org/10.1016/S1364-6826(02)00032-9.

    Article  Google Scholar 

  5. 5.

    [The effect of geomagnetic disturbances in capillary blood flow in ischemic heart disease patients] - Abstract - Europe PMC. Accessed May 17, 2020. https://europepmc.org/article/med/7495904

  6. 6.

    Ghione S, Mezzasalma L, Del Seppia C, Papi F. Do geomagnetic disturbances of solar origin affect arterial blood pressure? J Hum Hypertens. 1998;12(11):749–54. https://doi.org/10.1038/sj.jhh.1000708.

    CAS  Article  Google Scholar 

  7. 7.

    O’Connor RP, Persinger MA. Geophysical variables and behavior: LXXXII. A strong association between sudden infant death syndrome and increments of global geomagnetic activity—possible support for the melatonin hypothesis. Percept Mot Skills. 1997;84(2):395–402. https://doi.org/10.2466/pms.1997.84.2.395.

    Article  Google Scholar 

  8. 8.

    Sajedi SA, Abdollahi F. Geomagnetic disturbances may be environmental risk factor for multiple sclerosis: an ecological study of 111 locations in 24 countries. BMC Neurol. 2012;12:100. https://doi.org/10.1186/1471-2377-12-100.

    Article  Google Scholar 

  9. 9.

    Wing S, Rider LG, Johnson JR, Miller FW, Matteson EL, Crowson CS, Gabriel SE. Do solar cycles influence giant cell arteritis and rheumatoid arthritis incidence? BMJ Open. 2015;5(5). https://doi.org/10.1136/bmjopen-2014-006636.

  10. 10.

    Duarte-García A, Fang H, To CH, Magder LS, Petri M. Seasonal variation in the activity of systemic lupus erythematosus. J Rheumatol. 2012;39(7):1392–8. https://doi.org/10.3899/jrheum.111196.

    Article  Google Scholar 

  11. 11.

    Stojan G, Kvit A, Curriero FC, Petri M. A spatiotemporal analysis of organ-specific lupus flares in relation to atmospheric variables and fine particulate matter pollution. Arthritis Rheumatol Hoboken NJ Published online February 3. 2020. https://doi.org/10.1002/art.41217.

  12. 12.

    Stoupel E, Monselise Y, Lahav J. Changes in autoimmune markers of the anti-cardiolipin syndrome on days of extreme geomamagnetic activity. J Basic Clin Physiol Pharmacol. 2006;17(4):269–78. https://doi.org/10.1515/jbcpp.2006.17.4.269.

    CAS  Article  Google Scholar 

  13. 13.

    Hopkins PM, Cohort L. 1999 update. Rheum Dis Clin N Am. 2000;26(2):199–213 v.

    Article  Google Scholar 

  14. 14.

    Hochberg MC. Updating the American College of Rheumatology revised criteria for the classification of systemic lupus erythematosus. Arthritis Rheum. 1997;40(9):1725. https://doi.org/10.1002/1529-0131(199709)40:9&lt;1725:AID-ART29&gt;3.0.CO;2-Y.

    CAS  Article  Google Scholar 

  15. 15.

    Tan EM, Cohen AS, Fries JF, Masi AT, McShane DJ, Rothfield NF, Schaller JG, Talal N, Winchester RJ. The 1982 revised criteria for the classification of systemic lupus erythematosus. Arthritis Rheum. 1982;25(11):1271–7.

    CAS  Article  Google Scholar 

  16. 16.

    Petri M, Orbai A-M, Alarcón GS, Gordon C, Merrill JT, Fortin PR, Bruce IN, Isenberg D, Wallace DJ, Nived O, Sturfelt G, Ramsey-Goldman R, Bae S-C, Hanly JG, Sanchez-Guerrero J, Clarke A, Aranow C, Manzi S, Urowitz M, Gladman D, Kalunian K, Costner M, Werth VP, Zoma A, Bernatsky S, Ruiz-Irastorza G, Khamashta MA, Jacobsen S, Buyon JP, Maddison P, Dooley MA, van Vollenhoven RF, Ginzler E, Stoll T, Peschken C, Jorizzo JL, Callen JP, Lim SS, Fessler BJ, Inanc M, Kamen DL, Rahman A, Steinsson K, Franks AG, Sigler L, Hameed S, Fang H, Pham N, Brey R, Weisman MH, McGwin G, Magder LS. Derivation and Validation of Systemic Lupus International Collaborating Clinics Classification Criteria for Systemic Lupus Erythematosus. Arthritis Rheum. 2012;64(8):2677–86. https://doi.org/10.1002/art.34473.

    Article  Google Scholar 

  17. 17.

    Petri M, Genovese M, Engle E, Hochberg M. Definition, incidence, and clinical description of flare in systemic lupus erythematosus. A prospective cohort study. Arthritis Rheum. 1991;34(8):937–44. https://doi.org/10.1002/art.1780340802.

    CAS  Article  Google Scholar 

  18. 18.

    King JH, Papitashvili NE. Solar wind spatial scales in and comparisons of hourly Wind and ACE plasma and magnetic field data. J Geophys Res Space Physics. 2005;110(A2). https://doi.org/10.1029/2004JA010649.

  19. 19.

    Fourier Analysis of Time Series: An Introduction, 2nd Edition | Wiley. Accessed May 19, 2020. https://www.wiley.com/en-us/Fourier+Analysis+of+Time+Series%3A+An+Introduction%2C+2nd+Edition-p-9780471889489

  20. 20.

    Singh YP. Badruddin. Prominent short-, mid-, and long-term periodicities in solar and geomagnetic activity: wavelet analysis. Planet Space Sci. 2014;96:120–4. https://doi.org/10.1016/j.pss.2014.03.019.

    Article  Google Scholar 

  21. 21.

    Li Z, Yue J, Xiang Y, Chen J, Bian Y, Chen H. Multiresolution analysis of the relationship of solar activity, global temperatures, and global warming. Adv Meteorol. https://doi.org/10.1155/2018/2078057.

  22. 22.

    Hoffmann A, Levchenko A, Scott ML, Baltimore D. The IkappaB-NF-kappaB signaling module: temporal control and selective gene activation. Science. 2002;298(5596):1241–5. https://doi.org/10.1126/science.1071914.

    CAS  Article  Google Scholar 

  23. 23.

    Hirata H, Yoshiura S, Ohtsuka T, Bessho Y, Harada T, Yoshikawa K, Kageyama R. Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop. Science. 2002;298(5594):840–3. https://doi.org/10.1126/science.1074560.

    CAS  Article  Google Scholar 

  24. 24.

    Lev Bar-Or R, Maya R, Segel LA, Alon U, Levine AJ, Oren M. Generation of oscillations by the p53-Mdm2 feedback loop: a theoretical and experimental study. Proc Natl Acad Sci U S A. 2000;97(21):11250–5. https://doi.org/10.1073/pnas.210171597.

    CAS  Article  Google Scholar 

  25. 25.

    Stark J, Chan C, George AJT. Oscillations in the immune system. Immunol Rev. 2007;216(1):213–31. https://doi.org/10.1111/j.1600-065X.2007.00501.x.

    Article  Google Scholar 

  26. 26.

    Mandal S, Banerjee D. Sunspot sizes and the solar cycle: analysis using kodaikanal white-light digitized data. Astrophys J. 2016;830(2):L33. https://doi.org/10.3847/2041-8205/830/2/L33.

    Article  Google Scholar 

  27. 27.

    Le G-M, Cai Z, Wang H, Zhu Y. Solar cycle distribution of great geomagnetic storms. Astrophys Space Sci. 2011;339. https://doi.org/10.1007/s10509-011-0960-y.

  28. 28.

    Kirkby J. Cosmic rays and climate. Surv Geophys. 2007;28(5):333–75. https://doi.org/10.1007/s10712-008-9030-6.

    Article  Google Scholar 

  29. 29.

    Vieira CLZ, Janot-Pacheco E, Lage C, Pacini A, Koutrakis P, Cury PR, Shaodan H, Pereira LA, Saldiva PHN. Long-term association between the intensity of cosmic rays and mortality rates in the city of Sao Paulo. Environ Res Lett. 2018;13(2):024009. https://doi.org/10.1088/1748-9326/aaa27a.

    CAS  Article  Google Scholar 

  30. 30.

    Lane LJ, Nichols MH, Osborn HB. Time series analyses of global change data. Environ Pollut. 1994;83(1):63–8. https://doi.org/10.1016/0269-7491(94)90023-X.

    CAS  Article  Google Scholar 

  31. 31.

    Knox EG, Armstrong E, Lancashire R, Wall M, Haynes R. Heart attacks and geomagnetic activity. Nature. 1979;281(5732):564–5. https://doi.org/10.1038/281564a0.

    CAS  Article  Google Scholar 

  32. 32.

    Stoupel E, Kalediene R, Petrauskiene J, Starkuviene S, Abramson E, Israelevich P, Sulkes J. Twenty years study of solar, geomagnetic, cosmic ray activity links with monthly deaths number (n-850304). J Biomed Sci Eng. 2011;4(6):426–34. https://doi.org/10.4236/jbise.2011.46054.

    Article  Google Scholar 

  33. 33.

    Palmer SJ, Rycroft MJ, Cermack M. Solar and geomagnetic activity, extremely low frequency magnetic and electric fields and human health at the Earth’s surface. Surv Geophys. 2006;27(5):557–95. https://doi.org/10.1007/s10712-006-9010-7.

    Article  Google Scholar 

  34. 34.

    Zaporozhan V, Ponomarenko A. Mechanisms of geomagnetic field influence on gene expression using influenza as a model system: basics of physical epidemiology. Int J Environ Res Public Health. 2010;7(3):938–65. https://doi.org/10.3390/ijerph7030938.

    CAS  Article  Google Scholar 

  35. 35.

    Burch JB, Reif JS, Yost MG. Geomagnetic disturbances are associated with reduced nocturnal excretion of a melatonin metabolite in humans. Neurosci Lett. 1999;266(3):209–12. https://doi.org/10.1016/s0304-3940(99)00308-0.

    CAS  Article  Google Scholar 

  36. 36.

    Qu J. Is sunspot activity a factor in influenza pandemics? Rev Med Virol. 2016;26(5):309–13. https://doi.org/10.1002/rmv.1887.

    Article  Google Scholar 

  37. 37.

    Krylov VV. Biological effects related to geomagnetic activity and possible mechanisms. Bioelectromagnetics. 2017;38(7):497–510. https://doi.org/10.1002/bem.22062.

    Article  Google Scholar 

  38. 38.

    Baker KB, Wing S. A new magnetic coordinate system for conjugate studies at high latitudes. J Geophys Res Space Physics. 1989;94(A7):9139–43. https://doi.org/10.1029/JA094iA07p09139.

    Article  Google Scholar 

  39. 39.

    Rostoker G. Geomagnetic indices. Rev Geophys. 1972;10(4):935–50. https://doi.org/10.1029/RG010i004p00935.

    Article  Google Scholar 

  40. 40.

    Solar-Terrestrial Indices. Accessed May 17, 2020. https://ccmc.gsfc.nasa.gov/modelweb/solar/solar_index.html

  41. 41.

    Rodriguez JV, Sandberg I, Mewaldt RA, Daglis IA, Jiggens P. Validation of the effect of cross-calibrated GOES solar proton effective energies on derived integral fluxes by comparison with STEREO observations. Space Weather. 2017;15(2):290–309. https://doi.org/10.1002/2016SW001533.

    Article  Google Scholar 

Download references

Funding

Dr. George Stojan is supported by the Office of the Assistant Secretary of Defense for Health Affairs, through the Lupus Research Program under Award No. W81XWH19–1- 0793. Dr. George Stojan is a Jerome L. Greene Foundation Scholar. The Hopkins Lupus Cohort is supported by a grant from the National Institute of Health (NIH AR 43727 and 69572). This publication was also made possible by Grant Number UL1 RR 025005 from the National Center for Research Resources (NCRR), a component of the National Institutes of Health (NIH), and NIH Roadmap for Medical Research.

Author information

Affiliations

Authors

Contributions

G.S. contributed to the conception, design of the study, acquisition, analysis, and interpretation of data; and drafted the manuscript. F.G. contributed to the design of the study; acquisition, analysis, and interpretation of data, and revised the manuscript. M.P. contributed to the acquisition and interpretation of data and revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to George Stojan.

Ethics declarations

Ethics approval and consent to participate

The Hopkins Lupus Cohort has been approved by the Johns Hopkins University School of Medicine Institutional Review Board and complies with the Health Insurance Portability and Accountability Act. All patients gave written informed consent.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1. Supplementary section

The Ap index is a measure of the general level of geomagnetic activity over the globe for a given day. It is derived from measurements made at a number of stations world-wide of the variation of the geomagnetic field due to currents flowing in the earth's ionosphere and, to a lesser extent, in the earth's magnetosphere. The official values for Ap are calculated by the GeoForschungsZentrum Helmholtz Centre Potsdam (Germany) [39]. The sunspot number index R is a measure of the area of solar surface covered by spots. As the number of spots increases and their magnetic complexity grows, they become likely sources of large eruptive energy releases known as solar flares. The sunspot number index is also often called Wolf number in reference to the Swiss astronomer J. R. Wolf who introduced this index in 1848 [40]. The sun emits radio energy with slowly varying intensity. This radio flux, which originates from atmospheric layers high in the sun's chromosphere and low in its corona, changes gradually from day to day in response to the number of spot groups on the disk. Solar flux from the entire solar disk at a frequency of 2800 MHz has been recorded routinely by a radio telescope near Ottawa since February 1947 and is called the F10.7 index40. The AU (amplitude upper) index describes the disturbance level recorded by auroral zone magnetometers [40]. The radiation hazard from solar proton events (SPEs) has been characterized in terms of integral fluxes above selected threshold energies. We used high energy proton fluxes >60 MeV for our analysis [41].

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Stojan, G., Giammarino, F. & Petri, M. Systemic lupus Erythematosus and geomagnetic disturbances: a time series analysis. Environ Health 20, 28 (2021). https://doi.org/10.1186/s12940-021-00692-4

Download citation

Keywords

  • Systemic lupus erythematosus
  • Lupus
  • Environment
  • Geomagnetic disturbance
  • Time series
  • Disease activity