Journal Details
Format
Journal
eISSN
1891-5469
First Published
01 Jan 2010
Publication timeframe
1 time per year
Languages
English
Open Access

# Low error Kramers-Kronig estimations using symmetric extrapolation method

###### Page range:147 - 152
Journal Details
Format
Journal
eISSN
1891-5469
First Published
01 Jan 2010
Publication timeframe
1 time per year
Languages
English
Introduction

It is well known that the real and imaginary parts of certain complex function $R(ω)$ $R(\omega )$are related to each other by the KK equations. However, this is only possible if the complex function $R(ω)=R′(ω)+jR′′(ω)$ $R(\omega )={{R}^{\prime }}(\omega )+j{{R}^{\prime \prime }}(\omega )$represents the Fourier transformation of a linear and causal physical process. That is, if a stimulus $s(t)$ $s(t)$is applied to a specific linear physical system, its temporal response $r(t)$ $r(t)$is related to the stimulus applied to an ordinary differential equation with constant coefficients. So, if the stimulus is real, the response will also be real. That system can also be analyzedin the frequency domain. The link between both domains is given by the Fourier transformation:

$R(ω)=∫−∞∞r(t)e−jωtdt$ $$R(\omega )=\int_{-\infty }^{\infty }{r}(t){{e}^{-j\omega t}}dt$$

where $ω=2πf and f$ $\omega =2\pi f\,\,\text{and}\,f$is the frequency in Hz.

Kramers [1] and Kronig [2] developed the equations that bear their names. When working on linear optical spectroscopy, they can be applied to a large number of physical processes and are one of the most general relationships in electrodynamics. KK relationships are established between the real part of certain quantities that characterize the physical dispersion processes (dielectric permittivity, magnetic permeability, conductivity, etc.) and the imaginary part of these quantities that characterize the physical absorption processes.

Let $Z(ω)=Z′(ω)+jZ′′(ω)$ $Z(\omega )={{Z}^{\prime }}(\omega )+j{{Z}^{\prime \prime }}(\omega )$be the electrical impedance of a material medium between two metallic electrodes, where are real functions. If we know $Z′′(ω)$ ${{Z}^{\prime \prime }}(\omega )$and the value of $Z′(ω) for ω→∞, KK$ ${{Z}^{\prime }}(\omega )\,\,\text{for}\,\omega \to \infty ,\,\,\text{KK}$equations allow us to obtain $Z′(ω)$ ${{Z}^{\prime }}(\omega )$as shown in equation (2).

$Z′(ω)=Z′(∞)+(2π)∫0∞xZ′′(x)−ωZ′′(ω)x2−ω2dx$ $${{Z}^{\prime }}(\omega )={{Z}^{\prime }}(\infty )+\left( \frac{2}{\pi } \right)\int_{0}^{\infty }{\frac{x{{Z}^{\prime \prime }}(x)-\omega {{Z}^{\prime \prime }}(\omega )}{{{x}^{2}}-{{\omega }^{2}}}}dx$$

If instead of $Z′(∞),$ ${{Z}^{\prime }}(\infty ),$

we know $Z′(0),$ ${{Z}^{\prime }}(0),$then equation (3) allows us to obtain $Z′(ω)$ ${{Z}^{\prime }}(\omega )$for any values of the angular frequency ω.

$Z′(ω)=Z′(0)+(2ωπ)∫0∞(ωx)Z′′(x)−Z′′(ω)x2−ω2dx$ $${{Z}^{\prime }}(\omega )={{Z}^{\prime }}(0)+\left( \frac{2\omega }{\pi } \right)\int_{0}^{\infty }{\frac{\left( \frac{\omega }{x} \right){{Z}^{\prime \prime }}(x)-{{Z}^{\prime \prime }}(\omega )}{{{x}^{2}}-{{\omega }^{2}}}}dx$$

On the other hand, if we know Z′(ω), equation (4) allows us to obtain Z′′(ω).

$Z′′(ω)=−(2ωπ)∫0∞Z′(x)−Z′(ω)x2−ω2dx$ $${{Z}^{\prime \prime }}(\omega )=-\left( \frac{2\omega }{\pi } \right)\int_{0}^{\infty }{\frac{{{Z}^{\prime }}(x)-{{Z}^{\prime }}(\omega )}{{{x}^{2}}-{{\omega }^{2}}}}dx$$

