On the Ca2+ elevation in vascular endothelial cells due to inositol trisphosphate-sensitive store receptor activation: a data-driven modeling approach

Agonist-induced Ca2+ signalling is essential for the regulation of many vital functions in endothelial cells (ECs). A broad range of stimuli elevate the cytosolic Ca2+ concentration by promoting a pathway mediated by inositol 1,4,5 trisphosphate (IP3) which causes Ca2+ release from intracellular stores. Despite its importance, there are very few studies focusing on the quantification of such dynamics. Here, by using data from isolated ECs, we constructed and tested a minimalistic model that can be used for capturing the main features (averaged over a cell population) of the Ca2+ response to different IP3 stimulation levels. A suitable description of Ca2+-regulatory function of inositol 1,4,5 trisphosphate receptors (IP3Rs) is identified by comparing the different model variants against experimental mean population data. The same approach is used to numerically assess the relevance of cytosolic Ca2+ buffering, as well as Ca2+ store IP3-sensitivity in the overall cell dynamics. The variability in the dynamics’ features observed across the population can be explained (at least in part) through variation of certain model parameters (such as buffering capacity or Ca2+ store sensitivity to IP3). The results, in terms of experimental fitting and validation, support the proposed minimalistic model as a reference framework for the quantification of the EC Ca2+ dynamics induced by IP3Rs activation.


Introduction
Among the myriad of pathways regulating vascular function, Ca 2+ release from intracellular stores plays a fundamental role in the control of blood pressure and flow across the whole circulation [1,2]. The associated Ca 2+ elevation in the endothelial cytosol has indeed manifold effects on tone regulation including the enhancement of nitric oxide synthase (eNOS) activity [3] and potassium (SKCa/IKCa) channels activation which leads to the hyperpolarization of the plasma membrane [4].
Despite its involvement in different agonist stimulations [5,6], the contribution of Ryanodine Receptors (RyRs) located in the sarcoplasmic/endoplasmic reticulum (SR/ER) membrane to EC function is still not clear as their activation through caffeine does not seem to elevate cytosolic Ca 2+ [7,8]. On the other hand, Ca 2+ release from endoplasmic reticulum (ER) due to activation of IP3Rs plays a prominent role in several exogenous ligand-induced pathways. Among these, the Ca 2+ signalling induced by muscarinic cholinergic receptor M3 stimulation has been ascertained across various experimental settings [9,10,11,12]. Activation of muscarinic M3 receptors with acetylcholine (ACh) enhances the production rate of diacylglycerol and IP3, which is then followed by a cytosolic Ca 2+ elevation. The initial increase in cytosolic Ca 2+ is caused by IP3 binding to its receptor on ER which enables a store IP3mediated Ca 2+ release. The activation of a single IP3R on ER initiates the store Ca 2+ release as a localized Ca 2+ blip, which may expand to neighbouring IP3Rs generating a wider signal, defined as Ca 2+ puff. For rising IP3 concentration levels, this signal can augment in frequency and propagate to other IP3R clusters. The increase in Ca 2+ nearby the IP3R favours the ER Ca 2+ release through a positive feedback mechanism, defined as Ca 2+ -induced Ca 2+release (CICR) mechanism. Depending on the intracellular conditions and stimulation level, this signal may further propagate causing a transient global cell elevation in Ca 2+ , denoted as Ca 2+ wave, that can activate different types of membrane channels [13]. Whilst RyRs in the SR generate local Ca 2+ signalling events (denoted as 'sparks') allowing the communication with the juxtaposed plasma membrane, IP3Rs are capable of generating spatially fixed local release events (called 'pulsars') which project through the elastic lamina to the neighbouring smooth muscles membranes [14]. In pressurized mesenteric arteries, spontaneous Ca 2+ events appear to originate from ER IP3Rs even under resting conditions, and these appear to be regulated by surrounding smooth muscle cells via myoendothelial gap junctions [15]. The depletion of Ca 2+ in the intracellular stores triggers Ca 2+ influx through an capacitative entry mechanism (CCE) that serves for sustaining the cytosolic Ca 2+ elevation and ultimately enables ER Ca 2+ refilling [16,17,3]. This store-operated effect involves the activation of transient receptor potential (TRP) channels located on the EC plasma membrane, whose activation under agonist intervention was characterized in different studies [18,19,20,21,22,23,24]. Once the agonist effect vanishes, Ca 2+ extrusion and re-sequestration through Na + -Ca 2+ exchanger (NCX) and Ca 2+ ATPase pumps (PMCA) retake the cytosolic Ca 2+ concentration to its original basal level [25,26].
To operate as an effective sensor, the endothelium is endowed with different cell phenotypes [27], which may exhibit remarkably heterogeneous Ca 2+ responses to agents involving the IP3R pathway such as ACh [28,13]. However, Carter and Ogden [29] showed that the dependency between the maximum IP3Rs Ca 2+ flux and cytosolic IP3 concentration across a porcine aortic EC population can be well captured through an activation function which is regulated by the Ca 2+ concentration. Therefore the EC discrepancy in the response to ACh may not only be explained with different expression levels of in IP3R isoforms, but also due to other cellular components involved in Ca 2+ buffering, extrusion and store sequestration. Furthermore, this study well documented the inhibitory role of Ca 2+ on IP3Rs by considering a broad range (up to saturation level) of IP3 stimulations and showed that IP3Rs inhibition occurs independently of the level of IP3 stimulation. It is also worth noting that, apart from the study by Carter and Ogden [29], there is scarce quantitative information in the literature about the activation role of Ca 2+ on EC IP3Rs at low cytosolic Ca 2+ concentrations.
Depending on the level of detail required, Ca 2+ release via ER IP3R can be modelled by using different approaches including Markov models, gating models and modal models [30]. De Young-Keizer model developed a framework [31,32,33] for describing the kinetics of the IP3R channel release in the ER, which provided the basis for the development of stochastic models for studying CICR and Ca 2+ propagation phenomena (from puffs to waves) within the cytosol [34,35,36,37,38]. On the other hand, gating models provide a minimalistic description of the Ca 2+ release due to IP3Rs activation which is particularly suitable for the incorporation into whole-cell dynamics models [39,40,30].
With respect to vascular ECs, some relevant modelling studies were carried out towards the characterization of endothelial Ca 2+ response to agonist stimulation. Wiesner et al. [41] proposed a mathematical model for describing the calcium dynamics originating in endothelial after thrombin intervention. This work provided a detailed description of the kinetics involved from the membrane receptor to the induced Ca 2+ release. The modelling of the IP3R channels followed a previous study carried out on rat basophilic leukemia cells [42]. The work by Wiesner al. laid the foundation for the development of other models [43,44,45,46], which were used for investigating the effects of different agonists and shear stress on the EC Ca 2+ dynamics. Since the activation of IP3Rs by exogenous agonists can only occur via signalling mediators in the cytosol, the predictive capacity of the corresponding model is hindered by uncertainties associated with the multiple kinetics involved.
To the best of our knowledge, there is no established modelling framework that can quantitatively describe the endothelial Ca 2+ dynamics following the IP3-induced Ca 2+ release from ER. In this work, we leverage on the Ca 2+ dynamics data obtained by direct stimulation of IP3Rs by Carter and Ogden [29], which allow us to exclude the uncertainties associated with the kinetics regulating the ligand-induced IP3 cytosolic production from the model. Despite the heterogeneous behaviour observed in different stimulations involving the IP3R pathways [29,28,13], we hypothesize the existence of a 'reference' model able to capture the mean Ca 2+ transient features across an EC population. We test this hypothesis by assessing the performance (in terms of fitting and validation) of a novel, minimalistic Ca 2+ dynamics model (alongside some more complex variants) which uses a gating function for regulating the IP3-induced Ca 2+ release in the cytosol. We then further speculate that cellular variability may be explained as deviation in some of its parameters defining cellular functions such as Ca 2+ sequestration and cytosolic buffering. The identified model will ultimately enable us to the estimate physical quantities which are difficult to measure directly from the experiment.

