## Abstract

Cerebral blood flow regulation is based on a variety of different mechanisms, of which the relative regulatory role remains largely unknown. The cerebral regulatory system expresses two regulatory properties: cerebral autoregulation and neurovascular coupling. Since partly the same mechanisms play a role in cerebral autoregulation and neurovascular coupling, this study aimed to develop a physiologically based mathematical model of cerebral blood flow regulation combining these properties. A lumped parameter model of the P2 segment of the posterior cerebral artery and its distal vessels was constructed. Blood flow regulation is exerted at the arteriolar level by vascular smooth muscle and implements myogenic, shear stress based, neurogenic, and metabolic mechanisms. In eight healthy subjects, cerebral autoregulation and neurovascular coupling were challenged by squat-stand maneuvers and visual stimulation using a checkerboard pattern, respectively. Cerebral blood flow velocity was measured using transcranial Doppler, whereas blood pressure was measured by finger volume clamping. In seven subjects, the model proposed fits autoregulation and neurovascular coupling measurement data well. Myogenic regulation is found to dominate the autoregulatory response. Neurogenic regulation, although only implemented as a first-order mechanism, describes neurovascular coupling responses to a great extent. It is concluded that our single, integrated model of cerebral blood flow control may be used to identify the main mechanisms affecting cerebral blood flow regulation in individual subjects.

- cerebral blood flow regulation
- physiological model
- myogenic regulation
- neurogenic regulation
- visually evoked flow response

cerebral blood flow (cbf) regulation is based on a variety of different mechanisms, of which the relative regulatory role remains largely unknown. The cerebral regulatory system expresses two regulatory properties that work partly using the same mechanisms: cerebral autoregulation (CA), maintaining a constant CBF with changing blood pressure, and neurovascular coupling (NVC), adapting CBF to brain activity.

Dysfunction of CA or NVC is associated with several diseases (50, 65). NVC changes are observed, for example, with preeclampsia (41), type I diabetes (53), and Alzheimer's disease (22, 26, 55). CA changes occur, among others, with hypertension (5), long-lasting diabetes (29), and cerebral trauma (48). With the use of transcranial Doppler ultrasonography (TCD), these changes can be measured noninvasively. As such, investigation of NVC and CA using TCD on a patient-specific basis may serve as a valuable diagnostic tool. It is therefore that this study aimed to develop a mathematical model of CBF regulation to interpret CBF fluctuations. Both CA and NVC are accomplished by changing cerebrovascular resistance by vasoconstriction or vasodilation of arterioles and small arteries, induced by a change in vascular smooth muscle cell (VSMC) tone.

Four regulatory mechanisms are thought to influence VSMC tone (19). These all operate on a time scale larger than the heart beat. Myogenic regulation denotes the influence of circumferential mechanical stresses and/or strains on VSMC tone. Blood vessels respond to transmural pressure elevation with constriction and to pressure reduction with dilation, characterized by a time constant of the order of 2.5 s (35) to 10 s (28). Note that in Lanzarone et al. (35), myogenic regulation is characterized by a delay of 2.5 s, additional to the already mentioned time constant. Shear stress based regulation is initiated by the endothelium producing nitric oxide (NO) in response to an increased shear stress on the vessel wall. The NO produced causes VSMC relaxation, normalizing wall shear stress by increasing vessel diameter. Cerebral, shear stress based vasodilation has been shown at the arterial level (20) and, within a certain flow range, at the arteriolar level (47, 63). Shear stress based regulation is a slow process, which, in noncerebral vascular beds, is characterized by a time constant of the order of 60 s (6, 14). Contrary to the aforementioned mechanisms, neurogenic regulation requires some type of external, neural input. Astrocytes are thought to be able to sense neuronal activity and regulate CBF via their end feet (31). As such, smooth muscle tone of pial arteries and intraparenchymal arterioles is actively regulated by neural activity (34, 40), at a time constant of the order of 4 s (51). Metabolic regulation comprises many pathways. Kuschinsky (34) describes metabolic regulation as a mechanism correcting for a mismatch between the neurogenic flow increase and the flow demand. One pathway for metabolic regulation may lie in venular-arteriolar communication (25) of carbon dioxide (CO_{2}), as modeled by David and colleagues (2, 15, 16, 43, 44). Venous CO_{2} produced by metabolism can diffuse to the arterioles feeding the capillary bed, causing local vasodilation, governed by a time constant of the order of 20 s (61).

