Multifractal Analysis of Hemodynamic Behavior: Intraoperative Instability and Its Pharmacological Manipulation

• Fractal analysis is a nonlinear analytical method that can characterize fluctuating physiological processes

• Multifractal methods characterized intraoperative heart rate and blood pressure fluctuations in four patients undergoing aortic aneurysm repair and the effect of pharmacological intervention on underlying dynamics

• Increased signal complexity during physiological instability suggests recruitment of homeokinetic mechanisms to maintain normal physiology, with more complex fluctuations

• Pharmacological interventions that restore homeostasis or prevent the dynamic response to perturbation reduce fluctuation complexity because there are fewer active restorative homeokinetic processes

PHYSIOLOGIC instability is a common clinical problem both in the critically ill and in patients subjected to operative insult. It may result in abrupt and apparently random fluctuations in hemodynamic parameters with time. In health, homeostasis is traditionally thought of in terms of linear negative feedback loops serving to restore the system to stability. Within this paradigm, the violent and unpredictable swings seen in critical illness are usually considered to arise from a pathologic breakdown of homeostatic mechanisms, resulting in oscillation or sensitivity to external perturbation.

In the past decade there has been increasing interest in nonlinear systems as more realistic ways of describing natural processes. Such systems, where inputs and outputs are not directly proportional, are characterized by a very rich range of behavior that resembles the complexity observed in nature. Depending on the precise nature of the nonlinearity, nonlinear systems may be observed to be highly stable (sometimes with several possible behaviors), exhibit limit-cycle behavior, or be unstable. Instability may arise from extreme and persistent sensitivity to perturbation (popularly known as the “butterfly effect”) and may under certain circumstances lead to complex and apparently unpredictable de novo  behavior; a phenomenon which has come to be known as “chaos.” Such stereotypes provide a unified description that closely parallels what is observed under incipient and established physiologic instability. However, in this model, seemingly random fluctuations may emerge as a direct consequence of the underlying dynamics, rather than as a pathologic breakdown of homeostasis leading to the unwanted transmission of exogenous influences. In a nonlinear model, physiologic equilibrium is in fact achieved through an ensemble of dynamical fluctuating processes, a situation known as homeokinesis. It follows that these fluctuations may contain important information regarding the dynamical state of the system that is otherwise not apparent from traditional observations. Unfortunately, such systems are inherently difficult to study, and recent interest has been fuelled in no small part by advances in mathematical tools and the availability of affordable computing power. Nevertheless, the application of sophisticated analysis techniques has unmasked previously “hidden information” in time-series data of a wide variety of natural processes.1,–,7The observed fluctuations exhibit characteristics of the nonlinear processes that generate them, thus the underlying physical state of the system can be inferred.

Characterization of the underlying homeokinetic state from the fluctuation of physiologic signals may be of clinical importance.8It has been demonstrated that changes in the underlying dynamical state may be detectable before they propagate as instability,9suggesting that bedside black-box nonlinear analysis may be a possibility for tracking changes and providing an early warning of impending physiologic collapse. Furthermore, bedside detection and categorization of different “unstable states” has the potential to guide patient care and gauge response to therapy.

Fractal-based methods have been successfully applied to the study of a wide-range of nonlinear systems, including the prediction of stock markets trends,3weather patterns,7neuronal response in functional magnetic resonance imaging,4geophysics,5behavior of broadband internet traffic,6and highway congestion.9A fractal is a geometric structure constructed from an infinite set of increasingly small subunits. Each subunit is a scaled-down replica of the whole, producing a self-similar property across length scales, i.e. , across many levels of magnification.10Fractal theory is concerned with structural self-similarity over scales of space or time and techniques have been developed to quantify the statistical self-similarity of a structure, recognizing subunits that share such statistical properties across orders of scale.2,10The statistical notion of self-similarity is equally applicable in the time domain, permitting analysis of processes that are self-similar over many time scales.2,10Fractals provide a natural nonlinear metric for the characterization of fluctuating processes.

Early work in the field by Ivanov et al .1demonstrated that the human heart rate, in health, has a fractal temporal structure, and that the fractal properties change in congestive heart failure. Subsequent work by the same and other authors have studied fractality for a number of physiologic processes, including heart rate in age and disease,2stride length in health and chronic neurologic disease,2,11and middle cerebral artery blood flow in subarachnoid hemorrhage.12Importantly they have shown that fractal properties persist throughout episodes of rest and sleep and are not a result of superimposed complex physical or mental activity but instead reflect an emergent property of the underlying dynamical system.

