A quantitative analysis of cell-specific contributions and the role of anesthetics to the neurovascular coupling

The neurovascular coupling (NVC) connects neuronal activity to hemodynamic responses in the brain. This connection is the basis for the interpretation of functional magnetic resonance imaging data. Despite the central role of this coupling, we lack detailed knowledge about cell-specific contributions and our knowledge about NVC is mainly based on animal experiments performed during anesthesia. Anesthetics are known to affect neuronal excitability, but how this affects the vessel diameters is not known. Due to the high complexity of NVC data, mathematical modeling is needed for a meaningful analysis. However, neither the relevant neuronal subtypes nor the effects of anesthetics are covered by current models. Here, we present a mathematical model including GABAergic interneurons and pyramidal neurons, as well as the effect of an anesthetic agent. The model is consistent with data from optogenetic experiments from both awake and anesthetized animals, and it correctly predicts data from experiments with different pharmacological modulators. The analysis suggests that no downstream anesthetic effects are necessary if one of the GABAergic interneuron signaling pathways include a Michaelis-Menten expression. This is the first example of a quantitative model that includes both the cell-specific contributions and the effect of an anesthetic agent on the NVC.


Introduction
The temporal and spatial connection between the brain's neuronal activity and its hemodynamic responses, the neurovascular coupling (NVC), is a central aspect of brain function. This connection is the result of a complex interplay between neuronal cells and blood vessels which affects the cerebral blood flow (CBF), the cerebral blood volume (CBV), and the cerebral metabolic rate of oxygen (CMRO 2 ). It is essential to understand the underlying mechanisms of NVC to fully interpret neuroimaging, such as functional magnetic resonance imaging (fMRI) (Ogawa et al., 1990), which use hemodynamic changes as a proxy for neural activity (Logothetis et al., 2001;Shmuel et al., 2006).
To search for neuronal correlates of the fMRI signal, researchers rely on animal experiments, which normally are performed under anesthesia. Unfortunately, the commonly used anesthetics have a broad effect on brain cells, including changes in neuronal excitability, vascular reactivity, cerebral metabolism, and other physiological processes (Masamoto and Kanno, 2012;Gao et al., 2017). The quantitative effects of these changes vary between different anesthetics, and those detailed variations remain unknown. Nevertheless, the majority of anesthetics exert their effect by direct action on ion channels, like the ionotropic γ-aminobutyric acid (GABA)-activated GABA A and N-methyl-D-aspartate (NMDA)-activated glutamate receptors, or different types of K þ channels (Franks, 2008;Franks and Lieb, 1994) (Fig. 1A, upper arrow). However, how these effects propagate from altered ion-channel activities to other NVC aspects remain to be determined (Fig. 1A, lower arrow). This means that a precise and correct interpretation of the NVC and its alterations during anesthesia needs to be obtained for an exhaustive interpretation of fMRI data (Iadecola, 2017).
A basic chain of events, describing the NVC, in both awake and anesthetized animals, is slowly emerging. This chain of events (summarized in Fig. 1B) starts with the release of excitatory neurotransmitters, such as glutamate, and ends with the release of multiple vasoactive signaling molecules acting on pericytes and vascular smooth muscle cells (VSM) in the vessel walls (Attwell et al., 2010;Hillman, 2014). Activation of both excitatory and inhibitory neurons dilate cerebral vessels. More specifically, activation of excitatory pyramidal neurons release cyclooxygenase-2 (COX-2)-derived prostaglandin E2 (PGE 2 ), which in turn dilates arterioles (Lecrux et al., 2011;Lacroix et al., 2015) (Fig. 1B, grey triangle). Activation of inhibitory (i.e., GABAergic) interneurons can release nitric oxide (NO), a potent vasodilator, and different vasoactive neuropeptides with both vasodilative (vasoactive intestinal peptide, VIP) or vasocontractile (neuropeptide Y, NPY; somatostatin) effects (Fig. 1B, yellow oval). The type of vasoactive molecules released depends on which type of the interneuron is activated (Markram et al., 2004;Cauli et al., 2004;Cauli and Hamel, 2010). In addition, astrocytes release arachidonic acid-derived messengers and K þ ions, and affect CBF by acting on capillaries (Mishra et al., 2016) and arterioles (Filosa et al., 2006;Bazargani and Attwell, 2016;Mishra, 2017) (Fig. 1B, green rectangle). In summary, the NVC depends on a complex interplay between multiple cellular sources and targets, and the magnitude of each specific cellular contribution remains undetermined (Fig. 1A).
This complex interplay has been studied by the use of optogenetic (OG) stimulation, a technique where light-gated ion channels e.g., channelrhodopsins, are expressed in the cellular plasma membrane and opened by exposure to light of a specific wavelength. In a recent OG study, Uhlirova et al. (2016), investigated the cell type-specific contributions to the NVC by selectively expressing channelrhodopsin-2 (ChR2) channels in either pyramidal neurons or GABAergic interneurons in mice. They measured arteriolar diameter changes in vivo upon both light-induced and sensory-induced stimulation in awake and anesthetized animals ( Fig. 1C and D). They found that a light-induced activation of GABAergic interneurons evoked a biphasic (i.e., an overshoot followed by a post-peak undershoot) hemodynamic response (Fig. 1C). Compared to awake condition (Fig. 1C, red symbols), anesthesia caused a prolonged response, lengthening the complete response by a factor of~2.5 (Fig. 1C, black symbols). In contrast to these OG responses, sensory stimulation in awake or anesthetized animals evoked similar biphasic responses, with a larger post-peak undershoot in anesthetized animals (Fig. 1D, black and red symbols). In addition, Uhlirova and colleagues showed that NPY was responsible for the post-stimulus undershoot. This finding provides strong evidence for the active role of interneurons in facilitating and shaping the vascular response in the NVC. However, while such purely experimental research is essential, a quantitative systems-level understanding is still missing.
A quantitative systems-level understanding can be obtained using systems biology and mathematical modeling. A quantitative understanding is obtained if all remaining differences between simulations and experimental data can be viewed as measurement uncertainty (Cedersund and Roll, 2009;Cedersund, 2012). This stands in contrast to a qualitative understanding, which only requires that the qualitative features in experimental datasuch as the presence of an overshoot or a post-peak undershootis captured by the simulations. Systems biology models provide a systems-level understanding if the models are based on biological mechanisms that together produce the model simulations. In other words, systems biology is not a statistical black-box method, but a method to formally test mechanistic hypotheses using experimental data (Cedersund and Roll, 2009). Previously, various modeling approaches have been used for NVC analysis: e.g., multiple models describing phenomenological NVC mechanisms with a basis in fMRI (Buxton et al., Fig. 1. A: Open questions regarding the NVC and the effects of anesthesia. B: Simplified schematic illustration of the cellular pathways underlying the NVC. Neuronal activity activates GABAergic interneurons (yellow oval), pyramidal neurons (grey triangle), and astrocytes (green rectangle), by stimulating influx of Ca 2þ . Ca 2þ acts as a second messenger, triggering signaling pathways in the different cells. In GABAergic interneurons, Ca 2þ promotes nitric oxide (NO) through upregulation of nitric oxide synthase (NOS) as well as facilitating the release of different neuropeptides such as neuropeptide Y (NPY), vasoactive intestinal peptide (VIP), and somatostatin (SOM). In pyramidal neurons, Ca 2þ promotes the synthesis of the arachidonic acid (AA) via phospholipase A2 (PLA 2 ), which in turn is metabolized to prostaglandin E2 (PGE 2 ) via cyclooxygenase-2 (COX-2). In astrocytes, arachidonic acid is synthesized via phospholipase D2 (PLD 2 ) and its metabolites PGE 2 and epoxyeicosatrienoic acid (EET) are synthesized via cyclooxygenase-1 (COX-1) and cytochrome P450 (CYP) epoxygenase, respectively. Furthermore, extracellular K þ is increased upon neuronal activity. Together these vasoactive messengers act on arterioles and capillaries to modulate the vessel diameters, with the neuronal messengers acting primarily on the arterioles, and astrocytes acts on both arterioles and capillaries. C & D: Characteristics of the arteriolar diameter responses extracted from Uhlirova et al. (Uhlirova et al., 2016). Optogenetic (OG) (C) and sensory (D) stimuli-induced arteriolar diameter changes (y-axis) in awake and anesthetized animals (red and black symbols, respectively). The error bars depict the standard error of the mean (SEM). The bars in the lower left corner of each graph indicates the stimulus length for the OG light pulse of 450 ms (black) and 400 ms (red) and the sensory stimulation of 2 s (black) and 1 s (red) for awake and anesthetized animals, respectively. 1998; Mandeville et al., 1999;Vazquez et al., 2006;Havlicek et al., 2015;Kim et al., 2013;Kim and Ress, 2016) or optical imaging (Huppert et al., 2007(Huppert et al., , 2009. In recent years, mechanistic models, capturing the interplay between dilating and constrictive vasoactive substances produced as a result of intracellular signaling has gained traction (Yucel et al., 2009;Mathias et al., 2016Mathias et al., , 2018Lundengård et al., 2016;Sten et al., 2017). However, none of these models capture the precise mechanisms elucidated in Uhlirova et al. (2016). To our knowledge, no previous model has used OG data to either describe the interplay between GABAergic interneurons and pyramidal neurons in the NVC or their alterations during anesthesia.
Herein, we evaluate the experimental data presented in Uhlirova et al. (2016), using mechanistic systems biology modeling. We present a model capturing the mechanisms and interplay between GABAergic interneurons and pyramidal neurons with respect to data of arteriolar diameter changes (Fig. 2B). The model consists of three different levels: rapid neuronal activity, intracellular signaling pathways, and vascular responses. The model fits experimental data used for model training, featuring both awake and anesthetized animals. The model also correctly predicts experimental data not used for training, capturing the effect of different pharmacological inhibitors. In the model, the anesthetic agent acts on the rapid neuronal activity, and the OG-induced prolonged undershoot observed during anesthesia (Fig. 1C) can be explained if the vasoconstrictive effect of NPY follows Michaelis-Menten kinetics. In summary, our model provides a first example of a quantitative systems-level understanding of the cell-specific contributions to the NVC Fig. 2. Representation of the neurovascular system (A) and the model interaction graph (B). A: Two distinct types of cells are described: Pyramidal neurons (grey triangle, left) and GABAergic interneurons (further divided into NO-positive and NPY-positive; yellow oval shape, right). In pyramidal neurons, glutamate activates N-methyl-D-aspartate receptors (NMDAR) and α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic receptors (AMPAR), resulting in an influx of Na þ and Ca 2þ ions, which lead to depolarization. γ-aminobutyric acid (GABA) activates GABA A receptors (GABA A R) which prevents depolarization, thus counteracting the effect of glutamate. This balance between GABA and glutamate affect the membrane potential (V M ). Additionally, expressing a channelrhodopsin-2 (ChR2) ion channel in the cellular membrane allows for light-induced depolarization. Depolarization of a neuron opens voltage-gated calcium channels (VGCC), causing an additional influx of Ca 2þ ions. In pyramidal neurons, Ca 2þ activates phospholipase A2 (PLA2) which converts phospholipids (PL) into arachidonic acid (AA). In turn, AA is metabolized into prostaglandin E2 (PGE 2 ) by cyclooxygenase-2 (COX-2). Pyramidal neurons release glutamate during depolarization, which activate GABAergic interneurons through NMDAR and AMPAR activation. GABAergic interneurons release GABA when depolarized, inhibiting surrounding cells. In specific subtypes of GABAergic interneurons, Ca 2þ ions activate nitric oxide synthase (NOS) which triggers the production and subsequent release of nitric oxide (NO). Other subtypes of GABAergic interneurons synthetize and express vesicle-bound vasoactive peptides such as neuropeptide Y (NPY). The release of these peptides is facilitated by Ca 2þ ions. Vascular smooth muscle (VSM) cells (brown rectangle) enwrap arterioles (red cylinder) and regulates the arteriolar diameter. PGE 2 are transported over cellular membranes by prostaglandin E2 transporters (PGET) and promotes arteriolar dilation by relaxing VSM by activation of the prostaglandin EP 4 receptor and upregulation of cyclic adenosine monophosphate (cAMP) synthesis. NO diffuse freely over cellular membranes, and acts to increase the production of cyclic guanosine monophosphate (cGMP), which in turn promotes relaxation of VSM and arteriolar dilation. Lastly, NPY activates the G-protein coupled NPY Y1 receptor (NPY1R) expressed on VSM cells, which inhibits synthesis and thereby promote arteriolar constriction. B: Each box represents a state, while non-boxes represent variables. Mass-action kinetics are shown with fully drawn arrows, while non-mass action kinetics are shown in dashed lines. u 1 ; u 2 ; u 3 denotes the input signals used in the experimental setup. NPY vsm and NO vsm represents the effect of NPY and NO on the vascular smooth muscle cells. For the model evaluation, we postulate that the effect of the anesthetic can be replicated solely by altering the neuronal dynamics and input signal (brown arrows).
including their alterations by an anesthetic agent.

Model formulation
All models are formulated using ordinary differential equations (ODEs). The general notation of an ODE model is as follows: in which x is the vector of state variables, describing the concentration or amount of various substances, with respect to time; _ x represents the derivative of the states with respect to time; f and g are non-linear smooth functions; p is the vector of unknown constant parameters, here kinetic rate constants or scaling parameters; where u is the input signal corresponding to experimental data; where xðt 0 Þ contains the values of the states at the initial time point t ¼ t 0 , which are dependent on the parameters x 0 ðpÞ; and where b y are the simulated model outputs corresponding to the measured experimental signal, here the arteriolar diameter change. Note that x, u, and b y depend on t, but that the notation is dropped unless the time-dependence needs to be especially stated, as in Eq. (1b).

Model structure
In Fig. 2A, we present an NVC pathway. This pathway was simplified in the model, with a corresponding model interaction graph depicted in Fig. 2B. The interaction graph describes all of the reactions and interactions of the model. The stimulus u is given by the following equation where t on , t off are the times when the signal goes on and off, respectively.

Presynaptic activity and calcium influx
Glutamate and GABA are released from different types of neurons and bind to specific ion channel-coupled receptors located in the plasma membrane of post-synaptic neurons. The primary excitatory neurotransmitter, glutamate, bind to NMDA and α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors. Activation of these ion channel-coupled receptors open ion-conducting pores leading to an influx of Na þ and Ca 2þ ions, causing a depolarization of the neuron. This in turn opens voltage-gated calcium channels, allowing additional influx of Ca 2þ ions. In contrast, GABA prevents depolarization of neurons by binding to the ion channel-coupled GABA A receptor, which opens Cl À ion-conducting pores. Several general anesthetics, including α-chloralose used in Uhlirova et al. (2016), act as positive allosteric modulators on the GABA A receptor, potentiating the Cl À ion conductance by increasing the affinity of GABA (Garrett and Gan, 1998).
Pyramidal neurons are glutamatergic, meaning that glutamate is released upon depolarization, and acts on e.g. astrocytes and interneurons. GABAergic interneurons release GABA upon depolarization and act on pyramidal neurons or other interneurons. This forms a canonical neuronal circuit, where interneurons play a key role in regulating the activity of pyramidal neurons. Glutamate acts to depolarize neurons and to increase the intracellular Ca 2þ levels, while GABA acts to prevent neurons to depolarize. In the model, we represent these interplays between pyramidal and inhibitory interneuron activity and their effect on respective neuronal Ca 2þ levels as Where N NO ; N NPY ; and N Pyr describe the neuronal activity of the GABAergic NO and NPY interneurons, and the pyramidal neurons, respectively; Ca 2þ i represents intracellular calcium in respective neuron; k u;i is the input strength of the stimulation to respective neuron; k PF1 and k PF2 are the kinetic rate parameters governing pyramidal to GABAergic interneuron signaling; k INF1 and k INF2 are the kinetic rate parameters describing the negative feedback from GABAergic interneurons to pyramidal neurons; k IN1 and k IN2 are the kinetic rate parameters describing the negative feedback between GABAergic interneurons; sink N;i are kinetic rate parameters governing the degradation of the activity of respective neuron; k Ca is the basal inflow rate of Ca 2þ (k Ca ¼ 10 for all three neurons); sink Ca;i is the elimination rate of intracellular Ca 2þ and i ¼ (Pyr, NO, NPY), denotes the three different neuron types. To ensure reasonable neuronal dynamics, we limited the values of some parameters: sink Ca;i ; sink N;i > 1:28 s À1 in line with experimental observations by (Majewska et al., 2000).

Pyramidal neuron signaling
The rise in intracellular Ca 2þ levels in pyramidal neurons activates phospholipase A2 (PLA2) which leads, after some enzymatic steps, to an increase in intracellular arachidonic acid (AA) (Lacroix et al., 2015;Davis et al., 2004;Iadecola, 2017). This is described by the following equation where k PL and k COX are kinetic rate parameters. In pyramidal cells, AA is metabolized into PGE 2 through a COX-2 and PGE synthase rate-limiting reaction (Lecrux et al., 2011;Lacroix et al., 2015). PGE 2 promotes vasodilation by relaxing VSM cells and pericytes through activation of EP2 and EP4 receptors (Lacroix et al., 2015;Davis et al., 2004). In the model, this process is represented as Where PGE 2; vsm represents PGE 2 acting on the VSM cells; k PGE2 and sink PGE2 are kinetic rate parameters.

GABAergic interneuron signaling
The rise in intracellular Ca 2þ levels in GABAergic interneurons provoke the release of different vasoactive messengers. For instance, NO has previously been shown to be a potent vasodilator at the level of arterioles in both cortical slices and in vivo (Cauli et al., 2004;Mishra et al., 2016;Rancillac et al., 2006;Taniguchi, 2014). NO is released by specific NO-interneurons through a Ca 2þ dependent nitric oxide synthase (NOS) rate limiting reaction. Here, we represent this process as Where NO vsm represents NO acting on the VSM cells; k NOS , k NO and sink NO are kinetic rate parameters. Here, sink NO > 1 s À1 following experimental work by (Vaughn et al., 1998). The vasoconstrictive messenger responsible for the post-stimulus undershoot was identified by Uhlirova et al. as NPY. NPY is released by specific NPY producing subtypes of GABAergic interneurons and has previously been shown to induce vasoconstriction of vessels in slice preparations (Cauli et al., 2004;Kocharyan et al., 2008) and in vivo (Uhlirova et al., 2016). NPY acts by binding to the NPY receptor Y1, a G αi -protein coupled receptor which is expressed in VSM cells enwrapping arteries and arterioles (Bao et al., 1997;Abounader et al., 1995Abounader et al., , 1999. Activation of this receptor inhibits the synthesis of adenosine monophosphate (cAMP) and increase the intracellular Ca 2þ through an inositol-1,4,5-trisphosphate receptor (IP3R) dependent release from the sarcoplasmic reticulum (Tan et al., 2018). These processes lead to VSM contraction and a subsequent arteriolar vasoconstriction. In the model, these mechanisms are represented as: Where NPY vsm represents NPY acting on the VSM cells; k NPY and sink NPY are kinetic rate parameters. The conversion between intracellular NPY and NPY in the VSM, NPY vsm , is governed by Michaelis-Menten kinetics (see Eq.7a-b, v1) (Michaelis and Menten, 1913;Johnson and Goody, 2011), where K M is the Michaelis constant (the concentration of the substrate at which half max reaction rate is achieved) and V max is the maximal reaction rate. Michaelis-Menten kinetics was assumed for this reaction because initial tests with simpler mass-action kinetics could not produce a satisfactory agreement to the experimental data. The interpretation of Michaelis-Menten kinetics is that when ½NPY ≪ K M the reaction behaves like a first-order mass-action reaction, meaning that the reaction is linearly dependent on ½NPY (V max NPY). However, as ½NPY increases beyond K M , ½NPY ≫ K M , the reaction behaves like a zero-order reaction, meaning that the reaction is independent of ½NPY, and only depends on the maximal reaction rate V max . In between these two limits, ½NPY has a saturated effect on the reaction rate and when ½NPY equals K M , the reaction rate has reached half the maximum value.

Arteriolar diameter change
We express the arteriolar diameter as the per cent normalized weighted summation of the end states of the three vasoactive signaling pathways: Where k y1 , k y2 , k y3 are unknown weighting factors which determines the influence of each respective vasoactive pathway on the smooth muscle cells surrounding the arterioles.

Model evaluation
Once the model structure has been formulated ( Fig. 2B) and data has been collected, the parameters, p, need to be determined. This is commonly done by evaluating the negative log-likelihood function. Assuming independent, normally distributed additive measurement noise, the likelihood of observing data D given p is:

Experimental data
All experimental data were obtained from the study of Uhlirova et al. (2016), and no new data were acquired in our study. Specifically, arteriolar diameter responses measured using two-photon imaging were extracted from Figs. 5 and 6 in Uhlirova et al. (2016) using the online tool WebPlotDigitizer (Rohatgi). The original article state that all experimental procedures were performed in accordance with UCSD Institutional Animal Care and Use Committee guidelines. In brief, four different stimulation protocols were used in their study: three conducted under the administration of the anesthetic α-chloralose, and the fourth during awake condition. The time series from the original manuscript were shown with error bars depicting the standard error of the mean (SEM). Below, condensed descriptions of each stimulation protocol are presented, and we kindly refer the reader to the original publication for a complete detailed description. For the model training, we used all control conditions specified below, and for model validation we used the conditions with pharmacological applications.

OG stimulation of pyramidal neurons during anesthesia
Thy1-ChR2-YFP mice (n ¼ 2 subjects), with ChR2 selectively expressed in layer 5 pyramidal cells, were imaged (28 measurement locations along 13 arterioles within 40-380 μm range) during OG stimulation (single light pulse lasting 50-80 ms) under both control conditions (Fig. 3J, black error bars) and under the influence of both the AMPA/ kaniate receptor antagonist 6-cyano-7-nitroquinoxaline-2,3-dione (CNQX, 500 μM) and the NMDA receptor selective antagonist D-(-)-2-Amino-5-phosphonopentanoic acid (AP5, 200 μM) (Fig. 4F, black error bars). The responses were extracted from Fig. 5A in (Uhlirova et al., 2016) by sampling all error bars, resulting in a time resolution of 1 s between each data point. In the model, we simplified the stimulus by assuming a fixed length of 100 ms affecting the pyramidal neuron with a constant value (i.e. a continuous box-car function). The stimulus amplitude, k u; Pyr , was treated as a free parameter and was optimized to fit data.

OG stimulation of GABAergic interneurons during anesthesia
VGAT-ChR2(H134R)-EYFP mice (n ¼ 3 subjects), with ChR2 selectively expressed in GABA interneurons, were imaged (52 measurement locations along 25 arterioles within 50-590 μm range) during OG stimulation (one pair of light pulses separated by 130 ms, for a total of 450 ms duration) under both control conditions (Fig. 3H, black error bars) and during application of the neuropeptide Y receptor Y1 inhibitor BIBP-3226 (100 μM) (Fig. 4D, black error bars). The responses were extracted from Fig. 5B in (Uhlirova et al., 2016) by sampling all error bars, resulting in a time resolution of 1 s between each data point. In the model, we assumed a constant stimulus value for the whole duration of 450 ms, scaled by the amplitude factors k u;NPY and k u;NO for respective neuron which were treated as free parameters.

Sensory stimulation during anesthesia
Additionally, 4 wild-type mice underwent imaging ( In the model, we assumed a constant stimulus value for the whole duration of 2 s, scaled by k u;NPY , k u;NO and k u; Pyr affecting each respective neuron, and were treated as free parameters.

Sensory and OG stimulation during awake condition
Lastly, 4 wild-type and 3 VGAT-ChR2(H134R)-EYFP awake mice were imaged during sensory stimulation (1 s of air puffs, 100 ms puffs at 3 Hz) and OG stimulation (single light pulse lasting 150-400 ms), respectively ( Fig. 3B and D, black error bars). The responses were extracted from Fig. 6A and B in (Uhlirova et al., 2016) by sampling all error bars, resulting in a time resolution of 0.7 s between each data point. In the model, we assumed a constant stimulus value for the whole duration scaled by the amplitude parameters k u;i , with a length of 1 s and 400 ms for sensory and OG stimulation respectively.
Before being used in the model evaluation, SEM in all extracted time series with SEM measurements smaller than mean SEM was set to the mean SEM. This was done to avoid overfitting the model to a few extreme data points. Fig. 3. Model evaluation to arteriolar response data from awake (A-D) and anesthetized animals (E-J) for both optogenetic (OG) (C-D, G-H, I-J) and sensory (A-B, E-F) stimuli. A, C, E, G; I: Simple schematics illustrating the experimental design for each data set. Grey triangle: pyramidal neuron; yellow circle: GABAergic interneurons. u 1 ; u 2 ; u 3 indicate the active inputs for each experimental setup (see Fig. 2B). B, D, F, H, J: Model estimation paired with experimental data. For each graph, best estimated model simulation (solid red line) paired with model uncertainty (red shaded areas) compared to experimental data (black symbols, error bars depicting standard error of the mean). The stimulation length is indicated by the black bar in the lower left portion of each graph.

Model hypothesis formulation
We created a new model structure specified by the equations in Section 2.2. In this model, the effects of the anesthetic appear solely for the rates affecting the neuronal dynamics (N i ; i 2 fNO; NPY; Pyrg , see Eq. 3A-C). More specifically, we allow 12 parameters to change between awake and anesthetized animals: the three stimulus scaling parameters k u;i (Fig. 2B brown arrows); the three decay rate parameters sink N;i affecting N i (Fig. 2B, solid brown arrows) and finally the six neuronal interaction parameters k PF1 ; k PF2 ; k IN1 ; k IN2 ; k INF1 , and k INF2 (Fig. 2B, dashed brown arrows). Furthermore, we allow the stimulus scaling parameters k u;i to change between OG and sensory stimulation. Further details on the model, estimation of parameter values, model analysis, and all simulation files to reproduce the results are available as Supplementary Materials.

The model can describe all experimental data used for model training, both qualitatively and quantitatively
The model was trained on experimental data consisting of all control responses from the Uhlirova et al. study (Uhlirova et al., 2016) (Fig. 3  right panel (B, D, F, H, J), black symbols). In Fig. 3 left panel (A, C, E, G,   Fig. 4. Model predictions to drug perturbed arteriolar response data during anesthesia condition. Left: Simple schematics illustrating the experimental design for each data set. Yellow circle: GABAergic interneurons; Orange triangle: pyramidal neuron. u 1 ; u 2 ; u 3 indicate the active inputs for each experimental setup (see Fig. 2B). Right: Model predictions and experimental data. For each graph, model prediction uncertainty (red shaded areas) compared to experimental data (black symbols, error bars depicting standard error of the mean). The stimulus length is indicated by the black bar in the lower left portion of each graph. A & B: OG stimulation (400 ms) of GABAergic interneurons (B) and sensory stimulation (1 s) (C) during the presence of the NPY receptor Y1 antagonist BIBP. C: OG stimulation (100 ms) of pyramidal neurons in the presence of the glutamatergic signaling blockers AP5 and CNQX.

Fig. 5. Michaelis-Menten kinetics of NPY release
can explain prolonged undershoot seen in experimental data. A: simple schematic of the NPY signaling arm governed by the Michaelis-Menten reaction v1 described in Eq. 7a-b. Note that we in the following simulations use the full model described in section 2.2 and Fig. 2B. B: simulation of ½NPY for GABAergic OG stimulation for awake (red solid line) and anesthetized (black solid line) animals. The difference between the responses are highlighted by the hatched area. C: reaction rate v 1 (y-axis, normalized to the maximal reaction rate V max , black solid line) as a function of increasing ½NPY (x-axis). The achieved simulated maximal reaction rate for the different stimulus paradigms is shown in red (awake condition) and black (anesthesia condition) symbols. D: simulation of ½NPY vsm for GABAergic OG stimulation. For B & D, the change in respective state is shown on the y-axis, and time in seconds on the x-axis. I), a simple schematic describes the experimental setup. The model parameters were fitted to the experimental data, achieving a qualitative and quantitatively acceptable fit (Fig. 3 right panel, solid thick red lines), evaluated using a χ 2 -test (J lsq ðb p H1 Þ ¼ 42.4, cut-off: χ 2 (190 data points) ¼ 223.16 for α ¼ 0.05). Next, the uncertainty of the model was approximated using a Markov Chain Monte Carlo (MCMC) sampling algorithm. Posterior distributions of the parameters were sampled using 10 5 iterations. These distributions can be found in the supplementary material ( Fig. S1) along with the point-wise parameter estimate. In addition, all χ 2 acceptable parameters (J lsq ðpÞ < 223:16) found during the MCMC sampling were saved. Using the generated Markov chains and acceptable χ 2 parameters, a confidence interval: was drawn where the threshold Δ α ðχ 2 df Þ is the α quantile of the χ 2 df statistics (Kreutz et al., 2013;Maiwald et al., 2016). The significance was set to α ¼ 0:05 and degrees of freedom df H1 ¼ 46, equal to the number of estimated parameters. The resulting uncertainty is depicted in the form of shaded red areas in Fig. 3, right panel.
As seen in Fig. 3 (right panel), the best model fits (solid red lines) lie within experimental uncertainty (black symbols) in all five experiments. In addition, all quantitatively acceptable model simulations (red shaded areas) follows the same dynamic behavior as seen in the data. These results mean that the model is capable of describing experimental data used for model training, both qualitatively and quantitatively.  Fig. 2B for the complete model used for these simulations). The following states are displayed: the neuronal dynamics N NPY (A) and N NO (B); the calcium dynamics ½Ca NPY (C) and ½Ca NO (D); the first intracellular signaling state ½NPY(E) and ½NO(F); the vasoactive state of respective neuron, ½NPY vsm (G) and [NO vsm (H), and finally, the resulting vascular response (I). For each of these states, two graphs are shown depicting the dynamics for awake (right (left graph, black solid line) and anesthetized (right graph, red solid line) animals. In the middle, an interaction graph showing how these states are connected, and which reaction that are altered with anesthesia (brown arrows) is shown. All time courses shown in the graphs are normalized to their baseline value for easier comparisons of the dynamics and amplitude between awake and anesthetized animals.

The model can predict experimental data not used for training
To further evaluate the model, we used the previously generated Markov chains (see Section 3.1) to generate model predictions of experimental data not used for training. More specifically, all acceptable parameter sets within the confidence interval were used to predict the model response to a) sensory stimulation in the presence of NPY receptor Y1 antagonist BIBP (Fig. 4A and B); b) OG stimulation of GABAergic interneurons in the presence of the NPY receptor Y1 antagonist BIBP ( Fig. 4C and D), and, c) OG stimulation of pyramidal neurons during inhibition of glutamatergic signaling by administration of the AMPA and NMDA antagonists CNQX and AP5, respectively ( Fig. 4E and F). The corresponding experimental data used for these three pharmacological perturbations were extracted from Fig. 5 in Uhlirova et al. (2016). In the model, these inhibitions were carried out by: i) removing the transient effect of the NPY signaling arm on the arteriolar diameter change by setting S ¼ k y1 NO vsm ðtÞ þ k y2 PGE 2;vsm ðtÞ À k y3 NPY vsm ð0Þ (Fig. 4A-D), and, ii) removing feedback from pyramidal neurons to the GABAergic interneurons by setting k PF1 ¼ k PF2 ¼ 0, assuming a complete blockade of glutamatergic signaling (Fig. 4E and F). Fig. 4 shows a simple illustration of each experimental setup (left panel, A, C, E) and model prediction uncertainties (B, D, F: red shaded areas) together with corresponding experimental data (black symbols). The model predictions follow the same qualitative behavior as the experimental data, overlapping the majority of all error bars. For Fig. 4A-D, we assume a complete inhibition of the transient NPY effect on arteriolar diameter, removing the possibility of negative deflection in arteriolar diameter change and thereby truncating the simulated responses to be positive for all timepoints (grey shaded areas). Finally, there are many of the individual curves within the uncertainty area that passes a χ 2 -test, meaning that those predictions are also quantitatively correct.