Experimental conditions and stimulation
The mathematical model is defined by considering the experimental conditions adopted by Carter and Ogden [29]. In brief, isolated ECs from porcine thoracic aorta were stimulated with different concentration of IP3 in whole cell patch-clamp settings (holding membrane potential of 0 mV). Since the EC holding membrane potential was 0 mV, the contribution of transmembrane ionic currents driven by the potential gradient can be ignored in the cytosolic Ca 2+ concentration balance (membrane hyperpolarization elevates the cytosolic Ca 2+ by increasing Ca 2+ influx, as shown in [29]). Caged IP3 (P-4 and P-5 isomers) was released in the cytosol by means of photolysis lasting ∼ 1 ms. Therefore the cell stimulation consisted in a pulse increase of the average cytosolic IP3 concentration (i), whose dynamics can be modelled with a first-order kinetics reaction equation: where i0 is the initial cytosolic IP3 concentration after 1 ms pulse of near UV light, k id is the IP3 dissociation constant due to the interaction with other cytosolic components, and tr is the release time.

Ca 2+ store release
For the considered dynamics, the variation in free Ca 2+ concentration across the cytosol is not considered significant and therefore the currents between the cytosol and the ER and the extracellular space can be expressed as function of the average free cytosolic Ca 2+ concentration (cc). For the sake of brevity, from here onwards we omit the 'free' in front of cytosolic Ca 2+ , whilst the bounded cytosolic Ca 2+ is indicated as 'buffered' cytosolic Ca 2+ . The Ca 2+ dynamics triggered by the activation of the IP3-dependent store channels is the result of the simultaneous interaction between different channels, pumps and exchanger expressed at membrane, cytosolic and ER levels. Here we capture the key characteristics of IP3R-mediated Ca 2+ release by employing a minimalistic approach which incorporates a gating variable. The endogenous Ca 2+ release flux from ER due to IP3Rs activation (J ef , see Figure 1) is made dependent on the cytosolic levels of IP3 and Ca 2+ through the following expressions [39,40,30] where h is a gating variable representing the availability/state of IP3 receptors, k ef is the nominal ionic flux, τci is the time constant governing the channel inactivation, whist Kia, Kca, Kci and exponents nia, nca and nci are the coefficients of the corresponding Ca 2+ -dependent activation and inactivation functions. Here, due to the lack of experimental evidence, we assume that Kci and τci are i-independent and cc-independent, respectively. Eq. (2) implies Figure 1: Key components of the cellular Ca 2+ dynamics model. J ef : ER IP 3 -sensitive Ca 2+ release flux; J le : ER Ca 2+ leakage flux to cytosol; J se : ER Ca 2+ sequestration flux; J ex : Ca 2+ extrusion flux from cytosol; J bc : Ca 2+ buffering flux.
assuming the interaction between IP3 and IP3Rs instantaneous [47], and that the Ca 2+ release flux is proportional to the fraction of open channels, which is inline with previous work [41]. In the considered experiment, the timelength associated with the Ca 2+ store release is extremely short (∼ 1 s) and therefore it is reasonable to assume that the variation in Ca 2+ concentration in the ER (cs) is minimal and does not impact the Ca 2+ current from ER. Consequently, the membrane Ca 2+ influx associated with store-depletion can also be considered negligible. We highlight that the IP3-mediated Ca 2+ release is not made proportional to the difference in average Ca 2+ concentration between the store and the cytosol because the latter is not necessarily representative of the Ca 2+ gradient nearby the channels. The Ca 2+ leakage from the ER to the cytosol (see Figure 1) is accounted for as where k le indicates the permeability of the ER membrane while cs, as explained above, is assumed to be constant.