In literature, a vast number of CBF regulation models exist. However, these are either too detailed to be fitted patient specifically, e.g., Ref. 4, or purely descriptive, e.g., Ref. 54 (NVC) or Ref. 58 (CA), hampering physiological interpretation. Furthermore, models combining CA and NVC are scarce, a notable exception being the model by Payne (51), which, however, does not discriminate among all four aforementioned regulatory mechanisms. Therefore, in this study, a patient-specific model, describing all four mechanisms, has been developed and verified using NVC and CA data measured on the same subject.

## MATERIALS AND METHODS

### Model Definition

#### Lumped parameter model.

To study CBF regulation, a lumped parameter model was developed (Fig. 1), describing the P2 segment of the posterior cerebral artery (PCA) and its distal vessels, up to the large cerebral veins. This segment was chosen since it supplies the primary visual cortex with blood (1). The model consists of three lumped parts: the PCA, the arteriolar circulation and microcirculation, and the venous circulation. Fixed parameter values are specified in Table 1. Arteriolar parameters are taken from the modeled distal arterioles by Ursino and Lodi (61).

##### POSTERIOR CEREBRAL ARTERY.

The PCA is modeled by a transmission line model. Total PCA compliance is calculated assuming a thick-walled vessel:

##### VENOUS CIRCULATION.

Venous resistance is chosen such that the venous pressure drop is 14 mmHg (61) at baseline flow (*Q*_{bl}), leading to a resistance of 0.26 mmHg·ml^{−1}·min. Venous capacitance is 2.5·10^{−2} ml/mmHg, reflecting a venous RC-time (*R*_{v}*C*_{v}) of 0.4 s (61). All venous circulation lumped parameter model components are taken to be linear and constant. Relevant state equations are specified in the appendix.

##### ARTERIOLAR CIRCULATION.

Arteriolar resistance and capacitance are described by Laplace's law, cf. Ursino and Lodi (61). Total arteriolar vessel wall tension (*T*_{a}) is assumed to be the sum of elastic (*T*_{a,e}), muscular (*T*_{a,m}) and viscous (*T*_{a,v}) tensions:
*h*_{a}) in these equations is formulated by
*T*_{a,max} in *Eq. 6*, according to
*M*_{s} is a measure of smooth muscle activation, which varies from −1 to 1, representing maximum vasodilation and constriction, respectively.

