Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2720139

Formats

Article sections

Authors

Related links

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

Published in final edited form as:

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

PMCID: PMC2720139

NIHMSID: NIHMS117087

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

See other articles in PMC that cite the published article.

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.

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_{2} 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^{2} 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.

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_{2} and CO_{2} 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.

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).

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

$${V}_{s,j}-{V}_{s,j,us}={P}_{s,j}{C}_{s,j,0},$$

(1)

$${R}_{s,j}={R}_{s,j,0}{\left(\frac{{V}_{s,j,0}}{{V}_{s,j}}\right)}^{2},$$

(2)

and

$${q}_{s,j}=\frac{{P}_{s,j}-{P}_{s,k}}{{R}_{s,j}}.$$

(3)

Eq. (1) specifies the compliance constant *C _{s}*

$${V}_{s,j,us}={F}_{s,us}{V}_{s,j,0}.$$

(4)

Following [35], the compliance of the collapsible veins is increased 20 fold when the transmural pressure *P _{s}*

Calculating the system response with an ODE solver requires state equations that relate the time derivatives of the pressure
$\frac{d{P}_{s,j}}{dt}$ and volume
$\frac{d{V}_{s,j}}{dt}$ to the pressure *P _{s}*

$$\frac{d{P}_{s,j}}{dt}=\frac{1}{{C}_{s,j,0}}\frac{d{V}_{s,j}}{dt},$$

(5)

where the change in volume is the difference between blood outflow *q _{s}*

$$\frac{d{V}_{s,j}}{dt}={q}_{s,i}-{q}_{s,j}.$$

(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

$${C}_{s,j}={C}_{s,j,0}{\left(\frac{{P}_{s,j}}{{P}_{s,j,0}}\right)}^{{\scriptstyle \frac{1}{{\beta}_{s,j}}}-1},$$

(7)

where *β _{s}*

$${V}_{s,j}={P}_{s,j}{C}_{s,j}$$

(8)

is the pressure-volume relationship with time-varying compliance *C _{s}*

$$\frac{{dV}_{s,j}}{dt}={P}_{s,j}\frac{{dC}_{s,j}}{d{P}_{s,j}}\frac{{dP}_{s,j}}{dt}+{C}_{s,j}\frac{d{P}_{s,j}}{dt},$$

(9)

where

$$\frac{d{C}_{s,j}}{d{P}_{s,j}}=\frac{{C}_{s,j}}{{P}_{s,j,0}}\left(\frac{1}{{\beta}_{s,j}}-1\right){\left(\frac{{P}_{s,j}}{{P}_{s,j,0}}\right)}^{{\scriptstyle \frac{1}{{\beta}_{s,j}}}-2}.$$

(10)

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

$$\frac{d{P}_{s,j}}{dt}=\frac{d{V}_{s,j}}{dt}{\left({C}_{s,j}+\frac{{V}_{s,j}}{{C}_{s,j}}\frac{d{C}_{s,j}}{d{P}_{s,j}}\right)}^{-1}.$$

(11)

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

$$\frac{d{P}_{s,la}}{dt}=\frac{d{P}_{s,la,\mathit{input}}}{dt}$$

(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}*

$${q}_{s,la}=\frac{{P}_{s,la}-{P}_{s,sa}}{{R}_{s,la,0}{\left(\frac{{V}_{s,la,0}}{{P}_{s,la}{C}_{s,la,0}}\right)}^{2}}.$$

(13)

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

$$\frac{d{P}_{s,sa}}{dt}=\frac{{q}_{s,la}-{q}_{s,sa}}{{C}_{s,sa,0}}.$$

(14)

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

$$\begin{array}{c}{P}_{s,la}^{3}-\left({P}_{s,sa}-2\frac{{V}_{s,la,us}}{{C}_{s,la,0}}\right){P}_{s,la}^{2}-{V}_{s,la,us}\left(2\frac{{P}_{s,sa}}{{C}_{s,la,0}}-\frac{{V}_{s,la,us}}{{C}_{s,la,0}^{2}}\right){P}_{s,la}+\cdots \\ -{R}_{s,la,0}\frac{{V}_{s,la,0}^{2}}{{C}_{s,la,0}^{2}}\left(\frac{d{P}_{s,sa}}{dt}{C}_{s,sa,0}+{q}_{s,sa}+\frac{{P}_{s,sa}{V}_{s,la,us}^{2}}{{R}_{s,la,0}{V}_{s,la,0}^{2}}\right)=0,\end{array}$$

(15)

where unknown value of *P _{s}*

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

$${P}_{b,na}={P}_{s,la}.$$

(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

$$\frac{d{P}_{b,nv}}{dt}=\frac{{q}_{b,pv}+{q}_{b,{f}_{2}}-{q}_{b,nv}}{{C}_{b,nv,0}},$$

(17)

and is derived in the same manner as Eqs. (5) and (6) where *q _{b}*

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}*

$${V}_{b,j}-{V}_{b,j,us}=({P}_{b,j}-{P}_{b,ic}){C}_{b,j,0}.$$

(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}*

$${V}_{b,j}=({P}_{b,j}-{P}_{b,ic}){C}_{b,j},$$

(19)

where *j* = *c*, *v* and

$${C}_{b,j}={C}_{b,j,0}{\left(\frac{{P}_{b,j}-{P}_{b,ic}}{{P}_{b,j,0}-{P}_{b,ic,0}}\right)}^{{\scriptstyle \frac{1}{{\beta}_{b,j}}}-1}.$$

(20)

The state equation for *P _{b}*

$$\frac{d{V}_{b,j}}{dt}=({P}_{b,j}-{P}_{b,ic})\left(\frac{d{C}_{b,j}}{d{P}_{b,j}}\frac{d{P}_{b,j}}{dt}+\frac{d{C}_{b,j}}{d{P}_{b,ic}}\frac{d{P}_{b,ic}}{dt}\right)+{C}_{b,j}\left(\frac{d{P}_{b,j}}{dt}-\frac{d{P}_{b,ic}}{dt}\right),$$

(21)

where

$$\frac{d{C}_{b,j}}{d{P}_{b,j}}=\frac{{C}_{b,j}}{{P}_{b,j,0}-{P}_{b,ic,0}}\left(\frac{1}{{\beta}_{b,j}}-1\right){\left(\frac{{P}_{b,j}-{P}_{b,ic}}{{P}_{b,j,0}-{P}_{b,ic,0}}\right)}^{{\scriptstyle \frac{1}{{\beta}_{b,j}}}-2}$$

(22)

and

$$\frac{d{C}_{b,j}}{d{P}_{b,ic}}=-\frac{d{C}_{b,j}}{d{P}_{b,j}}.$$

(23)

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

$$\frac{d{P}_{b,j}}{dt}=\frac{d{V}_{b,j}}{dt}{\left({C}_{b,j}+\frac{{V}_{b,j}}{{C}_{b,j}}\frac{d{C}_{b,j}}{d{P}_{b,j}}\right)}^{-1}+\frac{d{P}_{b,ic}}{dt}.$$

(24)

The expressions for resistance *R _{b}*

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}*

