# Endothelial calcium dynamics elicited by ATP release from red blood cells

Ananta Kumar Nayak, Sovan Lal Das, Chaouqi Misbah

► To cite this version:

Ananta Kumar Nayak, Sovan Lal Das, Chaouqi Misbah. Endothelial calcium dynamics elicited by ATP release from red blood cells. Scientific Reports, 2024, 14 (1), pp.13550. 10.1038/s41598-024-63306-2. hal-04798161

HAL Id: hal-04798161

https://hal.science/hal-04798161v1

Submitted on 22 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Endothelial Calcium Dynamics Elicited by ATP Release from Red Blood Cells

Ananta Kumar Nayak<fnref n="1" />, Sovan Lal Das<fnref n="2" />, and Chaouqi Misbah<fnref n="1" />,*

## ABSTRACT

Red blood cells (RBCs) exhibit an interesting response to hydrodynamic flow, releasing adenosine triphosphate (ATP). Subsequently, these liberated ATP molecules initiate a crucial interaction with endothelial cells (ECs), thereby setting off a cascade involving the release of calcium ions ($Ca^{2+}$). $Ca^{2+}$ exerts control over a plethora of cellular functions, and acts as a mediator for dilation and contraction of blood vessel walls. This study focuses on the relationship between RBC dynamics and $Ca^{2+}$ dynamics, based on numerical simulations under Poiseuille flow within a linear two-dimensional channel. It is found that the concentration of ATP depends upon a variety of factors, including RBC density, channel width, and the vigor of the flow. The results of our investigation reveals several features. Firstly, the peak amplitude of $Ca^{2+}$ per EC escalates in direct proportion to the augmentation of RBC concentration. Secondly, increasing the flow strength induces a reduction in the time taken to reach the peak of $Ca^{2+}$ concentration, under the condition of a constant channel width. Additionally, when flow strength remains constant, an increase in channel width corresponds to an elevation in calcium peak amplitude, coupled with a decrease in peak time. This implies that $Ca^{2+}$ signals should transition from relatively unconstrained channels to more confined pathways within real vascular networks. This notion gains support from our examination of calcium propagation in a linear channel. In this scenario, the localized $Ca^{2+}$ release initiates a propagating wave that gradually encompasses the entire channel. Notably, our computed propagation speed agrees with observations.

## Introduction

Red blood cells (RBCs) are commonly known to transport oxygen to the tissues and to carry metabolic wastes from the tissues<fnref n="1" />. However, RBC also transport and release an important signaling molecule, adenosine triphosphate (ATP)<fnref n="2" />⁻⁴. ATP is an important second messenger, which helps in transmitting information from the surrounding of endothelial cell (EC) into cytoplasm of the EC. This is only possible by activating a cascade of reactions on plasma membrane (PM) of the EC⁵,⁶. A basic consequence of ATP is the release of sequestered $Ca^{2+}$ from the endoplasmic reticulum (ER) of EC into the cytoplasm⁷.

The endothelium, a confluent monolayer of ECs, forms the inner side of all vessel walls, so that it is constantly come in contact with the bloodstream. A basic function of the EC is to selectively allow the transportation of oxygen, nutrients, and ions to the tissue for a smooth functioning of organs. In addition, ECs synthesize one of the important vasodilators, namely the Nitric Oxide (NO), from the precursor L-arginine⁸,⁹. The production of NO is controlled by the endothelial nitric oxide synthase (eNOS) enzyme. However, the activation eNOS is regulated by the cytoplasmic $Ca^{2+}$ concentration⁴,¹⁰⁻¹². Thus, $Ca^{2+}$ is treated as an important ion and a second messenger that controls many cellular functions at different time scales¹³. $Ca^{2+}$ is also well known to regulate several upstream and downstream pathways such as PLC-dependent $IP_3$ production due to increase in the cytoplasmic $Ca^{2+}$ concentration¹⁴ and calmodulin-eNOS mediated NO synthesis¹⁵, respectively. Overall, $Ca^{2+}$ ion implicitly controls many cellular functions such as gene transcription, fertilization, and vasodilation¹³.

NO molecules produced due to the activation of eNOS by $Ca^{2+}$ rapidly diffuse into the neighboring ECs, the bloodstream, and the vascular smooth muscle cells (SMCs). This leads to the relaxation of SMCs¹⁶, reducing thereby the vascular resistance to blood flow by dilating the diameter of the vascular wall locally. In other words, NO helps in maintaining vascular tone within the vascular system¹⁷. In addition, NO in the bloodstream ensures several other functions such as prevention of the adherence and aggregation of platelets to the EC, and the adhesion of leukocyte and infiltration, as well as the proliferation of SMCs. Furthermore, any imbalance in NO production leads to disruption of vascular homeostasis. As a result, several life-threatening vascular diseases such as aortic aneurysms, hypertension, and atherosclerosis may arise¹⁷,¹⁸. An essential first step is to elucidate the early process leading to the rise of $Ca^{2+}$ in the cytoplasm due to ATP released from RBCs in the vascular vessels. This is the first objective of this paper.

The role of ATP in the rise of cytoplasmic $Ca^{2+}$ concentration is quite well understood. This is due to the activation of

<footnote n="1">¹Université Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France</footnote>

<footnote n="2">²Physical and Chemical Biology Laboratory, and Department of Mechanical Engineering, Indian Institute of Technology Palakkad, 678623 Palakkad, India</footnote>

<footnote>*chaouqi.misbah@univ-grenoble-alpes.fr</footnote>
purinergic receptors (P2Y₂) by ATP⁵,⁶. For example, binding of ATP to P2Y₂ receptors on the surface of ECs leads to inception of a second messenger called inositol trisphosphate (IP₃), which eventually helps in mobilization of sequestrated Ca²⁺ from ER¹⁹. The main source of ATP in the bloodstream originates from RBCs². These ATP molecules help in mobilization of ER Ca²⁺. However, another source of ATP originates from the ECs, directly due to shear stress, triggering thus Ca²⁺. ATP (from RBC and EC) activates G-protein (which is linked to purinergic receptors) leading to the production IP₃ molecules and then to the release of Ca²⁺ from the ER²⁰.

The typical ATP concentration inside RBC is in the range of mM²¹, whereas it is in the range of µM in the plasma²². This means that ATP released from RBCs may potentially exceed the need of ATP in the plasma, and this may become physiologically harmful²³. Therefore, ATP concentrations in the bloodstream must be maintained within the acceptable physiological range in the plasma (i.e., of the order of µM). This regulation is achieved thanks to an enzyme, ectonucleotidase, present on the surface of the EC, which hydrolyzes ATP into adenosine diphosphate (ADP) and eventually, ADP to adenosine monophosphate (AMP) and AMP to adenosine (AD). This means that there is a need to take into account the role of the EC in maintaining the physiological level of ATP in the plasma. Thus, we will consider here ATP release from both RBCs (the major source) and EC, as well as hydrolysis of ATP by the EC.

The elucidation of ATP pathways in the release by RBCs has improved thanks to in-vitro experiments². It has been reported that ATP released from RBCs follows two pathways: (1) the activation of Pannexin-1 (Px₁) hemichannels due to mechanical stress experienced by the membrane of RBC in flow and (2) the stretching of RBC membrane (due to topological defects in the spectrin-actin network) leads to activation of the cystic fibrosis transmembrane conductance regulator (CFTR), which ultimately up-regulates the Px₁ hemichannels. Zhang et al.²⁴ have recently proposed a phenomenological model based on the above experimental results². Their model accounted for the above two pathways and agreed with the experimental observations². However, the experimental results of Forsyth et al.² and the numerical study of Zhang et al.²⁴ were performed for a very dilute suspension of RBCs under a shear flow. The realistic blood flow situation in circulatory system is governed by the Poiseuille flow and the complex hydrodynamic interactions among RBCs, where the total volume of RBCs in a healthy situation may reach up to 45 % (designated as Hematocrit, *Hct*). Gou et al.²⁵ extended the study of Zhang et. al²⁴ to the Poiseuille flow in a two-dimensional (2-D) straight channel for different flow conditions as well as for different RBCs hematocrits. However, the role of endothelium (i.e., ATP release from ECs, hydrolysis of ATP by the ectonucleotidase enzymes, and the activation of Ca²⁺ signaling), were not studied. As a result, ATP concentrations in the channel increased with time without bound. Clearly, the role of the ECs must be incorporated into the model in order to maintain a stable ATP concentration. This is another goal of this study. Once this question is elucidated, we will be in a position to address the question of Ca²⁺ signaling in the endothelium.

Our next goal will be the analysis of calcium signaling due to ATP. In order to understand how the distribution of ATP on the surface of the endothelium affects the Ca²⁺ signaling in the endothelium, Plank et al.²⁶,²⁷ proposed a model which couples the main three variables, namely IP₃ (IP₃), the cytoplasmic Ca²⁺ (*Ca*<sub>c</sub>), and the ER Ca²⁺ (*Ca*<sub>s</sub>) with a given distributions of ATP on the endothelium surface. Although they considered the role of endothelium on the steady concentrations of ATP distribution in the 2-D channel, the source of ATP was not known *apriori*. In addition, the Ca²⁺ dynamics model used in these studies did not satisfy intracellular homeostasis. The intracellular homeostasis means that the concentrations of three important variables such as IP₃, *Ca*<sub>c</sub>, and *Ca*<sub>s</sub> must return to their initial concentrations in the presence of ATP. To fill this gap, we have recently developed a minimal model in which the intracellular homeostasis was achieved thanks to the phosphorylation of P2Y₂ receptors upon binding of ATP molecules²⁸.

