Sunday, June 30, 2019

DFA alpha1, HRV complexity and polarized training

An exciting aspect of heart rate variability assessment above and beyond the conventional HF and LF power analysis deals with the fractal complexity of the RR interbeat interval over a short time span.  To better understand what "fractal complexity" means, it may be helpful to first discuss the mathematical construct known as a fractal.  The description of a fractal is not easily defined - "The consensus is that theoretical fractals are infinitely self-similar, iterated, and detailed mathematical constructs having fractal dimensions, of which many examples have been formulated and studied in great depth".  One of my favorite examples is that of a coastline, no matter how much you zoom in on a map, the pattern is similar to the previous magnification.
Perhaps a picture (from an excellent review by the originators of the DFA concept) would be helpful to illustrate an example in biology:
You can see that as you enlarge either the area in space (on the left) or the ECG tracing in time (on the right), the result looks quite similar to the original.  They both have self similar, repeated and detailed structure.  When that structure becomes less detailed and complex, it may be due to disease states or other abnormalities.
For instance, below is another example of the increase of, or loss of complexity in certain cardiac diseases:

To take a look at the math and theory behind one computation of the randomness/complexity of the signal, here is a passage from the review:
An important methodologic challenge is how to detect and
quantify the scaling and correlation properties of physiologic
time series, which are typically not only irregular, but also
nonstationary (i.e., their statistical properties change with time).
To help deal with the ubiquitous biologic ‘‘complication’’ of
nonstationarity, we have introduced a modified rms analysis of
a random walk—termed detrended fluctuation analysis (DFA;
refs. 19 and 26–29). To illustrate the DFA algorithm briefly,
consider the cardiac interbeat series shown in Fig. 4A. First, the
original time series is integrated, and then divided into boxes of
equal length, n. For each box of length n, a least squares line
(representing the trend in that box) is fit to the data (Fig. 4B). For
a given box size n, the characteristic size of the fluctuations,
denoted by F(n), is then calculated as the rms deviation between
y(k) and its trend in each box. This computation is repeated over
all time scales (box sizes). Typically, F(n) will increase with box
size n. A linear relationship on a log-log graph indicates the
presence of scaling (self-similarity), such that fluctuations in
small boxes are related to the fluctuations in larger boxes in a
power-law fashion. The slope of the line relating log F(n) to log
n determines the fractal scaling exponent, (Fig. 4C). This
exponent provides a measure of the ‘‘roughness’’ of the original
time series: the larger the value of , the smoother the time
series. In this context, 1 f-like noise ( 1) can be interpreted
as a ‘‘compromise’’ between the complete unpredictability of
white noise ( 0.5) and the much smoother ‘‘landscape’’ of
Brownian noise ( 1.5; ref. 19).
The DFA alpha1 is a short term index of RR interbeat fractal complexity usually calculated from 4-12 beats.  
Here is the interpretation of values from the Kubios website:

Values of 1 to 1.5 represent noise with complexity, whereas values near .5 are less complex white noise. 

Here is another description of noise and correlation:

"Relationship Between Self-Similarity and Auto-Correlation Functions
The self-similarity parameter of an integrated time series is related to the more familiar auto-correlation function, tex2html_wrap_inline1157, of the original (non-integrated) signal. Briefly:
  • For white noise where the value at one instant is completely uncorrelated with any previous values, the integrated value, y(k), corresponds to a random walk and therefore tex2html_wrap_inline1161 [7, 17]. The auto-correlation function, tex2html_wrap_inline1157, is 0 for any tex2html_wrap_inline1165 (time-lag) not equal to 0.
  • Many natural phenomena are characterized by short-term correlations with a characteristic time scale, tex2html_wrap_inline1167, and an autocorrelation function, tex2html_wrap_inline1157 that decays exponentially, [i.e., tex2html_wrap_inline1171]. The initial slope of tex2html_wrap_inline1173 vs. tex2html_wrap_inline1175 may be different from 0.5, but tex2html_wrap_inline1001 will approach 0.5 for large window sizes.
  • An tex2html_wrap_inline1001 greater than 0.5 and less than or equal to 1.0 indicates persistent long-range power-law correlations, i.e., tex2html_wrap_inline1189. The relation between tex2html_wrap_inline1001 and tex2html_wrap_inline1193 is tex2html_wrap_inline1195. Note also that the power spectrum, S(f), of the original (non-integrated) signal is also of a power-law form, i.e., tex2html_wrap_inline1199. Because the power spectrum density is simply the Fourier transform of the autocorrelation function, tex2html_wrap_inline1201. The case of tex2html_wrap_inline1203 is a special one which has interested physicists and biologists for many years--it corresponds to 1/f noise (tex2html_wrap_inline1207).
  • When tex2html_wrap_inline1209, power-law anti-correlations are present such that large values are more likely to be followed by small values and vice versa [10].
  • When tex2html_wrap_inline1211, correlations exist but cease to be of a power-law form; tex2html_wrap_inline1213 indicates brown noise, the integration of white noise.