Arteriolar resistance (*R*_{a}) and volume (*V*_{a}) scale with *r*_{a}^{−4} and *r*_{a}^{2}, respectively:
*K*_{a,R} is chosen such that, for baseline arteriolar pressure and radius, *R*_{a} leads to a baseline flow (*Q*_{bl}). *K*_{a,V} ensures that baseline capacitance corresponds to an arteriolar RC-time (*R*_{a}*C*_{a}) of 1 s (51). Note that an explicit equation for arteriolar compliance is absent. Arteriolar compliance arises from the fact that *r*_{a} is pressure-dependent (via Laplace's law, cf. *Eq. 4*) and is therefore an intrinsic property of the model. Because of the Laplace formulation used in this model, *r*_{a} is chosen as the arteriolar state variable (contrary to *P*_{a}, which would have been an obvious choice if *C*_{a} had been a linear capacitance). The state equation is successively obtained by differentiating *Eq. 11* with respect to time. The resulting

#### Regulation.

CBF regulation is performed by varying *M*_{s} in *Eq. 9*. *M*_{s} is coupled to *M*_{s,l} via a sigmoidal function:
*M*_{s,1} is the summation of influences of myogenic (*x*_{myo}), shear stress based (*x*_{shear}), neurogenic (*x*_{neuro}), and metabolic (*x*_{meta}) regulation, multiplied by their respective gains (*G*_{myo}, *G*_{shear}, *G*_{neuro}, and *G*_{meta}):
*x*_{initial} is a parameter to eventually correct for steady-state flow deviations between model and measured flow. An overview of regulation is depicted in Fig. 2.

##### MYOGENIC REGULATION.

Activation of myogenic regulation (*A*_{myo}) is based on the deviation of the current arteriolar wall tension, *T*_{a}, from a reference tension, *T*_{myo,0}, and is normalized by division by a normalization tension, *T*_{myo,s} (7, 8, 56):
*T*_{myo,0} is chosen to be the tension at baseline pressure for *M*_{s} = 0. *T*_{myo,s} is a normalization tension for obtaining a dimensionless *A*_{myo} and is chosen to be 0.3 mmHg·cm. In all four regulatory mechanisms, the regulatory state (*x*_{regulation}) is coupled to the regulatory activation (*A*_{regulation}) via first-order behavior. For myogenic regulation, this equates to
_{myo} is the time constant governing myogenic regulation.

##### SHEAR STRESS BASED REGULATION.

Shear stress is caused by viscous friction of blood and the vessel wall. If viscosity-dominated (Poiseuille) flow is assumed, wall shear stress (τ_{w}) is quantified by
*Q*, *r*, and η are flow, vessel radius, and dynamic viscosity, respectively. From *Eq. 16*, it can be inferred that shear stress scales with
*K*_{shear} is chosen such that, at baseline, *A*_{shear} = 0. *x*_{shear} is coupled to *A*_{shear} via first-order time behavior, analogous to *Eq. 15*. τ_{shear} is the time constant governing this coupling and is fixed at 60 s, a typical time constant observed in other vascular beds (6, 14).

##### NEUROGENIC REGULATION.

Neurogenic regulation, together with metabolic regulation, is assumed to be determined by the amount of visual cortex activation. For the sake of simplicity, neurogenic activation (*A*_{neuro}) is assumed to change stepwise at eye opening, cf. Rosengarten et al. (54). Neurogenic regulation is then again described by a first-order response, analogous to *Eq. 15*.

##### METABOLIC REGULATION.

Metabolic regulation is assumed to be mediated by CO_{2} (15), a potent vasodilator (32). It is assumed that the set of arterioles feeding the capillary bed are in close proximity to the venous return, a concept reviewed by Hester and Hammer (25). Brain tissue CO_{2} concentration (*C*_{t,CO2}) is computed from the molar balance between CO_{2} production and CO_{2} disposal. *C*_{t,CO2} represents a molar concentration. Tissue CO_{2} is produced by metabolism (*M*_{CO2}) and disposed by blood perfusing the tissue. This is formulated in the following equation:
*V* is the brain volume perfused by one PCA and *C*_{v,CO2} and *C*_{a,CO2} are the venous and arterial CO_{2} concentrations, respectively. *V* represents the unilateral primary visual cortex volume and is chosen to be 5.99 cm^{3} (18). Visual cortex stimulation influences metabolism as follows:
*f*_{Q} is the ratio between open eye and closed eye flows at baseline. *n*_{QM} is a factor describing the ratio of procentual flow increase to procentual metabolism increase, which is chosen to be 2.2 (60).
^{−7} mol/s , calculated using *V* and a metabolism of 0.58 ml CO_{2}·kg^{−1}·s^{−1} (42). *C*_{a,CO2} is assumed to be constant at 20.65 mol/m^{3} (21). *C*_{v,CO2} is assumed to be in equilibrium with *C*_{t,CO2}, which is justified by the very fast diffusion of CO_{2} (24). Note that for two dissolved gases to be in equilibrium, their partial pressures (not their concentrations) are to be equal (12).

The relation between venous CO_{2} concentration and partial pressure is described by
^{3} and *Eq. 20* with data on brain tissue CO_{2} concentration and partial pressure measured in dog brains (30), a relationship between *C*_{t,CO2} and *C*_{v,CO2} can be computed (Fig. 3, circles). Linear fitting of the data points yields the following relationship between *C*_{v,CO2} and *C*_{t,CO2} (Fig. 3, line):
_{t,v} = 0.96 and β_{t,v} = 8.9 mol/m^{3}. *C*_{t,CO2} influences metabolic activation, *A*_{meta}, via
*Eq. 18*. *G*_{Ameta} is a scaling gain factor, which is chosen to be 0.59 m^{3}/mol. *x*_{meta} is coupled to *A*_{meta} via first-order dynamics, analogous to *Eq. 15*.

### Experiments

#### Data.

To characterize CA and NVC, measurements were performed on 11 volunteers. Experiments were conducted according to the guidelines formulated by the Institutional Review Board of Maastricht University Medical Center. All subjects gave written informed consent. The following data were acquired: *1*) finger blood pressure; *2*) bilateral TCD blood flow velocity in the left P2 segment of the posterior cerebral artery (PCA) and the right middle cerebral artery (MCA), respectively; *3*) a one-lead, three-electrode electrocardiogram (ECG); and *4*) an electro-oculogram (EOG). Finger blood pressure was acquired continuously using an ambulatory blood pressure monitor (PortaPres, TNO, Amsterdam, The Netherlands), based on the volume-clamping technique (52). TCD signals were acquired using two Doppler transducers mounted to a head frame and connected to a Multi Dop X4 (DWL Compumedics, Singen, Germany). ECG amplification was performed using an in-house developed ECG amplifier. An EOG was measured to determine the exact time of eye opening or closing of the subject. EOG amplification was performed using an in-house developed ECG amplifier. All signals were acquired using a personal computer equipped with a Keithley KPCI-3101 I/O card (Keithley Instruments, Solon, OH) at a sampling rate and bit depth of 250 Hz and 12 bit, respectively.

