Analysis of cardiac signals using spatial filling index and timefrequency domain
 Oliver Faust^{1}Email author,
 Rajendra Acharya U^{1},
 SM Krishnan^{2} and
 Lim Choo Min^{1}
https://doi.org/10.1186/1475925X330
© Faust et al; licensee BioMed Central Ltd. 2004
Received: 14 May 2004
Accepted: 10 September 2004
Published: 10 September 2004
Abstract
Background
Analysis of heart rate variation (HRV) has become a popular noninvasive tool for assessing the activities of the autonomic nervous system (ANS). HRV analysis is based on the concept that fast fluctuations may specifically reflect changes of sympathetic and vagal activity. It shows that the structure generating the signal is not simply linear, but also involves nonlinear contributions. These signals are essentially nonstationary; may contain indicators of current disease, or even warnings about impending diseases. The indicators may be present at all times or may occur at random in the time scale. However, to study and pinpoint abnormalities in voluminous data collected over several hours is strenuous and time consuming.
Methods
This paper presents the spatial filling index and timefrequency analysis of heart rate variability signal for disease identification. Renyi's entropy is evaluated for the signal in the WignerVille and Continuous Wavelet Transformation (CWT) domain.
Results
This Renyi's entropy gives lower 'p' value for scalogram than WignerVille distribution and also, the contours of scalogram visually show the features of the diseases. And in the timefrequency analysis, the Renyi's entropy gives better result for scalogram than the WignerVille distribution.
Conclusion
Spatial filling index and Renyi's entropy has distinct regions for various diseases with an accuracy of more than 95%.
Keywords
Background
Biosignals are essentially nonstationary signals; they display a fractal like selfsimilarity. They may contain indicators of current disease, or even warnings about impending diseases. The indicators may be present at all times or may occur at random – in the time scale. However, to (study and) pinpoint anomalies in voluminous data collected over several hours is strenuous and time consuming. Therefore, computer based analytical tools for indepth study and classification of data over day long intervals can be very useful in diagnostics.
Electrocardiography deals with the electrical activity of the heart. Monitored by placing sensors at defined positions on chest and limb extremities of the subject, electrocardiogram (ECG) is a record of the origin and propagation of the electric action potential through cardiac muscle. It is considered a representative signal of cardiac physiology, useful in diagnosing cardiac disorders. The state of cardiac health is generally reflected in the shape of ECG waveform and heart rate. It may contain important pointers to the nature of diseases afflicting the heart. However, biosignals being nonstationary signals, this reflection may occur at random in the time scale.
(That is, the disease symptoms may not show up all the time, but would manifest at certain irregular intervals during the day.) Therefore, for effective diagnostics, the study of ECG pattern and heart rate variability signal (instantaneous heart rate against time axis) may have to be carried out over several hours. HRV is a useful signal for understanding the status of the autonomic nervous system (ANS).
The interest in the analysis of heart rate variability (HRV), that is, the fluctuations of the heart beating in time, is not new. And much progress was achieved in this field with the advent of cheap and massive computational power, which provoked many recent advances.
HRV is a noninvasive measurement of cardiovascular autonomic regulation. Specifically, HRV is a measurement of the interaction between sympathetic and parasympathetic activity in autonomic functioning. There are two main approaches for analysis: time domain analysis of HRV [for standard deviation of normal to normal intervals (SDNN)]; and frequency domain analysis [for power spectrum density (PSD)]. The latter provides high frequency (parasympathetic activity) and low frequency (sympathetic and parasympathetic activity) and total power (sympathetic/parasympathetic balance) values [1–3]. Recent results on HRV signal analysis show that its dynamic behavior involves nonlinear components that also contribute in the signal generation and control [4]. The autonomic nervous system (ANS) modulates the cardiac pacemaker and provides beattobeat regulation of the cardiovascular rhythm. Application of wavelet transformation techniques to beattobeat heart rate variations (HRV) provides an important noninvasive tool for monitoring the autonomic nervous system functioning.
The cardiovascular system is a complex system that includes heart and vessels. ECG and HRV are two methods for study it. Hence, many attempts have been made to analyze these signals and extract information about the cardiovascular system. Most of the methods used are linear and it has been recognized that nonlinear methods may be more suitable for analyzing signals that originate from complex nonlinear living systems [5]. Recent developments in nonlinear analysis have provided various methods for the study of the complex cardiovascular system [6]. It is now generally recognized that many processes generated by the biological system can be described in an effective way by using the methods of nonlinear dynamics. The nonlinear dynamical techniques are based on the concept of chaos, which was first introduced with applications to complicated dynamical systems in meteorology [7]. Since then, it has been applied to medicine and biology [8, 9]. A particularly active area for the application of chaos theory has been cardiology [10, 11], where many aspects have been addressed including whether chaos can be used to represent healthy or diseased state [12].
A complex system like cardiovascular system can not be linear in nature and by considering it as a nonlinear system can lead to better understanding of the system dynamics. Recent studies have also stressed the importance of nonlinear techniques to study HRV in both health and disease. The progress made in the field using measures of chaos has attracted scientific community applying these tools in studying physiological systems, and HRV is no exception. There have been several methods of estimating invariants from nonlinear dynamical systems reported in the literature. Recently, Fell et al and Radhakrishna et al have tried the nonlinear analysis of ECG and HRV signals respectively [13, 14]. Also, Addison at al showed that coordinated mechanical activity in the heart during ventricular fibrillation may be made visible in the surface ECG using wavelet transform [15]. Rajendra et al, [16] have classified the HRV signals using Artificial Neural Networks (ANN) and Fuzzy equivalence relation. Recently, Renyi's entropy is used for texture analysis by Grigorescu et al [17]. Gokcay et al have applied Renyi's entropy to clustering and analyze the resulting staircase nature of the performance function that can be expected during learning [18]. In this work, different heart rate signals are analyzed using spatial filling index and time frequency techniques. Renyi's entropy is evaluated for the different cardiac abnormalities.
Methods
Number of subjects in various groups
Type  Normal  PVC  CHB  SSS  CHF  ISCDIL  AF  VF 