In this study, we will treat RBCs under flow as the first main source of ATP. This source will be shown to depend on the spatial organization of RBCs, which is affected by the channel width, the hematocrit and the flow strength. In addition, we will take into account the role of ECs (i.e., ATP release from ECs, hydrolysis of ATP by the ectonucleotidase enzymes, and the activation of Ca²⁺ signaling), which is an essential ingredient in order to achieve a steady pattern of ATP. Finally, we will use our recently developed model²⁸ for Ca²⁺ signaling, which takes into account the phosphorylation of P2Y₂ receptors upon binding of ATP molecules. This ingredient will be shown to ensure homeostasis for the three variables entering dynamics (IP₃, *Ca*<sub>c</sub>, and *Ca*<sub>s</sub>).

We have studied Ca²⁺ signaling as a function of hematocrit, channel size, and flow strength. Additionally, we have analyzed the possibility for Ca²⁺ propagation along the endothelium and determined the propagation speed. Our investigation reveals key findings: (1) The peak amplitude of Ca²⁺ in EC increases proportionally with RBC concentration. (2) Higher flow strength increases the magnitude of Ca²⁺ signals and reduces the time to reach the Ca²⁺ peak, for a given channel width. (3) Increasing channel width, while flow strength remains constant, elevates Ca²⁺ peak amplitude and reduces peak time. This suggests that Ca²⁺ signals have the potential to propagate from a region of higher Ca²⁺ (wide vessels) concentration to a region of lower Ca²⁺ (small vessel diameters) concentration within vascular networks at a speed of 5 µm/s. This study will offer a general basis for Ca²⁺ signaling in the presence of blood flow along ECs, which can serve for a more elaborate analysis in complex vascular networks.
Methods

We consider the 2-D model for the sake of numerical efficiency, as we have already shown that the 2-D model captures the essential features of the three-dimensional (3-D) RBC model (such as shapes under flow<sup>29</sup>, cross-streamline noninertial migration velocity<sup>30,31</sup>, inclination angle of cell with respect to flow direction<sup>32,33</sup>, multilobe shapes<sup>34</sup>, and ATP release under flow<sup>25</sup>). Note that for ATP release, a topic of direct interest here, an even semi-quantitative agreement between 2D and 3D is obtained<sup>25</sup>. Due to the 2-D assumption, the cell boundary is one-dimensional (1-D) incompressible boundary and can only undergo bending (like for pure phospholipid membranes). The fluid flow inside and outside the RBC is described by the Navier-Stokes equations for an incompressible flow

$$
\nabla \cdot \mathbf{v} = 0, \tag{1}
$$

$$
\rho \left( \frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} \right) = -\nabla p + \nabla \cdot (\mu \nabla \mathbf{v}) + \mathbf{F}, \quad (2)
$$

where ρ is the density of the incompressible fluid, p is the pressure, and μ is the dynamic viscosity of the fluid. In order to distinguish the dynamic viscosities of fluids in the interior and exterior of RBCs, we define a dimensionless number called the viscosity contrast, λ. λ is defined as the ratio between the dynamic viscosity of the interior, μ<sub>in</sub> to that of the exterior, μ<sub>ex</sub> fluids of the RBC. In this numerical study, we use the physiological value of the viscosity contrast and the dynamics viscosity of the fluid present inside the RBC i.e., λ = 6 and μ<sub>in</sub> = 0.012 Pa s<sup>2</sup>. The force, **F** (N/m<sup>3</sup>) in equation (2) is obtained from the deformation of the membrane due to flow. This force is distributed over fluids close to the RBC membrane using the immersed boundary method (IBM)<sup>35</sup>

$$
\mathbf{F} = \oint_{\Gamma} \mathbf{f}(s) \delta(\mathbf{x} - \mathbf{X}(s)) ds, \qquad (3)
$$

where $\mathbf{x} \in [0, L] \times [-W/2, W/2]$ ($L$ and $W$ are the length and width of the channel, respectively), $\Gamma$ represents the integration along the contour of a RBC (shown in Fig. 1), and $\mathbf{X}(s)$ is a set of Lagrangian points on the membrane in the x-y plane, $s$ is being the local arc length. This membrane force, $\mathbf{f}(s)$ is obtained by the functional derivative of the Helfrich bending energy for the vesicle model. The Helfrich bending energy per unit length (J/m), $E$ is defined as<sup>36</sup>

$$
E(\mathbf{X}(s)) = \frac{\kappa}{2} \oint_{\Gamma} c^2 ds + \oint_{\Gamma} \zeta ds, \quad (4)
$$

where $\kappa$ is the bending modulus of the membrane ($\kappa = 3 \times 10^{-19} \text{ J})^{37}$, $c$ is the local mean curvature of the membrane, and $\zeta$ is the Lagrangian multiplier, which preserves the constant membrane perimeter, $P$. The membrane force, $\mathbf{f}(s)$ (N/m$^2$) is obtained after taking the functional derivative of equation (4)$^{30}$ and is written as

$$
\mathbf{f}(s) = \frac{\delta E}{\delta \mathbf{X}(s)} = \kappa \left( \frac{\partial^2 c}{\partial s^2} + \frac{c^3}{2} \right) \mathbf{n} - \zeta c \mathbf{n} + \frac{\partial \zeta}{\partial s} \mathbf{t}, \quad (5)
$$

where **n** and **t** are the local unit outward normal and tangent vectors on the contour, respectively. The velocity of the Lagrangian points is obtained by interpolating the velocity of neighboring fluid nodes to the membrane using the IBM<sup>35</sup>, as follows

$$
\frac{\partial \mathbf{X}}{\partial t} = \mathbf{V}(s) = \int_{y} \int_{x} \delta(\mathbf{x} - \mathbf{X}(s)) \mathbf{v}(\mathbf{x}) dxdy. \qquad (6)
$$

Subsequently, the position of the Lagrangian points of the membrane is updated using Euler method. We define a dimensionless number called the reduced area, $v$, which is defined as the ratio between the vesicle area, $A$ and area of a circle having the same perimeter ($P$) i.e.,

$$
v = \frac{4\pi A}{p^2}
$$

where $A = \pi R_0^2$, $R_0$ is the characteristics radius of the RBC. In numerical simulation, we have taken, $v = 0.65$, which is the physiological value for a healthy RBC. We impose a pressure-driven Poiseuille flow, with a velocity field along the x-direction given by

$$
v_x = v_{\max} \left[ 1 - \left( \frac{y}{W/2} \right)^2 \right].
$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic shows: (a) RBC is moving in a straight two-dimensional channel (2-D) which is lined with a layer of endothelial cells known as endothelium; (b) the pathways by which RBC releases ATP molecules when it is subjected to shear stresses; (c) the activation of pannexin-1 (Px₁) hemichannels when the membrane shear stress exceeds the critical shear stress; and (d) up-regulation of ATP release by the cystic fibrosis transmembrane conductance regulator (CFTR) through Px₁ when the curvature change rate exceeds the critical curvature change rate.</figcaption></figure>

where $v_{max}$ is the maximum flow velocity at the center of the channel. We assume no slip boundary conditions for the channel walls. The immersed boundary lattice Boltzmann method (IBLBM)³⁸ is employed to numerically solve the equations (1)-(6) for the fluid velocity field.

Once the velocity field is obtained in the channel, it is used to find the distribution of ATP in the channel. The ATP released from RBCs due to shear stress and cell deformation, diffuses in the suspending fluid and is advected by the flow. Thus, the advection-diffusion equation for ATP is

$$ \frac{\partial c_{ATP}}{\partial t} + \mathbf{v} \cdot \nabla c_{ATP} = D_{ATP} \nabla^2 c_{ATP}, \qquad (7) $$

where $D_{ATP}$ is the diffusion coefficient of ATP in the plasma, and $c_{ATP}$ is the concentration of ATP in the plasma. Note that we only solve equation (7) for ATP concentration outside RBCs. The concentration of ATP in the cytoplasm of RBCs (~ 1 mM) is several folds larger than the ATP concentration in the plasma (~ 1 µM). Therefore, we consider the RBC as a large reservoir of ATP²⁵. The detailed model of ATP release is given in following sections.

## ATP release model

According to Zhang et al.²⁴, the release of ATP from RBCs follows two pathways (shown in Fig. 1): (1) the Px₁ hemichannels are activated when the local membrane shear stress exceeds a critical value and (2) the membrane deformation leads to activation
of the CFTR, which up-regulates the $Px_1$ hemichannels. These two ATP pathways are mathematically expressed as

$$p_{\sigma} = k_{\sigma} H\{\sigma_{mem}(\mathbf{X}(s)) - \sigma_{c}\}, \qquad (8)$$

and

$$p_{c} = \min[1 + k_{c}\{\dot{c}(\mathbf{X}(s)) - \dot{c}_{c}\}H\{\dot{c}(\mathbf{X}(s)) - \dot{c}_{c}\}, p_{c\max}], \qquad (9)$$

where $p_σ$ and $p_c$ are ATP release rates due to opening of $Px_1$ hemichannels and the CFTR amplification factor, respectively. $H(\cdot)$ is the Heaviside step function, $k_σ$ is a phenomenological coefficient, $\sigma_{mem}$ is the local membrane shear stress, and $\sigma_c$ is the critical shear stress required for opening of ATP channels. $k_{c}$ is the slope of local membrane curvature change rate, $\dot{c}$ is the local curvature change rate (quantifying the degree of cell deformation), $\dot{c}_c$ is the critical curvature change rate for up-regulation of $Px_1$, and $p_{c\max}$ is the maximum amount of CFTR. We define the flux on the outer surface of RBC²⁴ as