#### Maneuvers.

##### NEUROVASCULAR COUPLING.

Eleven visually-evoked flow responses (VEFRs) to a 60-s alternating checkerboard stimulus (0.5 s between alternations) followed by 60 s without stimulus were recorded. An audible beep indicated the beginning and end of each stimulus. During this experiment, the subject was sitting in a hospital bed, with the upper body tilted at an angle of ∼70°.

##### BASELINE.

Five-minute baseline pressure and blood flow velocity measurements were performed with (eyes open, alternating checkerboard stimulus) and without (eyes closed) visual stimulation.

##### CEREBRAL AUTOREGULATION.

CA was challenged by blood pressure variations induced by periodic squatting and standing of the subject, cf. Refs. 9, 57. The subject's eyes remained closed during the entire experiment. Maneuvers were conducted just after pressure calibrations (64), ensuring that pressure calibration occurred in the least dynamic period of the measurement cycle. Since pressure calibrations were separated by 70 heart beats, this led to an interval of ∼1 min between squat and stand movements, assuming a heart rate of 70 beats/min. In total, 11 squat-stand cycles were recorded.

### Data Processing and Model Implementation

Continuous blood pressure and flow velocity signals were converted to beat-to-beat averaged values. Data were inspected visually to check and if necessary correct for artifacts. Short periods of missing signal were interpolated to ensure a continuous, smooth recording. Further processing and cutting were performed using MATLAB R2012a (MathWorks, Natick, MA). Beat-to-beat data were resampled equidistantly at 5 Hz using cubic spline interpolation as previously described (23). CA and NVC responses were averaged for each subject. NVC response averaging was performed on an EOG time-locked basis.

The model defined above was implemented using MATLAB. The program code used containing the basic concept of the model has been made available as supplemental material on the *American Journal of Physiology: Heart and Circulatory Physiology* web site and was tested on both MATLAB R2012a and GNU Octave version 3.6.1 (18a). The reader is asked to use the model code cautiously, bearing the model's limitations (see discussion) in mind. In addition to the original model code, a CellML version of the model was prepared in the group of Beard et al. and is available at the CellML model repository (38a) and at the National Simulation Resource Physiome Project model repository (5a).

### Fitting

From the baseline measurement's closed eye period, average pressure and flow velocity were calculated. Average flow velocity is assumed to correspond to baseline flow, yielding a flow-to-flow velocity ratio. Physiologically, this ratio is dependent on the insonated vessel's cross-sectional area and flow profile and on the insonation angle. The baseline measurement was additionally used for calculating the ratio between average open and closed eye blood flow velocities (*f*_{Q}).

Measured finger blood pressure was used as model input. Model arterial blood flow (*Q*_{ai}) was fitted to measured arterial blood flow (*Q*_{ms}) by minimizing the root mean square error
*E*_{RMS}, *Eq. 23*).

To easily refer to the various measurements and maneuvers performed, the following terms are defined. For CA, an up maneuver refers to the subject's postural change from squatting to standing. The resulting fluctuations in physiological quantities are termed the up response. For NVC, an open maneuver refers to the subject opening her/his eyes, causing an open response.

