Toward Minimalistic Model of Cellular Volume Dynamics in Neurovascular Unit

The neurovascular unit (NVU) concept denotes cells and their communication mechanisms that autoregulate blood supply in the brain parenchyma. Over the past two decades, it has become clear that besides its primary function, NVU is involved in many important processes associated with maintaining brain health and that altering the proportion of the extracellular space plays a vital role in this. While biologists have studied the process of cells swelling or shrinking, the consequences of the NVU’s operation are not well understood. In addition to direct quantitative modeling of cellular processes in the NVU, there is room for developing a minimalistic mathematical description, similar to how computational neuroscience operates with very simple models of neurons, which, however, capture the main features of dynamics. In this work, we have developed a minimalistic model of cell volumes regulation in the NVU. We based our model on the FitzHugh–Nagumo model with noise excitation and supplemented it with a variable extracellular space volume. We show that such a model acquires new dynamic properties in comparison with the traditional neuron model. To validate our approach, we adjusted the parameters of the minimalistic model so that its behavior fits the dynamics computed using the high-dimensional quantitative and biophysically relevant model. The results show that our model correctly describes the change in cell volume and intercellular space in the NVU.


Introduction
The concept of the neurovascular unit (NVU) in neurophysiology was introduced to denote cells and their communication mechanisms involved in the autoregulation of blood supply to brain parenchyma neurons depending on their activity [1][2][3].
Intricate mechanisms of communication paths in NVU are a challenge for theoretical biology in terms of mathematical modeling. The complexity of this task lies in the fact that different types of cells are involved in the work of NVU: neurons, astrocytes, vascular smooth muscle cells, and endothelium, which communicate with each other in various ways. Nevertheless, to date, different mathematical models of processes in the NVU have been developed. In one of the first works on the topic [4], the description of the control pathway from a neuron to a blood vessel is proposed. This description is based on the nonlinear dependence of the NVU response from the potassium amount ejected by the astrocyte. In subsequent years, a group led by T. David proposed several models of increasing complexity, including almost all known mechanisms [5][6][7]. These models are computationally heavy and unsuitable for analyzing their dynamic properties since they contain about 60 ODEs and hundreds of control parameters.
Unlike synaptically connected neurons, cells in NVUs communicate mainly by releasing substances into the intercellular space, known as volume transmission [8][9][10][11]. Modeling studies that address this mechanism in neural models are also very relevant for analyzing the behavior of NVU. In particular, minimalistic models of a single neuron interacting with the intercellular space, the interaction of such neurons, and the mutual communication of neurons and astrocytes [12][13][14][15].
Recently, NVU studies received increased attention in connection with the glymphatic hypothesis-the proposed brain parenchyma drainage mechanism [16,17]. In particular, important differences were found in the NVU in awake and sleep states [18].
It has been shown that the efficiency of waste metabolites elimination is significantly facilitated during sleep, which is due to a change in the ratio between intracellular and extracellular volumes in NVU [19,20]. Thus, attention was drawn to how neurons and astrocytes change their volume following the current ion concentrations ratio inside and outside cells [21][22][23].
The mathematical modeling of the living cell volume autoregulation process based on the balance of osmotic forces has been known for a long time [24]. However, it is usually challenging to determine which of the ionic currents should be taken into account. Recently, regarding the neurovascular unit, it was shown that integrating the transmembrane flow of chlorine ions into the model can explain different patterns of volume regulation in neurons and astrocytes [25].
Therefore, unlike conventional mathematical models of neurons, the description of the NVU operation requires taking into account (1) volume transmission and (2) dynamic volume regulation. At the level of direct description of biophysical processes (transmembrane ion fluxes, osmotic forces), such models are multidimensional and dynamically intractable. For a better understanding of the possible dynamic consequences of such an extension of the model, a minimalistic description is appropriate, which would be based on the basic biophysical processes, but would represent them with the minimum number of dynamic variables and control parameters, similar to how the FitzHugh-Nagumo model [26,27] describes excitable neuron dynamics. Bearing in mind that biologically the brain parenchyma can be viewed as a 3D environment, consisting of many NVUs, the minimalistic model will allow the development of a class of models of active media specific for the brain.
In this work, we construct such a minimalistic model. In doing so, we (i) use a spiking neuron rather than a mean firing rate model; (ii) do not distinguish between the perisynaptic and the rest of the intercellular space; (iii) describe astrocytes only by potassium dynamics, whereas in a typical approach, a calcium activity is modeled. This allows us to reduce the description of the astrocyte contribution to a one-dimensional nonlinear ODE, assuming that the primary role is played by potassium uptake and that the fluxes of chlorine and sodium follow changes in the potassium balance (in reality, the situation is much more complicated, see, e.g., [28]).
Having formulated the model, we define reasonable intervals of dimensionless parameters. To do this, we rely on the results of quantitative modeling of a large physiologically based model [4][5][6][7]25] updated with our own improvements.
Further, we present the results of predicted dynamics of changes in the intercellular volume in NVU and evaluate them from the point of view of their qualitative agreement with the known physiological data. Finally, we discuss the resulting model as a whole, its limitations and perspectives in the Section 4.

