Parameter Sensitivity and Experimental Validation for Fractional-Order Dynamical Modeling of Neurovascular Coupling

Goal: Modeling neurovascular coupling is very important to understand brain functions, yet challenging due to the complexity of the involved phenomena. An alternative approach was recently proposed where the framework of fractional-order modeling is employed to characterize the complex phenomena underlying the neurovascular. Due to its nonlocal property, a fractional derivative is suitable for modeling delayed and power-law phenomena. Methods: In this study, we analyze and validate a fractional-order model, which characterizes the neurovascular coupling mechanism. To show the added value of the fractional-order parameters of the proposed model, we perform a parameter sensitivity analysis of the fractional model compared to its integer counterpart. Moreover, the model was validated using neural activity-CBF data related to both event and block design experiments that were acquired using electrophysiology and laser Doppler flowmetry recordings, respectively. Results: The validation results show the aptitude and flexibility of the fractional-order paradigm in fitting a more comprehensive range of well-shaped CBF response behaviors while maintaining a low model complexity. Comparison with the standard integer-order models shows the added value of the fractional-order parameters in capturing various key determinants of the cerebral hemody-namic response, e.g., post-stimulus undershoot. This investigation authenticates the ability and adaptability of the fractional-order framework to characterize a wider range of well-shaped cerebral blood flow responses while preserving low model complexity through a series of unconstrained and constrained optimizations. Conclusions: The analysis of the proposed fractional-order model demonstrates that the proposed framework yields a powerful tool for a flexible characterization of the neurovascular coupling mechanism.

Index Terms-Neurovascular coupling, Cerebral blood flow, Neural activity, Fractional-order calculus, fractional differentiation orders, Sensitivity analysis.
Impact Statement-The present study proposes a novel fractional-order framework for modeling neurovascular coupling. A parameter sensitivity analysis demonstrates the potential flexibility, and effectiveness of the fractionalorder paradigm in reconstructing the cerebral hemodynamics with manageable complexity; and a real experimental validation analysis demonstrates the ability of the model in modeling a wider range of well-shaped CBF responses.

