Is type 1 diabetes a chaotic phenomenon?

TL;DR

This study investigates whether blood glucose variations in type 1 diabetes patients exhibit chaotic behavior by analyzing nighttime glucose time series, estimating chaos indicators, and assessing predictability limits.

Contribution

It demonstrates that type 1 diabetes blood glucose dynamics may be chaotic, providing new insights into their complex behavior and potential predictability.

Findings

Type 1 diabetes blood glucose shows signs of chaos.

The maximal Lyapunov exponent indicates chaotic dynamics.

Predictability limit is approximately half of a 90-minute sleep cycle.

Abstract

A database of ten type 1 diabetes patients wearing a continuous glucose monitoring device has enabled to record their blood glucose continuous variations every minute all day long during fourteen consecutive days. These recordings represent, for each patient, a time series consisting of 1 value of glycaemia per minute during 24 hours and 14 days, i.e., 20,160 data point. Thus, while using numerical methods, these time series have been anonymously analyzed. Nevertheless, because of the stochastic inputs induced by daily activities of any human being, it has not been possible to discriminate chaos from noise. So, we have decided to keep only the 14 nights of these ten patients. Then, the determination of the time delay and embedding dimension according to the delay coordinate embedding method has allowed us to estimate for each patient the correlation dimension and the maximal Lyapunov exponent. This has led us to show that type 1 diabetes could indeed be a chaotic phenomenon. Once this result has been confirmed by the determinism test, we have computed the Lyapunov time and found that the limit of predictability of this phenomenon is nearly equal to half the 90-minutes sleep-dream cycle. We hope that our results will prove to be useful to characterize and predict blood glucose variations.

Tables

Table 1. Table 2: Data length for the 14 concatenated nights of each patient.

Patient	1	2	3	4	5	6	7	8	9	10	$⟨N⟩$
N	8243	8141	8226	6892	7387	8276	7843	8161	8199	8246	7961.4 $\pm$ 464.9

Table 2. Table 3: Fraction of false nearest neighbors as a function of the embedding dimension m 𝑚 m .

Patient	$\tau$	% fnn	$m$
1	{5, 10, 15, 20}	{0.0019 %, 0.0027 %, 0.0066 %, 0.0034 %}	{5,5,5,5}
2	{5, 10, 15, 20}	{0.0085 %, 0.0110 %, 0.0110 %, 0.0120 %}	{5,9,9,9}
3	{5, 10, 15, 20}	{0.0027 %, 0.0027 %, 0.0046 %, 0.0021 %}	{5,5,9,9}
4	{5, 10, 15, 20}	{0.0032 %, 0.0031 %, 0.0110 %, 0.0094 %}	{6,6,6,20}
5	{5, 10, 15, 20}	{0.0041 %, 0.0057 %, 0.0064 %, 0.0049 %}	{5,6,16,12}
6	{5, 10, 15, 20}	{0.0061 %, 0.0072 %, 0.0074 %, 0.0075 %}	{5,8,14,13}
7	{5, 10, 15, 20}	{0.0029 %, 0.0032 %, 0.0033 %, 0.0034 %}	{6,6,17,11}
8	{5, 10, 15, 20}	{0.0026 %, 0.0022 %, 0.0026 %, 0.0024 %}	{6,5,5,7}
9	{5, 10, 15, 20}	{0.0024 %, 0.0043 %, 0.0052 %, 0.0043 %}	{6,6,10,10}
10	{5, 10, 15, 20}	{0.0022 %, 0.0027 %, 0.0024 %, 0.0019 %}	{8,7,13,11}

Table 3. Table 4: Estimates of the correlation dimension for each patient.

Patient	$\tau$	Values of $r$ in the scaling region	Correlation Dimension $\nu$	$⟨\nu ⟩$
1	{5, 10, 15, 20}	{[5,9],[6,10],[2,9],[1,14]}	{1.72,1.71,1.71,1.74}	1.72
2	{5, 10, 15, 20}	{[16,17],[16,17],[15,16],[13,15]}	{1.19,1.22,1.19,1.17}	1.19
3	{5, 10, 15, 20}	{[10,12],[9,14],[4,7],[6,10]}	{2.42,2.66,2.69,2.65}	2.61
4	{5, 10, 15, 20}	{[16,20],[13,18],[13,17],[12,18]}	{1.32,1.34,1.40,1.43}	1.37
5	{5, 10, 15, 20}	{[18,19],[13,15],[13,18],[12,17]}	{1.36,1.39,1.36,1.37}	1.37
6	{5, 10, 15, 20}	{[13,16],[11,12],[10,12],[9,10]}	{1.38,1.47,1.48,1.47}	1.45
7	{5, 10, 15, 20}	{[15,17],[12,15],[11,14],[5,8]}	{1.89,1.91,1.90,1.89}	1.45
8	{5, 10, 15, 20}	{[8,10],[5,6],[7,8],[2,7]}	{2.19,2.16,2.23,2.18}	2.19
9	{5, 10, 15, 20}	{[12,22],[10,19],[10,12],[8,10]}	{1.78,1.79,1.79,1.81}	1.79
10	{5, 10, 15, 20}	{[12,14],[9,12],[4,8],[4,7]}	{2.26,2.25,2.27,2.30}	2.27