CA experimental data were used to estimate (fit) myogenic (0 s ≤ τ_{myo} ≤ 15 s and 0.1 ≤ *G*_{myo} ≤ 25) and shear stress based (−20 ≤ *G*_{shear} ≤ 0) regulation parameters, as well as a steady-state correction parameter (−1 ≤ *x*_{initial} ≤ 1). Neurogenic and metabolic regulations were switched off (*G*_{neuro} = 0 and *G*_{meta} = 0). Fitting was performed on averaged up responses. The up response was chosen for fitting since blood pressure fluctuations when standing up are larger than when squatting, yielding a larger regulatory stimulus. NVC experimental data were used to estimate (fit) neurogenic (1 s ≤ τ_{neuro} ≤ 20 s and −5 ≤ *G*_{neuro} ≤ 0) regulation parameters and again a steady-state correction parameter (−1 ≤ *x*_{initial} ≤ 1). Myogenic and shear stress based parameters were fixed to their values found using the CA experiment data. Metabolic regulation was not included in the fitting procedure, based on preliminary fits systematically yielding *G*_{meta} = 0. This choice is further elaborated in the discussion. The open response was chosen for fitting because this type of response had been already extensively studied [e.g., Rosengarten et al. (54)]. Preliminary studies showed that response averaging did not yield statistically different parameter values from fitting nonaveraged results for both CA and NVC (data not shown).

## RESULTS

### Measurements

#### Included subjects.

Eleven subjects were initially included in this study. One subject was excluded due to inability to locate the PCA. Two other subjects, in whom measured PCA and MCA flow velocity responses during NVC measurements were almost equal, were also excluded. The remaining eight subjects were used for further analysis. General subject properties are given in Table 2. All subjects were healthy, nonsmokers without any known cardiovascular disease. Five out of eight subjects wore glasses. This may be of influence, since, because of the TCD fixation mask, wearing of glasses is not possible. This could cause a slightly lower NVC response due to the fact that these subjects could not focus sharply on the screen.

#### Measurements.

Mean pressures and flow velocities for baseline measurements are shown in Table 3. Baseline blood pressures varied among subjects in a range of 28 mmHg. VEFRs ranged from 13 to 24%, which is low compared with other studies, e.g., Ref. 49. Flow increases in subjects wearing glasses were tested against flow increases in subjects not wearing glasses; groups were not statistically different (*P* = 0.82, Wilcoxon rank sum test). Velocities measured with eyes open differed significantly from velocities measured with eyes closed (*P* = 0.008, Wilcoxon paired signed rank test).

CA data are shown in Fig. 4, *A* and *B*. After a blood pressure drop caused by an up maneuver, PCA flow returned to its original value faster than blood pressure did, indicating autoregulatory activity. In the intersubject averaged up response, at *t* = 5 s, pressure was even still decreasing whereas PCA flow velocity was already increasing. NVC data are shown in Fig. 4, *C* and *D*. After eye opening, a visually evoked flow velocity increase was observed, which, in some subjects, showed an overshoot. In literature, for healthy subjects, an overshoot is reported [e.g., Martens et al. (41)]. Flow velocity responses in this study, however, only showed a modest overshoot for some of the subjects. After eye opening, between *t* = 0 s and *t* = 10 s, a small blood pressure increase was observed.

### Model Fitting

#### Representative subject.

Baseline data are used to find a flow velocity to flow conversion factor, and to find *f*_{Q}, which is the ratio of steady state open-eye and closed-eye flows. For a representative subject (*number 7*), time courses of measured and model pressures and flows, and of regulatory state variables are presented in Fig. 5. Figure 5*C* shows that for CA, myogenic regulation (*G*_{myo}·*x*_{myo}) has a large influence on *M*_{s} dynamics, whereas shear stress based regulation (*G*_{shear}·*x*_{shear}) has a much smaller amplitude. In NVC, myogenic regulation is of opposite sign than neurogenic regulation (*G*_{neuro}·*x*_{neuro}; Fig. 5*F*). Eye opening causes a small decrease in shear stress (visible as an increase in *G*_{shear}·*x*_{shear}, since *G*_{shear} < 0). Fitted parameter values for the representative subject are *G*_{myo} = 3.5, τ_{myo} = 6.5 s, *G*_{shear} = −2.5, *G*_{neuro} = −0.73, and τ_{neuro} = 10 s.