The Michaelis-Menten kinetics of NPY signaling removes amplitude differences but preserves duration differences
Having established validity of the model, we now examine mechanistic properties of the system. The most prominent characteristic of the experimental data is the prolonged post-peak undershoot evoked by OG stimulation in anesthetized animals ( Fig. 1C (black symbols); Fig. 3H-J) as compared to awake animals ( Fig. 1C (red symbols); Fig. 3D). The key feature in the model producing this difference is the saturated signal propagation in the NPY signaling pathway, possible via the Michaelis-Menten expression (v 1 in Fig. 5A; see Eq. 7a-b, v 1 ). Upon OG stimulation, The increase of ½NPY is larger in anesthetized animals than awake (Fig. 5B, compare black and lines). This difference in ½NPY response disappears in the next reaction, which is the conversion of NPY to NPY vsm (governed by v 1 ), because ½NPY is above K M , meaning that the reaction rate of v 1 has started to saturate (Fig. 5C). In other words, the increase from awake to anesthetized animals in ½NPY is large (>300% bigger; Fig. 5B), but the corresponding increase in v 1 is small (Fig. 5C, less than 20% bigger). This saturation implies that the dual difference between awake and anesthetized animals in ½NPY (Fig. 5B, higher amplitude and longer duration) is transformed to a basically single difference in ½NPY vsm (Fig. 5D, only longer duration).
This transformation is the key property in the NPY signaling leading to the prolonged post-peak undershoot when considered as a part of the entire system.