Minimalistic Model of Volume Regulation in NVU
In the detailed model of a neurovascular unit, the extracellular space is divided into several parts of various sizes, as shown in Figure 1a. The synaptic connection of two neurons occurs through a narrow synaptic cleft surrounded by the perisynaptic space, which, in turn, is part of the common extracellular space (ECS). The vascular wall cells are surrounded by a relatively isolated perivascular space, which is connected to the ECS by small (about 20 nm) gaps between the astrocytic endfeet.
In our minimalistic model, we do not consider the perivascular space, cells of the vascular wall, and the presynaptic neuron, assuming that an external random signal represents it. We consider the perisynaptic space and the rest of the extracellular space as a single ECS (since there is no clear anatomical border between them). As a result, our model consists of two cellular compartments (spiking neuron SN and astrocyte AST) and the extracellular volume ECS, as shown in Figure 1b.
In this panel, the simulated system is represented as a total unit volume, and the ratio of its three parts can change with the flow of water (blue). We believe that the primary marker of this process is potassium, which is released from the neuron during the generation of spikes, can accumulate in significant amounts in the extracellular space, which is hampered by the process of reuptake and buffering of excess potassium by astrocytes. Potassium pathways are shown in the figure with red arrows. Note that in reality, the processes occurring are much more complex and involve currents of other ions. To develop a minimalistic mathematical description, we believe that potassium represents the ongoing processes, while the dynamics of other ions correlate with its changes. Such an approximation is quite obvious for a neuron since each spike generation assumes a certain scenario of opening and closing ion channels and is less obvious (but in our opinion, acceptable) for an astrocyte that does not possess the properties of electrical excitability. This issue is explained below when describing the equations of the model.  The corresponding equation set is composed of six ODEs, so there are six state variables: x, y, z, u, w n , w a , where x and y are the variables of the 2D neuron model, z stands for the amount of potassium in extracellular space, u describes the potassium amount transported into the astrocyte, and w n , w a are the current volumes of neuron and astrocyte, respectively.
We consider the total volume of the simulated space to be constant and equal to 1; therefore, the volume of ECS is simply w e = 1 − w n − w a .
When modeling the cell dynamics, it is convenient to describe the ionic flows in terms of concentrations rather than amounts. In the current study, it is less suitable since volumes that contain ions are not constant. For this reason, we use the concentration value C z = z/w e , which is calculated at each integration step.
The task of obtaining a minimalistic model for a single neuron led to the popularity of the FitzHugh-Nagumo model (FHN) [26,27]. Its main feature is the presence of two variables (fast x and slow y), which are functional counterparts of the neuronal membrane potential and generalized recovery variable. While the FHN model was initially developed to describe the generation of an impulse by a neuron, it eventually became a widely used paradigm of excitable dynamics [29]. Therefore, when developing our model, we rely on previous uses of this model [29,30] and previously considered sets of control parameters. Note also that since the FHN model was written for dimensionless time and amplitudes, the model we built on its basis was also written in dimensionless quantities, retaining only qualitative but not quantitative correspondence to the described biophysical processes.