#### Intersubject parameter variability.

Fitted parameter values for all eight subjects are represented by the box plots in Fig. 6*A*. Fit quality is represented by the root mean square error, shown in Fig. 6*B*. Generally, errors range from 2 to 5%, indicating a good fit. In one subject, however, the error of the CA fit was >9%.

## DISCUSSION

### Main Findings

In this study, a mathematical model of CBF regulation was developed and verified, describing both CA and NVC. Modeled arteriolar wall tension is actively modulated by physiological mechanisms known to influence CBF. Maneuvers exciting both CA and NVC were used to trigger regulatory responses in eight healthy volunteers. In CA, myogenic regulation dominated the regulatory response. In NVC, a first-order neurogenic response described the VEFR to a great extent. Furthermore, due to the interaction of neurogenic and myogenic regulation, the typical flow overshoot observed in VEFRs can arise.

The model in this study mainly differs from the model by Payne (51) in two ways. First, Payne chooses to model CA based on a flow difference from a setpoint. The exact physiological mechanism (e.g., myogenic or shear stress based regulation) that is modeled, however, remains unspecified. In contrast to this, we choose to base our model on four mechanisms which elicit CA and NVC. Secondly, while Payne used a second-order neural model, we showed that a first-order neural model is sufficient to describe NVC dynamics.

### Cerebral Autoregulation Responses

In all subjects, myogenic regulation dynamics are of much larger amplitude than shear stress based regulation dynamics (e.g., Fig. 5, *C* and *F*). This implies that myogenic regulation plays a larger role in correcting for pressure-induced CBF fluctuations than shear stress based regulation does.

Myogenic regulation and shear stress based regulation are two counteracting mechanisms. An autoregulatory pressure drop leads to myogenic vasodilation, but also to diminished flow, causing less shear stress and therefore relative vasoconstriction. Shear stress regulation, therefore, “counteracts” CA, a concept also observed by Carlson et al. (7), Cornelissen et al. (13), and Liao and Kuo (38). It should be noted that this counteraction takes place at a different time scale (i.e., τ_{myo} ≈ 0.1τ_{shear}). Nonetheless, if CA is studied on a sufficiently large time scale (i.e., if static autoregulation is studied), in order for CA to work correctly, myogenic regulation should dominate shear stress based regulation.

Because of the relatively small contribution of shear stress based regulation to the CA response, fitting its time constant (τ_{shear}) proved difficult. It was therefore decided to fix τ_{shear} and omit it from the fitting procedure, decreasing fitting uncertainty of other parameters.

### Neurovascular Coupling Responses

It could be hypothesized that, in order for the VEFR to show an overshoot, a second-order neurogenic model is required, cf. Rosengarten et al. (54). However, model simulations incorporating only myogenic and neurogenic regulation (both of first order) clearly show an overshoot (Fig. 7). This overshoot occurs because the initial rise caused by neurogenic regulation is counteracted by a myogenic response. A similar observation is described by Trillenberg et al. (59) in which the NVC response is described by the sum of a pulse and a first-order response.

Metabolic regulation, as implemented in this model, has a typical effect. From fMRI studies (60), it is known that increased brain activity leads to larger CBF changes than metabolism changes. If, during a stimulus, CBF increases more than metabolism does, the CO_{2} concentration equilibrium will decrease, since CO_{2} removal (by CBF) changes to a greater extent than CO_{2} production (by metabolism) does (*Eq. 18*). Therefore, it can be concluded that metabolic regulation as implemented cannot be responsible for the NVC response on its own. However, because of the net vasoconstrictive effect occurring with an increased metabolism, it could be hypothesized that metabolic regulation is responsible for the overshoot in the VEFR generally observed, counteracting a fast, neurogenic flow increase. Indeed, such a response can be simulated (Fig. 8).

Metabolic regulation was disabled during NVC response fits. Preliminary NVC fits showed that *G*_{meta} was systematically fitted to 0. Therefore, to increase fitting performance and to reduce uncertainty in other fitted parameters, metabolic fitting was abandoned. This is conceivable from a physiological point of view, as metabolic regulation is sometimes described as a “back-up” mechanism, being effective only if neurogenic regulation is impaired or ineffective (34). Flow responses in our study were low compared with other studies, e.g., Panczel et al. (49), and may not have triggered a metabolic response.