The simplest fractal systems possess a self-similarity described by a single scaling parameter: the Hurst exponent.10Such systems are described as monofractal and several techniques for estimating the Hurst exponent have been developed, including detrended fluctuation analysis and rescaled-range analysis.10,13,14 

More recently,1it has been established that complex dynamical systems may instead result from a spectrum of processes with a range of different scaling parameters. Such systems with multiple scaling behaviors are called multifractal. These arise through several interacting processes, each with different self-similar behaviors acting in concert to produce the overall structure, or a single process whose self-similar statistical properties change within the timeframe under analysis. Parameterization of fractal properties, which can be thought of as a continuum from mono- to multifractal, provides an insight into the behavior and mechanisms of the underlying control systems.2,15 

Physiologic processes typically exhibit multifractal behavior. Rather than having a single scaling parameter, multifractal systems are described by a singularity spectrum. This defines the distribution of scaling parameters termed Hölder exponents (a generalization of the Hurst exponent). The Hölder exponent (h ) can be interpreted in terms of the statistical properties of the time-series, viz .

1. 0 ⩽h < 0.5: Antipersistent behavior. An increase at one time interval is more likely to be followed by a decrease and vice versa .

2. h = 0.5: Uncorrelated random walk. Increases and decreases are equally likely.

3. h > 0.5: Persistent behavior. An increase in one time period is more likely to be followed by a further increase in the next period.

A typical singularity spectrum is shown for illustration (fig. 1A). The spread of Hölder exponents depicted in the spectrum represent the multiple, different self-similar scaling behaviors in the underlying multifractal system. The y -axis charts the Hausdorff dimension, D(h).  The precise mathematical definition of Hausdorff dimension is rather involved; in essence it reveals the relative frequencies of the Hölder exponents, akin to a probability distribution. The Hölder exponent, h^m , with the greatest Hausdorff dimension is the most frequent and hence most dominant scaling behavior. In a manner common to other statistical distributions, the singularity spectrum is parameterized using measures of central-tendency and spread, conventionally h^m  and width-at-half-height (WHH ), the width of the curve at half the maximal Hausdorff dimension (fig. 1A). Such parameterization permits the recording and comparison of multifractal properties within or across datasets. Multifractal techniques can also quantify monofractal scaling behavior, thus providing a very general approach. Relatively monofractal systems are characterized by narrow singularity spectra. In contrast, strongly multifractal systems display a wider distribution of Hölder exponents (fig. 1).

This paper examines the fractal dynamics of heart rate (HR) and mean arterial blood pressure (MAP) in a small cohort of patients undergoing open aortic aneurysm repair. Patients undergoing open abdominal aortic aneurysm repair are subjected to a complex array of surgical and anesthetic perturbations that may overlap temporally and that are complicated to disentangle. Fractal behavior was characterized using a multifractal analysis technique. The application of these techniques to the acutely unstable patient has not previously been undertaken, neither has it been applied to the study of blood pressure physiology.

This study focused on the effect of intraoperative vasoconstrictor boluses on multifractal dynamics. Within this model, the periods before vasoconstrictor administration represent episodes of clinically significant intraoperative hypotension where homeokinetic systems should be maximally active and fluctuation behavior should be particularly complex. Comparison with the behavior immediately after vasoconstriction allows the effect of pharmacological normalization of blood pressure on system dynamical behavior to be studied.

Hemodynamic observational data from four nondiabetic patients undergoing elective open infrarenal aortic aneurysm repair and recruited under a previous study were retrospectively analyzed. The previous aneurysm study recruited a total of seven patients; three patients were excluded because it was impossible to extract blocks of hemodynamic data of sufficient length for analysis. Appropriate ethical review board approval for further analysis of the original, fully anonymous data were obtained (05/Q0108/470; Cambridge Central Research Ethics Committee, Cambridge, United Kingdom). Patients enrolled in the original study underwent written informed consent; further consent was not required for the study extension. A summary of patient demographics is presented in table 1.

Patients were prepared for surgery following usual local practice. Standard anesthetic monitoring was established. Peripheral intravenous access was obtained and the radial artery was cannulated with a 20-gauge catheter for invasive arterial pressure recording. Anesthesia was induced with fentanyl (1–1.5 μg/kg) and propofol (2–3 mg/kg). Following muscle relaxation with atracurium (0.6 mg/kg) the trachea was intubated and the lungs ventilated to normocapnia with a mixture of desflurane, oxygen, and air. A triple-lumen central venous catheter was sited in the right internal jugular vein and used, inter alia , for the intraoperative administration of cardiovascular drugs. Intraoperative analgesia was provided with remifentanil (0.1–0.2 μg⁻¹· kg⁻¹· min⁻¹). Epidural anesthesia was not employed intraoperatively.

