Initially, the geometrical shape of the modelled area must be determined. For this purpose, one needs to create a virtual (digital) phantom ( Ω ) which is in fact represent a combination of digital areas corresponding to biological prototype. A term “phantom” is used here as a description of the object where further numerical calculations will be fulfilled. It has a direct similarity with a physical imitation of the object during ionizing irradiation research. According to biological features of NVU described above (see Section 2.), all structures have been placed one after another (Figure 2A). A capillary lumen is represented as a short tube ( Ω I ) with a circular section and the length of L_capillary. Endothelial cells line the vessel intima, and they are regarded as a cylinder ( Ω E ) with the fixed thickness of h_end. The next cover is basal lamina which is also formalized as a thin cylinder ( Ω B L ) with the width of h_bl. For simplification, the line segment of a capillary without pericytes is used for modelling.

The rest part of NVU is approximated as a cylinder ( Ω I I ) surrounding all structures described above. This cylinder is supposed to consist of different cells. Astrocytes inhere in the first layer forming the end-feet contacting with basal lamina. In the further space neurons and astrocytes compose a heterogenic structure containing cells and interstitial fluid organized in sheets and tunnels. It was anatomical proved that diffusion distances to neurons and glial cell bodies for solutes and drugs are short because no brain cell is further than about 25 μm from a capillary (Abbott et al., 2010; Gould et al., 2017). Based on the indicated observation, the thickness of tissue-cylinder Ω I I is fixed on this value (L_surround). For each type of the considered area the whole surface will be indicated as ∂ Ω and the transverse and longitudinal surfaces as S ∂ Ω and H ∂ Ω respectively. The distances between neighbor end-feet are modelled by a division of H ∂ Ω I I i n n e r on two parts. One of them corresponds to the astrocytes processes with area S e n d − f e e t . Another one represents a cleft-area where a free diffusion occurs. Due to the system has a cylindrical symmetry the division of H ∂ Ω I I i n n e r should be made by intermittent circular stripes. The set of gemetrical parameters of the virtual phantom is represented in Table 1.

For many cases, the governing equation for the fluid flow is the non-steady-state Navier-Stokes equation for an incompressible flow without buoyancy effects (Huang et al., 2013; Iasiello et al., 2016) ρ ∂ u → ∂ t + ρ ( u → ⋅ ∇ ) u → = − ∇ p + ∇ ( μ ( | γ ˙ | ) ( ∇ u → + ( ∇ u → ) Τ ) ) ; (1a) γ ˙ = 2 ε = ∇ u → + ( ∇ u → ) Τ , | γ ˙ | = 2 ( ε : ε ) ; (1b) ∇ u → = 0 ; (1c) where u → is the velocity vector along the coordinate system; ε denotes strain rate tensor; p, t and ρ are pressure, time, and the fluid density respectively. There is no slip on the boundary of the capillary lumen. The specific feature of this Eq. 1a for the non-Newtonian fluid flow is the dependence of the dynamics viscosity μ on shear rate | γ ˙ | .

There are a lot of time-dependent models proposed to describe thixotropic and viscoelastic behavior of blood (Yilmaz and Yaar Gundogdu, 2008). They designate the dependence of μ ( | γ ˙ | ) directly. It was previously shown that the non-Newtonian nature of blood acts as a regulating factor to reduce the flow resistance and wall shear stress thereby considering shear thinning to have the most significant role in facilitating blood flow through stenotic vessels (Shukla et al., 1980). Moreover, blood is usually considered as a predominantly shear thinning fluid, especially under steady flow conditions (Fatahian et al., 2018). Thus, this property seems to have the most important non-Newtonian impact (Sochi, 2013). Having quickly considered a competitive analysis of the obtained modelled values in the literature, one is able to conclude that shear thinning is evaluated using the Carreau-Yasuda model, and yield stress is usually described by the Casson model (Yilmaz and Yaar Gundogdu, 2008; Sochi, 2013). As it was mentioned above shear thinning is the most important property in the description of hemorheology and hemodynamics.

