SimCVR consists of coupled fluid flow and gas exchange systems. The fluid flow system is a circuit-analog of nonlinear resistors and capacitors () 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 () 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).

Systemic fluid flow system

The circulatory model assumes that the pressure-volume relationship in each chamber is governed by compliance and that the pressure drop between chambers and the vascular resistance determines the blood flow. Following the nomenclature defined in , 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

*C*_{s}_{,} _{j}_{,0} in the systemic

*s*, vascular chamber

*j*, as the ratio of volume

*V*_{s}_{,} _{j} to transmural pressure

*P*_{s}_{,} _{j,} which is represented in this electrical analog as a capacitor (). The unstressed volume

*V*_{s}_{,} _{j}_{,}_{us} is the residual volume in the vessel that is present when the transmural pressure falls to zero. The unstressed volume is assumed to be constant and has been initialized as a fraction

*F*_{s}_{,}_{us} of the initial volume

Following [

35], the compliance of the collapsible veins is increased 20 fold when the transmural pressure

*P*_{s}_{,} _{j} is negative.

Eq. (2) follows from Poiseuille’s law for laminar flow in a tube in which resistance

*R*_{s}_{,} _{j} is inversely proportional to vessel diameter to the forth power.

Eq. (3) states that blood flow

*q*_{s}_{,} _{j} out of chamber

*j* equals the pressure drop from the current chamber

*j* to the subsequent chamber

*k* divided by the resistance

*R*_{s}_{,} _{j.}Calculating the system response with an ODE solver requires state equations that relate the time derivatives of the pressure

and volume

to the pressure

*P*_{s}_{,} _{j} and volume

*V*_{s}_{,} _{j} state variables. The state equations for

*i* =

*la*,

*v*,

*sv* and

*j* =

*sa*,

*sv*,

*lv* are obtained by solving

Eq. (1) for

*P*_{s}_{,} _{j} and then taking the time derivative

where the change in volume is the difference between blood outflow

*q*_{s}_{,} _{j} and inflow

*q*_{s}_{,}_{i} 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

where

*β*_{s}_{,} _{j} is the vascular reserve of the

*j*th compartment in the systemic

*s* region and

is the pressure-volume relationship with time-varying compliance

*C*_{s}_{,} _{j.} The equations for resistance and flow are the same as

Eqs. (2) and

(3) respectively. The state equations for the Windkessel compartments are obtained from the time-derivative of

Eq. (8) Solving

Eq. (9) for the derivative of pressure yields

Specifying inputs to the systemic system

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

because the input is directly connected to the large artery without resistance in the circuit model (). The time-varying pressure input

*P*_{s}_{,}_{la}_{,}_{input} must be specified by the user. An alternative to using the large artery pressure input of

Eq. (12) is to supply the derivative of blood pressure in the small artery

*sa* instead of the large artery

*la*. The advantage of specifying

is that it can be used to drive the model with the pressure variations measured in the peripheral circulation (e.g. with a finger cuff). However, the following additional calculations are required to compute the unknown value of

*P*_{s}_{,}_{la} from a known value of

*P*_{s}_{,}_{sa.} First

Eqs. (1),

(2) and

(3) are combined with

*j* =

*la* and

*k* =

*sa* to obtain

Next

Eqs. (5) and

(6) are combined with

*j* =

*sa*,

*i* =

*la* to obtain

Eqs. (13) and

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

where unknown value of

*P*_{s}_{,}_{la} is only real positive root of

Eq. (15).

Cerebral fluid flow system

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

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

*q*_{b}_{,} _{f}_{2} is CSF absorption into the neck vein (see ). The equations for resistance

*R*_{b}_{,} _{j} and flow

*q*_{b}_{,} _{j} are the same as

Eqs. (2) and

(3) with

*s* =

*b* and

*k* =

*a*,

*nv*.

Above the neck, the description of the cerebral compartments can be grouped into the cerebral capillary and venule, which are assumed to have Windkessel compliance, and the pial arteries and pial veins, which are assumed to have constant compliance. This distinction is made because prior studies of the Windkessel model have focused on the microvasculature [

36;

37] and not on larger vessels. The cerebral arteriole compliance is treated uniquely due to its involvement in dynamic cerebral autoregulation.