Ca 2+ sequestration, extrusion and buffering
Ca 2+ re-uptake into ER (Jse) occurs via SERCA pump while Ca 2+ extrusion from the cytosol (Jex) is operated by membrane proteins such as PMCA and NCX (see Figure 1). These can generally be described as Ca 2+ -activation functions [41,46,30] Jse = Vse c nse where Vse and Vex are the nominal fluxes and Kse, Kex, nse and nex are the activation function coefficients corresponding to the Ca 2+ sequestration and extrusion from the cytosol. More complex formulations for these pumps and exchangers can be found in literature [46,30]. Despite such pumps and exchangers playing a physiologically essential role in the Ca 2+ depletion from the cytosol, in the considered experiment their contribution to the cytosolic Ca 2+ balance is limited [29]. In ECs, the Ca 2+ extraction rate of the membrane pumps remains low also because their characteristic time for full activation is considerably higher than the considered time-frame [25,48]. Due to this, together with the fact that the functioning of such proteins in ECs is poorly characterized and distinguished in literature, we assume that all Ca 2+ removal mechanisms from the cytosol can be incorporated together under only one flux, whose activation depends solely on cytosolic Ca 2+ concentration (cc). To avoid introducing new notation, we refer to this contribution with the flux Jse introduced in Eq. (5) whilst the flux Jex introduced in Eq. (6) will be considered negligible for the cytosolic Ca 2+ balance. Different types of intracellular buffers, such as calmodulin and fluorescence dye, act by dampening the variation in free cytosolic Ca 2+ concentration. There are different ways for including this effect into the model. A classical approach treats the (averaged in space) buffered cytosolic Ca 2+ concentration (c b ) as an independent variable, whose rate of change is determined via [41,46] where k ba and k bd are, respectively, the association and dissociation rate constants for Ca 2+ binding to cytosolic proteins whilst BT is the concentration of the total buffering proteins present in the cytosol. However, in the considered experiments [29], the variation in cytosolic free Ca 2+ appeared proportional to the change in the total (free plus buffered) cytosolic Ca 2+ . Therefore, Ca 2+ buffering may be also approximated as an instantaneous process and simply accounted for by means of a parameter β cb representing the fraction of non-buffered Ca 2+ with respect to the total. In this case, the rate of change in the cytosolic Ca 2+ concentration can be expressed as where β cb can be assumed to be constant or, alternatively, defined through the following expression [49,50] in which k * bd is the corresponding Ca 2+ dissociation rate constant. If the Ca 2+ buffering is instead treated as a dynamic process, the ionic mass balance yields