Therefore, the Carreau model has been used for Eq. 1a in the present study. According to the Carreau model the dynamics viscosity is described by the following expression: μ ( | γ ˙ | ) = μ ∞ + ( μ 0 − μ ∞ ) [ 1 + ( λ γ ˙ ) 2 ] ( n − 1 ) / 2 (2) where μ ∞ connotes a constant limit value when blood is treated as a Newtonian fluid, μ 0 is the blood viscosity at a zero shear rate, λ is the time constant associated with the viscosity that changes with the shear rate, and n is an index parameter (Molla and Paul, 2012).

Although the exact value of p under a microcirculation level is controversial, the initial pressure in the capillary has been fixed at p 0 = 18.5 mmHg, and the difference of pressure Δ p between inlet and outlet surfaces assigns CBF. The alteration in CBF is modelled by an application of a step-like function f s h i f t d e c / i n c ( t ) , which shifts both p 0 and Δ p to the modified values. If one needs to consider a decrease of CBF, f s h i f t d e c ( t ) will fall down to the fixed value which is lower than 0.5. For the example of an increased CBF the same function ( f s h i f t i n c ( t ) ) raises up to the value 1.5. p i n l e t = ( p 0 + Δ p ) f s h i f t d e c / i n c ( t ) ; p o u t l e t = ( p 0 − Δ p ) f s h i f t d e c / i n c ( t ) ; f s h i f t d e c ( t ) = { 1 , t < t 0 ; 1 − δ ( t ) , t ≥ t 0 ; f s h i f t i n c ( t ) = { 1 , t < t 0 ; 1 + δ ( t ) , t ≥ t 0 ; ∀ t : δ ( t ) ∈ [ 0 , a ] ; 0.5 ≤ a < 1 ; (3)

The view of f s h i f t d e c / i n c ( t ) is represented in Figure 3. The considered time is in a second’s range with the time-point of shift t 0 = 2   s . Indeed, the changes in a blood flow are always appeared with a continues lag phase, and a shift time-point is only a mid-point of the transition time area. This effect is modelled by a time dependent function ( δ ( t ) ) which helps to evaluate alteration of the pressure in the capillary fluently. The range of a smoothed area near the time-point of shift is characterized by Δ t . It is taken as high as 50% of the considered area because the values are shifted smoothly from initial levels to final ones.

Forming of metabolic gradients near a capillary surface obviously composes of diffusion and advection in the interstitial space which is facilitated by convection of interstitial fluid (ISF) in the paravascular space (Jin et al., 2016). To describe these processes the form of tissue structure geometry should be set explicitly. In some case it has been done on the base of 3D reconstructions of neuropil from electron microscopy images (Kinney, 2009; Mishchenko et al., 2010; Kinney et al., 2013). However, a digital form of the phantom is sometimes created based on the anatomical information. For example, to consider the glymphatic mechanism by modeling diffusive and convective transport in brain ECS, the geometry of the microvascular lobule can be idealized as a hexagonal lattice with a venule at the center, surrounded by six arterioles (Jin et al., 2016).

If the parenchymal structure is not assigned plainly the velocity of fluid flow in brain ECS can be estimated using Darcy’s law, u → = − κ μ I S F ∇ p (4) where κ and μ I S F are permeability and dynamic viscosity of the interstitial fluid, respectively. The Eq. 4 is appropriate for description of convectional processes in Ω I I because it describes a flow in a porous media. This media includes conditions where ISF is incompressible fluid which is placed around the rigid bodies depicting neurons and glia cells. The set of the physical parameters used in the modelling of CBF and ECF is represented in Table 2.

Having obtained the velocity field u → using Eqs 1a, 4, one needs to consider the diffusion of glucose in the model. The equation governing the system describes several processes. They are diffusion, convection, and glucose consumption. The type of the equation depends on the considered area. For a blood stream, endothelium and the basal lamina large stain objects forming the obstacles are supposed to be insufficient. On the contrary, for parenchyma the most essential contribution is given by spatial heterogeneity of neurons and glia cells forming a typical porous media. Thus, the gradient of metabolites will be formed in the capillary lumen and two surrounding structures on the base of the following equation: ∂ c ( r → , t ) ∂ t = ∇ ⋅ ( D ⋅ ∇ c ( r → , t ) ) − u → ⋅ ∇ c ( r → , t ) + f c o n ( c ( r → , t ) , r → ) ; (5) where c is the volume concentration of glucose, D is a diffusion tensor and f c o n ( c ( r → , t ) , r → ) is the function describing the rate of metabolic reactions which consume the diffusing molecules. There is no convection in Ω E , B L . Thus, Eq. 5 will be applied in these areas with the zero second summand in the second member of equation. In the outside space diffusion is supplied by percolation of ISF into shears and tunnels of the parenchyma. It transforms Eq. 5 into the following form: ∂ ∂ t [ α I C F c ( r → , t ) ] = ∇ ⋅ ( ( D d + α I C F τ D ) ⋅ ∇ c ( r → , t ) ) − u → ⋅ ∇ c ( r → , t ) + f c o n ( c ( r → , t ) , r → ) ; (6) where α I C F is porosity and it is also called volume fraction of ISF and τ is interstitial tortuosity. Additionally, the dispersion tensor D d for ISF has been included. It depends also on the fluid’s average velocity u → in the void space, D ( d ) i j = 1 u ∑ k , l = 1 3 a i j k l u k u l ; u = ∑ k = 1 3 u k 2 ; (7) where a i j k l is a property, called dispersivity, of the porous medium only, and u k is the k^th component of the fluid’s average velocity vector u → (Bear, 1961). Since the early 60 s, almost all research on (solute) dispersion has been limited to isotropic porous media. For such media, the components of the dispersivity tensor have been shown to depend only on two material moduli, referred to as longitudinal ( D ( d ) L ) and transversal ( D ( d ) T ) dispersivities (Fel and Bear, 2010). In a common case the coefficients of transverse and longitudinal dispersion are non-linear functions of velocity: D ( d ) L = D c τ + α L u n ;       D ( d ) T = D c τ + α T u n ;       n ∈ [ 1,2 ] (8) where the coefficients α L and α T are the longitudinal and transverse dispersivities, respectively, of the porous medium in the direction of transport, D c is a diffusion coefficient of substance c in the solvent and n is an empirically constant (Delgado, 2007).

Having assumed longitudinal dispersivity is proportional to a linear size of the system (L_surround) (Pickens and Grisak, 1981) and taking into account the medium value of ISF in approximately 5 × 10^–7 m/s (Zhan et al., 2017), one can neglect the second summand in a right part of Eq. 7. Thus, the dispersion transforms into: D ( d ) L = D ( d ) T = D c τ ; (9)

Indeed, tensors in Eqs 5, 6 are the subject for experimental measurements. Despite the attempts to depict the true diffusion process, it was generally accepted that the gradients in a biological tissue be portrayed on a voxel scale of obtained images (Basser and Jones, 2002). In that case the physical diffusion coefficient has been replaced with a global, statistical parameter, the apparent diffusion coefficient (ADC) (Chanraud et al., 2010). The non-invasive observation of the water diffusion driven displacement distributions in vivo provides unique clues to the fine structural features and geometric organization of neural tissues, and the technique of brain tissue anisotropy measurement is actually based on a water diffusion magnetic resonance imaging (Le Bihan, 2003).

Taking together the mentioned above, the elements of diffusion tensors can be represented as production of the mean diffusion coefficient and anisotropy matrix. D = D c ⋅ ( σ x x σ x y σ x z σ y x σ y y σ y z σ z x σ z y σ z z ) ; (10)

Certainly, in a common case it is difficult to indicate the sharp orientation of a short capillary part to parenchyma anisotropy. Nevertheless, without loss of generality one can suppose that σ i j = δ i j ⋅ ξ i j where δ i j is Kronecker delta and ξ i j is a dimensionless value indicating the relative ratio between the directions. The physical parameters of diffusion in distinct parts of the considered virtual phantom ( Ω ) is represented in Table 3.

The third summand in the second member of Eqs 5, 6 is corresponds to consumption of glucose in the medium. As for many metabolites this process is mediated by an activity of enzymes transforming an initial substrate to the products of reaction. It should be noticed that the first step of the metabolic pathway for glucose oxidation is considered to be a unique reaction forming f c o n ( c ( r → , t ) , r → ) . Upon glucose has come to the cells, it is transformed into glucose-6-phosphate by hexokinase. It is generally considered that this first step of glycolysis has more flux control than its transmembrane transport. The activity of hexokinase is observed both in blood stream, endothelial cells and nervous parenchyma. Although the kinetic properties of such an enzyme is known well, the activity of the whole metabolic pathway should be taken into account. The rate equation of the hexokinase–phosphofructokinase system was proposed by Heinrich and Schuster (1996), but it should be extended by a term that accounts for the glucose effect according to classical Michaelis–Menten kinetics (Aubert and Costalat, 2005). The parameters f c o n ( c ( r → , t ) , r → ) = − ε G l c / A T P Ω I / E / I I c ( r → , t ) c ( r → , t ) + K G l c ;     ε G l c / A T P Ω I / E / I I =   k G l c Ω I / E / I I c A T P Ω I / E / I I [ 1 + ( c A T P Ω I / E / I I K I , A T P ) n H ] − 1           r → ∈ Ω I / E / I I     (11)

The Eq. 11 includes not only glucose concentration but also concentration of adenosine triphosphate (ATP) - c A T P Ω I / E / I I . The kinetic equation of steady-state hexokinase reaction depends on ATP concentration explicitly. However, due to the level of such a metabolite is essentially regulated in the living cells, it remains stable in a wide range of other parameters (Ainscow et al., 2002; Nartsissov and Mashkovtseva, 2006). Thereby, in the present study ATP concentration is considered as a constant and Eq. 11 is transformed to the classical simplified hyperbolic dependence. It is also supposed that endothelial cells and neuronal parenchyma have similar kinetic properties with respect to glucose consumption, but the level of hexokinase in blood cells is lower ε G l c / A T P Ω I ≪ ε G l c / A T P Ω E ≈ ε G l c / A T P Ω I I .

For any kind of partial differential equation problem, the boundary conditions are the most essential part which determines, in fact, quantitative characteristics of the solution. Moreover, the type of condition is related to the physical properties of the system. It is supposed that the capillary is essentially isolated from the external influence, and the main part of metabolites is coming with CBF. It means that on the end transverse surface of Ω I Danckwerts condition has been set up: n → ⋅ ( J → + u → ⋅ c ( r → , t ) ) = n → ⋅ ( u → ⋅ c ˜ ( t ) ) ;     r → ∈ S ″ ∂ Ω I   (12)

Another end transverse surface was validated with outflow condition: n → ⋅ ( D ⋅ ∇ c ( r → , t ) ) = 0 ;   r → ∈ S ′ ∂ Ω I (13)

Considering a convection field in ISF, the concentrations on H ∂ Ω I I o u t e r will be governed by Eq. 13 as well. Furthermore, appropriate end transverse surfaces of Ω E and Ω B L should be considered with an no flux condition: − n → ⋅ J → = 0   ;                               r → ∈ S ′ ∂ Ω E / B L ;     r → ∈ S ″ ∂ Ω E / B L (14)

If the virtual phantom length is long i.e., L c a p i l l a r y ≫ L s u r r o u n d and c 0 ( t ) ∼ 〈 c 〉 then Eq. 14 can be replaced on the fixed boundary concentration. c ( r → , t ) = c 0 ( t ) ;     r → ∈ S ′ ∂ Ω E / B L ;     r → ∈ S ″ ∂ Ω E / B L   (15)

It should be stressed that a peculiarity of the brain tissue causes an evident effect of ISF on medium convection properties. Under considered conditions, a flow takes the metabolites away from the blood vessel boundary and it turns mathematical expression of Eq. 14 equal to Eq. 13 for r → ∈ ( ∂ Ω I I \ H ∂ Ω I I i n n e r ) . The difference between 13) and 14) will appear when the external flow of interstitial liquid is directed inside the virtual phantom. This case corresponds to veins type of a blood vessel, and it causes the substitution of Eq. 13 with Eq. 12 including an appropriate value of   c ˜ ( t ) .

