HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 293/1/C277    most recent 00542.2006v2 00542.2006v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (1) Citing Articles via Google Scholar Google Scholar Articles by Silva, H. S. Articles by Tsoukias, N. M. Search for Related Content PubMed PubMed Citation Articles by Silva, H. S. Articles by Tsoukias, N. M.

VASCULAR BIOLOGY

A mathematical model of plasma membrane electrophysiology and calcium dynamics in vascular endothelial cells

Department of Biomedical Engineering, Florida International University, Miami, Florida

Submitted 22 October 2006 ; accepted in final form 20 March 2007

	ABSTRACT

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX: MODEL EQUATIONS GRANTS REFERENCES

 
Vascular endothelial cells (ECs) modulate smooth muscle cell (SMC) contractility, assisting in vascular tone regulation. Cytosolic Ca²⁺ concentration ([Ca²⁺]_i) and membrane potential (V_m) play important roles in this process by controlling EC-dependent vasoactive signals and intercellular communication. The present mathematical model integrates plasmalemma electrophysiology and Ca²⁺ dynamics to investigate EC responses to different stimuli and the controversial relationship between [Ca²⁺]_i and V_m. The model contains descriptions for the intracellular balance of major ionic species and the release of Ca²⁺ from intracellular stores. It also expands previous formulations by including more detailed transmembrane current descriptions. The model reproduces V_m responses to volume-regulated anion channel (VRAC) blockers and extracellular K⁺ concentration ([K⁺]_o) challenges, predicting 1) that V_m changes upon VRAC blockade are [K⁺]_o dependent and 2) a biphasic response of V_m to increasing [K⁺]_o. Simulations of agonist-induced Ca²⁺ mobilization replicate experiments under control and V_m hyperpolarization blockade conditions. They show that peak [Ca²⁺]_i is governed by store Ca²⁺ release while Ca²⁺ influx (and consequently V_m) impacts more the resting and plateau [Ca²⁺]_i. The V_m sensitivity of rest and plateau [Ca²⁺]_i is dictated by a [Ca²⁺]_i "buffering" system capable of masking the V_m-dependent transmembrane Ca²⁺ influx. The model predicts plasma membrane Ca²⁺-ATPase and Ca²⁺ permeability as main players in this process. The heterogeneous V_m impact on [Ca²⁺]_i may elucidate conflicting reports on how V_m influences EC Ca²⁺. The present study forms the basis for the development of multicellular EC-SMC models that can assist in understanding vascular autoregulation in health and disease.

microcirculation; vascular tone regulation; calcium influx pathway(s), plasma membrane Ca²⁺-ATPase

Mounting experimental evidence shows that EC V_m hyperpolarization itself, independently of its impact on EC Ca²⁺, is an important signal in resistance vessels. This electrical signal is generated by Ca²⁺-activated K⁺(K_Ca) currents and most likely comprises a key component of the EDHF response, which is transmitted to adjacent SMCs via myoendothelial gap junctions to cause vessel relaxation (14, 22, 52, 53). Another candidate for EDHF is simply K⁺ expelled by ECs into the vascular interstitial gap affecting K⁺-sensitive channels and transporters in both SMCs and ECs. The EC layer is most likely the main electrotonic current conduction route (for signals generated by either SMC or EC stimulation) along the vessel wall during conducted vasomotor responses (24, 31, 32, 65, 73). In addition, EC K_Ca channel inhibition can greatly affect arterial vasomotion (53), suggesting an important role of EC electrophysiology in this phenomenon.

Isolated ECs can be classified into two subtypes according to their resting V_m, namely, K-type and Cl-type (56). K-type EC resting V_m levels fall between –70 and –60 mV, which is close to the Nernst potential of K⁺ (E_K) and thus indicates a dominant K⁺ membrane conductance, mainly due to the inward rectifier K⁺ (Kir) channel. On the other hand, Cl-type EC potential at rest is usually between –40 and –10 mV, which is close to the Nernst potential for Cl^– (E_Cl) and suggests a Cl^– conductance dominance under resting conditions. Experiments on isolated vessels suggest that ECs under hypoosmotic stress (which activates volume-sensitive Cl^– conductance and brings V_m toward E_Cl) will not be able to hyperpolarize in response to raised extracellular K⁺ concentration ([K⁺]_o), via Kir activation, and thus cannot relax the artery (25). Experimental data on isolated ECs, however, have shown that Cl-type and not K-type cells hyperpolarize in response to K⁺ challenge (55). To understand the complex V_m role in EC behavior and vascular tone regulation, the theoretical understanding of EC electrophysiology must be enhanced.

Previous EC mathematical modeling efforts have made important contributions toward describing both Ca²⁺ signaling and plasmalemmal electrical activity. Two models by Wiesner et al. (77, 78) simulated the response of human umbilical vein ECs (HUVECs) to thrombin stimulation and fluid shear stress. Wong and Klassen (80) developed a model describing both [Ca²⁺]_i and electrical responses of vascular ECs to shear stress. Korngreen et al. (44) successfully simulated Ca²⁺ transients of electrically nonexcitable cells following agonist washout and intracellular Ca²⁺ store depletion. Schuster et al. (67) modeled the changes in V_m and [Ca²⁺]_i following bradykinin stimulation of coronary artery ECs.

However, ECs regulate Ca²⁺ entry and V_m by expressing an abundant and diverse collection of plasmalemmal ion channels (55), which are for the most part absent in previous EC models. In addition, considering the balance of the other major intracellular ionic species (i.e., Na⁺, K⁺, and Cl^–) is essential to modeling both single cell V_m behavior and cell-to-cell electrochemical coupling as Nernst potentials are concentration gradient dependent and gap junctions are nonselective, allowing the passage of small cytoplasmic ions. Consequently, there is a need to build on previous EC modeling work to provide a more biophysically detailed model of EC plasma membrane electrophysiology and further analyze the extent to which, if any, V_m affects EC intracellular Ca²⁺ dynamics.