Table 4. Table 5: Determinism factor κ 𝜅 \kappa for the 14 concatenated nights of each patient.

Patient	1	2	3	4	5	6	7	8	9	10	$⟨\kappa ⟩$
$\kappa$	0.84	0.84	0.82	0.82	0.85	0.82	0.83	0.79	0.85	0.82	0.82 $\pm$ 0.018

Table 5. Table 6: Estimates of the ⟨ M ​ L ​ E ⟩ delimited-⟨⟩ 𝑀 𝐿 𝐸 \langle MLE\rangle and of the ⟨ L ​ T ⟩ delimited-⟨⟩ 𝐿 𝑇 \langle LT\rangle for each patient.

Patient	$\tau$	$ϵ$	MLE	$⟨MLE⟩$	$⟨LT⟩$
1	{5, 10, 15, 20}	{3,2,2.5}	{0.0149,0.0146,0.0151,0.0145}	0.0147	68’
2	{5, 10, 15, 20}	{3,2,2.5}	{0.0338,0.0340,0.0338,0.0327}	0.0336	30’
3	{5, 10, 15, 20}	{3,2,2.5}	{0.0150,0.0146,0.0149,0.0148}	0.0148	68’
4	{5, 10, 15, 20}	{3,2,2.5}	{0.0151,0.0154,0.0142,0.0146}	0.0148	68’
5	{5, 10, 15, 20}	{3,2,2.5}	{0.0326,0.0323,0.0327,0.0320}	0.0324	31’
6	{5, 10, 15, 20}	{3,2,2.5}	{0.0352,0.0356,0.0352,0.0348}	0.0353	28’
7	{5, 10, 15, 20}	{3,2,2.5}	{0.0202,0.0203,0.0189,0.0196}	0.0198	51’
8	{5, 10, 15, 20}	{3,2,2.5}	{0.0177,0.0174,0.0180,0.0175}	0.0177	57’
9	{5, 10, 15, 20}	{3,2,2.5}	{0.0161,0.0160,0.0162,0.0160}	0.0161	62’
10	{5, 10, 15, 20}	{3,2,2.5}	{0.0163,0.0164,0.0161,0.0164}	0.0162	62’

Equations

$$p\left(i\right)=\left(x_{t-(m-1)\tau },\dots ,x_{i-2\tau },x_{i-\tau },x_t\right)$$

$$I\left(\tau \right)=\sum _h\sum _k{P}_{h,k}\left(\tau \right){\log }_2\frac{P_{h,k}\left(\tau \right)}{P_h{P}_k}.$$

$$R_i=\frac{|x_{i+\tau }-x_{j+\tau }|}{||p(i)-p(j)||}$$

$$C(r)=\frac{n(r)}{\left(\genfrac{}{}{0ex}{}{N}{2}\right)}.$$

$$dist\left(p(t),p(i);\delta \right)=|x_{t+\delta }-x_{i+\delta }|.$$

$$S\left(\delta \right)=\underset{t\in {\mathcal{T}}_{\delta }}{\mathrm{mean}}\left(\ln \left\{\underset{i\in {\mathcal{U}}_t}{\mathrm{mean}}\left[dist\left(p(t),p(i);\delta \right)\right]\right\}\right).$$

Full text

I Introduction

Diabetes is a chronic disease which affects more than two hundred millions of people in the world Wile et al. (2004). There are two main types of diabetes: type 1 and type 2. In the case of type 1 diabetes, the lack of insulin due to the destruction of insulin-producing beta cells in the pancreas leads to diabetes mellitus in which “insulin is required for survival” to prevent the development of ketoacidosis, coma and death. Type 2 diabetes is characterized by disorders of insulin action and insulin secretion, including insulin resistance Alberti & Zimmet (1998). This work only concerns type 1 diabetes. Insulin-dependent diabetes or type 1 diabetes is characterized by dramatic and recurrent variations in blood glucose levels. The effects of such variations are irregular (erratic) and unpredictable hyperglycemias, sometimes involving ketoacidosis, and sometimes serious hypoglycemias. Nowadays, the development of Continuous Glucose Monitoring (CGM) devices makes it possible to record blood glucose every minute during two weeks providing endocrinologists thousands of data in the form of time series. This has also led to prediction of continuous blood glucose variations based on computational methods such as support vector machine Hamdi et al. (2018).