Solving the Ca 2+ dynamics
The solution of the Ca 2+ dynamics can be obtained, depending on the adopted Ca 2+ buffering approach, by solving either the system of Eq. (3) and (8) for the variables h and cc or the system of Eq. (3), (7) and (10) for the variables h, cc and c b . The time evolution of IP3 concentration i is considered independent from cytosolic Ca 2+ concentration and therefore Eq. (1) is computed before the Ca 2+ dynamics. In this study we consider five Ca 2+ dynamics model variants, each of them governed by a set of parameters to be determined (see Table 1). The system of ODEs representing the Model variant Ca 2+ dynamics Eq. Parameters to be determined Baseline (3) and (8) k le , τ ci , k ef , β cb , V se , K se , n se , k id Civar (3) and (8) k le , τ ci , k ef , β cb , V se , K se , n se , k id , K ci , n ci Ci&cavar (3) and (8) k le , τ ci , k ef , β cb , V se , K se , n se , k id , K ci , n ci , K ca , n ca β cb var (3) and (8) (3), (7) and (10) k le , τ ci , k ef , k bd , V se , K se , n se , k id , B T , k ba Ca 2+ dynamics is integrated over time by using the function 'odeint' from the python library NumPy, with time-step ∆t=0.001 s (this value was set after appropriate time-convergence analysis).

Model parameter identification
Here we present the methodology for identifying the set of parameters for each model variant that best fit the mean EC population features. We devise an optimization problem in which the solution is the set of model parameters that minimizes the discrepancy (quantified by a cost function) between the features of the simulated Ca 2+ transient and the corresponding measurements across an EC population. A set of model parameters is accepted only if the error (see Section 2.2.1) with respect to the experimental mean population data is below a certain 'ad-hoc' threshold. Four features of the Ca 2+ concentration transient are considered for the fitting: the 'Peak', 'Max flux', 'Termination rate per unit of flux' (see Figure 2) and 'Flux ratio' obtained from two consecutive stimulations with interpulse interval ∆tinter (see Figure 4 in [29]). The first three features were recorded across the IP3 stimulation range 0.2-16.0 µM, whilst the latter feature was evaluated with a first i impulse at 0.8 µM, followed by a second impulse at 0.7 µM, as documented in [29]. The feature Flux ratio can be numerically evaluated as ratio between the values of the gating variable h obtained immediately (1 time-step ∆t) after the two consecutive stimulations (separated by a interpulse ∆tinter). According to the experimental study, 16.0 µM of IP3 is high enough to saturate all the corresponding store receptors, and therefore no IP3 stimulation beyond such level is considered. Details on the feature extraction from Ca 2+ concentration time-recordings are reported in the Appendix.

Cost function
The cost function Ltot to be minimized is made of two components (Ltot=L data +L sum con ). L data is the sum of all relative errors of the simulated features (ŷ) with respect to average values from experimental data (ȳ):  Table 2: Imposed constraints in the optimization procedure. ∆c c,fin is the difference in cytosolic Ca 2+ concentration between the final (t > 20.0 s) and initial time step. ∀t: at any time; ∀j: for any IP 3 stimulation level.
where κ is the penalty parameter (1.0e6) and ∆ycon defines the distance of the variable/feature with respect to the allowed interval boundary. Set of parameters that yield a 'First decline rate' (see Figure 2) lower than -5.0 µM/s are also discarded in the optimization process. According to measurements (see Figure 3 in [29]), this latter feature is sparse with respect to Peak, and therefore is used only for discarding undesirable solutions.

Optimization algorithm
The parameter set is identified by minimizing the cost function through the covariance matrix adaptation evolution strategy (CMA-ES) [51,52]. The search stops when either the cost function goes below a 'ad-hoc' threshold (chosen to be 3.5 after empirical evaluation) or the number of iterations employed by the search-algorithm go above a pre-set maximum number of iterations (500). Given the high number of model parameters and due to the stochastic nature of the search algorithm, the optimizer may lead to different optimal set of parameters for the same model variant. Therefore, for any model variant, multiple cases (1,000) are used for identifying the most likely space of parameters. Since any considered model parameter (ξ) can be only positive, we express it as ξ = ξ0e α , where ξ0 is a constant (kept fixed during the optimization process) whilst α is the term through which the parameter ξ is varied. At the start of the search α is set to zero for each parameter and then updated by the optimization algorithm.

Results
Here we target the identification of the parametric space for each model variant, and the assessment of their performance by employing the quantitative experimental data from the study Carter and Ogden [29]. For doing this, different model variants are considered (introduced from the simplest to the most complex). This allows to identify the model components which are essential for capturing the experimentally-observed EC behaviour. The model is then used for making predictions on how EC variability may impact the IP3-induced Ca 2+ dynamics.