The gradients of glucose will be always formed under circumstances of non-zero initial values of the metabolite. This essential feature of native compounds has a substantial difference with the same one to the toxic chemicals. Indeed, usual concentration of some mediator or a drug can be awfully close to zero because a high regulatory effect of the molecules needs to be appeared shortly. On the contrary, common participants in the metabolic pathways must be present in an essential amount in tissues all the time. It means that the initial concentrations of each metabolite are reasonably fixed at an average tissue level. c ( r → , t ) | t = 0 = 〈 c 〉 | I / E / B L / I I ;     r → ∈ Ω I / E / B L / I I (16)

In a common case the sharp values of initial concentrations are varied for different areas. However, the most diversity is observed between CBF and parenchyma.

The essential processes forming the gradients of metabolites near a capillary are the transport across the membranes. The transport of glucose is mediated by a well-known family of sodium-independent bi-directional facilitative transporters from the solute carrier 2 (SLC2) family of which 14 isoforms (GLUTs 1–14) are widely represented in endothelial cells, glia and neurons (Lizák et al., 2019). Despite variety of identified membrane carriers, GLUT1 and GLUT3 are the major glucose transporters in NVU (Patching, 2017). It is described by the usual equation of a passive transport. Under normal circumstances, as it was mentioned above brain glycolysis is not limited by glucose transport, but by phosphorylation of glucose to glucose-6-phosphate. Quantitative measurements suggest an asymmetric distribution of GLUT1 at the luminal and abluminal membranes and up to 40% of the GLUT1 protein may be sequestered within the cell cytoplasm at any given time (Patching, 2017). A number of other studies have quantified the relative amounts of GLUT1 in luminal and abluminal membranes and cytoplasm from humans and from other mammals with variable results (Deng and Yan, 2016).

