|Home | About | Journals | Submit | Contact Us | Français|
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 , 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  that can be used to grade autoregulatory function from noninvasively measured blood pressure and cerebral circulation  with applications in, for example, hypertension, stroke, and carotid sinus syndrome . 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 .
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  and Windkessel models , 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 . More recent studies offer sophisticated models for the relationship of cerebral blood flow to metabolic biochemistry and vascular smooth muscle  and the effects of cerebral autoregulation on neural activation . 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)  and in fMRI , 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) , autoregressive models , adaptive filters , and with Kalman filter in NIRS . 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 . CA assessment with TCD has been explored in hypertension , severe carotid stenosis , vertebral artery disease , critically ill infants  and acute ischemic stroke . 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  and a model of systemic circulation . Features have been added that account for the affinity of CO2 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  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 R2 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 O2 and CO2 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
Eq. (1) specifies the compliance constant Cs, j,0 in the systemic s, vascular chamber j, as the ratio of volume Vs, j to transmural pressure Ps, j, which is represented in this electrical analog as a capacitor (Fig. 1A). The unstressed volume Vs, 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 Fs,us of the initial volume
Following , the compliance of the collapsible veins is increased 20 fold when the transmural pressure Ps, j is negative. Eq. (2) follows from Poiseuille’s law for laminar flow in a tube in which resistance Rs, j is inversely proportional to vessel diameter to the forth power. Eq. (3) states that blood flow qs, j out of chamber j equals the pressure drop from the current chamber j to the subsequent chamber k divided by the resistance Rs, j.
Calculating the system response with an ODE solver requires state equations that relate the time derivatives of the pressure and volume to the pressure Ps, j and volume Vs, j state variables. The state equations for i = la, v, sv and j = sa, sv, lv are obtained by solving Eq. (1) for Ps, j and then taking the time derivative
where the change in volume is the difference between blood outflow qs, j and inflow qs,i
The systemic arteriole, capillary and venule compartments include Windkessel compliance to appropriately model nonlinear vascular properties . Using the notational substitutions i = sa, a, c; j = a, c, v; and k = c, v, sv, the compliance is
where βs, j is the vascular reserve of the jth compartment in the systemic s region and
is the pressure-volume relationship with time-varying compliance Cs, 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)
Solving Eq. (9) for the derivative of pressure yields
The state equation for chamber la is driven by the derivative of the blood pressure in the large artery
because the input is directly connected to the large artery without resistance in the circuit model (Fig. 1A). The time-varying pressure input Ps,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 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 Ps,la from a known value of Ps,sa. First Eqs. (1), (2) and (3) are combined with j = la and k = sa to obtain
where unknown value of Ps,la is only real positive root of Eq. (15).
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
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
and is derived in the same manner as Eqs. (5) and (6) where qb, f2 is CSF absorption into the neck vein (see Fig. 1A). The equations for resistance Rb, j and flow qb, 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 Cb, j,0 of the pial arteries pa and pial veins pv, the volume for j = pa, pv is
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 Cb, j and account for the intracranial pressure
where j = c, v and
The expressions for resistance Rb, j, flow qb, j and the volume derivatives 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 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 qb,sp and autoregulation gain Kaut on this feedback loop are influenced by the partial pressure of CO2 in the ISF. The following system equations that govern the CO2 effects are simplified versions of those presented in  because data simulations were found to operate in a linear regime of the more complex equations. The difference from baseline ISF CO2 is expressed as an autoregulation set point fraction
which correlates positively with the CBF set point
and negatively correlates with the autoregulation reactivity gain
The autoregulation state variable xb,aut is an auxiliary variable that will determine the arteriole compliance changes a few steps later. The dynamics of xb,aut are governed by its state equation
The variable cerebral arteriole compliance operates on a non-symmetric sigmoid with the bounds Cb,w and curvatures kb, p
The arteriole compliance is then described by the sigmoid curve
where exp(·) is the exponential function. Thus if from steady state, the cerebral arteriole flow qb,a falls below the CBF set point qb,sp, the autoregulation state xb,aut will become negative. This will cause the arteriole compliance Cb,a to increase towards Cb,a,0 +Cb,w1, which will dilate the arterioles and increase CBF. In the opposite scenario, Cb,a decreases towards Cb,a,0 − Cb,w1, 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 Vb,a, resistance Rb,a and flow qb,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 Pb,a in Eq. (19) with j = a, and then taking the time derivative to obtain
where qb, f1 is CSF formation from the arteriole (see Fig. 1A), and the compliance derivative is
The compliance input has been added to Eq. (32) to permit inclusion of an arteriolar compliance response to neural stimulus , which would be superimposed on the regulatory physiology. When this input is non-zero, the arteriole compliance Cb,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 . Following , the formation rate depends on the pressure gradient
and the absorption rate is unidirectional
The state equation for CSF accounts for the difference between the rates of formation and absorption
Now that all of the time-varying intracranial fluid volumes in the model are described, the intracranial pressure Pb,ic can be expressed in terms of the intracranial compliance Cb,ic and volume Vb,ic as
where the intracranial volume is the sum of all intracranial volumes
and the interstitial fluid Vb,ISF, intracellular fluid Vb,ICF and tissue volume Vb,tissue in the intracranial cavity are assumed to be known constants. Given the definition of Vb,ic in Eq. (37), the graphical representation of intracranial compliance Cb,ic as a single compartment in Fig. 1A is only symbolic. The state equation for intracranial volume Vb,ic is
where the inflow through the neck artery qb,na, outflow through the neck vein qb,nv, and CSF absorption qb, f2 are the only net effects on the volumes that comprise the intracranial volume Vb,ic in Eq. (37).
The remaining unknown is the intracranial compliance Cb,ic. In , 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  but the algebraic cancelation of Pb,ic is avoided
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.
The systemic gas exchange model tracks the bound and unbound O2 and CO2 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 CO2 in blood is the sum of bound and unbound O2
The bound O2 concentration depends on the percent oxygen saturation SO2, the hemoglobin concentration MHbT,blood, gas constant R, standard temperature TSTP and pressure PSTP
The percent oxygen saturation SO2 can be calculated from the partial pressure of oxygen PO2 using an analytic hemoglobin dissociation curve which depends on temperature T and pH  (included for convenience in Appendix Blood oxygen). The unbound O2 concentration simply depends on the partial pressure of oxygen PO2 and the gas solubility αO2
In the derivation of state equations that follows, it will be necessary to convert oxygen concentration CO2 state equations into oxygen partial pressure PO2 state equations. This conversion can be accomplished using the relation
The general expression for total CO2 concentration in blood is similar to the expression for oxygen in that both bound and unbound gas must be considered
The bound and unbound CO2 concentrations depend on the partial pressure of carbon dioxide PCO2, percent oxygen saturation SO2, hematocrit Hct, gas solubility αCO2, temperature T and pH  (see Blood carbon dioxide in Appendix). Also similar to the treatment of oxygen is the need for an effective solubility
The state equations for the systemic gas exchange can now be described beginning with the input to the partial pressure of gaseous species g = O2, CO2 to the systemic arterioles Ps,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 Qs,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 Cs,g,c for g = O2, CO2
where Vs,g,c is the gas volume and Vs,c is the blood volume. Differentiating Eq. (48) yields
where the derivative of gas volume Vs,g,c depends on the gas inflow from the arteriole, outflow from the capillary and diffusion with ISF
where Ds,L,g is the diffusion constant, Ps,g,c is the partial pressure in the capillary and Ps,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
A similar state equation governs the venule v chamber
The state equation in the intracellular chamber includes the effects of metabolism Qs,g and uses the solubility for unbound gas αg only because no hemoglobin is present
A similar expression is use for the interstitial compartment
Following the same procedure as with the systemic compartments, preliminary calculations are carried out to obtain the cerebral blood gas concentration Cb, g, j and effective solubility b, g, j. User inputs of the cerebral metabolic rates can model changing metabolic demands
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 qb, f1 and absorption qb, f2. The following state equations result for the arteriole
and intracellular compartments
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 MHbT, blood is constant. This assumption simplifies the conversion of time-varying blood oxygen saturation Sr,O2, j in the jth vascular chamber to time-varying concentration of oxyhemoglobin
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 Vr, j divided by the initial volume Vr, j,0. The effective hemoglobin concentration r,h, j for h = HbO, HbR is then
where Mr,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 Fr,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 Fr, 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
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 Vr, j and blood oxygen saturation Sr,O2, j as computed by the lumped parameter components of the SimCVR model.
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 Lr,λ 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 . 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 Lr,λ are already published [44; 45; 46].
In the second step, the calculated pathlengths Lr,λ are used in a first-order Rytov expansion  to model changes in optical density (ΔOD) that result from perturbations in the absorption coefficient Δμa,r,λ
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
where εHbO,λ and εHbR,λ are the tabulated extinction coefficients . Finally, the deviations in hemoglobin concentration are measured from an initial value
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.
The SimCVR model dynamics are first illustrated in simulation by showing how the model responds to step increases in ABP and arterial blood CO2. 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).
The following narrative qualitative describes the step response dynamics: The step increase in pressure in the systemic large artery Ps,la propagates through the neck arteries to the cerebral circulation where it causes an increase in cerebral arteriole pressure Pb,a and increased arteriole blood flow qb,a. When qb,a raises above the blood flow set point qb,sp, the autoregulation state xb,aut becomes positive [Eq. (28)]. The autoregulation state acts to decrease cerebral arteriole compliance Cb,a [Eq. (30)] because positive values of xb,aut pull Cb,a towards the lower bound [Eq. (29b)] of the sigmoid where Cb,a =Cb,a,0 − Cb,w. The effects of decreasing Cb,a are that cerebral arteriole volume Vb,a goes down [Eq. (18)] and arteriole resistance Rb,a goes up [Eq. (2)]. The increased Rb,a counteracts the higher Pb,a and lowers the arteriole flow qb,a [Eq. (3)] towards the set point qb,sp. This return of cerebral blood flow during a step increase in arterial pressure is the expected autoregulatory response.
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).
Relative to the pressure step response, the CO2 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 CO2 flows into the cerebral capillary chamber and diffuses into the ISF [Eq. (60)]. The autoregulation blood flow set point qb,sp increases with increased CO2 in the ISF Pb,CO2,ISF [Eqs. (25) and (26)]. Meanwhile, the autoregulation gain Kaut decreases [Eq. (27)]. Since the simulation started at steady state with all time derivatives equal to zero, the increased qb,sp results in a negative derivative of the autoregulation state xb,aut [Eq. (28)], which pulls its value negative. The lower Kaut just makes this happen more gradually. The autoregulation state acts on the cerebral arteriole compliance Cb,a [Eq. (30)] such that negative values of xb,aut raise Cb,a towards the upper bound [Eq. (29a)] of the sigmoid where Cb,a =Cb,a,0 + Cb,w. The effects of increasing Cb,a are that cerebral arteriole volume Vb,a goes up [Eq. (18)] and arteriole resistance Rb,a goes down [form of Eq. (2)]. Increasing Vb,a drives up the intracranial volume Vb,ic [Eq. (37)] and increases intracranial pressure Pb,ic [form of Eqs. (10) and (11)]. Considering the arteriole flow qb,a [form of Eq. (3)], decreasing Pb,a will lower qb,a and decreasing Rb,a will raise qb,a. In our simulation, qb,a drops slightly for about 15 seconds before increasing, resulting in an initial net outflow from the capillary. Now examining the capillary pressure Pb,c [Eqs. (24) and (22)], the initial negative derivative of capillary volume Vb,c would decrease Pb,c, except that it is dominated by the positive derivative of intracranial pressure Pb,ic. The increase in Pb,c drives up the capillary blood flow qb,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 O2 is also expected due to the interaction of CO2 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 (CMRO2) is constant.
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 . 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  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
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 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 Fr,bv and vascular compositions Fr, j separately for each detector by a constrained nonlinear optimization (FMINCON function in MATLAB®, MathWorks, Inc., Natick, MA) that minimized the objective function
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 Fr,bv and Fr, j act as regression coefficients that determine the linear combination of time-varying blood oxygen Sr,O2, j and time-varying blood volume Vr, j that best predict the measured optical density data. The optimization is constrained, however, because blood volume fractions Fr,bv and vascular composition fractions Fr, 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  and 5.2 ± 1.4 % in the gray matter of the brain . The arterial volume fraction in the brain is approximately 30% , 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
and the brain blood volume fraction
The vascular composition bounds for the large arteries and large veins (j = la, lv) were
and in the smaller vessels
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 Fs,bv was 4% ± 1% and for the cerebral blood volume fraction Fb,bv it was 7% ± 3%. The results for the vascular composition Fr, 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.
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 R2 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 R2 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 R2 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 R2 values were found in 9 out of 10 subjects for the 830 and 690 nm NIRS measurements. The higher R2 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.
Next we explored the possibility that the improved R2 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 , we chose to fit a 2nd order continuous-time system transfer function G of the form
where s is the Laplace variable, Kp 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.
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 R2 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 R2 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 R2 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)  and possibly also a flow-sensitive MRI technique like arterial spin labeling (ASL) .
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  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 , the anatomy of the cerebral vessels , and neurovascular coupling . 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 R2 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 R2 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
The calculation starts with a scale factor v, which accounts for changes in temperature, pH, and the partial pressure of CO2
Next the following empirical relationship is used to compute hemoglobin saturation
We have further differentiated Kelman’s equation with respect to the partial pressure of oxygen as needed for the effective solubility O2 in Eq. (45)
First and effective pH is calculated to account for differences in temperature T
Now the concentration ratio of oxyhemoglobin-bound CO2 to plasma CO2 is estimated as
and the concentration ratio of deoxyhemoglobin-bound CO2 to plasma CO2 is
The unbound CO2 concentration is computed with the Henderson-Hasselbalch equation
and the bound CO2 concentration accounts for hemoglobin binding with Dox and Dr
Finally, we will take the derivative with respect to the partial pressure of CO2
which is the effective solubility of carbon dioxide in the blood CO2 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 O2 and negative for CO2.
The initial partial pressures of O2 and CO2 in the systemic capillary are computed with a nonlinear search that minimizes the objective function
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
with g = O2, CO2. The partial pressure of O2 should monotonically decrease from artery to capillary to ISF to ICF, whereas the partial pressure of CO2 should monotonically increase along the same cascade. The diffusion constants can then be calculated as
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
where Fs,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 Fs,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
For j = pa, pv and k = a, nv we have
and for the Windkessel compartments, j = a, c, v and k = c, v, pv, the unstressed volume fraction Fb,us is neglected. The initial intracranial volume is
with the compliance
The resistances for the CSF flow are
For the cerebral gas exchange, there is again only one capillary chamber therefore with g = O2, CO2,
The partial pressure of O2 should monotonically decrease from artery to capillary to CSF to ISF to ICF, whereas the partial pressure of CO2 should monotonically increase along the same cascade. The initial partial pressures of O2 and CO2 in the cerebral artery are found with a nonlinear search that minimizes the objective function
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
and the diffusion constants are
Finally, the autoregulation state variable xb,aut should be initialized at zero.
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.