$$\psi = p_{\sigma}p_{c}, \qquad (10)$$

which postulates that no ATP is released from the RBC when $\sigma_{mem} < \sigma_c$ and is not up-regulated by cell deformation (equation (9)) when $\sigma_{mem} < \sigma_c$.

## Boundary conditions

The ATP release flux on the RBC membrane takes the form

$$D_{ATP} \frac{\partial c_{ATP}}{\partial \mathbf{n}} = -\psi, \qquad (11)$$

where $\mathbf{n}$ denotes the local normal vector pointing outward from the boundary to fluids, and the definition of $\psi$ is given in equation (10). Because of the strategic location of the endothelium (between the vascular lumen and tissues), it plays a critical role in maintaining the vascular tone. One of the important functions of EC is to maintain the concentration of ATP in the physiological range. This is achieved with help of ectonucleotidase enzymes. These enzymes act as a safeguard against unwanted activation of pathways or reactions with receptor due to the presence of abnormal ATP concentration²³. Another function of EC is to sense the fluid flow. It is then transmitted to the cytoplasm for synthesis and liberation of signaling molecules such as ATP³⁹. Several studies have been devoted⁴⁰⁻⁴² to the understanding the role shear stress on the rise of cytoplasmic $Ca^{2+}$ concentration. One of the mechanisms may be the autocrine signaling³⁹,⁴³,⁴⁴ i.e., EC releases ATP in response to the flow and these ATP molecules bind to the surface receptors of the same EC to initiate the $Ca^{2+}$ signaling reactions. These two functions of EC are mathematically expressed as the first and second terms in equation (12)⁴⁵. The ATP flux on the surface of endothelium is written as

$$D_{ATP} \frac{\partial c_{ATP}}{\partial \mathbf{n}} = \frac{V_{\max} c_{ATP}}{K_m} - S_{\max} \left[ 1 - \left( \frac{-\tau_w}{\tau_0} \right)^3 \right], \qquad (12)$$

where $V_{max}$ is the maximum enzymatic reaction velocity, $K_m$ is the Michaelis-Menten constant, $S_{max}$ is the maximum ATP release rate from the EC, $\tau_w$ and $\tau_0$ are the wall and reference shear stresses on the EC surface, respectively.

Once the velocity field is obtained by solving the Navier-Stokes equations (equations (1)-(2)), we use it to solve for ATP concentration field using the advection-diffusion equation (equation (7)). The advection-diffusion equation (equation (7)) is supplemented with boundary condition equations (11)-(12). Both fluid and ATP equations are solved numerically using an in-house lattice Boltzmann method (LBM) solver. The detailed parameter values for the ATP release model are given in Table 1, and these values are obtained from² and⁴⁵. The details of the numerical procedure can be found in²⁴,⁴⁶.

## Dimensionless numbers and some mathematical definitions

Dimensionless numbers are used to characterize the flow, the ATP distribution, and the ATP release level of RBCs.

1. The capillary number is defined as the ratio between the viscous stress of fluid and bending strength of the membrane

$$Ca = \frac{\mu_{ex} \gamma_{w} R_{0}^{3}}{\kappa}. \qquad (13)$$

where $\dot{\gamma}_w = 4v_{max}/W$ refers to the wall shear rate of the Poiseuille flow.

2. The confinement number is used to define the degree of confinement of a RBC by the gap of the channel. It is defined as

$$C_n = \frac{2R_0}{W}. \qquad (14)$$
3. To understand the role of ECs on the lumen ATP concentration, we define the spatial-averaged ATP concentration for a given time

$$ \bar{c}_{\text{ATP}} = \frac{1}{LW} \int_{0}^{L} \int_{-\frac{W}{2}}^{\frac{W}{2}} c_{\text{ATP}}(x, y, t) dy dx. \qquad (15) $$

This will allow us to follow the global ATP concentration as a function of time in order to examine the evolution towards a steady state for ATP. The time-averaged ATP concentration along the channel is defined as

$$ <c_{\text{ATP}}>_{\text{transverse}} = \frac{1}{T} \int_{0}^{T} \frac{1}{L} \int_{0}^{L} c_{\text{ATP}}(x, y, t) dx dt. \qquad (16) $$

This quantity will allow us to analyze the ATP spatial distribution in the y-direction, orthogonal to the flow direction. The time-averaged ATP concentration at the wall is defined as

$$ <c_{\text{ATP}}>_{\text{wall}} = \frac{1}{T} \int_{0}^{T} c_{\text{ATP}}(x,t) |_{\text{wall}} dt. \qquad (17) $$

This is an important quantity since it is defined at the wall where ATP reacts with EC.

4. In a 2-D channel, the hematocrit is defined as the area occupied by RBCs to the area of the channel

$$ Hct = \frac{N_c A}{LW}. \qquad (18) $$

<table><caption>**Table 1.** ATP release parameters</caption><thead><tr><th>Parameters</th><th>Description</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td>k<sub>σ</sub></td><td>Phenomenological constant</td><td>1.75 × 10<sup>2</sup> nM µm s<sup>-1</sup></td><td>25</td></tr><tr><td>σ<sub>c</sub></td><td>Critical shear stress</td><td>0.05 Pa</td><td>24</td></tr><tr><td>k&#x0307;<sub>c</sub></td><td>Slope of curvature change rate</td><td>6 × 10<sup>-3</sup> µm s</td><td>24</td></tr><tr><td>c&#x0307;<sub>c</sub></td><td>Critical curvature change rate</td><td>2 × 10<sup>2</sup> µm<sup>-1</sup>s<sup>-1</sup></td><td>24</td></tr><tr><td>p&#x0307;<sub>cmax</sub></td><td>Maximum amount of CFTR</td><td>2.5</td><td>24</td></tr><tr><td>D<sub>ATP</sub></td><td>ATP diffusion constant</td><td>2.36 × 10<sup>-10</sup> m<sup>2</sup>/s</td><td>24</td></tr><tr><td>V<sub>max</sub></td><td>Enzymatic reaction velocity</td><td>0.8 µmol m<sup>-2</sup> s<sup>-1</sup></td><td>45</td></tr><tr><td>K<sub>m</sub></td><td>Michaelis-menten constant</td><td>475 µM</td><td>45</td></tr><tr><td>S<sub>max</sub></td><td>Maximum ATP release rate</td><td>10<sup>-3</sup> nmol m<sup>-2</sup> s<sup>-1</sup></td><td>45</td></tr><tr><td>τ<sub>0</sub></td><td>Reference shear stress</td><td>1 Pa</td><td>45</td></tr></tbody></table>

## Ca<sup>2+</sup> dynamics model of endothelial cell signaling

Recently, we have developed a model for calcium that ensures return to homeostasis in the presence of agonists<sup>28</sup>. The schematics of the model are recalled in Fig. 2. The basic new element (in comparison to previous literature where homeostasis was not achieved) in that model is the introduction of ATP-receptors desensitization. The homeostatic model is described by 5 important variables<sup>28</sup>, namely unphosphorylated receptor ([R<sub>s</sub>]), phosphorylated receptor ([R<sub>p</sub>]), IP<sub>3</sub> ([IP<sub>3</sub>]), cytoplasmic Ca<sup>2+</sup> ([Ca<sub>c</sub>]), and store Ca<sup>2+</sup> ([Ca<sub>s</sub>]). The evolution of these variables are obtained based on balance of fluxes across the membrane (PM and ER) as well as based on the chemical kinetics, and read

$$ \frac{d[R_s]}{dt} = k_r[R_t] - \left( k_r + k_p \frac{[L]}{K_r + [L]} \right) [R_s] - k_r[R_p], \qquad (19) $$

$$ \frac{d[R_p]}{dt} = k_p \frac{[L][R_s]}{K_r + [L]} - k_e [R_p], \qquad (20) $$

$$ \frac{d[IP_3]}{dt} = q_{\text{gen}} - q_{\text{deg}}, \qquad (21) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic shows that red blood cells (RBCs) and endothelial cells (ECs) release adenosine triphosphate (ATP) molecules due to the application of parabolic-Poiseuille flow along x-direction. ATP is then transported to the wall-lined with ECs where ATP reacts with the surface receptors, in order to provoke the Ca²⁺ signaling. Here, both walls of the vascular lumen are lined with ECs. The activation process of Ca²⁺ signaling in an EC is shown in the enlarged Figure. This enlarged Figure depicts the activation of cascade of reactions and influx/efflux of Ca²⁺ fluxes, across the plasma membrane (PM) and the endoplasmic reticulum (ER) membrane, due to the binding of ATP molecules with the purinergic receptors (P2Y₂). Words starting with q represent fluxes and expression for these fluxes are reported in Ca²⁺ signaling subsection. The dynamics of P2Y₂ receptors is shown in the red colored enclosure. There are four states of a receptor: (a) unphosphorylated receptors, R, (b) ligand-receptor complexes, LR, (c) phosphorylated receptors, Rₚ, and (d) internalized receptors, Rᵢ. ATP concentration is denoted by L. k₁⁺, k₁⁻, kₚ, kₑ, and kᵣ are the chemical reaction constants. Abbreviations: inositol trisphosphate (IP₃), inositol disphosphate (IP₂), inositol tetrakisphosphate (IP₄), G-protein (G), phospholipase C (PLC), phosphatidylinositol 4,5-bisphosphate (PIP₂), and diacylglycerol (DAG).</figcaption></figure>