In the middle of the eighties, Wolf et al. [1985] proposed an algorithm allowing the “estimation of non-negative Lyapunov exponents from an experimental time series”. Thus, determination of Lyapunov exponents enabled, on the one hand, to decide whether the time series is chaotic or not and, on the other hand, to assess the Lyapunov time corresponding to the limit of predictability of the observed phenomenon. Since this algorithm was first published in Fortran code, many other versions have been developed in various languages such as C and C++ Kodba et al. (2005). Recently, this Fortran code has been implemented in MatLab by Wolf, as well as in Mathematica by Ruskeepää in 2014, see http://library.wolfram.com!infocenter/MathSource/8775/. This Mathematica software is used in the present work. The software is shortly explained in [Ruskeepää, 2015] and is also used in http://demonstrations.wolf’ram.com/ChaoticDataMaXimaILYapunovExponent/ (2017).

Long after the famous French mathematician, physicist and engineer, Henri Poincaré had discovered “deterministic chaos” in his seminal works concerning the motion of celestial bodies Poincaré (1892), many scientists searched for traces of chaotic behavior in various phenomena. The signature of chaos is the very well-known property of “Sensitive Dependence upon Initial Conditions” (SDIC) which makes the observed phenomenon unpredictable in long term. In the beginning of the sixties, Edward Norton Lorenz was the very first to identify such a feature in Meteorology Lorenz (1963). More than ten years later, Sir Robert May demonstrated the existence of “complex dynamics” (chaos) in ecological models May (1976). During the following decades, scientists highlighted chaotic behavior in human body. Indeed, according to Ives [2004] “innumerable, entwined (nonlinear) feedback loops regulate our internal processes, keeping us within the narrow bounds needed for survival. Despite this regulation, our systems are aperiodic and unpredictable in the long term.” Thus, a prime example of chaos was found in the brain Babloyantz et al. (1985) and then, in the beating of the heart di Bernardo et al. (1998); Signorini et al. (1998). Although type 1 diabetes is widely and intuitively considered by endocrinologists and clinicians as a chaotic phenomenon Goldstein & Little (1997); Berg & Sacks (2008); Little et al. (2011), this has not yet been established by numerical methods to our knowledge.

In this work, starting from a database of glucose from ten type 1 diabetes patients and while using well-known numerical algorithms with Mathematica, we give support to the conclusion that type 1 diabetes is a chaotic phenomenon and we provide the Lyapunov time corresponding to the limit of predictability of this phenomenon. These results will prove to be very useful to characterize and predict blood glucose variations. ”- This paper is organized as follows. In Section 2.1, we briefly present the three main continuous glucose monitoring devices and detail the main features of the one used in this study. In Section 2.2, we recall the method of time delay reconstruction, also referred to as delay reconstruction and phase space reconstruction Packard et al. (1980); Takens (1981), and the definitions of “time delay” and “embedding dimension”. To determine their proper values we use, respectively, “average mutual information” Fraser et al. (1986) and the method of “false nearest neighbors” proposed by Kennel et al. Kennel et al. (1992). We recall the method of correlation exponent Grassberger & Procaccia (1983) to estimate the correlation dimension and the method of local divergence rates Kantz (1994) to estimate the maximal Lyapunov exponent. In Section 3.1, we present, in some detail, all the results for the glucose data of patient 1. Then, in Section 3.2, we briefly summarize the results of all the ten patients. We also apply a direct test for determinism in a time series Kaplan et al. (1992) and the programs by Perc to state that type 1 diabetes is likely a chaotic phenomenon. Section 4 summarizes the results of the article i.e., that an estimation of the Lyapunov time which is nearly equal to half the 90-minutes sleep-dream cycle.

II Materials and Methods

II.1 Continuous Glucose Monitoring Systems