Spiking Neuron
In our model, we use the FitzHugh-Nagumo model (FHN) supplied with the volume balance equation. The resulted equation set reads: where x and y are the fast and slow FHN model variables, which are defined by the parameters ε x , and ε y . Variable w n stands for the volume of the neuron. Parameters a, b, c, are formal FHN parameters that define the location of the linear nullcline for the variable y. We take their values from a previous study in order to provide the excitable dynamics of this model. Parameter α x scales the depolarizing action of extracellular potassium concentration C z = z/w e , while C z0 = z 0 /w e0 denotes its value at rest. Figure 2 explains the selection of parameters and the action of Cz using the nullcline crafting. Separate fulfillment of the conditionsẋ = 0 andẏ = 0 gives two curves, the relative position and points of intersection of which give information about the dynamics on the [27] phase plane. For the selected parameter values (shown in the table below), there is a single nullcline intersection and therefore a single equilibrium state. This equilibrium state is stable if the intersection point is located to the left of the minimum of the x -nullcline, which has the shape of a cubic parabola. In our model, this condition is satisfied for small values of C z , which is clearly seen in panel (b) of the figure. As C z grows, the equilibrium point loses stability, and the neuron model goes into a spontaneous oscillation mode. This is very similar to the response of a real neuron to an increase in the extracellular potassium concentration [14]. We use this analogy in our model.
Cz=6 Cz=8 x y x y (a) (b) Another important term in (1) is the fluctuating external force ξ(t) scaled by D, which is introduced to simulate the normal operation regime of a neuron. Note, while the statistics of firing events of individual neuron fits well the model of shot noise (or Poisson noise), after summation of many synaptic inputs in the dendritic tree, the resulted process becomes similar to white Gaussian noise [31][32][33]. In our model, this term provides the random activation of the system.
In the equation for w n , term w eq n (C z ) sets the current equilibrium volume that is determined by the balance of osmotic forces (Donnan equilibrium). The shift of this balance is provided by the net effect changes of all the ionic concentrations [25]. For the sake of minimalistic modeling, we assume that w eq follows the changes of extracellular potassium concentration only, since the pattern of each neuronal spike includes the specific amounts of sodium, potassium, and chloride ions that crossed the cell membrane. Thus, we assume w eq n = 1/3 at-rest extracellular potassium concentration C z0 , which increases linearly with a rise in C z , w eq n = 1/3 + α n (C z − C z0 ).

Potassium Balance in Extracellular Space (ECS)
This part of the model qualitatively and in dimensionless variables describes the balance of the potassium amount (variable z), which is released by the neuron and then redistributed between neuronal reuptake J n (C z ) and transported into the astrocyte J a (C z , u/w a ). Such a simplified way of describing the process of potassium release by a neuron was proposed and substantiated in [12,30]. Here we have expanded the description of outflow pathways both back to neuron and to astrocyte.
We calculate the current volume of extracellular space w e simply from the volumes of the neuron and astrocyte (w n and w a , respectively) based on the fact that the total volume is constant. Thus, the equations for this model compartment read: where z describes the extracellular potassium amount; α z Ψ(x) describes the release of K + by the neuron per each spike: (Ψ(x) ≈ 1 during spike and Ψ(x) ≈ 0 otherwise). Potassium reuptake by the neuron is assumed to be the linear function of C z : J n (C z ) = g n (C z − C z0 ) and vanishes at C z = C z0 . It is a clear simplification of what actually happens there but can be accepted since the large deviations of C z from its rest value are handled mostly by the astrocytic reuptake; see below.