$${F}_{\mathit{aut},sp}=\frac{{P}_{b,C{O}_{2},\mathit{ISF}}-{P}_{b,C{O}_{2},\mathit{ISF},0}}{{P}_{b,C{O}_{2},\mathit{ISF},0}},$$

(25)

which correlates positively with the CBF set point

$${q}_{b,sp}={q}_{b,0}(1+{F}_{\mathit{aut},sp}),$$

(26)

and negatively correlates with the autoregulation reactivity gain

$${K}_{\mathit{aut}}={K}_{\mathit{aut},0}(1-{F}_{\mathit{aut},sp}).$$

(27)

The autoregulation state variable *x _{b}*

$$\frac{d{x}_{b,\mathit{aut}}}{dt}=\frac{{K}_{\mathit{aut}}}{{\tau}_{\mathit{aut}}}\left(\frac{{q}_{b,a}-{q}_{b,sp}}{{q}_{b,sp}}\right).$$

(28)

The variable cerebral arteriole compliance operates on a non-symmetric sigmoid with the bounds *C _{b}*

$$\begin{array}{c}{C}_{b,w}={C}_{b,{w}_{1}}\\ {k}_{b,p}={C}_{b,{w}_{1}}/2\end{array}\}\phantom{\rule{0.16667em}{0ex}}\text{when}\phantom{\rule{0.16667em}{0ex}}{x}_{b,\mathit{aut}}<0,$$

(29a)

$$\begin{array}{c}{C}_{b,w}={C}_{b,{w}_{2}}\\ {k}_{b,p}={C}_{b,{w}_{2}}/2\end{array}\}\phantom{\rule{0.16667em}{0ex}}\text{otherwise}.$$

(29b)

The arteriole compliance is then described by the sigmoid curve

$${C}_{b,a}=\left({C}_{b,a,0}-{C}_{b,w}+({C}_{b,a,0}+{C}_{b,w})exp\left(\frac{-{x}_{b,\mathit{aut}}}{{k}_{b,p}}\right)\right)/\left(1+exp\left(\frac{-{x}_{b,\mathit{aut}}}{{k}_{b,p}}\right)\right),$$

