Early detection of metabolic and energy disorders by thermal time series stochastic complexity analysis.

Maintenance of thermal homeostasis in rats fed a high-fat diet (HFD) is associated with changes in their thermal balance. The thermodynamic relationship between heat dissipation and energy storage is altered by the ingestion of high-energy diet content. Observation of thermal registers of core temperature behavior, in humans and rodents, permits identification of some characteristics of time series, such as autoreference and stationarity that fit adequately to a stochastic analysis. To identify this change, we used, for the first time, a stochastic autoregressive model, the concepts of which match those associated with physiological systems involved and applied in male HFD rats compared with their appropriate standard food intake age-matched male controls (n=7 per group). By analyzing a recorded temperature time series, we were able to identify when thermal homeostasis would be affected by a new diet. The autoregressive time series model (AR model) was used to predict the occurrence of thermal homeostasis, and this model proved to be very effective in distinguishing such a physiological disorder. Thus, we infer from the results of our study that maximum entropy distribution as a means for stochastic characterization of temperature time series registers may be established as an important and early tool to aid in the diagnosis and prevention of metabolic diseases due to their ability to detect small variations in thermal profile.

Introduction

A study performed at least 50 years ago demonstrated that reduced variability and linear responses of biological signals, by the statistical and mathematical analysis of intrauterus heart rate variability, would be a predictive value for future underlying disease manifestation (1). Since then, chaos theory has described elements manifesting behavior that is extremely sensitive to initial conditions, does not repeat itself, and yet is deterministic; a sequential series register of several biological parameters such as blood pressure measurements, brain electrical activity, and renal ion transport, and possibly temperature, presents nonlinear unpredictable chaotic behavior, a kind of order without periodicity, and several investigators have implemented studies trying to define and distinguish normal variability from pathological behavior profiles (2). Simultaneously, in the past few years, along with the development of descriptive statistics and computer tools, new models of variance analysis were proposed, and studies which were analyzed by interpretation of averages and deviations began to include new mathematical concepts such as geometric models, stochastic time series, chaos theory, and spectral analysis (3-6).

Observation of the thermal registers of core temperature behavior, in humans and rodents, permits identification of some characteristics of time series such as autoreference and immobility that fit adequately into a stochastic analysis. These characteristics are associated with the physiological nature of an internal thermal series, which is modulated by a complex superstructure that acts at the behavioral level up to molecular signaling (7,8). In this case, this regulatory superstructure reflects the importance of central thermal homeostasis for keeping temperature within a narrow variation range in mammalian species. In fact, among all the biological signals measured in humans, temperature is the one that has the smallest relative deviation (±0.5°C around the mean value of 36.7°C) (9). The hierarchical relevance of the central nervous system (CNS) to thermal control activities is explained by fine adjustments of maximum enzyme reaction activities to a given temperature in homeothermic animals.

In homeothermic animals, the inability to cope with the excess energy absorbed by a hypercaloric diet leads to the phenomenon known as lipotoxicity, a concept interpreted as the cost of keeping the temperature at a maximum but stable level, in order to promote increased fatty acid oxidation and heat release. To maintain thermal stability, animals develop obesity and impaired glucose tolerance, responses associated with reduced thermogenesis (10). Thus, we may hypothesize that the energy imbalance may reflect, at an earlier stage, a subtle change in the variability of the thermal series of these animals. However, from what we know so far, no studies have been performed using complexity stochastic analysis to define the temperature behavior in normal and pathological models. To identify this change, we have used, for the first time, a stochastic autoregressive model, the concepts of which match those associated with physiological systems involved (11,12) and applied in high-fat intake diet (HFD) rats compared with their appropriate standard food intake controls. The HFD model is extremely efficient for creating pathophysiological conditions such as hyperleptinemia, peripheral insulin resistance, diabetes mellitus, and obesity that lead to long-term metabolic and energy disorders. Blocking insulin and leptin pathways at the level of the hypothalamic nuclei promotes increased thermogenic activity in the body, leading to adjustments in the CNS to maintain temperature stability. Therefore, the purpose of the present study was two-fold: a) to identify a stochastic process associated with thermal time series involved with the HFD procedure, and b) from characterization of this stochastic process, to propose a model of an algebraic system as a tool for early and efficient detection of changes in physiological signals.