Astrocyte Potassium Balance and Volume Regulation
In contrast to most studies on astrocyte functions modeling, we do not deal with calcium dynamics, which is mainly represented by the intercellular processes and thus does not drive the cell volume. Furthermore, we do not consider the sodium transmembrane currents. While sodium channels play an important role in neuropeptige transmission [34], as well as in neuro-metabolic coupling [35], the resting membrane potential is mainly driven by potassium currents [36]. Since astrocytes are known as not electrically excitable cells, their response to the raised extracellular potassium is expected to be in the form of a one-dimensional model with a suitable set of nonlinearities. Specifically, we focus exclusively on the function of removing excess potassium and the associated cell swelling, as in [25], where ion pumps and potassium KIR-channels are predominantly used [37]. It is also known that astrocytic reuptake of potassium nonlinearly depends on its concentration in ECS and sharply increases when it reaches a certain threshold value [38].
In our model, we describe this complex process as a piecewise linear function J a , which includes the acceleration of the process as C z grows and its limitation as the maximum transport efficiency is achieved. Specifically, we define J a as follows: where g a and g a _max are the normal and maximal uptake rates, C z_max stands for the threshold potassium concentration, min() and max() take the minimal and maximal value of two arguments, respectively. Further, J a contributes to the rise of variable u, which describes the accumulation of potassium transported into the astrocyte: Adequate assignment of w eq a is not an easy task since, in reality, the change in cell volume is caused by a shift in Donnan equilibrium and can be created by different transmembrane ion currents. We tested both the most straightforward choice (10) in the form of linear dependence on concentration and a more realistic form (11), which takes into account the transfer of potassium into the astrocyte: While (10) describes the cell swelling only phenomenologically (swelling occurs after the considerable rise of extracellular potassium), expression (11) describes the balance of net amounts of ions across the astrocytic membrane. We discuss this issue in more detail below concerning the results of quantitative modeling.

Choice of Control Parameters
Although our model is extremely simplified and dimensionless, it refers to specific physiological mechanisms that prescribe control parameter selection. Generally, it can be achieved in two ways: (i) choose the parameters to achieve reasonable dynamics of the model from the viewpoint of qualitative description, and (ii) try to fit the parameters set to describe the specific quantitative data from some independent source. In the rest of the paper, we follow both ways.
Specifically, in the first two subsections of Results, we rely on a commonly used set of parameters in the FHN model, which gives us temporal characteristics such as typical spike duration and thus the model timescale. In addition, we use a physiologically relevant scale of potassium concentration in ECS (in our model, it is C z = z/w e ) since its values have been specially studied and are well known. All other parameters were selected in order to provide reasonable behavior. Say, volume changes should not be too fast, astrocytic uptake should not follow to single neuronal spike, etc.
The resulted "generic" parameter set is given below in Table 1: The alternative approach is described in the last subsection of the Section 3 and is based on comparing variable time courses with ones obtained from the qualitative and high-dimensional biophysical model. In fact, these computed data are obtained from the computation experiment and assumed to be physiologically realistic enough to serve for the validation of our model.

Excitable Regime Shows Local, but Not Global Stability
With the parameter values we used, the FHN model has a unique and globally stable state of equilibrium, while spike generation in response to a stimulus occurs due to the presence of a pseudo-orbit-an almost closed trajectory [39]. When exposed to noise, it is the irregular escape of the system from the equilibrium point to this trajectory that creates behavior that imitates the operation of a neuron in response to stimulation from other neurons in the network. In our model, however, taking into account the presence of extracellular space changes the dynamics. For a living neuron, the release of potassium into the extracellular space, in turn, changes the degree of excitability of the neuron itself. This is due to a shift in electrochemical equilibrium, which is usually described using the Nernst potentials [24]. In our model, this effect is mimicked by the α x (C z − C z0 ) term in Equation (1). As shown in [14], this method has an effect similar to changing the Nernst potential for potassium in the classical Hodgkin-Huxley model.
Note that since our model is initially driven by noise, it is not enough to obtain one or more solutions to conclude the characteristic behavior. Rather, one needs to obtain a whole ensemble of trajectories. Thus, we repeatedly ran the model at different noise intensities but with a blocked potassium intake by astrocytes, setting g a = 0 for this. The noise was turned on for half the integration time and turned off (D = 0) in the second half. The results are shown in Figure 3. As expected, under the influence of noise, the system generates a random number of spikes, from a single one, as shown in the top panel of Figure 3a, to several, or a large train of spikes (the second and third panels from the top, respectively). An important feature is that the generation of spikes was stopped with the noise turned off not for all such launches. The third panel from the top in Figure 3a shows that after turning off the noise (to the right of the dotted line), the neuron continues to produce spikes, turning into a pacemaker. This result shows the presence of bistability-the excitable regime coexists with self-sustained dynamics. Figure 3b shows an estimate of the neuron probability to enter this state depending on the intensity of the noise. This estimate was obtained based on a series of 100 runs of the program for each value of D and counting those cases where fluctuations persisted after turning off the noise. As one can see, as the noise increases, and the system completely goes into self-sustained mode. Note that if this were just a noise-induced system switching between two states, the probability of cessation of activity after turning off the noise would not become zero. However, the variable z of the model tends to accumulate with each spike, providing positive feedback, supporting further firing. This can be clearly seen in Figure 3b, where time courses of z for all three cases are shown together.
Thus, without the removal of potassium by astrocyte, the model demonstrates the instability of the excitable regime, which does not correspond to physiologically normal behavior.