(30)

where exp(·) is the exponential function. Thus if from steady state, the cerebral arteriole flow *q _{b}*

The arteriole compliance can now be incorporated into the arteriole state equations. The equations for cerebral arteriole volume *V _{b}*

$$\frac{d{P}_{b,a}}{dt}=\frac{{q}_{b,pa}-{q}_{b,{f}_{1}}-{q}_{b,a}}{{C}_{b,a}}-\frac{{V}_{b,a}}{{C}_{b,a}^{2}}\frac{d{C}_{b,a}}{dt}+\frac{d{P}_{b,ic}}{dt},$$

(31)

where *q _{b}*

$$\frac{d{C}_{b,a}}{dt}=-\frac{2{C}_{b,w}}{{k}_{b,p}}exp\left(\frac{-{x}_{b,\mathit{aut}}}{{k}_{b,p}}\right){\left(1+exp\left(\frac{-{x}_{b,\mathit{aut}}}{{k}_{b,p}}\right)\right)}^{-2}\frac{d{x}_{b,\mathit{aut}}}{dt}+\frac{d{C}_{b,a,\mathit{input}}}{dt}.$$

(32)

The compliance input
$\frac{d{C}_{b,a,\mathit{input}}}{dt}$ 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}*

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

$${q}_{b,{f}_{1}}=\frac{{P}_{b,a}-{P}_{b,ic}}{{R}_{b,{f}_{1},0}}$$

(33)

and the absorption rate is unidirectional

$${q}_{b,{f}_{2}}=\frac{{P}_{b,ic}-{P}_{b,nv}}{{R}_{b,{f}_{2},0}}\phantom{\rule{0.16667em}{0ex}}\text{when}\phantom{\rule{0.16667em}{0ex}}{P}_{b,ic}>{P}_{b,nv}\phantom{\rule{0.16667em}{0ex}}\text{and}$$

(34a)

$${q}_{b,{f}_{2}}=0\phantom{\rule{0.16667em}{0ex}}\text{otherwise.}$$

(34b)

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

$$\frac{d{V}_{b,\mathit{CSF}}}{dt}={q}_{b,{f}_{1}}-{q}_{b,{f}_{2}}.$$

(35)

Now that all of the time-varying intracranial fluid volumes in the model are described, the intracranial pressure *P _{b}*

$${P}_{b,ic}=\frac{{V}_{b,ic}}{{C}_{b,ic}},$$

(36)

where the intracranial volume is the sum of all intracranial volumes

$${V}_{b,ic}={V}_{b,pa}+{V}_{b,a}+{V}_{b,c}+{V}_{b,v}+{V}_{b,pv}+{V}_{b,\mathit{CSF}}+{V}_{b,\mathit{ISF}}+{V}_{b,\mathit{ICF}}+{V}_{b,\mathit{tissue}},$$

(37)

and the interstitial fluid *V _{b}*

$$\frac{d{V}_{b,ic}}{dt}={q}_{b,na}-{q}_{b,pv}-{q}_{b,{f}_{2}},$$

(38)

where the inflow through the neck artery *q _{b}*

The remaining unknown is the intracranial compliance *C _{b}*

$${C}_{b,ic}={C}_{b,ic,0}{\left(\frac{{P}_{b,ic,0}}{{P}_{b,ic}}\right)}^{{\scriptstyle \frac{{\beta}_{b,ic}-1}{{\beta}_{b,ic}}}}.$$

(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}*

The systemic gas exchange model tracks the bound and unbound O_{2} and CO_{2} 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}*

$${C}_{{O}_{2}}={C}_{\mathit{bound},{O}_{2}}+{C}_{un\mathit{bound},{O}_{2}}.$$

(40)

The bound O_{2} concentration depends on the percent oxygen saturation *S _{O}*

$${C}_{\mathit{bound},{O}_{2}}=\frac{{S}_{{O}_{2}}}{100}\frac{4R{T}_{\mathit{STP}}}{{P}_{\mathit{STP}}}{M}_{\mathit{HbT},\mathit{blood}}.$$

(41)

The percent oxygen saturation *S _{O}*

$${C}_{un\mathit{bound},{O}_{2}}={P}_{{O}_{2}}{\alpha}_{{O}_{2}}.$$

(42)

In the derivation of state equations that follows, it will be necessary to convert oxygen concentration *C _{O}*