The tex2html_wrap_inline1001 exponent can also be viewed as an indicator of the ``roughness'' of the original time series: the larger the value of tex2html_wrap_inline1001, the smoother the time series. In this context, 1/f noise can be interpreted as a compromise or ``trade-off'' between the complete unpredictability of white noise (very rough ``landscape'') and the much smoother landscape of Brownian noise [18]."

Although this is interesting academically, is there any practical value in DFA a1 during dynamic exercise?  

Gronwald and associates investigated the change in DFA a1 (RR fractal complexity) during exercise in young cyclists.   An incremental ramp test (20 w per 3 min) was performed until exhaustion.

Sixteen trained cyclists (age: 25.9 ± 3.8 years; height: 180.7 ± 6.1 cm; body mass: 77.4 ± 8.2 kg; body fat: 12.3 ± 3.4%) were recruited from local sports clubs. Eligible subjects had to fulfill the following inclusion criteria: non-smokers and at least 8 h cycling training per week within the last 6 months.

Methods (DFA a1 related):
In this study, we only computed the short term scaling exponent (window
width: 4 ≤ n ≤ 16 beats) due to the relatively short recording times for each physiological condition (see Goldberger et al., 2002; Hautala et al., 2003; Tulppo et al., 2001). DFA-alpha1 values reveal correlation properties of heart rate time series, i.e. type of noise: approx. 0.5 for uncorrelated white noise with random signals, approx. 1 for a mix
of uncorrelated and maximally correlated signal components with 1/f noise and approx. 1.5 for strongly correlated Brownian noise (Platisa & Gal, 2008).

Here is the entire DFA a1 response over the range of VO2 (load):
The authors felt there was a biphasic response, but the drastic reduction with moderate to high intensity is impressive.

Here is a more detailed look at the higher VO2 load range:
Values of .2 to .4 are reached at maximal VO2.

Their conclusions were as follows:
The gradual decrease of DFA-alpha1 with increasing
exercise intensity (with values of around 0.3 at exhaustion) indicates an intensity dependent change from strongly correlated to uncorrelated/stochastic or anti-correlated behaviour of the RR intervals (Platisa & Gal, 2008). This is consistent with previous studies reporting an almost linear reduction of complexity and approximation of the RR data structure towards a merely random signal.
Maybe in the future, non-linear analysis of HRV and especially DFA-alpha1 could be an easy detectable individual measure for programming and training control to get new insights into dose-response relationship at different exercise intensities and volumes. This could be an important issue for the application in health care, primary and secondary prevention as well as in sports therapy to avoid overload in different risk groups.
The concept of using this to avoid overload and over training will be discussed shortly, but the idea seems sound.  One of my original critiques of using muscle O2 as a marker of training overload (or effort readiness) was that there is no valid physiologic underpinning that would support this.  In contrast, the loss of RR interbeat complexity is a well researched field in cardiac disease and it's application could spill over to endurance sports.

The same group looked at DFA a1 in trained cyclists performing HIT intervals, particularly the trends during active recovery.
From the paper:
It was the aim of the present study to analyze the influence of short-term interval sessions combined with active recovery periods on standard time-domain measures and non-linear dynamics of HRV. The main objective was to determine whether active recovery is effective to prevent the loss of fractal organization in the short-term scaled characteristic of cardiac autonomic activity in endurance trained cyclists.
The HIT intervals were 60 second maximal efforts x 5 with 3 sets (in other words - intense). 

  • As expected, the lactate and heart rate rose during the HIT intervals.  
  • The DFA a1 did drop during the intervals as (expected) and rose during active recovery:

HR, Lactate:

Here is the DFA a1 response:

The summary of results:
  • I yellowed in the DFA a1 drops during the intervals but also notice how DFA a1 rose during recovery.

Conclusions (from the discussion):
  • The present data combining interval loads with active recovery periods show that the DFA-alpha1 values return to the level of the warm-up periods very quickly during active recovery and remain at nearly the same values until the end of the active recovery phase in endurance trained cyclists.
  • In addition, the results of the present study indicate that the assessment of DFA-alpha1 allows a distinction between different exercise intensities in short-term cycling interval sessions.
  • In the present study, a decrease in DFA-alpha1, which corresponded with increasing exercise intensity, verifies a demand-dependent change from strongly correlated to uncorrelated/stochastic or anti-correlated behavior of the RR-intervals and back to strongly correlated behavior in early recovery.

My data and how to "do it yourself"
Now that we have reviewed some published work on using DFA a1 during dynamic exercise, let's have a look at some of my data.

Going back to the ramp test I showed in the prior post, is the DFA a1 response similar?
According to the above, it should steadily decline with load.   

It certainly looks like it does:
DFA a1 in black, HR is blue (data from Hexoskin with artifacts about 1%).
  • This ramp was 30w/minute, a faster rise, short ramp time than the Gronwald study, but the trend is the same.
  • As heart rate rises (with power), the DFA a1 drops to a low of .2 which is quite similar to the ramp study of Gronwald.

MLSS vs MLSS-15w:
I was interested to see if DFA a1 could be used as a surrogate marker of MLSS as SmO2 desaturation or ventilation rate can be.  
The following tracing consists of 2 intervals, the first at about 15w below MLSS and the second just above it.

Here is the data profile, first with rectus femoris muscle O2 saturation then with costal (serratus) O2 saturation:

Rectus Femoris O2 sat:

Costal O2 sat:
  • Both RF and Costal muscles remain steady state during the 235w interval.
  • However they continue to desaturate through 6 minutes of 250w consistent with an exercise load above MLSS.

  • This breakpoint is also confirmed by the ventilation response (Hexoskin) that is stable at 235w but not at 250w:

Does DFA a1 show a breakpoint or different response between the two interval powers?
Heart rate in blue, DFA a1 in black:
  • The responses are very similar and cluster around a value of .35, definitely in the range associated with high intensity.
  • There does not seem to be a difference in DFA a1 between a load just below and above the MLSS.  The values are both in the uncorrelated, loss of complexity range.

I did repeat this and if anything the DFA a1 was suppressed more at the lower power.  Here is the tracing integrated into the entire session.  There was ample warm up.
  • The DFA a1 was still suppressed to a similar degree at the MLSS and MLSS-20w.

Are there breakpoints at lower power?
Can we find the power/heart rate upper limit for a low intensity exercise session?
Since there does not appear to be a defined step down of DFA a1 at the MLSS, perhaps there is one at a lower power, possibly near VT1.  In addition, as mentioned in the Gronwald paper discussion, it may be possible to utilize DFA a1 as an indicator of low load endurance training effort.  Knowing your "threshold of interbeat complexity loss" could allow you to do the large volumes of low intensity training recommended by many sources without undo cardiac stress.  To that end, less examine a longer interval ramp done at sub maximal power.
The following is a restricted low power, longer interval ramp I did looking for the DFA a1 shift into a less complex interbeat interval:
After a warm up, each segment was about 4 minutes with a 10 watt step up, the power and average heart rate were as follows:

DFA a1 results:
The top pane shows heart rate only, divided into the 4 minute intervals.
The lower pane contains the DFA a1 results (in black), the heart rate (in blue) as well as the label for watt power per interval.
A closer look at the last power zone of 218w shows that DFA a1 is already near maximal complexity loss at a value of .45.  Although continuing the ramp further may have driven the number further downward, it is now similar to that of the MLSS interval above:

  • The first two zones (152, 162w) do not appear to have a change in DFA a1 value to any significant degree. 
  • However, in the middle of the 172w interval, there is a mild drop.
  • The 183 watt interval seems to cause the largest per zone decrease in DFA a1.  
  • From here, it continues to drop during the 195w interval but seems to plateau at the 218w area to a nadir similar to that of a maximal effort or MLSS.
  • There is prompt recovery in DFA a1 after the last moderate intensity interval.
  • Whether or not the 183w zone represents a breakpoint change at or near VT1 is unknown.  
  • There certainly does seem to be a power level and heart rate which when exceeded causes DFA a1 to fall. 
  • It appears that even moderate power levels (sub MLSS) done for a long enough time can drive DFA a1 to near nadir values.