Shear stress based regulation is found to be of little influence. However, the “sign” of its contribution (either vasoconstrictive or vasodilative) in this model during VEFRs is remarkable. Under Poiseuille flow assumptions, shear stress (τ) scales with
*Eq. 16*); *Q* being flow and *r* being vessel radius. *Q* scales with 1/R and therefore with *r*^{4} (cf. *Eq. 3*); *R* being resistance. Combining these scalings gives τ ∝ *r*, i.e., shear stress increases linearly with vessel radius. After eye opening (and consequent vasodilation, implying an increase of *r*), however, shear stress is found to decrease in the model (Fig. 5*F*). This, at first sight unexpected observation, can be explained as follows. In the lumped parameter model, *R* is composed of a PCA, arteriolar, and venous part, of which only the arteriolar part is regulated. This leads to a smaller change in resistance and flow than is predicted by *r*^{4}. Shear stress, however, does scale with

### Limitations

The model developed in this study exerts CBF regulation only at one level in the vasculature, which is physiologically incorrect, since it is known that CBF regulation is heterogeneous along the cerebrovascular bed (17). Shear stress regulation is often described to operate in the larger arterioles or small arteries, whereas myogenic regulation is described to be effective downstream the vascular tree. If such a spatial division of flow regulation would be modeled, spatial heterogeneity of flow regulation could be studied in greater detail. However, in this study, for simplicity, CBF regulation was exerted at only one level in the model.

Arterial blood flow (*Q*_{ai}) is assumed to scale linearly with measured flow velocity in this study, assuming a constant vessel radius at the measurement location (the P2 segment of the PCA). Since the PCA is a distensible vessel, this assumption is formally invalid. In a future model, distensibility at the TCD site could be implemented by assuming a monoexponential pressure-radius relationship, analogous to Ursino and Lodi (61). PCA radius changes are, however, small, making the error caused by assuming a constant radius acceptable (see Ref. 33 and Ref. 36 for further discussion on this subject).

An important modeling assumption is that at measured baseline pressure and flow, regulatory influences of myogenic, shear stress based, neurogenic, and metabolic regulation are zero (i.e., *M*_{s} = 0), allowing an equal amount of vasodilation (*M*_{s} = −1) or vasoconstriction (*M*_{s} = 1) with respect to baseline conditions to take place. For one subject in this study, this yielded an insufficient regulation range, leading to a large root mean square error of 9.4% of the CA fit (upper bound in Fig. 6*B*). Physiologically, this assumption is not necessarily true: pressure and flow measured during baseline measurements do not necessarily represent “medium” VSMC tone. An alternative for this assumption would be to use reference values from literature to determine a regulatory reference. However, these are not subject-specific and discard intersubject differences in CA range and baseline pressure and flow. Additionally, the smooth muscle cell component of Laplace's law used in this study is based on parameter values for small pial arterioles (61). However, the smooth muscle cell component of the model developed represents pial arterioles of all sizes. Because of the bell-shaped relationship between vessel radius and smooth muscle tension used (*Eq. 6*), maximum vasoconstriction and vasodilation are highly dependent on the instantaneous vessel radius. More specifically, the amount of vasoconstriction or vasodilation corresponding to a certain value of *M*_{s} depends on the difference between *r*_{a} and *r*_{a,m}, being the instantaneous arteriolar vessel radius and the vessel radius at which the modeled smooth muscle cell develops maximal tone (*Eq. 6*), respectively. Therefore, the gains found for different regulatory mechanisms are to some extent dependent on the operating point in the smooth muscle radius-tension relationship. This problem could be overcome by either choosing *r*_{a,m} subject-specifically (probably using baseline measurement data) or by using another way of coupling regulatory influences to the lumped parameter model. Ursino et al. (62) and Payne (51), for example, regulate arterial compliance based on a sigmoidal, asymmetrical relationship.

In future studies, capnogram data should be studied and, if large fluctuations are found, an additional model input might be constructed modeling the vasodilatory effect of arterial CO_{2}, e.g., Ursino and Lodi (61).