In the original study, beat-to-beat cardiac output measurements were made using a LiDCO monitor (LiDCO Ltd., Cambridge, United Kingdom) slaved to the invasive arterial blood pressure recording system throughout the operative period from preintubation to discharge from theater recovery. In this study we only analyze the beat-to-beat records of HR and MAP from the LiDCO device. Surgical and anesthetic events (including the administration of vasoactive drugs) were recorded with a timestamp and a corresponding unique event marker on the LiDCO monitor.

Intraoperative data from each patient record were selected for further analysis. Data points from the connection of the LiDCO monitor up to and including tracheal intubation and subsequent transfer to the operating table were excluded because of the number of exogenous factors affecting arterial line transduction. Data points from aortic cross-clamp onward were also excluded because of the potential for aortic clamping to cause a profound change in cardiovascular physiology. The remaining data points, which were continuous and correspond to the middle of the intraoperative period, were analyzed in full.

Only minor preprocessing of the data were required before multifractal analysis and this was undertaken using custom-written software that provided visualization and data manipulation capabilities. The LiDCO monitor determines for each heartbeat whether the recorded values are “good” or “bad.”“Bad” data typically arises at times the system cannot elicit meaningful values from the arterial transducer, usually when the patient is being moved or the arterial line used for intraoperative blood sampling. Only 0.3–3% of the data were marked as “bad” for each patient and these regions were excluded.

Multifractal analysis of HR and MAP data were undertaken using the wavelet transform modulus maxima (WTMM) technique. An introductory explanation of the technique is given in Appendix 1 with full details of our implementation in Appendix 2. In summary, the data (fig. 2A) is initially subjected to a wavelet transform; a signal-processing technique used to explore structural properties at different length (or time) scales, akin to a mathematical zoom lens. The transform results in a matrix of coefficients (figs. 2B and C), whose maxima encode the self-similar behavior. Signals with self-similar scaling behavior show branching maxima patterns, appearing like an inverted leafless tree that can be seen throughout the pictorial representation, with smaller branches at fine scales and larger trunks at coarse scales. The WTMM technique characterizes the different branching patterns as singularity spectra (fig. 2D), enabling the self-similarity of different signals to be measured and compared.