$$\frac{d{C}_{{O}_{2}}}{dt}=\frac{d{C}_{{O}_{2}}}{d{P}_{{O}_{2}}}\frac{d{P}_{{O}_{2}}}{dt},$$

(43)

where

$$\frac{d{C}_{{O}_{2}}}{d{P}_{{O}_{2}}}=\frac{d{S}_{{O}_{2}}}{d{P}_{{O}_{2}}}\frac{4R{T}_{\mathit{STP}}}{100{P}_{\mathit{STP}}}{M}_{\mathit{HbT},\mathit{blood}}+{\alpha}_{{O}_{2}}$$

(44)

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

$${\stackrel{\sim}{\alpha}}_{{O}_{2}}=\frac{d{C}_{{O}_{2}}}{d{P}_{{O}_{2}}}.$$

(45)

The general expression for total CO_{2} concentration in blood is similar to the expression for oxygen in that both bound and unbound gas must be considered

$${C}_{C{O}_{2}}={C}_{\mathit{bound},C{O}_{2}}+{C}_{un\mathit{bound},C{O}_{2}}.$$

(46)

The bound and unbound CO_{2} concentrations depend on the partial pressure of carbon dioxide *P _{CO}*

$${\stackrel{\sim}{\alpha}}_{C{O}_{2}}=\frac{d{C}_{C{O}_{2}}}{d{P}_{C{O}_{2}}}.$$