$$ \frac{d[Ca_c]}{dt} = \beta (q_{rel} - q_{serca} + q_{leak} + q_{soc} + q_{in} - q_{pmca}), \qquad (22) $$

$$ \frac{d[Ca_s]}{dt} = v_r (q_{serca} - q_{rel} - q_{leak}), \qquad (23) $$

where [L] (or is equivalent to the c<sub>ATP</sub> at the boundary) is the ATP concentration on the endothelium surface and [R<sub>t</sub>] is the total receptor concentration. Furthermore, the fluxes given in equations (21)-(23) are defined below. The generation of IP₃
flux due to hydrolysis of the PIP₂ molecules by the PLC-enzymes is denoted by $q_{gen}$, and the degradation flux of IP₃ to inositol bisphosphate (IP₂) and inositol tetrakisphosphate (IP₄) due to the action of enzymes such as phosphates and kinases is denoted by $q_{deg}$. Other fluxes such as $q_{rel}$, $q_{serca}$, $q_{leak}$, $q_{soc}$, $q_{in}$, and $q_{pmca}$ are release of Ca²⁺ from ER, uptake by the ER from cytoplasm, passive leak from ER, influx of the extracellular Ca²⁺ to cytoplasm via the store-operated channel (SOC), passive leak of Ca²⁺ from the extracellular space through PM, and pumping of cytoplasmic Ca²⁺ to the extracellular space using PMCA, respectively. The value and meaning of each constant in equations (19)-(20) as well as in the fluxes of equations (21)-(23) are extracted from²⁸.

$$q_{gen} = \gamma \frac{[Ca_c]}{K_1 + [Ca_c]}, \qquad (24)$$

where $\gamma = \alpha[G] [(PIP_2)_t] / (N_A V_{EC})$, $N_A V_{EC}$ is the product of the Avogadro's number and the volume of an EC.

$$[G] = [G_t] \frac{\rho_r + \delta}{\rho_r + \delta + K_g} \equiv [G_t]\lambda, \qquad (25)$$

where $K_g = k_d/k_a$,

$$q_{deg} = k_2 [IP_3], \qquad (26)$$

$$q_{soc} = k_{soc} \frac{1}{\left(1 + \frac{[Ca_s]}{K_{soc}}\right)^n} [Ca_{ex}] ([Ca_s](0) - [Ca_s]), \qquad (27)$$

$$q_{pmca} = k_8 \frac{[Ca_c]^2}{K_4^2 + [Ca_c]^2}, \qquad (28)$$

$$q_{leak} = k_5 [Ca_s], \qquad (29)$$

$$q_{serca} = k_4 \frac{[Ca_c]^2}{K_3^2 + [Ca_c]^2}, \qquad (30)$$

$$q_{rel} = k_3 \frac{[IP_3]^3}{K_2^3 + [IP_3]^3} \frac{K_{i,Ca_c}^3}{[Ca_c]^3 + K_{i,Ca_c}^3} Ca_s, \qquad (31)$$

$$q_{in} = k_8 \frac{([Ca_c](0))^2}{K_4^2 + ([Ca_c](0))^2}, \qquad (32)$$

$$\beta = \left\{ 1 + \frac{[B_t]^{\frac{k_7}{k_6}}}{\left( \frac{k_7}{k_6} + [Ca_c] \right)^2} \right\}^{-1}, \qquad (33)$$

where $\beta$ is a buffer factor. The buffer factor is obtained using the rapid equilibrium approximation i.e., the binding and unbinding reactions of free Ca²⁺ and buffer proteins are governed by the fast-kinetics.
10(a) in Guglietta et al.⁵¹ and Fig. 2 in Guglietta et al.⁴⁹). Given the fact that our simulation of suspension requires about 2.5 s to reach a steady-state, we expect the corresponding time in the presence of membrane viscosity to be of about 4 s, significantly smaller than the time required for the development of calcium peak (which lies in the 10-30 s range). Therefore the separation of time scales (that corresponding to RBCs reaching steady-state, and that corresponding to calcium kinetics) seems to be reasonably legitimate. For this reason, we first solved for blood flow in order to determine the ATP distribution on the walls after RBCs have reached their steady-state configuration. We then inject this information into the calcium model. We repeated the same procedure for each set of parameters such as capillary number, confinement, and hematocrit.

Note that the magnitude of $Ca^{2+}$ concentrations in the blood lumen is several orders greater than that of the cytoplasmic $Ca^{2+}$ concentration⁵². As a consequence, the lumen $Ca^{2+}$ concentration is taken as a homogeneous and infinitely large reservoir of $Ca^{2+}$ (term in equation (27)).

## Enzymatic activity of nucleotidases

In addition to transmitting chemical or mechanical cues from the blood lumen into the cytoplasm of cell, ECs play a vital role in maintaining blood ATP concentration within the physiological range. As discussed in Subsection Boundary conditions, ectonucleotidase enzymes act as a safeguard against unwanted or abnormal activation of pathways or reactions with receptors due to uncontrolled ATP concentration in the vascular lumen²³. These enzymes hydrolyze ATP into ADP and eventually, ADP into AMP and AMP into AD. As a result, the potential harmful effects of ATP are avoided. A recent study²⁵ did not taken into account the role of ECs, so that ATP released by RBCs increased with time without bound. We include here the role of ECs in our model as a prescribed boundary condition (i.e., equation (12)). To understand the role of ECs in lumen ATP concentration, we have analyzed the behavior of the spatially-averaged ATP concentration with time. Interestingly, we found that ATP concentration reaches a steady-state in the lumen (Fig. 4). This result is exactly due to the balance of ATP release fluxes from RBCs and ECs with hydrolysis of ATP by the ectonucleotidase enzymes. We also demonstrate the behavior of ATP in the absence of endothelium intervention in that figure. Here, we have started the simulation with an initial ATP concentration higher than the steady ATP concentration.

## Wall ATP concentrations

Prior to coupling the response from the endothelium, i.e., studying $Ca^{2+}$ signaling, we first analyze how ATP concentrations vary along the top and bottom endothelium walls for different degrees of confinement, *Cn*, in a straight channel. There are almost no variation in ATP concentrations along both walls of the endothelium (Figs. 5(a)), as expected. The ATP concentration difference between the two walls is quite small. This is due to fact that after a long time has elapsed in the simulation, the RBCs in the channel occupy equal distances from the walls. As a consequence, the RBCs near to the two walls experience equal shear stress and ATP release. This leads practically to the same ATP concentration for both walls, as shown in Fig. 5(a). For this reason, in next sections we will only consider the maximum ATP concentration from the bottom wall, which will be used as an input for $Ca^{2+}$ signaling activation in ECs, for various *Cn*, *Ca*, and *Hct*.

To understand how the ATP concentration varies across the channel for different *Cn* with a given *Hct* and *Ca*, we represent the transverse time- and space-averaged ATP concentration, as shown in Fig. 5(b). We find that the ATP concentration is higher near the walls than that in the core region. This is owing to the facts: (1) ATP concentration is 0 in the area occupied by RBCs. Therefore, the overall ATP concentration in the core region is lower than that near the walls. (2) RBCs spread closer to the walls due to cell-cell hydrodynamic interactions. The high shear stress near the walls induces and increase of ATP release, resulting in a higher ATP concentration. This high concentration of ATP close to the walls is beneficial for calcium signaling in ECs, which we will discuss in upcoming sections. With varying *Cn*, we observe that the ATP concentration is higher in wider channels than in smaller width channels, as shown in Fig. 5(b). This can be explained by the following argument: the number of RBCs is inversely proportional to *Cn* for a constant *Hct*. Thus, there are more RBCs in the case where *Cn* = 0.2 than in the cases *Cn* = 0.4 and 0.8, leading to a greater accumulated ATP release by RBCs and consequently a higher ATP concentration. The ATP concentrations both near the walls and core regions are higher for the case *Cn* = 0.2 compared to the case *Cn* = 0.4 and 0.8. Similarly, the ATP concentrations both near the walls and core regions are higher in *Cn* = 0.4 than in *Cn* = 0.8; however, the differences are minimal.

## Effect of capillary number

We have carried out several simulations by choosing different *Ca* while fixing confinement (*Cn*) to 0.2. In real microvascular networks, RBCs experience different wall shear rates when traversing through the vessels of the same diameter. To emulate such a fact in the microvasculature, we have performed simulations by varying *Ca*. In numerical simulations, we are, in principle, in the Stokes limit (zero *Re*), but the lattice Boltzmann method (LBM) introduces inherently a finite Reynolds number, that must be kept small enough. The real *Re* is about $10^{-4}$, and a choice of this value in the LBM would make computation quite prohibitive. There is a consensus in the literature that fixing *Re* = 0.1 helps accelerate the simulation speed without compromising the accuracy of simulation results⁴⁷. Since we aim to keep the same *Re* for all simulations, we have fixed the
<figure><img src="image_4.png" /><figcaption>**Figure 4.** The evolution of bulk ATP concentration with time in the presence and absence of endothelium (*Ca* = 22.5, *Cn* = 0.4, and *Hct* = 16.7%).</figcaption></figure>

