Fractal Calculus Facilitates Rethinking ‘Hard Problems’: A New Research Paradigm

Abstract

This paper introduces a non-standard research technique to clarify how complex phenomena, such as those that are abundantly present in human physiology, can be faithfully described using fractal dynamical models with and without stochastic forces. This method for conducting research involves tracing the historical evolution of understanding an empirical medical process facilitated by the fractal-order calculus perspective. Herein, we trace the analysis of the time series for heart rate variability (HRV) developed for diagnosing the cardiovascular health of a patient. This is performed herein by introducing four (one empirical, which entails three theoretical fractal models) distinct but related fractal models, each one introduced to solve a particular problem arising from a fundamental defect in the previous model, but in generalizing a model at one stage to resolve the problem associated with the defect, another is invariably introduced by the replacement model. It is through the utilization of the fractal-order calculus that the necessity for rethinking how to systematically incorporate additional layers of complexity is revealed, ultimately resulting in a ‘complete’ description of its empirical dynamics in fractal terms.

1. Introduction

It is rare that the first pass at constructing the mathematical model of a complex process is error-free, which is the two-fold reflection of the difficulty associated with matching a sufficiently general mathematical infrastructure to the specific behavior of the complex phenomenon of interest. However, even after satisfactorily completing the construction of a fractal-order calculus (the terminology ‘fractal calculus’ is used throughout to emphasize the fact that the fractal descriptor means that ‘fractal geometry’ and ‘fractal statistics’ are vastly different from the usual terms of geometry and statistics. So too is the ‘fractal-order calculus’ (FOC) remarkably different from the integer-order calculus (IOC) and the integer-order differential equation (IODE)) dynamical model (FOCDM) of a complex phenomenon, an empirical property of that phenomenon may reveal itself to be incompatible with the FOCDM thus constructed. In response to the new information, typically one has to decide to either generate a new model from scratch or to generalize the existing FOCDM to incorporate the new information into the potential behavior modes of the existing FOCDM. Our intention here is to create a modeling strategy based on the second option, which is sufficiently robust that an empirical correction can be correctly handled through one or more modest tweaks of the initial FOCDM.

Herein, we do not perform any statistical processing of datasets but rather search elsewhere for explanations of any discrepancy between the dynamic models constructed in the literature and the making of reasonable medical interpretations of the statistical datasets. More specifically, rather than assuming the observed statistical discrepancies are due to data processing errors, we look for a deterministic source for the errors since the source of error lies either in the representation of those statistics or in the modeling of the physiological control of those statistics. We choose the latter route and accept without criticism the manner in which statistics have been treated in characterizing a ‘hard problem’ in physiology. The fractal-order calculus is used to construct and apply new models to examine the behavior and control of HRV datasets.

1.1. Fractal-Paradigm and the Fractal-Order Calculus

Herein, we address concerns regarding how complex phenomena, such as those that are apparently ubiquitous in physiology, can be faithfully described with stochastic nonlinear dynamic models using stochastic nonlinear FOCDM, or FOCSNDM if you wish. Since the closing decade of the last century, the anatomist/physiologist Ewald Weibel had been working to develop fractal geometry as a biological design principle by means of its scaling behavior to provide insights into the possibilities of the efficient genetic programming of the biological form [1]. He ultimately synthesized over two decades of research in his arguments on the design of biological fractal geometry into the general symmorphosis hypothesis, a name that he and C.E. Taylor had invented in 1981. The Symmorphosis, in Weibel’s book of the same name [2], posited that biological systems adhere to an ‘economy of design’ thereby pairwise matching their various structural and functional parameters. Over twelve years of the same period many of the empirical results obtained in this area of research were communicated in a series of workshops, using bottom-up theoretical FC approaches, and published in the four volumes: Factals in Biology and Medicine I [3], II [4], III [5] and IV [6].

One aspect of the fractal paradigm concerns itself with the 1/f-behavior of the self-similar statistics of a fractal time series [7,8] in which this 1/f-variability arises because the underlying process does not vary in “chronological time” but in an “intrinsic time”, which is a form of fractal time [9,10,11,12]. The fractal time has to do with the fact that the events are crucial, and therefore, in addition to their being renewal [13,14], the time interval statistics are IPL with an IPL index within the interval 1 < $\mu$ < 3 (see Section 7.1 for an overview). The fractal time suggests that the description of the evolution of the PDF is determined by the fractal-order probability calculus; see, e.g., [12,15,16].

A sketch of the formal solution to a Fractal–Order Kinetic Equation (FOKE) for the scaling PDF is given in Section 7.2 using renormalization group theory. The equation to be solved was first derived using chaos theory and the FOC by Zaslavsky [17] and subsequently rederived [18] using the continuous-time random walk model of Montroll and Weiss [19]. Consequently, the FOKE describing the space-time evolution for the scaling PDF was determined to be:

$$D_t^{\alpha }\left[P\left(x,t\right)\right]=K_{\beta }{D}_{\left|x\right|}^{\beta }\left[P(x,t)\right],$$

(1)

where the Caputo fractal-order derivative in time is operating on the left side of the equation, and the Riesz–Feller fractal-order derivative in space is operating on the right side of the equation. The FOKE is solved as an initial value problem for the FOPD, where $P_0\left(x\right)=P(x,t=0)=\delta \left(x\right);$ see Section 7.2.

Consequently, using the FOC approach to theory development is top-down rather than bottom-up, which is achieved by considering the high points in the historical record [20] of the relatively well-understood physiologic time series that constitutes heart rate variability (HRV), developed for diagnosing a patient’s cardiovascular health. The value of FOC in modeling such physiologic time series is displayed herein by examining the structure of four distinct but related nonlinear stochastic FOCDMs of HRV. Each model is introduced in turn to overcome a particular deficiency in a previous model, but resolving the defect at a given stage inadvertently introduces another defect that requires resolution to advance to any subsequent stage. It is through the utilization of the FOC that the necessity for rethinking how to systematically incorporate additional layers of interrelated complexity is revealed, ultimately resulting in a more complete understanding of HRV dynamics in general, both in terms of the mathematics used and the interpretation of those mathematics for bedside medicine.

The method we are proposing is an add-on to a familiar one in that it involves identifying the behavioral path taken by a complex phenomenon undergoing failure and requires intervention to reestablish health. Documenting the cause of a particular failure and determining how a more accurate modeling of the component that failed could have saved the network of which the component was a part by an appropriate intervention is the ’Holy Grail’ of every discipline whose purpose is related to the health and well-being of human beings and the groups to which they belong. Without question, the science and practice of medicine have developed this documentation procedure for the understanding of disease and injury into a fine art and have extended it into the establishment of protocols for new methods of rehabilitation as well. What we introduce into this time-proven research paradigm [20] is FOC mathematics for the proper description of what is known but which also has the flexibility to anticipate which of its FOCDM components has an increasing probability of failure and how to keep that FOCDM component failure from occurring, and failing that, to minimize the damage its failure entails and to do this all within the same FOC modeling paradigm. Therefore, let us begin by clarifying what is meant by a basic research paradigm.

Basic Research Paradigm

A basic research paradigm is a system of fundamental beliefs and agreements shared among scientists about how problems should be understood and addressed. Its framework guides the entire research process, including the choice of methods, and influences how researchers interpret their results. A paradigm can also be thought of as a worldview, involving shared understandings of reality detailing the philosophical underpinnings of a researcher’s approach and helping to improve the rigor of their projects and communication in research. For these reasons, any move to adopt new research paradigms ought to be treated with more than the usual scientific skepticism. Rahman [21] provides “An overview and critique” of the different research paradigms available to an investigator. Our main interest is in the practical circumstances requiring that a paradigm be replaced without loss of what was gained with the initial and subsequent FOCDM replacements along the lines of [22] and/or [23].

