Cerebral Flow Autoregulation (CFA) is the dynamic process by which cerebral blood flow is maintained within physiologically acceptable bounds during fluctuations of cerebral perfusion pressure. The distinction is made with “static” flow autoregulation under steady-state conditions of perfusion pressure, described by the celebrated “autoregulatory curve” with a homeostatic plateau. This paper studies the dynamic CFA during changes in perfusion pressure, which attains critical clinical importance in patients with stroke, traumatic brain injury and neurodegenerative disease with a cerebrovascular component. Mathematical and computational models have been used to advance our quantitative understanding of dynamic CFA and to elucidate the underlying physiological mechanisms by analyzing the relation between beat-to-beat data of mean arterial blood pressure (viewed as input) and mean cerebral blood flow velocity(viewed as output) of a putative CFA system. Although previous studies have shown that the dynamic CFA process is nonlinear, most modeling studies to date have been linear. It has also been shown that blood CO2 tension affects the CFA process. This paper presents a nonlinear modeling methodology that includes the dynamic effects of CO2 tension (or its surrogate, end-tidal CO2) as a second input and quantifies CFA from short data-records of healthy human subjects by use of the modeling concept of Principal Dynamic Modes (PDMs). The PDMs improve the robustness of the obtained nonlinear models and facilitate their physiological interpretation. The results demonstrate the importance of including the CO2 input in the dynamic CFA study and the utility of nonlinear models under hypercapnic or hypocapnic conditions.
The quantitative study of cerebral haemodynamics has been confounded by the fact that cerebral blood flow variations depend on multiple physiological factors. Variations in cerebral perfusion pressure (CPP) are viewed as a key factor affecting variations in cerebral blood flow, although many other physiological variables influence the latter. Among them, the effects of blood CO2 tension have been shown to be significant. Thus, we have included in our studies the dynamic effects of variations in blood CO2 tension (or its surrogate, the end-tidal CO2) as another variable modulating dynamic cerebral flow autoregulation (CFA). The latter is defined as the physiological process by which variations in cerebral blood flow are kept within physiologically acceptable bounds in the presence of variations in CPP (or systemic arterial blood pressure).
Since CFA is of vital importance for survival and the proper functioning of the brain, it has received considerable attention for many years starting with the pioneering work of Lassen [1Aaslid R, Lindegaard K F, Sorteberg W, and Nornes H, "“Cerebral autoregulation dynamics in humans,”", Stroke, vol. 20, pp. 45-52, 1989.] and Fog [2Bellapart J, and Fraser J F, "“Transcranial Doppler assessment of cerebral autoregulation,”", Ultrasound Med. Bio, vol. 35, pp. 883-893, 2009.]. Head injury and subarachnoid haemorrhage can impair CFA and worsen neurological outcome [3Busse R, Trogisch G, and Bassenge E, "“The role of endothelium in the control of vascular tone,”", Basic Res. Cardiol, vol. 80, pp. 475-490, 1985.-7Czosnyka M, Brady K, Reinhard M, Smielewski P, and Steiner L, "“Monitoring of cerebrovascular autoregulation: facts, myths, and missing links,”", Neurocrit. Care, vol. 10, pp. 373-386, 2009.]. Impaired CFA may also impede recovery from acute stroke [8Diehl R R, Linden D, Lucke D, and Berlit P, "“Spontaneous blood pressure oscillations and cerebral autoregulation,”", Clin. Auton. Res, vol. 8, pp. 7-12, 1998., 9Ekstrom-Jodal B, Haggendal E, Linder L E, and Nilsson N J, " “Cerebral blood flow autoregulation at high arterial pressures and different levels of carbon dioxide tension in dogs”", Eur. Neurol, vol. 6, pp. 6-10, 1971.] and lead to brain edema or hypertensive encephalopathy when CPP is persistently high, and hypoperfusion or cerebral ischemia when CPP is persistently low [10Enevoldsen E M, and Jensen F T, "“Autoregulation and CO2 responses of cerebral blood flow in patients with acute severe head injury,”", J. Neurosurg, vol. 48, pp. 689-703, 1978.]. Furthermore, changes in cerebrovascular resistance and compliance (i.e. frequency-dependent impedance or, its inverse, admittance) have been associated with hypertension, diabetes, atherosclerosis, dyslipidemia [8Diehl R R, Linden D, Lucke D, and Berlit P, "“Spontaneous blood pressure oscillations and cerebral autoregulation,”", Clin. Auton. Res, vol. 8, pp. 7-12, 1998., 11Fog M, "“Autoregulation of cerebral blood flow and its abolition by local hypoxia and-or trauma,”", Scand. J. Clin. Lab. Invest. Suppl, vol. 102, 1968.-14Giller C A, and Mueller M, "“Linearity and nonlinearity in cerebral hemodynamics,”", Med. Eng. Phys, vol. 25, pp. 633-646, 2003.] and other cardiovascular risk factors, including ageing and neuro-degenerative processes [15Haubrich C, Kruska W, Diehl R R, Moller-Hartmann W, and Klotzsch C, "“Dynamic autoregulation testing in patients with middle cerebral artery stenosis,”", Stroke, vol. 34, pp. 1881-1885, 2003., 16Hu K, Peng C K, Czosnyka M, Zhao P, and Novak V, " “Nonlinear assessment of cerebral autoregulation from spontaneous blood pressure and cerebral blood flow fluctuations,”", Cardiovasc. Eng, vol. 8, pp. 60-71, 2008.]. The extent to which these cerebrovascular changes are related to dysfunction of CFA remains an open question of great clinical importance [13Giller C A, Bowman B, Dyer H, Mootz L, and Krippner W, " “Cerebral arterial diameters during changes in blood pressure and carbon dioxide during craniotomy,”", Neurosurgery, vol. 32, pp. 737-741, 1993.].
CFA generally depends on the cerebrovascular characteristics and the proper function of the neurovascular unit, including intrinsic and autonomic neural activities [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983., 18Immink R V, van den Born B J, van Montfrans G A, Koopmans R P, Karemaker J M, and van Lieshout J J, "“Impaired cerebral autoregulation in patients with malignant hypertension,”", Circulation, vol. 110, pp. 2241-2245, 2004.]. Its assessment in human subjects requires direct or indirect measurements of CPP, which are often approximated by measurements of systemic arterial blood pressure and localized cerebral blood flow velocity using various imaging or Doppler modalities [3Busse R, Trogisch G, and Bassenge E, "“The role of endothelium in the control of vascular tone,”", Basic Res. Cardiol, vol. 80, pp. 475-490, 1985., 19Kalaria R N, "“Cerebral vessels in ageing and Alzheimer’s Disease,”", Pharmacol. Ther, vol. 72, pp. 193-214, 1996.]. For example, beat-to-beat changes in arterial pressure can be measured non-invasively via finger photoplethysmography that can be related to CPP by reference to intracranial pressure (if available). The cerebral blood flow velocity is typically measured at the middle cerebral artery via trans-cranial Doppler [20Kontos H A, Wei E P, Raper A J, and Patterson J L Jr, "“Local mechanism of CO2 action of cat pial arterioles,”", Stroke, vol. 8, pp. 226-229, 1977.-23Lassen N A, "“Cerebral blood flow and oxygen consumption in Man,”", Physiol. Rev, vol. 39, pp. 183-238, 1959.] and can be related to changes in flow by reference to the arterial lumen cross-sectional area [24Lindegaard K F, Lundar T, Wiberg J, Sjoberg D, Aaslid R, and Nornes H, "“Variations in middle cerebral artery blood flow investigated with noninvasive transcranial blood velocity measurements,”", Stroke, vol. 18, pp. 1025-1030, 1987.]. Assuming that the latter remains relatively constant through time, the measured changes in flow velocity are accepted as a surrogate for changes in flow [25Marmarelis P Z, and Marmarelis V Z, "Plenum, New York, 1978. Russian translation: Mir Press: Moscow, 1981", In: Analysis of Physiological Systems: The White-Noise Approach, Chinese translation: Academy of Sciences Press: Beijing, 1990., 26Marmarelis V Z, "“Coherence and apparent transfer function measurements for nonlinear physiological systems,”", Ann. Biomed. Eng, vol. 16, pp. 143-157, 1988.]. In the steady-state, the concept of “static” CFA has emerged in the form of the “autoregulatory curve” exhibiting a flow plateau demarcated by two flexion points (upper and lower limit of the plateau) [27Marmarelis V Z, "“Identification of nonlinear biological systems using Laguerre expansions of kernels,”", Ann. Biomed. Eng, vol. 21, pp. 573-589, 1993.]. However, this static curve reveals no temporal information on how blood flow responds to changes in arterial blood pressure. Thus, recent interest has focused on “dynamic” CFA that refers to the relation between continuous variations in blood pressure and flow. Several other physiological variables have been known to affect the cerebral vasculature directly – e.g. CO2 tension and hypoxia [28Marmarelis V Z, "“Modeling methodology for nonlinear physiological systems,”", Ann. Biomed. Eng, vol. 25, pp. 239-251, 1997., 29Marmarelis V Z, Chon K H, Chen Y M, Holstein-Rathlou N H, and Marsh D J, "“Nonlinear analysis of renal autoregulation under broadband forcing conditions,”", Ann. Biomed. Eng, vol. 21, pp. 591-603, 1993.], adenosine, hydrogen or potassium ions [30Marmarelis V Z, Chon K H, Holstein-Rathlou N H, and Marsh D J, "“Nonlinear analysis of renal autoregulation in rats using principal dynamic modes,”", Ann. Biomed. Eng, vol. 27, pp. 23-31, 1999.]. The vascular regulatory mechanisms of CFA are rather complicated and include the effects of the myogenic mechanism [31Marmarelis V Z, Nonlinear Dynamic Modeling of Physiological Systems, Wiley-Interscience & IEEE Press: Hoboken NJ, 2004.], flow-mediated endothelial mechanisms [32McCulloch J, and Edvinsson L, "“Cerebrovascular smooth muscle reactivity: A critical appraisal of in vitro and in situ techniques,”", J. Cereb. Blood Flow Metab, vol. 4, pp. 129-139, 1984.] and perivascular innervation [33Mitsis G D, Zhang R, Levine B D, and Marmarelis V Z, " “Modeling of nonlinear physiological systems with fast and slow dynamics II: Application to cerebral autoregulation,”", Ann. Biomed. Eng, vol. 30, pp. 555-565, 2002.-37Mraovitch S, Iadecola C, Ruggiero D A, and Reis D J, "“Widespread reductions in cerebral blood flow and metabolism elicited by electrical stimulation of the parabrachial nucleus in rat,”", Brain Res, vol. 341, pp. 283-296, 1985.]. A number of “autoregulatory indices” have been proposed in a clinical context for the assessment of CFA in human subjects that result either from static measures of cerebral flow and arterial pressure [16Hu K, Peng C K, Czosnyka M, Zhao P, and Novak V, " “Nonlinear assessment of cerebral autoregulation from spontaneous blood pressure and cerebral blood flow fluctuations,”", Cardiovasc. Eng, vol. 8, pp. 60-71, 2008., 38Muizelaar J P, Ward J D, Marmarou A, Newlon P G, and Wachi A, "“Cerebral blood flow and metabolism in severely head-injured children. Part 2: Autoregulation,”", J. Neurosurg, vol. 71, pp. 72-76, 1989., 39Newell D W, and Aaslid R, Transcranial Doppler, Raven Press: New York, 1992.] or dynamic measures [20Kontos H A, Wei E P, Raper A J, and Patterson J L Jr, "“Local mechanism of CO2 action of cat pial arterioles,”", Stroke, vol. 8, pp. 226-229, 1977., 40Newell D W, Aaslid R, Lam A, Mayberg T S, and Winn H R, " “Comparison of flow and velocity during dynamic autoregulation testing in humans”", Stroke, vol. 25, pp. 793-797, 1994.-42Panerai R B, Simpson D M, Deverson S T, Mahony P, Hayes P, and Evans D H, "“Multivariate dynamic analysis of cerebral blood flow regulation in humans,”", IEEE Trans. Biomed. Eng, vol. 47, pp. 419-423, 2000.].
Several studies of dynamic CFA have obtained the frequency-dependent cerebrovascular “impedance” (or its inverse “admittance”) through Fourier analysis of the time-series of beat-to-beat changes in mean arterial blood pressure (MABP) and mean cerebral blood flow velocity (MCBFV) over each R-R interval [41Panerai R B, Dawson S L, and Potter J F, "“Linear and nonlinear analysis of human dynamic cerebral autoregulation,”", Am. J. Physiol, vol. 277, pp. H1089-H1099, 1999., 43Panerai R B, Eames P J, and Potter J F, "“Variability of timedomain indices of dynamic cerebral autoregulation,”", Physiol. Meas, vol. 24, pp. 367-381, 2003., 44Panerai R B, Kerins V, Fan L, Yeoman P M, Hope T, and Evans D H, "“Association between dynamic cerebral autoregulation and mortality in severe head injury,”", Br. J. Neurosurg, vol. 18, no. 5, pp. 471-479, 2004.]. The model used in those studies is a putative linear system defined by MABP as the “input” and MCBFV as the “output”. Its admittance function showed high gain in the 0.07-0.15 Hz and 0.2-0.3 Hz frequency ranges and a reduction of gain for frequencies below 0.07 Hz, a fact interpreted as the effect of CFA over this frequency range [43Panerai R B, Eames P J, and Potter J F, "“Variability of timedomain indices of dynamic cerebral autoregulation,”", Physiol. Meas, vol. 24, pp. 367-381, 2003.]. It was also observed that the coherence function between changes in MABP and MCBFV was reduced below 0.07 Hz, a fact that was interpreted as indicative of intrinsic nonlinearities and/or nonstationarities in this frequency range [43Panerai R B, Eames P J, and Potter J F, "“Variability of timedomain indices of dynamic cerebral autoregulation,”", Physiol. Meas, vol. 24, pp. 367-381, 2003.]. Time-domain analysis using cross-correlation functions [45Panerai R B, Moody M, Eames P J, and Potter J F, "“Cerebral blood flow velocity during mental activation: interpretation with different models of the passive pressure-velocity relationship,”", J. Appl. Physiol, vol. 99, pp. 2352-2362, 2005.], Auto-Regressive Moving-Average models [46Panerai R B, "“Cerebral autoregulation: from models to clinical applications,”", Cardiovasc. Eng, vol. 8, pp. 42-59, 2008.] and Impulse Response Functions have also provided some useful insights [47Panerai R B, "“Transcranial Doppler for evaluation of cerebral autoregulation,”", Clin. Autonomic Res, vol. 19, pp. 197-211, 2009.].
The observation of low linear coherence and the fact that the static “autoregulatory curve” is nonlinear, provided the motivation for nonlinear modeling studies of dynamic CFA using the nonparametric input-output Volterra approach [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983., 47Panerai R B, "“Transcranial Doppler for evaluation of cerebral autoregulation,”", Clin. Autonomic Res, vol. 19, pp. 197-211, 2009.-50Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, and Hetzel A, "“Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations,”", Stroke, vol. 36, pp. 1684-1689, 2005.]. These studies quantified the intrinsic nonlinearities of CFA and showed that they reside primarily in the frequency range below 0.07 Hz, where reduced linear coherence was previously observed. The studies have also shown that an additional protagonist of CFA is the CO2 tension in the blood, measured indirectly through end-tidal CO2 (ETCO2), which modulates the pressure-flow relationship, especially in the low-frequency range from 0.02 to 0.08 Hz [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983., 49Reinhard M, Roth M, Muller T, Guschlbauer B, Timmer J, Czosnyka M, and Hetzel A, "“Effect of carotid endarterectomy or stenting on impairment of dynamic cerebral autoregulation,”", Stroke, vol. 35, pp. 1381-1387, 2004.]. The modeling of this modulatory action represents a quantification of CO2 vasomotor reactivity and its interactions with vascular impedance. The characteristics of these dynamic interactions were shown to be altered by orthostatic stress (simulated with low body negative pressure) and pharmaceutical interventions affecting autonomic neural activity [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983., 50Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, and Hetzel A, "“Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations,”", Stroke, vol. 36, pp. 1684-1689, 2005.]. It must be noted that these findings were made possible by the extension of Volterra-type dynamic nonlinear modeling to two inputs (MABP and ETCO2). These dynamic nonlinear models have the potential to explore whether the celebrated homeostatic mechanism of static CFA represented by the nonlinear “autoregulatory curve” can be extended to the broader notion of a “homeodynamic” mechanism pertinent to dynamic CFA. The ultimate goal is to elucidate the underlying physiological mechanisms in order to advance our understanding of cerebral blood flow autoregulation and improve clinical diagnostic and therapeutic procedures.
It is now evident that the underlying physiological complexity of CFA requires dynamic nonlinear modeling methods in a multivariate context, using at least two inputs: MABP and ETCO2. The key notion is that the dynamic interactions between the frequency-dependent vascular impedance/admittance and the CO2 vasomotor reactivity can be quantified through these nonlinear models. This represents a significant advance in our understanding and quantification of the CFA system – one that will naturally require a long maturation process – but it also presents a formidable challenge.
This paper presents a recent variant of efficient Volterra modeling that can be effective for this purpose. This variant utilizes the concept of Principal Dynamic Modes (PDMs) which has been developed in recent years and applied successfully to various physiological domains [51Seylaz J, Hara H, Pinard E, Mraovitch S, MacKenzie E T, and Edvinsson L, "“Effect of stimulation of the sphenopalatine ganglion on cortical blood flow in the rat,”", J. Cereb. Blood Flow Metab, vol. 8, pp. 875-878, 1988.]. The criterion for comparing the PDM-based approach with other approaches used before, is two-fold: (a) the predictive capability of the respective models, provided that the number of free parameters in each model is not more than 20% of the number of data samples (to prevent over-fitting); and (b) the interpretability of the obtained models. The PDM-based modeling methodology is effective with short data-records, thus suitable for nonstationarities in the time-scale of a few minutes (as in the case of CFA). It was applied to beat-to-beat time-series of MABP and MCBFV and breath-to-breath ETCO2 data from 12 normal subjects that were recorded over 6-min. The obtained models allow examination of the intrinsic nonlinearities of the CFA system in a manner that is expected to advance our understanding of the fundamental regulatory aspects of cerebral haemodynamics, including eventually the mechanisms of endothelial and neural processes, as well as the neuronal-vascular coupling relationship. The PDM-based model allows the reliable assessment of the effects of various physiological or pathophysiological conditions on CFA in a manner not possible heretofore. Thus, it has the potential to advance the clinical care of patients with cerebrovascular disease.
Twelve healthy subjects with mean (standard deviation) age of 37 (9) years voluntarily participated in this study and signed the Informed Consent Form that has been approved by the Institutional Review Boards of the University of Texas Southwestern Medical Center and Texas Health Presbyterian Hospital Dallas, where the data were collected as continuous recordings at the Institute for Exercise and Environmental Medicine. The arterial blood pressure was measured continuously and noninvasively with finger photoplethysmography (Finapres, Ohmeda, Colorado) and the cerebral blood flow velocity was measured in the Middle Cerebral Artery using a 2 MHz trans-cranial Doppler (TCD) probe (Multiflow, DWL, Germany) placed over the temporal window and fixed at constant angle with a custom-made holder [20Kontos H A, Wei E P, Raper A J, and Patterson J L Jr, "“Local mechanism of CO2 action of cat pial arterioles,”", Stroke, vol. 8, pp. 226-229, 1977., 23Lassen N A, "“Cerebral blood flow and oxygen consumption in Man,”", Physiol. Rev, vol. 39, pp. 183-238, 1959., 43Panerai R B, Eames P J, and Potter J F, "“Variability of timedomain indices of dynamic cerebral autoregulation,”", Physiol. Meas, vol. 24, pp. 367-381, 2003.]. The heart rate was monitored by a 12-lead ECG and the end-tidal CO2 tension was obtained via a nasal cannula using a mass spectrometer (Marquette Electronics). All experiments were performed in the morning in a quiet, environmentally controlled laboratory. After 20 minutes of supine rest, 6 minutes of recordings were made under resting conditions in supine position. These clinical measurements are reliable, reproducible, non-invasive, inexpensive, safe and comfortable for the subjects.
Continuous recordings of data were reduced to beat-to-beat time-series data by averaging the signals over each RR interval and re-sampling the beat-to-beat values at 2 Hz. The re-sampled data were subsequently de-meaned and de-trended by subtracting the 2-min non-causal moving-average. Finally, we clip outliers in the de-trended signal that are more than 3 root-mean-square values away from the zero mean. The latter is done in order to protect the modeling results from the effects of occasional outliers due to measurement errors. Fig. (1) shows illustrative preprocessed data, as well as the subtracted moving-average, for subject #10. It is evident that very low frequency variations (below 0.01 Hz) are present in the data but they will not be analyzed in this study, because dynamic CA has been traditionally studied in the frequency range 0.01-0.50 Hz.
Fig. (2) Schematic of the PDM-based model of the CFA system with three global PDMs and two pairs of cross-PDMs. |
In this paper, we employ the novel concept of Principal Dynamic Modes (PDMs) to obtain reliable dynamic nonlinear models of the causal relationship between two input signals, mean arterial blood pressure (MABP) and end-tidal CO2 tension (ETCO2) time-series and the output signal of mean cerebral blood flow velocity (MCBFV). The preprocessed time-series data over 6 min from 12 subjects in supine resting position were preprocessed as described above to remove low-frequency trends and clip extreme values due to experimental artifacts. The utility of the PDMs is found in the fact that they make the models more compact and allow their accurate and robust estimation from short data-records. The use of PDMs also facilitates the physiological interpretation of the obtained model, as will be discussed later. We briefly outline below the Volterra modeling approach and the proposed PDM-based variant of this approach. For the many mathematical and technical details of Volterra-type modeling, the reader is referred to the monograph [51Seylaz J, Hara H, Pinard E, Mraovitch S, MacKenzie E T, and Edvinsson L, "“Effect of stimulation of the sphenopalatine ganglion on cortical blood flow in the rat,”", J. Cereb. Blood Flow Metab, vol. 8, pp. 875-878, 1988.]. Some background information is given in Appendix I.
In the Volterra modeling approach, the output is expressed in terms of the input(s) by means of a functional expansion that represents the hierarchical nonlinear interactions among the inputs and among different parts of the input epochs as they affect the output (i.e. a 2^{nd} order Volterra functional between two inputs represents the quantitative pattern by which the two inputs interact in order to influence the output). The specific pattern of these nonlinear interactions for a given system is codified by the respective Volterra “kernels” (i.e. the 2^{nd} order self-kernel codifies the interactions between different parts of the epoch of the respective input and the 2^{nd} order cross-kernel codifies the interactions between two inputs). These kernels convolve the respective inputs to form the corresponding Volterra functional. The 1^{st} order kernel quantifies the linear dynamics of the system and the respective 1^{st} order functional represents the linear component of the Volterra model prediction. Thus, a linear model or a nonlinear model (incorporating a hierarchy of nonlinear functional defined by the respective kernels) can be built in each case as appropriate. It has been proven mathematically that the Volterra formulation is applicable to (almost) all nonlinear causal relationships between arbitrary input and output signals, with the exception of chaotic processes and non-dissipating nonlinear oscillations [51Seylaz J, Hara H, Pinard E, Mraovitch S, MacKenzie E T, and Edvinsson L, "“Effect of stimulation of the sphenopalatine ganglion on cortical blood flow in the rat,”", J. Cereb. Blood Flow Metab, vol. 8, pp. 875-878, 1988.-54Tiecks F P, Lam A M, Aaslid R, and Newell D W, " “Comparison of static and dynamic cerebral autoregulation measurements,”", Stroke, vol. 26, pp. 1014-1019, 1995.]. The importance of the latter notwithstanding, the Volterra model is general and rigorous formulation for all other cases which constitute the vast majority.
The conceptual and mathematical/computational appeal of the Volterra formulation is evident. However, practical problems arise from the fact that the high dimensionality of the high-order kernels makes their estimation difficult and their physiological interpretation daunting. For this reason, most applications to date have been limited to 2^{nd} order, with only a few attempts for 3^{rd} order modeling, including our previous work on renal autoregulation [55Tiecks F P, Lam A M, Matta B F, Strebel S, Douville C, and Newell D W, "“Effects of the valsalva maneuver on cerebral circulation in healthy adults. A transcranial Doppler Study,”", Stroke, vol. 26, pp. 1386-1392, 1995., 56van Beek A H, Claassen J A, Rikkert M G, and Jansen R W, " “Cerebral autoregulation: overview of current concepts and methodology with special focus on the elderly,”", J. Cereb. Blood Flow Metab, vol. 28, pp. 1071-1085, 2008.]. To mitigate this practical limitation, we have utilized kernel expansions on the Laguerre basis that simplify the kernel estimation problem [51Seylaz J, Hara H, Pinard E, Mraovitch S, MacKenzie E T, and Edvinsson L, "“Effect of stimulation of the sphenopalatine ganglion on cortical blood flow in the rat,”", J. Cereb. Blood Flow Metab, vol. 8, pp. 875-878, 1988., 57Volterr V, Theory of Functionals and of Integral and Integro- Differential Equations, Dover Publications: New York, 1930., 58Wiener N, Nonlinear Problems in Random Theory, MIT Press: Cambridge, Mass, 1958.]. Although the use of kernel expansions resulted in spectacular estimation benefits, it does not remove the “curse of dimensionality” associated with the multi-dimensional structure of high-order kernels. The key mathematical relations of this modeling approach are summarized in Appendix I.
We seek to overcome this fundamental practical limitation with the introduction of the PDM concept, which aims at identifying the minimum set of "basis functions" (that are distinct and characteristic for each system) capable of representing adequately the system dynamics (i.e. provide satisfactory expansions of the system kernels) and, at the same time, minimize the importance of cross-terms by separating the system nonlinearity into appropriate individual nonlinear functions cascaded with the PDMs to form collectively the system output. This is equivalent to splitting the multi-input static nonlinearity of the Wiener-Bose model, which was shown to be an equivalent formulation of the general Volterra model [52Steiner L A, Czosnyka M, Piechnik S K, Smielewski P, Chatfield D, Menon D K, and Pickard J D, "“Continuous monitoring of cerebrovascular pressure reactivity allows determination of optimal cerebral perfusion pressure in patients with traumatic brain injury,”", Crit. Care Med, vol. 30, pp. 733-738, 2002., 54Tiecks F P, Lam A M, Aaslid R, and Newell D W, " “Comparison of static and dynamic cerebral autoregulation measurements,”", Stroke, vol. 26, pp. 1014-1019, 1995.], into univariate nonlinear functions associated with each filter in the filter-bank used for kernel-expansion [51Seylaz J, Hara H, Pinard E, Mraovitch S, MacKenzie E T, and Edvinsson L, "“Effect of stimulation of the sphenopalatine ganglion on cortical blood flow in the rat,”", J. Cereb. Blood Flow Metab, vol. 8, pp. 875-878, 1988., 58Wiener N, Nonlinear Problems in Random Theory, MIT Press: Cambridge, Mass, 1958.]. This is accomplished by use of Singular Value Decomposition (SVD) of a properly constructed rectangular matrix that includes the significant eigenvectors of the estimated self-kernels, as described in detail below. Since the latter “separability” of the multi-input static nonlinearity of the Wiener-Bose model cannot be guaranteed in all cases, we include in the PDM-based model cross-terms that are properly selected to capture the most significant interactions among the PDM “channels”. The latter is achieved by means of SVD of another matrix that is constructed on the basis of the estimated cross-kernels and yields pairs of cross-PDMs. The structure of the PDM-based Volterra-Wiener model of the two-input/one-output CFA system is shown in Fig. (2), since three PDMs for each input and two pairs of cross-PDMs were found to be adequate in this case. Note that these PDMs are “global” in the sense that they are obtained from the data of all subjects in the reference group and represent a “functional coordinate system” (i.e. basis of functions) for efficient representation of all Volterra kernels of the CFA system, which provide a complete description of the dynamic characteristics of this system. The estimated static nonlinearities that are associated with each PDM channel and the coefficients of the cross-terms are subject-specific and can be used to characterize uniquely the CFA process for each subject.
The proposed methodology commences with the estimation of a 2^{nd} order Volterra model of the two-input CFA system using Laguerre expansions with the following parameters: L1= 5, alpha1= 0.5, L2= 3, alpha2= 0.6, where MABP is designated as input #1 and ETCO2 as input #2. This results in 45 free parameters (including a constant baseline term) for the two-input 2^{nd} order CFA Volterra model, which can be adequately supported (in terms of statistical estimation) by a minimum of 3 min of data. Then, we perform Singular Value Decomposition (SVD) of a rectangular matrix for each input composed of the 1^{st} order kernels of all subjects in the reference group (4 of the 12 subjects, randomly selected) and the three most significant eigenvectors of their 2^{nd} order self-kernels (resulting from eigen-decomposition) weighted by the respective eigenvalues in product with the root-mean-square (RMS) value of the respective input. We select the three most significant “singular vectors” of this matrix (corresponding to the three largest “singular values”) for each input as the global PDMs for the respective input. The pairs of global cross-PDMs are selected by means of SVD of a rectangular matrix that is composed of all the estimated 2^{nd} order cross-kernels of the reference group. The significant singular vectors from the “prior” and “posterior” SVD matrices (termed U and V in the standard SVD terminology) are selected as the pairs of “cross-PDMs” which form the cross-terms as products of the convolutions of these cross-PDMs with the respective inputs. It is important to note that the waveforms of the global PDMs and cross-PDMs were not affected significantly when different sets of 4 subjects (from the full set of 12 subjects) were randomly selected. This fact corroborates the premise of the existence of global PDMs and cross-PDMs for the CFA system and the potential of the PDM-based model to be generalizable in form (i.e. applicable) to all subjects.
In order to complete the PDM-based nonlinear model, we must further estimate the “associated nonlinear function” (ANF) for each global PDM, which is a static nonlinearity applied to the convolution of the input signal with the respective global PDM, as well as the appropriate coefficients of the cross-terms that account for the inter-modulation effects between the two inputs. The model output prediction is composed of the sum of all ANF outputs and cross-terms, along with a constant baseline value. The coefficients of six cubic ANFs (one for each of the six PDMs) and of four multiplicative cross-terms (representing all combinations of the two pairs of cross-PDMs) are estimated for the CFA PDM-based model via least-squares fitting of the output equation, along with the constant baseline value. The estimated coefficients of the ANFs and the cross-terms are distinct for each subject and can be used to quantify uniquely its CFA characteristics, offering a potential diagnostic tool in future clinical practice. In this study, it was found that cubic ANFs and four multiplicative cross-terms are adequate for the CFA system. The total number of free parameters for this two-input PDM-based model is 23 (including the constant baseline term), which implies a minimum required data-record of 2 min. By comparison, the Laguerre-based 2^{nd} order Volterra model (with L1=5, L2=3) has 45 free parameters and, if the self-kernels are extended to 3^{rd} order (to be comparable with the cubic ANFs of the PDM-based model), the number of free parameters rises to 90.
The presented PDM-based modeling methodology can be also applied to the single-input case (MABP-to-MCBFV), both linear and nonlinear, which has been studied extensively in the past because it yields measures of frequency-dependent “vascular admittance” (the inverse of vascular impedance). In this case, the number of free parameters is reduced to 6 for linear modeling (L=5), while for nonlinear modeling this number becomes 21 for the 2^{nd} order Laguerre-based model or 10 for the PDM-based model with 3 PDMs and cubic ANFs. A practical implication of the reduced number of free parameters is that shorter data-records can support the estimation task or, alternatively, greater estimation accuracy is achieved for the same data length. The results of PDM-based linear and nonlinear modeling of the CFA system are presented in the following section and are juxtaposed to results obtained using previously introduced methods, starting with the single-input (MABP-to-MCBFV) case.
We begin with linear single-input modeling of CFA using the Transfer Function measurement between MABP (input) and MCBFV (output) that was employed in previous studies [43Panerai R B, Eames P J, and Potter J F, "“Variability of timedomain indices of dynamic cerebral autoregulation,”", Physiol. Meas, vol. 24, pp. 367-381, 2003., 44Panerai R B, Kerins V, Fan L, Yeoman P M, Hope T, and Evans D H, "“Association between dynamic cerebral autoregulation and mortality in severe head injury,”", Br. J. Neurosurg, vol. 18, no. 5, pp. 471-479, 2004.], for the purpose of comparison with the proposed PDM-based approach. The Transfer Function is estimated as the ratio of the input-output cross-spectrum S_{fp} to the input spectrum S_{fp}, where p(t) denotes the input MABP signal and f(t) denotes the output MCBFV signal:
Note that the computed spectra are statistical estimates because they are based on finite data-records of random signals (i.e. they have nonzero statistical variance). The magnitude of the obtained Transfer Function estimate (Gain Function) represents a frequency-dependent measure of vascular admittance (the inverse of impedance). Fig. (3) shows the Gain Functions for the 12 subjects of this study: average waveform and standard-deviation bounds (left panel) and illustrative Gain Functions from 3 subjects (right panel). These plots demonstrate the high statistical variability of the Gain (or Transfer) Function measurement. This motivates the use of input-output modeling with Laguerre expansions that reduces this variability, as demonstrated below. The inverse Fourier Transform of the Transfer Function is the Impulse Response Function (IRF),h(τ) of the linear model of the MABP-to-MCBFV relationship in accordance with the convolutional relation:
Fig. (4) shows the obtained IRFs for the 12 subjects of this study, indicating considerable inter-subject variability and statistical estimation variance for each subject. This variability can be reduced significantly using Laguerre expansions for the estimation of the 1^{st} order (linear) Volterra model, as demonstrated below. Note that these IRF estimates are scaled via linear regression so that their prediction, based on (2), achieves minimum mean-square error. The IRF estimate, h(τ), can be used in conjunction with (2) to generate, through convolution with the input MABP signal p(t), linear model predictions of the output MCBFV signal f(t). The mean (standard deviation) of the prediction NMSE (Normalized Mean Square Error) of these linear models for the 12 subjects is: 63.94% (17.30%). The considerable error in the model predictions points to the prevailing low signal-to-noise ratio and possible system nonlinearities in the CFA system.
Before we extend our study to nonlinear modeling, we compare the linear model predictions based on the IRF estimates of Fig. (4) with the 1^{st} order kernels of the linear Volterra models using Laguerre expansions with L= 5, alpha= 0.5 [31Marmarelis V Z, Nonlinear Dynamic Modeling of Physiological Systems, Wiley-Interscience & IEEE Press: Hoboken NJ, 2004.], which are shown in Fig. (5) along with their frequency-domain counterparts (magnitude only). The latter are estimates of the “apparent” Gain Functions. The adjective “apparent” is used to point out that these are “linearized” model approximations of an actual nonlinear system. The phase of the Transfer Function may also be used to explore phase relations between the MABP and MCBFV signals (e.g. the phase lead of the flow signal relative to the pressure signal that is associated with the compliance of the vasculature in the Windkessel model). We note that the 1^{st} order Volterra kernel represents the best linear model of a nonlinear system and, therefore, it is not generally the same as the IRF which is derived under the assumption that the system is linear [59Zhang R, Zuckerman J H, Giller C A, and Levine B D, " “Transfer function analysis of dynamic cerebral autoregulation in Humans,”", Am. J. Physiol, vol. 274, pp. H233-H241, 1998.]. It is evident that the Laguerre-based 1^{st} order kernel estimates (and the corresponding apparent Gain Functions) have some differences but exhibit smaller statistical variance than their counterparts (IRFs) obtained via cross-spectral methods. We note that the zero-lag values of these kernels are significantly larger than their IRF counterparts for all subjects. The resulting mean (standard deviation) of the prediction NMSE of these Laguerre-based linear models for the 12 subjects is: 55.47% (15.42%) -- better than their IRF-based counterparts (55% versus 64%) and with smaller variability. We explore below further improvements in prediction accuracy using Laguerre-based nonlinear models.
Nonlinear modeling of CFA has been attempted previously by our group with Volterra models (for both the single-input and the two-input cases) following the method of Laguerre-Volterra networks which combine Laguerre kernel expansions with Volterra-equivalent networks through an iterative estimation procedure [33Mitsis G D, Zhang R, Levine B D, and Marmarelis V Z, " “Modeling of nonlinear physiological systems with fast and slow dynamics II: Application to cerebral autoregulation,”", Ann. Biomed. Eng, vol. 30, pp. 555-565, 2002.-36Mitsis G D, Zhang R, Levine B D, Tzanalaridou E, Katritsis D G, and Marmarelis V Z, "“Nonlinear analysis of autonomic control of cerebral hemodynamics,”", IEEE Eng. Med. Biol. Mag, vol. 28, pp. 54-62, 2009.]. The obtained results in the single-input (MABP) case showed significant prediction improvement for the nonlinear model, corroborating the importance of nonlinear modeling of CFA. Nonetheless, the complexity of the obtained models hinders their physiological interpretation and certain practical limitations persist in the convergence of the iterative approach employed by the Volterra-equivalent networks. These limitations are removed by the proposed PDM-based nonlinear modeling approach, because it is not iterative and reduces the number of free parameters in the nonlinear model. The advocated PDM modeling approach relies on a common set of “global” PDMs to offer a generalizable and (potentially) physiologically interpretable model, as shown below. The global PDMs are obtained from the data of a “reference group”, randomly selected among the available subjects (4 out of 12 available subjects in this case, with the remaining 8 subjects used as the “test set” for validation purposes). We note that the waveforms of the global PDMs did not change appreciably when different subjects were used in the randomly selected reference group.
Following the procedure for the estimation of the global PDMs for the single-input (MABP) case that was outlined in the Methods section, we obtain the three global PDMs shown in Fig. (6) both in the time-domain and the frequency-domain (FFT magnitude). We observe that the 1^{st} PDM exhibits a high-pass characteristic, consistent with the traditional Windkessel model and passive fluid mechanics of blood circulation. This PDM exhibits a blunt spectral peak ~ 0.25 Hz, probably related to the respiratory sinus arrhythmia (breathing cycle). The 2^{nd} PDM exhibits two spectral peaks around 0.03 Hz and 0.15 Hz, while the 3^{rd }PDM exhibits a single resonant peak around 0.05 Hz. The source of these resonant peaks ought to be explored in future studies, but it is likely to be related to myogenic mechanisms and the autonomic nervous system [36Mitsis G D, Zhang R, Levine B D, Tzanalaridou E, Katritsis D G, and Marmarelis V Z, "“Nonlinear analysis of autonomic control of cerebral hemodynamics,”", IEEE Eng. Med. Biol. Mag, vol. 28, pp. 54-62, 2009.]. Our previous studies have shown that most of the nonlinear characteristics of CFA are exhibited in the range 0.02 – 0.08 Hz [33Mitsis G D, Zhang R, Levine B D, and Marmarelis V Z, " “Modeling of nonlinear physiological systems with fast and slow dynamics II: Application to cerebral autoregulation,”", Ann. Biomed. Eng, vol. 30, pp. 555-565, 2002.-36Mitsis G D, Zhang R, Levine B D, Tzanalaridou E, Katritsis D G, and Marmarelis V Z, "“Nonlinear analysis of autonomic control of cerebral hemodynamics,”", IEEE Eng. Med. Biol. Mag, vol. 28, pp. 54-62, 2009.].
To complete the PDM model for the single-input case, we must now estimate the “associated nonlinear function” (ANF) for each global PDM, which is a static (cubic in this case) nonlinearity applied to the convolution of the input MABP signal with the respective global PDM. These ANFs are distinct for each subject and can be used to quantify the CFA characteristics in each patient, offering a potential diagnostic tool. As an illustrative example, Fig. (7) shows the three ANFs of the respective global PDMs for subject #4. It is evident that the ANF of the 1^{st} PDM is almost linear and much larger than the other two ANFs (i.e. its contribution to the model output is larger). We also observe that the positive slope of the 1^{st} PDM is consistent with passive fluid dynamics. However, the small negative slope of the 2^{nd} PDM is not consistent with passive fluid dynamics and suggests the presence of an active mechanism of counter-regulation. The ANF of the 3^{rd} PDM has a positive slope and exhibits a nonlinear characteristic that tends to emphasize the effects of blood pressure increases. The resulting PDM-based nonlinear model predictions have NMSE values with mean (standard deviation) of: 53.8% (14.7%) over all 12 subjects. We note that, although this is not a significant improvement in average prediction accuracy over the linear Laguerre-based models, these nonlinear models are based on the global PDMs of Fig. (6) and, therefore, have far fewer free parameters than their Volterra counterparts which have provided nonlinear model predictions of greater accuracy [48Paulson O B, Strandgaard S, and Edvinsson L, "“Cerebral autoregulation,”", Cerebrovasc. Brain Metab. Rev, vol. 2, pp. 161-192, 1990.].
We proceed with the estimation of the “global” PDMs when two inputs are used (MABP and ETCO2) for the same output (MCBFV). The procedure is the same as for the single-input case and involves the Laguerre-based 2^{nd} order self-kernel estimates for the respective input. The cross-kernels are analyzed separately to yield the cross-PDMs, as described in the Methods section. The two-input 2^{nd}-order Volterra model is obtained using the Laguerre expansion parameters: L1= 5, alpha1= 0.5, L2= 3, aplha2= 0.5. The obtained global PDMs for the MABP and ETCO2 inputs are shown in Fig. (8) in the time and frequency domains. It is encouraging that the global PDMs for the MABP input are almost identical with the single-input case, suggesting that the effects of the two inputs are segregated by the proposed methodology. The 1^{st} PDM of the ETCO2 input has a high-pass characteristic (like its MABP counterpart) and exhibits a blunt spectral peak ~0.15 Hz, possibly related to the autonomic nervous system (see Discussion). The 2^{nd} PDM of the ETCO2 input has a low-pass characteristic (half-max frequency of ~ 0.08 Hz), indicating integration of the input values over 6-8 sec. The 3^{rd }PDM of the ETCO2 input exhibits a single resonant peak around 0.07 Hz. The origin of these mechanisms (related to the action of chemoreceptors and autonomic control or endothelial processes) ought to be explored in future studies (see Discussion).
The obtained ANF coefficients for these global PDMs indicate the following:
Completion of the PDM-based two-input nonlinear model requires the inclusion of cross-terms that are derived from the cross-kernel estimates. Application of the procedure outlined in Methods yields the two cross-PDM pairs shown in Fig. (9) that account for inter-modulation effects between the two inputs (i.e. model terms that are products of the convolutions of the cross-PDM pairs with the respective inputs). The spectral characteristics of these cross-PDMs are similar to the global PDMs. The coefficients of the cross-terms in the PDM model are estimated, along with the coefficients of the cubic ANFs of the global PDMs (and the baseline value) through least-squares fitting to the output data. Note that these cross-PDMs can be expressed as linear combinations of the global PDMs for the respective input. The obtained coefficients for the cross-PDM terms indicate that their significance varies considerably among subjects, both in terms of magnitude and sign. The average (standard deviation) value for the prediction NMSE of the PDM-based two-input model was: 40.4% (14.4%).
In order to examine the functional characteristics of the two-input PDM nonlinear model, we compute the model response to step changes of the MABP input for fixed level of ETCO2. Two illustrative cases are shown in Fig. (10) for subject #1. The left panel shows the model response under hypercapnic conditions (ETCO2 elevated to +1 mmHg relative to baseline) for an initial step increase of MABP (+2 mmHg relative to baseline) for 20 sec, a subsequent drop to baseline for 20 sec, and another drop to -2 mmHg relative to baseline for another 20 sec (note that the settling time of this model is < 20 sec). The right panel shows the model response under hypocapnic conditions (ETCO2 level dropped by 1 mmHg relative to baseline) for the same sequence of MABP step changes. As expected, the PDM model correctly predicts the blood flow autoregulation for step changes in blood pressure (showing a return of MCBFV to baseline after a transient caused by the MABP step change) and an increase in steady-state MCBFV in response to elevated ETCO2 or a decrease in response to reduced ETCO2 (normal CO2 vasoreactivity).
In order to examine the nonlinear characteristics of the two-input CFA model, we compute the steady-state model response to step changes of the two inputs for various values within the physiological range. The resulting steady-state response surfaces for subjects #3 and #8 are shown in Fig. (11) for illustrative purposes. A homeostatic plateau is evident with respect to MABP step changes (demonstrating the effect of autoregulation). It is also evident that the CO2 vasoreactivity is supralinear.
In this study, we seek to advance our quantitative understanding of the process of cerebral flow autoregulation (CFA) by use of dynamic nonlinear models based on the concept of Principal Dynamic Modes (PDMs). CFA is defined as the dynamic process by which cerebral blood flow is maintained within physiologically acceptable bounds during fluctuations in cerebral perfusion pressure secondary to changes in arterial blood pressure. This study also addresses the dynamic modulatory effects of changes in blood CO2 tension on the CFA process. Our findings indicate that the effects of changes in CO2 tension on the pressure-flow relation (CO2 vasoreactivity) ought to be included in the PDM-based model. Note that pressure autoregulation and CO2 vasoreactivity are mediated by different physiological mechanisms.
The PDMs of the CFA model were estimated from beat-to-beat time-series data of two “inputs”: MABP (Mean Arterial Blood Pressure) and ETCO2 (End-Tidal CO2), and one “output”: MCBFV (Mean Cerebral Blood Flow Velocity). The obtained PDM-based models were validated through their predictive capability and are deemed capable of predicting the MCBFV “output” for any given pair of “inputs” MABP and ETCO2 within the dynamic range of the analyzed data – notwithstanding the effects of possible system nonstationarity. The memory of this input-output system was found to be less than 20 sec (i.e. the effects of a change in MABP or ETCO2 on MCBFV are expressed within 20 sec). Although there is considerable inter-subject variability in the CFA process, three “global” PDMs for each of the two inputs were found to constitute an adequate and consistent “functional coordinate system” for representing the dynamics of the CFA process for all 12 subjects in this study. The “Associated Nonlinear Functions” (ANFs) of the PDMs were estimated from the data and used to describe quantitatively the nonlinearities of the CFA process for each subject. The existence of such a set of global PDMs makes this formidable modeling task feasible in practice and facilitates the physiological interpretation of the modeling results. It is posited that the form of the global PDMs contains information of physiological importance regarding the underlying mechanisms of the CFA process. Future physiological studies should seek to examine which specific mechanisms give rise to these global PDM waveforms and their various features. There is a similar need for physiological interpretation of the form of the subject-specific ANFs in future studies, which will be essential for their prospective use in improved clinical diagnosis for a host of important pathologies such as stroke, traumatic brain injury or neurodegenerative diseases. The main findings of this study are summarized below.
Nonlinear modeling/analysis of the CFA process was attempted previously with Volterra kernels and the Laguerre-Volterra network [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983., 38Muizelaar J P, Ward J D, Marmarou A, Newlon P G, and Wachi A, "“Cerebral blood flow and metabolism in severely head-injured children. Part 2: Autoregulation,”", J. Neurosurg, vol. 71, pp. 72-76, 1989., 47Panerai R B, "“Transcranial Doppler for evaluation of cerebral autoregulation,”", Clin. Autonomic Res, vol. 19, pp. 197-211, 2009.-50Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, and Hetzel A, "“Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations,”", Stroke, vol. 36, pp. 1684-1689, 2005.]. It was shown that such nonlinear models reduced significantly the output (MCBFV) prediction error. However, these models retained formidable complexity and exhibited considerable inter-subject variability, both of which hindered their physiological interpretation and their broader acceptance by the research community. The PDM-based nonlinear models are functionally equivalent with the kernel-based Volterra model of the same order (from the output prediction point of view) but their complexity is significantly reduced in terms of model structure and number of free parameters. The PDM modeling approach allows estimation of higher order nonlinearities and facilitates the physiological interpretation of the obtained model. It also makes the model estimation more robust and accurate from short data-records – thus mitigating the effects of intrinsic nonstationarities in the CFA process [21Kontos H A, Wei E P, Raper A J, Rosenblum W I, Navari R M, and Patterson J L Jr, "“Role of tissue hypoxia in local regulation of cerebral microcirculation,”", Am. J. Physiol, vol. 234, pp. H582-H591, 1978., 46Panerai R B, "“Cerebral autoregulation: from models to clinical applications,”", Cardiovasc. Eng, vol. 8, pp. 42-59, 2008., 50Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, and Hetzel A, "“Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations,”", Stroke, vol. 36, pp. 1684-1689, 2005., 60Zhang R, Zuckerman J H, Iwasaki K, Wilson T E, Crandall C G, and Levine B D, "“Autonomic neural control of dynamic cerebral autoregulation in humans,”", Circulation, vol. 106, pp. 1814-1820, 2002.].
The obtained models demonstrate the existence of nonlinearities in the CFA process within the dynamic range of the analyzed data, especially in the frequency range below 0.08 Hz. These nonlinearities are relatively mild for the pressure-dependent components of the model, but somewhat stronger for the CO2-dependent components of the model that allow rapid flow increase in response to hypercapnia (as quantified by the large quadratic ANF coefficient for the 1^{st} PDM of the ETCO2 input). Cross-interaction terms were identified in the form of cross-PDMs that are equivalent to a Volterra cross-kernel. These terms describe the intermodulation of pressure and CO2 effects on cerebral blood flow. The Volterra approach (kernel-based or PDM-based) is the only known method that can quantify such cross-interactions.
The validity of the obtained PDM models was also demonstrated in the case of step changes in the MABP and ETCO2 inputs (see Fig. 10) where the physiologically expected response was predicted by the model – i.e. a return of MCBFV values to baseline after a step increase of MABP (pressure autoregulation) and an increase of MCBFV after a step increase of ETCO2 (CO2 vasoreactivity). We note the big difference between the transient and the steady-state responses that illustrates the clinical importance of understanding the dynamic CFA process (as opposed to static CFA).
The results of this study have confirmed the view that the effects of CO2 on cerebrovascular characteristics are of primary importance with respect to cerebral perfusion and, therefore, a CO2-dependent input should be included in the CFA model along with the pressure-dependent input. It was further shown that the CO2 effects on cerebral blood flow are more nonlinear than the pressure effects and are manifested mostly in the low frequency range from 0.02 to 0.08 Hz. The obtained PDM-based models show high inter-subject variability in both the CO2-dependent and the pressure-dependent component of the model.
The specific physiological mechanisms that determine the form of the obtained global PDMs and the subject-specific ANFs require future studies. We provide herein an initial interpretation of the observed spectral characteristics of the global PDMs. Specifically, the 1^{st} PDM of MABP exhibits a high-pass characteristic (half-max frequency of ~0.05 Hz), consistent with the traditional Windkessel model of fluid mechanics in the passive vasculature. This PDM has a blunt spectral peak ~ 0.25 Hz, probably related to the breathing cycle. Its ANF has a strong linear coefficient with consistent positive sign, suggestive of its fluid mechanical origin. The 2^{nd} PDM of MABP exhibits two spectral peaks around 0.03 Hz and 0.15 Hz, and the 3^{rd }PDM of MABP exhibits a resonant peak around 0.05 Hz. These resonant peaks may be related to the autonomic nervous system [50Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, and Hetzel A, "“Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations,”", Stroke, vol. 36, pp. 1684-1689, 2005., 61Zhang R, Witkowski S, Fu Q, Claassen J A, and Levine B D, " “Cerebral hemodynamics after short- and long-term reduction in blood pressure in mild and moderate hypertension,”", Hypertension, vol. 49, pp. 1149-1155, 2007.]. Although this question deserves careful examination in future studies, some initial indications in support of these hypotheses are found in our previous studies of Volterra modeling under conditions of ganglionic blockade of autonomic neural activity using trimethaphan and sympathetic activation during lower body negative pressure [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983., 50Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, and Hetzel A, "“Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations,”", Stroke, vol. 36, pp. 1684-1689, 2005.].
The interpretation of the PDM waveforms is even more challenging for the ETCO2 input, because less information is available about the specific dynamics of this process. Fig. (8) shows that the 1^{st} PDM of ETCO2 exhibits a high-pass characteristic with half-max frequency of about 0.05 Hz, akin to its MABP counterpart, but the blunt spectral peak is now around 0.15 Hz. The ANF of this 1^{st} PDM of ETCO2 has a large quadratic coefficient, indicating stronger response to hypercapnia than hypocapnia. The 2^{nd} PDM of ETCO2 exhibits a low-pass characteristic with bandwidth (half-max frequency) of about 0.08 Hz. This PDM integrates the ETCO2 input over ~8 sec and, therefore, it should be critical in cases of temporary hypercapnia or hypocapnia during hypoventilation or hyperventilation, respectively – a hypothesis supported by the fact that the sign of the first-degree (linear) coefficient of this 2^{nd} PDM is consistently positive. Note that our previous trimethaphan study [17Iadecola C, Nakai M, Mraovitch S, Ruggiero D A, Tucker L W, and Reis D J, "“Global increase in cerebral metabolism and blood flow produced by focal electrical stimulation of dorsal medullary reticular formation in rat,”", Brain Res, vol. 272, pp. 101-114, 1983.] showed suppression of the low-frequency gain for ETCO2, suggesting that this 2^{nd} PDM is influenced by autonomic neural activity. The importance of the 2^{nd} PDM of ETCO2 is underlined by the fact that most of the power of the ETCO2 signal resides within the high-gain low-frequency range of this PDM (see Figs. 1 and 8). The 3^{rd} PDM of ETCO2 exhibits a spectral peak around 0.08 Hz and the underlying physiological mechanism must be explored in future studies.
This point was made earlier but its importance invites further elaboration. The use of “global” PDMs in the advocated methodology implies that the CFA system dynamics can be properly constrained (in first approximation) for all subjects by the “functional coordinate system” defined by these global PDMs. This represents a “generalizability” of the CFA model that is potentially useful, because it reduces the model complexity and facilitates its physiological interpretation with potential clinical implications. Open questions remain with regard to the subject-specific form of the ANFs, where a cubic polynomial form has been initially imposed in order to explore this approach. This issue deserves more attention in future studies.
In addition to exploring the most appropriate form for the ANFs, some attention should be directed towards the use of PDM-based “linearized” models (i.e. the “best” linear approximations of the CFA system achieved in the PDM framework), since they may represent an efficient way to achieve improved clinical diagnosis without the “overhead” of a detailed PDM-based nonlinear model. Note that the “linearized” model is generally distinct from the linear component of a nonlinear Volterra or PDM model, since the high-order kernels generally contribute to the 1^{st} order kernel estimate of the “linearized” model [51Seylaz J, Hara H, Pinard E, Mraovitch S, MacKenzie E T, and Edvinsson L, "“Effect of stimulation of the sphenopalatine ganglion on cortical blood flow in the rat,”", J. Cereb. Blood Flow Metab, vol. 8, pp. 875-878, 1988.]. Specifically, we envision two CFA-related indices of potential clinical utility that are derived from these PDM-based “linearized” models: (1) the Effective Pressure Reactivity (EPRx) index that is defined as the steady-state change of MCBFV in response to a unit-step change of MABP (i.e. the integral of the 1^{st} order MABP-related kernel of the linearized CFA model); (2) the Effective CO2 Reactivity (ECRx) index that is defined as the steady-state change of MCBFV in response to a unit-step change of ETCO2 (i.e. the integral of the 1^{st} order ETCO2-related kernel of the linearized CFA model). These indices are subject-specific and can be easily computed in practice.
The PDM-based model can be related to commonly used clinical indices of CFA for clinical studies of traumatic brain injury, hypertension, diabetes and stroke. For instance: (a) the cerebrovascular resistance (CVR) index is obtained as the ratio of short-time averages of MCBFV over MABP or through computation of transfer functions [13Giller C A, Bowman B, Dyer H, Mootz L, and Krippner W, " “Cerebral arterial diameters during changes in blood pressure and carbon dioxide during craniotomy,”", Neurosurgery, vol. 32, pp. 737-741, 1993., 40Newell D W, Aaslid R, Lam A, Mayberg T S, and Winn H R, " “Comparison of flow and velocity during dynamic autoregulation testing in humans”", Stroke, vol. 25, pp. 793-797, 1994.]; (b) the correlation index (Mx) is computed from short-term correlations (averaged over 5-10 sec) between beat-to-beat time-series data of pressure and flow velocity [45Panerai R B, Moody M, Eames P J, and Potter J F, "“Cerebral blood flow velocity during mental activation: interpretation with different models of the passive pressure-velocity relationship,”", J. Appl. Physiol, vol. 99, pp. 2352-2362, 2005.]; (c) the “autoregulation index” (ARI) is computed from the integrated impulse response of an estimated second-order linear model [62Zhang R, Behbehani K, and Levine B D, "“Dynamic pressure–flow relationship of the cerebral circulation during acute increase in arterial pressure,”", J. Physiol, vol. 587, no. Pt 11, pp. 2567-2577, 2009.]; (d) the “pressure reactivity” index (PRx) is computed as the cumulative flow response to a step pressure change; (e) the “CO2 reactivity” index (CRx) is computed as the cumulative flow response to a step CO2 change CO2 [19Kalaria R N, "“Cerebral vessels in ageing and Alzheimer’s Disease,”", Pharmacol. Ther, vol. 72, pp. 193-214, 1996.]; and (f) the phase difference between oscillations in MABP and MCBFV has been used as another autoregulation index [42Panerai R B, Simpson D M, Deverson S T, Mahony P, Hayes P, and Evans D H, "“Multivariate dynamic analysis of cerebral blood flow regulation in humans,”", IEEE Trans. Biomed. Eng, vol. 47, pp. 419-423, 2000.]. These indices can be related to the PDM-based model as follows: (a) the CVR index via the Gain Function estimate (i.e. Fourier Transform magnitude of 1^{st} order MABP kernel); (b) the Mx index via averaging of the 1^{st} order MABP kernel over lags from 1 to 10 sec (because the input-output cross-correlation for a linearized model is equal to the convolution of the 1^{st} order kernel with the autocorrelation of the input ); (c) the ARI via the integral of the 1^{st} order MABP kernel of the “linearized” model fitted to a second-order differential equation (the Tiecks model); (d) the PRx index via the steady-state response to a step MABP input; (e) the CRx index via the steady-state response to a step ETCO2 input; and (e) the phase-difference index through the phase function of the Fourier Transform of the 1^{st} order MABP kernel of the “linearized” model.
Thus, the PDM-based model may unify all these clinical indices. Most importantly, the two-input PDM-based model of CFA provides critical clarity with respect to the contemporaneous effects of CO2 tension and arterial pressure on CFA. This allows computation of pressure-dependent indices from the PDM-based model for isocapnic conditions (steady CO2 tension), whereas this is not possible for previously defined indices – an important distinction for clinical practice, since the effect of CO2 on CFA is known to be significant. Likewise, the CO2-dependent indices from the PDM model (e.g. the PDM-based ECRx index) can be defined for isobaric conditions (steady arterial pressure), whereas this is not possible for other indices of CO2 vasomotor reactivity. Finally, we note that the PDM-based indices can be computed from spontaneous activity data and do not require specialized maneuvers (e.g. Valsalva) that may drive the CFA system outside its physiological range.
The potential importance of the presented modeling methods for advancing our understanding of the CFA process and its implications for medical science necessitate further studies with a larger population in order to confirm these initial results and explore the specific physiological mechanisms underlying the characteristics of the PDM-based model.
The general nonparametric Volterra model of a finite-memory dynamic nonlinear system with a single input x(t) and a single output y(t) is given by the functional series expansion [52Steiner L A, Czosnyka M, Piechnik S K, Smielewski P, Chatfield D, Menon D K, and Pickard J D, "“Continuous monitoring of cerebrovascular pressure reactivity allows determination of optimal cerebral perfusion pressure in patients with traumatic brain injury,”", Crit. Care Med, vol. 30, pp. 733-738, 2002.]:
In this context, the modeling task involves the estimation of the unknown Volterra kernels of the system {k_{n}} from given input-output data x(t) and y(t). It has been shown [27Marmarelis V Z, "“Identification of nonlinear biological systems using Laguerre expansions of kernels,”", Ann. Biomed. Eng, vol. 21, pp. 573-589, 1993.] that this task is facilitated immensely by Laguerre expansions of the kernels:
where {b_{j} (τ)} denote the orthonormal Laguerre function basis. Such kernel expansion yields the following nonlinear input-output relation which involves linearly the Laguerre expansion coefficients {c_{r}}:The fact that the Laguerre expansion coefficients enter linearly in the nonlinear Volterra model allows their estimation via least-squares regression (a relatively simple, robust and stable numerical procedure). Having estimated the Laguerre expansion coefficients, we can construct the Volterra kernel estimates using Equation (A2) and compute the model prediction for any given input using Equation (A1) or (A3). To date, Volterra models have been usually estimated up to 2^{nd} order.
The introduction of the concept of Principal Dynamic Modes (PDMs) has allowed the practical estimation of nonlinear models of higher order [31Marmarelis V Z, Nonlinear Dynamic Modeling of Physiological Systems, Wiley-Interscience & IEEE Press: Hoboken NJ, 2004.] as in the subject application. Briefly stated, the use of PDMs allows us to write the output equation as:
where {u_{p} (t)} are the PDM outputs (i.e. convolutions of the input with the respective PDM) and f_{p}[u_{p}]} are the static nonlinearities associated with each PDM (ANFs) which are typically given polynomial form (cubic in this application). This general modeling methodology has been extended to systems with multiple inputs and multiple outputs [31Marmarelis V Z, Nonlinear Dynamic Modeling of Physiological Systems, Wiley-Interscience & IEEE Press: Hoboken NJ, 2004.].