(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_{2}, CO_{2} to the systemic arterioles *P _{s}*

$${C}_{s,g,c}=\frac{{V}_{s,g,c}}{{V}_{s,c}},$$

(48)

where *V _{s}*

$$\frac{d{C}_{s,g,c}}{dt}=\frac{d{V}_{s,g,c}}{dt}\frac{1}{{V}_{s,c}}-\frac{{V}_{s,g,c}}{{V}_{s,c}^{2}}\frac{d{V}_{s,c}}{dt}$$

(49)

where the derivative of gas volume *V _{s}*

$$\frac{d{V}_{s,g,c}}{dt}={C}_{s,g,a}{q}_{s,a}-{C}_{s,g,c}{q}_{s,c}+{D}_{s,L,g}({P}_{s,g,\mathit{ISF}}-{P}_{s,g,c}),$$

(50)

where *D _{s}*

$$\frac{d{P}_{s,g,c}}{dt}=\frac{1}{{V}_{s,c}{\stackrel{\sim}{\alpha}}_{s,g,c}}\left({C}_{s,g,a}{q}_{s,a}-{C}_{s,g,c}{q}_{s,c}-{C}_{s,g,c}\frac{d{V}_{s,c}}{dt}+{D}_{s,L,g}({P}_{s,g,\mathit{ISF}}-{P}_{s,g,c})\right).$$

(51)

A similar state equation governs the venule *v* chamber

$$\frac{d{P}_{s,g,v}}{dt}=\frac{1}{{V}_{s,v}{\stackrel{\sim}{\alpha}}_{s,g,v}}\left({C}_{s,g,c}{q}_{s,c}-{C}_{s,g,v}{q}_{s,v}-{C}_{s,g,v}\frac{d{V}_{s,v}}{dt}\right).$$

(52)

The state equation in the intracellular chamber includes the effects of metabolism *Q _{s}*

$$\frac{d{P}_{s,g,\mathit{ICF}}}{dt}=\frac{1}{{V}_{s,\mathit{ICF}}{\alpha}_{g}}({D}_{s,M,g}({P}_{s,g,\mathit{ISF}}-{P}_{s,g,\mathit{ICF}})-{Q}_{s,g}).$$

(53)

A similar expression is use for the interstitial compartment

$$\frac{d{P}_{s,g,\mathit{ISF}}}{dt}=\frac{1}{{V}_{s,\mathit{ISF}}{\alpha}_{g}}({D}_{s,M,g}({P}_{s,g,\mathit{ICF}}-{P}_{s,g,\mathit{ISF}})+{D}_{s,L,g}({P}_{s,g,c}-{P}_{s,g,\mathit{ISF}})).$$

(54)

Following the same procedure as with the systemic compartments, preliminary calculations are carried out to obtain the cerebral blood gas concentration *C _{b}*

$$\frac{d{Q}_{b,g}}{dt}=\frac{d{Q}_{b,g,\mathit{input}}}{dt}.$$

(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}*

$$\frac{d{P}_{b,g,a}}{dt}=\frac{1}{{V}_{b,a}{\stackrel{\sim}{\alpha}}_{b,g,a}}\left({C}_{b,g,a}{q}_{b,pa}-{C}_{b,g,a}{q}_{b,a}-{q}_{b,{f}_{1}}{P}_{b,g,a}{\alpha}_{g}-{C}_{b,g,a}\frac{d{V}_{b,a}}{dt}\right),$$

(56)

capillary

$$\frac{d{P}_{b,g,c}}{dt}=\frac{1}{{V}_{b,c}{\stackrel{\sim}{\alpha}}_{b,g,c}}\left({C}_{b,g,a}{q}_{b,a}-{C}_{b,g,c}{q}_{b,c}-{C}_{b,g,c}\frac{d{V}_{b,c}}{dt}+{D}_{b,L,g}({P}_{b,g,\mathit{ISF}}-{P}_{b,g,c})\right),$$

(57)

venule

$$\frac{d{P}_{b,g,v}}{dt}=\frac{1}{{V}_{b,v}{\stackrel{\sim}{\alpha}}_{b,g,v}}\left({C}_{b,g,c}{q}_{b,c}-{C}_{b,g,v}{q}_{b,v}-{C}_{b,g,v}\frac{d{V}_{b,v}}{dt}\right),$$

(58)

cerebrospinal fluid

$$\frac{{dP}_{b,g,\mathit{CSF}}}{dt}=\frac{1}{{V}_{b,\mathit{CSF}}{\alpha}_{g}}({\alpha}_{g}({q}_{b,{f}_{1}}{P}_{b,g,a}-{q}_{b,{f}_{2}}{P}_{b,g,\mathit{CSF}})+{D}_{b,N,g}({P}_{b,g,\mathit{ISF}}-{P}_{b,g,\mathit{CSF}})),$$

(59)

interstitial

$$\begin{array}{l}\frac{d{P}_{b,g,\mathit{ISF}}}{dt}=\frac{1}{{V}_{b,\mathit{ISF}}{\alpha}_{g}}({D}_{b,L,g}({P}_{b,g,c}-{P}_{b,g,\mathit{ISF}})\cdots \\ +{D}_{b,M,g}({P}_{b,g,\mathit{ICF}}-{P}_{b,g,\mathit{ISF}})+{D}_{b,N,g}({P}_{b,g,\mathit{CSF}}-{P}_{b,g,\mathit{ISF}})),\end{array}$$

(60)

and intracellular compartments

$$\frac{d{P}_{b,g,\mathit{ICF}}}{dt}=\frac{1}{{V}_{b,\mathit{ICF}}{\alpha}_{g}}({D}_{b,M,g}({P}_{b,g,\mathit{ISF}}-{P}_{b,g,\mathit{ICF}})-{Q}_{b,g}).$$

(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.

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}*

$${M}_{r,\mathit{HbO},j}=\frac{{S}_{r,{O}_{2},j}}{100}{M}_{\mathit{HbT},\mathit{blood}}$$

(62)

and deoxyhemoglobin

$${M}_{r,\mathit{HbR},j}=\left(1-\frac{{S}_{r,{O}_{2},j}}{100}\right){M}_{\mathit{HbT},\mathit{blood}},$$

(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 *j*th chamber in the systemic or cerebral region (*r* = *s*, *b*) is the state variable *V _{r}*

$${\stackrel{\sim}{M}}_{r,h,j}=\frac{{V}_{r,j}}{{V}_{r,j,0}}{M}_{r,h,j},$$

(64)

where *M _{r}*

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}*

$${\stackrel{\sim}{M}}_{r,h}={F}_{r,bv}\sum _{j}{F}_{r,j}{\stackrel{\sim}{M}}_{r,h,j},$$

(65)

such that _{r}_{,}* _{h}* approximates the time-varying distributed concentration of hemoglobin species

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}*

In the second step, the calculated pathlengths *L _{r}*

$$\mathrm{\Delta}O{D}_{\lambda}=\mathrm{\Delta}{\mu}_{a,s,\lambda}{L}_{s,\lambda}+\mathrm{\Delta}{\mu}_{a,b,\lambda}{L}_{b,\lambda},$$

(66)

where *λ* is the wavelength of the NIRS measurement. Next the absorption coefficient perturbations Δ*μ _{a}*

$$\mathrm{\Delta}{\mu}_{a,r,\lambda}=\mathrm{\Delta}{\stackrel{\sim}{M}}_{r,\mathit{HbO}}{\epsilon}_{\mathit{HbO},\lambda}+\mathrm{\Delta}{\stackrel{\sim}{M}}_{r,\mathit{HbR}}{\epsilon}_{\mathit{HbR},\lambda},$$

(67)

where *ε _{HbO}*