In the last four decades, mathematical modeling of biological systems has contributed to the basic understanding of physiological behavior. For some organs (i.e., the heart), multiscale models have been developed that describe function at the macroscale level while incorporating mechanisms and events at the subcellular and molecular levels. A similar theoretical progress has not been paralleled in the vasculature. This mathematical model presents a first step toward this direction. The aims of this study were 1) to deliver a mathematical model that captures experimentally observed behavior of vascular ECs (and particularly ECs from rat mesenteric arterioles) and 2) to analyze how these cell responses emerge from the nonlinear interactions of individual cellular components.

	METHODS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX: MODEL EQUATIONS GRANTS REFERENCES

 
The present EC model can be divided into two components (Fig. 1). The first component represents the equivalent electrical circuit model of the EC plasma membrane (Membrane Electrophysiology), describing its electrophysiology (Fig. 1A). It was developed by identifying the main transmembrane currents affecting cell homeostasis and [Ca²⁺]_i dynamics in rat mesenteric ECs and other vascular ECs (1, 55, 72). The EC transmembrane ionic currents included in the model are Kir, volume-regulated anion channel (VRAC), small-conductance (SK_Ca) and intermediate-conductance K_Ca (IK_Ca), store-operated Ca²⁺ channel (SOC), nonselective cation channel (NSC), Ca²⁺-activated chloride channel (CaCC), Na⁺-K⁺-ATPase (NaK), plasma membrane Ca²⁺-ATPase (PMCA), Na⁺/Ca²⁺ exchanger (NCX), and Na⁺-K⁺-2Cl^– cotransport (NaKCl) system. Published experimental data were used to describe them mathematically using standard electrophysiological equations.

View larger version (40K): [in this window] [in a new window]   Fig. 1. Schematic diagram of model components and their interactions. The model can be divided into two components: plasmalemmal electrophysiology (A) and fluid compartment model (B). Model equations are given in the APPENDIX, and model parameters are given in Tables 1 and 2. Cm, membrance capacitance; ISOC, store-operated Ca2+ channel (SOC) current; INSC, nonselective cation channel (NSC) current; IVRAC, voltage-regulated anion channel (VRAC) current; ICACC, Ca2+-activated Cl– channel (CACC) current; IKir, inward rectifier K+ channel current; IKCa, Ca2+-activated K+channel (KCa) current; ISKCa, small-conductance KCa channel current; IIKCa, intermediate-conductance KCa channelcurrent; INaK, Na+-K+-ATPase (NaK) current; IPMCA, plasma membrane Ca2+-ATPase (PMCA) current; INCX, Na+/Ca2+ exchanger (NCX) current; ESOC, ENSC, ECl, and EK, Nernst potentials of SOC, NSC, Cl–, and K+, respectively; A, agonist; R, receptor; [Cl–]i and [Cl–]o, intracellular and extracellular Cl– concentration; [Na+]i and [Na+]o, intracellular and extracellular Na+ concentration; [K+]i and [K+]o, intracellular and extracellular K+ concentration; [Ca2+]i and [Ca2+]o, intracellular and extracellular Ca2+ concentration; Vm, membrane potential; CSQN, calsequestrin; ER, endoplasmic reticulum; SERCA, sarco(endo)plasmic reticulum Ca2+-ATPase; IP3, inositol (1,4,5)-trisphosphate; IP3R; IP3 receptor; PIP2, phosphatidylinositol (4,5)-bisphosphate; RyR, ryanodine receptor; CaM, calmodulin; NaKCl, Na+-K+-2Cl– cotransporter; EC, endothelial cell.

Model mathematical expressions are given in the APPENDIX, while model standard parameters and initial conditions are shown in Tables 1 and 2, respectively. Model components and parameter values were chosen to best describe rat mesenteric artery ECs. All abbreviations and symbol definitions can be found in the text and in Tables 1 and 2.

The EC plasma membrane electrical activity is modeled according to the classical Hodgkin-Huxley model as it is applied to a variety of single cell models, including ECs (37, 47, 67, 80, 83). The EC plasma membrane is thus regarded as a capacitor with its capacitance (C_m) shunted by several ionic currents that are described in detail in this section. Kirchhoff's current law describes dynamic V_m changes by:

(1)

where I_stim represents an external stimulation current and the other currents are carried via the aforesaid transmembrane ionic pathways and described below. The EC membrane expresses a cotransport system of Na⁺, K⁺, and 2 Cl^– (Fig. 1B), but since these ionic fluxes result in an electroneutral process, there is no net current generated and this is therefore excluded from Eq.1. NaKCl fluxes are taken into account, however, in the intracellular species balance equations (APPENDIX, I_NaKCl).

Kir channel current. Kir channels are considered key for resting cell V_m regulation, but their expression is rather heterogeneous among vascular EC types and some do not express it at all (55). I_Kir is presently modeled following Hodgkin-Huxley formalism for the instantaneous voltage-dependent gating of the channel's conductance (Eqs. A1–A3 in the APPENDIX) (50). Maximal Kir conductance is described by an increasing function of [K⁺]_o and fits well Kir conductance data from guinea pig coronary ECs (data not shown) (75). Equations A1–A3 were fitted to rat mesenteric EC Kir current-voltage (I-V) data in the –100 to +40 mV range (Fig. 2A, solid line) (17). Kir conductance and Kir velocity fitted values are in agreement with data from Ref. 75. The Kir I-V model predictions under physiological intracellular K⁺ concentration ([K⁺]_i) with control and elevated [K⁺]_o displayed the reported "crossover" effect (Fig. 2A, dashed and dotted curves) (34, 55).

View larger version (20K): [in this window] [in a new window]   Fig. 2. IKir and IKCa. A: whole cell current-voltage rescaled data (circles) (17), model fitting (solid line), and predictions (dashed and dotted lines) for IKirs at different [K+]i and [K+]o. B: linear (solid lines) model fittings to whole cell KCa rescaled data from porcine coronary artery ECs (10). IKCa data above +50 mV were excluded from the fittings ([K+]i = 125 mM and [K+]i = 5.4 mM).