Continuous Glucose Monitoring (CGM) devices began to be developed in the eighties. However, they became available for practical use only twenty years ago with the miniaturization and development of electronic sensors combined to growing storage capacities. These systems, which have been proven to reliably reflect glucose levels Bailey et al. (2015), replace henceforth the classical finger prick blood glucose readings by monitoring interstitial fluid (ISF) glucose levels continuously. Interstitial fluid is a thin layer of fluid that surrounds the cells of the tissus below the skin Djakour -Platonoff et al. (2003); Buckingham (2008); Benhamou et al. (2012). That’s the reason why there is a 5 to 10 minute delay in interstitial fluid glucose response to changes in blood glucose. This result of great importance will need to be compared to the Lyapunov time obtained in this study (see Sec. 4). Today, the three main manufacturers which propose devices with continuous glucose monitoring reading are Abbott (Freestyle Navigator II), Medtronic (MiniMed coupled with Veo-pump), and Novalab which offers the reader Dexcom G4 coupled to the insulin pump Animas Vibe. Whatever the system, it is composed of two parts. The glucose sensor which has a small, flexible tip that, inserted just under the skin, continuously measures the glucose concentration in the interstitial fluid, and stores data during several days. The blood glucose reader is used to scan the sensor and displays the current glucose reading based on the most recently updated glucose value during the latest hours of continuous glucose data. Some readers are also coupled to an external insulin pump (Paradigm VEOTM Medtronic or Animas VibeTM of Novalab). In this work, we have chosen to use the Abbott (Freestyle Navigator II) system because the blood glucose reader records the blood glucose variations continuously every minute all day long during fourteen consecutive days. This represents, for each patient, one value of glycaemia per minute during 24 hours and 14 days, i.e., 20,160 data. Thus, ten type 1 diabetes patients have accepted to provide us the recordings of their blood glucose during fourteen consecutive days so that they could be anonymously analyzed.

II.2 Methods

Following the works of Takens [1981], Sauer et al. [1991], and Abarbanel [1996], summarized in Kodba et al. [2005], we consider the reconstruction of the attractor in an $m$-dimensional phase space starting from the time series $\{x_{1},\ldots,x_{i},\ldots,x_{T}\}$ of blood glucose variations for each patient. Here $x_{i}$ denotes the glycaemia in $i^{th}$ minute. According to Takens [1981], the reconstructed attractor of the original system is given by the vector sequence

[TABLE]

where $\tau$ and $m$ are the time delay and the embedding dimension, respectively. Takens’ famous theorem states that for a large enough $m$, this procedure provides a one-to-one image of the original system. It follows that the attractor constructed according to equation (1) will have the same dimension and Lyapunov exponents as the original system. To reconstruct the attractor successfully, pertinent values of $\tau$ and $m$ have to be accurately determined.

II.2.1 Time Delay

Two criteria are to be taken into account for the estimation of the time delay $\tau$:

$\tau$ has to be large enough because the information we get from measuring the value of $x$ at time $i+\tau$ should be significantly different from $x$ at time $i$.

• -

$\tau$ must not be larger than the time in which the system loses memory of its initial state. This is important for chaotic systems, which are unpredictable and lose memory of the initial state.

Fraser et al. [1986] defined the mutual information between $x_{i}$ and $x_{i+\tau}$ as a suitable quantity for determining $\tau$. The mutual information between $x_{i}$ and $x_{i+\tau}$ measures the quantity of information according to the following expression

[TABLE]

Here $P_{h}$ and $P_{k}$ denote the estimated probabilities that supposes a value inside the $h^{th}$ bin and the $k^{th}$ bin, respectively and $P_{h,k}(\tau)$ is the estimated probability that $x_{i}$ is in $h^{th}$ bin and $x_{i+\tau}$ is in $k^{th}$ bin. The optimal choice for the time delay is given by the first minimum of $I(\tau)$, because $x_{i+\tau}$, then adds the largest quantity of information to the information at $x_{i}$.

II.2.2 Embedding Dimension

The method of false nearest neighbors proposed by Kennel et al. [1992] is used to determine the embedding dimension $m$. This method is based on the fact that points which are close in the reconstructed embedding space have to stay sufficiently close also during forward iteration. So the distance between two points of the reconstructed attractor cannot grow further as a threshold $R_{tr}$. Nevertheless, if an $i^{th}$ point has a close neighbor that does not fulfil this criterion, then this $i^{th}$ point is marked as having a false nearest neighbor. In order to calculate the false nearest neighbors, the following algorithm is used:

For a point $p(i)$ in the embedding space, we have to find a neighbor $p(j)$ for which $||p(i)-p(j)||<\varepsilon$, where $||\ldots||$ is the square norm and $\varepsilon$ is a small constant usually not larger than the standard deviation of data.

• -

The normalized distance $R_{i}$ between the points $p(i)$ and $p(j)$ is computed:

[TABLE]

• -

If $R_{i}$ is larger than the threshold $R_{tr}$, then $p(i)$ has a false nearest neighbor. According to Abarbanel et al. [1996], $R_{tr}=15$ has proven to be a good choice for most data sets.

After having computed (with the algorithms quoted above) the optimal time delay $\tau$ and the embedding dimension $m$, the attractor can be successfully reconstructed.

II.2.3 Correlation Dimension