flow speed (to have $Re = 0.1$) to the same value for all simulations, while varying the bending rigidity to explore different values of $Ca$. Fig. 6(a) shows the maximum time-averaged bottom wall ATP concentration. The concentration of ATP at the wall surface increases with an increase in $Hct$ concentration for each $Ca$. From Fig. 6(a), the following conclusions can be drawn: (1) the maximum ATP release level (ATP releasing capacity of a RBC to its maximum ATP releasing capacity), as reported by²⁵, does not play a vital role in deciding the quasi-steady ATP concentration, (2) the final ATP concentration depends on $Hct$ i.e., high $Hct$ means the higher cumulative ATP concentration and total ATP release level (3) $Ca = 9$ is unable to cause any significant change in the concentration of ATP for different $Hct$. This is due to fact that the critical shear stress required to activate $Px_1$ hemichannels on the RBC membrane is not attained. Consequently, the concentration of ATP does not change significantly, (4) furthermore upon increasing of $Ca$ (i.e., $Ca = 22.5$ and 45), more $Px_1$ hemichannels are opened due to the high shear stress experienced by RBC membrane. This membrane shear stress is higher than the critical shear stress required to activate these hemichannels. Further, the ATP released from ECs does not play a significant role to lumen ATP concentration for $Ca = 22.5$ and 45. And also, the contribution of CFTR on the up-regulation of $Px_1$ hemichannels is also negligible due to the small deformation of the RBC membrane, (5) with further increase in the flow strength i.e., $Ca$ ($Ca = 90$), almost all $Px_1$ hemichannels channels are activated as reported by Gou et al.²⁵. In the meantime, these channels are also up-regulated by CFTRs, which are activated due to the deformation of the membrane at $Ca = 90$. Furthermore, the shear stress experienced by the EC is sufficient to contribute to the lumen ATP concentration. Another important remark is that for high $Hct$ concentration and $Ca = 90$ more RBCs are pushed towards the wall due to cell-cell interaction, resulting in higher shear stress experienced by RBCs and higher ATP release, and (6) the ATP concentration at the wall is linear for the dilute concentrations (i.e., up to 10 %). After that regime, the wall ATP concentration becomes nonlinear with increasing in $Hct$.
<figure><img src="image_5.png" /><figcaption>**Figure 5.** ATP concentrations are reported for *Ca* = 90. (a) Time-averaged wall ATP concentrations for *Cn* = 0.8 and *Hct* = 33.4%; and (b) the time- and space-averaged ATP concentration across the width of the channel for various *Cn* and *Hct*.</figcaption></figure>

Finally, we provide ATP wall concentration as an input to solve Ca²⁺ signaling equations (19)-(23). We present the maximum cytoplasmic Ca²⁺ concentration (Ca²⁺ peak) in Fig. 6(b). The response of an EC, the cytoplasmic Ca²⁺ concentration, follows the same trend as the wall ATP concentration (see Fig. 6(a)). The cytoplasmic Ca²⁺ concentration rises above the homeostatic concentration, 0.1 µM for all *Ca* and *Hct*. However, the response for *Ca* = 9 is very negligible in comparison to other *Ca*. This is due to the fact that the low concentration of ATP on the surface of the endothelium wall is unable to produce IP₃ molecules in sufficient quantities by phosphorylating free receptors. As a consequence, the concentration of IP₃ molecules required to activate IP₃R (inositol trisphosphate receptor) channels present on the membrane of ER is not sufficient to mobilize Ca²⁺.

Effect of confinement

In this section, we study the effect of *Cn* on ATP and the cytoplasmic Ca²⁺ concentration for *Ca* = 45. The results obtained with increasing *Cn* are similar to the observations made in previous subsection. In addition, the following observations can be drawn from Fig. 7(a): (1) interestingly, the ATP concentrations for *Cn* = 0.4 and *Cn* = 0.8 are very close to each other. There are no significant differences in the results within the dilute regime (upto *Hct* = 10%). However, a difference is visible beyond the dilute regime. Gou et al.²⁵ have reported that the ATP release level per RBC for *Cn* = 0.8 is always higher than the ATP release level for *Cn* = 0.4 for all *Hct* concentrations. This is due to the fact that at *Cn* = 0.8, RBCs flow close to the wall (see Fig. 7(b)) like a train for different *Hct* concentrations. As result, the RBCs experience higher membrane shear stress, which helps increase the ATP release level. It has already been proven that the shear stress-dependent ATP release through *Px*₁ hemichannels predominates over the CFTR amplification factor. The CFTR amplification factor is associated with RBC deformation. This conclusion was drawn based on measurements of the averaged shear stress per RBC and the averaged rate of curvature change per RBC for different flow strengths. Specifically, it was observed that, particularly at high enough *Ca* levels, the average shear stress exceeds the critical shear stress required to activate *Px*₁ hemichannels. However, the average rate of local curvature change remains below the critical rate of curvature change in the majority of configurations²⁴,²⁵ (except when a cell hits a bifurcation, a geometry which is outside the present scope). Consequently, shear stress-dependent ATP release assumes greater significance than RBC deformation. However, for the same length of channel the number of RBCs in *Cn* = 0.4 is greater than in *Cn* = 0.8, resulting in a higher total ATP release level in *Cn* = 0.4 than the total ATP release level in *Cn* = 0.8; and (2) the rapid increase in ATP concentration is observed for *Cn* = 0.2. The reason is due to the development of the complex cell-cell interactions and dispersion of RBCs in the channel. As a result, RBCs experience higher shear stress and membrane deformation which results in greater ATP concentration in contrast to other values of *Cn* = 0.8 and 0.4. Fig. 7(c) shows that at *Cn* = 0.4, the RBCs are oriented in an order and close to the center of the channel. In contrast, at *Cn* = 0.2 the RBCs are randomly oriented and distributed in the channel, as shown in Fig. 7(d).

There is not a significant visible ATP gradient in the plasma for *Cn* = 0.8 and 0.4, as shown in Figs. 7(b)-(c). However, a visible ATP gradient does exist for *Cn* = 0.2, as shown in Fig. 7(d). The concentration of ATP between the spaces between
<figure><img src="image_6.png" /><figcaption>**Figure 6.** These results are obtained for $C_n = 0.2$ on the bottom wall of a 2-D channel: (a) maximum ATP concentrations on the bottom wall with varying *Hct* concentrations for various *Ca*; and (b) maximum cytoplasmic $Ca^{2+}$ concentrations (transient peaks) obtained from the maximum ATP concentrations for different *Ca*.</figcaption></figure>

RBCs is higher than in the cell-free layer (CFL) regions for $C_n = 0.2$. Due to cell-cell interactions, the released ATP from the core region does not diffuse freely into the CFL regions, where RBCs serve as obstacles. Consequently, the concentration of ATP in the free spaces between cells remains higher than in the CFL regions. At higher RBC concentrations, for $C_n = 0.4$ and 0.8, RBCs flow in double and single-file trains, respectively (shown in Figs. 7(b)-(c)). In both cases, the core regions are crowded, with minimal free spaces between cells, and the CFL layers are uniform. The ATP release is higher near the walls in these confinements, as cells are very close to the walls where shear stress is highest. The released ATP near the walls accumulates in the CFL regions and is balanced by the hydrolysis of enzymes present on the endothelial wall surface (shown in Fig. 5(b)), resulting in a stable ATP concentration. Therefore, due to the small CFL thickness and velocity in the wall vicinity, the ATP in the CFL regions can be assumed to be diffusion dominated, at least for $C_n = 0.4$ and 0.8.

The maximum cytoplasmic $Ca^{2+}$ concentration observed in an EC follows the same trend as ATP concentration (Fig. 7(e)). The maximum cytoplasmic $Ca^{2+}$ increases monotonously with increase in *Cn* and *Hct*. Moreover, an increase in the maximum cytoplasmic $Ca^{2+}$ concentration may lead to propagation of $Ca^{2+}$ signals from region of high cytoplasmic $Ca^{2+}$ concentration to lower cytoplasmic $Ca^{2+}$ concentrations. A real vascular network has different vessel sizes as well as different local hematocrits. According to the results of Fig. 7(e), it may happen that ATP concentration reaches the threshold value for calcium signaling in some locations (e.g. with high enough *Hct*, or large shear stresses). This naturally raises the question of calcium propagation in vascular network. It is hoped to investigate this matter in a future work. Here, we will first address the question of the transient peak time of the cytoplasmic $Ca^{2+}$ concentration, which will help clarify the coordinating cellular function for calcium propagation in a vascular network.

## Cytoplasmic $Ca^{2+}$ peak time

In order to understand the possibility of intercellular communication in the vascular wall, we present the transient peak time of the cytoplasmic $Ca^{2+}$. The response of an EC, the maximum cytoplasmic $Ca^{2+}$ concentration, is dose-dependent i.e., high ATP concentration leads to high cytoplasmic $Ca^{2+}$ concentration. The dependence between ATP and cytoplasmic $Ca^{2+}$ is sigmoidal²⁸ due to the fixed quantity of production of $IP_3$ molecules. Furthermore, the high ATP concentration leads to a decrease in the cytoplasmic $Ca^{2+}$ transient time. In contrast, the low ATP concentration leads to a slower response from an EC. This is elucidated in Fig. 8. Fig. 8(a) shows the dependence of the calcium peak time on *Hct* and *Ca*. The cytoplasmic transient time for a given *Ca* decreases with an increase in *Hct* concentration (see Fig. 8(a)). Similarly, the transient peak time decreases monotonically with an increase in *Ca*. The peak time difference between two different values of *Ca* and a given *Hct* concentration varies between 2-5 s. As a result calcium can propagate from one site to another, both downstream and upstream. This helps coordinating cellular functions, such as cellular proliferation, wound healing, vasodilation, etc. Since calcium propagation speed (in the range of $10\mu m/s^{53-55}$) is lower than the blood flow speed (in the range of mm/s in microcirculation), the vasodilation is controlled upstream only.

