A cerebrovascular response model for functional neuroimaging including dynamic cerebral autoregulation

doi:10.1016/j.mbs.2009.05.002

Math Biosci. Author manuscript; available in PMC 2010 August 1.

Published in final edited form as:

Math Biosci. 2009 August; 220(2): 102–117.

Published online 2009 May 13. doi: 10.1016/j.mbs.2009.05.002

PMCID: PMC2720139

NIHMSID: NIHMS117087

A cerebrovascular response model for functional neuroimaging including dynamic cerebral autoregulation

Solomon Gilbert Diamond,¹^* Katherine L. Perdue,¹ and David A. Boas²

¹ Thayer School of Engineering at Dartmouth, 8000 Cummings Hall, Hanover, NH, 03755, USA

² Massachusetts General Hospital, Athinoula A. Martinos Center for Biomedical Imaging, 13th Street, Building 149, Charlestown, MA, 02129, USA

^* Corresponding author. e-mail address: Solomon.G.Diamond/at/Dartmouth.edu

The publisher's final edited version of this article is available at Math Biosci

See other articles in PMC that cite the published article.

Publisher's Disclaimer

Abstract

Functional neuroimaging techniques such as functional magnetic resonance imaging (fMRI) and near-infrared spectroscopy (NIRS) can be used to isolate an evoked response to a stimulus from significant background physiological fluctuations. Data analysis approaches typically use averaging or linear regression to remove this physiological baseline with varying degrees of success. Biophysical model-based analysis of the functional hemodynamic response has also been advanced previously with the Balloon and Windkessel models. In the present work, a biophysical model of systemic and cerebral circulation and gas exchange is applied to resting state NIRS neuroimaging data from 10 human subjects. The model further includes dynamic cerebral autoregulation, which modulates the cerebral arteriole compliance to control cerebral blood flow. This biophysical model allows for prediction, from noninvasive blood pressure measurements, of the background hemodynamic fluctuations in the systemic and cerebral circulations. Significantly higher correlations with the NIRS data were found using the biophysical model predictions compared to blood pressure regression and compared to transfer function analysis (multifactor ANOVA, p<0.0001). This finding supports the further development and use of biophysical models for removing baseline activity in functional neuroimaging analysis. Future extensions of this work could model changes in cerebrovascular physiology that occur during development, aging and disease.

Keywords: Physiological modeling, dynamic cerebral autoregulation, NIRS, fMRI, Balloon model, Windkessel

Introduction

Models of blood circulation and gas exchange in the human body have advanced steadily over the past decade. System models in the literature span pulmonary [1], systemic [2; 3] and cerebral circulations [4; 5; 6] with inclusion of a wide range of regulatory physiology. The application focus for these models has been on predicting the dynamics of medically relevant events like carotid artery compression, hypotension, hypercapnia, and thigh cuff deflation. Significant progress has also been made in the development of a cerebral autoregulation index [7] that can be used to grade autoregulatory function from noninvasively measured blood pressure and cerebral circulation [8] with applications in, for example, hypertension, stroke, and carotid sinus syndrome [9]. The aim of the present work is to recast a nonlinear dynamic physiological model in the context of neuroimaging. The present approach differs from prior work in that the model dynamics are propagated through the relevant measurement biophysics for direct comparison with functional neuroimaging data. This approach will support the recognized need to evaluate model applicability in the clinical environment and to validate the reproducibility and normal ranges of model parameters [10].

Previous applications of physiological models in neuroimaging have largely focused on describing the functional hemodynamic response to brain activity [11; 12; 13; 14]. A great deal of this interest stems from the Balloon [15] and Windkessel models [16], which have helped to formulate physiologically-based interpretations of fMRI and optical brain imaging data. A limitation of these previous models is that they do not include the effects of cerebral autoregulation, which affects cerebral blood flow. This is important in neuroimaging because the resting state blood flow can confound the interpretation of blood oxygenation level-dependent (BOLD) fMRI [17]. More recent studies offer sophisticated models for the relationship of cerebral blood flow to metabolic biochemistry and vascular smooth muscle [18] and the effects of cerebral autoregulation on neural activation [19]. However, these recent studies used only qualitative comparisons with measurements from the human brain.

The present work focuses on capturing the temporal dynamics of spontaneous fluctuations in cerebral hemodynamics. These physiological fluctuations are apparent in near-infrared spectroscopy (NIRS) [20] and in fMRI [21], where they are usually considered signal clutter or noise. Previous studies have aimed to remove this signal clutter in fMRI with, for example, principal component analysis (PCA) and independent component analysis (ICA) [22], autoregressive models [21], adaptive filters [23], and with Kalman filter in NIRS [24]. Prior neuroimaging studies have not attempted to remove these fluctuations from the data with a physiological model.

An additional motivation to this work, besides reducing signal clutter in brain activation studies, is that model-based assessment of baseline cerebral physiology could serve as biomarkers of certain cerebrovascular pathologies such as impaired dynamic cerebral autoregulation (CA). A current technology for CA assessment is transcranial Doppler ultrasonography (TCD), which can measure cerebral blood flow velocity in certain large intracranial vessels but cannot perform imaging [25]. CA assessment with TCD has been explored in hypertension [26], severe carotid stenosis [27], vertebral artery disease [28], critically ill infants [29] and acute ischemic stroke [30]. Modeling the CA response in fMRI and NIRS data is an exciting new frontier. Some of the parameters in the present model may also be altered during the progression of neurodegerative diseases. For example, changes in vascular stiffness in patients with Alzheimer’s disease or vascular dementia [31; 32] could be represented in the model as a decreased compliance parameter. Another example is diabetes mellitus, which is known to directly affect cerebral vasculature and to cause cerebral blood flow (CBF) abnormalities and complications of the cerebral circulation [33; 34]. The detail included in the proposed model would allow for testing of specific hypotheses about the changes in physiology that accompany disease.

The present study builds on a previously published model of cerebral circulation [2] and a model of systemic circulation [4]. Features have been added that account for the affinity of CO₂ with hemoglobin, nonlinear resistances, Windkessel compliances, and initial parameter values that are derived from steady state relationships. The model is composed of seven systemic compartments, five cerebral compartments, and two neck compartments. These numbers of compartments were chosen to enable the simulated NIRS data to reflect the heterogeneous vascular composition of the scalp and the complex vascular structure of the brain where there are large arteries and veins on the surface, and then arteriolar and venous branches diving down into capillary beds in the cortex. As the modeled blood circulates through each chamber, the hemodynamics are successively filtered into the dynamics that comprise the NIRS data. The 48 system equations for the complete model are derived in this manuscript to make explicit all the necessary information that is needed to reproduce this work. This model will be referred to as SimCVR for Simulated Cerebro-Vascular Response.

Additional features were added to SimCVR that estimate NIRS time-series data from the circulation dynamics so that the model output can be directly compared with measurements from NIRS such as the CW5 optical imager [20] used in the present study. After describing the SimCVR model, it is used to predict the resting state fluctuations in NIRS data at 830 and 690 nm from 10 scalp locations on 10 different subjects. NIRS predictions were generated by driving the SimCVR model with finger-cuff pressure variations and then estimating the blood volume fraction under each NIRS probe. The SimCVR predictions had higher coefficient of determination R² values compared to direct blood pressure regression and compared to transfer function analysis. These results suggest that physiological models may help account for baseline physiology in functional brain imaging experiments and ultimately lead to physiological model-based data interpretation.

Methods

SimCVR consists of coupled fluid flow and gas exchange systems. The fluid flow system is a circuit-analog of nonlinear resistors and capacitors (Fig. 1A) that provide for resistance to fluid flow and compliance of the blood vessels. The circuit can be divided into systemic s and cerebral b regions. The systemic region accounts for blood circulation in the body, excluding the cerebral and pulmonary circulations. There are seven lumped vascular chambers in the systemic region representing the large arteries la, small arteries sa, arterioles a, capillaries c, venules v, small veins sv, and large veins lv. The cerebral circulation is connected through the neck arteries na and neck veins nv. The cerebral model contains five vascular chambers representing pial arteries pa, arterioles a, capillaries c, venules v, and pial veins pv. The cerebral region also accounts for cerebrospinal (CSF) circulation f and intracranial ic pressure. The gas exchange system (Fig. 1B) accounts for the partial pressures of O₂ and CO₂ blood gases g in all of the vascular chambers and, additionally, chambers representing CSF, interstitial fluid (ISF) and intracellular fluid (ICF). The gas transport accounts for convection with the blood flow between vascular chambers and diffusion from the capillaries to the metabolic action in ICF.