Baseline model
We start by considering the simplest model, defined as 'Baseline' model for which Ca 2+ buffering can be approximated as an instantaneous and Ca 2+ -independent process (β cb is constant). Furthermore, for this model, we assume that Ca 2+ affects the IP3Rs in the same manner as reported in [46] (we fix Kci=1.0 µM, nci=3.8, Kca=0.0 µM). The coefficients for the i-dependent activation function of the ER Ca 2+ flux are directly taken from [29] (we fix Kia=1.6 µM, nia=3.8). As anticipated earlier, the Ca 2+ concentration in the ER (cs) can be considered constant, equal to 3,000 µM. The model parameters (k le , τci, k ef , β cb , Vse, Kse, nse, k id ) are identified by solving the optimization problem described in Section 2.2. Around 15 % of the cases were discarded because their cost function did not satisfy the imposed error threshold within the maximum number of iterations of the optimizer (500). Figure 3 compares the curves of 10 different model parameter sets (obtained from 10 solutions of the optimization problem) and the solution obtained from the model with the geometric mean of (all) the parameters against the mean population experimental values. The geometric mean (for a given model variant) is obtained by averaging the variable term of the parameter α (and not the final parameter value ξ) across all the parameter sets obtained as solution of the optimization problem. Table 3 reports the initial/guess value and geometric mean for each parameter of the Baseline model, obtained by considering 847 cases/solutions. -  EC cytosolic Ca 2+ transient is known to depend on the Ca 2+ inactivation time constant (τci) and the IP3 dissociation rate (k id ), which for the Baseline model appear to fall in well defined range of values.
The predictive capacity of the Baseline model is also assessed by considering another Ca 2+ transient feature, the First decline rate, which reflects the rate of Ca 2+ removal from the cytosol after stimulation. Carter and Ogden [29] extracted this feature from experimental time recordings of Ca 2+ concentration for different IP3 stimulation levels and reported it as function of Ca 2+ Peak. This feature appears extremely heterogeneous across the EC population and therefore the extracted average trend can be considered less representative than the other population features. Figure  5 shows a comparison between the model predictions and experimentally reported values. Despite its simplicity, In the same figure we report the Pearson correlation coefficient between this feature and the model parameters, evaluated for each stimulation level. As expected, the parameter governing the Ca 2+ sequestration process Kse stands out, indicating that pumps with high sensitivity with respect to cytosolic Ca 2+ level (low Kse) are responsible for sharp Ca 2+ concentration drops. Altogether these results indicate that the Baseline model is able to capture the mean population features of the considered Ca 2+ dynamics, constituting a first reference model.

Ca 2+ effect on IP 3 Rs activation
The Ca 2+ dynamics under analysis is initiated by an IP3-induced stimulation of the intracellular Ca 2+ stores. Here we numerically investigate how the dependency of IP3-induced ER Ca 2+ release on cytosolic Ca 2+ affects the overall EC dynamics. For doing this we consider, with respect to the baseline model, the following variants i) 'variable' Ca 2+ inactivation component (Civar) and ii) 'variable' Ca 2+ inactivation and activation components (Ci&cavar). The parameters of these Ca 2+ -dependent functions are included into the model parameter set to be identified through the optimization algorithm (as in Section 3.1). For only a fraction of the cases (672/1000 for Civar, 642/1000 for Ci&cavar) a solution for the optimization problem was found (cost function below the threshold). The fittings for these two model variants against experimental curves and their parameter spaces are reported in the Appendix (see Figures 12 and  15, and Figures 13 and 16, respectively). With respect to the Baseline model, the variant Civar provides a slightly better estimate in terms of averaged First decline rate (see Figure 14). Table 4 reports the geometric mean (geometric standard deviation) of the parameter set cases across the considered model variants. For most of the simulated cases of this model variant, the sensitivity (Kci) of the inactivation function seems to be just below 1.0 µM (see Figures 13  and 16 in the Appendix), whilst the distribution of the exponent of the Hill equation (nci) appears to be more sparse.
On the contrary, a function for ER Ca 2+ release which accounts for both Ca 2+ -activation and -inactivation effects (Ci&cavar) does not improve the discrepancy with respect to experimental recordings (fitting and validation in the Appendix, respectively Figures 15 and 17). The Ca 2+ inactivation time constant (τci) does not significantly change across the considered models while the IP3 dissociation rate (k id ) appears to be reduced for the two more complex variants.
The Ca 2+ -dependency of the IP3-induced ER Ca 2+ release was experimentally assessed by Carter and Ogden by considering different pre-stimulation cytosolic Ca 2+ levels. The predictions from the Baseline model and the other variants (obtained from the previous fittings) are compared against these experimental data to identify the most appropriate relationship between IP3-induced ER Ca 2+ release and cytosolic Ca 2+ concentration (see Figure 6). The Baseline and Civar models, despite overestimating the 'Normalized max flux' for cytosolic Ca 2+ concentration below    Figure 6: Role of Ca 2+ on IP 3 Rs function: experimental vs numerical data. In the experiment by Carter and Ogden [29], the cytosolic Ca 2+ concentration was increased by membrane hyperpolarization before agonist intervention. The 'Normalized max flux' was evaluated as a ratio between the Max flux obtained at the imposed (through membrane hyperpolarization) cytosolic Ca 2+ concentration and the Max flux obtained at basal (before membrane hyperpolarization) cytosolic Ca 2+ concentration.
1.0 µM, are able to capture the main trend of the experimental data. On the other hand, the introduction of a biphasic function accounting for the Ca 2+ activation component (Ci&cavar) leads to very similar results. We conclude that the inclusion of the Ca 2+ -dependent activation component into ER Ca 2+ release model does not improve its predictive capacity (and therefore it is discarded in the following) while the Baseline model can be considered representative of the model Civar.