The correlation dimension is a generalization of the usual integer-valued dimension. The method used is called the method of correlation exponent. Suppose that we have data $x_{1},x_{2},\ldots$. Prepare the $m$-dimensional delay coordinates $p(t)=(x_{t-(m-1)\tau},\ldots,x_{t-2\tau},x_{t-\tau},x_{t})$. Let $N$ be the number of these delayed points $p(t)$. Define $n(r)$ as the number of pairs $(p(i),p(j))$, $i<j$, whose distance is smaller than $r$. Define then $C(r)$ to be the corresponding relative frequency:

[TABLE]

Thus, $C(r)$, also called the correlation sum, represents the probability that a randomly chosen pair of points in the reconstructed phase space will be less than a distance $r$ apart. When $N$ approaches infinity, $C(r)$ is called the correlation integral. It is a measure of the spatial correlation between data points. It can be shown that when $N$ is large and $r$ small, $C(r)$ behaves approximately like a power function: $C(r)\approx\alpha r^{\nu}$, that is, $C(r)$ is proportional to $r^{\nu}$. Here, $\nu$ is called the correlation exponent. Because $\log C(r)\approx\log a+\nu\log r$, we see that $\log C(r)$ behaves linearly as a function of $\log r$. Furthermore, $\nu$ is the slope of $\log C(r)$ versus $\log r$. To estimate the correlation dimension, we estimate $\nu$ for increasing $m$. If the correlation exponent saturates with increasing $m$, the system is considered to be chaotic and the saturation value is called the correlation dimension of the attractor. Thus, in estimating the correlation dimension, we look for an interval for $r$ where $\log C(r)$ versus $\log r$ behaves approximately linearly. Such an interval is called a scaling region. If the slopes of the curves in the scaling region are approximately constant for a range of values of $m$, this constant is an estimate of the correlation dimension. The saturation of the correlation exponent for increasing $m$ is a necessary condition for accepting an estimate as the correlation dimension Galka (2000). For chaotic data, the correlation dimension gives a fractional dimension for the strange attractor Grassberger & Procaccia (1983). Note that the presence of noise in the data makes the estimation of the correlation dimension much more uncertain. Noise may increase the estimate and may cause nonconvergence of the estimate for increasing $m$ Galka (2000). Define the noise level to be the ratio of the standard deviation of the noise to the standard deviation of the noise-free data. A small noise level of 2% will typically make the estimation of the correlation dimension impossible! Thus, nonlinear noise reduction is important to noisy data. For a stochastic system (in contrast to a chaotic system), the correlation exponent increases without bound when the embedding dimension is increased. Thus, the correlation exponent can be used to distinguish between low-dimensional dynamics and randomness. In the calculation of the correlation sum, we encounter one problem, namely dynamic correlation. To avoid the dynamic correlation, only pairs of points whose time distance is large enough are used. To this aim, Provenzale et al. [1992] have proposed the use of so-called space-time-separation plot. Once such time distance has been evaluated, the correlation sums are computed and so, the correlation dimension can be estimated.

II.2.4 Determinism Test

In order to confirm that type 1 diabetes is a chaotic phenomenon we have used the determinism test proposed by Kaplan et al. [1992]. The determinism test is of great important because it enables us to distinguish between deterministic chaos and irregular random behavior, which often resembles chaos. A program by Perc [2006] enables to compute the determinism factor $\kappa$ which is substantially smaller than 1 for a system with a stochastic component while it is close to 1 for a system with a chaotic component. According to Kodba et al. [2005] “the deterministic signature is still good enough to be preserved, so that the chaotic appearance of the reconstructed attractor cannot be attributed to stochastic influences.”

II.2.5 Maximal Lyapunov Exponent

The maximal Lyapunov exponent is a characteristic of the dynamical system and quantifies the strength of chaos. Indeed, in chaotic systems nearby trajectories diverge exponentially fast. The properly averaged exponent of this divergence is called the maximal Lyapunov exponent. It describes the average rate at which predictability of the system is lost. More than twenty years ago, Kantz [1994] presented the method of local divergence rates to calculate the maximal Lyapunov exponent (see also Kantz and Schreiber [2004]). This method estimates local divergence rates of nearby trajectories. Let the data be $x_{1},x_{2},\ldots,x_{T}$. Choose a delay time $\tau$ and an embedding dimension $m\geqslant 2$ and prepare the $m$-dimensional delay coordinates $p(t)=(x_{t-(m-1)\tau},\ldots,x_{t-2\tau},x_{t-\tau},x_{t})$, $t=1+(m-1)\tau,\ldots,T$. Choose an $\epsilon>0$, to be used later to find points of the data in an $\epsilon$-neighborhood. Choose an integer $\delta$ (typically, $\delta$ is between 0 and, say, 20). Choose a $t$ from the set $\mathcal{T}_{\delta}=\{1+(m-1)\tau,\ldots,T-\delta\}$. Find all points $p(i)$ that are in the $\epsilon$-neighborhood of $p(t)$; let $\mathcal{U}_{t}$ be the set of these indices $i$. Calculate, for all $i\in\mathcal{U}_{t}$, the distance between $p(t)$ and $p(i)$ after the relative time $\delta$; this distance is defined as