The model explains the prolonged post-peak undershoot in anesthesia conditions by a prolonged NPY signaling
To fully understand how the model explains the anesthesia-induced change in vascular dynamics, let us now consider the case of GABAergic OG stimulation for awake and anesthetized animals for the entire model. Fig. 6 shows the dynamics of the NPY and NO neurons, from the neuronal signaling to the release of their respective vasoactive agents. In the middle of the figure, a cut-out of the full interaction graph (Fig. 2B) is shown, highlighting the reactions changed between awake and anesthetized animals (Fig. 6, brown arrows).
Let us first follow the signal propagation in the left NPY branch (Fig. 6A, C, E, G). First, the anesthesia-induced changes (brown arrows) implies a difference in amplitude for the rapid neuronal response (Fig. 6A, N NPY Þ: This difference is then transformed to a difference in both amplitude and duration for the much longer ½NPY response (Fig. 6E).
Let us now follow the signal propagation in the right NO branch (Fig. 6B, D, F, H). As can be seen, there is here a smaller difference in amplitude, which stays relatively small throughout the entire signaling cascade.
Let us finally consider how these two branches combine to form the final response (Fig. 6I). The combined effect of the two branches is obtained by taking the difference between them: ½NPY vsm is vasoconstrictive and ½NO vsm is vasodilating. Since [NO vsm is slightly lower in anesthesia, but with a preserved timing, and [NPY vsm has a longer duration but roughly the same amplitude, the diameter dynamics has an approximately the same peak amplitude but a longer post-peak undershoot.
This simple analysis illustrates how straightforward model simulations can help unravel the mechanism by which our validated model explains the various facets of data.

Discussion
A model-based approach is necessary to fully understand the complexity of the NVC, but no previously existing model describes neither the relevant neuronal subtypes nor the effects of anesthetics (Fig. 1A). Herein, we provide a first systems biology model for the crosstalk between different GABAergic interneurons and pyramidal neurons. The model can quantitatively explain experimental data used for model training (Fig. 3) and correctly predict experimental data not used for training (Fig. 4), both for awake ( Fig. 3A-D) and anesthetized ( Fig. 3E-J; Fig. 4A-F) animals. The analysis shows that the prolonged post-peak undershoot in anesthesia can be explained solely by excitability changes since amplitude changes can be translated to downstream signal duration changes via a Michaelis-Menten expression (Figs. 5 and 6).

Previous models of the NVC
Numerous earlier studies have done mathematical modeling of arteriolar diameter changes e.g., convection-diffusion models based on: 1) transient blood oxygen level dependent (BOLD) signals (Kim et al., 2013;Kim and Ress, 2016;Zheng et al., 2005), 2) two-photon optical measurements (Huppert et al., 2007(Huppert et al., , 2009 or, 3) direct measurement of arteriolar diameter in mice (Barrett et al., 2012). Recently, spatio-temporal models such as vascular anatomical network (VAN) models (Fang et al., 2008;Boas et al., 2008;Gould et al., 2017) were shown capable of simulating the entire vascular network using two-photon imaging data and associated BOLD responses (Gagnon et al., 2015). However, these models do not include any biochemical mechanisms governing the intracellular processes underlying vasodilation and vasoconstriction. Some models exist which include such processes (Yucel et al., 2009;Mathias et al., 2016Mathias et al., , 2018Dormanns et al., 2015;Blanchard et al., 2016), but none of these models describe the molecular mechanism of NPY elucidated in Uhlirova et al. (2016); the difference between sensory and OG stimulation; the crosstalk between GABAergic interneurons and pyramidal neurons; or the effects of anesthetics.

The role of anesthetics in modulating the NVC
Previous experimental research has characterized different anesthetics and their effects on evoked hemodynamic responses, for instance by looking at changes in response dynamics. Multiple studies have shown that 1) the amplitude of the evoked hemodynamic response (BOLD, CBV, CBF) are suppressed during anesthesia (Peeters et al., 2001;Martin et al., 2006;Goense and Logothetis, 2008;Pisauro et al., 2013;Aksenov et al., 2015;Drew et al., 2011); 2) in some cases, the observed hemodynamic response is slightly delayed compared to awake condition (Martin et al., 2006;Pisauro et al., 2013), and, 3) anesthetics induce changes in large scale neuronal networks, as seen in functional connectivity (Williams et al., 2010;Paasonen et al., 2018). The used anesthetic in Uhlirova et al. (2016), α-chloralose, is commonly used in animal research on the NVC because it preserves a robust hemodynamic response to sensory stimulation (Lindauer et al., 1993;Ueki et al., 1992). Garrett and Gan (1998) has reported that α-chloralose potentiates the GABA-induced GABA A current by increasing the affinity of GABA, thereby acting as a positive allosteric modulator on the GABA A receptor complex. Other studies have shown that α-chloralose reduces baseline CBF values (Masamoto and Kanno, 2012), increases blood pH, and decreases heart rate (Low et al., 2016), which could affect the measured BOLD responses.
In the used experimental data (Fig 1C; Fig. 3D, H, J), a clear difference can be seen between awake and anesthetized animals. For OG stimulation, the response during the awake condition (Fig. 3D) exhibits a faster post-peak undershoot phase compared to the anesthesia condition (Fig. 3H). For the sensory stimulation, the post-peak undershoot is reduced in awake animals (Fig. 3B, F). These changes are driven either by changes in neural processing or changes in vascular reactivity, or both. In the model, we have assumed that the anesthetic only affects the neuronal processing, allowing the parameters which governs the rapid neuronal dynamics to change between awake and anesthetized animals (see section 3.1). This assumption is sufficient to explain the used experimental data. This is possible because amplitude changes on the GABAergic interneuron activities (N NPY ) is translated into first amplitude and duration changes (for [NPY) and then into only duration changes for v 1 and [NPY vsm via the Michaelis-Menten expression for v 1 (Fig. 5; Fig. 6A, C, E, G). Since the GABAergic NO interneurons only display a slightly smaller amplitude (Fig. 6B, D, F, H), the effect of adding the two GABAergic interneuron contributions together is a slightly decreased dilation and a clearly prolonged post-peak undershoot (Fig. 6I). Importantly, our analysis does not rule out that the observed effects of the anesthetic could be due direct effects of anesthetics on the arteriolar VSM or other vascular segments, since this possibility was not investigated with the model. However, a recent study by Vanlandewijck et al. (2018) reports that arteriolar VSM lack GABA A receptors in mice.

Other mechanisms facilitating the neurovascular coupling
The astrocytic contributions to the NVC in vivo is still unclear. The end-feet of astrocytes envelope cortical blood vessels, making them a potential mediator of the NVC. Indeed, previous studies have shown that astrocytes respond to neuronal activity via activation of glutamate receptors mGluR 1 and mGluR 5 , triggering a rise in intracellular Ca 2þ concentration through an inositol-1,4,5-trisphosphate receptor (IP3R) mediated mechanism and culminates in the release of vasoactive gliotransmitters such as prostaglandins, epoxyeicosatrienoic acids, or K þ ions (Bazargani and Attwell, 2016). However, this explanation has recently been challenged, as Sun et al. (2013) reported that adult human and mouse astrocytes lack expression of mGluR 1 and mGluR 5 and other studies have reported that IP3R knockout mice retain unaltered NVC in vivo (Bonder and McCarthy, 2014). A recent study by Mishra et al. (2016) however, provided evidence for a new mechanism, where exogenous release of ATP from post-synaptic neurons triggers astrocytic Ca 2þ influx through the ATP receptor P2X 1 with a subsequent capillary dilation as a result of cyclooxygenase-1 mediated PGE 2 release.
A recent study by Longden et al. (2017) shows that extracellular K þ released by neuronal activity opens inward-rectifier K þ (K IR 2.1) channels in capillary endothelial cells, inducing retrograde hyperpolarization which rapidly propagates to upstream arterioles, causing arteriolar dilation. Following this, Rungta et al. (2018) report a complete spatio-temporal vascular response from juxta-capillaries back to the pial arteriole. They show that local synaptic activation triggers a synchronous Ca 2þ drop in mural cells along the upstream vascular tree. However, the parenchymal and first order arterioles dilate rapidly, preceding capillary dilation. Together, these studies show that the hemodynamic response can be initiated at the capillary level, spreading quickly through electrical signals propagating upstream the vascular tree, but arterioles form the primary unit for the onset of vascular responses. In the model, we do not account for the mechanism of retrograde hyperpolarization, and instead describe the onset of arteriolar dilation as dependent on local vasoactive messengers.

Sensitivity analysis versus prediction and parameter uncertainties
A sensitivity analysis is typically done to examine the result of local variations around a single point in the parameter space. Since we do not have a uniquely determined parameter point, but rather a long list of acceptable parameter values, such a point-wise local analysis is of limited value. We therefore believe that a global approach, which considers all acceptable parameters and basically ignores non-acceptable parameters, is the most appropriate analysis. MCMC is such a global approach. We use MCMC both to provide an estimate of model uncertainty for predictions and the uncertainty of each individual parameter. This approximates a sensitivity analysis, as it quantifies the impact of parameter changes on predictions: this impact is seen by the shaded regions in connection to each model simulation (for instance, Fig. 3). These uncertainty ranges should be considered in combination with the uncertainty of each individual parameter, which is given by the parameter distributions in Figure S1. This analysis does not exactly translate to a parameter specific sensitivity analysis, but rather provides an overall assessment of the relationship between the uncertainties in parameters and predictions considered as a whole. Nevertheless, given a specific parameter, a wide distribution implies that the parameter has a low sensitivity, and a narrow distribution implies that the parameter has a high sensitivity (when allowing for compensations of the other parameters to maintain a good agreement with data).

Assumptions and limitations
The present study has a number of limitations. The model structure is built upon underlying assumptions limiting the strength and scope of the conclusions that can be drawn. Nevertheless, such assumptions are necessary in order to develop a comprehensible model. In the model, we assume a simple relationship between GABAergic interneurons and pyramidal cells, where pyramidal neurons excite the GABAergic interneurons, which in turn act to inhibit the pyramidal neuron in a phenomenological manner (Eq. (3)). Here, a possible extension would be to include Hodgkin-Huxley description of ion-channel kinetics (Hodgkin and Huxley, 1952), describing the activity of each respective neuron as fluctuations of the membrane potential arising from opening and closing of ion channels due to glutamatergic and GABAergic signaling. Furthermore, the observed effect of disinhibition, where GABAergic interneurons are inhibited leading to excitation is not included. Other parts of the model structure lacking physiological details are the action of VSM cells and vessel elasticity of the arterioles. These mechanisms are simply represented by intermediary states (e.g., Eq. (6b)) in our model, which do not accurately capture the mechanisms of smooth muscle contraction or relaxation. In addition, our model does not include spatial dimensions or structural information regarding the cortical vessels, captured in other models (Boas et al., 2008;Aquino et al., 2014). Lastly, throughout the model, the amount of each state is specified in arbitrary units, which means that we describe biological values such as concentrations, scaled by unknown scaling constants. Note that this implies that the model still can be used to predict and simulate key observations that are possible to observe in experimental data, for example, the response time, relative amplitudes between different stimulation strengths, whether there is a post-peak undershoot, etc.

Conclusions
We have provided a first quantitative systems-level understanding of GABAergic interneurons and pyramidal neurons in the NVC, in awake and anesthetized animals, in response to both sensory and OG stimulations. The model can qualitatively predict the effect of different pharmacological inhibitors (Fig. 4), and explain a key characteristic seen in experimental data: the post-peak undershoot is prolonged during anesthesia (Figs. 5 and 6). This explanation shows how our quantitative systems-level model can be used to understand the mechanistic underpinnings of complicated experimental observations in the NVC.

Funding
(The journal requires that funding is visible on the cover page during the review process).
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: SS and ME received funding from the Swedish Research Council