Accounting for the intracranial *ic* pressure and linear compliance *C*_{b}_{,} _{j}_{,0} of the pial arteries *pa* and pial veins *pv*, the volume for *j* = *pa*, *pv* is

The state equations for the pial vessels are obtained from

Eq. (5) with

*s* =

*b; i* =

*na*,

*v*; and

*j* =

*pa*,

*pv*. The cerebral capillary

*c* and venule

*v* compartments include Windkessel compliance

*C*_{b}_{,} _{j} and account for the intracranial pressure

The state equation for

*P*_{b}_{,} _{j,} derived from

Eqs. (19) and

(20), further assumes that the intracranial pressure

*P*_{b}_{,}_{ic} may be time varying

Combining and solving

Eqs. (21) and

(23) for the derivative of pressure results in

The expressions for resistance

*R*_{b}_{,} _{j,} flow

*q*_{b}_{,} _{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

*q*_{b}_{,}_{sp} and autoregulation gain

*K*_{aut} on this feedback loop are influenced by the partial pressure of CO

_{2} in the ISF. The following system equations that govern the CO

_{2} effects are simplified versions of those presented in [

4] because data simulations were found to operate in a linear regime of the more complex equations. The difference from baseline ISF CO

_{2} 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 *x*_{b}_{,}_{aut} is an auxiliary variable that will determine the arteriole compliance changes a few steps later. The dynamics of *x*_{b}_{,}_{aut} are governed by its state equation

The variable cerebral arteriole compliance operates on a non-symmetric sigmoid with the bounds *C*_{b}_{,}_{w} and curvatures *k*_{b}_{,} _{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

*q*_{b}_{,}_{a} falls below the CBF set point

*q*_{b}_{,}_{sp,} the autoregulation state

*x*_{b}_{,}_{aut} will become negative. This will cause the arteriole compliance

*C*_{b}_{,}_{a} to increase towards

*C*_{b}_{,}_{a}_{,0} +

*C*_{b}_{,}_{w}_{1}, which will dilate the arterioles and increase CBF. In the opposite scenario,

*C*_{b}_{,}_{a} decreases towards

*C*_{b}_{,}_{a}_{,0} −

*C*_{b}_{,}_{w}_{1}, which constricts the arterioles and decreases CBF.

The arteriole compliance can now be incorporated into the arteriole state equations. The equations for cerebral arteriole volume

*V*_{b}_{,}_{a}, resistance

*R*_{b}_{,}_{a} and flow

*q*_{b}_{,}_{a} are the same as

Eqs. (1),

(2) and

(3) with

*s* =

*b; j* =

*a*; and

*k* =

*c*. The arteriole state equation is obtained from solving for

*P*_{b}_{,}_{a} in

Eq. (19) with

*j* =

*a*, and then taking the time derivative to obtain

where

*q*_{b}_{,} _{f}_{1} is CSF formation from the arteriole (see ), 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 [

14], which would be superimposed on the regulatory physiology. When this input is non-zero, the arteriole compliance

*C*_{b}_{,}_{a} must be treated as a state variable that is solved by integration and

Eq. (30) would not be used.

In addition to the time-varying fluid flow in the cerebral vascular compartments, the CSF volume in the intracranial space also varies in time because of continuous formation by the choroid plexus and absorption by the arachnoid villi [

38]. Following [

4], the formation rate depends on the pressure gradient

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 *P*_{b}_{,}_{ic} can be expressed in terms of the intracranial compliance *C*_{b}_{,}_{ic} and volume *V*_{b}_{,}_{ic} as

where the intracranial volume is the sum of all intracranial volumes

and the interstitial fluid

*V*_{b}_{,}_{ISF,} intracellular fluid

*V*_{b}_{,}_{ICF} and tissue volume

*V*_{b}_{,}_{tissue} in the intracranial cavity are assumed to be known constants. Given the definition of

*V*_{b}_{,}_{ic} in

Eq. (37), the graphical representation of intracranial compliance

*C*_{b}_{,}_{ic} as a single compartment in is only symbolic. The state equation for intracranial volume

*V*_{b}_{,}_{ic} is

where the inflow through the neck artery

*q*_{b}_{,}_{na}, outflow through the neck vein

*q*_{b}_{,}_{nv}, and CSF absorption

*q*_{b}_{,} _{f}_{2} are the only net effects on the volumes that comprise the intracranial volume

*V*_{b}_{,}_{ic} in

Eq. (37).

The remaining unknown is the intracranial compliance

*C*_{b}_{,}_{ic.} In [

4], the intracranial compliance is defined as inversely proportional to the intracranial pressure. Substituting this inverse proportionality into

Eq. (36), however, results in a constant intracranial volume, which implies zero intracranial compliance. Given this inconsistency, we opted instead to use the Windkessel compliance of

Eq. (7) with

*s* =

*b* and

*j* =

*ic*. Rearranging

Eq. (7) illustrates that larger values of

*β*_{b}_{,}_{ic} result in a compliance expression that is nearly equivalent to the one in [

4] but the algebraic cancelation of

*P*_{b}_{,}_{ic} is avoided

The intracranial pressure state equations that follow from

Eqs. (36) and

(39) are the same as

Eqs. (10) and

(11) with

*s* =

*b* and

*j* =

*ic*.

A value of

*β*_{b}_{,}_{ic} can be chosen based Monroe-Kellie doctrine [

39;

40], which states that the cranium is rigid and any volume change in one of three constituents (blood, CSF, and brain tissue) must be compensated by change in another constituent, which implies that

*β*_{b}_{,}_{ic} is infinitely large. Monroe-Kellie also permits, however, the ability of the CSF to be displaced into the spinal canal, and some compliance in the falx cerebri between the hemispheres and the tentorium between the hemispheres and the cerebellum, which implies that

*β*_{b}_{,}_{ic} is not infinite. We have, for the moment, set

*β*_{b}_{,}_{ic} =5 based on the empirical observation that this value predicts NIRS data reasonably well and we leave further exploration of this parameter to future studies that include intracranial pressure recordings.

Systemic gas exchange

The systemic gas exchange model tracks the bound and unbound O_{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}_{2} in blood is the sum of bound and unbound O_{2}

The bound O_{2} concentration depends on the percent oxygen saturation *S*_{O}_{2}, the hemoglobin concentration *M*_{HbT}_{,}_{blood}, gas constant *R*, standard temperature *T*_{STP} and pressure *P*_{STP}

The percent oxygen saturation

*S*_{O}_{2} can be calculated from the partial pressure of oxygen

*P*_{O}_{2} using an analytic hemoglobin dissociation curve which depends on temperature

*T* and pH [

41] (included for convenience in

Appendix *Blood oxygen*). The unbound O

_{2} concentration simply depends on the partial pressure of oxygen

*P*_{O}_{2} and the gas solubility

*α*_{O}_{2}In the derivation of state equations that follows, it will be necessary to convert oxygen concentration *C*_{O}_{2} state equations into oxygen partial pressure *P*_{O}_{2} state equations. This conversion can be accomplished using the relation

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]

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

The bound and unbound CO

_{2} concentrations depend on the partial pressure of carbon dioxide

*P*_{CO}_{2}, percent oxygen saturation

*S*_{O}_{2}, hematocrit

*Hct*, gas solubility

*α*_{CO}_{2}, temperature

*T* and pH [

42] (see

*Blood carbon dioxide* in

Appendix). Also similar to the treatment of oxygen is the need for an effective solubility

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}_{,}_{g,a,} which is specified as an input in the form of