[TABLE]

Thus, the distance is the absolute value of the difference between the scalars $x_{t+\delta}$ and $x_{i+\delta}$. Calculate then the logarithm of the mean of these distances. Do this for all $t\in\mathcal{T}_{\delta}$, and then calculate the mean of the calculated values, that is, the $\mathcal{T}_{\delta}$-mean of the logarithms of the $\mathcal{U}_{t}$-means of the distances. In this way, we get a value that we denote $S(\delta)$. A formula for $S(\delta)$ is as follows:

[TABLE]

It is a measure of the local divergence rate of nearby trajectories. Repeat these calculations for various values of $\delta$, such as $\delta=0,\ldots,20$. In this way we get the values of $S$ at discrete points. For chaotic dynamical systems, the $S(\delta)$ function behaves as follows:

For small values of $\delta$, $S(\delta)$ may behave a little irregularly.

• -

For intermediate values of $\delta$, $S(\delta)$ increases linearly. Such values of $\delta$ define a so-called scaling region. The slope of the linear growth is an estimate of the maximal Lyapunov exponent. A positive Lyapunov exponent is an indicator of chaotic behavior.

• -

For larger value of $\delta$, $S(\delta)$ tends towards a constant.

In the following section we have used the programs mentioned above for the determination of the time delay $\tau$, the embedding dimension $m$ and of course for the estimation of the correlation dimension, the maximal Lyapunov exponent as well as the Lyapunov time, i.e., the inverse of the maximal Lyapunov exponent.

III Results

We applied the method of reconstruction of the attractor in an $m$-dimensional space phase (recalled in the previous section) to the time series of blood glucose variations for each of our ten patients. In the first place, we considered the whole time series including days and nights. Nevertheless, because of the stochastic inputs induced by daily activities of any human being (work, eat, walk, sport,…), it has not been possible to discriminate chaos from noise. So, in the second place, we have decided to keep only the nights of these ten patients starting from 10:00 p.m. to 08:00 a.m. Therefore, for each patient, we have one value of glycaemia per minute during 10 hours, i.e., 600 data. But, numerical analysis methods imply the time series to have a sufficient length, typically comprising several thousand points or more. Consequently, we have concatenated the 14 nights of each patient in order to obtain for each of them nearly 8,400 data. Then, we applied the programs developed by Ruskeepää and Perc to each of our ten patients. In the first subsection, we detail the procedure for patient 1. In subsection 3.2 we will only give the results for all the ten patients.

III.1 Patient 1

III.1.1 Blood Glucose Data

For the 14 concatenated nights of this first patient, starting from 10:00 p.m. to 08:00 a.m., we have $N=8243$ data (see Fig. 1).

III.1.2 Time Delay

According to subsection 2.2.1, the optimal choice for the time delay is given by the first minimum of the average mutual information $I(\tau)$. For our data, there is no point of minimum (see Fig. 2). So, in this case, we can consider the value of $\tau$ at which the first substantial decrease of the average mutual information stops and a slower decrease begins. Thus, we can choose for the time delay say the value $\tau=5$. We could also have chosen as the time delay another value like $10,15$, or $20$. Consequently, we will use these four values of the time delay in the following to show that the results of our computations are invariant with respect to $\tau$.

III.1.3 Embedding Dimension

According to subsection 2.2.2, we use the method of false nearest neighbor to compute the embedding dimension. However, since we consider four possible values of the time delay, we have plotted in Fig. 3 the fraction of false nearest neighbors for $\tau=5,10,15,20$ as a function of the embedding dimension $m$ (the lowest curve is for $\tau=5$ and the highest curve for $\tau=20$). It appears that the fraction of false nearest neighbors does not drop to zero with any value of $m$ (whatever the value of $\tau$) but with $m=5$ the fraction is almost zero.

These results are summarized in Table I. So, for patient 1 we can choose the embedding dimension $m=5$.

III.1.4 Correlation Dimension

According to subsection 2.2.3, we use the method of correlation exponent. First, to avoid the dynamic correlation we use space-time-separation plot proposed by Provenzale et al. [1992]. The plot on Fig. 4 indicates that if points are at least, say, 500 steps from each other, they will no more have dynamic correlation.

We calculate correlation sums by first using $\tau=5$, $m=1,\ldots,12$, minimum time distant 500, and space distances $r=10^{d}$, where $d$ ranges from 0 to 2.5 so that we have a total of 25 values of $d$ (see Fig. 5).