Ca²⁺-activated Cl^– channel current. The physiological significance of CaCC in vascular ECs remains controversial, but CaCCs can function in feedback control loops for Ca²⁺ entry into the cell via membrane depolarization, antagonizing the action of K_Ca channels (33, 62). The model I_CaCC (Eqs. A6–A9) was based on and fitted to rescaled HUVEC data (Fig. 3A) (27). Calcium-dependent activation of CaCC was captured by a Hill function (Eq. A6). High positive V_m enhanced Ca²⁺-dependent activation, producing outward rectification of I_CaCC (Fig. 3A) (27, 33, 62). This V_m-dependent activation was modeled using a Boltzmann function with a large positive half-activation V_m (Eq. A8). The time constant () of the voltage activation was fitted with a Gaussian function (Eq. A9) (data not shown).

View larger version (19K): [in this window] [in a new window]   Fig. 3. ICaCC and ISOC. A: whole cell current-voltage rescaled CaCC data (circles and triangles) (27) and model fits (solid and dashed lines). B: negative of ISOC data (circles) from mouse aortic ECs (29) and model fit (solid line) versus [Ca2+]o. The inset shows predicted Na+ and Ca2+ components of ISOC (ISOC,Na and ISOC,Ca; dashed and dashed-dotted lines).

View larger version (29K): [in this window] [in a new window]   Fig. 7. Overall EC model control resting condition [time (t) < 20 s] and Vm, intracellular species, transmembrane current, and IP3-sensitive store current control dynamic responses to 100 s of agonist stimulation (at t = 20 s). Plots display the dynamic behavior of Vm (A); ISOC,Na and ISOC,Ca (B); IPMCA, INCX, and INaK (C); [Ca2+]i and IP3 responses (D); K+, Na+, and Ca2+ components of INSC (INSC,K, INSC,Na, and INSC,Ca; E); IIP3R and IP3-sensitive store ISERCA (ISERCA,IS) and leak currents (Ileak,IS) (F); intracellular ([Na+]i, [K+]i, and [Cl–]i) and IP3-sensitive store ionic concentration ([Ca2+]IS) changes (G); IKCa and IKir (H); and ICaCC, IVRAC, and Cl– component of the NaKCl system current (INaKCl_Cl) (I). IP3 and [Na+]i values are rescaled by a factor of 10 (10x) and [Ca2+]IS values are rescaled by a factor of 50 (50x) for plotting convenience.

NSC current. NSCs permeable to Ca²⁺ (39, 40) constitute Ca²⁺ entry pathways other than store operated in ECs (55). A Ca²⁺-permeable NSC is present in rat mesenteric ECs and has been characterized at the single channel level. It is permeable to major intracellular cations (i.e., Ca²⁺, Na⁺, and K⁺), displays PKG sensitivity, and plays a key function in the EDHF response of small mesenteric arterioles to KT-5823 (a potent PKG inhibitor) (23). However, since the study (23) did not measure whole cell I_NSC and its gating mechanisms were not thoroughly analyzed, these data were not used in the present I_NSC formulation. The model I_NSC contributes to the balance of Na⁺, K⁺, and Ca²⁺ and provides an alternative pathway for Ca²⁺ entry into the cytosol, as observed experimentally (55). The formulation developed here to model I_NSC (Eqs. A16–A20) is based on rescaled I-V data from a constitutively open NSC found in rabbit aorta ECs (63). This study (63) found a basally active I_NSC having a permeability ratio of 1:0.40:0.18 for K⁺:Na⁺:Ca²⁺, respectively, and this current was the major resting V_m determinant in these cells, which lack Kir channels. Like in SOC, external Ca²⁺also had a blocking effect on I_NSC Na⁺ permeability (P_NSC,Na), whose description fits a Hill relationship similar to the one used to model P_SOC,Na (Eq. A20). I_NSC was described by the sum of three GHK equations (Eq. A19), one for each permeant species, namely, Na⁺ (I_NSC,Na), K⁺ (I_NSC,K), and Ca²⁺ (I_NSC,Ca). GHK has been used before to model rat brain endothelial NSC (18). Additionally, I_NSC is not gated by store depletion or [Ca²⁺]_i, as evidenced in calf pulmonary artery ECs (CPAECs) (55). Model I_NSC,Na is the major inward NSC current (data not shown), as found in EC experiments (63).

NCX current. The NCX is an electrogenic transporter of Na⁺ and Ca²⁺ with a stoichiometry of 3:1, respectively, as found in different cell types including ECs (6, 20, 30, 77, 83). Although the presence of NCX in ECs has been confirmed, its involvement in the modulation of Ca²⁺ signaling is controversial (55). In CPAECs, NCX contributes to cytosolic Ca²⁺ removal as a low-affinity system due to the high [Ca²⁺]_i needed for significant activation, and, in effect, this exchanger counteracts large and rapid rises in EC [Ca²⁺]_i during early stages of stimulation (68). In cultured bovine pulmonary artery ECs, NCX was found to be voltage dependent but not responsible for resting [Ca²⁺]_i maintenance (64). The model's I_NCX formulation (Eqs. A21–A23) is based on previously reported equations (30, 47, 76). The I_NCX expression given in Ref. 47 has fewer parameters from the one given in Ref. 76, and it was adapted in this study. A Hill-type instantaneous [Ca²⁺]_i-dependent activation term is also included in this I_NCX description (Eq. A21), as used in modeling NCX in cardiac myocytes and pancreatic -cells (30, 76). I_NCX equations were fitted to rescaled data from mouse pancreatic -cells (30).