$$\mathrm{\Delta}{\stackrel{\sim}{M}}_{r,h}={\stackrel{\sim}{M}}_{r,h}-{\stackrel{\sim}{M}}_{r,h,0},$$

(68)

where _{r}_{,}* _{h}* is obtained from Eq. (65) and

The SimCVR model dynamics are first illustrated in simulation by showing how the model responds to step increases in ABP and arterial blood CO_{2}. 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.

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).

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}*

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).

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_{2} 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_{2} flows into the cerebral capillary chamber and diffuses into the ISF [Eq. (60)]. The autoregulation blood flow set point *q _{b}*

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

$$\mathrm{\Delta}O{D}_{\lambda ,\mathit{data}}=-log\left(\frac{{I}_{\lambda}}{{I}_{\lambda ,0}}\right),$$

(69)

where *I _{λ}*

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 $\frac{d{P}_{s,sa,\mathit{input}}}{dt}$ 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}*

$$f=\frac{\sum {\left(\mathrm{\Delta}O{D}_{{\lambda}_{1},\mathit{model}}-\mathrm{\Delta}O{D}_{{\lambda}_{1},\mathit{data}}\right)}^{2}}{\sum {\left(\mathrm{\Delta}O{D}_{{\lambda}_{1},\mathit{data}}\right)}^{2}}+\frac{\sum {\left(\mathrm{\Delta}O{D}_{{\lambda}_{2},\mathit{model}}-\mathrm{\Delta}O{D}_{{\lambda}_{2},\mathit{data}}\right)}^{2}}{\sum {\left(\mathrm{\Delta}O{D}_{{\lambda}_{2},\mathit{data}}\right)}^{2}}.$$

(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 *λ*_{1} = 890 nm and *λ*_{2} = 690 nm. The summations in Eq. (70) are over time and the volume fractions *F _{r}*

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

$$0.01<{F}_{s,bv}<0.06,$$

(71)

and the brain blood volume fraction

$$0.02<{F}_{b,bv}<\mathrm{0.12.}$$

(72)

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

$$0.01<{F}_{s,j}<0.10,$$

(73)

and in the smaller vessels

$$0.05<{F}_{r,j}<\mathrm{0.50.}$$

(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}*

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.

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^{2} 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^{2} 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^{2} 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^{2} values were found in 9 out of 10 subjects for the 830 and 690 nm NIRS measurements. The higher R^{2} 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.

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^{2} values from the 5-minute runs are reported **...**

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

Next we explored the possibility that the improved R^{2} 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

$$G(s)={K}_{p}\frac{1+{\tau}_{z}s}{1+2\zeta {\tau}_{w}s+{\tau}_{w}^{2}{s}^{2}},$$

(75)

where *s* is the Laplace variable, *K _{p}* is the static gain,

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

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^{2} 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.

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^{2} 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.

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^{2} 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].

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.

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^{2} 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^{2} 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.

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.

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.

The Kelman equation for oxygemoglobin dissociation uses the coefficients

$$\begin{array}{l}{a}_{1}=-8532.2289\\ {a}_{2}=2121.4010\\ {a}_{3}=-67.073989\\ {a}_{4}=935960.87\\ {a}_{5}=-31346.258\\ {a}_{6}=2396.1674\\ {a}_{7}=-\mathrm{67.104406.}\end{array}$$

(76)

The calculation starts with a scale factor *v*, which accounts for changes in temperature, pH, and the partial pressure of CO_{2}

$$v={10}^{[0.024(37-T)+0.40(pH-7.4)+0.06({log}_{10}(40)-{log}_{10}({P}_{C{O}_{2}})\left)\right]}.$$

(77)

Next the following empirical relationship is used to compute hemoglobin saturation

$${S}_{{O}_{2}}=100\frac{{a}_{1}{P}_{{O}_{2}}v+{a}_{2}{P}_{{O}_{2}}^{2}{v}^{2}+{a}_{3}{P}_{{O}_{2}}^{3}{v}^{3}+{P}_{{O}_{2}}^{4}{v}^{4}}{{a}_{4}+{a}_{5}{P}_{{O}_{2}}v+{a}_{6}{P}_{{O}_{2}}^{2}{v}^{2}+{a}_{7}{P}_{{O}_{2}}^{3}{v}^{3}+{P}_{{O}_{2}}^{4}{v}^{4}}.$$

(78)

We have further differentiated Kelman’s equation with respect to the partial pressure of oxygen as needed for the effective solubility _{O}_{2} in Eq. (45)