Eq. (12). The blood gas inputs effectively occur at the arteriole chamber because gas exchange is assumed not to occur in any of the arterial chambers. The user may also supply inputs that change the systemic metabolic rates

*Q*_{s}_{,}_{g} in the form of

Eq. (12) to capture, for example, the effects of a stimulant drug. Moving from the arteriole

*a* chamber to the capillary

*c* chamber will require accounting for the systemic

*s* blood gas concentration

*C*_{s}_{,}_{g}_{,}_{c} for

*g* = O

_{2}, CO

_{2}
where

*V*_{s}_{,}_{g}_{,}_{c} is the gas volume and

*V*_{s}_{,}_{c} is the blood volume. Differentiating

Eq. (48) yields

where the derivative of gas volume

*V*_{s}_{,}_{g}_{,}_{c} depends on the gas inflow from the arteriole, outflow from the capillary and diffusion with ISF

where

*D*_{s}_{,}_{L}_{,}_{g} is the diffusion constant,

*P*_{s}_{,}_{g,c} is the partial pressure in the capillary and

*P*_{s}_{,}_{g}_{,}_{ISF} is the partial pressure in the ISF. The state equation in the capillary is now obtained by combining

Eqs. (48),

(49) and

(50) and dividing by the effective solubility [

Eq. (45) or

(47) as appropriate] to arrive at partial pressure in terms of gas species