NaK current. The NaK pump has been found in various EC types under a multitude of experimental protocols, and it represents the main transmembrane ionic pump present in virtually all mammalian cell membranes (1, 70). This pump has a coupling ratio of 3Na⁺:2K⁺, and its inhibition causes a 5 to 10 mV depolarization in cultured ECs (1, 70). NaK activity maintains intracellular Na⁺ concentration ([Na⁺]_i) and [K⁺]_i low and high, respectively, thus keeping physiological Na⁺ and K⁺ gradients across the plasmalemma (70). NaK activity has been described mathematically in detail, taking into account the metabolic state of the cell (ATP levels), which is suitable for modeling studies concerned with the effects of ischemia on the functionality of cardiac myocytes (70). However, a simpler I_NaK formulation is applied as the effects of low ATP in EC behavior were not pursued at present. The I_NaK formula used here (Eq. A24) was based on previous models (47, 83) and fitted to rescaled NaK I-V data from human embryonic kidney cells expressing the rat brain NaK ₁-subunit (84).

NaKCl cotransport flux. The NaKCl mechanism is electroneutral, pivotally involved in the regulation of resting cell volume and tissue permeability, and believed to constitute the major influx pathway for K⁺ in vascular ECs, even surpassing the NaK contribution (1, 41, 58). It also largely contributes to EC Cl^– transport (21). In cultured bovine corneal ECs, the NaKCl function needs extracellular Na⁺, K⁺ and Cl^–, and it is modulated by the cytosolic Cl^– concentration ([Cl^–]_i), extracellular bicarbonate, protein kinases, and the cytoskeleton state (21). The mathematical formulation of I_NaKCl (Eqs. A25 and A26) is based on the work by Strieter et al. (71), which described cotransport fluxes using nonequilibrium thermodynamics. This NaKCl description requires the specification of only one parameter, namely, the cotransport coefficient (L) (71), which was adjusted to the present model (Model Parameter Estimation). Model I_NaKCl follows an inherent 1:1:2 stoichiometry for Na⁺:K⁺:Cl^–, respectively, as shown experimentally.

PMCA current. The function of the PMCA is to extrude Ca²⁺ to the extracellular medium (77). In both ECs and rabbit atrial cells, the PMCA is considered a low-capacity high-affinity transport system (47, 66). In CPAECs, the PMCA function was experimentally observed to be high affinity and [Ca²⁺]_i, time, and CaM regulated as well as linked to NCX activity via [Ca²⁺]_i, comprising an intricate dynamic compensatory mechanism (68). The PMCA is responsible for determining resting [Ca²⁺]_i levels and counteracting moderate perturbations in [Ca²⁺]_i, having a higher affinity for Ca²⁺than NCX (68). The I_PMCA formulation (Eq. A27) is based on previous descriptions (47, 83). The average Hill coefficient and [Ca²⁺]_i for half activation of PMCA current values were calculated from cultured CPAE cell data (68). _PMCA was estimated by rescaling the value reported for a rabbit atrial cell model (47).

Fluid Compartment Model

In this EC model component, the balance of intracellular chemical species takes place (Eqs. A33–A41). The model assumes that Na⁺, K⁺, and Cl^– can occupy the total intracellular volume while Ca²⁺ is restricted to discrete cellular compartmental volumes, namely, the cytosolic volume available to free Ca²⁺ and the IP₃-sensitive store volume. Another major model assumption is that chemical species are homogeneously distributed throughout the intracellular environment (67).

Intracellular Ca²⁺ store. The model contains an intracellular Ca²⁺store mainly composed by the ER (74). It is sensitive to IP₃ and considered the principal Ca²⁺-release site in ECs (55) via IP₃R activation, releasing Ca²⁺ from the store into the cytosol once IP₃ binds IP₃R (74). Additionally, EC experimental data have suggested that ryanodine-sensitive stores coexist with IP₃-sensitive ones and that the former plays an essential role in shaping EC Ca²⁺ signals (36, 60). However, the ryanodine-sensitive stores are not implemented in the model as they lack a quantitative description level in ECs.

IP₃ generation and IP₃R current. Upon agonist binding to cell surface receptors, a G protein activates PLC, which in turn hydrolyzes phosphatidylinositol (4,5)-bisphosphate to give IP₃ and diacylglycerol (not shown) (Fig. 1B). IP₃ dynamics are described in the model (Eqs. A40 and A41) like in Ref. 43. Thus, agonist stimulation in the present model is simulated by a step increase in the steady-state IP₃ generation rate constant, which causes an exponential increase in IP₃ generation (Eq. A41) with time constant _IP3, and the slower IP₃ breakdown occurs via a first-order reaction with rate constant k_DIP3 (Eq. A40). This approach allows the model to be easily adapted to simulate EC responses to different types of agonists (e.g., thrombin or ACh). The I_IP3R formulation selected here (Eqs. A29 and A30) is based on two receptor models found in Refs. 28 and 80. These IP₃R models were combined and modified based on experimental data from porcine aortic ECs (11). The data suggest that IP₃-mediated receptor activation follows a Hill function with a coefficient of 3.8, suggesting approximately four binding sites for IP₃, which reflects the channel's tetrameric nature (74). The data gathered from porcine ECs also underscore an inhibitory effect of [Ca²⁺]_i on I_IP3R. Ca²⁺ inhibits IP₃R according to a modified version of the inactivation function (P_i,IP3R) from the IP₃R model in Ref. 28. The P_i,IP3R function was fitted to porcine EC data (11) (Fig. 4). Although cytosolic Ca²⁺ can also activate IP₃R (5) and this effect is included in other IP₃R models (28, 80), such activation is absent in vascular ECs (11) and thus not considered here.

View larger version (12K): [in this window] [in a new window]   Fig. 4. IP3R inactivation probability (Pi,IP3R) model fit (solid line) as a function of [Ca2+]i to experimental data (circles) from porcine aortic ECs (11). The intracellular IP3 level used experimentally, as well as in the model fit, was 1.4 µM.