Material and Methods

Animals

General guidelines established by the Brazilian College of Animal Experimentation (COBEA) and approved by the Institutional Ethics Committee (CEEA/UNICAMP 1697-1) were followed throughout the investigation. Our local colonies originated from a breeding stock supplied by CEMIB/UNICAMP, Campinas, SP, Brazil. Experiments were conducted on age-matched, male offspring of sibling-mated Sprague Dawley rats (n=7 for each group). At 6 weeks of age (180.5±11.7 g body weight), the animals were housed in individual cages and maintained under controlled temperature (22°C) and lighting conditions (7:00 am to 7:00 pm), with free access to tap water, and were randomly distributed into two dietary groups: a control group fed with pelleted standard rodent laboratory chow (standard diet, SD; Nuvital, Brazil) or a high-fat diet group (HFD) fed a high-fat diet for 4 weeks, up to 10 weeks of age. The standard chow diet contained 11.9% kcal as fat and a total of 2.9 kcal/g, and the high-fat diet contained 58.3% kcal as fat and a total of 5.44 kcal/g. All male rats were weighed weekly and had food and water intake measured daily throughout the experiment. After an adjustment Subcue period of 48 h, the 10-week-old rats were anesthetized with a mixture of ketamine (75 mg/kg body weight, ip) and xylazine (10 mg/kg body weight, ip), and a thermal probe (Data Logger, USA) was placed inside the abdominal cavity. The recording of internal temperature was scheduled to start at 1 week after the surgical proceedings, and data were recorded for 6 h in both the SD and HFD groups, with a 1-min period between recordings. During the experiments, the animals were kept isolated in an environment of low external stimuli and monitored with respect to the acceptance of diet and vitality. At the end of the assay, a lethal dose of thiopental was administered and data were transferred to the program Subcue-Data Logger Mathematica 7.0 (Wolfram Research, Inc., USA).

Assumptions

The chosen time interval for the recorded temperature of each pair of animals submitted to analysis was between 0:00 and 6:00 am, in SD rats and HFD animals. For mammals, although internal temperature is considered to be an inherently stationary time series with variability fixed at about 0.5°C, we reinforced this physiological condition by subtracting the mean value of the entire series from each individual measurement (13), to neutralize the action of hormonal circadian rhythms, which could lead to a certain degree of seasonality in the series (averaging). According to regulatory dynamic characteristics of the thermal control, the mechanisms of heat production and heat dissipation have to work in synergy (7).

Based on this thermal profile, we consider animal temperature as a source of information. The characterization of a source of information comes from establishment of the stochastic process. On the other hand, from the information theory point of view, every source has an associated measure quantifying its output. In our case, the associated measure is the entropy (rate). As a consequence, one of the goals in source coding is to maximize the source entropy under an appropriate set of constraints.

In this article, we focus on the characterization of the stochastic process associated with animal temperature, which will maximize a statistical measure (the entropy rate) under some constraints (the autocorrelation function) (14).

In this direction, we assume a stationary stochastic process represents animal temperature x = {X(k)}. A source is said to be stationary if given: X(1),X(2),…,X(n), and X(1 + I),X(2 + I),…,X(n + I) then, P[X(1) = x 1,X(2) = x 2,…,X(n) = x] = P[X(1 + I) = x 1,P(2 + I) = x 2,…,P(n + I) = x] for any integer n, I ≥ 1, therefore χ has a constant variance and its autocorrelation function depends only on I–(I + n).

For stationary processes, it is possible to determine the entropy rate of the source {X(k)} (14). The entropy rate of an information source is defined aswhen the limit exists. On the other hand, we also assume the source is ergodic, that is, the time average exhibited by only one realization of the source output is equal to the ensemble average with probability one. Under these conditions, the Shannon-McMillan-Breiman theorem (15) guarantees convergence of the entropy rate.

From these facts, the problem we have to solve may be stated as that of determining maximum entropy distributions, as follows.where P is the set of all probabilities or probability density functions αi are known constants, r(x) = x the i th moment of the random variable X, and p(x) denotes either the probability of the discrete random variable X or the probability density function of the continuous random variable X. The summation symbol is used either when X is a discrete random variable or when X is a continuous random variable. In this latter case, identification with the integral symbol should be clear. As a consequence, differential entropy h(x) is used instead of entropy rate H(x).

