A Nitric Oxide–Modulated Variable-Order Fractional Maxwell Viscoelastic Model of Cerebral Vascular Walls

It is well known that the mechanical behavior of arterial walls plays an important role in the pathogenesis of vascular diseases. Most studies existing in the literature focus on the mechanical interactions between the blood flow and wall’s deformations. However, in the brain, the smaller vessels experience not only oscillatory forces due to the pulsatile blood flow but also structural and morphological changes controlled by the surrounding brain cells. In this study, the mechanical deformation of the cerebral arterial wall caused by the pulsatile blood flow and the dynamics of the neuronal nitric oxide (NO) is investigated. NO is a small diffusive gaseous molecule produced by the endothelial cells and neurons, which is involved in the regulation of cerebral blood flow and pressure. The cerebral vessel is assumed to be a hollow axial symmetric cylinder whose wall thickness is much smaller than the cylinder’s radius and longitudinal length is much less than the propagating wavelength. The wall is an isotropic, homogeneous linear viscoelastic material described by an NO-modulated variable-order fractional Maxwell model. A fractional telegraph equation is obtained for the axial component of the displacement. Patterns of wall’s deformation are investigated through numerical simulations. The results suggest that a significantly decreased inactivation of the neuronal NO may cause a reduction in the shear stress at the blood-vessel interface, which could lead to a decrease in the production of shear-induced endothelial NO and neurovascular disease.


INTRODUCTION
Cerebral vasculature plays a critical role in brain's metabolism and neurovascular conditions. The literature abounds with studies of cerebral blood flow and its interactions with the vasculature and brain cells [reviews of models and computer simulations can be found in chapter 4 of Drapaca and Sivaloganathan (2019)]. Recent advancements in technology have allowed researchers to gain invaluable knowledge about the intricate chemo-mechanical connections among neurons, glial, vascular, and blood cells. It is now acknowledged that the neurons and glial cells control the cerebral blood flow through the chemo-mechanical activation of the cells within the vascular wall (Attwell et al., 2010).
One of the many particles that facilitate the chemo-mechanical communications among the brain cells, blood, and the vascular wall is the nitric oxide (NO). A small diffusive gaseous molecule, the cerebral NO is mainly produced by shear-induced mechanotransduction at the blood-vessel interface (Sriram et al., 2016) and by synthesis reactions within the endothelial cells of the blood vessels and neurons (Forstermann and Sessa, 2012). NO diffuses and is removed from brain through some specialized chemical processes (Palacios-Callender et al., 2007;Unitt et al., 2010;Santos et al., 2011;Santos et al., 2012;Helms et al., 2016). In its role as a neuro-glial-vascular messenger, NO controls the cerebral blood flow and the release of neurotransmitters (Huang, 1999;Iadecola, 2004;Attwell et al., 2010;Contestabile et al., 2012;Iadecola, 2017). The regulation of the blood flow in brain is achieved through vasomotor mechanisms in which both the neuronal and endothelial NO are involved (Cockcroft, 2005;Metea and Newman, 2006;Attwell et al., 2010;Atochin and Huang, 2011;Petzold and Murthy, 2011;Contestabile et al., 2012;Lourenco et al., 2014;Haselden et al., 2020). However, throughout the entire cardiovascular system, the endothelial NO usually acts as a vasodilator (Attwell et al., 2010;Schuler et al., 2014). While the cerebral NO activity causes the local vasodilatation of downstream cerebral vessels at the neuroglial-vascular unit site, the NO diffusion within the vascular wall relaxes the smooth muscle cells leading to the so-called remote vasodilation due to the propagation of muscle's relaxation to the upstream arteries via the intercellular communications among endothelial and smooth muscle cells (Freed and Gutterman, 2017) facilitated by vascular gap junctions (Iadecola, 2004;Iadecola, 2017). Impaired cerebral NO production and/or decay can signal the presence of a neurovascular disease (Parker and Parks, 1995;Maurer et al., 2000;Wilkinson et al., 2004;Unitt et al., 2010;Santos et al., 2011;Haselden et al., 2020).
Although the NO-modulated vasodilation contributes to the mechanical deformation of the vascular walls of intracerebral vessels, existing mathematical models and corresponding numerical simulations of blood flow interacting with the deformable vascular wall do not account for it since their focus is the mechanics of big arteries in the presence of the pulsatile blood flow. The arteries are usually modeled as anisotropic, incompressible, nonlinear elastic materials whose constitutive stress-strain relationships may also incorporate collagen fibers' orientations and waviness, and/or the activation of smooth muscle cells (see Ebrahimi, 2009;Holzapfel and Ogden, 2010;Kim and Wagenseil, 2014;Espinosa et al., 2018 and references within). For instance, muscle activation has been modeled using 1) the (original or modified) Hill model (Hill, 1938), 2) an elastic strain-energy function dependent on the concentration of free intracellular calcium (Rachev and Hayashi, 1999), 3) a strain-energy function dependent on the chemical kinetics of the smooth muscle (Stalhand et al., 2008), and 4) a lumped Hodgkin-Huxley-like electrical circuit of the smooth muscle cell membrane coupled with a fluid compartment model describing the mass balances of considered ions and a contractile kinematics model regulated by intracellular calcium (Yang et al., 2003). Given the viscoelastic behavior of constituent cells (Kasza et al., 2007) of the blood vessels (and biological tissues, in general), the vascular wall has also been modeled as a viscoelastic material (Toth et al., 1998;Orosz et al., 1999;Holzapfel et al., 2002;Hodis and Zamir, 2008;Ebrahimi, 2009). Since these models do not incorporate the NO influence on the vascular wall, they are not able to predict neurovascular pathologies in which NO plays a critical role. Furthermore, the coupling of most of the above models with the cerebral NO dynamics will probably increase the already large number of model parameters that are practically impossible to find in a living brain with minimally invasive procedures using present-day technologies.
The aim of this study is to propose a novel mechanical model for cerebral arterioles that accounts for changes in the wall's mechanical behavior due to the NO-activated vascular cells and has few parameters. Orosz et al. (1999) used stress-relaxation measurements in cerebral arteries of human cadavers to show that four-and five-element generalized Maxwell viscoelastic models provide more accurate descriptions of the vascular wall mechanics than the two-element Maxwell viscoelastic model. One way to obtain a linear viscoelastic constitutive law with fewer parameters that perform as well as (or better) a springdashpot model with many elements is to use an integral formulation instead of the differential operator representation commonly associated with rheological (spring-dashpot) models. Integral constitutive laws admit equivalent differential formulations only for certain expressions of their relaxation functions (Gurtin and Sternberg, 1965;Drapaca et al., 2007). It was observed experimentally that relaxation functions represented as power functions of negative fractional exponents accurately describe the fading memory of many viscoelastic materials, including polymers and soft biological tissues (Nutting, 1921;Gemant, 1935;Gemant, 1936;Scott Blair and Coppen, 1939;Scott Blair and Coppen, 1942;Nutting, 1943;Guttinger, 1966;Caputo and Mainardi, 1971;Bagley and Torvik, 1983a;Bagley and Torvik, 1983b;Koeller, 1984;Torvik and Bagley, 1984;Suki et al., 1994;Mainardi, 2010). Integral constitutive laws with decaying fraction power-law relaxation functions have equivalent differential representations like those used in classic rheological models where integer-order derivatives were replaced by fractionalorder derivatives (convolutions between decaying fraction power functions and integer-order derivatives). These fractional viscoelastic models are causal at zero time (Bagley and Torvik, 1983b;Torvik and Bagley, 1984) and can be derived from molecular theories that incorporate the molecular complexity of various polymers (Bagley and Torvik, 1983b;Suki et al., 1994). In this context, the fractional order of the strain history models the "contribution of the long-chain molecules to the macroscopic stress" (Suki et al., 1994). Also, Suki et al. (1994) noticed the structural similarities between the lung tissue and some polymers and showed that a fractional viscoelastic model can successfully predict the viscoelastic behavior of lung tissue. Given that lungs are highly vascularized, it is reasonable to assume that a fractional viscoelastic model could be employed in studies of blood vessels mechanics.
Therefore, this study proposes a new NO-modulated variableorder fractional Maxwell viscoelastic model for the cerebral vascular wall and investigates its predictions through numerical simulations. The suggested constitutive equation is a generalization of the (constant order) fractional Maxwell viscoelastic model (see, for instance, Mainardi, 2010). The classic (first-order) Maxwell viscoelastic model has been previously used in the literature to model vascular walls (see, for instance, Hodis and Zamir, 2008), and therefore, this is the model which was chosen to be generalized here. According to Lorenzo and Hartley (2002), variable-order fractional operators are suitable in modeling, for instance, the fading memory characteristic of viscoelastic materials and the order memory which records the order in which memories are recalled [a recent review of applications of variable-order fractional operators can be found in Patnaik et al. (2020)]. In the model proposed here, the neuronal NO dynamics controls the order memory. With only three parameters, the proposed model can account for vessel's mechanics and the fact that the order memory of chemo-mechanical events is essential to the proper functionality of the constituent cells of the vascular wall. Lastly, the vascular wall is assumed to be isotropic and homogeneous.
As in Hodis and Zamir (2008), it is further assumed that the cerebral vessel is a hollow axial symmetric cylinder whose wall thickness is much smaller than the cylinder's radius, and longitudinal length is much less than the propagating wavelength. These assumptions together with the new constitutive equation led to a variable-order fractional telegraph equation for the unknown axial displacement of the vessel's wall. By assuming that the blood-vessel interface was exposed to the pulsatile blood flow and the vessel-tissue interface was tethered (Hodis and Zamir, 2008), an initialboundary value problem was obtained and solved using the separation of variables method and finite difference schemes.
The assumption of a tethered vessel-tissue interface is supported by experiments performed using high-resolution ultrasonic scanning that shows significantly larger longitudinal displacements of the wall's middle layers than the corresponding displacements of the wall's outer layers (Persson et al., 2003;Cinthio et al., 2006). Lastly, the variable fractional order was taken to be proportional to the concentration of neuronal NO. Two cases of neuronal NO synthesis (NOS) were considered: 1) a stepwise activation of NO with drastically decreased NO inactivation and 2) a dynamic activation of NO (Hall and Garthwaite, 2006). Case 1 may not be physiological because the main neuronal NOS activation matches the short pulsatile Ca 2+ -calmodulin binding activation in dendritic spines (Hall and Garthwaite, 2006), and thus, this case may model impaired neuronal NO dynamics. Indeed, the temporal variation of NO resembles a stepwise function when the NO inactivation is drastically reduced by either brain ischemia (Santos et al., 2011) or the inhibition of cytochrome c oxidase (COX) (Palacios-Callender et al., 2007;Unitt et al., 2010). Brains of Alzheimer disease (AD) patients may suffer from a COX inhibition and thus an accumulation of NO (primarily in the temporal cortex and hippocampus) because a selective defect of COX causing a significant reduction in COX activity was found in AD brains (Parker and Parks, 1995;Maurer et al., 2000). Numerical simulations generated in MATLAB show smaller displacements within the vascular wall in the case of stepwise NOS than in the case of dynamic NOS. Also, the shear stress at the bloodvessel interface is smaller in the case of stepwise neuronal NOS. Since the production of shear-induced endothelial NO is proportional to the shear stress at the inner boundary of the vessel's wall (Sriram et al., 2016), it follows that less endothelial NO will be produced in this case which could ultimately lead to neurovascular disease. Thus, the model could be used as a complementary clinical tool for early detection of disease and intervention (Alzheimer, for instance).
The study is structured as follows. The proposed mathematical model and the initial-boundary value problem under investigation are presented in Mathematical Model. The corresponding semi-analytic solution is given in Semi-Analytic Solution. Numerical simulations are shown in Results, which is followed by the last section containing a discussion of the results and final conclusions.

MATHEMATICAL MODEL
The vessel is modeled as a hollow horizontal axial symmetric circular cylinder of radius a and thickness h ( Figure 1). As in Hodis and Zamir (2008), it is assumed that the vessel's wall is made of a homogeneous, linear viscoelastic material that is tethered at the outer boundary (the wall-tissue interface) and exposed to oscillations caused by the heart pulsations at its inner boundary (the blood-wall interface). In addition, it is assumed that the spatial variations of the oscillations can be neglected in a first approximation because they are of the same order of magnitude as the propagating wavelength (which is approximately 10 m at a frequency of 1 Hz in the systemic circulation) while the vessel's length is much smaller than the wavelength. Thus, only the time variations of the oscillations at the inner boundary of the vessel's wall will contribute to the wall's deformation. If (r, θ, x) are the cylindrical coordinates, then the modeling assumptions made so far reduce the equation of motion which is relevant to the work presented here to the following equation (Hodis and Zamir, 2008): FIGURE 1 | A schematic of the geometric domain occupied by a vessel of radius a and wall thickness h. The model assumes a hollow horizontal axial symmetric circular cylinder whose wall has a fixed outer boundary and an oscillatory inner boundary.
for (r, t) ∈ (a, a + h)×[0, T]. In Eq. 1, ρ is the mass density, ξ(r, t) is the axial displacement of a material point of the vessel's wall, and σ rx (r, t) is the shear stress of a wall's point. Lastly, because the wall thickness h is usually much smaller than the radius a, the curvature term σ rx /r is smaller than the other terms of Eq. 1, and thus, it can be neglected. Therefore, Eq. 1 reduces to (Hodis and Zamir, 2008): In this study, it is assumed that the mechanical behavior of the vessel's wall is described by the constitutive equation of a variable-order fractional Maxwell linear viscoelastic material: where, based on the model's assumptions, the infinitesimal strain ε rx reduces to: The constitutive Eq. 3 has three physical parameters: E, the modulus of elasticity, μ α , and α(t), 0 < α(t) ≤ 1. Without loss of generality, it is assumed further that ε rx (r, 0 +) 0 and σ rx (r, 0 + ) 0 (Bazhlekova and Bazhlekov, 2017). Lastly, the variable-order fractional derivative used in Eq. 3 is as follows (see for instance Ramirez and Coimbra, 2010;Moghaddam and Machado, 2017): where m ∈ 1, 2, 3 . . . , is an arbitrary, continuously m th -order differentiable function whose m th -order derivative is is a continuously differentiable function, and Γ(z) From Eq. 5 and the zero initial conditions satisfied by ε rx and σ rx , it follows that Eq. 3 becomes the fractional Maxwell model for α(t) a constant in the interval (0, 1) and reduces to the classic Maxwell model for α(t) 1, ∀t > 0. Note that in the classic Maxwell model parameter μ α becomes the viscosity denoted by μ. When the fractional order is a constant α ∈ (0, 1), properties of the (Caputo) fractional derivative Eq. 5 and the Laplace transform can be combined to obtain the following expressions for the creep compliance J and relaxation modulus G (Mainardi, 2010;Bazhlekova and Bazhlekov, 2017): where g E/μ α and the Mittag-Leffler function is, by Figure 2 shows plots of dimensionless material functions EJ and G/E given by Eq. 7 versus the dimensionless time t t/T for various values of α. Thus, the material functions for a variable order α(t) ∈ (0, 1) will combine the behavior of multiple material functions shown in Figure 2.
Replacing Eq. 4 in Eq. 3, applying operator z −α(t) t to Eq. 3, and using Eq. 6 give:  Differentiating Eq. 8 with respect to the spatial variable r and using Eq. 2 yield: Lastly, by replacing Eq. 9 in Eq. 2 the following equation is obtained: Eq. 10 is a variable-order fractional telegraph equation (or twoterm time-variable fractional diffusion-wave equation).
It is assumed that the variable fractional order α(t) ∈ (0, 1) models the temporal effects of the neuronal NO on the vessel's wall. At this incipient stage, the shape of α(t) is considered to look like the temporal profile of the concentration of neuronal NO. Hall and Garthwaite (2006) proposed two models of neuronal NO dynamics: 1) stepwise activation of NO and 2) dynamic activation of NO. Both models include the NO inactivation. The equation of the first model is as follows: whose approximate analytic solution for zero initial condition is (Mehala and Rajendran, 2014) as follows: The equation of the second model is as follows: whose approximate analytic solution for zero initial condition is (Mehala and Rajendran, 2014) as follows: Here, [NO] is the concentration of NO, ν 1 is the constant rate of NOS, V max is the maximum rate at saturating concentration of NO, K m is the concentration of NO at which the reaction rate is V max 2 , k 1 and k 2 are constant kinetic parameters, and κ V max /K m . The values in Table 1 were chosen such that the numerical solutions to the above differential equations for [NO] and the corresponding approximate analytic solutions agree (Mehala and Rajendran, 2014). Figures 3A,B show plots of the NO concentrations given by Eqs 13, 14, respectively, while Figures 3C,D show two proposed profiles for the variable order α(t) ∈ (0, 1) which are obtained by simply multiplying the NO concentrations shown in Figures 3A,B by a factor of 10 2 . The scaling factor was chosen as follows. As seen in Figures 3A,B, the maximum concentrations of NO are in the low picomolar rage, while physiological values are in the low nanomolar range (Hall and Garthwaite, 2006). A mere multiplication of the NO concentrations by 10 2 brings the maximum values within the physiological range for the NO concentrations. Since these scaled concentrations are still less than 1 nM , they can be used as expressions for α(t). Another advantage of this scaling factor is that significant differences are observed between the displacements and the shear stresses at the blood-vessel interface corresponding to the two profiles of α(t) shown in Figures 3C,D (see later).
The aim of the study is to find a semi-analytic solution to the initial-boundary value problem (10-12) for two variable orders given by scaled Eqs 13, 14. The first step is to formulate a corresponding non-dimensional problem. By introducing the dimensionless quantities: Due to a lack of experimental data, the value of T was fixed first, and then, the values of f and g were chosen such that a semi-analytic solution could be found for the initial boundary value problem Eqs 15-17. The normalized value of the wall's density was chosen for mathematical simplicity.

Drapaca
NO-Modulated Viscoelastic Model of Vessels the initial-boundary value problem (10-12) transforms into the following dimensionless form: In Eq. 15, f E/(ρh 2 ). The solution for problem (15-17) is presented in the next section.

SEMI-ANALYTIC SOLUTION
A classic approach of solving initial-boundary value problems for partial differential equations is used. Look for a solution to problem (15-17) of the form: such that V(0, t) ξ(0, t) sin(ωT t), V(1, t) ξ(1, t) 0 and z 2 V z r 2 0. It is straightforward to solve for V( r, t) and find that: Combining Eqs 18, 19 yields: ξ r, t (1 − r) sin ωT t + W r, t .
By replacing expression (29) in problem (27, 28), the following problem for the unknown function Z n ( t), n ∈ 1, 2, 3 . . .} { is obtained: where B gT α( t) , C n (nπT) 2 f , D n 2nπT 2 f , and F n 2ωT (−1) n+1 +1 nπ . Once the solution to the initial value problem (30, 31) is found, Eqs 20, 26, and 29 can be combined to find the solution to the original problem (15-17): Thus, the last step is finding a solution to problem (30, 31). If the fractional order is a constant α ∈ (0, 1), then a closedform solution of problem (30, 31) can be obtained by using either the Laplace transform method (Gorenflo et al., 2014) or an operational method proposed by Luchko and Garenflo (1999) (it is worth noticing here that the Laplace transform method belongs to the class of operational methods). Not only that this analytic solution is cumbersome and difficult to implement in a computer program but also may not be generalizable to the case of variable fractional order. Therefore, looking for a numerical solution to problem (Eqs 30, 31) is a better approach. Eq. 30 is a linear multiterm fractional differential equation which for a constant fractional order α ∈ (0, 1), can be solved numerically using the elegant explicit, implicit, and predictor-corrector methods involving product-integration rules proposed by Garrappa (2018). For the variable-order fractional differential Eq. 30, the finite difference schemes given by Moghaddam and Machado (2017) can be used. Thus, in this study, the second-order derivative of Eq. 20 will be approximated by a center difference scheme, and the variable-order fractional derivative will be approximated by the forward difference scheme proposed in theorem 3.1 of Moghaddam and Machado (2017). For the sake of completeness, the numerical discretization of Eq. 30 is presented further.

RESULTS
Numerical scheme Eq. 35 was implemented in MATLAB, and plots of the distributions of dimensionless displacements within the vessel's wall (Eq. 32) and the dimensionless shear stress at the blood-vessel interface (Eq. 37) were generated for two expressions of the variable fractional order: the stepwise NOS and the dynamic NOS. The values of the parameters used in the numerical simulations are given in Table 1. Numerical simulations used a step size Δ t 0.00025 and 20 terms of the Fourier series in Eq. 32. Figure 4 shows temporal profiles of the dimensionless axial displacement ξ( r, t) at various radial locations inside the vascular wall and for constant fractional orders. These plots were generated using MATLAB's function MT_FDE_PI12_PC developed by Garrappa which is an implementation of the predictor-corrector scheme proposed by Garrappa (2018). The same results are obtained with the forward scheme Eq. 35 when α is constant. This agreement of numerical solutions validates scheme Eq. 35. The plots of Figure 4 show that while the amplitudes of the oscillations are approaching zero as r → 1 for all values of α, the decaying in  amplitude decreases faster as α increases. The solutions for constant fractional orders provide a first glimpse into the shapes of the displacements for variable fractional orders since these are expected to be a combination of the behaviors of the solutions for constant α. Figures 5, 6 show spatio-temporal variations of the dimensionless axial displacement ξ( r, t) for the variable fractional order α(t) corresponding to the case of decaying dynamic NOS ( Figure 5) and to the case of stepwise NOS ( Figure 6). The amplitude of the oscillations appears to be bigger in the dynamic NOS case than in the stepwise NOS case. This is confirmed in Figure 7 which shows temporal profiles of the dimensionless axial displacement at various radial locations inside the vascular wall which was extracted from the surfaces presented in Figures 5, 6. Since the results in Figure 7 resemble those in Figure 4, it can be concluded that the solutions for variable fractional orders are plausible. According to Figures 3C,D, the maximum value of α(t) for the case of dynamic NOS is 0.15 which is smaller than the almost constant value of 0.5 that α(t) in the case of stepwise NOS. The results in Figure 7 confirmed what was already known from Figure 4, namely that the amplitude of oscillations corresponding to the higher values of α(t) is smaller than the one corresponding to lower values of α(t).
To better understand the differences noticed between the oscillatory patterns corresponding to the two cases of stepwise NOS and dynamic NOS, the dimensionless shear stress at the blood-vessel interface was calculated using Eq. 37 and MATLAB's built-in function gradient. The results of this comparison are shown in Figure 8. The amplitude of the shear stress is bigger in the case of dynamic NOS which was expected given that the amplitude of the axial displacement is also bigger in this case. By varying the scaling factor of order α(t) corresponding to the case of stepwise neuronal NOS, the  difference between the shear stresses at the inner boundary of the vessel wall corresponding to the two cases 1) vanishes when the fractional order of the stepwise neuronal NOS gets smaller and 2) is almost unchanged when the fractional order of the stepwise neuronal NOS approaches 1. The possible implications of this finding will be discussed in the next section.

DISCUSSION
The main contribution of this study is modeling the wall of a cerebral blood vessel using an NO-modulated variable-order fractional Maxwell viscoelastic model. The variable fractional order is assumed to be proportional to the neuronal NO dynamics, and thus, the order memory introduced by this choice guides the pattern of the fading memory of this viscoelastic material. Two cases of NOS are considered: a stepwise activation of NO and a dynamic activation of NO. Following the approach of Hodis and Zamir (2008), a variable-order fractional telegraph equation for the axial displacement of the wall was obtained which was solved under the assumptions that the outer boundary of the vascular wall was tethered, and the inner boundary of the wall was exposed to the pulsatile blood flow. Numerical simulations were created in MATLAB using numerical scheme (Eq. 35) and Garrappa's function MT_FDE_PI12_PC (Garrappa, 2018) for constant fractional orders. The function MT_FDE_PI12_PC was used to validate scheme (Eq. 35) in the case of constant fractional orders. The main finding of these simulations is that a significantly decreased inactivation of the neuronal NO causes a reduction in the shear stress at the blood-vessel interface which could lead to a decrease in the production of shear-induced endothelial NO and ultimately to neurovascular disease.
Testing the computer code for various values of the parameters f , g, and T which are present in the coefficients of Eq. 15 showed that both approaches, scheme (Eq. 35) and MT_FDE_PI12_PC, were sensitive to these parameters for variable and constant fractional orders. These numerical schemes are stable only for narrow ranges of f , g, and T. These parameters couple mechanical parameters E, μ α , wall thickness h, and the characteristic time T. Thus, if these narrow ranges of f , g, and T are subsets of their corresponding physiological ranges, this stability issue could be disregarded. However, for E 6 × 10 5 N/m 2 , μ α 1.5 × 10 5 Kg/(m × s 2−α ) derived from Hodis and Zamir (2008) for arteries and (α 1), a wall thickness h 10 −3 m measured in human intracranial arteries in vivo (Yuan et al., 2021), and a more accurate value of the wall density ρ ≈ 1.1 × 10 3 Kg m 3 (IT'IS Foundation, 2021), the following values are obtained: f 5.45 × 10 8 s −2 and g 4 s −α . These values, together with the chosen T 10 s, make Eq. 15 stiff since the order of magnitude of the coefficient of the righthand side term is much bigger than the magnitudes of the coefficients of the left-hand-side terms. For a smaller wall thickness which is more appropriate for an intracerebral arteriole which experiences NO-modulated vasodilation, the order of magnitude of the coefficient of the right-hand-side term of Eq. 15 becomes even bigger. The proposed numerical schemes are unstable for a stiff Eq. 15, and a mere decrease of the step size Δ t does not resolve this issue. Finding algorithms for stiff variable-order fractional differential equations is an open problem in numerical analysis. Nevertheless, for a characteristic time T of order of ms which may be more suitable for NO dynamics in vivo (Hall and Garthwaite, 2006), a possibly lower value of E corresponding to cerebral arterioles (Medical Physiology, 2021), and a more realistic value of the wall density, the numerical stiffness can be avoided, and thus, the presented results will hold.
Emerging imaging techniques could help validate the proposed model in animal models and thus make the model relevant to clinical applications. For instance, a multimodal in vivo magnetic resonance (MR)/electron paramagnetic resonance (EPR) spectroscopy/fluorometry could be used to visualize NO production and spatio-temporal distribution (Sharma et al., 2014). Also, intravascular optical coherence tomography could be used for the in vivo real-time estimation of the vascular stiffness (Potlov et al., 2020). By combining these imaging techniques and the high-resolution ultrasonic scanning of Persson et al. (2003) and Cinthio et al. (2006) for the visualization of wall's longitudinal displacements, the cerebral NO dynamics and vascular wall mechanics may be investigated simultaneously. This approach can be used to estimate α(t), f , g, and T and validate the mathematical model proposed here. If a significant decrease in the neuronal NO inactivation is observed, then the model's prediction could suggest the use of a preventive therapy [such as NO inhalation proposed by Terpolilli et al. (2016)] to reduce imminent brain damage.
Until the above-mentioned imaging techniques are adapted and approved for clinical applications, animal models can be used to find healthy physiological ranges for the model's parameters and investigate the potential of using these parameters as biomarkers. For example, animal models of cerebral ischemia have shown that the decreasing amount of endothelial NO will act as a protective agent for a few minutes after the injury, while the amount of neuronal NO will increase causing neuronal injury (Huang, 1999;Wei et al., 1999). If the production of shear-induced endothelial NO mediated by the blood flow ceases to happen following ischemia, then the stepwise concentration of neuronal NO is a consequence of the significantly decreased NO inactivation caused by ischemia (Santos et al., 2011) and thus, according to the prediction made by the proposed model, may contribute to the reduction of the localized endothelial NO produced by the lower shear stress in the vessel wall at its inner boundary. Thus, the model's prediction may explain the interplay between the endothelial and neuronal NO seen in cerebral ischemia. Lastly, model's parameters estimated using the multimodal in vivo imaging techniques mentioned above could suggest the presence of cerebral ischemia and thus be used as a complementary diagnostic tool. Improving the mathematical model may provide more sensitive biomarkers and helps avoiding the numerical stiffness issue mentioned earlier. According to Iadecola (2004), the NOmodulated local vasodilation of the intracerebral arterioles and capillaries propagates upstream in the vascular network which causes an increase of blood flow in the upstream arteries that leads to increased shear stress at the blood-vessel interface and thus an increase in the amount of shear-induced endothelial NO and a further flow-mediated vasodilation. The proposed model represents the blood flow as an inner boundary condition, and thus, this mechanism of global production of shear-induced endothelial NO is not accounted for. Thus, the decrease in the local production of the shear-induced endothelial NO due to a stepwise neuronal NOS predicted by the model proposed in this study might not be significant enough to cause adverse effects if the cerebral blood flow is intact. Coupling the deformation of the vessel's wall and the blood flow and incorporating in the model the mechanism described above should provide a better prediction of the amount of shear-induced endothelial NO and the possible role that the neuronal NO dynamics may play in this process. Introducing more detailed information about the complex, multi-layered structure of the vascular wall could also enhance the accuracy of model's predictions. Lastly, the full three-dimensional fluid-structure interaction problem should be formulated and solved since the simplifying assumptions made by Hodis and Zamir (2008) and implemented here that reduced the model to a onedimensional problem might not be valid for cerebral arterioles given the overly complex geometry of the cerebral vascular network (Reina-De La Torre et al., 1998). The parameters of this enhanced mathematical model could be found either from in vitro/in situ measurements or from in vivo, real-time observations using medical imaging techniques.
In conclusion, the study proposes a novel NO-modulated variable-order fractional Maxwell viscoelastic model of the cerebral arterioles and investigates the effects of the neuronal NOS on the mechanical behavior of the vessel's wall. Numerical simulations show how neuronal NO dynamics influence the deformation and shear stress within the vascular wall. A generalization of this model to a three-dimensional geometry and the incorporation of the blood flow into the model should provide a better understanding of the coupling between NO dynamics and mechanical damage and their combined role in neurovascular diseases.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material; further inquiries can be directed to the corresponding author.