Ca²⁺ buffering. Both cytosolic and ER Ca²⁺ buffering are included in the present model, unlike in previous EC modeling efforts. The rate of [Ca²⁺]_i buffering by CaM (Eq. A35) was based on the HUVEC model (77). The buffering mechanism inside the IP₃-sensitive store is provided by CSQN kinetics, which is described mathematically by the rapid buffering approximation (Eq. A36) as applied in Ref. 79. Mitochondrial buffering was excluded from the model since this organelle does not modulate Ca²⁺ signals in cultured CPAECs (42).

Sensitivity Analysis

Sensitivity analysis of the EC model was performed utilizing the Latin hypercube sampling (LHS) method and multiple regression techniques as previously described (8, 19). Parameters with significant uncertainty and/or contribution to model's behavior were chosen for the present analysis (Table 3). A variation range of 20% was assigned to parameters obtained from rat mesenteric ECs. Uncertainties of 40% and 60% were assigned to parameters obtained from EC types other than rat mesentery and cell types different from ECs, respectively. The highest uncertainty range (80%) was applied to parameters whose experimental method of determination did not allow an absolute measurement of the parameter. Although extracellular ionic concentrations can be accurately determined, uncertainty ranges of 20% were assigned to them to simulate variation in experimental media preparations in different studies. Fifty simulations were performed by selecting parameter values randomly and without replacement. Sensitivity of each output to selected model parameters is quantified by the sensitivity index (_i,j) as calculated from linear least-squares multiple regression of Eq. 2:

(2)

where Y_i^k is the value of the ith output at the kth run and X_j^k is the change in percentage of the central value of the jth parameter. Confidence intervals for the sensitivity indexes were calculated from an "elliptically shaped" joint 100(1 – )% confidence region for all the sensitivity indexes to determine whether they were statistically different from zero using the F-statistic (26). Five key simulation outputs were selected for analysis based on their impact on EC physiology: 1) resting V_m (V_m,rest), 2) resting [Ca²⁺]_i levels ([Ca²⁺]_i,rest), 3) maximum V_m hyperpolarization (V_m,max), 4) peak [Ca²⁺]_i level achieved ([Ca²⁺]_i,peak), and 5) [Ca²⁺]_i level at the end of stimulation ([Ca²⁺]_i,end).

All the parameters used in model simulations along with their reference sources are defined and shown in Table 1. A C_m of 14 pF was assumed for the model as this value falls within the range of capacitances (10–25 pF) reported in the literature for cultured vascular ECs (1). A capacitance to surface area ratio of 1 µF/cm² was assumed for the estimation of EC plasma membrane surface area. Although parameter optimization was not performed to fit integrated cellular responses, a few parameters had to be modified from their initial (literature derived) values as explained below. Adjustments in parameters such as channel conductances and permeabilities are justified when values are based on experimental data from cells other than rat mesenteric ECs, given that channel expression can significantly vary among different cell types or even in the same cell type at different cell cycle stages (56).The parameters describing [Ca²⁺]_i-dependent activation of K_Ca channels were selected from Refs. 2 and 81 to render them dormant at rest, as reported in rat mesenteric ECs (16, 54). Moreover, both channel conductances were increased, aiming to achieve physiologically comparable agonist-induced rat mesenteric EC hyperpolarization (54). The CaCC conductance was decreased 10-fold to allow significant agonist-induced hyperpolarization, as normally observed in rat mesenteric arteries. The [Ca²⁺]_IS necessary for about half-maximal SOC activation was assumed to be 50% of the total content of the model IP₃ store, as reported experimentally (69). I_SOC parameters (Eq. A15) were modified from their initial values to achieve a better agreement between simulations and experimentally recorded [Ca²⁺]_i transients in rat mesenteric ECs (Figs. 2C and 6C in Ref. 54). P_NSC,K, _NaK, and L_NaKCl were also adjusted to achieve physiological intracellular ionic concentrations and V_m at rest. _IP3 was reduced because the ACh-mediated [Ca²⁺]_i transient in rat mesenteric ECs reached peak 40 s earlier than the reference HUVEC responses, which were caused by thrombin (77). The maximum current via IP₃R was changed due to modifications in the original IP₃R mathematical formulations (28, 80).

Numerical Methods

The model contains a total of 10 nonlinear differential equations, which were coded in FORTRAN 90 and solved numerically using Gear's backward differentiation formula method for stiff differential equation systems (IMSL Numerical Library routine). The maximum time step was 4 ms, and tolerance for convergence was 0.0005.

	RESULTS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX: MODEL EQUATIONS GRANTS REFERENCES

 
Resting Conditions

EC model initial conditions (Table 2) were obtained by letting the system reach equilibrium using the parameters shown in Table 1 according to adjustments discussed in Model Parameter Estimation. The simulated EC had a resting V_m of approximately –20 mV, which is close to the value reported (–19.2 mV) for rat mesenteric ECs in Ref. 54. Moreover, the membrane resistance, calculated by injecting a known transmembrane current (I_stim = ±1 pA) and observing the magnitude of V_m change, was close to 2 G, which is in the range (1–5 G) reported for ECs (56). Blockade of NaK under resting conditions caused a 4 mV depolarization in V_m, which is close to the range (5–10 mV) recorded from ECs after [K⁺]_o-free solution or ouabain-induced NaK inhibition protocols (1). Resting model [Ca²⁺]_IS was physiological (2.25 mM) (74).

Resting V_m Modulation

Figure 5 presents model responses for different inhibition levels of VRAC. Blockade of VRAC was simulated by multiplying the G_VRAC parameter by 0.0, 0.25, and 0.50 for 100%, 75%, and 50% levels of inhibition, respectively. Inhibition of VRAC was performed for 120 s starting at three different time points, namely, time = 100, 320, and 680 s to mimic the experimental protocol in Fig. 11D from Ref. 55. As Fig. 5 shows, the model was able to replicate the experimentally observed EC V_m hyperpolarization response when a compound such as NPPB or mibefradil blocked G_VRAC (Fig. 11D in Ref. 55). Model outputs rapidly reached a stable V_m condition, and the magnitude of the hyperpolarization for 100% inhibition of VRAC was similar to the experiments (Fig. 11D in Ref. 55).