It is true that when new experimental results consistently contradict established theory, the replacement theory must satisfactorily explain both the pre-existing as well as the new experimental results, and this may require abandoning the existing research paradigm. The necessity of quantum mechanics to resolve the empirically-based wave/particle-duality paradox changed our view of the physical universe as well as how the world of microbiology can dominate disease through such elusive things as molecules, germs and microbes, being extreme examples of what may occur. Einstein’s theory of relativity was no less fundamental in its influence on physics but its influence on our understanding of the everyday world in which we live only affects those aspects outside our daily experience, involving things whose speed approaches that of light.

The kind of paradigmatic change we present herein typically lacks the drama associated with these examples in that our purpose is to examine and process empirical data associated with the failure of a network to better understand the who, what, where, why and when (5 Ws) of the failure. The author maintains that this was not typically performed in the past and has only emerged piecemeal in the study of failure modes in complex systems, for example, through the joining together of complexity theory with network science [24,25,26].

Herein, we offer a practical strategy for using what is learned in the early stages of a 5 W empirical history describing the failure of an organ network (ON) within a network-of-organ-networks (NoONs) to formulate ways, in the most extreme case, to keep the NoONs from failing when the ON of interest fails, e.g., by replacing the failing ON with a viable synthetic one. The exemplar used herein concerns the human body as the NoON and the cardiovascular system as the failing ON, as determined using HRV time series. However, the form of what is revealed is not discipline-specific and could readily be adapted to, for example, a failing organization-network within a society of such interacting and interdependent networks; see. e.g., [25].

1.2. Brief History of HRV as a ‘Hard Problem’

The history of using the time intervals between successive heartbeats to monitor an individual’s state of health is a long and interesting one [27]. Heart rate variability (HRV) is the variation in rate between a succession of RR intervals, which has become a popular clinical and investigational tool since its 1996 endorsement by the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology [28]. Over the past half-century, one of the central topics in physiological signal analysis has been the characterization of HRV time series, which, when properly processed, has evolved into a vital non-invasive indicator of the performance of the cardiovascular system. It is well established that the cardiovascular system forms a complex structure in which nonlinear control mechanisms shape the cardiac activity on the scales of milliseconds, seconds, minutes, hours, and days [29,30]. The decrease in both statistical and spectral indices of HRV has been linked to poor outcomes in patients with, for example, myocardial infarction [31], coronary artery disease [32] and heart failure [33], while distinct peaks in a HRV power spectral density (PSD) have been attributed to autonomic, endocrine and vagal modulation [34].

The spectral analysis of HRV time series has, since its acceptance in 1996, revealed the complexity of cardiovascular regulation in the non-normal statistical nature of heart-rate time series [35]. The fractal fluctuations in HRV time series entail important physiologic advantages associated with the ability of the cardiovascular system to respond to unpredicted stresses and stimuli, since nonlinear responsiveness offers greater flexibility than linear reactivity [36].

One of the first studies focusing on the scaling properties of HRV time series was conducted by Peng et al. [34] and revealed that the PSD of the time increments between successive beats scale as an inverse power law (IPL) in frequency:

$$S_p\left(f\right)\propto f^{-\mu }.$$

(2)

The scale-free behavior of the PSD implies that the HRV time series lacks a characteristic time scale, which has been interpreted to mean that the fluctuations are not completely random but reflect the influence of various processes occurrring at different time scales in the control and regulation of the cardiovascular ON.

The fact that so much has been published concerning the empirical evidence addressing HRV’s diagnostic utility is a clear indication of the need to understand these limitations even outside the HRV area of study. These limitations are the boundaries of the applications of a useful model and are interesting to even the larger medical community since such limitations may provide insight into potential limitations within other distinct research areas and suggest an expansion or contraction of their own boundaries. Investigators may be seeking better processing techniques, where by better we mean capable of achieving the same level of accuracy using only a fraction of the data needed with today’s methods, or perhaps a more sophisticated method for localizing the peak of the R component to the QRS-complex [27], the source of the HRV time series. These things are only briefly mentioned to emphasize that the final interpretation of HRV time series has not been reached, and it remains a ‘hard problem’ worthy of reexamination. Chalmers [37] coined the term ‘hard problem of consciousness’ to distinguish the totality of consciousness from the easy problems that are amenable to reductive logic. Herein, we classify nearly every fundamental problem remaining in science as a ’hard’ problem in the sense that they typically arise out of empirical paradox [38]. How successful we are in justifying the utilization of fractal calculus as the facilitator for this reexamination, the reader will have to judge.

1.3. Outline of Modeling Synthesis Using 4 Models

In Section 2, the first analysis of the HRV time series conducted by Peng et al. [34], which revealed the statistical fluctuations in heart-rate data to have a Lévy stable probability density function (PDF), is reviewed, yielding the zeroth FOC model (FOCM-0), which is based on the processing of observed HRV time series. In Section 3, a theoretical derivation of the empirical Lévy PDF as the first-order FOC model (FOCM-1) using a fractal-order Fokker–Planck equation (FOFPE) is reviewed, as are the reasons why such a strict empirical result is improbable. We introduce recent empirical evidence in Section 4, explaining how the severity of a pathology increases the tails of the HRV PDFs prompting various proposed methods for obtaining a tempered (reduced) Lévy PDF with finite second moments as a second-order FOC model (FOCM-2), which includes a modest formal change in the FOFPE. The procedure is mathematically elegant but requires the fluctuating signal from the autonomous nervous system (ANS) to be self-regulating and was therefore discarded. Instead, we thought it more reasonable to hypothesize the existence of a cardiovascular control mechanism for a third-order FOC model (FOCM-3), which in Section 5 produces a slightly different FOFPE whose steady-state solution is a tempered Lévy PDF characterizing HRV time series in healthy as well as in certain diseased individuals. In Section 6, we draw some generic conclusions regarding the generality of the obtained results using the new method. Section 7 contains supplemental material clarifying fundamental theory, data processing, etc. Section 7.1 consists of a brief review of the basis for crucial event (CE) time series (CETS) and the reasons why such physiologic time series contain information about being multifractal and why that is important. Section 7.2 contains a micro-review for solving a simple fractal-order kinetic equation (FOKE). Section 7.3 is a review of what goes into the diffusion entropy data processing technique, the DEA. Section 7.4 contains a discussion of the quality of a Gaussian fit to HRV data.

2. Empirical Lévy Statistics for HRV: FOCM-0

The HRV PSD is determined by the correlations captured in the spectral analysis, characterizing temporal ordering of beat-to-beat differences, $I\left(n\right)=B(n+1)-B\left(n\right)$, where $B\left(n\right)$ denotes the length of the time interval for the beat number n. Peng et al. [34] graphed the histogram for the differences in the beat-to-beat intervals $I\left(n\right)$ and found that this is a stationary Lévy process, such as [see Section 7.2] shown in Figure 1. It was determined by them that the $B\left(n\right)$ are anti-correlated for typical healthy individuals using a random walk model. Correspondingly, this correlation was determined to vanish in patients with heart disease, but its vanishing was determined to not influence the PDF of the heartbeat increments. Both healthy and diseased individuals were described by Lévy-stable PDFs without being able to statistically distinguish between them, thereby demonstrating that the distribution of inter-beat fluctuations $\left\{I\right(n\left)\right\}$ does not uniquely characterize healthy HRV time series since it is their time ordering that has accurately captured the fractal nature of the beat-to-beat interval time series.