Fig. 1

SimCVR model for (A) blood circulation and (B) gas exchange. The nomenclature is defined in Tables 1 and and22.

One of the significant challenges with lumped parameter modeling in physiology is parameter estimation. This issue has been addressed in the circulation model by requiring that the user specify only a subset of the parameters such as the steady-state volumes and pressures. Remaining parameters such as resistances and compliances are calculated based on steady-state relationships. Similar parameter constraints are exploited in the gas exchange system where the user must specify the steady-state blood gas concentrations in the fluid chambers but the diffusion constants are calculated. Exploiting the steady-state parameter relationships makes the task of lumped parameter estimation easier because established physiological quantities such as the steady-state pressure drops through the circulatory system are specified and abstract lumped parameters like effective lumped compliance are calculated. The parameter calculation methodology is detailed in the Appendix section: Parameter relationships at steady state.

In the following model description, the governing state equations are presented in a format that can be implemented directly in an ordinary differential equation (ODE) solver such as the ODE113 function in Matlab® (MathWorks, Inc., Natick, MA).

Systemic fluid flow system

The circulatory model assumes that the pressure-volume relationship in each chamber is governed by compliance and that the pressure drop between chambers and the vascular resistance determines the blood flow. Following the nomenclature defined in Table 1, the blood volume V, vascular resistance R, and blood flow q in vasculature chambers j = la, sa, sv, lv and k = sa, a, lv, ra, are given by

Table 1

Subscript nomenclature

(1)

(2)

and

(3)

Eq. (1) specifies the compliance constant C_s_, _j_,0 in the systemic s, vascular chamber j, as the ratio of volume V_s_, _j to transmural pressure P_s_, _j, which is represented in this electrical analog as a capacitor (Fig. 1A). The unstressed volume V_s_, _j_,_us is the residual volume in the vessel that is present when the transmural pressure falls to zero. The unstressed volume is assumed to be constant and has been initialized as a fraction F_s_,_us of the initial volume

(4)

Following [35], the compliance of the collapsible veins is increased 20 fold when the transmural pressure P_s_, _j is negative. Eq. (2) follows from Poiseuille’s law for laminar flow in a tube in which resistance R_s_, _j is inversely proportional to vessel diameter to the forth power. Eq. (3) states that blood flow q_s_, _j out of chamber j equals the pressure drop from the current chamber j to the subsequent chamber k divided by the resistance R_s_, _j.

Calculating the system response with an ODE solver requires state equations that relate the time derivatives of the pressure equation M5

and volume equation M6

to the pressure P_s_, _j and volume V_s_, _j state variables. The state equations for i = la, v, sv and j = sa, sv, lv are obtained by solving Eq. (1) for P_s_, _j and then taking the time derivative

(5)

where the change in volume is the difference between blood outflow q_s_, _j and inflow q_s_,_i

(6)

The systemic arteriole, capillary and venule compartments include Windkessel compliance to appropriately model nonlinear vascular properties [16]. Using the notational substitutions i = sa, a, c; j = a, c, v; and k = c, v, sv, the compliance is

(7)

where β_s_, _j is the vascular reserve of the jth compartment in the systemic s region and

(8)

is the pressure-volume relationship with time-varying compliance C_s_, _j. The equations for resistance and flow are the same as Eqs. (2) and (3) respectively. The state equations for the Windkessel compartments are obtained from the time-derivative of Eq. (8)

(9)

where

(10)

Solving Eq. (9) for the derivative of pressure yields

(11)

Specifying inputs to the systemic system

The state equation for chamber la is driven by the derivative of the blood pressure in the large artery

(12)

because the input is directly connected to the large artery without resistance in the circuit model (Fig. 1A). The time-varying pressure input P_s_,_la_,_input must be specified by the user. An alternative to using the large artery pressure input of Eq. (12) is to supply the derivative of blood pressure in the small artery sa instead of the large artery la. The advantage of specifying equation M15

is that it can be used to drive the model with the pressure variations measured in the peripheral circulation (e.g. with a finger cuff). However, the following additional calculations are required to compute the unknown value of P_s_,_la from a known value of P_s_,_sa. First Eqs. (1), (2) and (3) are combined with j = la and k = sa to obtain

(13)

Next Eqs. (5) and (6) are combined with j = sa, i = la to obtain

(14)

Eqs. (13) and (14) are then combined and simplified into a cubic equation

(15)

where unknown value of P_s_,_la is only real positive root of Eq. (15).

Cerebral fluid flow system

The cerebral circulation model connects to the systemic circulation model through vascular compartments of the neck. The neck artery na is assumed to be directly coupled to the systemic large artery la such that their pressures are matched

(16)

The equations for volume, resistance and flow for the arteries and veins of the neck have the same form as the Eqs. (1), (2) and (3) with the exchange of b subscripts for s and the substitutions j = na, nv; and k = pa, lv. The state equation for the neck vein nv is

(17)

and is derived in the same manner as Eqs. (5) and (6) where q_b_, _f_₂ is CSF absorption into the neck vein (see Fig. 1A). The equations for resistance R_b_, _j and flow q_b_, _j are the same as Eqs. (2) and (3) with s = b and k = a, nv.

Above the neck, the description of the cerebral compartments can be grouped into the cerebral capillary and venule, which are assumed to have Windkessel compliance, and the pial arteries and pial veins, which are assumed to have constant compliance. This distinction is made because prior studies of the Windkessel model have focused on the microvasculature [36; 37] and not on larger vessels. The cerebral arteriole compliance is treated uniquely due to its involvement in dynamic cerebral autoregulation.

Accounting for the intracranial ic pressure and linear compliance C_b_, _j_,0 of the pial arteries pa and pial veins pv, the volume for j = pa, pv is

(18)

The state equations for the pial vessels are obtained from Eq. (5) with s = b; i = na, v; and j = pa, pv. The cerebral capillary c and venule v compartments include Windkessel compliance C_b_, _j and account for the intracranial pressure

(19)

where j = c, v and

(20)

The state equation for P_b_, _j, derived from Eqs. (19) and (20), further assumes that the intracranial pressure P_b_,_ic may be time varying

(21)

where

(22)

and

(23)

Combining and solving Eqs. (21) and (23) for the derivative of pressure results in

(24)

The expressions for resistance R_b_, _j, flow q_b_, _j and the volume derivatives equation M28

are obtained from Eqs. (2), (3) and (6) with s = b; i = a, c; j = c, v; and k = v, pv. An expression for the time derivative of intracranial pressure equation M29

will be introduced after describing the rest of the cerebral vascular and fluid compartments.

The remaining vascular compartment to be described is the cerebral arteriole. In our model, dynamic cerebral autoregulation of CBF operates via feedback to the cerebral arteriole compliance. This autoregulation-compliance model will first be described and then incorporated into the arteriole state equations.

The CBF set point q_b_,_sp and autoregulation gain K_aut on this feedback loop are influenced by the partial pressure of CO₂ in the ISF. The following system equations that govern the CO₂ effects are simplified versions of those presented in [4] because data simulations were found to operate in a linear regime of the more complex equations. The difference from baseline ISF CO₂ is expressed as an autoregulation set point fraction

(25)

which correlates positively with the CBF set point

(26)

and negatively correlates with the autoregulation reactivity gain

(27)

The autoregulation state variable x_b_,_aut is an auxiliary variable that will determine the arteriole compliance changes a few steps later. The dynamics of x_b_,_aut are governed by its state equation

(28)

The variable cerebral arteriole compliance operates on a non-symmetric sigmoid with the bounds C_b_,_w and curvatures k_b_, _p

(29a)

(29b)

The arteriole compliance is then described by the sigmoid curve

(30)

where exp(·) is the exponential function. Thus if from steady state, the cerebral arteriole flow q_b_,_a falls below the CBF set point q_b_,_sp, the autoregulation state x_b_,_aut will become negative. This will cause the arteriole compliance C_b_,_a to increase towards C_b_,_a_,0 +C_b_,_w_₁, which will dilate the arterioles and increase CBF. In the opposite scenario, C_b_,_a decreases towards C_b_,_a_,0 − C_b_,_w_₁, which constricts the arterioles and decreases CBF.

The arteriole compliance can now be incorporated into the arteriole state equations. The equations for cerebral arteriole volume V_b_,_a, resistance R_b_,_a and flow q_b_,_a are the same as Eqs. (1), (2) and (3) with s = b; j = a; and k = c. The arteriole state equation is obtained from solving for P_b_,_a in Eq. (19) with j = a, and then taking the time derivative to obtain