Astrocyte Action Provides Global Stability
Further, we increased the contribution of the astrocyte by increasing the parameters g a , at D = 0.07, for which the self-sustained dynamics was 100%, having also performed 100 program runs for each value of g a . As one can see in Figure 4b, even small values of the g a block self-sustained firing, returning stability to the excitable mode.  The increased contribution from astrocytic potassium reuptake rapidly prevents self-sustained neuronal firing. All displayed quantities are dimensionless. Figure 4a shows the representative time courses of the model variables for g a = 0.0005, which allow us to analyze its operation. As the top panel shows, a neuron actively generates spikes when noise is on. At the same time, both Cz (the second panel from the top) and the concentration of potassium in the astrocyte, calculated as u/w a , the lower panel, increase. Both the neuron and the astrocyte react to this activity with an increase in volume, and as a result, the intercellular volume w e decreases significantly. This, in turn, affects the value of C z . The second panel from the top shows two curves, one (blue)-for C z = z/w e , and the second-for C z = z/w e0 where w e0 = 1.3 is the initial value of the intercellular volume, as if it did not change. As can be seen, the volume change significantly increases C z , which provides positive feedback, supporting the further neuron activity. However, when the noise is turned off, any activity stops, and the system returns to its original state since astrocytic reuptake provided the stability of the excitable regime.
Thus, we can conclude that taking into account the change in cell volume generates another positive feedback loop in the model since a physiologically normal and expected decrease in the intercellular volume in itself increases the concentration of potassium, which destabilizes the state of the neuron. Note, the intensity of potassium reuptake by the astrocyte becomes a key factor stabilizing the model in the excitable mode.