Role of cytosolic Ca 2+ buffering on Ca 2+ elevation
Here we investigate the role of cytosolic Ca 2+ buffering on the model's dynamics. Two variants, with respect to the Baseline model, are considered: i) the buffering parameter β cb depends via Eq. (9) on cytosolic Ca 2+ concentration (β cb var) and ii) buffering considered as a dynamic process (Dynbuff), with the buffered cytosolic Ca 2+ concentration obtained by solving Eq. (7). Note that in the latter case the cell dynamics is obtained by solving Eq. (7) in conjunction with Eq.   Table 5: Geometric mean of the parameter set cases for the considered model variants.
non-constant buffering effect (associated with the two variants) reduces the necessary Ca 2+ sequestration and the maximum store Ca 2+ efflux (Vse and k ef are significantly smaller for the two variants). Figure 7 shows the Ca capacity of the cytosolic proteins reduces the Ca 2+ Peak but it does not significantly alter the secondary recovery phase (after 4 s). Furthermore, the transient of the gating variable h is similar across the model variants, indicating that the (modelled) store Ca 2+ -dependent inactivation is not remarkably affected by the adopted buffering approach. In the Baseline model, buffering plays a more limited action in compensating the Ca 2+ increase after IP3Rs activation. In this case, indeed, the Ca 2+ sequestration flux (Jse) has a much higher peak than in the buffering variants. These results suggest that, upon IP3Rs stimulation, two cells with remarkably different Ca 2+ sequestration capacity can still exhibit a very similar Ca 2+ response, as long as there is a cytosolic Ca 2+ buffering compensation.
Despite the Ca 2+ store IP3-sensitivity function for porcine aortic EC being well-defined in [29], it is important to assess its role on the overall cell dynamics. We therefore numerically explore, for all the buffering approaches, how the variability in IP3Rs' sensitivity to IP3 may affect the EC Ca 2+ response to different IP3 stimulation levels (see Figure 8). As expected, an increase in IP3Rs' sensitivity to IP3 (corresponding to a Kia reduction) makes the cell more responsive to low concentration of stimulating agent. The Ca 2+ peaks obtained with IP3 stimuli 1.4 µM and 10.0 µM are relatively similar (∼30% difference) for Kia=1.0 µM, whilst they are significantly different (>100% difference) for higher Kia. Also in this case, no substantial changes in the Ca 2+ transients are recorded between the different buffering approaches. Given the results of the comparison with experimental data, the reported variants (with more sophisticated buffering approach) may also serve as reference model. However, the increased model complexity (larger

Role of IP 3 dissociation and Ca 2+ inactivation time constants on Ca 2+ dynamics
The IP3 dissociation constant (k id ) and the Ca 2+ inactivation time constant (τci) are expected to play an important role on the cell dynamics. The impact of the former on the Ca 2+ transient is assessed by means of Figure 9 (for different stimulation levels and buffering approaches). Surprisingly, the IP3 rate of dissociation plays an important role on the Ca 2+ response amplitude (Peak) only upon low IP3 stimulation levels. For higher IP3 stimuli, the k id impacts only the second part of the Ca 2+ decline/recovery phase whilst the previous part of the transient seems unaffected. Upon high stimulation levels, cytosolic Ca 2+ concentration decline faster for high k id (0.38 s −1 ) which is, from a qualitative point of view, more inline with the experimentally recorded single-cell traces (see Figure 2B-C in [29]) than for the intermediate k id case (0.19 s −1 ). Values of k id in the range 0.15-0.19 s −1 were identified through the parameter optimization procedure as geometric mean for the Baseline model and the other variants as they guarantee a better fit over the whole range of IP3 stimulation. However, in case the experimentally-derived mean population Peak [29] was overestimated for low IP3 stimulation levels (or a more complex IP3 dynamics takes place), a higher IP3 dissociation constant would represent a more appropriate choice. While the impact of k id on the Ca 2+ transient seems the same across all buffering approaches, the same cannot be said regarding the parameter τci (see Figure 10). For the Baseline model and upon medium-high stimulation levels, the Ca 2+ inactivation time constant can significantly modulate the amplitude of the signal but the other two variants appear almost insensitive with respect to this parameter. Interestingly, for different buffering approaches, almost the same geometric mean of τci was found (∼0.045 s). The Ca 2+ -dependent gating variable (h) sharply drops after the stimulation in a similar fashion across all buffering approaches (see Figure 11 in the Appendix). The discrepancy in τci-sensitivity between the Baseline and the other variants may be explained as result of the significant difference in the nominal ionic flux (k ef ). Since the magnitude of k ef for the Baseline model is significantly higher than for the variants, a variation in IP3Rs inactivation time will have a much bigger impact on the cytosolic Ca 2+ balance.

Discussion and concluding remarks
Through this work we introduced a minimalistic modelling framework for describing the EC Ca 2+ dynamics induced by IP3Rs activation. For doing this, we hypothesized the existence of a reference model able to the capture the average Ca 2+ transient features from the cell population used in the study by Carter and Ogden [29]. This experimental protocol enabled the characterization of the Ca 2+ efflux following the IP3Rs activation across a range of stimulation levels. The proposed structure of Baseline model, which considers the cell variables homogeneously distributed in space, allowed a more than satisfactory fitting against four features of the time-recorded Ca 2+ traces. Furthermore, the validation against the measured First decline rate and Normalized max flux corroborated its predictive capacity. These results together indicate that IP3-induced EC dynamics, despite involving diffusive phenomena across the cytosol, can be captured with a good level of accuracy by a lumped modelling approach. Model variants led to similar results under most of the aspects, indicating that the minimalistic Baseline model can constitute a reference framework for capturing the key features of the considered Ca 2+ dynamics. Whilst variations in the description of IP3Rs Ca 2+ -dependent inactivation (Civar) can reduce the discrepancy with respect to measurements, the inclusion of the IP3Rs Ca 2+ activation effect (Ci&cavar) does not substantially improve the IP3Rs Ca 2+ -dependent function in terms of accuracy. Supported by findings on other types of cells, more complex mathematical constructs could be used in future for describing the bidirectional link between IP3, Ca 2+ and Ca 2+ store release. Mak et al. [53] investigated the IP3 dependence of the Ca 2+ sensitivity (of the ER receptors) at single channel level in isolated Xenopus oocyte nuclei by employing a biphasic Hill equation. Whilst the Ca 2+ -activation effect appeared to be independent from i, the dissociation constant describing the Ca 2+ inhibitory role increased with higher i. The IP3-dependent effect on the Ca 2+ inhibitory role may be accounted for by defining the parameter Kci as dependent (via activation function) on IP3 concentration i. The inclusion of more sophisticated Ca 2+ buffering did not significantly alter the shape of Ca 2+ response, even though the Ca 2+ sequestration function is reduced. The proposed framework was meant to provide a detailed description of the Ca 2+ efflux due to IP3Rs activation by minimizing the model complexity. For this reason, a very simplified description of the cytosolic Ca 2+ removal mechanisms was adopted. Their description accounts for a simple dependency on cytosolic Ca 2 + concentration via activation function and does not take into consideration the delay in the activation of such cell components nor the interaction with other ions. This is expected to affect, at least in part, the shape of the second phase of the Ca 2+ recovery transient. Considering two separate pumps with different Ca 2+ sensitivity (Jex different from zero) in the model identification did not lead to any significant improvement (results not shown). Since the activation of such cell components affects in a minimal way the Ca 2+ efflux following the IP3Rs activation, its detailed analysis is considered beyond the scope of this work.
Through the proposed approach, we were also able to explore parameters characterizing the IP3-induced store Ca 2+ release. The IP3 dissociation constant (k id ) within the EC cytosol obtained from the parameter optimization process appeared to be significantly lower than what estimated/used in other cell models (1.0 s −1 in [54], 1.25 s −1 in [50], 2.0 s −1 in [41]). Here we showed that k id affects the peak of the Ca 2+ transient only for low stimulation levels. As expected, higher values of IP3 dissociation rate are associated to a faster cytosolic Ca 2+ decline, but in this case the discrepancy with the experimental data for low IP3 stimuli increases. Despite the evidence for the regulatory role of Ca 2+ and PLC on the cytosolic IP3 dynamics [54], this feedback mechanism is not mentioned in the considered experiment [29] and its role remains uncertain. Furthermore, the manipulation of the system by introducing fluorescent dye and/or altering the electrolytic balance may also have a significant impact on the IP3 dissociation rate. To justify the reported results, one could speculate that the lower IP3 dissociation constant (with respect to other cell types) may be due to the fact that, since IP3 represents in ECs the primary pathway for agonist stimulation, it needs to be metabolized by cytosolic proteins at slower rate to enable its signalling at low concentrations. This might be the case, as well as a more complex IP3 dynamics taking place, and therefore future studies are necessary for addressing this question. This study also led to a first estimate of the Ca 2+ -inactivation time constant associated to IP3Rs, which appeared to be consistent across all the model variants.
In the considered settings, the membrane potential was kept equal to 0 mV throughout the experiment with the aim to isolate the store Ca 2+ efflux. Behringer and Segal [55] investigated the controversial role of membrane potential on Ca 2+ response to ACh in endothelial tubes. They showed that, under submaximal activation of muscarinic receptors, membrane potential modulates Ca 2+ influx through the plasma membrane channels in accord with the electro-chemical driving force. On the other hand, previous experimental evidence [56] indicated that, for isolated ECs from small resistance vessels, membrane potential does not seem to play a significant role on sustained muscarinic agonist-induced cytosolic Ca 2+ elevation. In pressurized rat mesenteric arteries, Ca 2+ influx through TRPV4 channels plays a primary role in EC cytosolic Ca 2+ elevation [55,13]. However, despite their activation is associated with mechanisms of Ca 2+ -induced Ca 2+ -release at IP3Rs, these do not seem to significantly contribute to ACh-evoked relaxations [13].The complexity and uncertainty associated with the role of membrane potential in the regulation of endothelial cytosolic Ca 2+ motivated our choice of considering the simplified setting by Carter and Ogden [29] where membrane potential was kept constant equal to 0. Despite IP3Rs activity depends primarily on cytosolic IP3 and Ca 2+ levels, it may be affected by interaction with other organelles such as mitochondria [57] and other factors such as oxidative stress [58]. Furthermore, Wilson et al. suggested that also pressure, by altering the cell geometry, can affect the Ca 2+ diffusive environment near IP3 receptor micro-domains, limiting the IP3-mediated Ca 2+ signals as pressure increases. Therefore, future studies investigating the interaction between EC Ca 2+ store release and other cytosolic and membrane components are warranted.
This work lays the foundation for future experimental and modelling studies targeting the analysis and quantification of Ca 2+ dynamics in ECs and its effects on blood vessel contractility and permeability. The proposed model could indeed be combined with other models representing SMCs' activity [59,60] for quantifying the interaction between functionally different vascular cells, the resulting wall stress level and lumen deformation. These theoretical findings may also contribute to the development of new models for studying the blood-brain barrier function [61]. Furthermore, this work might help in the analysis of Ca 2+ intercellular communications and Ca 2+ wave propagation across cellular networks such as in [62]. To conclude, the proposed minimalistic theoretical model aims to be a reference in-silico tool for deciphering vascular function.

Appendix
Ca 2+ transient features extraction

Experimental data
Here we report the method for extracting the statistics from the datasets reported in [29]. Digitization of data from plot images was carried out by employing the open-source tool 'WebPlotDigitizer'.
Max flux The population data reporting the maximum Ca 2+ release from stores against IP3 (i) concentration was fitted by Carter and Ogden with the following Hill equation (see Figure 3B in [29]): where k4=200.0 µM/s, nia=3.8 and Kia=1.6 µM. This relationship was obtained by considering IP3 stimulation up to 16.0 µM. Beyond this concentration, the max Ca 2+ release response recorded was extremely variable among cells and therefore was not included.
Termination rate per unit of flux Carter and Ogden fitted the measurements of Termination rate per unit flux against IP3 concentration (i) with a single site binding curve (see Figure 5C in [29]) where k6=0.5 µM −1 k7=0.38 µM and k8=0.1 µ M −1 .
First decline rate Here the dependency between the first rate of decline and the Ca 2+ Peak (Pc) was established by fitting with a quadratic function the data reported in Figure 3A in [29] y = k12P 2 c + k13Pc + k14, where k12=-0.003 µM −1 s −1 , k13=-0.127 s −1 and k14=0.218 µM/s. Role of time constant on Ca 2+ -dependent gating variable

Simulated data
In Figure 11 we report how the Ca 2+ inactivation time constant τci affects the transient of the Ca 2+ -dependent gating variable for different buffering approaches.

Model variants fitting, parameter space and validation
Below the plots of the experimental fitting, parameter space and experimental validation for all the considered model variants (except for the Baseline) are reported.