It turns out that a scaling region is between the fifth and ninth value of $r$: there (between the two vertical lines, that is, for $r$ from somewhat below 3 to somewhat below 7) the $C(r)$-curves in Fig. 5 behave approximately linearly and, as Fig. 6 shows, the slopes of the $C(r)$-curves in the scaling region approach a saturation value, approximately 1.72. Indeed, for $m=5,6,7$, and 8 the slopes are near to each other: 1.717, 1.720,1.723, and 1.725. The mean of these slopes is 1.721. This leads, for $\tau=5$, to an estimate of 1.721 for the correlation dimension of the attractor.

These calculations for the correlation dimension was done by first using a Mathematica program to calculate the correlation sums (for a fixed value of $\tau$) and then applying a dynamic Mathematica interface to find a scaling region and the corresponding estimate of the correlation dimension. This dynamic interface contains sliders to manually adjust the start and end points of the possible scaling region and shows the slopes of the correlation sums, dynamically reflecting the current choice of the scaling region. So, to find a scaling region the sliders are adjusted until a promising region is found where the correlation sums behave approximately linearly and the slopes in the scaling region saturate with increasing $m$. If, in addition, the invariance with respect to $\tau$ can be verified, this gives support to the estimated correlation dimension.

III.1.5 Determinism Test

While using Perc’s program an estimation of the determinism factor is $\kappa=0.84$.

III.1.6 Maximal Lyapunov Exponent

To estimate the maximal Lyapunov exponent, we use the method of local divergence rates. First, we choose the delay time $\tau=5$ and set $\epsilon=3$. The embedding dimension $m$ varies from 1 to 20. We obtain the values of $S(\delta)$ plotted in Fig. 7.

It turns out that a scaling region is for $\delta=1,\dots,6$: there the $S(\delta)$-curves in Fig. 7 behave approximately linearly (for larger values of $m$) and, as Fig. 8 shows, the slopes of the $S(\delta)$-curves in the scaling region are approximately a constant for $m=12,\ldots,17$; the mean of these slopes is 0.0150. This leads, for $\tau=5$, to an estimate of 0.0150 for the maximal Lyapunov exponent.

The maximal Lyapunov exponent is invariant under the embedding procedure and also under different values of $\epsilon$. Fig. 8 shows, for $\tau=5$ and $\epsilon=3$, the invariance with respect to the embedding dimension $m$ (the slopes are approximately constant for $m=12,\ldots,17$). For $\tau=5$, we also considered the values $\epsilon=2.5$ and 2, and we got the estimates 0.0141 and 0.0156. The estimates for the three values of $\epsilon$ are near to each other, supporting, for $\tau=5$, the invariance with respect to $\epsilon$. The mean of the three estimates is 0.0149.

We also estimated the maximal Lyapunov exponent for $\tau=10,15$, and 20. For each $\tau$, we considered the values 3, 2.5, and 2 for $\epsilon$. For each $\tau$, the estimate was invariant with respect to $\epsilon$. The means of the estimates for $\tau=5,10,15$, and 20, were 0.0149, 0.0146, 0.0151, and 0.0145. These values are near to each other, supporting the invariance with respect to $\tau$. The mean of the four estimates is 0.0147 (here, we used the original, more accurate estimates). Thus, the final estimate of the maximal Lyapunov exponent is 0.0147. The Lyapunov time is then $1/0.0147\approx 68minutes$.

These calculations for the maximal Lyapunov exponent done by first using a Mathematica program to calculate the values of $S(\delta)$ (for a fixed value of $\tau$ and $\epsilon$) and then applying a dynamic Mathematica interface to find a scaling region and the corresponding estimate of the maximal Lyapunov exponent.

In the next subsection, we summarize the main results for all the patients.

III.2 Patient 1 to 10

III.2.1 Data Length

The lengths of the data for the ten patients (see Table II) varied between 6892 and 8276 with a mean of 7961.4. For each patient, the data contains values of glucose for 14 concatenated nights. The lengths of the data vary somewhat because some values of the glucose have not been recorded by the glucose sensor. The lengths of the data turned out to be large enough for the numerical programs to work properly.

III.2.2 Average Mutual Information: determination of the time delay $\tau$

For the time delay $\tau$, the situation for all the ten patients was similar: the average mutual information did not have a local point of minimum. As examples of the time delay, we considered the values $\tau=5,10,15$, and 20 when we estimated the embedding dimension, the correlation dimension, and the maximal Lyapunov exponent.

III.2.3 False Nearest Neighbors: determination of the embedding dimension $m$

The embedding dimension $m$, estimated by the method of false nearest neighbors, varied between 5 and 20 for the four values of $\tau$. The fraction of false nearest neighbors did not drop exactly to zero but decreased to a small enough level (typically below 0.01 but above 0.001). The results are presented in Table III.