(31)

where q_b_, _f_₁ is CSF formation from the arteriole (see Fig. 1A), and the compliance derivative is

(32)

The compliance input equation M39

has been added to Eq. (32) to permit inclusion of an arteriolar compliance response to neural stimulus [14], which would be superimposed on the regulatory physiology. When this input is non-zero, the arteriole compliance C_b_,_a must be treated as a state variable that is solved by integration and Eq. (30) would not be used.

In addition to the time-varying fluid flow in the cerebral vascular compartments, the CSF volume in the intracranial space also varies in time because of continuous formation by the choroid plexus and absorption by the arachnoid villi [38]. Following [4], the formation rate depends on the pressure gradient

(33)

and the absorption rate is unidirectional

(34a)

(34b)

The state equation for CSF accounts for the difference between the rates of formation and absorption

(35)

Now that all of the time-varying intracranial fluid volumes in the model are described, the intracranial pressure P_b_,_ic can be expressed in terms of the intracranial compliance C_b_,_ic and volume V_b_,_ic as

(36)

where the intracranial volume is the sum of all intracranial volumes

(37)

and the interstitial fluid V_b_,_ISF, intracellular fluid V_b_,_ICF and tissue volume V_b_,_tissue in the intracranial cavity are assumed to be known constants. Given the definition of V_b_,_ic in Eq. (37), the graphical representation of intracranial compliance C_b_,_ic as a single compartment in Fig. 1A is only symbolic. The state equation for intracranial volume V_b_,_ic is

(38)

where the inflow through the neck artery q_b_,_na, outflow through the neck vein q_b_,_nv, and CSF absorption q_b_, _f_₂ are the only net effects on the volumes that comprise the intracranial volume V_b_,_ic in Eq. (37).

The remaining unknown is the intracranial compliance C_b_,_ic. In [4], the intracranial compliance is defined as inversely proportional to the intracranial pressure. Substituting this inverse proportionality into Eq. (36), however, results in a constant intracranial volume, which implies zero intracranial compliance. Given this inconsistency, we opted instead to use the Windkessel compliance of Eq. (7) with s = b and j = ic. Rearranging Eq. (7) illustrates that larger values of β_b_,_ic result in a compliance expression that is nearly equivalent to the one in [4] but the algebraic cancelation of P_b_,_ic is avoided

(39)

The intracranial pressure state equations that follow from Eqs. (36) and (39) are the same as Eqs. (10) and (11) with s = b and j = ic.

A value of β_b_,_ic can be chosen based Monroe-Kellie doctrine [39; 40], which states that the cranium is rigid and any volume change in one of three constituents (blood, CSF, and brain tissue) must be compensated by change in another constituent, which implies that β_b_,_ic is infinitely large. Monroe-Kellie also permits, however, the ability of the CSF to be displaced into the spinal canal, and some compliance in the falx cerebri between the hemispheres and the tentorium between the hemispheres and the cerebellum, which implies that β_b_,_ic is not infinite. We have, for the moment, set β_b_,_ic =5 based on the empirical observation that this value predicts NIRS data reasonably well and we leave further exploration of this parameter to future studies that include intracranial pressure recordings.

Systemic gas exchange

The systemic gas exchange model tracks the bound and unbound O₂ and CO₂ in the vascular chambers and tissue compartments and accounts for concentration changes that occur from fluid flow, diffusion and cellular metabolism. First, some general expressions for the calculation of blood gas concentration are explained and then applied to the specific chambers in the systemic system. Next, the state equations are obtained for the partial pressure of each gas species and these gas exchange state equations are coupled to the circulatory state equations that were already presented.

In general, the total oxygen concentration C_O_₂ in blood is the sum of bound and unbound O₂

(40)

The bound O₂ concentration depends on the percent oxygen saturation S_O_₂, the hemoglobin concentration M_HbT_,_blood, gas constant R, standard temperature T_STP and pressure P_STP

(41)

The percent oxygen saturation S_O_₂ can be calculated from the partial pressure of oxygen P_O_₂ using an analytic hemoglobin dissociation curve which depends on temperature T and pH [41] (included for convenience in Appendix Blood oxygen). The unbound O₂ concentration simply depends on the partial pressure of oxygen P_O_₂ and the gas solubility α_O_₂

(42)

In the derivation of state equations that follows, it will be necessary to convert oxygen concentration C_O_₂ state equations into oxygen partial pressure P_O_₂ state equations. This conversion can be accomplished using the relation

(43)

where

(44)

is derived from Eqs. (40), (41) and (42) and is referred to as an effective solubility _O_₂ since it accounts for both bound and unbound oxygen [2]

(45)

The general expression for total CO₂ concentration in blood is similar to the expression for oxygen in that both bound and unbound gas must be considered

(46)

The bound and unbound CO₂ concentrations depend on the partial pressure of carbon dioxide P_CO_₂, percent oxygen saturation S_O_₂, hematocrit Hct, gas solubility α_CO_₂, temperature T and pH [42] (see Blood carbon dioxide in Appendix). Also similar to the treatment of oxygen is the need for an effective solubility

(47)

The state equations for the systemic gas exchange can now be described beginning with the input to the partial pressure of gaseous species g = O₂, CO₂ to the systemic arterioles P_s_,_g,a, which is specified as an input in the form of Eq. (12). The blood gas inputs effectively occur at the arteriole chamber because gas exchange is assumed not to occur in any of the arterial chambers. The user may also supply inputs that change the systemic metabolic rates Q_s_,_g in the form of Eq. (12) to capture, for example, the effects of a stimulant drug. Moving from the arteriole a chamber to the capillary c chamber will require accounting for the systemic s blood gas concentration C_s_,_g_,_c for g = O₂, CO₂

(48)

where V_s_,_g_,_c is the gas volume and V_s_,_c is the blood volume. Differentiating Eq. (48) yields

(49)

where the derivative of gas volume V_s_,_g_,_c depends on the gas inflow from the arteriole, outflow from the capillary and diffusion with ISF

(50)

where D_s_,_L_,_g is the diffusion constant, P_s_,_g,c is the partial pressure in the capillary and P_s_,_g_,_ISF is the partial pressure in the ISF. The state equation in the capillary is now obtained by combining Eqs. (48), (49) and (50) and dividing by the effective solubility [Eq. (45) or (47) as appropriate] to arrive at partial pressure in terms of gas species g

(51)

A similar state equation governs the venule v chamber

(52)

The state equation in the intracellular chamber includes the effects of metabolism Q_s_,_g and uses the solubility for unbound gas α_g only because no hemoglobin is present

(53)

A similar expression is use for the interstitial compartment

(54)

Cerebral gas exchange

Following the same procedure as with the systemic compartments, preliminary calculations are carried out to obtain the cerebral blood gas concentration C_b_, _g_, _j and effective solubility _b_, _g_, _j. User inputs of the cerebral metabolic rates can model changing metabolic demands

(55)

The state equations for the cerebral gas exchange system follow a derivation similar that for Eq. (51). Where needed, additional terms account for gas transport with CSF formation q_b_, _f_₁ and absorption q_b_, _f_₂. The following state equations result for the arteriole

(56)

capillary

(57)

venule

(58)

cerebrospinal fluid

(59)

interstitial

(60)

and intracellular compartments

(61)

The lumped parameter model of circulation and gas exchange is now fully described. Additional steps are required, however, to relate the state variables in the lumped model such as blood volume and oxygen saturation in individual vascular chambers to measurements that are made with NIRS from spatially distributed physical tissue.

Hemoglobin concentration changes

The purpose of this section is to relate the state variables in the lumped parameter model to spatially distributed concentration changes in oxy- and deoxyhemoglobin in tissue. Establishing this relationship requires knowledge about the local tissue such as the blood volume fraction, arterial blood volume fraction and so forth. Given the composition of vascular chambers in the tissue and the temporal dynamics in oxy- and deoxyhemoglobin in those chambers, it is possible to simulate NIRS measurements and the final goal of the model will be accomplished. The following four assumptions are used to relate the state variables of the SimCVR model to the spatially distributed hemodynamics that are measured by NIRS.

The first assumption is that the concentration of hemoglobin in the blood M_HbT_, _blood is constant. This assumption simplifies the conversion of time-varying blood oxygen saturation S_r_,_O_₂, _j in the jth vascular chamber to time-varying concentration of oxyhemoglobin

(62)

and deoxyhemoglobin

(63)