I did go back to the data of Gronwald and it is similar to my own:
  • I yellowed in the DFA a1 area of sudden drop (on the bottom line), along with the corresponding load (50 to 58% max) and lactate change (1.24 to 1.57 to 2.07).  Although VT1 was not mentioned specifically, I wonder if it occurs around this point in the ramp.
  • Regardless, it seems that DFA a1 may be a marker for the transition from light to more moderate intensities.
Doing this on your own:
A measure of DFA a1 drop can be shown with the "free" version of Kubios as well.  You will not have the nice time varying graph, but you still can get the DFA a1 values at a particular time and span.  Make sure you are using a good quality chest strap (Polar H10) and have either put cream or water on the conductors.
You should do a ramp starting at a low power level (120watts/HR-110) then every 4 minutes increase by 10 watts or 5 BPM up until your MLSS/FTP is reached.  
Save the .fit file
Download and install Kubios and Matlab runtime.
Open the Kubios software.
Go to File-open and select the .fit file that contains the ramp data you have done.

The Kubios window will look like this:
Notice several areas
1 - The start and length of the analysis window (circled in red in upper pane #1).  The software will calculate the HRV over the time chosen.  In this example the window is in the first minute of the recording.
2 - The Nonlinear tab (circled as #2)- make sure you have pressed it since the DFA a1 will be there.
3 - The DFA a1 result - (#3 circled)

Next - Move the analysis window to the end of each ramp interval(we need to get the DFA a1 value for each power segment).  This is done by dragging the blue zone to the right:
Drag the blue analysis area on the top pane to the end of your first interval (make sure it does not overlap into the next interval).  If you are using a Garmin product, the lap buttons will serve as interval markers (circled in red).
Write down the DFA a1 with the corresponding watt/heart rate of each interval.  Repeat the drag of the blue area on the top pane for each power segment.  Make sure you place the analysis window at the end of the interval.
You can also get the mean heart rate on each power segment by looking at the Time Domain tab:

Make a simple graph of the DFA a1 vs watt or heart rate (the values may be different than the time varying tab since the analysis window is longer).  You can experiment with longer or shorter windows but they are usually 45 to 60 seconds in length.
Look for a drop off:


DFA a1 during intervals near VO2 max?
The DFA a1 response at near maximal levels should be predictable and result in values near .5 (white noise - uncorrelated) or lower.  However, it will also be interesting to see the recovery at different power levels.  The following is an example of DFA a1 response during a fast start 3 minute maximal aerobic power (MAP) interval, followed by pace just below MLSS (well above active recovery levels) then finally active recovery at about 170w:

  • The MAP interval drives the DFA a1 to a nadir near .3, which is maintained as the cycling power continues at near MLSS.  
  • It does not rise until power drops to about half of the VO2 max equivilant.  This is also the same power level corresponding to the transition to a more or less complex RR interbeat interval shown on the ramp above.
  • As discussed above, DFA a1 could be used as an index of "duration of high stress" especially in sports that don't have power meters.  Since some exercise sessions are made up of complex segments of high, moderate and low power, the total time of DFA a1 suppression could be used as an index of accumulated high intensity effort.
Further evidence as a marker of recovery:
If DFA a1 signifies interbeat complexity loss during cardiovascular stress, it should quickly vary unless some disease process is present or the person is not acclimated to intense exercise.  To illustrate this I did a near maximal 3 minute interval followed by a short recovery, another 45 second 300 watt section then a final recovery.  This piece of cycling was at the end of a 50 mile, 3 hour ride in 95F humid conditions, so stress levels should have been high.
Here is the summary:
  1. Maximal 3 minute interval (red arrow down showing drop in DFA a1)
  2. Brief recovery at 130w for about 45 seconds (DFA a1 rise - red arrow now up)
  3. Another 45 second burst of 300w (DFA a1 drop - red arrow down)
  4. Final recovery at 130w with rise of DFA a1 above 1.2.
The DFA a1 value certainly follows the segments of effort and recovery.  The red arrows illustrate how well this tracks each time load rises and falls (similar to the Gronwald HIT study):
Heart rate in blue, DFA a1 in black

DFA a1 as a marker of mild to moderate intensity - A practical example:
Conversely, can we use this index as a guarantee that we are staying within a true, low stress recovery mode of exercise?  For instance, the day after a challenging interval session or race, you may want to still get some training in but need to stay out of an intensity that could affect your recovery.  What better way than making sure DFA a1 stays above the more intense zones.
Here is an example.  
From several prior DFA a1 session reviews, I already had an idea what power/HR (170w with HR below 130) would be "safe".  This is what the DFA a1 looked like:

  • DFA a1 stays above 1 throughout the ride.  Hopefully this will allow optimal recovery to occur.
Why a need for alternate "method of training-intensity quantification"
According to the recent study, the training distribution in a group of athletes was significantly different depending on how each zone was created.  
The training zones were demarcated by the running speed and heart rate corresponding to the gas exchange threshold (GET) and the respiratory compensation threshold (RCT) and fixed RPE values (training zone 1: 0-4, training zone 2: 5-6, training zone 3: 7-10). The TID (total and percentage of time spent in each training zone) of all running training sessions were quantified using continuous running speed and heart rate monitoring, and RPE.
With the Result:
The results from the present study suggest that the method of training intensity quantification affects the computation of the TID in highly-trained middle-distance runners. In particular, the total training time spent in three pre-defined training zones demarcated by RPE provided significantly greater training time in training zones two and three compared to heart rate and running speed time in zone calculations. Furthermore, a greater training time was spent with a heart rate in zone 2 compared to running speed despite a similar training time in zone 1, while heart rate tended to accumulate less time in training time in zone 3 compared to running speed. These findings are important for athletes, coaches and sport science practitioners who monitor the TID of their training programs.

  • Defining the training zone is not as simple as one may think!
  • According to the study, a huge variation of time spent in zone 2 was a result of the different ways of quantifying the target areas.

Examples from other athletes:
The relationships noted above are consistent with the published studies reviewed.  However, I wanted to look at other individual tracings to make sure my data was not skewed on the basis of age, genetics or type of training.
I was fortunate enough to get some data from my friend the XCSkier who is a former member of a national cross country ski team and is still in superb condition.  He provided some ramp and recovery sessions.

The ramp:
He did a lactate test (running) with intervals of speed/heart rate below and above the MLSS.  
Here is the interval sequence, with lactate values after interval 2, 3 and 4: 
Heart rate in blue, DFA a1 in black. 
  • I drew a green line through DFA a1 at the .4 level.  
  • There does not seem to be a breakpoint in DFA a1 at the MLSS (lactate 3.5 to 6.7).  In addition, even though the heart rate is near max, the DFA a1 does not drop any lower than it it at just below the MLSS.  This is very similar to my results.
  • However, there may be a breakpoint at the lactate 1.3 to 3.5 transition.  Whether this represents a shift near VT1 is unclear.

He also sent me a recovery skiing session:

  • For the most part, the DFA a1 stays above .8 (red line).  Only on a few occasions does it get into the .4 zone, which are areas associated with higher heart rates (the green circles with low DFA a1 corresponded with hills). 

Potential use of DFA a1 to enforce polarized training
Staying out of the moderate zone:
Over the past several years multiple opinions have been given on what an ideal training intensity distribution should consist of.  Without getting into specifics, most would agree that the bulk of one's effort should be in the "easy" zone.  The question is, what should that load cutoff be?  Clearly, efforts at or above MLSS are going to be stressful and can't be considered as part of the easy majority of training.  We have already seen from the above that DFA a1 can be used as an index of the "stress of training" and seems sensitive to even moderate loads.
I looked at a previous, 20 minute sustained effort of mine to see if DFA a1 could be helpful in distinguishing the "no man's land" of the light-middle zone.  This is generally a moderate pace below MLSS, but above VT1 (talk test).  In polarized training, this is usually not a pace that much time should be spent doing.  Perhaps a reason to avoid this area is that it is still stressful (from the HRV standpoint) and could lead to overtraining with excessive time in that zone.

This session consisted of a fast start (1 min at 350w), followed by 9 minutes at just above MLSS then finally 10 minutes at 210w which is close to 20% below my MLSS:
I pasted the Power profile in the top pane.
The bottom pane is heart rate (blue) and DFA a1 (black).
The heart rate was recorded with a Hexoskin and artifacts were generally below 3%.

  • During both the high power 1 minute, the 9 minute MLSS+ and 10 minute MLSS -20%, DFA a1 stays below .5 (red circle).  
  • The second red circle on the right was during a Wingate 60 second effort (also with DFA a1 below .5).
  • There is rapid recovery of DFA a1 after both the 20 minute and 1 minute intervals.
  • It is possible that even at moderate intensities, DFA a1 will remain suppressed indicating loss of interbeat complexity and cardiac stress.  The suppression in this case was the same as the MLSS and maximal 1 minute pace.  This situation may need to be taken into account with regards to training distribution balance.

The Astroskin:
Lastly, a look at the Astroskin data from a short ramp test, a strenuous walk uphill and a weight training session that I did with my son a few weeks back.  He does regular strength training, walking but no cycling.

Cycling Ramp:
The following data was from a ramp on a bike trainer, each interval was 4 minutes with the last spike in HR (blue) being a 10 second maximum of about 700 watts.
The top pane is heart rate in blue.
The bottom pane is heart rate (blue) and DFA a1 (black).

  • The DFA a1 does drop with increasing effort and seems to have a breakpoint near the end of the 180 watt section.  It did drop further at 200 watts and ended with a 10 sec maximum sprint.
  • There is quick recovery in the DFA a1 to over 1.2 within 1 minute of stopping the ramp. 
  • One notable observation was the zero percent artifact through the entire session!  The Astroskin data quality is superb.
Practical value:
According to both the cycling ramp (top figure) and walk (below it), it seems that DFA a1 begins to drop in earnest at a heart rate over 120:

  • Therefore in this particular individual, large volume, low intensity training should take place below this heart rate.

DFA a1 during fast walking uphill:
On the same day as the ramp above, my son did a weight training session followed by a walk home.  The walk was up one of those extreme San Francisco hills.  Toward the end we both were going as fast as possible.  I went back and looked at the DFA a1, which was quite interesting.
The top pane is heart rate in blue.
The bottom pane is heart rate (blue) and DFA a1 (black).

  • At mild/moderate intensity there is some DFA a1 reduction, near .65.
  • However at max intensity DFA a1 drops to .5 indicating loss of signal complexity.  This is associated with near max heart rate.
  • There is prompt correction upwards during recovery indicating return of normal homeostasis.
Weight training session:
I also recorded a weight training session using the Astroskin.  This involved a complex, intense routine with kettle bells as described in the post.
The artifact percent was higher than usual (1-10%) probably from the twisting body motion.

Here is a DFA a1 analysis of the 4 sets:
The top pane is heart rate in blue.
The bottom pane is heart rate (blue) and DFA a1 (black). Green rectangles were drawn around the active sets.
  • With progressive higher intensity, heart rate rose and DFA a1 dropped.  The lowest DFA a1 was close to .5.
  • There was good recovery of interbeat complexity during between set rests.
  • Similar to the Gronwald ramp study, heart rate was correlated to DFA a1 reduction.

Summary and Conclusions:

  • DFA a1 is a measure of the fractal complexity of the ECG RR interbeat series.  Values near .5 indicate loss of signal complexity equating to non correlated white noise.  The index does not depend on long recording intervals and can be used during dynamic exercise.
  • The degree of DFA a1 complexity loss corresponds to the exercise load and heart rate elevation.  Maximal loads will reduce the value/complexity more than lighter loads.
  • Under normal circumstances and health status, DFA a1 suppression should reverse relatively quickly when the exercise is stopped.
  • Although not seen with the data above, failure to return to a more complex pattern (DFA a1 rise) after significant efforts could be a sign of systemic disease or maladaptation.
  • There is a general "dose response" curve of DFA a1 change with exercise load, but it may or may not behave in a linear step wise fashion.
  • There does not seem to be a significant change in interbeat complexity just below and above the MLSS.
  • However, there may be some sort of breakpoint well below the MLSS (perhaps VT1) that could serve as a demarcation of mild to moderate intensity.
  • Using DFA a1 trends to avoid exercising in moderate intensity zones could potentially optimize a polarized training regime.


No comments:

Post a Comment