Fig. 8(b) shows that for a given *Hct* concentration at $C_n = 0.4$ and 0.8, there is no significant difference in the cytoplasmic
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** (a) Maximum ATP concentration on the bottom wall by varying *Hct* concentration for various *Cn*; (b) spatial concentration of ATP at *Cn* = 0.8, *Ca* = 45, *Hct* = 33.4%, and *t* = 18.56 s; (c) spatial concentration of ATP at *Cn* = 0.4, *Ca* = 45, *Hct* = 43.4%, and *t* = 26.9 s; (d) spatial concentration of ATP at *Cn* = 0.2, *Ca* = 45, *Hct* = 45.9%, and *t* = 25.87 s; and (e) the maximum cytoplasmic Ca²⁺ concentration (transient peak) obtained from the corresponding maximum ATP concentration for various *Cn*. All simulations are carried out for *Ca* = 45. The Color bars are removed for *Cn* = 0.8 and 0.4 due to the absence of visible ATP concentration gradients in plasma.</figcaption></figure>

Ca²⁺ peak time. The cytoplasmic Ca²⁺ peak time at *Cn* = 0.2 is much faster than the transient peak time at *Cn* = 0.4 and 0.8. This is due to the fact that there exists higher ATP concentration at *Cn* = 0.2 (see Fig. 7(a)). Therefore, there is a possibility of signal transfer from large width channel (weakly confined channel) to small width channels (strongly confined channels). In experiments with the vascular network stimulating a region with an agonist concentration⁵⁶⁻⁵⁸, they observed propagation of Ca²⁺ to the downstream smaller vessel diameters. Our findings would be useful for understanding Ca²⁺ propagation both in-vitro and in-vivo vessel networks.

Calcium wave propagation

Intercellular Ca²⁺ waves are responsible for coordinated functions of the cell activities such as proliferation or division of cells and transmission of vasodilation signals to upstream and downstream in blood vessel⁵⁴. To demonstrate the propagation of a Ca²⁺ wave from an EC stimulated by a high ATP concentration to its neighboring ECs, we perform a simple study by considering a group of horizontally aligned ECs. These ECs communicate with each other via gap junctions, which allow Ca²⁺ and IP₃ to transfer to their neighboring ECs⁵⁵. This is one of several mechanisms through which Ca²⁺ and IP₃ diffuse in vascular wall. However, there are other mechanisms through which Ca²⁺ propagation occurs in vascular wall. The propagation of action potential through ECs and SMCs from the stimulated site (high ATP concentration region)⁵⁹ to the neighboring ECs and SMCs. The release of ATP from an EC into the extracellular space due to the mechanical stimulation and then diffusion
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The cytoplasmic Ca²⁺ peak in an EC: (a) the peak time obtained by varying *Hct* concentration for different *Ca* and *Cn* = 0.2; and (b) the peak time obtained by varying *Hct* for different *Cn* and *Ca* = 45.</figcaption></figure>

of ATP results in the propagation of Ca²⁺ waves over several distances in vascular wall⁶⁰,⁶¹. To investigate intercellular

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Schematic shows that a group of ECs are horizontally aligned and the middle EC is stimulated to a ATP concentration. Following stimulation, IP₃ and cytoplasmic Ca²⁺ diffuse through gap junctions into the neighboring ECs. IP₃ molecules then release sequestered Ca²⁺ from neighboring ECs. The propagation of Ca²⁺ ions continues until IP₃ molecules are able to activate IP₃R (inositol trisphosphate receptor) present on ER membrane. The width of the arrow indicates the strength of propagation of IP₃ molecules and cytoplasmic Ca²⁺ ions into the neighboring ECs.</figcaption></figure>

Ca²⁺ waves in the current model, we begin with the same model equations given in Ca²⁺ signaling subsection, and we additionally include the spatial variation of IP₃ and the cytoplasmic Ca²⁺ concentrations by incorporating the diffusion of IP₃ and cytoplasmic Ca²⁺ in equations (21)-(22). So, the new modified equations (equations (21)-(22)) become

$$ \frac{\partial [IP_3]}{\partial t} = D_{IP_3} \frac{\partial^2 [IP_3]}{\partial x^2} + q_{\text{gen}} - q_{\text{deg}}, \qquad (34) $$

$$ \frac{\partial [Ca_c]}{\partial t} = D_{ca} \frac{\partial^2 [Ca_c]}{\partial x^2} + \beta (q_{\text{rel}} - q_{\text{serca}} + q_{\text{leak}} + q_{\text{soc}} + q_{\text{in}} - q_{\text{pmca}}), \qquad (35) $$

where $D_{IP_3}$ is the diffusion coefficient of IP₃ (i.e., 280 µm²/s⁵⁵), and $D_{ca}$ is the effective cytoplasmic diffusion coefficient of Ca²⁺ (i.e., 20 µm²/s⁵⁵). The boundary conditions at the EC-EC interface are

$$ D_{IP_3} \frac{\partial [IP_3]}{\partial x} = P_{IP_3} ([IP_3]_+ - [IP_3]_-), \qquad (36) $$
$$D_{ca} \frac{\partial [Ca_c]}{\partial x} = P_{ca} ([Ca_c]_+ - [Ca_c]_-), \qquad (37)$$

where $P_{IP_3}$ and $P_{ca}$ are the permeability rate constants for $IP_3$ and cytoplasmic $Ca^{2+}$, which are estimated to 0.05 µm/s and 5.0 µm/s, respectively$^{55}$. $[IP_3]_+$ and $[IP_3]_-$ are $IP_3$ concentrations on either sides of the boundary, and $[Ca_c]_+$ and $[Ca_c]_-$ are cytoplasmic $Ca^{2+}$ concentrations on either sides of the boundary. The concentration variation in the lateral direction is smaller

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Evolution of cytoplasmic $Ca^{2+}$ concentration for the stimulated EC and also ECs on the right hand side of that EC. The stimulated EC is treated with 1 µM concentration of ATP.</figcaption></figure>

than that in the longitudinal direction, implying that molecular diffusion occurs more rapidly in the lateral direction. For the purpose of demonstration, we limit our investigation to a cluster of 9 endothelial cells (ECs). Specifically, there are 4 ECs on each side of the stimulated EC. For instance, we locally stimulate the middle EC (indicated by the downward arrow in Fig. 9) with a concentration of 1 µM ATP. This localized stimulation triggers signal propagation on both sides of the stimulated EC. The stimulated EC is designated as “EC0”, and the ECs to its right are labeled as “EC1”, “EC2”, “EC3”, and “EC4”. Similarly, the ECs to its left are named “EC5”, “EC6”, “EC7”, and “EC8”. We maintain the constants as provided in$^{28}$, with the exception of altering the values of $K_r$ and $k_2$ to 0.5 µM and 0.02 s$^{-1}$ respectively, to assess their impact on the speed of $Ca^{2+}$ wave propagation.

Because of the symmetry of the problem, we present the evolution of cytoplasmic $Ca^{2+}$ (shown in Fig. 10) for ECs on the right hand of the stimulated EC (“EC0”) such as “EC1”, “EC2”, “EC3”, and “EC4”. As we move away from “EC0”, we find variations in cytoplasmic $Ca^{2+}$ peak magnitude as well as in peak time in the neighboring ECs. The increase in $Ca^{2+}$ concentration in ECs neighboring to “EC0” is due to the fact that $IP_3$ molecules diffuse from “EC0” to “EC1” and eventually to
other ECs. As a consequence, the action of these molecules on IP₃R channels results in release of Ca²⁺ from the ER. In Fig. 10, we find that the Ca²⁺ peak in “EC1” occurs approximately 2 s after that in “EC0”. Therefore, in the average sense, the Ca²⁺ signal moves at a speed of 5 µm/s. As we move further away from “EC1”, the speed of the Ca²⁺ signal decays. This is due to the fact that the concentration of IP₃ decreases as it moves away from “EC0” because of passive diffusion of IP₃ molecules between ECs (equation (36)) and the degradation of its concentration in ECs (third term in equation (34)). Consequently, a decline in the cytoplasmic Ca²⁺ peak concentration is observed (shown in Fig. 10).

In experiments⁵³⁻⁵⁵,⁶², it has been reported that the Ca²⁺ wave can propagate at a speed of 10 µm/s or more, which compares well with the above order of magnitude. However, an important remark is in order. We have not observed a Ca²⁺ wave front movement in our simulations, strictly speaking, in contrast to experiments⁵³⁻⁵⁵,⁶². It is likely that for wave front propagation to occur, an active generation of Ca²⁺ or IP₃ would be required in the neighboring cells while Ca²⁺ or IP₃ spreads from the stimulated cell through cell-cell junction⁵⁴,⁶²,⁶³. This effect will be thoroughly investigated in the future in order to have a complete picture of wave propagation.

## Conclusions

In this study, we have successfully integrated the dynamic behavior of RBCs within a parabolic flow with the intricate dynamics of calcium ions (Ca²⁺). While the primary function of RBCs is to facilitate gas exchange between lungs and tissues, they also play a crucial secondary role. RBCs release ATP as a signaling molecule in response to externally imposed flow. These ATP molecules act as messengers, transmitting signals from the external environment to the ECs, initiating the release of ubiquitous Ca²⁺ ions. The presence of Ca²⁺ ions plays a pivotal role in regulating numerous enzymatic functions. Notably, Ca²⁺ governs the production of NO, a vasodilator crucial for maintaining vascular tone and overall cellular function. Consequently, this intricate system occurring in the vascular wall plays a vital role in both cellular regulation and vascular homeostasis.