If Z′′(ω) has the formula:

$Z′′(ω)=a ω1+b ω2$ $${{Z}^{\prime \prime }}(\omega )=\frac{a\,\omega }{1+b\,{{\omega }^{2}}}$$

The integrals (2) and (3) give as a result:

$∫0txZ′′(x)−ωZ′′(ω)x2−ω2 dx=a arctan(b t)b(1+b ω2)$ $$\int_{0}^{t}{\frac{x{{Z}^{\prime \prime }}(x)-\omega {{Z}^{\prime \prime }}(\omega )}{{{x}^{2}}-{{\omega }^{2}}}\,}dx=\frac{a\,arctan(\sqrt{b}\,t)}{\sqrt{b}\left( 1+b\,{{\omega }^{2}} \right)}$$ $∫0t(ωx)Z′′(x)−Z′′(ω)x2−ω2 dx =−ab ω arctan(b t)(1+b ω2)$ \begin{align}& \int_{0}^{t}{\frac{\left( \frac{\omega }{x} \right){{Z}^{\prime \prime }}(x)-{{Z}^{\prime \prime }}(\omega )}{{{x}^{2}}-{{\omega }^{2}}}}\,dx \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=-\frac{a\sqrt{b}\,\omega \,arctan(\sqrt{b}\,t)}{\left( 1+b\,{{\omega }^{2}} \right)} \\ \end{align}

If now Z′(ω) is expressed in terms of ω as shown in equation (8)

$Z′(ω)=c+d1+b ω2$ $${{Z}^{\prime }}(\omega )=c+\frac{d}{1+b\,{{\omega }^{2}}}$$

The integral (4) gives as a result:

$∫0tZ′(x)−Z′(ω)x2−ω2 dx=−db arctan(b t)(1+b ω2)$ $$\int_{0}^{t}{\frac{{{Z}^{\prime }}(x)-{{Z}^{\prime }}(\omega )}{{{x}^{2}}-{{\omega }^{2}}}}\,dx=-\frac{d\sqrt{b}\,arctan(\sqrt{b}\,t)}{\left( 1+b\,{{\omega }^{2}} \right)}$$

In many situations, especially when conducting experiments, a limited number of measured values of Z′ and Z′′ are available in a narrow range of frequencies. This limited frequency range avoids knowing ω →∞ or ω → 0, limiting the validity of the KK equations. The limitation may have different origins, for example, it may be purely technological, such as the Solartron 12508W or HP 4291 impedance analyzers, which have ranges from 10 μHz to 65 kHz and 1 MHz to 1.8 GHz respectively.

A frequency range is also limited when separating a phenomenon of interest from others, reducing the range of measured frequencies, as in the monitoring of α, β and δ dispersions of biological suspensions. For example, all the information needed to characterize the β dispersion of a yeast suspension from an industrial fermenter is in the range of 0.1 MHz to 10 MHz [3]. Other dispersions outside this range can also be observed, but they are not of practical interest. The necessary information to monitor yeast growth is only in the mentioned range.

In the case of yeasts, real and imaginary components for ω → ∞ and ω → 0 can be obtained, simply by assuming that the values outside the working range remain constant. This leads us to assume that there are measurements

between 0 and ∞. This is not the case when impedance is used to monitor electrochemical sensors. Here, the ranges of frequency to observe the behavior of the electrode-electrolyte interfaces must include very low frequencies, difficult to measure experimentally.

For example, to complete the sigmoid curve of interface impedance modelled with the Randles circuit, it is necessary to measure frequencies lower than 0.1 Hz [4]. This is impractical because it is very difficult to keep stable measurements for enough time to complete a cycle at very low frequencies, where a cycle can take several minutes.

All these limitations of information that Z′ and Z′′ give us, play an important role when using KK equations to verify that the system of measurement is working correctly.

Several works about extrapolating experimental results below the frequency at which the imaginary part of the impedance presents the maximum have already been published.

These methodologies are based on the symmetry of the imaginary impedance around the frequency of the maximum and are limited to systems having one time constant. Kendig and Mansfeld [5], for example, obtained the polarization resistance from the following equation:

$Rp=Z′(0)−Z′(∞)=(4ln(10)π)| ∫−∞log(ωm)Z′′(ω) dlog10(ω) |$ \begin{align}& {{R}_{p}}={{Z}^{\prime }}(0)-{{Z}^{\prime }}(\infty ) \\ & =\left( \frac{4\ln (10)}{\pi } \right)\left| \int_{-\infty }^{\log \left( {{\omega }_{m}} \right)}{{{Z}^{\prime \prime }}}(\omega )\,d{{\log }_{10}}(\omega ) \right| \\ \end{align}

Macdonald and Urquidi-Macdonald [6] presented a method of experimental data extrapolation consisting of polynomial fit. The authors evaluated the integral in equation (11) by segments, adjusting the experimental data of Z′′ versus frequency with a fifth-order polynomial, given by equation (12), using the least squares technique.

$R p ≈ 2 π ∫ x min x max Z ′ ′ ( x ) x d x$ \begin{align}& {{R}_{p}}\approx \left( \frac{2}{\pi } \right)\int_{{{x}_{\min }}}^{{{x}_{\max }}}{\frac{{{Z}^{\prime \prime }}(x)}{x}}\,dx\end{align} $Z ′ ′ ( x ) = ∑ i = 0 5 a i x i$ \begin{align}{{Z}^{\prime \prime }}(x)=\sum\limits_{i=0}^{5}{{{a}_{i}}}{{x}^{i}} \end{align}

Equation (11) follows from equations (2) and (3).

The segments on which the integral is evaluated are chosen to coincide with the sign of Z′′ changes or with changes in the gradient of Z′′, This methodology of experimental data extrapolation is also presented in ref [7].

Esteban and Orazem [8] proposed an approach that avoids the inconveniences derived from extrapolations with polynomials and also does not require making assumptions about asymptotical behaviors. Authors use an algorithm to determine the functions of the real and imaginary impedance components in a region of ω for which no experimental data are available. Then, for each frequency, the polynomial functions of K and M order given by

equations (13) and (14) are replaced in equations (4) and (2) respectively.

$Z r ( ω ) = ∑ k = 0 K a k ( 1 ) ( log ⁡ ω ) k$ \begin{align}& {{Z}_{r}}(\omega )=\sum\limits_{k=0}^{K}{a_{k}^{(1)}}{{(\log \omega )}^{k}}\end{align} $Z i ( ω ) = ∑ m = 0 M a k ( 2 ) ( log ⁡ ω ) m$ \begin{align}&{{Z}_{i}}(\omega )=\sum\limits_{m=0}^{M}{a_{k}^{(2)}}{{(\log \omega )}^{m}} \end{align}

Due to the complicated form of the integrands, the integrals were numerically solved.

Another approach shown in the literature is based on the expected asymptotical behavior of a typical electrochemical system [9].

That is, to extrapolate the Z′′(ω) value to ω = 0, it is assumed that the function is proportional to ω, when ω → 0, This is consistent with the behaviour in frequency of Z′′(ω) in a Randles circuit.

For its part, the real part of the impedance tends to a constant value, which is the sum of the electrolytic solution resistance and the charge transfer resistance. The last resistance models the non-linearity of the electrode-electrolyte interface.

Finally, it should be mentioned that the lack of agreement between experimental data and the corresponding KK transformations can be attributed to two factors: on the one hand, the questions associated with the non-linearity or non-causality of the system and, on the other hand, the questions related to the low precision of the measurement systems at a very low frequency.

In order to overcome those limitations, we propose to extrapolate data outside of the available range of measurement, significantly improving the results of applying the KK transformation. The method, which is based on the symmetry of Z′′(ω), can be applied in impedance measurements of the electrode-electrolyte interface (EEI), when the expected behavior responds to the Randles model.

Materials and methods
Electrochemical cell

In order to illustrate the application of the KK transformation and validate the method, measurements were made using a tripolar cell in an electrolytic solution.