where for the systemic region (r = s) we have j = la, sa, a, c, v, sv, lv, and for the cerebral region (r = b), j = pa, a, c, v, pv. Relating blood oxygen saturation to hemoglobin concentration captures the dynamics of the systemic and cerebral gas exchange models.

Next the blood volume dynamics from the circulation model are incorporated with the introduction of a second assumption, which is that a fractional change in blood volume causes a fractional change in effective hemoglobin concentration. Specifically, the fractional volume in the jth chamber in the systemic or cerebral region (r = s, b) is the state variable V_r_, _j divided by the initial volume V_r_, _j_,0. The effective hemoglobin concentration _r_,_h_, _j for h = HbO, HbR is then

(64)

where M_r_,_h_, _j is the hemoglobin concentration from Eq. (62) or (63) as appropriate. The effective hemoglobin concentration _r_,_h_, _j now includes the dynamics from the lumped gas exchange and lumped circulatory models but does not yet reflect the composition of the spatially distributed tissues.

Third and forth assumptions are now needed to incorporate tissue composition and complete the relationship between lumped parameters and distributed tissue properties in a local region. The third assumption is that an initial blood volume fraction F_r_,_bv can be used as a weighting factor in the systemic and cerebral regions (r = s, b) to account for the local blood-tissue composition. The forth assumption is that initial volume fractions of the vascular chambers F_r_, _j can account for the relative composition of the vascular chambers that comprises the local blood content of the tissue. Combining the third and forth assumptions relates the effective hemoglobin concentration in the vascular chambers _r_,_h_, _j to the spatially distributed hemoglobin concentration in the tissue

(65)

such that _r_,_h approximates the time-varying distributed concentration of hemoglobin species h in the systemic and cerebral regions (r = s, b). Since the volume fractions F are held constant, all the temporal variation in _r_,_h comes from dynamics in blood volume V_r_, _j and blood oxygen saturation S_r_,_O_₂, _j as computed by the lumped parameter components of the SimCVR model.

Simulating near-infrared spectroscopy measurements

Based on the model development thus far, the spatially distributed hemoglobin concentrations in the systemic and cerebral tissue regions are known. The present task is now to simulate NIRS measurements from these known hemoglobin concentration changes in the scalp (i.e. systemic region) and brain (i.e. cerebral region). Simulating NIRS will require three steps.

First the effective pathlength L_r_,_λ of the near-infrared light in each region r and each wavelength λ is estimated. For this step we used a Monte Carlo simulation of photon migration in the tissue, which is a numerical solution to radiative transport [43]. The Monte Carlo method allows for spatially varying optical properties and complex geometries. The resulting effective pathlength distributions characterize the spatial sensitivity for each source-detector pair to each tissue region. Details of how this Monte Carlo method is used to obtain pathlengths L_r_,_λ are already published [44; 45; 46].

In the second step, the calculated pathlengths L_r_,_λ are used in a first-order Rytov expansion [47] to model changes in optical density (ΔOD) that result from perturbations in the absorption coefficient Δμ_a_,_r_,_λ

(66)

where λ is the wavelength of the NIRS measurement. Next the absorption coefficient perturbations Δμ_a_,_r_,_λ are related to perturbations in the concentration of oxy- and deoxyhemoglobin

(67)

where ε_HbO_,_λ and ε_HbR_,_λ are the tabulated extinction coefficients [48]. Finally, the deviations in hemoglobin concentration are measured from an initial value

(68)

where _r_,_h is obtained from Eq. (65) and _r_,_h_,0 is calculated from steady state, which completes the relationship between NIRS measurements and the lumped model states of SimCVR.

Results

The SimCVR model dynamics are first illustrated in simulation by showing how the model responds to step increases in ABP and arterial blood CO₂. Next we describe the procedure and results of fitting the SimCVR model to NIRS experimental data from 10 human subjects. The complete SimCVR model is also compared to a blood pressure regression, a second order transfer function model, and simplified versions of SimCVR to assess the benefit of using the complete model.

Simulated Step Increase in Arterial Blood Pressure

We simulated a 10 mmHg increase in the blood pressure of the systemic large artery. The increase was modeled as one-half period of a raised cosine with a 1 s transition to the final value. The circulation and gas exchange dynamics are shown in the capillary chambers (Fig. 2).

Fig. 2

SimCVR circulation and blood gas response to a simulated step 10 mmHg increase in arterial blood pressure shown at the capillary level.

The following narrative qualitative describes the step response dynamics: The step increase in pressure in the systemic large artery P_s_,_la propagates through the neck arteries to the cerebral circulation where it causes an increase in cerebral arteriole pressure P_b_,_a and increased arteriole blood flow q_b_,_a. When q_b_,_a raises above the blood flow set point q_b_,_sp, the autoregulation state x_b_,_aut becomes positive [Eq. (28)]. The autoregulation state acts to decrease cerebral arteriole compliance C_b_,_a [Eq. (30)] because positive values of x_b_,_aut pull C_b_,_a towards the lower bound [Eq. (29b)] of the sigmoid where C_b_,_a =C_b_,_a_,0 − C_b_,_w. The effects of decreasing C_b_,_a are that cerebral arteriole volume V_b_,_a goes down [Eq. (18)] and arteriole resistance R_b_,_a goes up [Eq. (2)]. The increased R_b_,_a counteracts the higher P_b_,_a and lowers the arteriole flow q_b_,_a [Eq. (3)] towards the set point q_b_,_sp. This return of cerebral blood flow during a step increase in arterial pressure is the expected autoregulatory response.

Simulated Step Increase in Arterial Blood Carbon Dioxide

We also simulated the response to a 10 mmHg an increase in the partial pressure of blood carbon dioxide. The increase was again modeled as one-half period of a raised cosine but with a 30 s transition to the final value. The circulation and gas exchange dynamics are again shown in the capillary chambers (Fig. 3).

Fig. 3

SimCVR response to a simulated step 10 mmHg increase in blood carbon dioxide shown at the capillary level.

Relative to the pressure step response, the CO₂ response dynamics follow a somewhat more complicated sequence of events that arise from competing factors in the system equations as illustrated in the following narrative: The elevated arterial blood CO₂ flows into the cerebral capillary chamber and diffuses into the ISF [Eq. (60)]. The autoregulation blood flow set point q_b_,_sp increases with increased CO₂ in the ISF P_b_,_CO_₂,_ISF [Eqs. (25) and (26)]. Meanwhile, the autoregulation gain K_aut decreases [Eq. (27)]. Since the simulation started at steady state with all time derivatives equal to zero, the increased q_b_,_sp results in a negative derivative of the autoregulation state x_b_,_aut [Eq. (28)], which pulls its value negative. The lower K_aut just makes this happen more gradually. The autoregulation state acts on the cerebral arteriole compliance C_b_,_a [Eq. (30)] such that negative values of x_b_,_aut raise C_b_,_a towards the upper bound [Eq. (29a)] of the sigmoid where C_b_,_a =C_b_,_a_,0 + C_b_,_w. The effects of increasing C_b_,_a are that cerebral arteriole volume V_b_,_a goes up [Eq. (18)] and arteriole resistance R_b_,_a goes down [form of Eq. (2)]. Increasing V_b_,_a drives up the intracranial volume V_b_,_ic [Eq. (37)] and increases intracranial pressure P_b_,_ic [form of Eqs. (10) and (11)]. Considering the arteriole flow q_b_,_a [form of Eq. (3)], decreasing P_b_,_a will lower q_b_,_a and decreasing R_b_,_a will raise q_b_,_a. In our simulation, q_b_,_a drops slightly for about 15 seconds before increasing, resulting in an initial net outflow from the capillary. Now examining the capillary pressure P_b_,_c [Eqs. (24) and (22)], the initial negative derivative of capillary volume V_b_,_c would decrease P_b_,_c, except that it is dominated by the positive derivative of intracranial pressure P_b_,_ic. The increase in P_b_,_c drives up the capillary blood flow q_b_,_c until a new steady state is reached, which is the expected cerebral blood flow outcome in a hypercapnic condition. The small increase in capillary O₂ is also expected due to the interaction of CO₂ with oxyhemoglobin dissociation (see Blood oxygen in Appendix), and because the increased CBF will decrease the oxygen extraction fraction (OEF) since cerebral metabolic rate of oxygen (CMRO₂) is constant.

Fitting the Model to Near-Infrared Spectroscopy Data

The model fitting approach was to drive the circulation and gas exchange components of the SimCVR model with finger cuff measurements from human subjects and then to fit for the blood volume fractions that result in a best between simulated and actual NIRS data.

Baseline optical data was collected from approximately 12 healthy subjects, including 11 males and 1 female. The subjects’ mean age was 35 years old, with a standard deviation of 12 years as previously described [20]. The subjects were instructed to lie quietly in a darkened room and to breathe freely during repeated 5-minute runs. The raw optical measurements were collected with a continuous wave DOT instrument, demodulated and down sampled to 10 Hz. We concurrently monitored blood pressure variations with a custom-made finger cuff. The optical and blood pressure measurements were band-pass-filtered in a forward then reverse direction for zero phase distortion with a 4th order IIR Butterworth filter [49] with a pass band of 0.05 to 1.5 Hz. The intensity measurements were converted to a change in optical density for the measurement wavelengths λ = 830 and 690 nm

(69)

where I_λ_,0 is an initial value for the measurements, which we took as the median intensity.

Given the large quantity of data (50 detectors × 2 wavelengths × 12 subjects × 3 runs × 5 minutes), and the sensitivity of time-series modeling to noise and motion artifacts, we pruned the data set by the following procedure. Of the 12 subjects, 2 were dropped from the analysis due to more frequent motion artifacts in the blood pressure measurements as judged by visual inspection. For the remaining subjects, runs and detectors, we computed the power spectrum of the optical data and used the magnitude of the cardiac power relative to the background signal as a quality metric. We selected the runs with the highest cardiac power from each of the remaining 10 subjects. Then for each of the runs we selected the 10 detectors with the highest cardiac power. The resulting data subset contained 10 detectors, 2 wavelengths, and 10 subjects, each with 1 run of 5-minute duration and a concurrent finger-cuff pressure recording.

The finger cuff measurements were converted to inputs for the SimCVR model in the form of systemic small artery pressure derivative equation M78

as in Eq. (12) by taking finite difference derivative of the finger cuff pressure. The finite difference derivative of the pressure signal introduced negligible noise amplification because the attenuation of the band pass filter described above exceeds 120 dB at 4 Hz. The finger cuff measurements obtained from the human subjects were uncalibrated. Although this is a source of uncertainty in the analysis that follows, the fitting of blood volume fraction effectively compensates to a degree. The gain on the finger cuff measurements was manually adjusted after data collection to obtain physiologically plausible systolic and diastolic blood pressure fluctuations spanning a range of 80 to 120 mmHg.

We estimated the blood volume fractions F_r_,_bv and vascular compositions F_r_, _j separately for each detector by a constrained nonlinear optimization (FMINCON function in MATLAB®, MathWorks, Inc., Natick, MA) that minimized the objective function

(70)

The objective function f of Eq. (70) is a joint minimization of the sum-squared difference between the simulated and measured data at the wavelengths λ₁ = 890 nm and λ₂ = 690 nm. The summations in Eq. (70) are over time and the volume fractions F_r_,_bv and F_r_, _j act as regression coefficients that determine the linear combination of time-varying blood oxygen S_r_,_O_₂, _j and time-varying blood volume V_r_, _j that best predict the measured optical density data. The optimization is constrained, however, because blood volume fractions F_r_,_bv and vascular composition fractions F_r_, _j can be bounded to physiologically plausible ranges. We selected the FMINCON function because of it supports the use of constraints while minimizing an objective function.

Constraints on the blood volume fractions in the scalp and brain were guided by values from the literature. Normal blood volume fractions are approximately 4% in the skin [50] and 5.2 ± 1.4 % in the gray matter of the brain [51]. The arterial volume fraction in the brain is approximately 30% [52], meaning that capillary plus venous blood volume fraction is in the range of 70%. With these values as a general guideline and recognizing that the actual blood volume and composition may vary significantly depending on exactly where the NIRS detectors are placed relative to vascular structures, we choose the following bounds happens for the scalp blood volume fraction

(71)

and the brain blood volume fraction

(72)

The vascular composition bounds for the large arteries and large veins (j = la, lv) were

(73)

and in the smaller vessels

(74)

for r = s; j = sa, a, c, v, sv and r = b; j = pa, a, c, v, pv. The vascular composition fractions were also subject to the constraint that they must sum to 1 in the systemic and cerebral systems. By performing this fitting procedure separately on each detector, we accommodate local anatomical variations in the tissue composition. The mean and standard deviation of the fit values for the systemic blood volume fraction F_s_,_bv was 4% ± 1% and for the cerebral blood volume fraction F_b_,_bv it was 7% ± 3%. The results for the vascular composition F_r_, _j are summarized in Fig. 4. Future studies that include anatomical MRI could set these values explicitly. All other parameters in the SimCVR model were left constant across all subjects with the values reported in Table 3.

Fig. 4

Fit values for blood volume fraction (bvf) and the distribution of vascular compartments. The error bars show the standard deviation across measurement location after correcting for inter-subject variability.

Table 3

Parameter values

Control conditions for comparison with model fitting

Our primary control condition for comparison with the SimCVR model fit was to directly rescale the finger-cuff measurements to maximize the coefficient of determination R² with the NIRS data (BP fit). A short time series of the SimCVR model fit and BP fit is shown in Fig. 5 for the 830 and 690 nm NIRS measurements. The high frequency oscillations that are apparent in the NIRS data are the cardiac fluctuations; the lower frequency fluctuations are likely related to respiration, vasomotion and other factors. In this short time series, the R² values for the model fit are higher than for the BP fit, which is the average trend that we observed and summarize in Fig. 6. Our purpose in showing this short time series, however, is to qualitatively illustrate that the R² improvement in the SimCVR model fit reflects better agreement of the relative amplitude and phase of the high and low frequency physiological fluctuations. From this perspective, the SimCVR model is a nonlinear filter of the finger-cuff measurements that results in improved correlations with the baseline NIRS data. The results summary (Fig. 6) shows that improved R² values were found in 9 out of 10 subjects for the 830 and 690 nm NIRS measurements. The higher R² values appear consistently higher at 830 nm relative to 690 nm. This wavelength dependence is potentially related to the dominance of HbO absorption at 830 nm and HbR absorption at 690 nm.

Fig. 5

Example of data fit results from subject 7, detector 5 for the finger-cuff pressure (top row) and SimCVR model (bottom row) and for measurements at 830 nm (left column) and 690 nm measurements (right column). The R² values from the 5-minute runs are reported (more ...)

Fig. 6

Within-subject comparison of R² values for the full SimCVR model and finger-cuff fits. We performed 2-way ANOVA on Fisher z transformed R² values to test significance. Differences at the level p < 0.0001 are indicated by **. The p < 0.01 (more ...)

Next we explored the possibility that the improved R² values from the SimCVR fit relative to the BP fit were due to the time-dynamics of the system that could be captured with a simple transfer function (TF fit). Based on published methods for transfer function analysis on cerebral autoregulation [53], we chose to fit a 2^nd order continuous-time system transfer function G of the form

(75)

where s is the Laplace variable, K_p is the static gain, τ_z is the time constant of the zero, τ_w is the resonance time constant, and ζ is the damping factor. We fit the transfer function parameters in Eq. (75) using prediction error minimization (PEM function in MATLAB® System Identification Toolbox). The overall results of the TF fit (Fig. 7) show significant improvement over the BP fit. For the 830 nm data, the SimCVR model fit is superior to the TF fit at the level of p < 0.0001, which suggests that the SimCVR captures temporal-dynamics beyond the complexity of the transfer function. For the 690 nm data, however, the TF fit and SimCVR fit performed equally well, suggesting that perhaps the transfer function sufficiently captures HbR dynamics, which dominate the NIRS measurements at this wavelength.

Fig. 7

Overall comparison of R² values for the full SimCVR model and control conditions. We performed 3-way ANOVA on Fisher z transformed R² values to test significance. Differences at the level p < 0.0001 are marked **. The error bars show the standard (more ...)

Two additional control conditions were used to explore the potential for over-fitting with the SimCVR model. A third control condition (Systemic only) forces all the time derivatives in the cerebral region to zero, which effectively turns off the cerebral region so that all the measurement dynamics are fit only with the systemic (i.e. scalp) compartment. The forth condition turns off the systemic region leaving the dynamics to be fit only by the cerebral compartment (Cerebral only). The overall R² values from all four control conditions (TF fit, BP fit, Systemic Only, Cerebral only) are compared with SimCVR in Fig. 7. We found for 830 nm NIRS measurements that the systemic only fit was better than the cerebral only. The SimCVR model fit was better than any of the four control conditions at 830 and better than all but the TF fit at 690 nm (multifactor ANOVA, p<0.0001), which supports our hypothesis that the SimCVR model captures some salient features of the systemic and cerebral physiology.

Discussion

We have presented an integrated systemic and cerebral circulation and gas exchange model with a number of improvements to prior published work. The improvements include explicit inclusion of bound and unbound gases, nonlinear resistances, Windkessel compliances, and a pressure input in the small systemic artery. We have also used the steady-state relationships to calculate many of the required parameters and included a complete description of the state equations. The model was first illustrated in simulation and then used to predict NIRS data. Statistical comparisons were made between the SimCVR model predictions, blood pressure regression and transfer function analysis.

Our principal findings are the following: i) The SimCVR model simulation appears to reproduce the main characteristics of cerebral autoregulation (Fig. 2), where the cerebral blood flow returns to its baseline value on the order of seconds after a step increase in systemic arterial blood pressure. The effects of this regulatory physiology are also reflected through the coupling between the flow and gas exchange systems (Fig. 3). ii) The R² goodness-of-fit to NIRS data was better with the SimCVR model than with blood pressure regression, transfer function analysis, or reduced versions of the model. iii) The time-series predictions from the SimCVR model appear to more closely match the relative amplitude and phase of the cardiac and respiratory fluctuations in the NIRS data compared to blood pressure alone (Fig. 5).

Although many issues remain to be addressed in this model with regards to validation, un-modeled phenomena, and the potential of over-fitting, this version of SimCVR is a significant advancement towards physiological model-based interpretation of neuroimaging data. In the following discussion, we comment on the issues of validation, model limitations, and parameter fitting. We conclude by mentioning some of the potential research and clinical applications of this work.

Validation issues

The data that we used to compare model predictions with observed human physiology includes a high degree of spatial averaging due to the photon scattering in the NIRS biophysics. The variations in the measurements arise from a mixture of scalp and brain physiology and, from a modeling perspective, some unknown weighting of the vascular chambers. Given these limitations, our outcome measure is an aggregate outcome from all of these regions rather than a validation of the relative attenuation and phase of the vascular chambers separately. The bottom line is that additional data will be required to validate the model.

Although the current data comparison that we present falls short of validation, the improvement in R² fits between the data and the model with the whole model versus control conditions (Fig. 7) is consistent with the hypothesis that the SimCVR model captures some salient features of real human physiology.

Another validation issue is with the pressure input to the small arteries in the model [Eqs. (12) through (15)]. Instead of acquiring a Finapres monitor (Finapres Medical Systems, Amsterdam, Netherlands), we built our own custom finger cuff pressure monitor. The Finapres uses online servo control of the air pressure in the finger cuff and an assumed finger to aorta transfer function to estimate arterial pressure. Our finger cuff uses manually pressurized water, which was then fixed and included no transfer function for preprocessing. Lacking the online servo pressure control increases the level of uncertainly in our pressure readings and we are, for now, tolerating this error source. Rather than use a separate finger-aorta transfer function, we have used the SimCVR model with a small artery pressure derivative input. While we hypothesize that using the systemic circulation model for this transformation is a sound approach, our large artery pressure estimates are not validated at present. Although the large artery pressure estimates are a source of uncertainty, the model analysis and results rely predominately on the temporal dynamics of the finger cuff data and not on the absolute scale.

Future validation studies will ideally include measurement of intracranial pressure and high-resolution anatomical scans that segment the brain anatomy including vascular composition. It would also be desirable to obtain measures of blood flow velocity with transcranial Doppler (TCD) [54] and possibly also a flow-sensitive MRI technique like arterial spin labeling (ASL) [55].

Limitations of the model

The challenge of building a mathematical model is to balance its predictive capability and its simplicity to enable interpretation of the physical system dynamics. There must also be a balance between the scope of questions posed and comprehensiveness of the model. In this work we have chosen 48 states and 92 lumped parameters that represent millions of physical elements and dynamic relationships. Although by modeling standards this is moderately complex, the physical system is vastly more complex. The complexity of SimCVR compared to other models allows for testing of specific hypotheses about brain physiology. In addition, SimCVR serves as a bridge between the simplified transfer function model [53] and the millions of elements that comprise physiological reality.

In the current version of this model there are still some entire categories of physical phenomena that have been ignored such as vasomotion [56], the anatomy of the cerebral vessels [57], and neurovascular coupling [58]. This fact of missing features is apparent even when the model predictions appear to be working very well. For example in Fig. 5, there is a dip in the observed NIRS measurement in the 3–4 second time range that is completely missed by both the blood pressure and SimCVR predictions. It is possible that this dip arose from local vasomotion rather than a systemic blood pressure fluctuation that passed through the autoregulated cerebral circulation.

Perhaps the most significant limitation of the current SimCVR model is the use of lumped parameters to simulate spatially distributed hemodynamics. This limitation would become particularly apparent if applied to voxel-wise fMRI data. As a whole-brain lumped model, SimCVR would only be useful for capturing whole-brain signal dynamics in fMRI. Additional representation of the cerebral vascular structure would be required to describe regional differences.

Parameter fitting

We limited our parameter fitting to the initial blood volume fractions and vascular compositions for the systemic and cerebral regions. The number of fit parameters totaled 8 for the systemic region and 6 for the cerebral region. The parameter values were also bounded to physiologically plausible ranges; they also did not very in time, and did not alter the dynamics of the flow or gas exchange systems. The resulting fit values vary over a fairly wide range (Fig. 4) but have, on average, sufficiently reasonable values for the present crude level of validation.

The ratio of fit parameters to data points was 14 to 6000. We did not, however, attempt to assess the actual degrees of freedom in light of the temporal autocorrelation, correlation between blood pressure and NIRS measurements and the nonlinear, bounded parameter search. Instead we tested the possibility for the model fit to have a lower R² value compared to four control conditions. We also tried altering some steady-state parameter values non-physiological values and re-fit the simulations to the data. The R² values from SimCVR dropped significantly (results not shown). Since there are infinite possible combinations of parameter values for the dynamic system, we chose to illustrate the a similar decrease in the goodness-of-fit by shutting off first the cerebral loop and then the systemic loop prior to fitting the blood volume fractions and vascular compositions (Fig. 7). These findings help build confidence that we are not over-fitting.

A limitation of the present parameter fitting approach in which only the blood volume and vascular composition fractions are fit is that there is unlikely to be any one-size-fits-all dynamic model for all human physiology. Individual differences can be significant and then there could also be changes during development, aging, disease and so forth. However, we were not confident in the present study that available data would support unique parameter fits with 92 dynamic parameters plus the 14 we fit already. A move towards fitting more of the parameters will require more data sources beyond the optical brain and finger cuff pressures. Anatomical MRI, for example, could help to fix the arterial and venous volume fractions in the brain regions that are measured optically. Multimodality studies involving concurrent measurements with arterial spin labeling (ASL) MRI and NIRS [59; 60] would help to resolve uncertainties in the blood flow and volume dynamics.

Potential research and clinical applications

The application of the present model to NIRS data is a significant step forward in our longer-term aim of model-based neuroimaging analysis. Although further validation studies are required, we can already begin to use this model to examine the contributions of baseline physiology to both NIRS and, with some addition to the forward model, fMRI as well. We can also now pursue pilot studies in clinical data interpretation by applying restricted and selective parameter fits to neuroimaging data from different patient populations such as stroke, cerebrovascular disease, coma, and traumatic brain injury (TBI). In the future, we aim to fit model parameters such as vascular compliance as potential biomarkers of cerebrovascular disease processes with the hope that the measured dynamics provide sufficient sensitivity. Perhaps the most immediate application, however, is to simulate ground truth for analysis algorithm development in the functional neuroimaging field. This application could be particularly useful for comparing time-series analysis methods of the baseline physiology. The simulated ground truth data, if combined with a multi-compartment Windkessel model of the functional activation and spatially mapped, could also be used to compare different tomography and functional analysis methods.

Table 2

Glossary of terms

Acknowledgments

The authors wish to thank M. A. Francescini for the data, T. J. Huppert for consultation, and S. R. Arridge for an early conceptualization of this project. Funding was provided by NIH T32-CA09502, P41-RR14075, R01-EB002482, R01-EB001954, and the MIND Institute. KLP received support from an NSF Graduate Research Fellowship.

Appendix

Blood oxygen

The Kelman equation for oxygemoglobin dissociation uses the coefficients

(76)

The calculation starts with a scale factor v, which accounts for changes in temperature, pH, and the partial pressure of CO₂

(77)

Next the following empirical relationship is used to compute hemoglobin saturation

(78)

We have further differentiated Kelman’s equation with respect to the partial pressure of oxygen as needed for the effective solubility _O_₂ in Eq. (45)

(79)

Blood carbon dioxide

First and effective pH is calculated to account for differences in temperature T

(80)

Now the concentration ratio of oxyhemoglobin-bound CO₂ to plasma CO₂ is estimated as

(81)

and the concentration ratio of deoxyhemoglobin-bound CO₂ to plasma CO₂ is

(82)

The unbound CO₂ concentration is computed with the Henderson-Hasselbalch equation

(83)

and the bound CO₂ concentration accounts for hemoglobin binding with Dox and Dr

(84)

Finally, we will take the derivative with respect to the partial pressure of CO₂

(85)

which is the effective solubility of carbon dioxide in the blood _CO_₂ in Eq. (47).

Parameter relationships at steady state

One of our objectives with this model was to obtain self-consistency among the parameters and system equations at steady state. Achieving this objective requires specifying a set of constraints on the parameters. The general approach to determining these constraints is to set all the time derivatives in the state equations equal to zero. The specific form of these constraint equations requires a somewhat judicious choice of which parameters to set and which to calculate. Although our fixed parameter set is not unique, it follows a sequence that we found to be logical and that permits direct algebraic calculation of all but two parameter pairs, which must be estimated with nonlinear searches. By also following simple rules on the monotonic variation in fluid and gas pressure throughout the system, our approach will always yield positive quantities for the diffusion constants and cerebral metabolic rates that are positive for O₂ and negative for CO₂.

The initial partial pressures of O₂ and CO₂ in the systemic capillary are computed with a nonlinear search that minimizes the objective function

(86)

where the gas concentrations are calculated from the oxyhemoglobin dissociation, blood oxygen and blood carbon dioxide functions described in the Appendix. Since there is only one capillary chamber for the systemic system

(87)

with g = O₂, CO₂. The partial pressure of O₂ should monotonically decrease from artery to capillary to ISF to ICF, whereas the partial pressure of CO₂ should monotonically increase along the same cascade. The diffusion constants can then be calculated as

(88)

(89)

In the systemic fluid flow system, the initial blood pressures should monotonically decrease through the system and all the initial compliance and resistances are calculated from the initial values of pressure, volume and flow

(90)

(91)

where F_s_,_us is unstressed volume fraction, j = la, sa, sv, lv and k = sa, a, lv, ra. For the Windkessel compliance compartments, j = a, c, v, and k = c, v, sv, the term F_s_,_us is neglected.

In the cerebral fluid flow system, the initial blood pressures again decrease monotonically through the system with intracranial ic pressure also greater than the neck vein nv pressure. The initial compliances and resistances are, in the neck arteries and veins

(92)

(93)

(94)

(95)

For j = pa, pv and k = a, nv we have

(96)

(97)

and for the Windkessel compartments, j = a, c, v and k = c, v, pv, the unstressed volume fraction F_b_,_us is neglected. The initial intracranial volume is

(98)

with the compliance

(99)

The resistances for the CSF flow are

(100)

(101)

For the cerebral gas exchange, there is again only one capillary chamber therefore with g = O₂, CO₂,

(102)

The partial pressure of O₂ should monotonically decrease from artery to capillary to CSF to ISF to ICF, whereas the partial pressure of CO₂ should monotonically increase along the same cascade. The initial partial pressures of O₂ and CO₂ in the cerebral artery are found with a nonlinear search that minimizes the objective function

(103)

where the gas concentrations are calculated from the Oxyhemoglobin dissociation, Blood oxygen and Blood carbon dioxide functions described in the Appendix. The metabolic rates that balance the system at steady state are

(104)

and the diffusion constants are

(105)

(106)

(107)

Finally, the autoregulation state variable x_b_,_aut should be initialized at zero.

Footnotes

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

1. Liu CH, Niranjan SC, Clark JW, Jr, San KY, Zwischenberger JB, Bidani A. Airway mechanics, gas exchange, and blood flow in a nonlinear model of the normal human lung. J Appl Physiol. 1998;84:1447–69. [PubMed]

2. Lu K, Clark JW, Jr, Ghorbel FH, Ware DL, Bidani A. A human cardiopulmonary system model applied to the analysis of the Valsalva maneuver. Am J Physiol Heart Circ Physiol. 2001;281:H2661–79. [PubMed]

3. Mukkamala R, Cohen RJ. A forward model-based validation of cardiovascular system identification. American Journal of Physiology-Heart and Circulatory Physiology. 2001;281:H2714–H2730. [PubMed]

4. Lu K, Clark JW, Ghorbel FH, Robertson CS, Ware DL, Zwischenberger JB, Bidani A. Cerebral autoregulation and gas exchange studied using a human cardiopulmonary model. American Journal of Physiology-Heart and Circulatory Physiology. 2004;286:H584–H601. [PubMed]

5. Olufsen MS, Nadim A, Lipsitz LA. Dynamics of cerebral blood flow regulation explained using a lumped parameter model. Am J Physiol Regul Integr Comp Physiol. 2002;282:R611–22. [PubMed]

6. Ursino M, Lodi CA. Interaction among autoregulation, CO2 reactivity, and intracranial pressure: a mathematical model. Am J Physiol. 1998;274:H1715–28. [PubMed]

7. Tiecks FP, Lam AM, Aaslid R, Newell DW. Comparison of static and dynamic cerebral autoregulation measurements. Stroke. 1995;26:1014–1019. [PubMed]

8. Panerai RB, White RP, Markus HS, Evans DH. Grading of cerebral dynamic autoregulation from spontaneous fluctuations in arterial blood pressure. Stroke. 1998;29:2341–2346. [PubMed]

9. van Beek AHEA, Claassen JAHR, Rikkert MGMO, Jansen RWMM. Cerebral autoregulation: an overview of current concepts and methodology with special focus on the elderly. Journal of Cerebral Blood Flow and Metabolism. 2008;28:1071–1085. [PubMed]

10. Panerai RB. Cerebral autoregulation: From models to clinical applications. Cardiovascular Engineering. 2008;8:42–59. [PubMed]

11. Buxton RB, Uludag K, Dubowitz DJ, Liu TT. Modeling the hemodynamic response to brain activation. Neuroimage. 2004;23(Suppl 1):S220–33. [PubMed]

12. Mintun MA, Lundstrom BN, Snyder AZ, Vlassenko AG, Shulman GL, Raichle ME. Blood flow and oxygen delivery to human brain during functional activity: theoretical modeling and experimental data. Proc Natl Acad Sci U S A. 2001;98:6859–64. [PubMed]

13. Fantini S. A haemodynamic model for the physiological interpretation of in vivo measurements of the concentration and oxygen saturation of haemoglobin. Phys Med Biol. 2002;47:N249–57. [PubMed]

14. Behzadi Y, Liu TT. An arteriolar compliance model of the cerebral blood flow response to neural stimulus. Neuroimage. 2005;25:1100–11. [PubMed]

15. Buxton RB, Wong EC, Frank LR. Dynamics of blood flow and oxygenation changes during brain activation: The balloon model. Magnetic Resonance in Medicine. 1998;39:855–864. [PubMed]

16. Mandeville JB, Marota JJA, Ayata C, Zaharchuk G, Moskowitz MA, Rosen BR, Weisskoff RM. Evidence of a cerebrovascular postarteriole windkessel with delayed compliance. Journal of Cerebral Blood Flow and Metabolism. 1999;19:679–689. [PubMed]

17. Stefanovic B, Warnking JM, Rylander KM, Pike GB. The effect of global cerebral vasodilation on focal activation hemodynamics. Neuroimage. 2006;30:726–34. [PubMed]

18. Banaji M, Tachtsidis I, Delpy D, Baigent S. A physiological model of cerebral blood flow control. Math Biosci. 2005;194:125–73. [PubMed]

19. Payne SJ. A model of the interaction between autoregulation and neural activation in the brain. Math Biosci. 2006;204:260–81. [PubMed]

20. Franceschini MA, Joseph DK, Huppert TJ, Diamond SG, Boas DA. Diffuse optical imaging of the whole head. Journal of Biomedical Optics. 2006 Accepted for publication.

21. Lund TE, Madsen KH, Sidaros K, Luo WL, Nichols TE. Non-white noise in fMRI: does modelling have an impact? Neuroimage. 2006;29:54–66. [PubMed]

22. Thomas CG, Harshman RA, Menon RS. Noise reduction in BOLD-based fMRI using component analysis. Neuroimage. 2002;17:1521–37. [PubMed]

23. Deckers RH, van Gelderen P, Ries M, Barret O, Duyn JH, Ikonomidou VN, Fukunaga M, Glover GH, de Zwart JA. An adaptive filter for suppression of cardiac and respiratory noise in MRI time series data. Neuroimage. 2006;33:1072–81. [PMC free article] [PubMed]

24. Diamond SG, Huppert TJ, Kolehmainen V, Franceschini MA, Kaipio JP, Arridge SR, Boas DA. Dynamic physiological modeling for functional diffuse optical tomography. Neuroimage. 2006;30:88–101. [PMC free article] [PubMed]

25. Sloan MA, Alexandrov AV, Tegeler CH, Spencer MP, Caplan LR, Feldmann E, Wechsler LR, Newell DW, Gomez CR, Babikian VL, Lefkowitz D, Goldman RS, Armon C, Hsu CY, Goodin DS. Assessment: Transcranial Doppler ultrasonography - Report of the Therapeutics and Technology Assessment Subcommittee of the American Academy of Neurology. Neurology. 2004;62:1468–1481. [PubMed]

26. Fu CH, Yang CCH, Kuo TBJ. Effects of different classes of antihypertensive drugs on cerebral hemodynamics in elderly hypertensive patients. American Journal of Hypertension. 2005;18:1621–1625. [PubMed]

27. Reinhard M, Muller T, Guschlbauer B, Timmer J, Hetzel A. Transfer function analysis for clinical evaluation of dynamic cerebral autoregulation- a comparison between spontaneous and respiratory-induced oscillations. Physiological Measurement. 2003;24:27–43. [PubMed]

28. Haubrich C, Kohnke A, Diehl RR, Moller-Hartmann W, Klotzsch C. Impact of vertebral artery disease on dynamic cerebral autoregulation. Acta Neurologica Scandinavica. 2005;112:309–316. [PubMed]

29. Ramos EG, Simpson DM, Panerai RB, Nadal J, Lopes JMA, Evans DH. Objective selection of signals for assessment of cerebral blood flow autoregulation in neonates. Physiological Measurement. 2006;27:35–49. [PubMed]

30. Reinhard M, Roth M, Guschlbauer B, Harloff A, Timmer J, Czosnyka M, Hetzel A. Dynamic cerebral autoregulation in acute ischemic stroke assessed from spontaneous blood pressure fluctuations. Stroke. 2005;36:1684–1689. [PubMed]

31. Hanon O, Haulon S, Lenoir H, Seux ML, Rigaud AS, Safar M, Girerd X, Forette F. Relationship between arterial stiffness and cognitive function in elderly subjects with complaints of memory loss. Stroke. 2005;36:2193–2197. [PubMed]

32. Qiu CX, Winblad B, Viitanen M, Fratiglioni L. Pulse pressure and risk of Alzheimer disease in persons aged 75 years and older - A community-based, longitudinal study. Stroke. 2003;34:594–599. [PubMed]

33. Fulesdi B, Limburg M, Bereczki D, Michels RPJ, Neuwirth G, Legemate D, Valikovics A, Csiba L. Impairment of cerebrovascular reactivity in long-term type 1 diabetes. Diabetes. 1997;46:1840–1845. [PubMed]

34. Fulesdi B, Limburg M, Bereczki D, Kaplar M, Molnar C, Kappelmayer J, Neuwirth G, Csiba L. Cerebrovascular reactivity and reserve capacity in type II diabetes mellitus. Journal of Diabetes and Its Complications. 1999;13:191–199. [PubMed]

35. Snyder MF, Rideout VC. Computer simulation studies of venous circulation. IEEE Transactions on Biomedical Engineering. 1969;Bm16:325-&. [PubMed]

36. Huppert TJ, Allen MS, Benav H, Jones PB, Boas DA. A multicompartment vascular model for inferring baseline and functional changes in cerebral oxygen metabolism and arterial dilation. Journal of Cerebral Blood Flow and Metabolism. 2007;27:1262–1279. [PMC free article] [PubMed]

37. Boas DA, Jones SR, Devor A, Huppert TJ, Dale AM. Avascular anatomical network model of the spatio-temporal response to brain activation. Neuroimage. 2008;40:1116–1129. [PMC free article] [PubMed]

38. Mountcastle VB. Medical physiology. C. V. Mosby Co; St. Louis: 1980.

39. Monro J. Observations on the Structure and Functions of the Nervous System. W. Creech; Edinburgh: 1783.

40. Kellie G. An account …, with some reflections on the pathology of the brain. Edinburgh Med Chir Soc Trans. 1824;1:84–169.

41. Kelman GR. Digital Computer Subroutine for Conversion of Oxygen Tension into Saturation. Journal of Applied Physiology. 1966;21:1375-&. [PubMed]

42. Kelman GR. Digital computer procedure for the conversion of PCO2 into blood CO2 content. Respir Physiol. 1967;3:111–5. [PubMed]

43. Boas DA, Culver JP, Stott JJ, Dunn AK. Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head. Optics Express. 2002;10:159–170. [PubMed]

44. Boas DA, Dale AM. Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function. Applied Optics. 2005;44:1957–1968. [PubMed]

45. Diamond SG, Huppert TJ, Kolehmainen V, Franceschini MA, Kaipio JP, Arridge SR, Boas DA. Physiological system identification with the Kalman filter in diffuse optical tomography. Med Image Comput Comput Assist Interv Int Conf Med Image Comput Comput Assist Interv. 2005;8:649–56. [PubMed]

46. Strangman G, Franceschini MA, Boas DA. Factors affecting the accuracy of near-infrared spectroscopy concentration calculations for focal changes in oxygenation parameters. Neuroimage. 2003;18:865–879. [PubMed]

47. Ntziachristos V, Chance B, Yodh AG. Differential diffuse optical tomography. Optics Express. 1999;5:230–242. [PubMed]

48. Prahl S. Tabulated molar extinction coefficient for hemoglobin in water. Oregon Medical Laser Center; 2005.

49. Rabiner LR, Gold B. Theory and application of digital signal processing. Prentice-Hall; Englewood Cliffs, N.J: 1975.

50. Meglinski IV, Matcher SJ. Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions. Physiological Measurement. 2002;23:741–753. [PubMed]

51. Leenders KL, Perani D, Lammertsma AA, Heather JD, Buckingham P, Healy MJR, Gibbs JM, Wise RJS, Hatazawa J, Herold S, Beaney RP, Brooks DJ, Spinks T, Rhodes C, Frackowiak RSJ, Jones T. Cerebral blood-flow, blood-volume and oxygen utilization - normal values and effect of age. Brain. 1990;113:27–47. [PubMed]

52. Ito H, Kanno I, Iida H, Hatazawa J, Shimosegawa E, Tamura H, Okudera T. Arterial fraction of cerebral blood volume in humans measured by positron emission tomography. Annals of Nuclear Medicine. 2001;15:111–116. [PubMed]

53. Payne SJ, Tarassenko L. Combined transfer function analysis and modelling of cerebral autoregulation. Annals of Biomedical Engineering. 2006;34:847–858. [PubMed]

54. Aaslid R, Lindegaard KF, Sorteberg W, Nornes H. Cerebral autoregulation dynamics in humans. Stroke. 1989;20:45–52. [PubMed]

55. Detre JA, Zhang W, Roberts DA, Silva AC, Williams DS, Grandis DJ, Koretsky AP, Leigh JS. Tissue specific perfusion imaging using arterial spin labeling. NMR Biomed. 1994;7:75–82. [PubMed]

56. Nilsson H, Aalkjaer C. Vasomotion: mechanisms and physiological importance. Mol Interv. 2003;3:79–89. 51. [PubMed]

57. Duvernoy HM, Delon S, Vannson JL. Cortical blood vessels of the human brain. Brain Res Bull. 1981;7:519–79. [PubMed]

58. D’Esposito M, Deouell LY, Gazzaley A. Alterations in the BOLD fMRI signal with ageing and disease: a challenge for neuroimaging. Nat Rev Neurosci. 2003;4:863–72. [PubMed]

59. Huppert TJ, Hoge RD, Diamond SG, Franceschini MA, Boas DA. A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans. Neuroimage. 2006;29:368–82. [PMC free article] [PubMed]

60. Hoge RD, Franceschini MA, Covolan RJ, Huppert T, Mandeville JB, Boas DA. Simultaneous recording of task-induced changes in blood oxygenation, volume, and flow using diffuse optical imaging and arterial spin-labeling MRI. Neuroimage. 2005;25:701–7. [PubMed]