I. INTRODUCTION
C HANGES in neural activity lead to changes in local cerebral blood flow (CBF) and energy metabolism. The complex relationship between neural activity and cerebral blood flow, referred to as neurovascular coupling (NVC), is a subject of intensive investigation. Understanding factors and mechanisms that orchestrate this relationship will improve our understanding of the physiological underpinnings of measurements from Functional Magnetic Resonance Imaging (fMRI) [1]. Investigating NVC in humans has become a possibility with the development of neuroimaging techniques that measure local hemodynamics, including CBF. Thus, quantitative models have been proposed to describe the different mechanisms linking transient neural activity to the changes in CBF [1]- [3]. NVC models can be arranged in a broad spectrum ranging from simple models (with fewer details and fewer parameters) to more complicated models involving many biophysical parameters and complex relationships [4]. An example of such models includes the simple model developed in Friston et al. [5], [6] and which takes a stimulus waveform and outputs a CBF response shape. The initial goal of developing this model was to fill in the neurovascular compartment to the well-known Balloon Model to predict Blood Oxygen Level Dependent (BOLD) (measured with functional Magnetic Resonance Imaging (fMRI)) This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ VOLUME 3, 2022 69 responses given neural activity. A second order differential equation is used to describe CBF changes given an input stimulus. Another example is the model developed in Buxton et al. [7] which consists of a simple neural adaptation model and a linear convolution of neural activity with a flow response function.
While exploring biological processes underlying NVC using more complicated models is desired, their applicability is limited due to several mathematical constraints. For example, accurate and reliable estimation of model parameters is more difficult in complex models due to the higher number of parameters. Furthermore, parameters of those models are harder to interpret because of the lower model identifiability in complex models [4]. The difficulty of parameter estimation combined with the low identifiability and interpretability may encourage using simpler models with fewer parameters. Simpler NVC models, however, remain limited in their flexibility to fit the dynamics of CBF responses derived from specific experimental conditions or pathological populations. Deneux et al. [8], for example, showed that two of the dynamic linear NVC models (namely Friston's model and Buxton's model) along with their nonlinear variations were unable to capture different dynamics of CBF responses at various stimulation lengths. Linear variations appear to fail to fit the amplitude variations of different responses. While they underfit responses to short stimulations, they overfit responses to longer stimulations. Nonlinear variations also fail to account for temporal dynamics of the CBF responses especially those for shorter stimulations [8]. Those results call for developing more flexible models that can fit experimental data without compromising model simplicity.
In the last decades, non-integer differentiation, the so-called fractional-order differential calculus, became a popular tool for characterizing real-world physical systems and complex behaviors from various fields such as biology, control, electronics, and economics [9]- [13]. The long-memory and spatial dependence phenomena inherent to the fractional-order systems present unique and attractive feature which raises exciting opportunities. These opportunities include an accurate representation of power-law behavioral phenomena, [14], [15]. For instance, the power-law behavior has been demonstrated in describing human soft tissues viscoelasticity and characterizing the elastic vascular arteries. In-vivo and in-vitro experimental studies have pointed that fractional-order calculus-based approaches are more decent to precisely represent the hemodynamic; the viscoelasticity properties of soft collagenous tissues in the vascular bed; the aortic blood flow; red blood cell (RBC) membrane mechanical properties and the heart valve cusp [16]- [19]. Consistent with the ability of fractional-order models to fit temporal dynamics, Belkhatir et al. in [20] have shown that the fractionalorder model can yield better fit to Blood-Oxygenation-Level-Dependant (BOLD) signals measured with fMRI when compared with the original integer-order NVC model proposed by Friston et al. in [5]. Fractional calculus has been used as a powerful tool to understand better the dynamic processes that span spatiotemporal scales. Essentially, the fractional-order model is a continuous-time model with high flexibility to fit high-order dynamics and complex nonlinear phenomena. Based on this study, fractional calculus seems to be a suitable approach for NVC modeling. In fact, one of the most important properties of fractional-order derivatives is that they depend on the entire history of a function, not only the value of the function at the evaluated point. This property, called non-locality or memory effect, is relevant for modeling systems that exhibit temporal dynamics and delays such as CBF response since the model's response at any given time depends on the whole history of the CBF response.
The presented work aims to extend the previous work [20] by studying the parameter sensitivity analysis of the fractionalorder dynamical model of NVC and validated this model using both synthetic and real CBF data. The contributions of this paper can be summarised in two main parts: r In the first part, the mathematical model is illustrated through a series of numerical simulations, which demonstrates the fractional-order model's ability to fit a wider range of well-shaped CBF responses that cannot be captured with the standard models. In addition, using an extensive parameter sensitivity analysis, we study the effect of the fractional differentiation order and the model's parameter on the CBF response.
r In the second part, the proposed model has been applied and validated using experimental CBF data obtained from [21]. Moreover, we evaluated the performance of the model and compared it to the integer-order model. This paper is organized as follows. In Section II, we will recall some basic concepts from the fractional-order derivatives and the fractional-order neurovascular coupling mode.In addition, this section presents the parameter sensitivity analysis of the NVC fractional-order model along with the adopted method to fit the real CBF data. Section III discusses the obtained results and provides some future directions on the use of the model for analyzing the cerebral hemodynamic.

A. Background 1) Fractional-Order Calculus:
The concept of fractional calculus is very old and goes back to the seventeenth century. Fractional calculus is defined as a generalization of the integerorder integration and differentiation operators to the non-integer order. Because of its interesting properties of non-locality and memory, the interest on fractional derivatives has grown in many fields of engineering and science. Examples of real-life applications include but not restricted to: viscoelastic, diffusive, biomedical and biological systems [22]- [27].
The continuous fractional integro-differential operator D α t , where α and t are the limits of the operation, is defined as follows For fractional derivative, several definitions exist in the literature [28], [30]. In this work, we consider the generalized Riemann-Liouville (RL) definition which is proposed in [29] and recalled in definition 1. This definition is more appropriate for mathematical analysis. 1) Definition 1:: Equitation [30] The initialized RL fractional derivative of order α ∈ (0, 1) of a function g, denoted 0 D α t g(t), is given by: where 0 d α t g(t) and Ψ(g i , α, −a, 0, t) are the uninitialized α th order RL derivative and the initialization function, respectively. They are given as follows: For numerical implementation the definition of Grunwald-Letnikov (GL) given in definition 2 is used [28]. 1) Definition 2:: Equation [28] The Grunwald-Letnikov derivative of order α of a function g, denoted D α t g(t), is given by: where h > 0 is the time step, c are the binomial coefficients recursively computed using the following formula,