$$\begin{array}{l}\frac{d{S}_{{O}_{2}}}{d{P}_{{O}_{2}}}=100\frac{{a}_{1}v+2{a}_{2}{P}_{{O}_{2}}{v}^{2}+3{a}_{3}{P}_{{O}_{2}}^{2}{v}^{3}+4{P}_{{O}_{2}}^{3}{v}^{4}}{{a}_{4}+{a}_{5}{P}_{{O}_{2}}v+{a}_{6}{P}_{{O}_{2}}^{2}{v}^{2}+{a}_{7}{P}_{{O}_{2}}^{3}{v}^{3}+{P}_{{O}_{2}}^{4}{v}^{4}}\cdots \\ -100\frac{({a}_{1}{P}_{{O}_{2}}v+{a}_{2}{P}_{{O}_{2}}^{2}{v}^{2}+{a}_{3}{P}_{{O}_{2}}^{3}{v}^{3}+{P}_{{O}_{2}}^{4}{v}^{4})({a}_{5}v+2{a}_{6}{P}_{{O}_{2}}{v}^{2}+3{a}_{7}{P}_{{O}_{2}}^{2}{v}^{3}+4{P}_{{O}_{2}}^{3}{v}^{4})}{{({a}_{4}+{a}_{5}{P}_{{O}_{2}}v+{a}_{6}{P}_{{O}_{2}}^{2}{v}^{2}+{a}_{7}{P}_{{O}_{2}}^{3}{v}^{3}+{P}_{{O}_{2}}^{4}{v}^{4})}^{2}}.\end{array}$$

(79)

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

$$pK=6.086+0.042(7.4-pH)+(38-T)(0.00472+0.00139(7.4-pH)).$$

(80)

Now the concentration ratio of oxyhemoglobin-bound CO_{2} to plasma CO_{2} is estimated as

$$\mathit{Dox}=0.590+0.2913(7.4-pH)-0.0844{(7.4-pH)}^{2},$$

(81)

and the concentration ratio of deoxyhemoglobin-bound CO_{2} to plasma CO_{2} is

$$Dr=0.664+0.2275(7.4-pH)-0.0938{(7.4-pH)}^{2}.$$

(82)

The unbound CO_{2} concentration is computed with the Henderson-Hasselbalch equation

$${C}_{un\mathit{bound},C{O}_{2}}=(1-\mathit{Hct}){\alpha}_{C{O}_{2}}{P}_{C{O}_{2}}(1+{10}^{pH-pK}),$$

(83)

and the bound CO_{2} concentration accounts for hemoglobin binding with *Dox* and *Dr*

$${C}_{\mathit{bound},C{O}_{2}}={C}_{un\mathit{bound},C{O}_{2}}\frac{\mathit{Hct}}{1-\mathit{Hct}}\left(\mathit{Dox}\frac{{S}_{{O}_{2}}}{100}+Dr\left(1-\frac{{S}_{{O}_{2}}}{100}\right)\right).$$

(84)

Finally, we will take the derivative with respect to the partial pressure of CO_{2}

$$\frac{d{C}_{C{O}_{2}}}{d{P}_{C{O}_{2}}}=(1-\mathit{Hct}){\alpha}_{C{O}_{2}}(1+{10}^{pH-pK})\left(1+\frac{\mathit{Hct}}{1-\mathit{Hct}}\left(\mathit{Dox}\frac{{S}_{{O}_{2}}}{100}+Dr\left(1-\frac{{S}_{{O}_{2}}}{100}\right)\right)\right),$$

(85)

which is the effective solubility of carbon dioxide in the blood _{CO}_{2} in Eq. (47).

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_{2} and negative for CO_{2}.

The initial partial pressures of O_{2} and CO_{2} in the systemic capillary are computed with a nonlinear search that minimizes the objective function

$$f={({q}_{s,0}({C}_{s,{O}_{2},a,0}-{C}_{s,{O}_{2},c,0})-{Q}_{s,{O}_{2},0})}^{2}+{({q}_{s,0}({C}_{s,C{O}_{2},a,0}-{C}_{s,C{O}_{2},c,0})-{Q}_{s,C{O}_{2},0})}^{2},$$