Autoregressive model

From the previous subsection, taking the earlier problem, we are faced with finding the stochastic process {X(i)} associated with animal temperature, which maximizes the entropy rate under the constraint of the autocorrelation function, that isfor k = 1,2,…,κ and all i.

The solution to this problem is the k th order Gauss-Markov process (14) of the formwhere Z(i)'s are independent, identically distributed normal random variables with zero mean and variance σ2, denoted by N(0,σ2), a1,a2,...,ak, σ2 are unknown coefficients chosen to satisfy the autocorrelation function, and k is the order of the Markov process.

In a practical problem, we usually know a sample sequence X(1),X(2),…X(n) from which the autocorrelation R(0),R(1),…,R(κ) may be calculated. Therefore, the question is: knowing the process has the form in Equation 2, is it possible to choose ak values to satisfy the constraints (Equation 1)? The answer to this question is yes. The procedure to achieve this goal, after some algebraic manipulations (in Equation 2), is given by

Equations 3 and 4 are known as Yule-Walker equations. This set of equations has a unique solution, and from them it is possible to determine the a values from the known R(k). The next question is, how many correlation lags should we consider? That is, what is the optimum value of κ? The answer to this question is to make use of the Levinson-Durbin algorithm (16,17). Formally, from the Yule-Walker equations, a system of recurrent linear equations, the aforementioned algorithm may find the coefficients. In addition to this, an adjustment is needed between the set of experimental data and the algorithm (model) used to represent it. Hence, this adjustment is achieved by use of the selection method based on the information criteria (18,19). The idea is to select memory orders smaller or larger than the true memory order in order to balance them. For instance, if we consider Equation 5, the term may be interpreted as a penalty for selecting memory orders smaller, and the term can be interpreted as a penalty for selecting memory orders larger than the true value. Therefore, minimizing these two quantities leads to the true memory order. One of these information criteria is known as the Bayesian Information Criterion (BIC). In the case of thermal time series, the memory order is determined for the model satisfyingwhere k is the memory order of the AR model, n is the sample length, and is the maximum likelihood estimator of the noise variance given by , where X(j) is the original series and is the series generated by the model and is the coefficient of the correlation from the original series.

System theory

One could ask if a system theory approach may be employed to provide the same information as that obtained from rat temperature sample sequences. In this direction, consider a discrete time system with input Z(n), impulse response H(n) and output X(n). Let Z(z), H(z) and X(z) be the corresponding z-transform, that is, Z(z) = Σ Z(n)z −, H(z) = Σ H(n)z −, and X(z) = Σ X(n)z −. From the previous subsection, we know that the stochastic process is a k th order Gauss-Markov process and, consequently, an autoregressive process. Under a system theory point of view, an autoregressive process is such that the system has only poles, that is,Equivalently,

The next step is to determine the k th coefficient of the previous polynomial equality, Equation 7. Let n + j = k in the second term on the right-hand side of Equation 7. This implies that n − j = k. As a consequence, the k th coefficient in Equation 7 is given by . Therefore,Note that this is an autoregressive system with order κ, denoted by AR(κ).

Spectral analysis