2) Fractional Neurovascular Coupling Model:
The proposed fractional neurovascular coupling can be formally written as: where f is the CBF, is the neural efficacy, k s is signal decay and k f is the feedback term, q 1 and q 2 are the fractional differentiation orders which range between 0 and 1. Note that when both fractional orders are set to 1, the model represents the original model, which we refer to as integer-order model, proposed by Friston et al. [5]. Fractional dynamics are only present when any fractionalorder is set to a value less than 1. We refer to it as fractional-order model. This note is the reason why the integer-order model is a special case of the fractional-order model where both fractional parameters are set to a value of one. According to the integer-order model (where q 1 = q 2 = 1 in (6)), an increase in neural activity u(t) leads to an increase in the flow-inducing signal s (which is assumed to control CBF, f , at the arteriole level). The flow inducing signal s, then, leads to an increase in CBF, f . This integer-order NVC model is simple. However, it fails to account for more complex effects that arise from the fractional dynamics underlying the CBF response [8]. Hence, two new fractional differentiation order (namely, q 1 and q 2 ) are introduced to fully model and account for the fractional properties of the CBF responses.

B. Characterizing the Unique Contribution of the Fractional Parameters
In the first analysis, we ask whether the fractional-order model can generate unique well-shaped CBF responses that cannot be produced with integer-order models, how those contributions change the shape of the CBF response. More formally, this analysis aims mainly to exclude any equivalence that may exist between the integer-order model and the fractional-order model. In other words, we ask whether we can match any output of the fractional-order model (that has two more parameters) by solely tuning parameters of the integer-order model? If the two models are equivalent, then there exists a set of values for k f and k s in the integer-order model that can match any output of the fractional-order model. We can directly investigate this claim by running an optimization problem that minimizes signal dissimilarity between the two models. The minimization only optimizes the parameters k f and k s of both models, by fixing q 1 and q 2 to constant values lying between 0 and 1.
By fixing q 1 and q 2 of fractional differentiation order model at some constant(s) ranging between 0 and 1, we can compare the obtained optimization results (i.e. the value of the function at optimal values) across the range of q 1 and q 2 . If the two models are essentially equivalent, then signal dissimilarity should not change across all values used for q 1 and q 2 . However, a change in signal dissimilarity across the range of q 1 and q 2 indicates non-equivalence.
The cost function we used to calculate signal dissimilarity is the L 1 -norm of the difference between the two signals. It can be formally described as follows: where f q 1 <1,q 2 <1 denotes the CBF computed using the fractional-order model and f q 1 =1,q 2 =1 denotes the CBF computed using the integer-order model. Each model has a separate set of four parameters. The fractional-order model has: q 1 (where q 1 <1), q 2 (where q 2 <1), k f and k s . Similarly, the integer-order model has: q 1 (where q 1 = 1), q 2 (where q 2 = 1), k f and k s . All k f , k s , k f and k s are free to vary while q 1 and q 2 are fixed at some predetermined values. For the integer-order model, q 1 and q 2 are always fixed at a value of 1 as have been noted before. However, in the fractional-order model, we use different combinations of values to test for the contribution of fractional differentiation parameters. The goal, then, is to find some combination of the four parameters that best minimize their dissimilarity, under some values for q 1 and q 2 in the fractional-order model.
We carried out two series of optimizations. The first one was assumption-and bounds-free where we set the lower bounds of the four parameters (i.e., k f , k s , k f and k s ) to zero and the upper bound to infinity. The aim was to numerically prove that the fractional-order model can fit CBF flow response regardless the values (or the upper bound) of the four parameters. As this leads to unrealistic CBF responses, we defined permissible ranges for the four parameters and re-ran a constrained optimization. The selection of this range (or the upper bound) followed a visual inspection of the CBF response shapes using different values for the upper bound (e.g. 100, 50, 10 and 2). More importantly, it was also guided by previous literature on the upper limits of k f and k s (i.e. the range values in [8], [31]).

C. Parameters Sensitivity Analysis of NVC Fractional-Order Model
After showing that the fractional-order model gives raise to unique contributions to CBF response, we conduct a sensitivity analysis to study how the parameters of the fractional-order model affect or control the CBF response above and beyond the key parameters in integer-order model. To this end, we followed the sensitivity analysis approach conducted in [32]. Using a CBF response with a reference signal of fixed parameter values (shown in Table I), we slowly manipulate each parameter (and pairs of parameters) to quantify the deviation in output behavior by computing the L 1 -norm of the difference between the reference CBF and the manipulated CBF. Then we normalize the result by the norm of the reference signal to give more interpretable values.

D. Fitting Fractional-Order Model to Experimental CBF Data
In this section, we fit the proposed fractional-order model to real CBF data obtained from another study [21]. The inputoutput data used in this paper were acquired using electrophysiology and laser Doppler flowmetry (LDF) recordings, respectively. The data acquisition of the CSD data and CBF data are described briefly in the following and for more details we refer the reader to [21], [33], [34].
Hooded Lister rats with weight's range 200 − 300 g were used. The animals were prepared in a way to meet certain predefined specifications. After locating the whisker barrel cortex region, electrophysiology and LDF probes were placed based on the alignment of optical imaging maps with images of the cortical surface. The inserted electrophysiology electrodes were coupled to a data acquisition device (Medusa Bioamp, TDT, FL) with a custom written Matlab interface. Field potential recording was sampled at 6103.5 Hz with 16-bit resolution. To avoid the intrinsic spatial ambiguity which is inherent in the electrophysiology data, current source density (CSD) is used [33], [35]. The CSD profiles were given to us with a sampling time between each CSD point of 200 ms. To allow concurrent measurements of CBF, the LDF probe was positioned over the active area, adjacent to the electrodes. An LDF spectrometer including a low-pass filter was used to analyze the signal from the LDF probe with minimized errors due to the measurements noise. The CBF changes recorded with LPF were normalized to the baseline CBF which is collected for 8 s period before the onset of each trial. The CBF data acquired each 33.33 ms which correspond to an LDF with sampling rate of 30 Hz.
Regarding the experimental paradigm, it consists of conditioning block of stimulation followed by a probing block of stimulation per trial. For each trial, two blocks of stimuli were used. The first conditioning block has three different durations (2, 8, 16 s) which are followed by the probing block of 1 s duration. These two blocks are separated by 7 different time gaps (0.6, 1, 2, 3, 4, 6, 8 s). Therefore, there were 21 types of stimuli paradigms run for each animal. Last, the data were animal averaged. The averaged CSD and CBF data recorded for the 21 types of paradigms are shown in Fig. 4 (blue lines).
Using the CSD data as input to the model and the CBF data as its output, we evaluate how the model will fit the real-data over a certain range of parameters. To get those fittings, we first resampled the CSD input time points to the rate of CBF output data. Then, we used optimization function 'patternsearch' (in Matlab) to find the best set of parameters that can minimize the L 1 norm difference between the model outputs and the real CBF data, given the same input. The procedure was repeated for each subject, in each condition x stimulation combination. We then visualized the average fit for all subjects (while showing the standard deviation as a grey area), in each stimulation x gap combination.

III. RESULTS
The effects of how parameter change the CBF response have already been discussed elsewhere (see [5] for k s and k f , [20] for fractional orders q 1 and q 2 ]. Fig. 1 recalls the different shapes of CBF response in a range of values for each parameter. As shown in Fig. 1 A and 1 B, decreasing q 1 reduces signal width whereas decreasing q 2 increases the signal width and hence slows down signal decay. In most cases, decreasing q 2 also maintains a longer post-stimulus undershoot (associated with the slower decay) as well as delayed negative peak. On the other hand, q 1 only exhibits the negative undershoot in the first few fractions (approximately up to q 1 = 0.8) after which CBF response returns to baseline rapidly with no apparent negative undershoot. In general, while q 1 has more control over the positive aspects of the signal (i.e. overshoot width), q 2 has more control over the negative aspects of the signal like signal decay and undershoot duration. Positive signal amplitude is a common feature that both q 1 and q 2 effectively can change. Figs. 1 C and 1D shows the effect of k s and k f , holding both fractional orders q 1 and q 2 at 1 (hence, it represents the integer-order model). Decreasing k s increases signal oscillations while increasing k f eliminates the CBF undershoot.  Fig. 2 shows the results of the two series of optimizations described in Methods section. More specifically, signal dissimilarity increases as either q 1 and q 2 decreases. If the two models are equivalent, then a perfect match would be obtained resulting in an error of zero. However, increasing dissimilarity as a function of q 1 and q 2 indicates the noticeable effect of the new fractional parameters on the CBF response that cannot be obtained by tuning k f and k s of the integer-order model. More specifically, as we can see in Fig. 2 C, when varying q 1 while keeping q 2 at 1, we can see that CBF undershoots still matched, but overshoot amplitude is increasing due to the effect of the fractional parameter q 1 . We note that the actual values of error (or dissimilarity between the two models) are not of interest, but rather how they change as a function of fractional parameters.

B. Sensitivity Analysis for the Parameters of the Fractional-Order Model
Each subplot of Fig. 3 corresponds to the relative error of the CBF signal when one (or two) of the parameters takes values in a grid around the reference value within a specific range (shown in Table I). A unique extreme in the neighborhood of the reference value of the parameter is observed for both k s (Fig. 3 A) and k f (Fig. 3 B). Notably, the relative error is asymmetric around the observed global minimum. This asymmetry indicates the CBF response becomes less sensitive to changes in those parameters as the value of either parameter increases above the reference value (see Fig. 3 A and 3 B). Hence, initial guesses should always be taken less than the expected values of those two parameters [32].
Similarly, Fig. 3 C -3 F shows the relative error in CBF in the case of varying two parameters (again, while keeping other parameters at their reference values). Although a global minimum is shown around the reference values, those figures show a correlation between k s and both fractional orders which indicate that k s can, to some extent, "undo" the effects of both fractional orders. Concerning estimating the fractional parameters in light of those results, it can be argued that estimation of both k s and k f should always start with values less than their expected values as model dynamics are slower (less sensitivity) to big variations in values greater than the nominal values. For q 1 and q 2 , initial guesses are best set to 1 as CBF is less sensitive to variations in fractional orders when they approach zero. Colors represent the L1-norm of the difference between the two CBF outputs at the optimal values. In the constrained optimization, we used range from zero to 2. Signals of both models within the black rectangle are visualized in (C).

C. Fitting Fractional-Order Model to Experimental Data
Both integer-order and fractional-order models have been fitted to real CBF experimental data. Paired t-test of errors for each subject x condition combination show a significant difference between the errors derived from the two models (t(230) = −5.68, p-value < 1e−7; Mean f ractional = 60.706 < Mean integer = 62.913). Fig. 4 illustrate the results of fitting the fractional model to experimental CBF data. Each subplot contains a combination of stimulation x gap parameters that were used in the experiment. Model fits (in blue) show a very good fit for short stimulation paradigm (2 s) and a moderate fit for the 8 s stimulation paradigm. However, it captures less well the dynamics of the long (16 s) stimulation paradigm. Specifically, the model is only able to fit the second peak while missing other dynamics involved. Therefore, those results suggest that the fractional-order model has a moderately higher flexibility in fitting experimental data when compared with integer-order model due to the fact it generalizes the integer-order model.

IV. DISCUSSION
In this work, we have shown that the fractional equations have big potential in describing and characterizing a wide range of cerebral hemodynamic responses than the standard integer-order  models. The fractional-order model's ability to cover different response shapes is due to its flexibility and compliance offered by the two extra fractional order parameters, namely the fractional differentiation orders. Fractional parameters seem to characterize the CBF response in distinct ways by controlling the overshoots and the undershoots observed in the real CBF signal. Our sensitivity analysis clearly shows that these two parameters have different contributions in characterizing the cerebral hemodynamic determinants. While the variation of the positive overshoot of the CBF signal is sensitive to the value of the fractional differentiation order q 1 , the negative part of the CBF signal, namely the signal decay and undershoot, are more sensitive to the value of the second fractional differentiation order q 2 .
In this study, we noticed that the model fails to fit longer stimulation paradigms. When the framework involves highly nonlinear dynamics, this fact may limit the model's applicability to the long block design experiments of blocks longer than 8 seconds.
It is worth mentioning the strong correlation between k s and both q 1 and q 2 (see Fig. 3 C and 3 E). These correlations indicate that, within a specific interval, k s parameter can "undo" the effect of both fractional-order parameters. This collinearity translates into a difficulty of accurately estimating k s in the fractionalorder model.
Physiological interpretations of the fractional-order parameters q 1 and q 2 may be challenging. Similar to other NVC models, the model is meant to be descriptive to fit a wider range of experimental data that previous models cannot account for (see [4] for a discussion). Although these descriptive models have less physiologically interpretable parameters, they are extremely useful for comparison between different groups and conditions.

V. CONCLUSION
Modeling the NVC mechanism is undertaking considerable development. It is essential for a better assessment of BOLD fMRI data and also a necessary step towards understanding the physiology behind the complex phenomena involved and finding biomarkers that represent the key features observed in measured CBF profiles. Current NVC models still lack the flexibility to fit a wider range of observed experimental data. In this paper, the framework of fractional calculus is used to model the CBF response to neural activity. A fractional-order oscillator is proposed based on the well-accepted and known minimal model proposed by Friston et al.. Through an optimization scheme that compares the original integer-order model and the proposed fractional one, we showed that the added fractional parameters provide a unique contribution in describing the CBF that can not be captured using the integer-order model's parameters. Moreover, we assessed how sensitive CBF measure is to changes in the parameters of the model. Furthermore, using real neural activity-CBF data, the fractional model has proven capable of fitting wider CBF responses to both event and block design input paradigms. Although fractional model parameters are harder to interpret physiologically, they offer a great opportunity to compare groups or conditions. This paper does not deal with the estimation problem of the model's unknown parameters and fractional orders. This task may be the subject of a forthcoming study where model-based estimation techniques for fractional systems, for instance, the modulating function based non-asymptotic technique [36], [37], will be used to estimate the unknown parameters.