View larger version (24K): [in this window] [in a new window]   Fig. 5. Effect of blockade of VRAC on resting model Vm dynamics. Simulated VRAC blocking protocols for 50%, 75%, and 100% reductions in VRAC conductance (GVRAC) at 3 distinct periods for 120 s each are shown.

View larger version (17K): [in this window] [in a new window]   Fig. 6. Effects of raising extracellular K+ on resting model Vm. A: model Vm dynamics in different simulated experimental protocols, namely, control, blocked (100%) and doubled (2x) Kir conductance (GKir), and VRAC blockade (100%). B: simulation results showing steady-state Vm values at various extracellular K+ levels in control, Kir, and VRAC inhibitory conditions (100%).

View larger version (24K): [in this window] [in a new window]   Fig. 8. Model agonist-induced [Ca2+]i (A), total Ca2+ influx currents (B), Vm (C), and Ca2+ entry driving force (D) changes before and after blocking the effectiveness of KCa channels by either direct conductance block (SKCa and IKCa conductances set to 0) or increasing extracellular K+ to 35 mM.

Model Responses to Agonist Stimulation

Figure 7 depicts the overall dynamics of the main components of the EC model, namely, V_m, ionic currents, and intracellular species concentrations, upon short-term agonist stimulation. The agonist-induced EC response was simulated by changing the steady-state IP₃ generation rate constant from 0 to 5.5 x 10^–8 mM/ms at time = 20 s and maintaining for 100 s. Thus, as IP₃ slowly approached a constant value (Fig. 7D), [Ca²⁺]_i levels increased (Fig. 7D) prior to the activation of K_Ca channels (Fig. 7H) and before V_m began to hyperpolarize (Fig. 7A). As [Ca²⁺]_i approached its peak, I_NCX and I_PMCA (Fig. 7C) and I_SERCA (Fig. 7F) achieved maximum activity, which decayed as cytosolic Ca²⁺ decreased, thus reflecting their Ca²⁺ dependence. The depolarizing I_SOC (Fig. 7B) was enhanced during the response as V_m hyperpolarized and to a small extent by IP₃-sensitive store depletion, with I_SOC,Ca magnitude being larger than I_SOC,Na at rest and during stimulation. I_NSCs (Fig. 7E) also exerted a depolarizing effect on V_m as I_NSC,K decreased, and both I_NSC,Ca and I_NSC,Na became more negative, following changes in electrochemical gradients. The largest hyperpolarizing currents came from the activity of K_Ca and Kir channels (Fig. 7H). I_CaCC and I_VRAC (Fig. 7I), on the other hand, were mostly responsible for V_m depolarization (more so than NSC and SOC). Given the activity of ionic channels, [K⁺]_i and [Cl^–]_i levels decreased while [Na⁺]_i increased during agonist stimulation (Fig. 7G). NaK activity (Fig. 7C) was minimal as V_m reached maximum hyperpolarization, which hampered normal NaK function by hindering Na⁺ efflux, but increased following the rise in [Na⁺]_i and V_m repolarization. The major intracellular Ca²⁺ current came from IP₃R release of Ca²⁺ (Fig. 7F), which peaked around the same time as [Ca²⁺]_i. The Cl^– component of the cotransporter (I_NaKCl,Cl) slightly increased during stimulation (Fig. 7I). Moreover, as [Ca²⁺]_i reached peak, [Ca²⁺]_IS continued to decrease, suggesting a dynamic uncoupling of intracellular Ca²⁺ in this phase of the transient (Fig. 7, D and G).

Figure 8, A and C, illustrates long-term model Ca²⁺ and V_m dynamics, respectively, under agonist stimulation (at time = 20 s) as performed in Fig. 7. Figure 8 also shows how simulated transmembrane Ca²⁺ influx (via NSC and SOC; Fig. 8B) and its driving force (Fig. 8D) change during long-term stimulation. The model predicted that [Ca²⁺]_i changes are followed by variations in V_m (Fig. 8, A and C), and both responses displayed a similar transient pattern, as reported in experiments (55, 72). These simulations replicate characteristic Ca²⁺ and V_m responses experimentally recorded from rat mesenteric ECs under ACh stimulation (Figs. 2D and 6C in Ref. 54) at least for the first 100 s of stimulation reported in the study. Blockade of V_m hyperpolarization, by either direct K_Ca inhibition (i.e., SK_Ca conductance and IK_Ca conductance set to 0) or raised [K⁺]_o from 5 to 35 mM, caused a small reduction in [Ca²⁺]_i during the transient but did not affect the shape and dynamics of the agonist-induced [Ca²⁺]_i signal (Fig. 8A), as found experimentally (Fig. 6C in Ref. 54). Once K_Ca channels were inhibited, model V_m was unable to hyperpolarize, remaining essentially constant (Fig. 8C), and thus [Ca²⁺]_i alterations could no longer affect V_m significantly (Fig. 8, A and C). High [K⁺]_o was not as effective in reducing the [Ca²⁺]_i transient as direct K_Ca blockade, which agrees with experimental data (Fig. 6C in Ref. 54). V_m hyperpolarization increased [Ca²⁺]_i. The relative effect of V_m, however, was not homogeneous throughout the transient, becoming more significant as the later part of the transient approached (i.e., [Ca²⁺]_i increased 24% and 64% at peak and plateau, respectively, from hyperpolarization block to control conditions). Control transmembrane Ca²⁺ currents comprised 45% of the total Ca²⁺ flux (mainly NSC, SOC, and net store release) at 100 s after the initiation of the stimulation but comprised 87% of the total at 1,000 s. Agonist-stimulated Ca²⁺ influx currents increased more than twofold from hyperpolarization block to control conditions (Fig. 8B), although the maximal Ca²⁺ entry driving force change was only 35 mV (Fig. 8D), which was just 20% of the total resting driving force. Activation of SOC by store depletion became significant in the latter part of the Ca²⁺ transient (Fig. 8B).