The time series analysis of HRV performed by Peng et al. [34] indicated that the statistics of the fluctuations in the HRV are Lévy stable, which was obtained using a random walk model. This observation refuted the widely held belief that cardiac activity observed in healthy individuals is normal in a statistical sense. They demonstrated that rather than following a bell-shaped curve of the Gaussian type, the statistics of heartbeat intervals are characterized by heavy tails (see Figure 1), signifying that deviations from the average are far more likely than expected from a Gaussian process, as visually highlighted in Figure 1b. This is the zeroth order FOC model (FOCM-0) since it implies the need for Lévy statistics, thereby entailing the FOC, but that calculus was not used in its empirical determination by Peng et al. and was, in fact, first derived as a solution to an FOC differential equation by West and Seshadri [39] only a decade earlier.

However, such functional dependence shares an undesirable property, which is also observed in the HRV PSD, it being a diverging second moment. This divergence of the HRV fluctuations implies that the heartbeats of all healthy individuals would be asymptotically unstable, a conclusion that is medically unacceptable for a healthy HRV time series. A number of suggestions have been made to overcome this unrealistic divergence of the HRV, which we subsequently touch on. Herein, we adopt some ideas from control theory adapted to the fractional calculus as a means of resolving the objections raised with regard to these early suggestions.

3. Theoretical Lévy Statistics for HRV: FOCM-1

The empirical Lévy PDF was at the center of the analysis of an HRV time series conducted by Peng et al. [34], which suggested to us [40] that the FOC might provide a strategy for going beyond the random walk interpretation of the HRV time series proposed by these groundbreaking investigators. Consequently, West and Trualska [40] introduced a Langevin equation to replace the random walk argument of Peng et al. to determine the fractal form of the diffusion equation used to find the Lévy statistics:

$${\displaystyle \frac{dI\left(t\right)}{dt}}=\xi \left(t\right).$$

(3)

Here, $\xi \left(t\right)$ is the non-Gaussian random force produced by the autonomous nervous system (ANS), which manifests the conflict between the time series produced by the sympathetic system (increases HR) and that produced by the parasympathetic system (decreases HR), which was empirically determined to be an infinitely divisible IPL distribution using characteristic functions.

The probability that the dynamic variable $I\left(t\right)$ lies in the phase space interval ($\mathcal{I},\mathcal{I}+d\mathcal{I}$) at time t is $P(\mathcal{I},t)d\mathcal{I}$, where the PDF is defined as a solution to a fractal-order Fokker–Planck equation (FOFPE) [41]:

$${\displaystyle \frac{\partial P\left(\mathcal{I},t\right)}{\partial t}}=K_{\beta }{\partial }_{\left|\mathcal{I}\right|}^{\beta }\left[P\left(\mathcal{I},t\right)\right],$$

(4)

where we have introduced the symmetric Reisz–Feller fractal-order derivative (FOD) ${\partial }_{\left|\mathcal{I}\right|}^{\beta }\left[·\right],$ whose Fourier transform is $-{\left|\kappa \right|}^{\beta }$ with $\kappa$ the Fourier variable conjugates to the phase space variable $\mathcal{I}$, and $K_{\beta }$ is a ‘diffusion’ constant. Equation (4) is a fractional diffusion equation for the evolution for the PDF in the heartbeat intervals, whose solution is the Lévy distribution given by the integral [41]:

$$P(\mathcal{I},t,\beta )=\underset{0}{\stackrel{\infty }{\int }}{\displaystyle \frac{d\kappa }{\pi }}cos\left[\kappa \mathcal{I}\right]exp\left[-t{K}_{\beta }{\kappa }^{\beta }\right].$$

(5)

The Lévy PDF given by Equation (5) agrees with the histograms of the HRV data depicted in Figure 1 for a fixed time, with the Lévy index $\beta =\gamma$ $=1.7$ and a positive definite diffusion coefficient $K_{\beta }>0$. Note that, in each case, the value for the numerically generated Lévy random force is obtained for a healthy individual, which is provided by the ANS drivers.

The FOFPE is a model FOC equation of evolution assumed to govern the behavior of the PDF for the dynamics of healthy heart rate fluctuations. However its value for the present study lies in its forming the first link in the chain of logic connecting this unsatisfactory model with one, which allows for the inclusion of the cardiac control mechanisms, which satisfy all our present-day modeling concerns.

FOCM-0 and FOCM-1 Are Both Unacceptable

Peng et al. [34] find that the statistics of healthy and diseased (dilated cardiomyopathy) conditions are the same, that being Lévy stable with index $\gamma =1.7$; however, the spectra for the two cases are quite different. The PSD $S_p\left(f\right)$ is given by the square of the Fourier transform of $I\left(n\right)$ yielding $S_p\left(f\right)\propto 1/f^{\beta }$, where $\beta =2\delta -1$ and the mean-square level of the interbeat fluctuations increases as $n^{2\delta }$. Here, again, $\delta =0.5$ corresponds to Brownian motion, so that $\beta =0$ indicates the absence of correlations in the time series $I\left(n\right)$ (“white noise”). See the table and the discussion of the relations among the scaling indices in Section 7.1.

They [34] observed that for a diseased dataset, the PSD IPL index $\beta$ is approximately zero in the low-frequency regime confirming that the $I\left(n\right)$ are not correlated over long times. On the other hand, they also observed that for the healthy dataset, the PSD IPL index $\beta$ is approximately equal to one, indicating a long-time correlation in the interbeat interval differences. The anti-correlated properties of $I\left(n\right)$ are consistent with a nonlinear feedback system that “kicks” the heart rate away from fluctuation extremes. This tendency operates on a wide range of time scales, not on a beat-to-beat basis. The conclusion is that the different scaling patterns must be a consequence of the ordering of the differences, rather than their statistics, i.e., in the correlations produced by the underlying dynamics.

The first step in the logic chain has, therefore, been the introduction of the IPL fluctuations in a Langevin equation giving rise to the FOFPE Equation (4) and the Lévy PDF solution given by Equation (5). However, this model of the HRV statistics is also unacceptable due to the divergence of the second moment in the Lévy model of any phenomena including an HRV time series. However, it does capture what was known at that time, which leads to our formulation of the next link in the chain of FOC modeling logic, which is how to temper the large-scale fluctuations and thereby retain stability of the HRV time series by introducing a reasonable mechanism to reduce the tails of the Lévy PDF while simultaneously maintaining the proper scaling in the central region of the empirical PDF.

4. Tempered Fluctuations: FOCM-2

Figure 2 depicts a number of representative exemplar datasets of HRV time series for a healthy control group (top row) and two subsequent rows with pathological HRV time series [see Section 7.4 for details]. The value of the parameter ${\lambda }_{25}$ in the upper right-hand corner in each of the final panels of the three rows provides a measure of the deviation of the time series of interest from the best Gaussian fit to the depicted data [42]. Scanning down the figure, it is clear that the larger the value of ${\lambda }_{25}$, the less the tails of the fitted Lévy PDFs conform to those required for healthy function. Notice that the red curve in the bottom row flares out to $\pm 10$ standard deviations indicating that this person’s heart spends a large fraction of its time in the off or relaxed position.

We hypothesize a physiologic feedback mechanism that exponentially suppresses very large fluctuations at both positive and negative extremes. The feedback mechanism in FOCM-2 is equivalent to the so-called tempering mechanism previously used by [22] to explain the three PDFs characterizing a cohort group of 670 post-AMI (acute myocardium infarction) patients using 24 h Holter monitor data to construct heartbeat interval data from the HRV time series depicted in the figure from [22] (see next section). It has been argued that a modified FOFPE yields the tempered Lévy PDF and can be obtained by the direct inclusion of a dissipation rate in the Reisz–Feller fractal-order derivative term in Equation (4) to obtain FOCM-2:

$${\displaystyle \frac{\partial P\left(\mathcal{I},t\right)}{\partial t}}=K_{\beta }{({\partial }_{\left|\mathcal{I}\right|}+\lambda )}^{\beta }\left[P\left(\mathcal{I},t\right)\right].$$

(6)

The exact solution to this equation is given by an exponentially truncated Lévy PDF; see Section 7.6.1 of [18] for details:

$$P\left(\mathcal{I},t;\beta ,\lambda \right)=e^{-\lambda \left|\mathcal{I}\right|}{L}_{\beta }\left(\mathcal{I},K_{\beta }t\right),$$

(7)

where $L_{\beta }\left(\mathcal{I},K_{\beta }t\right)$ is the Lévy PDF given by Equation (5). The completely scaled form of the PDF is determined in terms of the probability for the scaled variable ${\mathcal{I}}_s=\mathcal{I}{\left(K_{\beta }t\right)}^{-1/\beta }$:

$$P\left(\mathcal{I},t;\beta ,\lambda \right)d\mathcal{I}=e^{-{\lambda }_s\left|{\mathcal{I}}_s\right|}{L}_{\beta }\left({\mathcal{I}}_s\right)d{\mathcal{I}}_s,$$

(8)

where we have used the infinite divisibility of the PDF to obtain the exact form of the solution to Equation (6) given by Equation (8).

Why Is FOCM-2 Unacceptable?

Figure 3 depicts the empirical result that those surviving and those that have died in the cohort group of 670 patients due to non-cardiac causes have essentially the same distribution due to the values of the hypothesized cardiac control mechanisms, which are significantly larger than those who suffer cardiac deaths. Said simply, those who suffer cardiac deaths do so because the control mechanism that suppresses the larger fluctuations has failed. This elegant mathematical result has a fatal drawback which is, of course, that the physiological source mechanism for the dissipation parameter ${\lambda }_{25}$ in Equation (6) resides in the ANS model, if it existed (it does not), and not in the model’s sinoatrial pacemaker cell. FOCM-2 would be the preferred model if the fluctuations and dissipation were generated by the same process, but that is not the physiological case presented here.

We need a control mechanism that “kicks” the heart rate fluctuations away from its extremes as so colorfully explained by Peng et al. [34], whereas FOCM-2 suppresses the extreme fluctuations before they have a chance to disrupt the HRV time series, which is not the situation observed by Peng et. al., and for that reason, the model must be rejected in favor of an explicit nonlinear control process that can explain the results recorded in both Figure 2 and Figure 3 equally well.

5. HRV Control Hypothesis: FOCM-3

As the third and last step in our logic chain for the HRV example, we hypothesize the existence of cardiovascular control mechanisms, which leads to a tempered Lévy PDF, characterizing CETS in healthy individuals using an FOC model. This control mechanism suppresses the larger fluctuations of those in the Lévy PDF tail, which persists in various cardiac disorders. The pathophysiology of the HRV PDF being Lévy stable, therefore, results from this physiological control process, suppressing the extreme fluctuations in the HRV, and differs from the exponential tempering of the Lévy PDF found in the literature, e.g., see [24,40,43]. Here, we again model the HRV time series to be determined by a nonlinear Langevin equation, where the nonlinear deterministic force corresponds to the control process acting to suppress the fluctuations exceeding a given physiological level beyond which the ON becomes unstable. The corresponding Fokker–Planck equation takes on a fractional form, whose steady-state solution is a tempered Lévy distribution. West and Turalska [40] captured the dynamics of cardiovascular-like systems in the framework of an FOC and their results are adapted to the needs of the present FOC model.

The HRV time series is a consequence of numerous control mechanisms that shape the cardiac activity, and the most fundamental component of the HRV is formed by the firing of the pacemaker cells in the heart’s sinoatrial node [44]. This periodic action is continuously driven by the two competing branches of the ANS with sympathetic stimulation increasing the firing rate and parasympathetic action decreasing that rate [34,36,45], as stated earlier. Now we address the significance of the heart’s primary pacemaker, the sinoatrial node, in controlling the balance in the competition between these two branches of the ANS. The importance of the response of the sinoatrial node to autonomic signaling of heart rate and HRV has been emphasized by [40,46].

IPL spectra have been observed in a number of dynamical systems with chaotic solutions; see, for example, Reichl [47]. Goldberger and West [36] suggested that the observed PSD of the HRV may be a consequence of such dynamics and that one may interpret the modifications in the IPL PSD as being indicative of pathology and, therefore, may be of diagnostic and prognostic value. Loss of complexity in the HRV has already been described in numerous settings, including multiple sclerosis [48], diabetes mellitus with neuropathy [47], fetal distress [49], bed-rest deconditioning [50], aging [51] and in certain patients at risk for sudden cardiac death [52,53]. Presumably, the more severe pathologies will be associated with the greatest loss of spectral power, which have been referred to as losses of spectral reserve [36,54].

West and Turalska [40] point out that there is no a priori reason to restrict the control analysis of HRV dynamics to linear negative feedback. In physical systems, the fluctuations and linear dissipation have a common origin and consequently are tied together by Einstein’s fluctuation–dissipation theorem [55], which implies an underlying linear response mechanism. In the life sciences, however, the significant random force and dominant determinisic force that appear together in the Langevin equation have different sources, and the deterministic feedback must be adaptive to regulate the random force while simultaneously responding to it. They [40] assumed in the case of the HRV time series that the deterministic force $F\left(I\right)$ is generated by the collective behavior of the sinoatrial pacemaker cells, and the sum of the competing contributions from the sympathetic and parasympathtic components of the ANS is $\xi \left(t\right)$ are written as a fractal-order nonlinear Langevin equation (FONLE):

$${\displaystyle \frac{dI\left(t\right)}{dt}}=F\left(I\right)+\xi \left(t\right).$$

(9)

This argument was used to determine the FOFPE, which when the deterministic force is taken to be linear, the solution to the fractal partial differential equation was determined by West and Seshadri [39] to yield a Lévy PDF solution, as discussed with regard to the FOCM-1.

A generic version of Equation (9) was analyzed by Chechkin et al. [56] wherein the deterministic term was given by the nonlinear force written as $F\left(I\right)=$$-\lambda \left(I\right)I:$

$${\displaystyle \frac{dI\left(t\right)}{dt}}=-\lambda \left(I\right)I+\xi \left(t\right),$$

(10)

Equation (10) reduces to the linear Langevin equation of the last section with ${\lambda }_0\equiv \lambda \left(0\right)$, and $\xi \left(t\right)$ is again assumed to be a renewal process with index $\beta$ [14]. In general, the feedback term is a polynomial, but [40] restricted the analysis to the highest-order non-vanishing term, since this term asymptotically dominates the control process and wrote $\lambda \left(I\right)\approx {\lambda }_0+{\lambda }_{2n}{I}^{2n}.$ But here, this replacement is considered to be unsatisfactory because the nonlinear form of the dissipation coefficient $\lambda \left(I\right)$ was selected solely to overcome the divergence of the variance $〈I^2〉.$

Herein, we acknowledge the weakness of the model put forward in [40] and replace the deterministic force used therein, in keeping with the overall strategy of learning from previous defects in modeling the FONLE. The error in the earlier work is corrected for the present model by replacing the deterministic force in the FONLE by:

$$F\left(I\right)\approx -\left({\lambda }_0 +\ge {}_{\nu }\left|I\right|^{2\nu }\right)I$$

(11)

where the power-law index $\nu$ herein can be non-integer, and the absolute value of $I\left(t\right)$ ensures that the deterministic force is symmetric. Note that Equation (11) is equivalent to our replacement used in [40] with ${\left|I\right|}^{2\nu }=I^{2\nu }$, where $\nu \ge 0.$ We also mention the recent work of Mackay [57], who analyzed non-equilbrium noise in a linear V-shaped potential whose solutions overlapped those studied by [40,56].