To delve into the intricate relationship between RBC dynamics and biochemical signaling, we conducted a series of numerical simulations. By systematically varying parameters such as channel width (*Cn*), flow strength (*Ca*), and RBC concentrations (*Hct*) in a straight channel, we sought to uncover the underlying mechanisms. Employing the immersed boundary lattice Boltzmann method (IBLM), we solved fluid-RBC interactions and ATP release mechanisms. Our investigation addressed a range of unanswered questions, including the impact of maximum ATP release levels on wall ATP concentration, the variation of ATP concentrations on vascular walls, the interplay between flow strength and confinement, and the potential existence of Ca²⁺ wave propagation within a vascular network.

The main outcomes of the current study are summarized as follows: (1) Ectonucleotidase enzymes on the blood vessel walls play a critical role in preventing an abnormal rise in plasma ATP concentration and the activation of biochemical pathways in ECs. (2) There is no variation in steady ATP concentration along the flow direction, as well as minimal difference in its concentration between the bottom and top walls. (3) For a given channel length, hematocrit (*Hct*), and capillary number (*Ca*), the number of RBCs is higher in wider channels than in smaller channels with changes in *Cn*. Consequently, wider channels have a higher ATP concentration due to the increased total ATP release level by RBCs. (4) With an increase in *Ca* for a given *Cn*, changes in *Hct* lead to a monotonic increase in ATP concentration, resulting in corresponding monotonic rises in cytoplasmic Ca²⁺ peak concentrations. This response is linear in dilute RBC concentrations; however, it becomes nonlinear in concentrated RBC concentrations due to the increase in cell-cell hydrodynamic interactions pushing RBCs towards the channel walls. This leads to an increase in ATP release due to higher membrane stress experienced by cells near the walls. (5) Similarly, keeping *Ca* constant, an increase in *Hct* leads to an increase in ATP concentration for a given *Cn*. The ATP concentration in wider channels is higher than that in smaller channels. This is due to the fact that a greater number of RBCs are flowing in wider channels compared to smaller channels, resulting in a higher total ATP release level from RBCs. The response of cytoplasmic Ca²⁺ peak concentrations follows a similar trend to the ATP response, with larger ATP concentrations leading to higher cytoplasmic Ca²⁺ peak concentrations. (6) The cytoplasmic Ca²⁺ peak time is shorter with higher ATP concentration, indicating a faster response. The difference in relative peak times between different magnitudes of ATP concentrations at various regions of the channel leads to the propagation of Ca²⁺ waves at a speed of 5 µm/s upstream and downstream of the channel. This is demonstrated by locally stimulating an EC within a confluent layer of ECs with ATP concentration, resulting in the release of Ca²⁺ in neighboring ECs. This is mainly attributed to the diffusion of IP₃ molecules from a stimulated cell to neighboring cells, causing the release of sequestered Ca²⁺ from the ER. The passive propagation of the Ca²⁺ wave continues until IP₃ molecules successfully activate IP₃R channels present on the ER membrane. The propagation of the Ca²⁺ wave upstream from a region in a vascular network helps coordinate cellular functions and control the flux of RBCs to fulfill the metabolic needs required for cellular functions.

While our study was conducted within a straight channel, it establishes a solid foundation for understanding RBC dynamics in the context of biochemical signaling under different (patho)physiological conditions. Our work can be extended to more complex vascular networks, where RBCs experience varying shear rates, local variations in RBC concentrations, ATP concentrations across different branches, and the propagation of Ca²⁺ signals throughout upstream and downstream vessel
branches. These expansions would significantly contribute to our understanding of localized blood flow regulation in real vasculature. In current study, we limited ourselves to 2D, which, at least, should provide qualitative results. this constitute an essential step in this broad area of research which is still in its nascent stage. This is also motivated by the limited availability of quantitative experimental data. For example, experiments regarding Ca²⁺ signaling of endothelial cells (ECs) in the presence of blood flow are largely missing. The present study can be viewed as guide for future experiments before dealing with a full advanced 3D model.

Moreover, beyond ATP release, RBCs are also recognized as potential scavengers of NO, a vasodilator⁶⁴. Therefore, a comprehensive investigation is required to unveil the intricate interplay between factors such as the cell-free layer (CFL), ATP concentration, RBC concentration, Ca²⁺ concentration, and their combined effects on the bioavailability of NO within the vascular wall. Such endeavors would not only advance our scientific knowledge but also provide insights into the multifaceted mechanisms governing vascular health and function. Furthermore, incorporating membrane viscosity into the RBC membrane model influences the dynamics of RBCs, including deformation, suspension viscosity, and particle stress⁴⁸,⁴⁹. It would be interesting to explore the effect of membrane viscosity on ATP release from RBCs under various conditions and to quantitatively validate it with experimental results in future studies²

## References

1. Amiri, F. A. & Zhang, J. Oxygen transport across tank-treading red blood cell: individual and joint roles of flow convection and oxygen-hemoglobin reaction. *Microvasc. Res.* **145**, 104447 (2023).

2. Forsyth, A. M., Wan, J., Owrutsky, P. D., Abkarian, M. & Stone, H. A. Multiscale approach to link red blood cell dynamics, shear viscosity, and ATP release. *Proc. Natl. Acad. Sci.* **108**, 10986–10991 (2011).

3. Lohman, A. W., Billaud, M. & Isakson, B. E. Mechanisms of ATP release and signalling in the blood vessel wall. *Cardiovasc. research* **95**, 269–280 (2012).

4. Xu, S. et al. New insights into shear stress-induced endothelial signalling and barrier function: cell-free fluid versus blood flow. *Cardiovasc. research* **113**, 508–518 (2017).

5. Pirotton, S., Raspé, E., Demolle, D., Erneux, C. & Boeynaems, J.-M. Involvement of inositol 1, 4, 5-trisphosphate and calcium in the action of adenine nucleotides on aortic endothelial cells. *J. Biol. Chem.* **262**, 17461–17466 (1987).

6. Song, Z., Vijayaraghavan, S. & Sladek, C. D. Atp increases intracellular calcium in supraoptic neurons by activation of both P2X and P2Y purinergic receptors. *Am. J. Physiol. Integr. Comp. Physiol.* **292**, R423–R431 (2007).

7. Carter, T., Newton, J., Jacob, R. & Pearson, J. Homologous desensitization of ATP-mediated elevations in cytoplasmic calcium and prostacyclin release in human endothelial cells does not involve protein kinase c. *Biochem. journal* **272**, 217–221 (1990).

8. Zhao, Y., Vanhoutte, P. M. & Leung, S. W. Vascular nitric oxide: Beyond eNOS. *J. pharmacological sciences* **129**, 83–94 (2015).

9. Sriram, K., Laughlin, J. G., Rangamani, P. & Tartakovsky, D. M. Shear-induced nitric oxide production by endothelial cells. *Biophys. journal* **111**, 208–221 (2016).

10. Adams, D. J., Barakeh, J., Laskey, R. & Van Breemen, C. Ion channels and regulation of intracellular calcium in vascular endothelial cells. *The FASEB J.* **3**, 2389–2400 (1989).

11. Hutcheson, I. R. & Griffith, T. M. Central role of intracellular calcium stores in acute flow-and agonist-evoked endothelial nitric oxide release. *Br. journal pharmacology* **122**, 117–125 (1997).

12. Kuchan, M. & Frangos, J. Role of calcium and calmodulin in flow-induced nitric oxide production in endothelial cells. *Am. J. Physiol. Physiol.* **266**, C628–C636 (1994).

13. Berridge, M. J., Bootman, M. D. & Roderick, H. L. Calcium signalling: dynamics, homeostasis and remodelling. *Nat. reviews Mol. cell biology* **4**, 517–529 (2003).

14. Bartlett, P. J., Cloete, I., Sneyd, J. & Thomas, A. P. IP₃-dependent Ca²⁺ oscillations switch into a dual oscillator mechanism in the presence of PLC-linked hormones. *IsScience* **23**, 101062 (2020).

15. Blatter, L. A. et al. Simultaneous measurements of Ca²⁺ and nitric oxide in bradykinin-stimulated vascular endothelial cells. *Circ. Res.* **76**, 922–924 (1995).

16. Marrelli, S. P. Mechanisms of endothelial P2Y₁- and P2Y₂-mediated vasodilatation involve differential [Ca²⁺]ᵢ responses. *Am. J. Physiol. Circ. Physiol.* **281**, H1759–H1766 (2001).

17. Davignon, J. & Ganz, P. Role of endothelial dysfunction in atherosclerosis. *Circulation* **109**, III–27 (2004).
18. Kolluru, G. K., Siamwala, J. H. & Chatterjee, S. eNOS phosphorylation in health and disease. *Biochimie* **92**, 1186–1198 (2010).

19. Carter, T. D. & Ogden, D. Kinetics of intracellular calcium release by inositol 1, 4, 5-trisphosphate and extracellular ATP in porcine cultured aortic endothelial cells. *Proc. Royal Soc. London. Ser. B: Biol. Sci.* **250**, 235–241 (1992).

20. Chachisvilis, M., Zhang, Y.-L. & Frangos, J. A. G protein-coupled receptors sense fluid shear stress in endothelial cells. *Proc. Natl. Acad. Sci.* **103**, 15463–15468 (2006).

21. Xu, Z., Dou, W., Wang, C. & Sun, Y. Stiffness and atp recovery of stored red blood cells in serum. *Microsystems & Nanoeng.* **5**, 51 (2019).

22. Gorman, M. W., Feigl, E. O. & Buffington, C. W. Human plasma ATP concentration. *Clin. chemistry* **53**, 318–325 (2007).