III.2.4 Correlation Dimension

The correlation dimension (see Table IV), estimated by the method of correlation exponent, varied between 1.19 and 2.61, the mean being 1.74. For each patient, the correlation dimension was estimated with the four values of the time delay $\tau$ (the values of the correlation dimension given in Table IV are the means of the four estimates got by using the four values of $\tau$).

III.2.5 Determinism Test

The estimation of the determinism factor $\kappa$ for the 14 concatenated nights of each of the 10 patients is presented in Table V.

Another program developed by BenSaïda [2015] has enabled to confirm the presence of a chaotic dynamics, as opposed to stochastic dynamics, in the 14 concatenated nights of each of our 10 patients.

III.2.6 Maximal Lyapunov Exponent and Lyapunov Time

The maximal Lyapunov exponent (see Table VI), estimated by the method of local divergence rates, varied between 0.0147 and 0.0353, the mean being 0.0215. For each patient, the exponent was estimated with the four values of $\tau$ and, for each $\tau$, the exponent was estimated with $\epsilon$ equal to 2, 2.5, and 3 (the values given in Table VI are the means of the resulting 12 estimates). The Lyapunov time (see Table VI) varied between 28 and 68 minutes, the mean being 52.5 minutes. The inverse of the final estimate of the Maximal Lyapunov exponent, 0.0215, is 46.5 minutes.

IV Discussion

Starting from a database often type 1 diabetes patients wearing a continuous glucose monitoring device, we have used numerical methods to analyze the continuous variations of the blood glucose every minute all day long during fourteen consecutive days. Nevertheless, because of the stochastic inputs induced by daily activities of any human being (work, eat, walk, sport, etc.), it has not been possible to discriminate chaos from noise in the whole time series including days and nights. So, we decided to keep only the fourteen nights of these ten recordings of the blood glucose, from 10:00 p.m. to 08:00 a.m. We then concatenated these recordings over fourteen nights and so we obtained, for each patient, nearly 8,400 data values; this amount of data was large enough for the computations.

We used the average mutual information to investigate the time delay $\tau$ for each of our ten patients. As the average mutual information did not have a local point of minimum, in later computations we considered, as examples, the values 5, 10, 15, and 20 for $\tau$. We were able to show that the estimates of the correlation dimension and the maximal Lyapunov exponent were invariant with respect to these values of $\tau$. The method of false nearest neighbors led to a determination of the embedding dimension $m$, a typical value being between, say, 5 and 12, the value varying with $\tau$ and with the patient in question.

The test of determinism provided a mean value 0.83 of the determinism factor for the ten patients. This value is close to one and so supports the chaoticity of type 1 diabetes.

With the four values of $\tau$, we applied the method of correlation exponent to the ten time series. In each case, the correlation sum saturated with increasing $m$ and was then approximately constant in a scaling region, that is, for some values of $m$. This behavior supports the chaoticity of our data. The correlation dimension was, on the average, 1.74.

Finally, we estimated the maximal Lyapunov exponent with the method of local divergence rates. For each of the ten patients, we estimated the exponent for the four values 5, 10, 15, and 20 of $\tau$ and, for each $\tau$, for the three values 2, 2.5, and 3 of $\epsilon$. The estimates turned out to be invariant with respect to $\tau$ and $\epsilon$. For each patient, the mean of these twelve estimates is given in Table VI. The overall mean of these ten estimates is 0.0215. The corresponding estimate of the Lyapunov time is 46.5 minutes.

According to Feinberg and Floyd [1979], “The notion that there exists in man a basic rest-activity cycle with a period of 90 min has been accepted by many investigators. This hypothesis was introduced by Kleitman (1963), who based it upon early observations (e.g. Dement & Kleitman, 1957) that REM sleep in adults recurred at about 90-min intervals.” More than ten years before this publication. Hartman [1968] proved that “the mean cycle length for 15 normal adult subjects in our laboratory, studied for eight or more nights each, was 95.8 minutes, with a standard deviation of 8.7 minutes.” Thus, the “sleep-dream”, also called rapid eye movement (REM)-sleep, recurred in adults at about 90 minutes intervals and is linked to real time, i.e., to the circadian cycle. Since we have found an estimate of the Lypaunov time equal to 46.5 minutes, which corresponds nearly to half the sleep cycle or rest-activity cycle, it seems to indicate that some correlations could be highlighted between blood glucose variations in type 1 diabetes patients and their REM-sleep. This would be subject of a new protocol of research in which type 1 diabetes patients wearing a continuous glucose monitoring device will undergo a polysomnography.

Acknowledgements.

Matjaž Perc was supported by the Slovenian Research Agency (Grants J1-7009 and P5-0027).

V Figures

VI Tables