*g* A similar state equation governs the venule *v* chamber

The state equation in the intracellular chamber includes the effects of metabolism *Q*_{s}_{,}_{g} and uses the solubility for unbound gas *α*_{g} only because no hemoglobin is present

A similar expression is use for the interstitial compartment

Cerebral gas exchange

Following the same procedure as with the systemic compartments, preliminary calculations are carried out to obtain the cerebral blood gas concentration

*C*_{b}_{,} _{g}_{,} _{j} and effective solubility

_{b}_{,} _{g}_{,} _{j}. User inputs of the cerebral metabolic rates can model changing metabolic demands

The state equations for the cerebral gas exchange system follow a derivation similar that for

Eq. (51). Where needed, additional terms account for gas transport with CSF formation

*q*_{b}_{,} _{f}_{1} and absorption

*q*_{b}_{,} _{f}_{2}. 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.

Hemoglobin concentration changes

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

The first assumption is that the concentration of hemoglobin in the blood *M*_{HbT}_{,} _{blood} is constant. This assumption simplifies the conversion of time-varying blood oxygen saturation *S*_{r}_{,}_{O}_{2,} _{j} in the *j*th 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

*j*th chamber in the systemic or cerebral region (

*r* =

*s*,

*b*) is the state variable

*V*_{r}_{,} _{j} divided by the initial volume

*V*_{r}_{,} _{j}_{,0.} The effective hemoglobin concentration

_{r}_{,}_{h}_{,} _{j} for

*h* =

*HbO*,

*HbR* is then

where

*M*_{r}_{,}_{h}_{,} _{j} is the hemoglobin concentration from

Eq. (62) or

(63) as appropriate. The effective hemoglobin concentration

_{r}_{,}_{h}_{,} _{j} now includes the dynamics from the lumped gas exchange and lumped circulatory models but does not yet reflect the composition of the spatially distributed tissues.

Third and forth assumptions are now needed to incorporate tissue composition and complete the relationship between lumped parameters and distributed tissue properties in a local region. The third assumption is that an initial blood volume fraction

*F*_{r}_{,}_{bv} can be used as a weighting factor in the systemic and cerebral regions (

*r* =

*s*,

*b*) to account for the local blood-tissue composition. The forth assumption is that initial volume fractions of the vascular chambers

*F*_{r}_{,} _{j} can account for the relative composition of the vascular chambers that comprises the local blood content of the tissue. Combining the third and forth assumptions relates the effective hemoglobin concentration in the vascular chambers

_{r}_{,}_{h}_{,} _{j} to the spatially distributed hemoglobin concentration in the tissue

such that

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

*h* in the systemic and cerebral regions (

*r* =

*s*,

*b*). Since the volume fractions

*F* are held constant, all the temporal variation in

_{r}_{,}_{h} comes from dynamics in blood volume

*V*_{r}_{,} _{j} and blood oxygen saturation

*S*_{r}_{,}_{O}_{2,} _{j} as computed by the lumped parameter components of the SimCVR model.

Simulating near-infrared spectroscopy measurements

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

First the effective pathlength

*L*_{r}_{,}_{λ} of the near-infrared light in each region

*r* and each wavelength

*λ* is estimated. For this step we used a Monte Carlo simulation of photon migration in the tissue, which is a numerical solution to radiative transport [

43]. The Monte Carlo method allows for spatially varying optical properties and complex geometries. The resulting effective pathlength distributions characterize the spatial sensitivity for each source-detector pair to each tissue region. Details of how this Monte Carlo method is used to obtain pathlengths

*L*_{r}_{,}_{λ} are already published [

44;

45;

46].

In the second step, the calculated pathlengths

*L*_{r}_{,}_{λ} are used in a first-order Rytov expansion [

47] to model changes in optical density (ΔOD) that result from perturbations in the absorption coefficient Δ

*μ*_{a}_{,}_{r}_{,}_{λ}
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 [

48]. 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.