(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

$${P}_{s,g,v}={P}_{s,g,c},$$

(87)

with *g* = O_{2}, CO_{2}. The partial pressure of O_{2} should monotonically decrease from artery to capillary to ISF to ICF, whereas the partial pressure of CO_{2} should monotonically increase along the same cascade. The diffusion constants can then be calculated as

$${D}_{s,L,g}={Q}_{s,g,0}/({P}_{s,g,c,0}-{P}_{s,g,\mathit{ISF},0}),$$

(88)

$${D}_{s,M,g}={Q}_{s,g,0}/({P}_{s,g,\mathit{ISF},0}-{P}_{s,g,\mathit{ICF},0}).$$

(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

$${C}_{s,j,0}={V}_{s,j,0}(1-{F}_{s,us})/{P}_{s,j,0},$$

(90)

$${R}_{s,j,0}=({P}_{s,j,0}-{P}_{s,k,0})/{q}_{s,0},$$

(91)

where *F _{s}*

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

$${C}_{b,na,0}={V}_{b,na,0}(1-{F}_{b,us})/{P}_{s,la,0},$$

(92)

$${R}_{b,na,0}=({P}_{s,la,0}-{P}_{b,pa,0})/({q}_{b,0}+{q}_{b,f,0}),$$

(93)

$${C}_{b,nv,0}={V}_{b,nv,0}(1-{F}_{b,us})/{P}_{b,nv,0},$$

(94)

$${R}_{b,nv,0}=({P}_{b,nv,0}-{P}_{s,lv,0})/({q}_{b,0}+{q}_{b,f,0}).$$

(95)

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

$${C}_{b,j,0}={V}_{b,j,0}(1-{F}_{b,us})/({P}_{b,j,0}-{P}_{b,ic,0}),$$

(96)

$${R}_{b,j,0}=({P}_{b,j,0}-{P}_{b,k,0})/({q}_{b,0}+{q}_{b,f,0}),$$

(97)

and for the Windkessel compartments, *j* = *a*, *c*, *v* and *k* = *c*, *v*, *pv*, the unstressed volume fraction *F _{b}*

$$\begin{array}{l}{V}_{b,ic,0}={V}_{b,pa,0}+{V}_{b,a,0}+{V}_{b,c,0}+{V}_{b,v,0}+{V}_{b,pv,0}\cdots \\ +{V}_{b,\mathit{CSF},0}+{V}_{b,\mathit{ISF},0}+{V}_{b,\mathit{ICF},0}+{V}_{b,\mathit{tissue},0}\end{array}$$

(98)

with the compliance

$${C}_{b,ic,0}={V}_{b,ic,0}/{P}_{b,ic,0}.$$

(99)

The resistances for the CSF flow are

$${R}_{b,{f}_{1},0}=({P}_{b,a,0}-{P}_{b,ic,0})/{q}_{b,f,0},$$

(100)

$${R}_{b,{f}_{2},0}=({P}_{b,ic,0}-{P}_{b,nv,0})/{q}_{b,f,0}.$$

(101)

For the cerebral gas exchange, there is again only one capillary chamber therefore with *g* = O_{2}, CO_{2},

$${P}_{b,g,v}={P}_{b,g,c}.$$

(102)

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

$$\begin{array}{l}f={({C}_{s,{O}_{2},a,0}({q}_{b,0}+{q}_{b,f,0})-{C}_{b,{O}_{2},a,0}{q}_{b,0}-{\alpha}_{{O}_{2}}{P}_{b,{O}_{2},a,0}{q}_{b,f,0})}^{2}\cdots \\ +{({C}_{s,C{O}_{2},a,0}({q}_{b,0}+{q}_{b,f,0})-{C}_{b,C{O}_{2},a,0}{q}_{b,0}-{\alpha}_{C{O}_{2}}{P}_{b,C{O}_{2},a,0}{q}_{b,f,0})}^{2},\end{array}$$

(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

$${Q}_{b,g,0}={C}_{s,g,a,0}({q}_{b,0}+{q}_{b,f,0})-{C}_{b,g,c,0}{q}_{b,0}-{\alpha}_{g}{P}_{b,g,\mathit{CSF},0}{q}_{b,f,0}$$

(104)

and the diffusion constants are

$${D}_{b,L,g}={q}_{b,0}({C}_{b,g,a,0}-{C}_{b,g,c,0})/({P}_{b,g,c,0}-{P}_{b,g,\mathit{ISF},0}),$$

(105)

$${D}_{b,M,g}={Q}_{b,g,0}/({P}_{b,g,\mathit{ISF},0}-{P}_{b,g,\mathit{ICF},0}),$$

(106)

$${D}_{b,N,g}={\alpha}_{g}{q}_{b,f,0}({P}_{b,g,a,0}-{P}_{b,g,\mathit{CSF},0})/({P}_{b,g,\mathit{CSF},0}-{P}_{b,g,\mathit{ISF},0}).$$

(107)

Finally, the autoregulation state variable *x _{b}*

**Publisher's Disclaimer: **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.

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. [PMC free article] [PubMed]

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |