Poiseuille Flow of a Non-Local Non-Newtonian Fluid with Wall Slip: A First Step in Modeling Cerebral Microaneurysms

Cerebral aneurysms and microaneurysms are abnormal vascular dilatations with high risk of rupture. An aneurysmal rupture could cause permanent disability and even death. Finding and treating aneurysms before their rupture is very difficult since symptoms can be easily attributed mistakenly to other common brain diseases. Mathematical models could highlight possible mechanisms of aneurismal development and suggest specialized biomarkers for aneurysms. Existing mathematical models of intracranial aneurysms focus on mechanical interactions between blood flow and arteries. However, these models cannot apply to microaneurysms since the anatomy and physiology at the length scale of cerebral microcirculation are different. In this paper we propose a mechanism for the formation of microaneurysms that involves the chemo-mechanical coupling of blood and endothelial and neuroglial cells. We model the blood as a non-local non-Newtonian incompressible fluid and solve analytically the Poiseuille flow of such a fluid through an axi-symmetric circular rigid and impermeable pipe in the presence of wall slip. The spatial derivatives of the proposed generalization of the rate of deformation tensor are expressed using Caputo fractional derivatives. The wall slip is represented by the classic Navier law and a generalization of this law involving fractional derivatives. Numerical simulations suggest that hypertension could contribute to microaneurysmal formation.


Introduction
Cerebral aneurysms are abnormal swellings of the vasculature with high rupture risk. Although most cerebral aneurysms develop near bifurcations of the large arteries [Liebeskind(2016)], microaneurysms can also form on arterioles located either in the retinas of patients with diabetes [Aguilar(2003), ] or deeper inside the brain tissue [Rosenblum(1977), Rosenblum(2008)]. Aneurysmal ruptures could cause life threatening intracranial hemorrhages or, in the case of diabetic retinopathy, could lead to blindness. With the exception of retinal microaneurysms, finding cerebral aneurysms before they rupture is very hard because either there are no symptoms or the symptoms can be linked to other more common brain diseases. Effective biomarkers for cerebral aneurysms do not exist yet [Hussain(2015)]. Mathematical models incorporating experimental and clinical observations of aneurysms can help us better understand aneurysmal pathophysiology, assess therapeutical success, identify biomarkers for the risk evaluation of aneurysmal formation and rupture, and develop pharmacological remedies for the inhibition of aneurysmal growth [Liang(2012)].
Existing models of intracranial aneurysms are macroscopic and study the roles played by blood hemodynamics, arterial wall biomechanics and geometry, blood-wall mechanical interactions, mechanobiology, and the intracranial environment in the onset, evolution and rupture of aneurysms. With few exceptions, the blood is considered an incompressible Newtonian fluid in laminar flow through an artery of idealized geometry and having impermeable thin walls modeled as either rigid or viscoelastic solids. Recently, computer simulations of cerebral hemodynamics involving patientspecific geometries reconstructed from medical images and computational fluid dynamics (CFD) software have been performed and some geometrical and hemodynamical parameters were proposed as biomarkers for aneurysmal rupture prediction [Kulcsár(2011), Jou(2008), Meng(2014), Cebral(2005), Oubel (2007), Xu(2011), Ngoepe(2016), Grinberg(2013)]. A CFD model using an incompressible non-Newtonian constitutive law for the blood flowing with no-slip through a rigid-walled vessel was given in [Cebral(2005)], while CFD models of incompressible Newtonian blood circulating with no-slip through vessels with deformable walls were presented in [Oubel (2007), Xu(2011)]. The role played by blood clotting in aneurysmal growth and rupture was studied in [Grinberg(2013)] using a complex atomistic-continuum model of hemodynamics, and [Ngoepe(2016)] using mixture theory and diffusion-reaction equations for some blood proteins responsible for clotting formation and growth. A multi-particle collison dynamics approach was used in [Paudel(2011)] to simulate blood flow through an axi-symmetric circular cylinder with slip, and the mechanical response of a deformable arterial wall to mechanical forces exerted on it by the exterior surrounding structures was studied in [Hodis(2009), Hodis(2011)]. However, the mechanical parameters of these models lack predictive robustness and thus they cannot be used as biomarkers for aneurysms. In addition, the models require knowledge of numerous parameters which are not only difficult if not impossible to measure in vivo but also their relevance and established physical meanings might be obscured by the complex structure and dynamics possessed by this living system. This severely restricts the applicability of these models in clinical practice at this time. Lastly, the above mentioned mathematical models cannot be used in studying microaneurysms since different brain structures and chemo-mechanical interactions are relevant at this length scale. However, mathematical models for microaneurysms do not appear to exist yet although with the increasing prevalence of diabetes throughout the world [Bhavsar(2010)] such models could be critical not only for understading the underlying mechanisms of retinal microaneurysms and diabetic retinopathy but also in developing effective treatments.
In this paper we propose the following mechanism for the formation of microaneurysms. A chemo-mechanical entanglement of blood and endothelial cells at the blood -vascular wall interface that might be caused by a chemical imbalance and/or the forward and backward moving waves reflected from a bifurcation site [Zamir(2016)] will cause slip of the blood at the wall. The slip will produce deviations from the expected shear stress level at the wall which will make the red blood cells release the energy-making enzyme called ATP [Fisher(2003)] that will activate a chemical response from the endothelial cells of the wall which will lastly trigger neuronal-induced release of vasodilators from astrocytes whose endfeets continuously probe and control the amount of mechanical forces on the arterial wall and facilitate complex mechanotrasduction processes [Iadecola(2007), Metea(2006), MacVicar(2015, Gordon(2007), Petzold(2011)]. According to [MacVicar(2015), Gordon(2007)], the two astrocyte endfeet use Ca 2+ -mediated chemical signals to control vasoconstriction (one endfoot) and vasodilation (the other endfoot). We hypothesize that slip could cause an excessive localized vasodilation which could lead to the formation of a microaneurysm. If some vasoconstriction is still present during this abnormal vasodilation, then the microaneurysm might evolve without rupture. However, if the vasoconstriction is missing, we conjecture that there is a high risk of microaneurysmal rupture. If this mechanism of microaneurysmal formation, growth and rupture is possible in practice then biomarkers for microaneurysms could most likely be found among chemical species located outside an arteriole at risk of microaneurysmal development. In figure 1 we show a schematic of the structures that are involved in the formation of microaneurysms.

Risk of Microaneurysm Formation
Astrocyte Endfeet Figure 1: Schematic of the chemo-mechanical coupling among blood, cerebral arterioles and astrocytes. Slip at the blood-arteriole's wall interface could trigger an excessive vasodilation of the wall that is controlled by the astrocyte endfeet. This process could further lead to the formation of a microaneurysm. Brain cells and vasculature are immersed in cerebrospinal fluid (CSF), a water-like fluid produced by the brain which facilitates the intracellular and extracellular transport of molecules and ions within the brain.
A mathematical model capable of predicting the formation, growth and risk of rupture of a microaneurysm needs to incorporate the wellknown non-Newtonian behavior of blood when flowing through smaller vessels [Fung(1993), Fisher(2003, Tomaiuolo(2016)], the deformability of the vascular wall [Monson(2003)], and relevant mechanotransduction processes taking place among red blood cells, endothelial cells and astrocytes [MacVicar(2015), Gordon(2007)]. The first step in building such a model is presented in this paper. Here we study the fully developed steady flow of blood through an arteriole due to an externally imposed pressure gradient and in the absence of body forces under the following assumptions: 1). blood is an incompressible non-Newtonian fluid with a non-local stress-deformation rate relationship, 2). the arteriole is a horizontal axisymmetric circular pipe with rigid and impermeable walls, and 3) slip at the blood-vascular wall interface. In the constitutive equation of blood a generalized rate of deformation tensor is defined using non-local fractional order Caputo spatial derivatives of order α, m−1 < α < m, m = 1, 2, 3.... This is an extension of our model in [Drapaca(2012)]. The physical parameter α is seen as a measure of the long range interactions among fluid's particles during flow and it could vary with time, location, temperature, pressure, shear rate, and/or concentration of particles. In this paper we assume that α is constant. As is customary in fluid mechanics, the Cauchy stress tensor is a linear function of this rate of deformation tensor where the corresponding coefficient of proportionality µ that relates these two tensors is a physical parameter that depends on a characteristic length scale of the problem and the fractional order α. In this paper we consider µ to be constant. We assume further that the flow has no slip at the wall until the shear stress at the wall reaches a critical value after which slip happens. The slip at the blood-vascular wall interface is modeled using the classic Navier slip condition and a new non-local generalization of Navier slip condition that involves a fractional order Caputo derivative. The parameter α in this generalized slip condition models the non-local entanglement and aggregation of the blood and endothelial cells. Thus we use fractional order derivatives in space to link cells from different structures. For now the fractional derivatives model the macroscopic behavior of aggregated blood cells and the macroscopic slip of entangled blood and endothelial cells. In our future work we plan to use the same approach involving spatial fractional order derivatives to account for the neuroglial effects on the arterioles, as well.
Since at this first modeling step we do not consider the role played by the neuroglial dynamics in the vasodilation of arterioles, we show numerical simulations for arterioles and arteries for two characteristic slip lengths and various values of α ∈ (0, 2). The results show that the velocity profiles and volumetric flow rates differ for the two slip lengths and for two values of α chosen such that one is less than 1 and the other is greater than 1. The velocity profiles and flow rates corresponding to the two slip conditions are relatively similar for the same slip length and value of α, except for the case when the radius of the pipe is small (arteriole) and the slip length is big. In this case the velocity at the wall and the volumetric flow rate increase dramatically with a small increase in the pressure gradient. This suggests that if the blood pressure is consistently above a critical shear stress level responsible for slip which may happen for instance in the case of hypertension (caused by diabetes or not) then the produced blood flow rates will be persistently high. This repeatedly high flow rate will not only cause (mainly mechanically-driven) vasodilation but also activate, through mechanotransduction, neuroglial processes involved in the release of vasodilators that will further increase vasodilation which ultimately could lead to the formation of microaneurysms. Our findings are in agreement with some clinical observations [Rosenblum(1977), Rosenblum(2008)]. We emphasize that, to the best of our knowledge, this is the first time that such a mechanism for microaneurysms has been proposed in the literature. Its experimental and/or clinical validation will be the aim of future work.
In conclusion, the contributions of this paper are as follows: 1). a mechanism for the formation, growth and rupture of microaneurysms that involve the coupled dynamics of blood, arterioles and neuroglia was proposed, 2). a non-local constitutive law involving spatial fractional order Caputo derivatives and only two physical parameters (α and µ) was proposed to describe the non-Newtonian behavior of blood, 3). a non-local slip condition at the blood-vascular wall interface was introduced, 4). analytic solutions to a Poiseuille flow through an axi-symmetric circular rigid and impermeable pipe with wall slip were proposed, and 5). numerical simulations indicated that hypertension might contribute to the formation of microaneurysms which agrees with some clinical studies.
From a physical point of view, the fractional order α is a measure of the non-local interactions of fluid particles during flow and could also control their deformability and aggregation patterns. The other physical parameter of the model, µ, reduces to the apparent viscosity of the fluid when α = 1.
The structure of the paper is as follows. Our mathematical model and corresponding analytic solutions are presented in section 2. In section 3 we show and discuss some numerical simulations for blood flowing through an artery and an arteriole for various values of α and two characteristic slip lengths. The paper ends with a section of conclusions and future work.
The reasons for using fractional order integro-differential operators in equations (2) and (3) are threefold. In general, it is very difficult to use the physics of a problem to create non-local operators especially for materials with complex microstructural components that act and interact on multiple length and time scales. Fractional order operators model well spatio-temporal uncertainties in dynamical systems [West(2015)] which is desirable in modeling for instance chemo-mechanical behavior of living biological tissues. Lastly, using these operators is mathematically convenient because we can apply known concepts from fractional calculus. Under certain mathematical assumptions the operators used in (2) and (3) could perhaps be fit into the general framework for non-local calculus proposed in [Du(2013)].
According to [Drapaca(2012)], the region of integration H = [a1, b1] × [a2, b2] × [a3, b3] ⊂ R 3 respresents the region of influence of X and contains all the material points involved in the deformation of X into x. The limits ai, bi, i = 1, 2, 3 can be functions of X, t and α α α. H can be determined by the physics of interactions between particles, by the geometry of the domain occupied by the body under observation, or by curve fitting to experimental data. The multiple length scales introduced by H vary during deformation and are related through power laws of orders contained in the parameter matrix α α α(t). The components of matrix α α α(t) can be regarded as measures of the dynamic deformation of a body with evolving microstructure.
Definition 2: The deformation gradient of order α α α(t), with t > 0 and αIi ∈ R, I, i = 1, 2, 3 is either: or: Definition 2 is an adaptation of definition 8 in [Drapaca(2012)]. We notice that the deformation gradient of order α α α(t) given by formulas (4) and (5) resembles the definitions of the left-sided and right-sided Riemann-Liouville fractional derivatives of order αIi(t), I, i = 1, 2, 3. Using some properties of fractional order derivatives it can be shown that inversions can be defined only for the deformations (1) and (3). The very long spatial memory of each material point during a deformation (2) is lost and thus an inversion of this motion is not possible.
Other strain and strain rate tensors of order α α α(t) can be further obtained by combining definitions known from continuum mechanics and formulas (4) or (5). For instance, the velocity gradient tensor of order α α α(t) is given by: where v is the velocity field associated to the deformation given in definition 1. It follows that the rate of deformation tensor of order α α α(t) is D α α α(t) = L α(t) α(t) α(t) + L T α(t) α(t) α(t) /2. In particular, if H is independent of t and α α α(t) and α α α is a constant matrix, then the velocity fields v = ∂x ∂t corresponding to the three cases in definition 1 are given by formulas (1)- (3) where χ χ χ is replaced by ∂χ χ χ ∂t . In this paper we will work under these simplifying assumptions and use the rate of deformation tensor of order α α α to model the non-local non-Newtonian behavior of blood.

Poiseuille Flow
We present now the mathematical model of the three-dimensional fully developed steady laminar flow of an incompressible non-local non-Newtonian fluid through an horizontal circular pipe with rigid and impermeable walls. The flow is axi-symmetric and is driven by an externally imposed pressure gradient. For simplicity, body forces are neglected. In cylindrical coordinates (r, θ, z), the components of fluid's velocity are: and therefore the equation of continuity: is identically satisfied. The equilibrium equations involving the components of σ σ σ(r, z), the Cauchy stress tensor of the fluid, in cylindrical coor-dinates are: The only non-zero components of σ σ σ are given by the following constitutive formulas: where p is the hydrostatic pressure of the fluid, and µ and α ∈ (m − 1, m), m = 1, 2, 3... are the two physical parameters of the model. In the constitutive formula (10), the shear stress is assumed to be proportional to the shear rate of order α introduced earlier.
Here the region of influece is H = [0, r] where the lower limit represents the axis of symmetry of the pipe r = 0. We assume that where R is the constant radius of the pipe, and introduce an extra term under the integral in (10) which is not present in definition 1 so that we can impose non-zero boundary conditions. In this case the integral in (10) is equal to the leftsided Caputo fractional derivative of order α [Mainardi(2013)] which we denote here by D α r w(r). By definition, the left-sided Caputo fractional derivative of order α is [Mainardi(2013)]: where m ∈ {1, 2, 3....}. In formula (11) the definition of the Riemann fractional integral operator of order α denoted by J α was also given [Mainardi(2013)].
For simplicity, we assume further that 1 : In particular, if α = 1 equation (10) becomes the well-known constitutive law of an incompressible viscous Newtonian fluid where µ is fluid's apparent viscosity.
To get an intuition into the physical meaning of parameter µ, we perform now a dimensional analysis for equation (10). We denote by M, L, and T the characteristic mass, length, and time scales, respectively. Then, the physical units of the quantities in equation (10) are: 1 We expect that for more complex flows other boundary conditions for d k w dr k (0 + ) could be specified.
It follows then from (10) that the physical unit of µ is: Formula (14) suggests that µ is a parameter that could for instance be represented as: where µ0 has the physical unit of apparent viscosity: [µ0] = M LT . Another representation of µ could be: where x is a preferred direction of fluid's particles aggregation during flow which may not coincide with the flow direction, and J α x is the Riemann fractional integral operator of order α.
Let assume now that µ is given by (15) and µ0 and L are constant. Then, according to (10), the shear stress varies linearly with the shear rate of fractional order α, but not with respect to the classical shear rate unless of course α = 1. A simple calculation can tell us how the shear rate (10) varries with the classical shear rate. If the classical shear rate dw dr is linear in the variable r, then w = r 2 /2 and [Mainardi(2013)]: In figure 2 we show the variations of the shear stress obtained by replacing (17) into (10) with the classical shear rate for µ0 = 1, L = 0.01 (with generic physical units) and α ∈ {0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75}. The shear thinning observed experimentally for non-Newtonian fluids including the blood flow through small vessels is obtained for 1 < α < 2 (right plot of figure 2). This suggests that the fractional order α is a parameter that provides an intinsic linkage between flow and the continuous evolution and rearrangement of fluid's microstructure during flow. When α = 1 the information about this coupling is lost and other ad-hoc assumptions are then to be made about the dependency of the apparent viscosity on the shear rate so that shear-thinning can be predicted by the constitutive law. Constitutive laws of non-Newtonian fluids that are obtained for α = 1 in the manner described here involve usually more physical parameters that need to be found experimentally. On the other hand, our constitutive model (10) has only two physical parameters, µ and α, which is a desirable feature especially when trying to observe blood flow in a living body in vivo with minimally invasive procedures.
In this paper, we assume that µ is constant 2 . Then, by replacing equations (9) and (10) into the system of equations (8) we obtain: Equations (18) and (19) suggest that dp dz = C < 0 is a given constant.
Integrating once equation (19) gives the following expression for the shear rate: where c1 is a constant of integration. One boundary condition requires the fractional shear rate (20) to be bounded at r = 0. This implies c1 = 0 and thus: We apply the operator J α to (21) and use the following properties of fractional integrals and derivatives [Mainardi(2013)]: to get: To obtain a physically consistent velocity, the second term of (22) will be found from the boundary condition at the wall r = R. We denote by τw = σrz|r=R the shear stress at the wall and by τc > 0 the critical shear stress at which slip starts. Inspired by [Ferras(2012), Ferras(2012), Kaoullas(2013)], we impose the following boundary condition: where the prescribed slip velocity v slip will be specified shortly.

Two Slip Conditions
We provide now two expressions for the slip velocity v slip . The first expression is the classic Navier slip condition [Ferras(2012), Ferras(2012), Kaoullas(2013)]: where l is a characteristic slip length and the negative sign suggests the friction between the fluid motion and the wall. The second expression that we propose here is a novel generalization of the Navier slip condition: where the power law l α is needed for dimensional consistency. We notice that if α = 1 then formula (26) reduces to (25). The slip condition (26) models the non-local entanglement and aggregation of the wall and fluid particles.
By substituting formula (24) into expression (25) and respectively formula (21) into expression (26) the following formulas for v slip are obtained: and respectively: From formulas (10) and (21) we have: τw = CR 2 and thus we define τc = C R 2 with C > 0 a constant. We now replace the expressions for τw and τc, and formulas (27), (28) into (24) and obtain the following formulas for the fluid velocity: corresponding to the Navier slip condition, and: corresponding to the generalized Navier slip condition.
Lastly, we calculate the volumetric flow rate which, in our case, is given by: By using formula (29) we get: and by using formula (30):

Non-Dimensionalization
We introduce the following non-dimensional quantities: The above formulas were adapted to our case from [Kaoullas(2013)]. The quantitiesw andQ given by (33) are dimensionless because of formula (14).
Then the dimensionless representations of formulas (29) and (30) are: and respectively: The dimensionless representations of formulas (31) and (32) are: and respectively: Formulas (34)-(37) are used in the results section.

Results
We present the effects of the pipe radius R, characteristic slip length l, and fractional order α on the fluid's velocity and volumetric flow rate.
Since the application we are interested in is cerebral blood flow and in our mathematical model we did not incorporate yet possible neuroglial effects on the vasculature, we show results for flow through an artery and an arteriole. We used R = 0.25 cm for the radius of a cerebral artery [Jain(1964)], and R = 0.001 cm for the radius of a cerebral arteriole which we estimated from the chosen arterial radius by applying Murray's law 3 to a 24-level branching tree of symmetric bifurcations [Zamir(2016)]. We used two values for the characteristic slip length: 0.0002 cm, and 0.02 cm, which are approximately 10 times less and respectively 10 times more than the diameter of an endothelial cell (about 20 microns). The results for a characteristic slip length of 0.002 cm are similar to those for l = 0.0002 cm and therefore we will not present them. The fractional order α varies from 0.1 to 1.7, and we used the values of 0.5 and 1.7 to compare the flows corresponding to the classic Navier slip condition and its generalization proposed in this paper.
In figure 3 we show the dimensionless velocity of the fluid in no slip flow. We notice that the velocity profile for α = 0.9 is close to the classic Poiseuille flow of a Newtonian fluid, while the profile for α = 1.7 resembles blood's velocity observed experimentally [Santisakultarm(2012), Zhong(2011)].  shows dimensionless velocity profiles for l = 0.0002 cm for two values of α, 0.5 and 1.7. More slip is observed in the velocity profiles satisfying the generalized Navier slip condition (bottom plot, left column of figure  4) than in the velocity profiles corresponding to the classic Navier slip condition (top plot, left column of figure 4). The first column of figure 5 shows that for l = 0.0002 cm the velocity profiles for the two slip conditions are close for α = 0.5 (top plot, left column of figure 5) and α = 1.7 (bottom plot, left column of figure 5). However, for the characterisitc slip length l = 0.02 cm big slips are observed for both slip conditions (right columns of figures 4 and 5). Not only that there are big differences between the velocity profiles for α = 0.5 and α = 1.7 for the Navier slip condition (top plot, right column of figure 4) and for the generalized Navier condition (bottom plot, right column of figure 4), but also there are big differences between the velocity profiles corresponding to two different slip conditions and the same value of α. For α = 0.5 the velocity profile satisfying the classic Navier slip condition is about 4 times bigger than the velocity corresponding to the generalized Navier slip condition (top plot, right column of figure 5), while for α = 1.7 the opposite is observed but with a much larger difference between the profiles (bottom plot, right column of figure 5). Similar differences between the classic Navier slip and generalized Navier slip conditions are observed in the velocity profiles at the pipe wall (figure 6) and the corresponding dimensionless volumetric flow rates (figure 7) with varyingC. When R = 0.25 cm andC = 2, the velocities corresponding to the classic Navier slip and to the generalized Navier slip conditions increase with decreasing α for l = 0.0002 cm and for l = 0.02 cm (figure 8), and more slip is observed for the velocity profile satisfying the generalized Navier slip condition when l = 0.02 cm (bottom plot, right column of figure 8). In fact, the difference between the velocity profiles for the two slip conditions is either very close to zero when l = 0.0002 cm (left column  Figure 6: Non-dimensional velocities calculated at the wall using formulas (34) and (35) for R = 0.001 cm. The square symbol is used for the flow with classic Navier slip, and the diamond stands for the flow satisfying the generalized Navier slip condition. The top row corresponds to α = 0, 5, while the bottom row corresponds to α = 1.7. The left column is for l = 0.0002 cm, and the right column is for l = 0.02 cm.  (36) and (37) for R = 0.001 cm. The square symbol is used for the flow with classic Navier slip, and the diamond stands for the flow satisfying the generalized Navier slip condition. The top row corresponds to α = 0, 5, while the bottom row corresponds to α = 1.7. The left column is for l = 0.0002 cm, and the right column is for l = 0.02 cm.  (34) and (35) forC = 2, R = 0.25 cm and two values of α: α = 0.5 and α = 1.7. The top row corresponds to the classic Navier slip condition, while the bottom row represents the speed for the generalized Navier slip condition. The left column is for l = 0.0002 cm, and the right column is for l = 0.02 cm.
In figure 12 we summarize our results presented above. Significant differences between the classic Navier slip and the generalized Navier slip conditions are observed only for R = 0.001 cm and l = 0.02 cm (bold letters in figure 12). The results show that a moderate increase of the dimensionless pressure gradientC causes a dramatic increase in the dimensionless volumetric flow rateQ which suggests that a persistent blood pressure above the critical shear stress when slips starts which could represent the presence of hypertension (by itself or due to diabetes) might contribute to the formation of microaneurysms since this big increase iñ Q asC increases could physically be achieved only by the dilatation of the arterioles. The validity of this statement will be investigated in our further work when we will consider a pipe with deformable walls. We also notice that both slip conditions could cause big volumetric flow rates but for different values of α. It is possible that variations of α during flow could change the slip condition at the wall. Since we envision α as a measure of entanglement not only among blood cells but also among blood and endothelial cells, α could be changed by chemo-mechanical signaling between glial cells and these entangled arterial cells. We plan to explore this aspect in the near future.

Conclusions
In this paper we proposed a mechanism for the formation, growth and rupture of microaneurysms that involve the coupled dynamics of blood, arterioles and neuroglia. The slip at the blood-vascular wall interface is seen as the initiator of this chemo-mechanical mechanism. A novel non-local constitutive law for the blood modeled as an incompressible non-Newtonian fluid is proposed and the Poiseuille flow of such a fluid  (34) and (35) forC = 2, R = 0.25 cm. The square symbol is used for the flow with classic Navier slip, and the diamond stands for the flow satisfying the generalized Navier slip condition. The top row corresponds to α = 0, 5, while the bottom row corresponds to α = 1.7. The left column is for l = 0.0002 cm, and the right column is for l = 0.02 cm.  Figure 10: Non-dimensional velocities calculated at the wall using formulas (34) and (35) for R = 0.25 cm. The square symbol is used for the flow with classic Navier slip, and the diamond stands for the flow satisfying the generalized Navier slip condition. The top row corresponds to α = 0, 5, while the bottom row corresponds to α = 1.7. The left column is for l = 0.0002 cm, and the right column is for l = 0.02 cm.  Figure 11: Non-dimensional volumetric flow rates calculated using formulas (36) and (37) for R = 0.25 cm. The square symbol is used for the flow with classic Navier slip, and the diamond stands for the flow satisfying the generalized Navier slip condition. The top row corresponds to α = 0, 5, while the bottom row corresponds to α = 1.7. The left column is for l = 0.0002 cm, and the right column is for l = 0.02 cm.
through an axi-symmetric circular pipe with rigid and impermeable walls is solved analytically. Two slip conditions are considered: the classical Navier slip and a new generalized slip condition. Both, the constitutive law and the generalized slip conditon, use spatial fractional order Caputo derivatives to express aggregation and entanglement of either only blood cells or blood and endothelial cells. Our numerical simulations suggest that hypertension could contribute to the formation of microaneurysms which agrees with some clinical observations. In our future work we plan to incorporate the deformability of the pipe due to blood flow and neuroglial chemical signaling.   Figure 12: Review of the results. Subscript 'C' stands for the classic Navier slip condition, while subscript 'F' represents the results corresponding to the generalized Navier slip condition.