The existence of the transmembrane transporters in endothelial cells and astrocytes end-feet forms an irregularity in a diffusion process. Indeed, having reached the surfaces of different areas, the molecules of glucose collide with an obstacle. There is no free diffusion available when they pass through H ∂ Ω E i n n e r and H ∂ Ω E o u t e r . For H ∂ Ω I I i n n e r there is a combination of the space with a free diffusion—the cleft between end-feet, and the membrane of astrocytes where the process of glucose entrance into the cell is the same as in endothelium. The phenomena described above form the internal flux irregularities in Ω which are introduced by the dysconnectivity functions.

These functions includes both kinetic properties of membrane transporters and their capacity on the surfaces H ∂ Ω E i n n e r , H ∂ Ω E o u t e r and H ∂ Ω I I i n n e r . It should be remarked that even though only one type of the transporter is considered but the formed fluxes are variated. Since the affinity of the protein to glucose is unchangeable under fixed conditions an obvious reason for such a diversity is the amount of GLUT incorporated into the membrane. A fluctuation of the typical value of GLUT content causes the quantitative multiplicity of the glucose transport in endothelium and end-feet. The combination of different GLUT yields the transmembrane flux of the metabolite. For each transporter, the rate equation is a typical hyperbolic function with maximal rate ( V max k ) and affinity ( K m k ).