Several frequency sweeping experiments were carried out at constant overpotential and using electrodes polished with sandpaper #180, The three-electrode electrochemical cell used (Fig. 1) is composed of a working electrode (WE) and an Ag/AgCl reference electrode (Re1). An AISI 304 stainless steel concave counter-electrode (CE) 85 mm was also used. The WE is a solid cylinder (1.5 cm long) made also of AISI 304 stainless steel with only 1 cm2 of its transversal section exposed; the rest was insulated with Grilon. The CE area was larger than the working electrode to minimize its impedance. The sampling frequency was such that the measurement points were equally spaced in a logarithmic scale taking 5 points per decade. The integration time of the measurements was 1 cycle for frequencies lower than 6.28 rad/s and 16 cycles for the rest.

In every case, the potential was stabilized in an open circuit until the voltage shift was lower than 0.05 mV/seg. Then, the equivalent series resistance and reactance measurements were carried out. The overpotentials used were 10 mV. The electrolyte solution used was NaCl 0.9%.

Electrochemical measurements were performed with a Solartron 12508W Impedance Analyzer composed of a Solartron 1287 Electrochemical Analyzer and a Solartron 1250 Frequency Response Analyzer, commanded by the software provided by the manufacturer (ZPlot®, Scribner Associates Incorporated).

The symmetric data extrapolation method

The method consists of expanding the amount of data available in a range of frequencies greater than that measured, extrapolating data from the available measurements.

The method requires two conditions for its application: (i) the curves of Z′(ω) and Z′′(ω) must be symmetric respect to the central frequency, a condition that is fulfilled by dispersions with a single time constant, and (ii) the frequency range of the measurements must cover at least the central frequency of the dispersion until the percentage change of Z′(ω) is below 0.1, Extrapolation allows us to reduce the difference between the curves of the measured and the calculated data using the KK transformation.

Extrapolation is performed at low frequency, from the lowest limit of frequency ωmin up to the central frequency of the dispersion ωp, taking into account that Z′(ω) − Z′(∞) is symmetric respect to ωp , This frequency ωp is also the frequency in which Z′′(ω) is maximum. ωp , often called the peak frequency.

In this work, we used $ωp$ ${{\omega }_{p}}$to locate the point of $Z′(ω),$ ${{Z}^{\prime }}(\omega ),$It is where its concavity changes, as schematically shown in Figure 4.

Figure 2 graphically shows how to extrapolate frequency and amplitude of $Z′(ω)$ ${{Z}^{\prime }}(\omega )$for points before $ωp.$ ${{\omega }_{p}}.$The corresponding equations are:

$ω 1 = 2 ω p − ω 2$ \begin{align}{{\omega }_{1}}=2{{\omega }_{p}}-{{\omega }_{2}}\end{align} $Z ′ ω 1 = 2 Z ′ ω p − Z ′ ω 2$ \begin{align} {{Z}^{\prime }}\left( {{\omega }_{1}} \right)=2{{Z}^{\prime }}\left( {{\omega }_{p}} \right)-{{Z}^{\prime }}\left( {{\omega }_{2}} \right)\end{align}
Ethical approval

The conducted research is not related to either human or animal use.

Results

Figure 3 shows the plots of $Z′(ω) and Z′′(ω)$ ${{Z}^{\prime }}(\omega )\,\,\text{and }\,{{Z}^{\prime \prime }}(\omega )$versus frequency for AISI 304 stainless steel in NaCl 0.9% solution between 0.628 and 4.084 ∙ 105 (rad/s).

Equation (4) solved numerically is used to obtain $Z′′(ω)$ ${{Z}^{\prime \prime }}(\omega )$from $Z′(ω)$ ${{Z}^{\prime }}(\omega )$as shown in Figure 4, A difference of more than 20% between the measured and calculated data at $ωp=$ ${{\omega }_{p}}=$$2.04 (rad/s) is observed. Zp′′=Z′′(ωp)$ $2.04\,(rad/s)\,is\,observed.\,\,Z_{p}^{\prime \prime }={{Z}^{\prime \prime }}\left( {{\omega }_{p}} \right)$is the maximum of $Z′′ at ω=ωp.$ ${{Z}^{\prime \prime }}\,\text{at}\,\omega ={{\omega }_{p}}.$

Figure 5 shows the expanded dataset of $Z′(ω)$ ${{Z}^{\prime }}(\omega )$after applying the proposed extrapolation method. In it, the interpolated frequency range is indicated with vertical lines $( ωmin to ωp )$ $\left( {{\omega }_{min}} \right.\,to\,\,\left. {{\omega }_{p}} \right)$.

Figure 6 shows the result of applying the KK equation (4) again, but this time to the expanded dataset.

It should be noted that the agreement between maximum of $Z′′(ω)$ ${{Z}^{\prime \prime }}(\omega )$experimental data and the result obtained using the KK equation (4) depends on the lower frequency limit $ωmin.$ ${{\omega }_{min}}.$To analyse this aspect in more detail, Figure 7 shows the percentage change of $Zp′′ (ΔZp′′%)$ $Z_{p}^{\prime \prime }\ \ \ \left( \Delta Z_{p}^{\prime \prime }\text{%} \right)$versus ωmin , $ΔZp′′%$ $\Delta Z_{p}^{\prime \prime }%$is calculated with equation (17).

$ΔZp′′%=[ Zp′′(ωmin)−Zp′′(exp) ] ⋅100/Zp′′(exp)$ \begin{align}& \Delta Z_{p}^{\prime \prime }\text{%}=\left[ Z_{p}^{\prime \prime }\left( {{\omega }_{\min }} \right)-Z_{p}^{\prime \prime }(\exp ) \right] \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\cdot 100/Z_{p}^{\prime \prime }(\exp ) \\ \end{align}
Discussion and conclusions

Ideally, equations (2) through (4) of KK allow us to obtain $Z′(ω) from Z′′(ω)$ ${{Z}^{\prime }}(\omega )\,\text{from}\,{{Z}^{\prime \prime }}(\omega )$and vice versa. In practice, the range of measurement frequencies must be enough to obtain an acceptable error within the required application. Normally, electrochemical measurement systems cover the range of high frequencies $(>ωp)$ $\left( >{{\omega }_{p}} \right)$quite well and the measurements are quite precise. For example, Solartron 1250 in combination with the 1287 electrochemical interface allows performing electrochemical impedance measurements with an error less than 0.3% between 10 μHz and 10 kHz. However, in the low frequencies range $(<ωp),$ $\left( <{{\omega }_{p}} \right),$it is very difficult, since it requires a long measurement time and it is difficult to keep the stability of experimental conditions. From the analysis presented, it is evident that the quality of the results obtained by application of the KK relations, $Z′′(ω) from Z′(ω)$ ${{Z}^{\prime \prime }}(\omega )\,\,\text{from}\,\,{{Z}^{\prime }}(\omega )$in the example presented, depends on the quality of the available experimental dataset and measurement errors. The extrapolation method of missing data presented in this work, as in several published works [3, 5, 8, 10], assume symmetry of the real and imaginary components of the impedance $Z(ω)$ $Z(\omega )$in the frequency domain and a system with a single time constant.

These restrictions are fulfilled in many physical and biophysical systems for certain frequency ranges.

However, unlike published works, the method presented does not use polynomial fit or assumptions about asymptotic behaviors. For frequencies lower than ωp, the method replicates the behaviour observed for $ω>ωp.$ $\omega >{{\omega }_{p}}.$This simple fact, described in Figure 2, allows us to extrapolate enough points to comply with the requirements needed to apply the KK equations.

A plot of measured data and its corresponding KK transformations using equation (4) are shown in Figure 4, It can be observed a mismatch due to insufficient measurement range.

Figure 7 shows how the theoretical adjustment depends on the amplitude of the extrapolated frequency range, by representing the percentage change of $Zp′′$ $Z_{p}^{\prime \prime }$as a function of the lower limit of the extrapolated dataset $ωmin$ ${{\omega }_{\min }}$It can be seen that the worst case (red point) corresponds to an error of 24% when using only measured data. On the other hand, when extrapolating data at frequencies lower than $ωmin=$ ${{\omega }_{min}}=$$10−2rad/s$ ${{10}^{-2}}rad/s$(blue point) the error dropped below 1%.

In Figures 5 and 6 data were extrapolated up to $ωmin=$ ${{\omega }_{min}}=$$10−5rad/s$ ${{10}^{-5}}rad/s$, achieving an error less than 0.2%.

In conclusion, we present a simple method to extrapolate impedance data sets from linear, causal systems and with a single time constant. The new method allows us to apply the KK equations with a minimum estimation error of the real or imaginary parts of impedance.