The spectral density function of an autoregressive process as described in Equation 8 may be obtained as follows: From Equation 8 we havemultiplying both sides by z−κ, and summing both sides over k and noticing that the left-hand side is the z-transform of X, we have X(z) = Σ a 1 X(k −1)z − k + … + Σ a(k − k)z − + Σ Z(k)z −. Changing variables leads to X(z) = Σ(a 1 z −)X(m)− + … + Σ a − X(m)z − + Σ(k)z −. Notice that the terms a −, for 1 ≤ j ≤ 1, are constants, yielding X(z) = (a 1 z −)(Σ X(m)z − + … + (a −)(Σ X(m)z −) + Σ Z(k)z −.

The left-hand side is the z-transform of X(.) and also of the terms in parenthesis on the right-hand side. Hence, X(z){1 − a 1 z −1 −a 2 −2 −... −a −} = Z(z).

Since E|Z(i)Z(j)| = σ2δ, it follows that where z = e. Since |X(z)|2 is just half the power spectrum, then the spectral density function is given byHence, Equation 10 is the maximum entropy spectral density subject to the constraints R(0),R(1),…,R(k).

Blood pressure measurement

Systolic blood pressure (SBP) was measured in conscious rats at 5 and 8 weeks of treatment (SD or HFD), employing an indirect tail-cuff method using an electrosphygmomanometer combined with a pneumatic pulse transducer/amplifier (BpMonWin Monitor Version 1.33, IITC Life Science, USA). This indirect approach allowed repeated measurements with a close correlation (correlation coefficient=0.975), compared to direct intra-arterial recordings (20,21). The mean of three consecutive readings represented the blood pressure.

Glucose tolerance test (GTT)

GTTs were performed in 10-week-old SD and HFD rats, after 12 h of fasting in order to determine changes in insulin sensitivity. Eleven rats from independent litters were tested. To establish basal values of glucose and insulin, blood samples were taken by lancing the tail vein before a glucose challenge (time 0). They then received a single bolus of 1 g/kg glucose ip. Blood samples were taken from the tail vein at 0, 15, 30, 60, 90, and 120 min. Glucose from whole blood was measured with a glucometer (MediSense/Optium, Abbott, USA). Serum was separated (50 μL) and kept at −20°C for measurement of insulin levels by radioimmunoassay. The incremental area under the glucose tolerance curve (AUC) was calculated as the integrated AUC above the basal value (time 0) over the 120-min sampling period using Prism 4 for Windows.

Biochemical analysis

Plasma and urine sodium and potassium concentrations were obtained at 4 weeks of treatment and measured by flame photometry (B262; Micronal, Brazil), while creatinine concentrations were determined spectrophotometrically (model 143, Instrumentation Laboratories, USA). The parameters glucose, albumin, globulin, total protein, chloride, magnesium, calcium, phosphorus, HDL, LDL cholesterol, and triglycerides were also collected at the 4th week in control and post-HFD treatment (SD and HFD) and measured using enzyme immunoassay kits with a Modular Analytic P Biochemistry Analyzer (Roche Diagnostics, Germany). Enzyme-linked immunosorbent assays for plasma insulin were used to determine rat/mouse insulin (EZRMI-13K; Linco Research, Millipore, USA), according to the manufactureŕ’s protocols.

Data presentation and statistics

All numerical results are reported as means±SD of the indicated number of experiments. Diagnostic checking was performed by analysis of the correlation of residues between the model and the samples (19,22). As shown in the Results section, the correlation of residues falls inside the boundary ±1.95/ (23). The AUC of the spectral density function was considered as the index of variability, and we have used it to establish the statistical difference between SD and HFD values. The calculi of these areas respected the same interval defined along the axis of frequency (w). Comparisons involving only two means within or between groups were carried out using the Student t-test. The level of significance was set at P≤0.05.

Results

Experimental model data

Table 1 shows the serum parameter analysis for animals from age-matched SD and HFD groups. There were no significant differences among serum biochemical parameters in SD rats compared with the age-matched HFD group, except for magnesium plasma levels, which were significantly lower in HFD compared to SD rats (P=0.0086). Arterial blood pressure (in mmHg) levels were similar in both experimental groups (SD: 137±4.24 vs 131.75±5.8 mmHg in HFD, P>0.05).

GTT

GTTs were performed to verify the effect of a high-fat diet on glucose tolerance, compared with the normal-diet group. Our study shows that, on the 4th week, the HFD group showed similar basal glycemia after overnight fasting, when compared to the SD group. Otherwise, after glucose intraperitoneal loading, the SD and HFD groups achieved similar plasma glucose concentrations at 30, 60, 90, and 120 min expressed by the incremental AUC. The measurements of plasma insulin levels during the same period (expressed in ng/mL) were also similar in SD and HFD animals.

Stochastic complexity analysis

From the previous subsections, the autoregressive model (19,22,23) may be seen as the best model to characterize temperature as the biological parameter. This model takes into consideration the inherent memory (past data) of the biological system to generate or predict future data. In mammals, adjustments in temperature are made in accordance with information observed previously. Through self-regulating mechanisms, animals are capable of avoiding significant thermal fluctuations. Estimation of the order of equations was performed via a partial autocorrelation function, , where P(k) is the autocorrelation matrix and P *(k) is the matrix P(k) with the last column being substituted by the autocorrelation vector; ||.|| denotes the determinant and the a's denote coefficients of the AR model, by use of the Levinson-Durbin algorithm. For an autoregressive process with memory κ, denoted by AR (κ), ϕ ≠ 0, for k κ, and ϕ = 0, for k κ. Note in both panels of Figure 1 that the greatest value of the partial autocorrelation function is achieved when k = 1. Best-fit equations were also chosen according to the minimum value obtained by BIC. All steps were performed with use of the software Mathematica 7.0 (16-18).

Figure 1

Partial auto-correlation function for the standard diet (panel A) compared to the high-fat diet (panel B) as a function of κ(lag).

The total number of samples was the same for all animals and achieved 361 consecutive valid measurements, covering a 6-h period from 0:00 to 6:00 am (Figure 2). These measurements were organized as follows: T = {t 1,t 2,…,t 361}, X = (X(t 1),…,X(t 361)) = (X 1,…,X 361), HFD : T = {t′ 1,…,t′ 361}, and X = (X(t′1),…,X(t′361)) = (X′1,…X′361)

Figure 2

Temperature series X for the standard diet (panel A) compared to X for the high-fat diet (panel B).

Figure 2 exemplifies the time averaging series, which were the sources used for derivation of the process for the autoregressive model. The best-fit model is achieved with the first order AR model, AR (Equation 1), given bywhere a 1 denotes the first order autoregressive coefficient, and σ2 denotes the variance of white noise. AR models derived from experimental data, Equations 12 and 13, originated from the Levinson-Durbin method from temperature samples collected from SD rats at T and and from HFD rats at T : and

Figure 3 shows the cross-correlation of residues, denoted by (the difference between correlation of the original and AR series) of the temperature sample data with those of Equation 8, up to κ = 30. As can be seen, corresponding cross-correlation values of residues κ as a function falls between the bounds: , in agreement with the model and its data. Since the coefficients in AR (Equation 1) models are less than 1 (see Equations 4, 5, 12, and 13), it follows that the associated stochastic process is stationary.

Figure 3

Representative illustration of cross-correlation of the residues for the standard diet group.

Now, if the autoregressive process has κ = 1 and a 1 = ρ, then fromby using recursively Equation 16, we arrive atFrom Equation 6, X(n) is the convolution of H(n) by Z(n), that is,By comparing Equation 18 with Equation 19, it follows that the impulse response H(n) is given byConsequently, the frequency response is

Therefore, the previous two Equations, 19 and 20, specify the time and frequency responses, respectively (see Figures 4 and 5) of the corresponding system.

Figure 4

System time response for the standard diet group.

Figure 5

Representative illustration of the system frequency response for the standard diet group.

Note that the spectral density function under a system theory point of view may be obtained as |X(z)|2 = X(x).(z), where denotes the complex conjugate of X leading to that of Equation 10. From the AR (1) model, the spectral density function is given bywhere for SD rats a 1 = 0.98 and , and for HFD a 1 = 0.96 and . The spectral density function (Equation 21), for the SD and HFD groups, is depicted in Figure 6. As can be seen, the signal energy content is located at low frequencies (39.78×10−4 rad/min) after the animals underwent the HFD for 4 weeks compared to appropriate controls. Thus, the mean AUC of this specific region to demonstrate the significant difference between the thermal profile of SD rats compared to HFD rats showed 0.04±0.0014 and 0.018±0.0028, respectively, with P<0.001.

Figure 6

Spectral density function Sw of standard diet (panels A and C) and high-fat diet rats (panels B and D).

Discussion

In the present study, we evaluated, for the first time, stochastic complexity analysis as a tool to detect, early and efficiently, presumable alterations in thermal-series profile registers in a specific model of metabolic and energy disorder. As previously defined, the temperature register obeys a nonlinear and randomized model (24). Thus, by a chaotic solution to a deterministic equation, we mean a solution whose outcome is very sensitive to initial conditions and whose evolution through phase space appears to be quite random. Taking into account the time course series of core temperature measured in control (SD for rodents) and experimental groups of animals (HFD-treated rats), as a source of information, the characterization of thermal profile comes from establishment of the stochastic process. On the other hand, from the information theory point of view, every source has an associated measure quantifying its output. In this case study, the associated temperature time course registers were treated through maximum entropy distribution as a means for stochastic characterization of the thermal series. In this way, in the current study, we focus on the characterization of thermal series maximized by a specific stochastic process (autoregressive model) under appropriate sets of constraints (the autocorrelation function) to produce the best variability analysis (spectral analysis).

The HFD model is extremely efficient for creating pathophysiological conditions such as hyperleptinemia, peripheral insulin resistance, diabetes mellitus, and obesity that lead to long-term metabolic and energy disorders. However, in this study no difference was observed in glycemic levels, blood pressure, and body mass measurement after 4 weeks of HFD intake compared to age-matched animals treated with standard chow (Table 1). On the other hand, the present study was able to establish, by autoregressive analysis, a consistent and reliable analytical method that has permitted us to characterize temperature as a referential biological parameter in homeothermic HFD-treated rats. Thus, the current temperature data analyzed by spectral density function has shown a distinct behavior when compared to the SD control and HFD-treated groups, respectively, as illustrated in Figure 6. Moreover, when statistical tests were conducted on the spectral analysis for these same series of temperature, we observed a significant decrease in the variability of values found, mainly, in the low frequency spectrum (Figure 6). This variability corresponded to one cycle every 4 h or, in the international system of notation, ∼39.78×10-4 rad/min or 0.25 cycle/h, after the animal underwent the HFD for 4 weeks, with a different circhoral rhythm to maintain the core and peripheral temperatures on small-range control in mammals (25). This reduction was observed by parametric analysis of the autoregressive model.

The role of the CNS has been demonstrated in the control of thermal and metabolic homeostasis in mammals (9,24). It is well known that thermal control pathways involve central and peripheral afferent nerve stimuli arriving into the hypothalamic preoptic area and efferent neural and humoral signaling control on thermogenesis and heat loss. Thus, studies have demonstrated that, after only 1 day of HFD intake, it is already possible to observe in rats an attenuated activity of both the insulin and leptin pathways signaling at the level of the hypothalamic nucleus (26-28). These hormones, acting either through stimulation of uncoupling proteins or via increased sympathetic activity, promote a stimulated thermogenic activity in those major organs of metabolic activity (29). Conversely, this metabolic hyperactivity leads to adjustments in the CNS to maintain stability of body temperature, necessary for the proper functioning of enzyme activities of all the chemical reactions associated with basal metabolism (30). Since this occurs in homeothermic mammals, these adjustments of body temperature are supposedly made in accordance with information data already generated and memorized. We hypothesize that it takes into account the inherent memory of biological parameters that predict data or estimate generation of future behavior profiles. In the present study, through self-regulating mechanisms, HFD-treated animals were capable of avoiding significant thermal fluctuations when compared to the SD group. Here, we may suppose that there are effective and fine-driven response sensors that are modulated by neural thermal control centers located in the hypothalamus, with regard to maintaining in a stationary condition the body temperature time series (Figure 2).

As previously demonstrated by hemodynamic studies, when analyzing heart rate time registers, reduction in variability of physiological parameters could be a marker for a deleterious prognosis in the long run (31,32), although a few studies have demonstrated some relationship between thermal profile and mortality (1,32). In the specific case of heart rate, statistical analysis of variability, by itself, was found to be very sensitive but less specific when applied in practice (1,31,32). In the current study, we accessed stochastic complexity analysis as being highly sensitive and specific for detecting discrete and subtle variations in thermal series in a specific metabolic-induced disorder. Data emerging from a stochastic analysis model of sequential temperature registers seem to be highly relevant, since there is a hierarchical priority that keeps thermal homeostasis under a narrow range of control to reach an adequate internal milieu for cells and life survival, compared with other biological parameters. We may state that mathematical analysis is an efficient tool for early prediction of the expression of pathophysiological syndromes even before appraisal of metabolic and clinical disorders. Thus, we may infer from the results of our study that thermal series analysis by a stochastic model could be proposed as a sensitive technique for clinical and pathophysiological use to predict disorders in metabolic and energy homeostasis.

Conclusion

In conclusion, maximum entropy distribution as a means for stochastic characterization of temperature time series registers may be established as an important and early tool to aid in diagnosis and prevention of metabolic diseases due to its high ability for detecting small variations in thermal profiles. Additional studies are being performed with other tools to establish the precise correlation between biological and mathematical concepts that involve representation and more specific analysis of thermal series profiles in experimental models in normal or disease states.