In this study, an alternating checkerboard stimulus is used. However, preliminary experiments in our laboratory showed that a more complex visual stimulus may enlarge the VEFR and its overshoot, which confirms similar findings by Panczel et al. (49). A larger flow response may increase signal-to-noise ratio and trigger metabolic regulation, which could improve fitting of metabolic regulatory mechanisms.

The neurogenic time constants (τ_{neuro}) found using our model are larger than those reported by e.g., Payne (51) (8 s vs. ∼4 s). However, VEFRs in their study clearly show an overshoot, whereas measured VEFRs in our study show a small overshoot at most. If an overshoot is absent, the neurogenic time constant will be fitted to a value approximately equal to the myogenic time constant, ensuring that when neurogenic and myogenic regulation are combined, no overshoot arises (first-order behavior). Therefore, it is expected that the larger overshoot occurring with a more complex stimulus would cause the fitted τ_{neuro} to decrease.

The model proposed can be used to interpret measurement data in different patient groups with specific pathological conditions. Parameter values obtained in this way can be compared with literature, either quantitatively or qualitatively, providing for another model validation step. A parameter study comparing characteristic CBF responses commonly observed in clinical studies is currently ongoing (Martens et al., in preparation).

In conclusion, our developed, integrated model of cerebral blood flow regulation describes both cerebral autoregulation and neurovascular coupling responses well. As such, it may be used to identify the main mechanisms affecting cerebral blood flow regulation in individual subjects, owing to the fact that the four regulatory mechanisms involved are explicitly modeled.

## GRANTS

This study was partially supported by a Kootstra talent fellowship awarded to B. Spronck by Maastricht University, Maastricht, The Netherlands.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: conception and design of research: B.S., E.G.H.J.M., E.D.G., and F.N.v.d.V. Performed experiments: B.S., E.G.H.J.M., and E.D.G. Analyzed data: B.S., E.G.H.J.M., E.D.G., and F.N.v.d.V. Interpreted results of experiments: B.S., E.G.H.J.M., E.D.G., and F.N.v.d.V. Prepared figures: B.S. Drafted manuscript: B.S. Edited and revised manuscript: E.G.H.J.M., E.D.G., and F.N.v.d.V. Approved final version of manuscript: B.S., E.G.H.J.M., E.D.G., and F.N.v.d.V.

## ACKNOWLEDGMENTS

D. A. Beard and C. T. Thompson are gratefully acknowledged for translating the MATLAB program code into CellML.

## APPENDIX

#### State Equations

To facilitate a straightforward reproduction of the mathematical model, all required lumped parameter model state equations are summarized in this appendix. Model simulation requires solving a set of nine state equations: two for the PCA (*Eqs. A2* and *A3*), one for the arteriolar circulation (*Eq. A4*), and one for the venous circulation (*Eq. A5*). In addition, five regulation equations must be solved: *Eq. 15* (myogenic regulation) and its analogs for the other three regulatory mechanisms; and a brain tissue CO_{2} concentration equation (*Eq. 18*).

#### Posterior Cerebral Artery

The PCA model involves three state equations:
*P*_{a} in *Eq. A3* directly follows from *r*_{a} via Laplace's law (*Eq. 4*).
*Eq. A4*).

*Equation A1* concerns the proximal PCA compliance. Since *P*_{ai} is prescribed, this state equation is not formally required to solve the set of model equations. However, it can be used to calculate *Q*_{ai} after simulation.

#### Arteriolar Circulation

Although the arteriolar compliance is depicted as a normal capacitance in Fig. 1, it is described indirectly by Laplace's law. Hence,
*Eq. 11* to time, yielding an expression for
*Q*_{a}) is found using Ohm's law for the distal arteriolar and proximal venous resistance, combined with Laplace's law (*Eq. 4*). Note that the viscous term of this equation (*T*_{a,v}; *Eq. 7*) has been explicitly filled in in *Eq. A4* to properly account for its

#### Venous Circulation

The venous circulation is represented by a single state equation:

- Copyright © 2012 the American Physiological Society

## REFERENCES

- 1.↵
- 2.↵
- 3.
- 4.↵
- 5.↵
- 5a.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 18a.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.
- 38.↵
- 38a.↵
- 39.
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.
- 46.
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