A superposition of the transport activity is defined as the following expression: − n → ⋅ ( J → + u → ⋅ c ( r → , t ) ) H ∂ Ω E i n n e r / H ∂ Ω E o u t e r / H ∂ Ω I I i n n e r = ∑ k f G L U T ( k ) ( r → ) ⋅ V max k ⋅ c ( r → , t ) K m k + c ( r → , t ) ; | r → ∈ S e n d − f e e t       − n → ⋅ [ ( J → + u → ⋅ c ( r → , t ) ) H ∂ Ω I I i n n e r − ( J → + u → ⋅ c ( r → , t ) ) H ∂ Ω B L o u t e r ] = 0 ;       r → ∈ H ∂ Ω I I i n n e r \ S e n d − f e e t (17)

The function f G L U T ( k ) ( r → ) in Eq. 17 describes a spatial distribution of GLUTs in the compartment’s membranes. It is validated in H ∂ Ω E i n n e r , H ∂ Ω E o u t e r and H ∂ Ω I I i n n e r . The first and the second surfaces contain the transporters in the endothelium cells and the last one comprises the GLUTs in end-feet of astrocytes. Generally, it is supposed that glucose transport in both barriers before parenchyma is mediated by GLUT1 isoenzymes (Weiler et al., 2017). According to this assumption the right part of Eq. 17 is transformed to ∑ k f G L U T ( k ) ( r → ) ⋅ V max k ⋅ c ( r → , t ) K m k + c ( r → , t ) = N G L U T 1 H Ω E i n n e r / H Ω E o u t e r / H Ω I I i n n e r ⋅ k c a t 1 ⋅ c ( r → , t ) N A ⋅ ( K m 1 + c ( r → , t ) ) (18) where N G L U T 1 H Ω E i n n e r / H Ω E o u t e r / H Ω I I i n n e r is the number of carriers per μm²; k c a t 1 is the number of GLUT1 turnovers and N A is the Avogadro constant.