Fitting Model to Quantitative Data
For a minimalistic and non-quantitative model, it is always a challenge to show that it is related to the dynamics of a real living system and captures important features of its behavior. For neuronal models, the argument is the generated signal itself-more or less similar to what is directly recorded in the experiment when recording the cell's electrical activity. In our case, the task is more complicated since not all model variables are available experimental data, and those available are usually obtained in different experiments under different conditions and cannot be joined directly. To solve this problem, we compare our model not with experimental data directly but with the computed behavior of a high-dimensional and multi-parameter model, some parts and modifications of which were justified and tested with a number of publications. The interested reader can find further details (relevant bibliography and current version of our version of such model ) in Supplementary Materials.
In this section, we present the results of the model adjustment for the results of a specific simulation run of the large model. We compare variable time courses obtained from both models by presenting them as a percentage of change from the initial level (for volumes and concentrations) or by comparing the shape (for neural spikes). Specifically, we compared the results for two cases-(i) stimulation with current from the 100th second to the first spike ( To fit the behavior of the minimalistic model to the scenarios above from the large model, we selected the parameters and modified the functions to: (i) have comparable time scales for both the dynamics of the neuron and the simulation protocol; (ii) scaled the responses of the variables so as to have a response of the same order of magnitude; (iii) reshaped the astrocytic uptake and equilibrium volume functions to bring the observed response closer to the computed quantitative data.
To do this, we completed the following: We adjusted the response of the neuron volume, without pretending to be particularly realistic since our description does not reflect the features of several biophysical mechanisms, panel (c); • Repeated the process for changing astrocyte volume, panel (d); • The noise level was adjusted to obtain a similar neuron activation process.
As a result, we obtained a pretty good match, in our opinion. In particular, we found that it is very important to adequately set the astrocytic uptake function since volume change has at least two phases: if the extracellular potassium concentration increases and the astrocyte has not yet responded to it, its volume will decrease over time until it starts to grow. This is a physical effect, but there is no data on its observation in the experiment since it is difficult to observe. As a result, we adjusted the parameters to minimize this effect, although it is slightly visible in panel (d) of the figure.
The resulting set of control parameters for the minimalistic model is shown in Table 2: With the manipulations above, we have shown that the dynamics of a simple dimensionless model do not contradict the behavior of a more complex quantitative model. However, that does not mean that we can make quantitative predictions based on a simple model, especially at another set of control parameters.

Discussion
Summing up, we have proposed a minimalistic description of the processes regulating extracellular and intracellular volumes in a neurovascular unit. We would like to stress that our goal was not to obtain quantitative predictions but to describe the dynamic properties of the process. Our model shows reasonable agreement with both basic physiological facts and the behavior of a much more complex high-dimensional quantitative model.
With a full dimension of six ODEs, our model is not easy to interpret dynamically. However, it can be made more transparent and tractable if one (i) assumes that the cell volumes are constant and maintains their initial values, and (ii) blocks the reuptake of potassium by astrocytes. Then, the model will become three-dimensional, in the variables x, y, z and nullcline planes as follows: These surfaces are shown in Figure 6. Here, the equilibrium state EP is represented by a magenta-filled point located at the intersection of all three planes. Region A, corresponding to the growth of z, is located under the light shaded part of the z-nullcline. In region B, z decreases with time. In the location of the x and y nullcline, one can recognize the analogy with Figure 2.
When a neuron is at a resting state, the trajectory sits in the vicinity of EP. Being ejected from this area, the trajectory follows the arrows in the figure and travels around the cubic nullcline, similar to the conventional FHN model. However, differently from FHN, there is the net change of z during this excursion, so the new noise-triggered firing is facilitated immediately after the trajectory approaches the intersection of green and red surfaces. This is the dynamical mechanism underlying noise-supported firing shown in Figure 3.
The decrease in extracellular volume w e observed in Figures 3 and 4 should further facilitate this effect, since it increases the concentration C z = z/w e at the same z values. In turn, the potassium reuptake by astrocyte counteracts the z rise and prevents the selfsustained firing, as shown in Figure 3.
From the consideration above, we can conclude that the extension of the neuron model by including the extracellular space leads to instability of the excitable regime. However, this instability is blocked by the action of the astrocyte. This fits well-known physiological data.
Earlier publications (e.g., [12]) have shown the similar structure of phase space in a similar (but not the same) model that implies the existence of attractors, other than the equilibrium point at neuron rest state. The modeling hint derived from this fact is that the astrocyte plays an important role in the NVU model, and models without it should be treated with caution. However, not every simple equation of potassium uptake by astrocytes is suitable.
The reason is that the growth of potassium only, induced by a neuron, leads to a decrease in volume-both of the neuron itself and of astrocytes.
The biophysically correct modeling approach is to take into account the fluxes of several ions, in particular, chlorine ions, as it was performed in the work [25].
Then, while maintaining electroneutrality, the astrocyte demonstrates a significant flow of ions inward, which leads to its swelling. In our minimalistic model, we have achieved the desired behavior by composing the potassium uptake function from three parts, which provide a realistic response of astrocytic volume.
The natural limitation of the proposed model is that our description of potassium flows is not closed since it is removed from the astrocyte to some remote storage, and from there, it is replenished in the neuron. Such a solution, however, seems to be inevitable in a l umped model. In reality, there is a process of redistributing potassium from the astrocyte in several directions, both into the perivascular space, where it affects the cells of the blood vessel wall, and laterally along the network of astrocytes, moving to the area of neighboring neurovascular units. Describing each of these processes is beyond the scope of this work and deserves the development of specialized models.

Conflicts of Interest:
The authors declare no conflict of interest.