To better examine the impact of V_m on [Ca²⁺]_i, voltage-clamp simulations were performed under resting and agonist stimulation conditions. Figure 9 depicts [Ca²⁺]_i as a function of clamped V_m (i.e., the V_m was set to a given value and [Ca²⁺]_i was estimated at rest and following stimulation). The resting, peak, and plateau [Ca²⁺]_i in Fig. 9 were normalized with respect to their corresponding values after being clamped at –20 mV. (Note that –20 mV is the V_m under control conditions or during agonist stimulation with K_Ca channels blocked.) Plateau and resting Ca²⁺ levels were more sensitive to V_m (and consequently to transmembrane Ca²⁺ influx) than peak Ca²⁺ levels, with resting levels a little more sensitive than plateau levels.

View larger version (17K): [in this window] [in a new window]   Fig. 9. Model Vm sensitivity of resting, peak, and plateau EC [Ca2+]i levels. To determine the sensitivity of resting [Ca2+]i, the cell was clamped at –20 mV followed by a corresponding voltage step for 1,000 s to reach a new steady state. To determine the sensitivity of peak and plateau [Ca2+]i, the voltage step was applied simultaneously with agonist stimulation (steady-state IP3 generation rate constant from 0 to 5.5 x 10–8 mM/ms). Data were normalized by the rest, peak, and plateau [Ca2+]i at –20 mV clamp.

The outputs of all 50 simulations as well as the model control response are shown in Fig. 10. Only one of the simulations, indicated by the arrow, exhibited extraneous behavior, and it was discarded from the further analysis as an outlier. The remaining 49 simulations were used for the application of the LHS method. Table 4 presents semirelative sensitivity coefficients calculated according to Eq. 2 for the five selected outputs. V_m,rest and V_m,max were significantly sensitive to 4 and 6 parameters, respectively, whereas [Ca²⁺]_i,rest, [Ca²⁺]_i,peak, and [Ca²⁺]_i,end were significantly affected by 2, 10, and 6 of the parameters, respectively. _PMCA was the parameter that had the most significant impact on model predictions, and [Ca²⁺]_i,peak was the output with the largest sensitivity on model parameters.

View larger version (44K): [in this window] [in a new window]   Fig. 10. EC model control responses and outputs for 50 simulations of agonist stimulation performed for sensitivity analysis using the Latin hypercube sampling method. A: intracellular Ca2+ transient profiles. B: Vm dynamic outputs. The outlier curve (arrow) was not considered for statistical analysis.

	DISCUSSION

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX: MODEL EQUATIONS GRANTS REFERENCES

Model Verification

The model's performance was assessed by its ability to reproduce established features of rat mesenteric EC behavior while utilizing physiological parameter values. The model's resting values for V_m, for the four intracellular ionic concentrations, and for store Ca²⁺ agree with experimental data. Whole cell conductance was also within the physiological range, and this provides an overall validation for the descriptions of transmembrane currents. V_m responses to NaK blockade, to VRAC blockade (Fig. 5), extracellular K⁺ challenge (Fig. 6), and agonist stimulation (Fig. 7A) also validate plasma membrane electrophysiology. Intracellular Ca²⁺-handling machinery was tested against experimental data with agonist-induced stimulation (Fig. 8A). Most importantly, proper model behavior in these scenarios was accomplished by adjusting (within physiological ranges) a limited number of unknown parameters (i.e., SK_Ca conductance, IK_Ca conductance, CaCC conductance, K_SOC,CaIS, _NaK, P_NSC,K, L_NaKCl, _IP3, and _IP3R).

Resting V_m Modulation

Experimental data from bovine pulmonary artery ECs (Fig. 11D in Ref. 55) have shown that VRAC is crucial for establishing resting EC V_m. The three levels of simulated VRAC blockade (Fig. 5) can be interpreted as different levels or types of inhibition (i.e., distinct blocker concentrations or degrees of hyperosmotic stress) or dissimilar levels of VRAC expression in ECs. Model simulations concur with VRAC inhibition experiments (25, 55). The model also supports the notion that different levels of VRAC and Kir channel expression or activation can lead to EC switching from K type (resting V_m close to E_K) to Cl type (resting V_m near E_Cl), or vice versa, as proposed to occur throughout the cell cycle (56) or using channel inhibition protocols (Fig. 7b in Ref. 25). New experimental studies are required to verify model predictions (Fig. 5) and confirm basal VRAC activity in rat mesenteric ECs as found in other ECs (55, 72).

The EC model was used to investigate two related hypotheses: one regarding the ability of ECs to hyperpolarize in response to raised [K⁺]_o under certain conditions (25) and another suggesting that K⁺ is the EDHF (22, 45) (Fig. 6). The former hypothesis was formulated based on rat mesenteric artery experiments showing that these vessels could relax to increased [K⁺]_o (from 5.88 to 10.58 mM) only if the VRAC in ECs was blocked. To explain these results, it was postulated that ECs in intact vessels can only hyperpolarize to raised [K⁺]_o if they are switched from Cl- to K-type cells by VRAC inhibition, using either channel antagonists or hyperosmotic media. However, isolated EC experiments (Fig. 8, C and D, in Ref. 55) have demonstrated that when EC V_m is close to E_Cl (Cl type), the cell responds to raised [K⁺]_o (from 6 to 12 mM) by hyperpolarizing. This V_m behavior was captured by the model (Fig. 6A) and is attributed to the crossover effect of the Kir I-V relationship. Figure 6A also shows that a K-type EC (i.e., after VRAC block) depolarizes upon [K⁺]_o challenge (from 5 to 12 mM). The model V_m (Fig. 6A, solid and dotted lines) achieved a new resting value (during [K⁺]_o stimulation) faster than what was observed experimentally (Fig. 8C in Ref. 55). This may attributed to a slower increase of [K⁺]_o in the experiments relative to the simulations. Despite these minor differences, there is an overall agreement between simulated EC V_m behavior and experiments in isolated cells.