The WTMM analysis was implemented in Matlab v7.5.0 (MathWorks, Natick, MA) with Wavelet Toolbox v4.1 (MathWorks), and the support tools were coded in the Java programming language, Java v1.6 (Oracle Corporation, Redwood Shores, CA). The Matlab scripts and Java support programs were designed and written in-house and the Matlab code is available for download (see Supplemental Digital Content 1–6, https://links.lww.com/ALN/A875, a collection of files containing the Matlab code, usage instructions, and open-source license agreement). Full details of the design, implementation, and optimization decisions can be found in Appendix 2.

The remaining data were divided into windows of 256 consecutive beats and the HR and MAP were subject to fractal analysis using the WTMM technique with maximum scale 48 beats and partition function fitting over scales 5–32 beats. A definition of these parameters and discussion on their optimization based upon our experimental experience can be found in Appendix 2. The corresponding spectra were parameterized in terms of dominant Hölder exponent, h^m , and WHH , as previously described (see Introduction and fig. 1). The spectral parameters for heart rate (h^mHR , WHH^HR  )  and MAP (h^mMAP , WHH^MAP ) were recorded in operation order and paired with the logbook of operative events for subsequent interpretation and analysis. Episodes of hypotension where α-1 agonist metaraminol (0.5 mg IV bolus) had been administered at the anesthetist's discretion were identified. Multifractal analysis was performed for the periods immediately before and after vasoconstrictor administration.

The WTMM technique is subject to several technically challenging issues regarding wavelet choice and the selection and optimization of parameters. Importantly the particular choice of wavelet and optimization strategies have effect only on estimation of the absolute spectral parameters; provided that these choices remain static throughout the analysis, the trends in spectral parameters (which we are interested in in this study) are preserved. Interested readers are referred to Appendix 2 and Oświęcimka et al.  for further detail.16 

The Hölder exponents, h^mMAP  and h^mHR  before and after metaraminol bolus are represented as mean ± SE. These are further subdivided upon whether the patient received treatment with atropine. A two-tailed, paired Student t  test is used to compare per-patient mean values of h^mMAP  and h^mHR  before and after metaraminol bolus (per-patient mean values are used to eliminate errors because of multiple comparisons). Robust hypothesis testing on h^mHR  to compare patients based upon administration of atropine is difficult, as independence between measurements from the same patient cannot be assured. However, results of a two-tailed paired Wilcoxon signed-rank test (with 1,000-round Monte Carlo jitter modification for tied ranks) is presented under the assumption of independence for reference.

Statistical analysis is similarly presented for WHH^HR  and WHH^MAP . Further details on the quality of the WHH  results and the methodological issues surrounding WHH  estimation are included in the Discussion.

All statistical tests were performed using the R statistical package v2.13.1 (The R Foundation for Statistical Computing, Vienna, Austria) and two-tailed P < 0.05 was considered to be statistically significant.

Between the four patients, a total of nine events of suitable data quality were obtained where no other confounding anesthetic (e.g. , fluid administration) or surgical perturbation was present. Summary statistics for the spectral parameters calculated from these events are detailed in tables 2and 3.

Figure 3A shows the effect of vasopressor boluses on the dominant Hölder exponent for MAP, h^mMAP . Immediately after the bolus, a statistically significant increase in h^mMAP  is seen (table 2), reflecting a change to a more correlated physiologic process with greater long-term memory behavior. During this time the system exhibits a short period of correlated behavior (Hölder exponents more than 1) or near- correlated behavior (Hölder exponents tending toward 1). Subsequently h^mMAP  decreases to a value greater than its starting point, suggesting a more correlated underlying behavior persists. Metaraminol causes almost pure α-1 adrenergic mediated vasoconstriction, changing an under-filled vasculature into a much “tighter” system with fewer homeokinetic mechanisms in action. In other words, less homeokinetic complexity is present. Our observation of altered fractal behavior in response to pharmacological perturbation indicates a profound alteration in the nonlinear character of the ensemble of interacting physiologic homeokinetic mechanisms. Thus, cardiovascular pharmacology allows the clinician not only to manipulate traditional physiologic parameters, such as blood pressure or HR, but also to fundamentally control the stability of the system as a whole by modifying the total activity/number of homeokinetic mechanisms in operation.

The effect of a metaraminol bolus on the dominant Hölder exponent for heart rate, h^mHR  is shown in figure 3B and behaves differently to that for blood pressure. Metaraminol has essentially no direct chronotropic effect and this may explain why the dynamical nature of heart rate homeostasis is unaffected. Statistical nonsignificance is demonstrated in table 2.

Measures of the degree of multifractality for MAP and HR (WHH^MAP  and WHH^HR ) varied throughout the intraoperative period and were altered by metaraminol but without discernable pattern (table 3).

In some cases, it had been necessary for the attending anesthetist to give atropine (600 μg, IV bolus) intraoperatively. Figure 3is marked based upon whether the patient had received atropine, with summary statistics shown in table 2. Atropine tended to make the heart rate fluctuations (fig. 3B) more positively correlated in behavior (higher h^mHR ). This is clearly evident even before metaraminol is given. Atropine antagonizes the background parasympathetic tone and prevents further modulation of parasympathetic effects on the cardiac electrical systems. It could be considered as removing a degree of freedom from the nonlinear heart rate feedback, reducing the homeokinetic complexity and tending to make the system more correlated in its behavior. Following a metaraminol bolus, the dominant correlated behavior is sustained. In contrast, the effect of atropine is not clearly seen in MAP fractality (fig. 3A; identical marking), suggesting that heart rate variation is not a key homeokinetic mechanism in the determination of MAP.

The retrospective analysis performed in this study presents its own challenges, particularly in providing assurance that the effects demonstrated are not artifact. We consider in turn a number of factors that we show do not influence fractal structure.

The effects of ectopic beats and arterial line artifacts on fractal structure were evaluated using a technique similar to Ivanov et al. .1Beat-to-beat HR was converted to interheart beat interval and intervals of greater than three standard deviations were corrected through linear interpolation, including the point either side of the outlier. A similar approach was adopted for mean arterial blood pressure. The interbeat difference in mean arterial blood pressure was calculated and differences of greater than three standard deviations were corrected by linear interpolation. Overall approximately 3 or 4% of data points were interpolated by this method. Comparison of the fractal spectra for preprocessed and unprocessed data sets demonstrated that no significant changes occurred after preprocessing, except in blocks where many errant data points clustered. A number of the 256 data point blocks (12% of total blocks) contained a large proportion of interpolated points (up to 21% points were interpolated) and despite showing multifractal behavior in the unprocessed analysis they demonstrated monofractal behavior after linear interpolation. An unstable patient with frequent marked fluctuations in HR or blood pressure would be heavily interpolated by this technique, removing the fractal structures of interest. We therefore concluded that interpolation at the preprocessing stage was unnecessary as it did not provide a significant increase in the quality of results, neither were artifacts or ectopics found to be responsible for the vast majority of the detected fractal structure.

The analysis of mean arterial blood pressure in ventilated patients is subject to the effects of changes in intrathoracic pressure with ventilation (respiratory swing). We compared multifractal analysis of the beat-to-beat mean arterial blood pressure recording with a detrended version. Blood pressure was detrended to remove respiratory swing by fitting a piecewise cubic hermite interpolating polynomial to the maxima identified through a peak-finding algorithm, and subtracting these polynomials from the original signal. Piecewise cubic hermite interpolating polynomials are cubic splines that are not guaranteed to be continuous in derivatives greater than the first order, thus permitting a tight fit to the data. Subsequent multifractal analysis of detrended MAP demonstrated no significant change in the multifractal spectra. We therefore concluded that respiratory detrending was unnecessary.

Finally, to exclude analytical artifacts, we tested the validity of analyzing blocks of a fixed number of heart beats (256 consecutive beats) with blocks of fixed time duration, because this may have led to inconsistent results when pharmacological interventions with chronotropic effect were used intraoperatively. The patients in figure 3B have heart rates that vary from 31 beats/min to 80 beats/min, thus two 256-beat blocks could vary by as much as 5 min in duration. We repeated our analysis with fixed length blocks of 3.2-min duration. No significant changes were seen, with the h^mHR  response similarly partitioning on prior administration of atropine. Limiting the blocks to a fixed time period includes more or fewer points than previously, hence small differences in the fractal estimation are to be expected. In practice, the length of each block could be fixed at any length. Our choice of 3.2 min was restricted by the timing and clustering of interventions during the operative period and was selected to ensure that the blocks were long enough to encompass the same behavior of interest as blocks in figure 3B.

No obvious difference in behavior was observed in patients who were perioperatively β-blocked. Furthermore, similar behavior was seen following both aortic clamping and unclamping. This may seem rather surprising because it is reasonable that large physiologic consequences might be expected from aortic cross-clamping. In practice, the authors see very few and, at most, only short-lived, traditional hemodynamic consequences in patients undergoing elective aortic aneurysm repair, and hence it is likely that this study is inadequately powered to detect more subtle changes in physiologic behavior.

To our knowledge this is the first application of multifractal analysis to study the dynamics of heart rate and blood pressure under acute intraoperative physiologic perturbation. Goldberger et al.  demonstrated that the resting heart rate dynamics of young, healthy subjects show multifractal structure, suggesting that the classic understanding of “homeostasis” is too simple and underlying control systems are paradoxically “far from equilibrium.”2Periods of physical stress including vigorous exercise and critical illness demonstrate the great adaptability of homeostatic mechanisms; such dynamic scaling behavior may be conferred by the fractal properties of the physiologic feedback cascades. However, the self-similar scaling demonstrated in HR behavior degenerates with advancing age and chronic heart failure. Loss of fractal scaling ability is implicated in the reduced adaptability and consequent susceptibility of the elderly and comorbid at times of physiologic stress.2 

We demonstrate that fractal properties change during the intraoperative period and physiologic mechanisms change fractal properties in two different ways. First, the process can change its dominant scaling (or self-similar) property, becoming more or less persistent in behavior. Second, control systems can become more or less multifractal, hence displaying greater or fewer self-similar behaviors. We postulate that these changes are the result of homeostatic adaption during times of increased stress in critical illness and are modulated by a combination of endogenous and exogenous factors. Furthermore, we provide the first demonstration that physiologic multifractal behavior can be acutely modified by pharmacological intervention.

Situations of intraoperative hypotension represent a homeostatic challenge. In a dynamic model, as the system moves further from equilibrium, more/stronger homeokinetic processes are recruited. Within this paradigm, physiology is maintained but the fluctuations must also increase in complexity. Conversely, restoration of physiology by therapeutic intervention may be expected to reduce fluctuation complexity, as fewer restorative homeokinetic processes are active. Similarly, interventions that limit the possible dynamic response of the system to perturbation (for example, pharmacological vagolysis with an antimuscarinic, which limits the heart's ability to respond with changes in HR) may similarly reduce fluctuation complexity.

A key result is the observation of a fundamental difference in response of h^mMAP  and h^mHR  to vasoconstrictor therapy. Vasoconstriction led to an increase in h^mMAP  (fig. 3A) but not in the case of h^mHR  (fig. 3B). This suggests that vascular tone is an important homeokinetic mechanism in the control of MAP but not of heart rate, which seems reasonable. Conversely, administration of atropine is expected to eliminate HR variability from the ensemble of homeokinetic processes. Figure 3B demonstrates a consistently higher value of h^mHR  for those patients pretreated with atropine. This implies a more correlated and less complex signal, consistent with a reduction in homeokinetic degrees of freedom. This partitioning is not evident for MAP data (fig. 3A), suggesting that pharmacological “paralysis” of HR variability does not much affect control and stability of MAP, which seems intuitively reasonable.

Goldberger et al.  2demonstrated that physiologic processes possess fractal structure and the fractal properties change in chronic disease, reflecting a reduction in homeokinetic repertoire in the failing organism. It seems plausible that the fractal properties of the same systems change or degenerate after acute insults. External perturbations such as changes in mechanical ventilation, surgical insults, and pharmacological interventions can have a marked effect on vascular physiology, and these alone might explain some of the detected changes in fractal structure. Analysis of how fractal properties change in response to a perturbation can also give useful insight into the underlying system. An extreme example occurs when nonlinear systems are on the verge of instability. A small external perturbation can have a destabilizing effect, leading to profound changes in behavior that perhaps in other situations would not occur.

Oświęcimka et al.  17demonstrated that multifractal analysis using short blocks of data produced errors in estimation of the multifractal spectrum tails, hence leading to errors in the calculation of WHH.  The short time intervals used in this study and the inherent inaccuracies in WHH  estimation are the most probable reason why trends are not seen in WHH  during periods of instability and pharmacological manipulation.

Our values of h^mHR  Hölder exponents are rather larger than those reported previously.1,2We attribute this effect to the comparatively short blocks of length 28we analyzed, much shorter than those used in other work.1,16,17This will tend to give larger, more correlated, estimates for h^mHR . As demonstrated by Oświęcimka et al.  17errors in the estimation of h^m  are reasonably preserved with decreasing block length, and despite our absolute values being larger than previous works, we believe the overall trends demonstrated in this study are physiologically meaningful. Our measurements are therefore necessarily distinct from those of others1,2in that they represent short-range scaling behavior and we cannot draw conclusions over longer time scales in our strongly perturbed systems where exogenous influences are likely to be causing significant changes of state. However, we suggest that it is the short-range scaling behavior that is of most clinical relevance in describing hemodynamic instability, which occurs over short time scales.

The authors acknowledge that a retrospective study of this size can only draw limited conclusions but believe that it does provide a basis for future research avenues. We have chosen to study acute changes in intraoperative physiology. However, it seems reasonable to assume that similar considerations should apply to acute, critical illness. Recording of real-time physiologic data are commonplace in the intensive care setting and this suggests the possibility of using the technique not simply as a means of analyzing physiologic behavior, but also as a clinical tool for detecting and characterizing incipient instability and homeostatic activity and breakdown.

The field of nonlinear analysis, of which fractal analysis is only a small part, may provide insight into human physiology beyond that which can be derived directly from classic measures. Traditional approaches to homeostasis have concentrated on the interplay between linear feedback mechanisms. The application of techniques such as fractal dynamics presupposes a shift in paradigm to one in which the physiology is regarded holistically as a nonlinear “black box” and the patient's behavior is defined in terms of an abstract “state space,”18where external perturbations and clinical interventions are shown to move a patient from one state to another. Fractal properties are one measure of this state space and are shown in this study to change in association with external perturbations in a seemingly predictable way, with vasopressor agents and atropine producing more correlated and less complex fractal properties of blood pressure and HR control mechanisms, respectively. We suggest that studying the effect of perturbations to fractal systems is a valid and powerful technique for analyzing the underlying dynamical properties, giving information about current state and potential future behavior. Complex behavior is an emergent property of both the perturbing inputs and their interplay with nonlinearities: Many clinicians involved in the management of critically ill patients will recall patients who have responded unexpectedly to an intervention that should have improved their illness. An understanding of the state space and the ability in real-time to determine where in the state space a given patient lies may provide a powerful tool in guiding clinicians to the most appropriate package of interventions to restore the patient to a stable state.

The authors thank Professor J. Suckling, Ph.D., Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom, for helpful discussions and suggestions.

The wavelet transform of a signal is a mathematical method that permits exploration of its properties at different length (or time) scales, akin to a mathematical zoom lens. The transform calls upon an underlying function, a wavelet, that itself can be thought to behave like a mathematical measuring tape. Many different types of wavelets exist and the particular choice of wavelet is application-dependent. The wavelet most frequently used in physiologic applications is the third-order derivative of the Gaussian.

In a similar fashion to the Fourier transform, which transforms a signal into frequency components, the wavelet transform decomposes a signal into its constituent wavelets. Analogously to the Fourier transform, the wavelet transform represents the synthesis of the original signal from superimposed copies of the wavelet varying only in scale and position. It is this decomposition that provides insight into the properties of hidden structures such as self-similarity within the original signal across length or time scales.

The wavelet transform applied to a Cantor set (a mathematically defined data series that is apparently irregular but has known self-similarity properties, fig. 4A) is shown as a two-dimensional map (fig. 4B) and three-dimensional surface (fig. 4C), as a pedagogical example. The transform reveals an underlying structure within the signal. Branching patterns akin to an inverted leafless tree can be seen throughout the pictorial representation, with smaller branches at finer scales and larger trunks at coarser scales. The recursive, branching structure is a representation of the self-similarity and it is this pattern that encodes the hidden fractal information of interest. The wavelet transform modulus maxima technique is one way to quantify and characterize the branching pattern, allowing the self-similarity of different signals to be measured and compared. In the absence of self-similarity (fractality) the branching pattern is lost.

Characterization of the branching pattern is performed through the detection of maxima (peaks) in the surface plot (fig. 4C). Subsequent computational analysis of the maxima across length (or time) scales produces a singularity spectrum. The singularity spectrum for the Cantor set is a multifractal spectrum with h^m  centered at 0.5 (fig. 4D).

A step-by-step guide detailing the implementation of the wavelet transform modulus maxima is contained below:

1. Represent the signal or data of interest as a continuous vector f[Combining Macron] , e.g. , intraheart beat intervals.

2. To minimize edge artifacts, pad the data f[Combining Macron]  of size n  with a further n  values; with n/2  pad values at the start of the data and n/2  pad values at the end of the data. The pad at the start should take the value f[Combining Macron]  (1) and at the end of the data the value f[Combining Macron]  (n ).

3. Calculate the absolute wavelet transform W[Combining Macron]^a  (f[Combining Macron] ,s ) by convolution of the vector with a third-order Gaussian wavelet over scales s = 1…scale_max . The parameter scale_max  is described later.

   The third-order Gaussian wavelet is defined by:

   The wavelet transform acts as a mathematical “zoom lens,” analyzing the signal according to length scale. The choice of wavelet derivative is important. A wavelet derivative of order n  will remove all polynomial trends from the data of degree n-1 , thus a third-order Gaussian kernel will remove all polynomial trends up to quadratic degree. A third-order Gaussian is sufficient for most purposes.

4. Remove the padding from the beginning and end of the coefficient matrix. The wavelet transform coefficient matrix is smoothed before the next step; a suitable approach is described later.

5. The set L  of individual lines of local maxima, lⁱ , are found by iterative search through W[Combining Macron]^a , and represent the paths along the peaks in the wavelet surface plot, figure 2C. Each maxima line, lⁱ εL , is itself a vector containing the point at each scale that lies on the maxima line.

6. The supremum is found for every point on each maxima line, defined in pseudocode as:

   In words, the supremum of the line at scale s is the largest wavelet transform coefficient that lies on the maxima line from scale 1 up to and including the coefficient at scale s . Thus the supremums along each line form a monotonically increasing sequence.

7. Define a partition function over scales s = 1…scale_max  and lines lⁱ εL  as:

   and construct a double-log plot of Z  (q ,s ) against scale s  for each order q = qmin…qmax  in step sizes of qstep , resulting in one line on the plot for each q  value. The partitioning orders qmin , qmax , and the step size qstep  are described later.

8. The self-similarity exponents are derived from the linear relation:

   where the gradient of the straight lines for individual q -lines in the double-log plot defines τ(q ). Straight line fitting only occurs between the scale ranges fit_min  to fit_max , discussed below.

9. Compute the singularity spectrum via  a Legendre transform:

   and plot D(h ) against h .

10. Fit a polynomial to the plot of D  (h ) against h  (using least squares regression) and find the maxima h^m  and width-at-half-height using analytical methods. The choice of polynomial order is discussed below.

The choice of tuning parameters and the design of key algorithm routines is critical to producing an automated implementation of the algorithm.

The continuous wavelet transform function introduces numerical noise into the coefficient matrix, particularly at small scales, which complicates the discovery of maxima lines. The noise is an artifact resulting from the discrete representation of the continuous wavelet function and its convolution with the quantized, discrete representation of the continuous physiologic signal. To overcome difficulties in maxima finding, a triangular smoothing function is applied to the coefficient matrix along each scale, with the width of the smoothing proportional to the scale, i.e. , a wider smoothing function at larger (coarser) scales.

The scale_max  parameter dictates the largest scale over which the algorithm looks for self-similar behavior. This is necessarily smaller than the amount of data; it would not be worthwhile searching for self-similarity between two different 3-min blocks in a sample of only 4-min duration. Experimentally an upper limit of scale_max = 48 beats was found to be realistic when working with data of length 256 beats. Scales greater than 48 beats added little to the analysis other than extra noise, and in turn resulted in a longer computation time.

The discovery of maxima lines within the wavelet coefficient matrix is a difficult task because of coefficient noise (even after the coefficient smoothing described above). Ideally the set of all possible connecting maxima lines would be determined and then the solution with the smallest residual distance between line segments chosen. Unfortunately this is numerically intractable. A close approximation is instead found using a breadth-first algorithm. The algorithm starts at the coarsest scale and for each maxima at the coarse scale computes the distances to the maxima at the next finest scale. The distances are then ordered (shortest first) and the maxima with the shortest distances are joined, ensuring that branching only occurs from coarser to finer scales. This procedure is iterated over decreasing scales until the finest scale is reached. To finish, some cleansing is performed. First, any maxima lines that do not run from the finest to the coarsest scale are analyzed. If any of these lines start and terminate at adjacent scales and are separated by no more than a predefined distance, they are joined to make a longer maxima line. Empirically this identifies the majority of maxima lines that are missed by the initial approximation technique. Second, lines that do not start at the finest scale are pruned from the set of maxima lines, as all singularities must be connected to the finest scale.

A key part of the algorithm is the calculation of the partition function in step 7. To ensure the resulting singularity spectrum is well defined, a large number of moments of order q  are calculated at the partitioning stage. Many implementations of the wavelet transform modulus maxima technique use a range of q 's from −5 to 5 (with a step-size of 1); here a range from q_min = −8 to q_max = +8 with a step-size q_step = 0.15842 (resulting in 101 steps between −8 and 8) is adopted. The use of smoothing and other techniques, which minimize the effect of noise, resulted in an implementation that could operate over such a large number of moments.

Occasionally it was necessary to limit the range of moments used. Monofractal datasets that result in very narrow singularity spectra required q_min  to be increased to 1. Partitioning over these datasets produced very small numbers after raising each partition to a negative power; the numbers were sufficiently small that they could not be accurately represented by the computer's floating-point mechanisms. The errors introduced by the inaccurate representation were large compared with the values being stored, resulting in errors in the singularity spectrum, producing a nonsense spectrum. The restriction of q_min  to 1 solved this issue, but it results in only half of the narrow monofractal spectrum being calculated under these conditions (fig. 1C).

The double-log plot of Z  (q, s ) against scale s  results in 101 different curves (one for each moment q ) and each of these curves is linear over a limited region of interest. The location of the linear region varies and is dependent upon scale_max  and the size of the original dataset. Analyses for a block length of 256 beats using the data in this study revealed that the linear region fell in the range of scales fit_min = 5 to fit_max = 32 beats. Linear regression lines were fitted for each of the 101 different curves in the double-log plot between fit_min  and fit_max. 

A sixth-order polynomial was fitted to points defining the singularity spectra; this provided the best fit over a wide range of resulting spectra. On the rare occasion that an incorrect fit was produced by the regression algorithm, this was identified by eye and corrected before calculation of the spectral parameters. The spectral parameters (h^m  and width-at-half-height) were calculated directly from the polynomial regression equation using analytical techniques.