Number of datasets  60  60  20  20  40  20  35  45 
Each dataset is taken consists of more than10,000 samples and the sampling frequency of the data is 360 Hz. The interval between two successive QRS complexes is defined as the RR interval (t_{rr} seconds) and the heart rate (beats per minute) is given as:
HR = 60/t_{rr} (1)
Spatial Filling Index
Let the signal be represented by the coordinates of a point X(k) in phase space. Then the dynamical behavior of the signal is reconstructed by succession of these points X(k) in the phase space. Phase space reconstructions are based on the analysis of dynamic systems by delay maps. The vectors X(k) in the multidimensional phase space are constructed by time delayed values of the time series, which determine the coordinates of the phase space plot.
X(k) = {x (k), x((k + τ), ...,x (k + (E1)τ)} for k = 1,2,...,N  (E  1)τ (1)
where X(k) is one point of the trajectory in the phase space at time k, x(k + τ) are the coordinates in the phase space corresponding to the time delayed values of the time series, τ is the time delay between the points of the time series considered and E is the embedding dimension, which is the number of coordinates of the phase space plot. The attributes of the reconstructed phase space plot depend on the choice of value of τ. One way to choose τ is to take it as the time it takes the autocorrelation function of the data to decay to 1/e [21]. Another method is to take the first minimum in the graph of average mutual information [22], which appears to be better since it considers the nonlinear structure in the signal. It has been established using this method that the value of 7 for τ is the best choice for ECG signals and 5 for HRV signals [23].
From the given signal x (1), x(2), ..., x (N), a matrix A _{ E }is obtained as
where E is the number of dimensions and M is related to N by the equation:
M = N  (E 1)τ (3)
By plotting column 2 of matrix A against column 1 (for the case E = 2), the phase space plot for two dimensions is obtained.
Similarly, the first three columns of matrix A _{3} represent a phase space plot in three dimensions. Now, a normalized matrix B _{ E }is obtained by dividing each element of A _{ E }by x _{max} where
x _{max} = max x(k) 1 ≤ k ≤ N (5)
The matrix B _{2} (in two dimensions) is hence represented as
In two dimensions, the phase space plot corresponding to the normalized matrix spans from 1 to +1 on either axis. The phase space area is now divided into small square areas of size {R × R R ∈ Real, 2/R ∈ Integer}. Then the number of grids in the normalized phase space is n = 2/R. A matrix C is now obtained with its elements c(i,j) equal to the number of phase space points falling in a grid g(i,j). The matrix C is called the phase space matrix and its elements are divided by m, where
This division yields P(i,j), the probability of a phase space point falling in a grid g(i,j). A matrix Q is now formed by squaring each element of P to get q(i,j) as the elements of Q. The sum of elements of matrix Q is calculated as
The spatial filling index η is defined as:
η = s / n^{2} (9)
Now the value of η is used to quantify the degree of variability in the test signals.
TimeFrequency analysis
There are three common approaches to generating the timefrequency (TF) plots. These are the short Time Fourier Transform; the WignerVille based bilinear distributions and the Continuous Wavelet Transform. In this investigation the latter two were used.
WignerVille analysis
The WignerVille distribution (WVD) is defined as:
where z(t) is the analytic signal and h(τ) is a window function. The results where obtained using a Hamming window. This window attenuates the interferences oscillating perpendicularly to the frequency axis. The WVD satisfies a large number of desirable mathematical properties. In particular, the WVD is always realvalued; it preserves time and frequency shifts and satisfies the marginal properties. Moreover, the WVD conserves the Energy of the signal. We obtain the Energy (E _{ x }) by integrating the WVD of z all over the time frequency plane:
With the Energy conservation property the WVD can be interpreted in terms of probability density: expression (10) is the Fourier transform of an acceptable form of characteristic function for the distribution of energy. Therefore, the WVD can be used to obtain the information content of a signal; this thought is further extended in Section 4.3.
Continuous Time Wavelet Transform (CWT) analysis
A 'wavelet' implies a small wave of finite duration and finite energy, which is correlated with the signal to obtain the wavelet coefficients [24]. The reference wavelet is known as the mother wavelet, and the coefficients are evaluated for the entire range of dilation and translation factors [25]. Initially the mother wavelet is shifted (translated) continually along the time scale for evaluating the set of coefficients at all instants of time. In the next phase, the wavelet is dilated for a different width – also normalized to contain the same amount of energy as the mother wavelet – and the process is repeated for the entire signal. The wavelet coefficients are real numbers usually shown by the intensity of a chosen color, against a two dimensional plane with yaxis representing the dilation (scaling factor) of the wavelet, and the xaxis, its translation (shift along the time axis). Thus the wavelet transform plot (scalogram) can be seen as a color pattern against a two dimensional plane. In the CWT the wavelet coefficients are evaluated for infinitesimally small shifts of translation as well as scale factors. That is, the color intensity of each pixel in the scalogram is separately evaluated, and the resulting pattern contains information about the size and location of the 'event' occurring in the time domain [26, 27]. Since the dilated wavelet is normalized to contain the same amount of energy as the mother wavelet; the scalogram representation of even high frequency, low energy 'events' occurring in the time scale are more conspicuous than in the Fourier Transform. Thus the color patterns in the scalogram can be useful in highlighting the abnormalities specific to different types of disease. MATLAB version 6.1 is used to plot the various scalogram plots.
For a given wavelet Ψ _{ a,b }(t), the coefficients are evaluated using Eq. (12):
Just like the WVD, the CWT representation preserves also the energy of the signal. The total energy (E _{ x }) is obtained by integrating over all scale and translation factors:
The scalogram patterns thus obtained also depend on the wavelet chosen for analysis. Biosignals usually exhibit self similarity patterns in their distribution, and a wavelet which is akin to its fractal shape would yield the best results in terms of clarity and distinction of patterns. In the present work, the analysis is based on the Morlet wavelet. This wavelet gives good result compared to all the other wavelets.
The Morlet wavelet function is given by:
Renyi's Entropy (RE)
The previous sections detailed WVD and CWT as two methods to represent a signal in the timefrequency domain. This section is concerned with the interpretation of the timefrequency representation. The signals represent measurements taken from patients being either normal or suffering from different vascular diseases. The goal is to find a measure which allows classifying the different signals according to the medical conditions.
One interesting information that one may obtain from the timefrequency representation is the number of elementary signals present in the current observation. This leads to the following question: How much separation between two elementary signals must one achieve in order to be able to conclude that there are two signals present rather then one?
A solution to this problem is given by applying an information measure to a timefrequency distribution of a signal. This can be done, because CWT and WVD preserve the energy of the signal.
Unfortunately, the well known Shannon information can not be applied to the timefrequency representation of a signal, because it contains negative values. One information measure, which allows negative values in the distribution, is Renyi's entropy. This information measure was used to analyze the timefrequency representation of the measurement data.
Renyi's entropy definition is derived from his proposed theory of means [28]
where φ(.) is a continuous and strictly monotonic function subclass of KolmogorovNagumo functions. To satisfy the constraints of an information measure
I(p _{ k }) any information measure
Simplifying the above relation, we have
The third order Renyi's entropy (α = 3) is calculated from the WVD timefrequency representations as follows:
similarly, the third order Renyi's entropy is calculated form the CWT as follows:
OneWay Analysis of Variance (ANOVA)
The purpose of oneway ANOVA is to find out whether data from several groups have a common mean. That is, to determine whether the groups are actually different in the measured characteristic.
Oneway ANOVA is a simple special case of the linear model. The oneway ANOVA form of the model is where:
y _{ ij }= α _{.j }+ ε _{ ij }

y _{ ij }is a matrix of observations in which each column represents a different group.

α _{.j }is a matrix whose columns are the group means. (The "dot j" notation means that applies to all rows of the jth column. That is, the value α _{ ij }is the same for all i.)

ε _{ ij }is a matrix of random disturbances.
The model posits that the columns of y are a constant plus a random disturbance. You want to know if the constants are all the same.
Results
Results for various cardiac abnormalities
Type  SSS  PVC  CHB  NORMAL  CHF  AF  ISCDIL  VF  pvalue 

η Phase Space  1.56 ± 0.08  3.77 ± 11.82  7.07 ± 0.57  6.76 ± 1.61  7.71 ± 0.20  2.26 ± 0.05  7.44 ± 1.24  3.77 ± 8.78  0.00005 
3.38 ± 2.74  4.79 ± 1.46  5.65 ± 0.61  4.00 ± 1.52  2.10 ± 0.61  3.31 ± 2.45  4.07 ± 2.78  4.44 ± 0.96  0.021  
2.84 ± 1.89  2.25 ± 1.06  3.01 ± 0.29  1.67 ± 0.84  1.78 ± 0.82  2.04 ± 1.37  2.15 ± 0.57  3.29 ± 0.33  0.001 
For the time frequency plots the normalized frequency is shown over the hart rate values. It is not useful to state an absolute frequency, because such a value is not relevant for the cardiac system under observation. Moreover, the relative frequency representation allows comparing the time frequency analysis results form varying observation intervals. As example, the observation interval for the VF data is significantly shorter as for the rest of the data, but still the results can be compared.
Discussion
The resulting phase space plots for various types of disease are shown in Figure 1(a),2(a),3(a),4(a),5(a),6(a),7(a),8(a). In SSS – III (Sick Sinus Syndrome – III, BradycardiaTachycardia) there is a continuous variation of heart rate between Bradycardia and Tachycardia. The phase space plot spreads over a larger area (Figure 1(a)). In the Ectopic beat abnormality; there would be a sudden impulsive jump in the heart rate. This may be due to a PrematureVentricular beat in the ECG signal. This is indicated as a sudden spike in the phase space plot (Figure 2(a)). In Complete Heart Block (CHB) cases, as the atrioventricular node fails to send electrical signals rhythmically to the ventricles, the heart rate remains low. The phase space plot reduces almost to a point, indicating very little change with time (Figure 3(a)). For Normal cases, the phase space plot looks like a cluster of points (Figure 4(a)). In the Congestive heart failure (CHF), the heart rate variation is lower and hence the phase space plot spread in a very small area (Figure 5(a)). In the Atrial Fibrillation (AF), heart rate signal records highly erratic variability; this is depicted as scattering of points in the phase space plot (Figure 6(a)). In the case of Ischemic/Dilated cardiomyopathy, the ventricles are unable to pump out blood to the normal degree. Here the heart rate variation is very small. And hence the phase space plot will be almost a point (Figure 7(a)). And its phase space plot resembles that of Normal class. Finally, in VF, the heart rate variation is high and hence the phase space plot is randomly distributed (Figure 8(a)).
The contour plots of scalogram and WignerVille distribution plot for the different abnormalities are shown in figures 1(b),1(c),2(b),2(c),3(b),3(c),4(b),4(c),5(b),5(c),6(b),6(c),7(b),7(c),8(b),8(c) respectively. In the contour plot of scalogram (Figure 1(b)), for SSS, there is clear indication of variation of high frequency and low frequency in the form of irregular circles at these frequencies. In PVC (Figure 2(b)), a irregular circle is shown at high frequency indicating the spike of the signal. These irregular circles or contours are at low frequencies for CHB (Figure 3(b)), indicating smaller RR variation. In normal case (Figure 4(b)) these contours are in the middle frequency due to variation in the RR interval. In CHF (Figure 5(b)) and Ischemic/Dialted cardiomyopathy (Figure 7(b)), the RR variation is extremely low. Hence the contours are aligned at the low frequency. In AF (Figure 6(b)), due to very high RR variation are shown as irregular contours at high frequency. For VF, this RR variation is slightly low and as result the contours are aligned at the middle of the contour plot (Figure 8(b)).
The contour plots of the WignerVille distribution does not indicate as clearly as contour plot of scalogram for various cardiac diseases.
The spatial filling index decreases or increases from the normal class for the abnormal subjects in different ranges (Table 2) depending on the RR variation. This value decreases for the abnormalities of high RR variation and increases for CHB, CHF and Ishemic/Dilated cardiomyoapathy, which has low RR variation. This parameter has excellent 'p' value for various classes (0.00005). The Renyi's entropy has high value for cardiac abnormalities like Ischemic/Dilated cardiomyopathy, CHB, VF and it decreases for Normal, PVC, AF, SSS and CHF. This RE gives lower 'p' value for scalogram than WignerVille distribution and also, the contours of scalogram visually show the features of the diseases. Hence, in the timefrequency analysis, the Renyi's entropy gives better result for scalogram than the WignerVille distribution.
Conclusion
Considering heart as a nonlinear complex system and processing various cardiovascular signals (HRV) seems to provide very useful information for detection of abnormalities in the condition of the heart that is not available by conventional means. In this paper, a phase space and the timefrequency analysis of these cardiac signals using spatial filling index and Renyi's entropy has been proposed for detecting cardiac dysfunction. The ANOVA test was used to compare the different analyzing methods. The Renyi's entropy gives better result for the scalogram than the WignerVille distribution. The evaluation of the proposed technique on a larger data set will improve the efficacy of the technique. It is left as future work to compare the different methods with more sophisticated statistical methods, such as post hoc comparisons. It is hoped that the graphical representation along with its corresponding analytical index and Renyi's entropy proposed here will find potential applications in computer analysis of cardiac patients' status in intensive care units.
Declarations
Authors’ Affiliations
References
 Myers G, Martin GJ, Magid NM, Barnett PS, Schaad JW, Weiss JS, Lesch M, Singer DH: Power Spectral Analysis of Sudden Cardiac Death: Comparison to Other Methods. IEEE Transactions on Biomedical Engineering 1986, 1149–1156.Google Scholar
 Akselrod S, Godo D, Ubel FA, Shannon DC, Baeger AC, Cohen RJ: Power spectrum analysis of heart rate fluctuation: a quantitative probe of beattobeat cardiovascular control. Science 1981, 220–222.Google Scholar
 Saul P, Yukata AL, Berger RD, Leonard SL, Wilson CS, Cohen RJ: Assessment of Autonomic Regulation in Chronic Congestive Heart Failure by Heart Rate Spectral Analysis. American Journal of Cardiology 1988, 1292–1299. 10.1016/00029149(88)911721Google Scholar
 Iyengar N, Peng CK, Goldberger AL, Lipsitz LA: Agerelated alterations in the fractal scaling of cardiac interbeat interval dynamics. American journal of Physiology 1996, 271: R1078–1084.Google Scholar
 Abarbanel HI: Analysis of Observed Chaotic Data Springer Verlag, Berlin 1996.View ArticleGoogle Scholar
 Karrakchou M, Rheymer KV, Vesin JM, Pruvot E, Kunt M: Improving cardiovascular monitoring through modern techniques. IEEE Engineering in Medicine and Biology Magazine 1996, 68–78. 10.1109/51.537062Google Scholar
 May RM: Simple mathematical models with very complicated dynamics. Nature 1976, 261: 459–467.View ArticleGoogle Scholar
 Tsonis PA, Tsonis A: Chaos: Principles and implications in biology. Computer Applications in Biosciences 1989, 5: 27–32.Google Scholar
 Eberhart RC: Chaos theory for biomedical engineer. IEEE Engineering in Medicine and Biology Magazine 1989, 41–45. 10.1109/51.35577Google Scholar
 Glass L, Zeng WZ: Complex bifurcations and chaos in simple theoretical models of cardiac oscillations. Annals of the New York Academy of Sciences 1990, 316–327.Google Scholar
 Denton TA, Diamond GA: Fascinating rhythms: a primer on chaos theory and its applications to cardiology. American Heart Journal 1990, 1419–1440. 10.1016/00028703(90)90258YGoogle Scholar
 Lipsitz LA, Goldberger A: Loss of complexity and aging. Journal of the American Medical Association 1992, 1806–1809. 10.1001/jama.267.13.1806Google Scholar
 Fell J, Mann K, Roschke J, Gopinath MS: Nonlinear analysis of continuous ECG during sleep I reconstruction. Biological Cybernetics 2000, 477–483. 10.1007/s004220050600Google Scholar
 Radhakrishna RA, Vikram Kumar Yergani , Narayana Dutt D, Vedavathy TS: Characterizing Chaos in heart rate variability time series of panic disorder patients. In Proceedings of ICBME, Biovision Bangalore, India 2001, 163–167.Google Scholar
 Addison PS, Watson JN, Clegg Gareth R, Steen PA, Robertson CE: Finding Coordinated Atrial Activity During Ventricular Fibrillation Using Wavelet Decomposition. IEEE Engineering In Medicine and Biology Magazine 2002, 58–61. 10.1109/51.993194Google Scholar
 Rajendra Acharya U, Pubbanna Bhat SS, Iyengar , Ashok Rao , Summet Dua : Classification of heart rate using artificial neural network and fuzzy equivalence relation. Pattern Recognition 2003, 61–68. 10.1016/S00313203(02)000638Google Scholar
 Grigorescu SE, Petkov N: Texture analysis using Renyi's generalized entropies. In IEEE International Conference on Image Processing 2003, 241–244.Google Scholar
 Gokcay E, Principe JC: A new clustering evaluation function using Renyi's information potential. In International Conference on Acoustics, Speech, and Signal Processing (ICASSP) 2000.Google Scholar
 Ichimaru Y, Moody GB: Development of the polysomnographic database on CDROM. Psychiatry and Clinical Neurosciences 1999, 175–177. 10.1046/j.14401819.1999.00527.xGoogle Scholar
 Jiapu P, Tompkins WJ: Real Time QRS Detector algorithm. IEEE Transactions on Biomedical Engineering 1985, 230–23.Google Scholar
 Bangham JA, Marshall S: Image and signal processing with mathematical morphology. Electronics & Communication Engineering Journal 1998, 117–128.Google Scholar
 Mallat S: Characterization of signals from multiscale edges. IEEE Transactions on Pattern Analysis and Machine Intelligence 1992, 1019–1033.Google Scholar
 Falconer K: Fractal Geometry J Wiley and Sons, New York 1990.Google Scholar
 Vetterli M, Herley C: Wavelets and filter banks: theory and design. IEEE Transactions on Signal Processing 1992, 2207–2232. 10.1109/78.157221Google Scholar
 Raghuveer MR, Bopardikar AS: Wavelet Transforms Introduction to theory and applications Addison Wesley Longman Inc 1998.Google Scholar
 Vetterli M, Kovacevic J: Wavelets and Subband Coding Englewood Cliffs, NJ: PrenticeHall 1995.Google Scholar
 Grossmann A, KronlandMartinet R, Morlet J: Reading and understanding continuous wavelet transforms. In Wavelets: TimeFrequency Methods and Phase Space (Edited by: Combes JM, Grossmann A). Tchamitchian Ph 1990, 2–20.View ArticleGoogle Scholar
 Waheed K, Salam FM: A DataDerived Quadratic Independence Measure for Adaptive Blind Source Recovery in Practical Applications. In 45th IEEE International Midwest Symposium on Circuits and Systems, Tulsa, Oklahoma August 4–7 2002, 473–476.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.