The model suggests that hyperpolarization and depolarization can occur in response to K⁺ increase, under both control and VRAC block conditions, depending on the initial and final K⁺ concentration (Fig. 6B). Although the model predicts depolarization after VRAC blockade for the particular K⁺ concentration range examined in Ref. 25, simulations do not directly translate to EC behavior in intact vessels as coupling to SMC may alter EC electrophysiology. Interestingly, Fig. 6B shows that a large V_m hyperpolarization is possible also under VRAC block, when the initial value of [K⁺]_o is low (36 mV, from 1 to 3 mM [K⁺]_o). Under control or Kir block conditions (i.e., solid and dotted lines in Fig. 6B), resting V_m at low starting [K⁺]_o (e.g., 0.5 mM) is rather hyperpolarized (i.e., about –40 mV). This is attributed to a very negative E_K. As the [K⁺]_o increases, the increase in E_K drives V_m to more depolarized values. As [K⁺]_o increases between 5 and 8 mM, V_m hyperpolarizes, when Kir are active, due to the crossover effect. When VRAC is blocked (i.e., dashed-dotted line in Fig. 6B) and at very low [K⁺]_o, the model predicts a significant depolarization as a result of a decrease in I_NaKCl that leads to the establishment of a larger transmembrane Na⁺ gradient and an increase in I_NSC,Na.

Assuming that a similar behavior as predicted in Fig. 6B occurs in intact vessels, the model suggests that by inhibiting VRAC, the EC V_m hyperpolarization response to [K⁺]_o gets larger (for a particular [K⁺]_o range), which may explain why rat mesenteric artery relaxations are only observed after this channel is blocked (25). Additionally, experimental setup differences could create different levels of VRAC activation and thus dramatically affect the predicted profiles (Fig. 6B). This may explain in part the variability observed in the response of rat small mesenteric arteries to [K⁺]_o (25). This issue merits further experimental and theoretical (via coupled EC-SMC models) investigations. Larger hyperpolarizing responses in the model might be also possible with higher values for resting [K⁺]_i.

Regarding the K⁺ as EDHF hypothesis, model simulations (Fig. 6) support the idea that Kir channels on ECs may play an important role in the relaxation of rat mesenteric arteries upon [K⁺]_o stimulation and thus in part of the EDHF response (23, 45). As SMCs and ECs in rat mesenteric arteries are electrically coupled via myoendothelial gap junctions, EC V_m hyperpolarization caused by [K⁺]_o-induced Kir activation is then able to be transmitted to the adjacent SMCs and thus relax the vessel (23, 25, 45, 65). However, simulations are based on an isolated EC model and may not correspond to the in vivo EC-SMC coupling situation, which might shift model response profiles (Fig. 6B) and may explain why maximum artery relaxations (23, 45) are observed at [K⁺]_o levels higher than predicted for maximum EC V_m hyperpolarization. Additionally, elevated [K⁺]_o may activate NaK in SMCs, causing further SMC V_m hyperpolarization and thus affect relaxation profiles (22). Therefore, theoretical studies using EC-SMC integrated models are necessary to further analyze the role of EC Kir channels and K⁺ in the EDHF response.

Model Responses to Agonist Stimulation

Following the examination of EC plasma membrane electrophysiological behavior, the next task was to investigate if the model was able to replicate EC responses to agonist stimulation. The overall dynamic agonist-induced responses of the main model components are displayed in Fig. 7. All EC model components performed their expected roles in the agonist-induced response, in keeping with literature (1, 55, 61). Model behavior can be explained in the following way. IP₃ activates IP₃R, releasing store Ca²⁺ into the cytosol, which in turn activates K_Ca channels and thus causes V_m hyperpolarization, increasing the driving force for Ca²⁺ entry. A significant increase in I_Kir during hyperpolarization (a result of the negative slope of the Kir I-V relationship; Fig. 2A, dashed line) acts as a positive feedback to hyperpolarization. NCX, PMCA, and SERCA work to counteract this sudden elevation in [Ca²⁺]_i, shaping the Ca²⁺ transient peak. As [Ca²⁺]_i reaches a peak, [Ca²⁺]_IS continues to decrease during agonist stimulation, and such intracellular Ca²⁺ dynamic uncoupling has been observed in CHO cells stimulated with the IP₃-generating agonist ATP (7). As IP₃-sensitive stores are emptied and V_m hyperpolarizes, I_SOCs are activated and contribute to the sustained [Ca²⁺]_i rise. I_NSCs, on the other hand, tend to offset the hyperpolarizing effect of K_Ca activation by passively following their electrochemical gradients. The overall depolarizing effect on V_m exerted by NSC during agonist stimulation, as I_NSC,K decreases and both I_NSC,Ca and I_NSC,Na become more negative (Fig. 7E), has been observed experimentally (56). Moreover, the [Ca²⁺]_i increase and V_m hyperpolarization activate CaCCs, while I_VRACs increase by V_m hyperpolarization alone (Fig. 7I). These two Cl^– conductances are the main players in the V_m repolarization process, which decreases the Ca²⁺ influx driving force. As a result of the agonist-induced processes described above, the levels of [K⁺]_i, [Cl^–]_i, and [Na⁺]_i are altered (Fig. 7G). During stimulation, NCX, NaK, and NaKCl activities tend to compensate intracellular ionic species changes and thus restore cellular homeostasis, thus reflecting the importance of these model components.

Figure 8 illustrates that the model reproduces dynamic [Ca²⁺]_i and V_m signals recorded from ACh-stimulated rat mesenteric ECs under both control and K_Ca blockade (54)