H. Kramers, La diffusion de la lumière par les atomes, Atti del Congresso Internazionale dei Fisici, vol 2 (Bologna: Zanichelli Editore), (1927) 545.Kramers H. La diffusion de la lumière par les atomes, Atti del Congresso Internazionale dei Fisici, vol 2 (Bologna: Zanichelli Editore) 1927 545Search in Google Scholar

L. Kronig, On the theory of dispersion of x-rays, J. Opt. Soc. Am., 12 (1926) 547-557, 10.1364/JOSA.12.000547Kronig L. On the theory of dispersion of x-rays J. Opt. Soc. Am 12 1926 547 557 10.1364/JOSA.12.000547

Yardley JE, Kell DB, Barrett J, Davey CL, On-line, real-time measurements of cellular biomass using dielectric spectroscopy, Biotechnol Genetic Eng Rev 17 (2000) 3-35, 10.1080/02648725.2000.10647986Yardley JE Kell DB Barrett J Davey CL On-line, real-time measurements of cellular biomass using dielectric spectroscopy Biotechnol Genetic Eng Rev 17 2000 3 35 10.1080/02648725.2000.10647986

C.C.M. Martinez, E.F. Treo, R.E. Madrid, C.C. Felice, Evaluation of chronoimpedance technique as transduction method for a carbon paste/glucose oxidase based glucose biosensor, Biosens. Bioelectron. 26 (2010) 1239-1244, 10.1016/j.bios.2010.06.033Martinez C.C.M. Treo E.F. Madrid R.E. Felice C.C. Evaluation of chronoimpedance technique as transduction method for a carbon paste/glucose oxidase based glucose biosensor, Biosens Bioelectron 26 2010 1239 1244 10.1016/j.bios.2010.06.033

M. Kendig and E Mansfeld, Corrosion rates from impedance measurements: an improved approach for rapid automatic analysis, Corrosion, 39 (1983), 466-467, 10.5006/1.3581908Kendig M. and Mansfeld E Corrosion rates from impedance measurements: an improved approach for rapid automatic analysis Corrosion 39 1983 466467 10.5006/1.3581908

D. D. Macdonald and M. Urquidi-Macdonald, Application of Kramers-Kronig transforms in the analysis of electrochemical systems, J. Electrochem. Soc., 132 (1985), 2316-2318, 10.1149/1.2113570Macdonald D. D. and Urquidi-Macdonald M. Application of Kramers-Kronig transforms in the analysis of electrochemical systems J. Electrochem. Soc 132 1985 23162318 10.1149/1.2113570

M. Urquidi-Macdonald, S. Real, and D. D. Macdonald, Applications of Kramers-Kronig Transforms in the analysis of electrochemical impedance DATA-III. Stability and linearity, Electrochim. Acta, 35 (1990), 1559-1566, 10.1016/0013-4686(90)80010-LUrquidi-Macdonald M. Real S. and Macdonald D. D. Applications of Kramers-Kronig Transforms in the analysis of electrochemical impedance DATA-III Stability and linearity, Electrochim. Acta 35 1990 15591566 10.1016/0013-4686(90)80010-L

J. M. Esteban and M. E. Orazem, On the application of the Kramers-Kronig relations to evaluate the consistency of electrochemical impedance data, J. Electrochem. Soc., 138 (1991), 67-75, 10.1149/1.2085580Esteban J. M. and Orazem M. E. On the application of the Kramers-Kronig relations to evaluate the consistency of electrochemical impedance data J. Electrochem. Soc 138 1991 6775 10.1149/1.2085580

C. Haili, M.S. Thesis, University of California, Berkeley, CA (1987).Haili C. Thesis M.S. University of California Berkeley, CA 1987Search in Google Scholar

P. Agarwal, M.E. Orazem and L.H. Garcia-Rubio, Application of measurement models to impedance spectroscopy: III. Evaluation of Consistency with the Kramers-Kronig Relations, J. Electrochem. Soc., 142 (1995), 4159-4168, 10.1149/1.2048479Agarwal P. Orazem M.E. and Garcia-Rubio L.H. Application of measurement models to impedance spectroscopy: III. Evaluation of Consistency with the Kramers-Kronig Relations J. Electrochem. Soc 142 1995 41594168 10.1149/1.2048479

Recommended articles from Trend MD