23. Burnstock, G. & Ralevic, V. Purinergic signaling and blood vessels in health and disease. *Pharmacol. reviews* **66**, 102–192 (2014).

24. Zhang, H., Shen, Z., Hogan, B., Barakat, A. I. & Misbah, C. ATP release by red blood cells under flow: model and simulations. *Biophys. journal* **115**, 2218–2229 (2018).

25. Gou, Z., Zhang, H., Abbasi, M. & Misbah, C. Red blood cells under flow show maximal ATP release for specific hematocrit. *Biophys. J.* **120**, 4819–4831 (2021).

26. Plank, M. J., Wall, D. J. & David, T. Atherosclerosis and calcium signalling in endothelial cells. *Prog. biophysics molecular biology* **91**, 287–313 (2006).

27. Plank, M., Wall, D. & David, T. The role of endothelial calcium and nitric oxide in the localisation of atherosclerosis. *Math. biosciences* **207**, 26–39 (2007).

28. Nayak, A. K., Gou, Z., Das, S. L., Barakat, A. I. & Misbah, C. Mathematical modeling of intracellular calcium in presence of receptor: a homeostatic model for endothelial cell. *Biomech. Model. Mechanobiol.* **22**, 217–232 (2023).

29. Kaoui, B. et al. Complexity of vesicle microcirculation. *Phys. Rev. E* **84**, 041906 (2011).

30. Kaoui, B., Ristow, G., Cantat, I., Misbah, C. & Zimmermann, W. Lateral migration of a two-dimensional vesicle in unbounded Poiseuille flow. *Phys. Rev. E* **77**, 021903 (2008).

31. Coupier, G., Kaoui, B., Podgorski, T. & Misbah, C. Noninertial lateral migration of vesicles in bounded Poiseuille flow. *Phys. Fluids* **20**, 111702 (2008).

32. Beaucourt, J., Rioual, F., Séon, T., Biben, T. & Misbah, C. Steady to unsteady dynamics of a vesicle in a flow. *Phys. Rev. E* **69**, 011906 (2004).

33. Kantsler, V. & Steinberg, V. Transition to tumbling and two regimes of tumbling motion of a vesicle in shear flow. *Phys. review letters* **96**, 036001 (2006).

34. Abbasi, M. et al. Dynamics and rheology of a single two-dimensional multilobe vesicle in a confined geometry. *Phys. Rev. Fluids* **7**, 093603 (2022).

35. Peskin, C. S. The immersed boundary method. *Acta numerica* **11**, 479–517 (2002).

36. Helfrich, W. Elastic properties of lipid bilayers: theory and possible experiments. *Zeitschrift für Naturforschung C* **28**, 693–703 (1973).

37. Scheffer, L., Bitler, A., Ben-Jacob, E. & Korenstein, R. Atomic force pulling: probing the local elasticity of the cell membrane. *Eur. Biophys. J.* **30**, 83–90 (2001).

38. Shen, Z., Farutin, A., Thiébaud, M. & Misbah, C. Interaction and rheology of vesicle suspensions in confined shear flow. *Phys. Rev. Fluids* **2**, 103101 (2017).

39. Burnstock, G. Release of vasoactive substances from endothelial cells by shear stress and purinergic mechanosensory transduction. *The J. Anat.* **194**, 335–342 (1999).

40. Ando, J., Komatsuda, T. & Kamiya, A. Cytoplasmic calcium response to fluid shear stress in cultured vascular endothelial cells. *In vitro cellular & developmental biology* **24**, 871–877 (1988).

41. Shen, J., Luscinskas, F. W., Connolly, A., Dewey Jr, C. F. & Gimbrone Jr, M. Fluid shear stress modulates cytosolic free calcium in vascular endothelial cells. *Am. J. Physiol. Physiol.* **262**, C384–C390 (1992).

42. Geiger, R. V., Berk, B. C., Alexander, R. W. & Nerem, R. M. Flow-induced calcium transients in single endothelial cells: spatial and temporal analysis. *Am. J. Physiol. Physiol.* **262**, C1411–C1417 (1992).
43. Yamamoto, K. et al. Endogenously released ATP mediates shear stress-induced Ca²⁺ influx into pulmonary artery endothelial cells. *Am. J. Physiol. Circ. Physiol.* **285**, H793–H803 (2003).

44. Ando, J. & Yamamoto, K. Flow detection and calcium signalling in vascular endothelial cells. *Cardiovasc. research* **99**, 260–268 (2013).

45. John, K. & Barakat, A. I. Modulation of ATP/ADP concentration at the endothelial surface by shear stress: effect of flow-induced ATP release. *Annals biomedical engineering* **29**, 740–751 (2001).

46. Zhang, H. & Misbah, C. Lattice boltzmann simulation of advection-diffusion of chemicals and applications to blood flow. *Comput. & Fluids* **187**, 46–59 (2019).

47. Kaoui, B. & Harting, J. Two-dimensional lattice boltzmann simulations of vesicles with viscosity contrast. *Rheol. Acta* **55**, 465–475 (2016).

48. Gürbüz, A., Pak, O. S., Taylor, M., Sivaselvan, M. V. & Sachs, F. Effects of membrane viscoelasticity on the red blood cell dynamics in a microcapillary. *Biophys. J.* **122**, 2230–2241 (2023).

49. Guglietta, F., Pelusi, F., Sega, M., Aouane, O. & Harting, J. Suspensions of viscoelastic capsules: effect of membrane viscosity on transient dynamics. *J. Fluid Mech.* **971**, A13 (2023).

50. Prado, G., Farutin, A., Misbah, C. & Bureau, L. Viscoelastic transient of confined red blood cells. *Biophys. journal* **108**, 2126–2136 (2015).

51. Guglietta, F., Behr, M., Biferale, L., Falcucci, G. & Sbragaglia, M. On the effects of membrane viscosity on transient red blood cell dynamics. *Royal Soc. Chem.* **16**, 6191–6205 (2020).

52. Bisello, A. & Friedman, P. A. PTH and PTHrP actions on kidney and bone. *Princ. bone biology* 665–712 (2008).

53. Halidi, N., Boittin, F.-X., Bény, J.-L. & Meister, J.-J. Propagation of fast and slow intercellular Ca²⁺ waves in primary cultured arterial smooth muscle cells. *Cell calcium* **50**, 459–467 (2011).

54. Leybaert, L. & Sanderson, M. J. Intercellular Ca²⁺ waves: mechanisms and function. *Physiol. reviews* **92**, 1359–1392 (2012).

55. Sera, T. et al. Three-dimensional model of intracellular and intercellular Ca²⁺ waves propagation in endothelial cells. *Biochem. biophysical research communications* **505**, 781–786 (2018).

56. Uhrenholt, T. R., Domeier, T. L. & Segal, S. S. Propagation of calcium waves along endothelium of hamster feed arteries. *Am. J. Physiol. Circ. Physiol.* **292**, H1634–H1640 (2007).

57. Seppey, D., Sauser, R., Koenigsberger, M., Bény, J.-L. & Meister, J.-J. Intercellular calcium waves are associated with the propagation of vasomotion along arterial strips. *Am. J. Physiol. Circ. Physiol.* **298**, H488–H496 (2010).

58. Bagher, P., Davis, M. J. & Segal, S. S. Visualizing calcium responses to acetylcholine convection along endothelium of arteriolar networks in Cx40BAC-GCaMP2 transgenic mice. *Am. J. Physiol. Circ. Physiol.* **301**, H794–H802 (2011).

59. Domeier, T. L. & Segal, S. S. Electromechanical and pharmacomechanical signalling pathways for conducted vasodilatation along endothelium of hamster feed arteries. *The J. physiology* **579**, 175–186 (2007).

60. Edwards, J. R. & Gibson, W. G. A model for Ca²⁺ waves in networks of glial cells incorporating both intercellular and extracellular communication pathways. *J. theoretical biology* **263**, 45–58 (2010).

61. Yamamoto, K. et al. Visualization of flow-induced atp release and triggering of Ca²⁺ waves at caveolae in vascular endothelial cells. *J. cell science* **124**, 3477–3483 (2011).

62. Long, J., Junkin, M., Wong, P. K., Hoying, J. & Deymier, P. Calcium wave propagation in networks of endothelial cells: model-based theoretical and experimental study. *PLoS computational biology* **8**, e1002847 (2012).

63. Schuster, S., Marhl, M. & Höfer, T. Modelling of simple and complex calcium oscillations: From single-cell responses to intercellular signalling. *Eur. J. Biochem.* **269**, 1333–1355 (2002).

64. Wei, Y., Mu, L., Tang, Y., Shen, Z. & He, Y. Computational analysis of nitric oxide biotransport in a microvessel influenced by red blood cells. *Microvasc. research* **125**, 103878 (2019).

## Acknowledgements

We acknowledge the financial support from CNES (Centre National d'Etudes Spatiales) and for having access to data, and the French-German University Programme "Living Fluids" (Grant CFDA-Q1-14). The simulations were performed on the Cactus cluster of the CIMENT infrastructure, which is supported by the Rhône-Alpes region (Grant No. CPER07 13 CIRA).
**Author contributions statement**

A.K.N has contributed to the code and performed simulations and analysis of data. C.M. and S.L.D has designed the research topic and plans. C.M has also involved in funding acquisition, resources, and project administration. All the authors have contributed to the article writing and interpretations.

**Data availability**

All outputs studied during current work are included in this published paper.

**Competing interests**

The author(s) declare no competing interests.