The technique presented in Section 7.2 enables the replacement of the FONLE containing the lowest-order nonlinear term in the sinoatrial node dynamics given by Equation (11) with a corresponding FOFPE that is expressed in the present situation as:

$${\displaystyle \frac{\partial P\left(\mathcal{I},t\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial \mathcal{I}}}\left[({\lambda }_0 +\ge {}_{\nu }\left|\mathcal{I}\right|^{2\nu })\mathcal{I}P(\mathcal{I},t)\right]+K_{\beta }{\partial }_{\left|\mathcal{I}\right|}^{\beta }\left[P\left(\mathcal{I},t\right)\right].$$

(12)

Note that this equation is the same ‘form’ of FOFPE studied by Chechkin et al. [56] where the ‘drift’ term ${\lambda }_{0}I$ is analogous to the harmonic Ornstein–Uhlenbeck potential. The next integer-order contribution in their analysis was ${\lambda }_1{I}^3$ and would be analogous to a quartic potential. It is clear that the quartic potential suppresses the divergent second moment, but the fourth moment still diverges and requires the next-order dissipation to quell that divergence.

Of course, we cannot solve the full time-dependent FOFPE given by Equation (12) in general, but it is not necessary to know the exact solution to obtain the information we need to test our present control hypothesis. The asymptotic behavior of the steady-state solution is sufficient for that purpose. The asymptotic steady-state solution to the FOFPE is given by setting Equation (12) to zero and retaining only the lowest-order nonlinear term to obtain:

$${\lambda }_{\nu }{\mathcal{I}}^{2\nu +1}{P}_{ss}\left(\mathcal{I}\right)=K_{\beta }{\displaystyle \frac{d}{d\mathcal{I}}}\underset{-\infty }{\stackrel{\mathcal{I}}{\int }}{\displaystyle \frac{P_{ss}\left({\mathcal{I}}^{\prime }\right)d{\mathcal{I}}^{\prime }}{{\left|\mathcal{I}-{\mathcal{I}}^{\prime }\right|}^{\beta -1}}}.$$

(13)

A formal argument using the ansatz $P_{ss}\left(\mathcal{I}\right)\backsim C/{\left|\mathcal{I}\right|}^{\gamma },\gamma >0$ yields the normalized solution for the approximated steady state to be the IPL PDF (see Section 7.2 for details) [40,56]:

$$P_{ss}\left(\mathcal{I}\right)={\displaystyle \frac{N_{\gamma }}{{\left|\mathcal{I}\right|}^{\gamma }}};\gamma =\beta +2\nu +1,$$

(14)

where $N_{\gamma }$ is the normalization constant, and the equation is valid for $\mathcal{I}\to \pm \infty$ due to symmetry. Here, we paraphrase the conclusion reached by Chechkin et al. [56] that for all the values of the IPL index greater than the critical value $\nu >{\nu }_{cr}=1$, the variance $〈{\mathcal{I}}^2〉$ is finite, and consequently, a nonlinearity whose highest power $\nu$ exceeds or is equal to the crucial value ${\nu }_{cr}=1$ counterbalances the information overload injected by the complexity of the driving ANS fluctuations to the information exchange with the sinoatrial node.

In Figure 4, three steady-state IPL PDFs are calculated using the asymptotic solution to the FOFPE at the two extremes of the interval $2.7\le \gamma \le 4.7$, with two of them being degenerate when $\nu =0$. The steady-state scaling index $\gamma$ of the PDF is fit to the numerical solution of the steady-state FOFPE (see Section 7.2). In the numerical simulations, $\beta =1.7$ is used to drive the dynamics of the sinoatrial node modeled as a dynamic network near a tipping point. It is clear from Figure 4 that the deterministic force is selected to have the formal steady-state solution that captures the three IPL forms predicted by the steady-state solution to the FOFPE with either $\gamma$ $=$ with $\nu =0$ (pathological) or $\gamma$ $=4.7$ with $\nu =1$ (healthy). It is clear that the extreme values of the control index $\nu$ act like a switch for turning the physiologic control of HRV fluctuations on or off indicating that there is a significant potential range of $\nu$ over which clinicians may be able to intervene and restore control of the HRV to an appropriate healthy level.

Why Is FOCM-3 Acceptable?

The model is acceptable because it provides a reasonable separation into healthy and diseased individuals, as shown in Figure 4, as well as further discrimination between those who have died due to heart disease from those who have passed due to other causes, yielding the interesting partitoning into the three groupings of patients recorded in Figure 3.

It is well known that for an IPL index $\gamma \ge 3$, the second moment $〈{\mathcal{I}}^2〉$ is finite as a consequence of the central limit theorem holding asymptotically. Consequently, for $\nu =1$, the statistics of $\mathcal{I}$ have a finite second moment asymptotically, and the steady-state PDF asymptotically transitions from Lévy stable to Gaussian. This transition explains the result obtained by Kiyono et al. [43] without the need to exponentially temper the random force in the FOLLE, relying solely on the sinoatrial node dynamics to temper the extreme fluctuations resulting from the competition between the two branches of the ANS.

Figure 4. The numerical integration of Equation (12) is used to determine the steady-state PDF with the empirical Lévy index = 1.7. The computational results are fit extremely well by the analytic steady-state IPL PDF given by Equation (14) with IPL index $\gamma =\beta +2\nu +1$ for $v=0$ : ${\lambda }_0={\lambda }_{\nu =0}=0$ (black squares); ${\lambda }_0=1, {\lambda }_{\nu =0}=0$ (blue circles): and for $v=0$ : ${\lambda }_0=1.0, {\lambda }_1=0.1$ (green diamonds). The dashed curves denote the analytic results where the healthy individual $v=1$ has a control index of $\gamma =4.7$ (green dashed line), and a person with a totally inoperable control index $v=0$ is $\gamma$ = 2.7 (black dashed line). Note the modest deviation in the numerical calculation away from the IPL due to the linear dissipation. Adapted from [40] with permission.

Figure 4. The numerical integration of Equation (12) is used to determine the steady-state PDF with the empirical Lévy index = 1.7. The computational results are fit extremely well by the analytic steady-state IPL PDF given by Equation (14) with IPL index $\gamma =\beta +2\nu +1$ for $v=0$ : ${\lambda }_0={\lambda }_{\nu =0}=0$ (black squares); ${\lambda }_0=1, {\lambda }_{\nu =0}=0$ (blue circles): and for $v=0$ : ${\lambda }_0=1.0, {\lambda }_1=0.1$ (green diamonds). The dashed curves denote the analytic results where the healthy individual $v=1$ has a control index of $\gamma =4.7$ (green dashed line), and a person with a totally inoperable control index $v=0$ is $\gamma$ = 2.7 (black dashed line). Note the modest deviation in the numerical calculation away from the IPL due to the linear dissipation. Adapted from [40] with permission.

6. Discussion and Conclusions

There are a number of take-home lessons to be derived from the presentation, starting with never being satisfied with the early FOC models developed to explain a complex process, which is not to say the argument or model should not be published. However, no matter how elegant the mathematics, it is important to repeatedly revisit until you have answered all your questions and those of your colleagues. Most importantly, make notes of the changes in your thinking as the levels of FOC modeling unfold. Therefore, let us quickly review what was learned at each step of the modeling process presented herein.

The analysis of HRV data by Peng et al. [34] indicated that the statistics of the changes in heartbeat intervals are Lévy stable. The Lévy statistics have the property that the second moment of a time series with such statistics diverges, consequently implying that the heartbeat of every healthy person would be asymptotically unstable. A number of investigators have argued that this instability is readily quenched by introducing a mechanism to dampen the largest fluctuations produced by the Lévy driver, thereby producing a tempered Lévy PDF; see, e.g., [22,58].

FOCM-0 verified the processing of the empirical HRV time series for the healthy state of the variability of the heartbeat intervals as a Lévy PDF, but ultimately, it is unsatisfactory because the Lévy PDF has an unphysiologic divergence of the second moment. The FOCM-0 is too simple.

The scaling of the statistics for the central or Lévy part of the tempered PDF is the same for both the healthy and CHF patients in Figure 2. Consequently, those who suffer from a serious cardiological condition have greater heart rate variability, but their HRV is certainly less complex than those whose HRV variability, which is restricted by a nonlinear control process. Thus, interpreting complexity to be proportional to the second moment (variability) in the HRV context is too simplistic.

FOCM-1 was designed to provide a formal mathematical description of the Lévy PDF empirically observed by [34] in FOCM-0. The natural calculus to obtain this PDF is the FOC but the FOFPE given in Equation (4) requires additional modifications to incorporate the tempering mechanism in the dynamics and thereby avoid the second-moment divergence.

The data shown in Figure 2 suggest the use of exponentially tempered Lévy statistics to explain the difference between healthy individuals and CHF patients. The tempering control mechanism would suppress the largest excursions in time, which persist in the severe CHF groups. The pathophysiology of the HRV PDFs being Lévy stable would then be the result of the suppression of a physiological control mechanism, as suggested in [40]. However, the form of this earlier hypothesized control mechanism would need to operate at the level of the involuntary nervous system and not be a separate and independent balancing mechanism, such as the one proposed herein and supported by the empirical evidence.

FOMC-2 uses the infinite divisibility of the PDF to obtain Equation (6), which has an elegant tempered Lévy PDF solution. However, it does not realize the proper control mechanism to suppress the heart rate at the extremes and instead suppresses the extreme fluctuations themselves before they can disrupt the HRV time series.

However, there is the additional scaling of the HRV in the tempered Lévy PDF produced by the criticality of the healthy sinoatrial node dynamics. Criticality, as manifest in the nonlinear term in the new FONLE given by Equation (11), reduces the span of the fluctuations from those that persist in severe CHF conditions to those present in more manageable CHF cases. This second scaling, which is produced by the proposed control mechanism, certainly adds to the overall complexity of the process, by introducing a second scaling mechanism. It is the weakening and ultimate loss of this control modeled by the control parameter $\nu \to 0$ that lifts the restriction on the span of the fluctuations in HRV time series, resulting in certain cardiac pathologies. The source of the physiological control mechanism in the FONLE given by Equation (11) is based on the sinoatrial node dynamics being poised at criticality [40] and is inserted into Equation (9) to create the full nonlinear control mechanism in the FOFPE Equation (12).

They [40] noted further that the sinoatrial node is comprised of a network of interacting pacemaker cells, the collective behavior of which provides insight into critical phase transitions in non-physical networks. The model they chose to mimic the complexity of the sinoatrial node was the Decision Making Model [59], which, near criticality, reduces Equation (10) to a cubic Langevin equation, which although qualitatively correct, like the tempered Lévy PDF discussed earlier, errs in quantitative detail for this particular application.

FOCM-3 solved this remaining problem by introducing a nonlinear negative feedback mechanism that is self-regulating and is asymptotically consistent with the hypothesis that “disease is the loss of complexity”, made by [36].

Let us use our new-found insight to explain the three curves in Figure 3. Here, we separate the cohort groups into the living and the dead. When the hypothesized control network is on and working properly, i.e., when $\nu =1,$ the largest fluctuations in the HRV time series are suppressed, thereby producing a tempered Lévy PDF. This would account for the survivor curve consisting, as it does of only the living. This first cohort group consists of those alive all having $\nu =1$, and the second cohort group consists of those who are all dead. The latter group can be usefully further subdivided into those who died by extreme cardiovascular malfunction and those who are equally unlucky but whose hearts are fine up to the point where they no longer beat, but the cause of death must be sought outside the cardiovascular network. It is this latter subdivsion that accounts for the cardic ($\nu =0$) and non-cardiac ($\nu =1$) deaths.

Thus, we may tentatively conclude that the FOC provides a sufficiently general infrastructure to build successively more complex models of real physiological networks, as unseen tasks being performed in the background of healthy networks are revealed when they fail to properly respond to the signals from other networks due to illness or injury ($IoI$). This opens the door to a general theory of rehabilitation using crucial events [60]. A correct model of the diagnostic for $IoI$ as in the case of the weakening of the control mechanism (${\nu }_C$ $<1$) often suggests how the $IoI$ can be healed through the introduction of an aritificial control mechanism (${\nu }_A$ $<1-{\nu }_C)$ to temper the extreme fluctuations such that the transform of a dynamic variable $\nu ={\nu }_C+{\nu }_A\approx 1$, thereby restoring the heart’s rehabilitated HRV to a healthy scaling value.

7. Methods

7.1. On Crucial Events (CEs)

We have hypothesized [12] that the ONs making up a NoON such, as the human body, each generates a time series $X\left(t\right)$ consisting of a sequence of CEs modulated by a multifractal dimension (MFD) index $\delta \left(t\right)$, such that by scaling the time, we obtain the scaling relation $X\left(\lambda t\right)={\lambda }^{\delta \left(t\right)}X\left(t\right).$ Empirically, the phenomenon of a sequence of crucial events (CEs) is based on time series consisting of random discrete events with renewal statistics, which enable two interacting organ-networks (ONs) to detect and quantify the complexity of the other ON and to minimize the difference in their relative complexity so as to become synchronized. The quantification of complexity is given by the MFD, which becomes equal in each ON time series (ONTS) even though they may be operating on largely different time scales, such as the brain, heart and lungs; see [61,62,63] for a discussion of the synchronization of the time-depenent scaling index $\delta \left(t\right)$. The discrete events in such time series have been named crucial because they emerge from the nonlinear complex sub-dynamics to determine the efficiency of information exchange among the members of the NoONs.

A crucial event time series (CETS) is generated by a renewal event time series (RETS) with an IPL waiting-time PDF $\psi \left(\tau \right)\propto {\tau }^{-\mu }$, with an IPL index in the domain 1 < $\mu$ < 3. Asymptotically, the generated CETS describes an ergodic process for the IPL index $\mu$ in the range 2 < $\mu$ < 3 with a finite average waiting time. Ergodic is the technical name for statistical processes for which averages taken over long time series are equal to those taken over PDFs. The understanding of complexity in many-body physics is formally understood by assuming the ‘ergodic hypothesis’, dating back over a century to the time of Boltzmann.

The function $F\left(y\right)$ of the scaled variable $y=x/t^{\delta }$ in a scaled PDF in Equation (21) is unknown in general but is well known for certain values of the scaling parameter. This generic form of the solution is obtained from an equation of evolution of the PDF that has an $\alpha —$order derivative in time and a $\beta —$order derivative in space (see Equation (1)) such that the solution to the fractal-order diffusion equation (FODE) can be solved using renormalization group theory, as shown in FOCM-2 to obtain the scaling PDF with the scaling index $\delta =\alpha /\beta$ [12]. The ordinary diffusion equation has the integer-order indices $\alpha =1$ and $\beta =2$, which yields the well-known scaling index $\delta =1/2$, in which case the unknown function $F\left(y\right)$ becomes Gaussian with the scaled variable $y=x/\sqrt{t}$. Note that the fractal dimension $\mathcal{D}$ for this process is obtained from $\mathcal{D}=2-\delta =1.5$, which is, therefore, a monofractal process that is completely random, which is to say, it is simple diffusion with no memory. Consequently, scaling PDFs with $\delta \ne 1/2$ are anomalous in that their fractal dimension and $D\ne 1.5$, and they contain various forms of memory.

If the anomaly has a completely temporal ($\alpha$ $<1)$ origin we have for the spatial derivative index $\beta =2$ for a spatially homogeneous process and $\delta <1/2$ thereby producing the fractal dimension $\mathcal{D}>1.5$ indicating the existence of a memory-like process. The more the scaling index deviates from 1/2, the more the fractal dimension deviates from 1.5, and the stronger the anomaly. If the time deviation is produced by a fractal time series with an IPL waiting-time PDF, when these time intervals are statistically independent of one another, it constitutes a renewal time series. The renewal time series is ergodic for $2<\mu <3$, with a finite average waiting time, and is non-ergodic for $1<\mu <2$, with a diverging average waiting time.

Suppose the trajectory $X\left(t\right)$ crosses a known level at a specific time and we want to know how long we must wait to recross that same level. Given that the waiting-time PDF has the IPL form, we denote the generic IPL index for the waiting-time PDF $\mu$ by the symbol ${\mu }_{\mathcal{D}}.$ From the parameter relations in Figure 5, it is clear that it is possible to prove that the IPL index is equal to the fractal dimension, so we obtain [63,64]:

$${\mu }_{\mathcal{D}}=\mathcal{D}=2-\delta ,$$

(15)

which is possible to establish using the probability of crossing and recrossing any fixed value of the diffusion trajectory, which has been shown to be a renewal process and is consequently a CETS.

The PSD $S_p\left(f\right)$ for a CETS is also IPL in terms of the frequency f:

$$S_p\left(f\right)\propto f^{-\beta },$$

(16)

indicating $1/f—$variability. This IPL PSD is the Fourier transform of a slowly decaying auto-correlation function indicative of a fractal time series having global or local self-similarity. The IPL indices $\mu$ and $\beta$ are interrelated for the class of fractal time series of interest in medicine such that the fractal dimension given by Equation (15) is related to the IPL PDF index defined by Equation (16), and the IPL index for the PSD is [12]:

$$\beta =3-\mu .$$

(17)

Thus, true $1/f—$noise at $\beta =1$ arises only at the boundary $\mu =2$ between the ergodic and non-ergodic variability domains.

Figure 5 is copied from West et al. [60] and provides easy reference to all the derived relations among the IPL scaling indices for the PSD index $\beta$, the IPL waiting-time PDF index $\mu$ and the scaling index $\delta$ for the scaled variable $X\left(t\right)$ as well as the scaled PDF.

The work in [65,66] has demonstrated reduced multifractality associated with impaired parasympathetic control and the pathological condition of CHF, respectively. This corresponds to a reduction in long-range correlations present in the HRV of healthy individuals, an outcome related to a loss of dynamical properties of the cardiovascular system. The control mechanism introduced herein is presented in terms of a FOFPE [16] given by Equation (12), which leads to investigations of HRV through its PDF distribution. The FOFPE is a theoretical model that describes the evolution of the PDF of the physiologic process, and its solution (when it can be obtained) captures the full dynamics of an ensemble of time series [67].

The adaptive feedback control process produced by the collective behavior of the pacemaker cells, postulated herein, is physiologically simpler than the exponential tempering suggested earlier and, applying Occam’s razor, is the more likely explanation. The nonlinear term is a collective effect within the sinoatrial node dynamics and produces a tempered response to the stochastic excitations. This tempering of the extreme HRV excursions is pathologically suppressed within the severe CHF group. The pathophysiology of the HRV PDF being Lévy stable is then the result of a suppressed control process, which can no longer determine the maximum size of the interbeat intervals, resulting in losing control of the disease [68].

7.2. Scaling Solution to Fractal-Order Kinetic Equation

Zaslavsky et al. [15,17] applied the renormalization group transformation R to the network dynamics such that the scaling properties of the incremental changes are: $R:\delta x={\lambda }_X\delta x;R:\delta t={\lambda }_T\delta t,$ which apply, after some averaging, to a restricted space-time domain, and the scaling parameters are (${\lambda }_X;{\lambda }_T$). They continue with the observation that the fractal-order kinetic equation (FOKE) given by Equation (1) is invariant under the renormalization group transformation:

$$R:D_t^{\alpha }\left[P\left(x,t\right)\right]=K_{\beta }R:D_{\left|x\right|}^{\beta }\left[P(x,t)\right]$$

which implies that the transformed fractal equation satisfies the scaling behavior:

$${\lambda }_T^{\alpha }{D}_t^{\alpha }\left[P\left(x,t\right)\right]=K_{\beta }{\lambda }_X^{\beta }{D}_{\left|x\right|}^{\beta }\left[P(x,t)\right].$$

(18)

The lowest-order renormalization group solution is given by equating the renormalization parameters raised to their respective powers: ${\lambda }_T^{\alpha }={\lambda }_X^{\alpha }$. The solution to the renormalized FOKE is given in terms of the Fourier transform of the PDF, which is the characteristic function $\varphi (k,t)$, expressed in terms of the Mittag–Leffler function $E_{\beta }(.)$ as:

$$\varphi \left(k,t\right)=E_{\beta }\left(-K_{\beta }{\left|k\right|}^{\beta }{t}^{\alpha }\right).$$

(19)

Consequently, taking the inverse Fourier transform of this characteristic function and expressing the Mittag–Leffler function as an infinite series, after some algebra [17,18], results in the scaling solution for the PDF:

$$P\left(x{\lambda }_T^{\alpha /\beta },{\lambda }_{T}t\right)=P\left(x,t\right)/{\lambda }_T^{\alpha /\beta }.$$

(20)

Selecting the ratio of the scaling parameters to be the exponent of the remaining scaling parameter: $\delta =\alpha /\beta ,$ and choosing the time scale of interest to reduce the dependence of the PDF on the left side of the equality in Equation (20) to be dependent on the single scaled variable $x/t^{\delta }$ using ${\lambda }_T=1/t,$ enables us to rewrite Equation (20) in the form of the scaling PDF given by

$$P(x,t)={\displaystyle \frac{1}{t^{\delta }}}F\left({\displaystyle \frac{x}{t^{\delta }}}\right).$$

(21)

Thus, the solution to the FOKE yields the scaling PDF along with the DEA data processing technique to interpret the datasets from the brain, heart and lungs in the text as CETS with synchronized levels of complexity, as measured by their respective multifractal dimensions.

Steady-State Solution

Chechkin et al. [56] argue that the behavior of the steady-state PDF can be determined from an analysis of the contributions to the integral in the fractal derivative in Equation (12) in the vicinity of the pole:

$${\lambda }_{2\nu }{\mathcal{I}}^{2\nu +1}{P}_{ss}\left(\mathcal{I}\right)-C_{\beta }{\displaystyle \frac{d}{d\mathcal{I}}}\underset{-\infty }{\stackrel{\mathcal{I}}{\int }}{\displaystyle \frac{P_{ss}\left({\mathcal{I}}^{\prime }\right)d{\mathcal{I}}^{\prime }}{{\left|\mathcal{I}-{\mathcal{I}}^{\prime }\right|}^{\beta -1}}}=0.$$

(22)

In particular, the dominant contribution to the integral arises from the region to the left of the pole $\mathcal{I}$ = ${\mathcal{I}}^{\prime }$, and consequently, the operator may be approximated by the expression shown below. The integral term may be expressed exactly as:

$$\underset{-\infty }{\stackrel{\mathcal{I}}{\int }}{\displaystyle \frac{P_{ss}\left({\mathcal{I}}^{\prime }\right)d{\mathcal{I}}^{\prime }}{{\left|\mathcal{I}-{\mathcal{I}}^{\prime }\right|}^{\beta -1}}}={\mathcal{I}}^{1-\beta }\underset{-\infty }{\stackrel{\mathcal{I}}{\int }}{\displaystyle \frac{P_{ss}\left({\mathcal{I}}^{\prime }\right)d{\mathcal{I}}^{\prime }}{{\left|1-\frac{{\mathcal{I}}^{\prime }}{\mathcal{I}}\right|}^{\beta -1}}}$$

so that expanding the numerator and noting:

$$\underset{-\infty }{\stackrel{\mathcal{I}}{\int }}{P}_{ss}\left({\mathcal{I}}^{\prime }\right)d{\mathcal{I}}^{\prime }\approx 1$$

yields the lowest-order approximation to the numerator expansion:

$$\underset{-\infty }{\stackrel{\mathcal{I}}{\int }}{\displaystyle \frac{P_{ss}\left({\mathcal{I}}^{\prime }\right)d{\mathcal{I}}^{\prime }}{{\left|\mathcal{I}-{\mathcal{I}}^{\prime }\right|}^{\beta -1}}}\backsim {\mathcal{I}}^{1-\beta }.$$

Inserting this expression for the integral into Equation (22) and taking the derivative yields, after some algebra, the asymptotic symmetric steady-state PDF given by Equation (14).

7.3. Diffusion Entropy Analysis (DEA)

In this section, we indicate how to determine the consistency of the empirical IPL scaling indices with the theory presented in the previous section and summarized from the appendices in [60,62]. This consistency is accomplished using the DEA to process the recorded time series generated by the members of the triad of ONs consisting of the heart, lungs and brain. The main text in the above two papers contains the results of applying this method to the 66 ONs of the triad. The sequential steps in the DEA technique referenced to Figure 6 are:

(1).

The raw data of each channel are projected onto the interval [0,1] by normalizing each time series by the total time interval of the dataset, thereby enabling the processing of each time series to be directly compared.

(2).

Divide the normalized data profile into parallel stripes of the size of 0.01 (panel a, ECG data).

(3).

Extract events by defining them as unit amplitude pulses if the signal at that time is in a different stripe with respect to its previous value (panel b) and zero if it remains in the same stripe.

(4).

Create a diffusion trajectory (panel c) using the time series of the extracted events from step 3.

(5).

Determine the statistics of a single diffusion trajectory by selecting a window size w and partitioning the diffusion trajectory into many pieces, each starting from an event.

(6).

Initiating all the trajectories from an event enables their being shifted to start from a common origin (panel d).

(7).

Finally, we evaluate the ensemble distribution of histograms at a given time (panel (e)) since the events are statistically independent.

To make the statistics of the single time series diffusion trajectory correspond to what was conducted in the DEA processing of the data, we picked a window $\mathit{w}=\mathit{t}$ starting from a common origin event (panel d) and evaluated the distribution of trajectories at time t (panel e). Denoting the time-dependence PDF for different windows as $\mathit{P}(\mathit{x},\mathit{w})$, we can define the SW-entropy as:

$$S\left(w\right)=-\int dxP(x,w){log}_{2}P(x,w).$$

(23)

Assuming that $\mathit{P}(\mathit{x},\mathit{w})$ is the PDF corresponding to window size $\mathit{w}$, we can define diffusion entropy using the SW-entropy as being the information contained in the time series. Using the scaling PDF, without knowing the $F(.)$ PDF, the deviation of the SW-entropy from its reference state defined by the unknown function is:

$$\Delta S\left(w\right)=S\left(w\right)-S_{ref}=\delta {log}_{2}w.$$

(24)

Consequently, if a graph of the SW-entropy for an empirical process versus the logarithm of the time yields a straight line, we interpret the positive slope to be the scaling index $\delta$, as in the main text.

It is clear that the entropy given by Equation (24) is a theortical quantity being determined by the form of the scaling PDF. On the other hand, it is to be used to determined if an empirical time series can be obtained that has the scaling statistical behavior given by Equation (23) using DEA. The latter is an experimental question that can only be answered by going over multiple values of w to determine which value range yields reasonable values for the scaling index. Of course, it may be true that no values of w give the proper slope that can be interpreted as a scaling index in which case the theory does not apply.

The healthy living network hosts healthy ONs by maintaining equal strength in transferring information among ONs in a NoONs but gives rise to time series whose multifractal fluctuations contain control information that guides both the internal behavior in intra-ON and external information exchange during inter-ONs within the NoONs. In fact, the health of the human body is determined by the multifractal dimensions (MFDs) $D_j\left(t\right)$ or equivalently by the scaling indices ${\delta }_j\left(t\right)$ due to the generic relation between the two quantities: $D_j\left(t\right)=2-{\delta }_j\left(t\right)$. The scaling index has the ideal value of 1 as the condition for the health of the human body. It is a singular condition corresponding to the largest possible scaling and has been used in the analysis of time series as a diagnostic indicator separating patients who have congestive heart failure from those that are healthy [7].

7.4. Non-Gaussian Index

The analysis of non-Gaussianity of HRV starts from integrating and detrending (third-order polynomial for the trend was used herein) R-R interval time series [43,58]. After the detrending procedure, heart rate increments at scale s are defined as ${\Delta }_{s}B\left(t\right)=B^{*}(t+s)-B^{*}\left(t\right)$, where $B^{*}\left(t\right)$ is a deviation of the HRV from the polynomial trend. Then, the PDF of differences ${\Delta }_{s}B\left(t\right)$ normalized by the SD is evaluated. Finally, the non-Gaussianity index ${\lambda }_s$ is estimated as:

$${\lambda }_s=\sqrt{{\displaystyle \frac{2}{q\left(q-2\right)}}ln\left({\displaystyle \frac{\sqrt{\pi 〈{\left|{\Delta }_{s}B\right|}^q〉}}{2^{q/2}}}\right)-ln\Gamma \left({\displaystyle \frac{1+q}{2}}\right),}$$

(25)

where $〈|{\delta }_s{B|}^q〉$ denotes an estimated value of the $q$-th order absolute moment of $|{\Delta }_{s}B|$. Following [22], we calculate ${\lambda }_s$ based on the 1/4-th-order moment ($q=0.25$), and we set the scale s at $25s$. The parameter ${\lambda }_s$ was originally introduced as a parameter of a log-normal cascade model, describing non-Gaussian distributions in the study of intermittency of hydrodynamic turbulence [42]:

$$P_{\lambda ,{\sigma }_0}\left(x\right)=\stackrel{\infty }{\underset{0}{\int }}{\displaystyle \frac{dln\sigma }{\sqrt{2\pi }\lambda }}exp\left(-{\displaystyle \frac{{ln}^2\left(\sigma /{\sigma }_0\right)}{2{\lambda }^2}}\right)\times {\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}exp\left(-{\displaystyle \frac{x^2}{2{\sigma }^2}}\right).$$

(26)

The value of ${\lambda }_s$ close to zero corresponds to a PDF close to a Gaussian distribution, while a larger value of ${\lambda }_s$ means that the observed PDF has fatter tails and a sharper peak in comparison with the Gaussian distribution.

The values (mean SD) of listed HRV indices for the group of healthy subjects and two groups of patients with CHF are presented in Figure 2. All indices overlap with values reported in the literature for comparable datasets [22,35]. In line with previous studies [22], congestive heart failure leads to a decrease in values of numerous HRV indices. We observe that the power of the HRV in LF and HF frequency bands, the LF/HF ratio, or the scaling coefficients estimated through DFA are affected by the pathological condition. This corresponds to a change in the sympatho-vagal balance, where the CHF leads to a dominance of the parasympathetic branch of the ANS. Reduced DFA scaling coefficients reflect diminished correlations present in HRV time series, reflecting reduced complexity of cardiovascular control mechanisms.