Finally, the kinetic constants and the parameters of transporters densities can be found in Table 4.

The model is created on the base of COMSOL Multiphysics ver. 5.5. All physical processes described above are included into structural geometry of NVU in appropriate way (Figure 2B). As it was mentioned above, the size of the blood vessel is taken as R₀ = 7 μm. The set of concentration values is the following   c ˜ ( t ) = 5   m M , 〈 c 〉 | E / B L / I I = 1   m M , 〈 c 〉 | I = 5   m M . The calculations have been evaluated using Intel^® Core ™ i9-7960X CPU 4.40 GHz and AMD Ryzen Threadripper 3990 × 64-Core Processor 4.3 GHz. The velocities fields and concentration gradients were obtained using the finite element method (FEM) applying User-Controlled Extra-Fine Meshes with minimum angle between boundaries in 240° and elemental size scaling factor 0.1. The governing equations of ISF (Eq. 4) and CBF (Eq. 1a) are solved to generate a flow field, and the velocities of both ISF and CBF are implemented into the reaction-diffusion equation with convection for transient simulations of the glucose gradients in blood flow (Eq. 5) and surrounding nervous tissue (Eq. 6), respectively. A mean calculation time was approximately 2.25 h for glucose estimations. The changes of CBF were simulated by a step-like function with a jump of pressures at t₀ = 2 s (Figure 3).

The calculations have been made for the time range from t = 0 to t = 4 s. The considered area in a time scale is chosen as 1.5   s ≤ t ≤ 2.5   s because the calculated values of blood flow velocities and glucose concentrations are changed in the smoothing area 1.75   s ≤ t ≤ 2.25   s only. They eighter decrease or increase the base value p 0 and Δ p with smoothing size of transitional zone in Δt = 0.5 s. The structure of NVU will be characterized by the ratio which describes the density of end-feet covering basal lamina. The parameter is calculated on the base of H ∂ Ω I I i n n e r division and it can be written as the following expression: δ e n d − f e e t = S e n d − f e e t S t o t a l ⋅ 100 % . (19)

The results of the modelling are both the velocity field in Ω I , I I and spatial-time distribution of the glucose concentration in all parts of the virtual phantom. If one needs to compare the results of modelling to the experimental data, the values are represented as a mean ± SD.