- A mathematical model of the cytosolic-free calcium response in endothelial cells to fluid shear stress

Intracellular transients in Ca2+ are a typical cellular response to many different stimuli, including shear stress (1). Flow enhances the delivery of blood-borne agonists to the cell surface (2). These agonists in turn mediate the release of calcium from intracellular stores into the cytosol. Shear stress increases the permeability of the cell membrane to extracellular calcium, thereby increasing the cytosolic concentration (3, 4). A third transduction mechanism is a direct shear stress effect upon the generation of inositol trisphosphate, which could mediate the release of calcium from intracellular stores (5). This last would operate via a putative mechanoreceptor behaving as a biological "strain gauge." It is likely that all three pathways act in concert to elicit the calcium transient.

The response of cytosolic calcium to a fluid shear stress has been illuminated to a large degree by experimental work. It is therefore natural to propose mathematical models that account for this behavior. Several models have been proposed that relate cytosolic-free calcium dynamics to stimulation by agonists (6, 7). The concentration profiles of agonists in flow chambers containing cell monolayers have been described (8). There has been at least one report of a mathematical model linking calcium transients to mechanical stimulation (9). There have been no reports of a model encompassing the signal recognition and transduction cascade from flow to cytosolic calcium.

This work presents a model that integrates previous models on agonist stimulation with a new model of shear stress-facilitated calcium ion transport into cells. It builds upon the earlier lumped parameter model of intracellular calcium dynamics by Wiesner et al. (6) by adding flow effects. In this earlier model, a balance on cytosolic-free calcium ion consists of four fluxes.

<FR><NU>dCa<SUB>c</SUB></NU><DE>dt</DE></FR>=q<SUB>s</SUB>&plus;q<SUB>in</SUB>−q<SUB>b</SUB>−q<SUB>out</SUB>. [ 1 ]

In this expression, Cac is the cytosolic-free calcium concentration, qs is the net flux released from intracellular stores into the cytosol in response to stimulation by agonists, qin is the rate of influx of extracellular calcium into the cytosol, qb is the rate of buffering of calcium ion by soluble cytosolic proteins, and qout is the rate of extrusion and exchange of calcium to the extracellular space. The complete calcium dynamics model, with its parameter values and detailed equations, is given in ref. 6.

An illustration of the flow effects in a capillary tube flow chamber is given in Fig. 1. A confluent monolayer of endothelial cells has been cultured on the bottom surface. The chamber is perfused with media in steady laminar flow with known concentrations of the agonist thrombin and extracellular calcium. Agonist is transported to the cell surface by coupled convection and diffusion. Flow of the perfusate exerts a viscous shear force upon the monolayer. The effect of flow upon mobilization of calcium from intracellular stores is described through its influence on the surface concentration of agonist, cs. This is an effect of shear rate. Shear stress itself affects extracellular calcium influx, qin. This is described via a model of mechanosensitive ion channels. Direct mechanical effects triggering the release of calcium from intracellular stores are not included. 

Mass Transfer of Agonist in the Capillary Flow Chamber

The coupled convection and diffusion of agonist to the endothelium may be described by a species balance on thrombin in the perfusate.

<FR><NU>∂c</NU><DE>∂t</DE></FR>&plus;v <FR><NU>∂c</NU><DE>∂z</DE></FR>=<UP>D</UP><FENCE><FR><NU>∂<SUP>2</SUP>c</NU><DE>∂x<SUP>2</SUP></DE></FR>&plus;<FR><NU>∂<SUP>2</SUP>c</NU><DE>∂y<SUP>2</SUP></DE></FR></FENCE>. [ 2 ]

In this equation, v is the z-component of velocity, t is time, c is the local concentration of agonist, and D is the diffusivity.

The species balance (Eq. 2) is subject to the conditions that no agonist is present initially (c[x, y, z, 0] = 0), and that the inlet agonist concentration is fixed (c[x, y, 0, t] = co). Two other boundary conditions are required. There is no agonist consumption at cell free surfaces.

<FR><NU>∂c</NU><DE>∂y</DE></FR><FENCE><SUB><SUB>y<UP>=&plus;</UP>b</SUB></SUB>=<FR><NU>∂c</NU><DE>∂x</DE></FR></FENCE><SUB><SUB>x<UP>=&plus;</UP>a</SUB></SUB>=0. [ 3 ]

In the experiment, the cells cover only the bottom surface at y = -b. At this surface, the agonist binds to cell surface receptors. The rate of binding depends upon the surface concentration of agonist. Therefore, the consumption of agonist at y = -b is

S<UP>D</UP> <FR><NU>∂c</NU><DE>∂y</DE></FR><FENCE><SUB><SUB>y<UP>=−</UP>b</SUB></SUB>=<FR><NU>k<SUB>on</SUB>R<SUB>u</SUB>c<SUB>s</SUB>−k<SUB>off</SUB>R<SUB>b</SUB></NU><DE>N<SUB>A</SUB></DE></FR></FENCE>. [ 4 ]

Here Rb is the number of bound receptors, Ru is the number of unbound receptors, and cs is the concentration of agonist at the cell surface. kon, koff, and kd are on, off, and degradation rate constants, respectively. NA is Avogadro's number, and the quantity S is the en face area of the cell.

The steady-state velocity profile in a rectangular duct undergoing a step increase in pressure gradient is established almost immediately (10). Therefore, an analytical steady-state expression for the profile is used (11). Differentiation of this profile yields the shear stress as a function of position and perfusion rate.

 &tgr;(X, Q)= [ 5 ]
&tgr;(X, Q)= 

In this equation, the dimensionless coordinates X = x/a and Y = y/b have been introduced. The quantity gamma = b/a is the aspect ratio of the channel and Q is the perfusion rate. tau is the shear stress and eta is viscosity of the perfusate.

Using the bulk agonist concentration, c0, shear stress (tau), and flow chamber geometry as inputs, the species continuity balance is solved to obtain the surface concentration of agonist, cs, which in turn drives the activation of cell surface receptors.

<FR><NU>dR<SUB>b</SUB></NU><DE>dt</DE></FR>=k<SUB>on</SUB>R<SUB>u</SUB>c<SUB>s</SUB>−(k<SUB>off</SUB>&plus;k<SUB>d</SUB>)R<SUB>b</SUB>. [ 6 ]

The number of bound receptors governs the generation of inositol trisphosphate, which in turn drives the release of calcium from intracellular stores, qs. Thus, the calcium transient is linked dynamically to the applied shear stress, tau.

Influx Pathway

Shear stress-activated calcium channels have been identified in human umbilical vein endothelial cells (HUVECs) (3) and distal renal tubule cells (4). Receptor-operated channels are also expressed in endothelial cells, and are opened by the binding of extracellular ligands such as thrombin. The increase in calcium influx due to receptor-operated channels (2.53 µM/sec, ref. 12) is quantitatively smaller than that due to shear stress (138 µM/sec, ref. 3). Therefore, the transduction mechanism between agonists and calcium influx is not included in this model. Calcium-activated calcium channels and voltage-gated ion channels are not expressed in cultured endothelial cells (13, 14).

Passive influx of calcium ions has also been established in endothelial cells (12). Shear stress augments this passive leak by inducing local transient defects in the structure of the plasma membrane. The increase in passive calcium transport as a function of shear stress has been measured in liposomes (15). The increase at a shear stress of 3.6 N/m2 is five orders of magnitude less than the basal passive leak of calcium into endothelial cells. Therefore, increased influx of calcium due to shear stress-induced defects in membrane structure is neglected in this study.

The steady-state flux of an ion across the membrane is described by the Goldman-Hodgkin-Katz equation (16), which depends in part upon the membrane permeability coefficient P. In this model, the value of the permeability coefficient is reflected by the fraction of ion channels that are in the open state. From Kirchhoff's law for parallel conductors and the relationship between the membrane permeability coefficient and membrane conductance (16), it can be shown that the permeability is linearly related to the open channel fraction fo and the permeability of the membrane when all channels are open, Pmax: i.e., P = foPmax.

The application of shear stress induces a uniaxial tension gradient in the membrane and deforms it. Part of the applied load is borne by the membrane and part is borne by underlying cellular structures. The fraction borne by submembranous structures is denoted by varepsilon. The deformation of the membrane imparts a strain energy density that depends upon the applied shear stress W(tau). The fraction of open channels has a Boltzmann dependence upon the level of gating energy in the membrane (17). The gating energy in turn can be related to the strain energy density in the membrane. A fraction of the elastic strain energy stored within the membrane is considered to gate the shear stress-activated channel. This fraction is designated as fe. Thus, the fraction of channels in the open state dependsupon the strain energy density as indicated in the following equation:

f<SUB>o</SUB>(&tgr;)=<FR><NU>1</NU><DE>1&plus;&agr;∗ <UP>exp</UP> <FR><NU><UP>−</UP>f<SUB>e</SUB>W(&tgr;)</NU><DE>kTN</DE></FR></DE></FR>. [ 7 ]

Here N is area channel density, k is the Boltzmann constant, T is the absolute temperature, and alpha is a measure of the probability that a channel is in the open state in the no-load case.

The strain energy density in the membrane is now related to the applied shear stress. The equations of mechanical equilibrium are used in conjunction with a constitutive relationship to develop this dependence. Employing a simplifying analysis, Fung and Liu (18) demonstrated that the equations of mechanical equilibrium for the endothelial cell membrane reduce to:

<FR><NU>∂T<SUB>z</SUB></NU><DE>∂z</DE></FR>&plus;(1−&epsiv;)&tgr;=0 [ 8 ]


T<SUB>x</SUB>=0. [ 9 ]

The tensions in the direction of the flow and orthogonal to flow in the plane of the membrane are denoted by Tz and Tx, respectively. In Eq. 8, the tension gradient in the direction of flow is related to the load fraction borne by the membrane, (1 - varepsilon)tau. The quantity varepsilon is the fraction of the stress borne by the interior of the cell, which is assumed to have a solid-like content. The tension in the direction perpendicular to flow in the plane of the membrane is taken as zero by virtue of the tension field hypothesis. This hypothesis states that the membrane cannot support compressive stresses in its own plane because its thickness is small compared with its length and width. Considering the stresses acting in the plane of the endothelial membrane, Fung and Liu argue that the stress orthogonal to the direction of flow is small and, thus, the simplifying assumption of the stress in the membrane being a uniaxial tension field in the direction of flow is valid. This same assumption is used here.

A well-known constitutive relation for cell membranes in pure shear deformation under uniaxial tension is given by Evans and Skalak (19).

T<SUB>z</SUB>=&mgr;(&lgr;<SUB>z</SUB><SUP>2</SUP>−&lgr;<SUB>z</SUB><SUP><UP>−</UP>2</SUP>). [ 10 ]

The quantity lambdaz is the stretch ratio in the flow direction. The elastic constant µ is known as the membrane shear modulus.

Skalak et al. (20) give a strain energy density function, W, for a two-dimensional membrane. Integration of that expression and Eq. 8 yields the relationship between strain energy density and the applied shear stress.

W(&tgr;)=<FR><NU><FENCE>(1−&epsiv;)&tgr;L&plus;<RAD><RCD>16&mgr;<SUP>2</SUP>&plus;&tgr;<SUP>2</SUP>L<SUP>2</SUP>(&epsiv;<SUP>2</SUP>−2&epsiv;&plus;1)</RCD></RAD>−4&mgr;</FENCE></NU><DE>8<FENCE>(1−&epsiv;)&tgr;L&plus;<RAD><RCD>16&mgr;<SUP>2</SUP>&plus;&tgr;<SUP>2</SUP>L<SUP>2</SUP>(&epsiv;<SUP>2</SUP>−2&epsiv;&plus;1)</RCD></RAD></FENCE></DE></FR><SUP>2</SUP>, [ 11 ]

where L is the length of the cell in the flow direction.

The foregoing treatment assumes a steady state. In Schwarz et al. (3), the calcium current is not established immediately, but rather over a period of approx70 sec. The dynamics are associated with changes in membrane potential, kinetics of ion channel activation, and viscoelasticity of the membrane. The membrane potential hyperpolarizes, i.e., becomes more negative in response to shear stress (21). Membrane potential depends upon transmembrane currents of potassium, sodium, and chloride ions. For the sake of simplicity, an empirical relationship between shear stress and membrane potential, Deltaphi(t, tau), is used herein. The temporal aspect of membrane permeability, P(t, tau), is also described empirically. The expressions for Deltaphi(t, tau) and P(t, tau) are given in ref. 22.

Estimation of Parameters in the Model

The model derived in the preceding section describes the effects of flow upon the surface concentration of agonist, cs, and on the influx of extracellular calcium, qin. The constants required to calculate the effect of flow on cs are readily available from handbooks and the experimental setup. The effect of flow upon calcium influx is not as well characterized.

The influx of calcium as a function of the applied fluid shear stress has been measured at only one level of shear stress, 1.0 N/m2, in endothelial cells (3). In the absence of experimental data, a prediction of the magnitude of the flow-induced calcium influx is offered. The load fraction, varepsilon, is given a value of 0.90. This initial estimate is based upon the micropipette experiments of Sato et al. (23), in which the main load bearing component of the cell was found to be the cytoskeleton. With varepsilon fixed a priori, values for fe, alpha, and Pmax are obtained by satisfaction of the following three criteria: qin(tau = 0) = 5.33 µM/sec, qin(tau = 3) = 10 µM/sec, qin,max = 23 µM/sec.

The value of qin(tau = 0) is that for passive leak of calcium into the cytosol (12). The value of qin(tau = 3) was selected to obtain a value for the plateau level of calcium corresponding to an agonist-free experiment (24). The maximum value of calcium influx, qin,max, was chosen because influxes above this value removed the experimentally observed biphasic character of the calcium response to agonists when used in Eq. 1. All parameter values are tabulated in Table 1, along with their respective sources. The equations of the model were coded in FORTRAN and integrated numerically on a CRAY-EL computer.


The calcium influx, qin, is illustrated as a function of shear stress in Fig. 2. The calculation assumes constant intracellular calcium concentration, extracellular calcium concentration, and transmembrane potential. A sigmoidal relationship is seen. For a value of varepsilon = 0.90, there is a threshold value of approx0.5 N/m2 below which the influx is negligible. In the range of shear stress from approx0.5 N/m2 to approx6 N/m2, the flux is shear-sensitive. For shear stresses higher than the shear-sensitive range, the influx is maximal. This corresponds to the situation in which all ion channels are open. 

 The effect of varepsilon, the fraction of the load borne by the cytoplasm, is also shown in Fig. 2. The threshold value of shear stress and the shear-dependent range are altered. As more of the stress is borne by the membrane, the threshold level of influx decreases. Also, the shear-sensitive range narrows.

The response of cytosolic-free calcium is shown in Fig. 3. In this calculation, the perfusate is agonist free (cs = 0) and contains physiologic calcium (Caex = 1.5 mM). Note that no peak in the transient is observed when agonists are absent. The plateau level increases with increasing shear stress. This is due to enhanced transplasmalemmal calcium influx. Both the lack of a peak and the graded response to the level of shear stress in serum-free media was observed in bovine aortic endothelial cells (25). 

The agonist pathway is now examined by setting Caex to zero, and the bulk agonist concentration, co, to 1 unit per ml. When an endothelial cell monolayer is exposed to a step change in flow, a concentration profile is established over time by coupled convection and diffusion. The surface concentration of agonist reaches its steady-state value after more than 100 sec (simulation not shown).

In Fig. 4, the effects of mass transfer limitations and receptor dynamics are illustrated. The fastest curve shows the response of the calcium dynamics model to a step change in receptor activity. The peak of this curve is in excess of 18 mM, a nonphysiologic limiting case. The second curve shows the response of the calcium dynamics model coupled to receptor activation. The stimulus here is a step change in surface-agonist concentration. The slowest curve shows the response of the complete agonist model to a step change in shear stress to a level of 0.1 N/m2. This curve contains the effects of mass transfer limitations, receptor dynamics, and intracellular calcium dynamics. As more dynamic lags are included in the simulation, the time to peak increases, and the peak value is attenuated. It can be seen that receptor dynamics and the development of the concentration profile introduce significant delays. 

Given that mass transfer limitations can govern the response, it is of interest to vary the flow rate parametrically. In Fig. 5 one sees an acceleration in the response as the shear stress is increased from 0.01 N/m2 to 3 N/m2. This is due to enhanced mass transfer of agonist to the surface at higher flow rates. Noteworthy here is the fact that there is little difference between the responses at 0.5 N/m2 and 3 N/m2. This indicates that the agonist-mediated calcium mobilization process is controlled by the development of the concentration profile in the perfusate at low shears (2). On the other hand, the response is controlled by the intrinsic kinetics of receptor activation at shear stresses above 3 N/m2. This modulation of the agonist response by the level of shear stress has been observed experimentally (24). In contrast to the no-agonist simulation, the peak here is very prominent. Also, there is no elevated plateau in the absence of extracellular calcium. 

Synergistic effects of thrombin and physiologic extracellular calcium are shown in Fig. 6. The calcium transient is calculated for shear stresses from 0.1 N/m2 to 5 N/m2. Both the peak value and the plateau level increase at elevated shear stresses. The effects of shear stress become pronounced only at values higher than 1 N/m2. This is above the threshold level for which calcium influx becomes appreciable. At shear stresses below 1 N/m2, the primary modulation of the response is by mass transfer of agonist to the cell. This modulation is minor and becomes kinetically limited between 0.1 and 0.5 N/m2. It is only when calcium influx becomes large that shear stress begins to significantly alter the peak value and the plateau level. 
To assess the validity of the model, the simulation was compared with a benchmark experiment in HUVECs. Here, a HUVEC monolayer in a square capillary tube was exposed to a sudden onset of flow corresponding to a shear stress of 3 N/m2. The perfusate contained 1 unit per ml (approx9 nM) of alpha-thrombin and physiologic calcium. Calcium concentrations were measured with fluorescent microscopy. Details of the experimental protocol are given by Helmlinger et al. (24).

The model prediction is compared with the data in Fig. 6. The simulation for tau = 3 N/m2 reasonably predicts the calcium level for most of the time course, particularly in the peak region. The recovery phase is not reproduced as accurately as the peak region. The simulation goes to an elevated plateau. This is due to the sustained increase in qin. The plateau region is examined in more detail in the Discussion.


The agreement with the experiment was judged adequate for a first step in the modeling of the intracellular calcium transient, in that the model reproduces many experimental findings. A biphasic response is predicted. The transient peak is due to agonist-stimulated release from intracellular stores, whereas the sustained plateau is attributable to the influx of intracellular calcium. This prediction agrees with the conclusions of Hallam et al. (30). The model reproduces the experimental finding that flow and agonist together elicit a stronger calcium transient than either alone (24). This phenomenon is manifested not only in agonist experiments, but also in the use of serum-supplemented media, which may contain calcium-mobilizing agonists. The simulations clarify the relative contributions of the two pathways to the recognition and transduction of the fluid mechanical signal. This model predicts that shear stress modulates the calcium transient primarily by the influx pathway, and this modulation occurs at shear stress levels above 1 N/m2. Calcium is one of the most exquisite second messengers in a cell. Given the wide variety of flow patterns in the circulation, the influence of calcium-dependent processes can vary widely as well.

This modeling study also offers new insights into the individual mechanisms by which shear stress elicits a calcium transient. Receptor dynamics and development of concentration boundary layer both appear to be rate-limiting processes in comparison to intracellular calcium dynamics. Resistance to diffusion of agonist is the controlling step at shear stresses below 0.1 N/m2. The calcium transient moves from transport-limited to kinetically limited regime in the range of 0.1 N/m2 to 0.5 N/m2. The kinetic limitation here lies in the receptor dynamics.

The mechanosensitive ion channel model predicts a sigmoidal dependence of calcium influx upon shear stress. Many cellular responses to shear stress are sigmoidal in nature. Schwarz et al. (31) found a sigmoidal dependence of the amplitude of calcium transients upon shear stress. Tseng et al. (32) likewise noted that microtubule-associated protein kinase is stimulated in a sigmoidal fashion by fluid shear stress. The sigmoidal dependence of calcium influx arises from the Boltzmann dependence of channel open probability upon membrane strain energy density.

As a first approximation to estimating membrane stresses and membrane strain energy density, a simple uniaxial tension field employed by Fung and Liu (18) has been used. The key to this analysis is the assumption that the stress perpendicular to the flow direction in the plane of the membrane is zero. The endothelial surface, however, is not one which is flat but rather a topography of "hills" and "valleys." For such a surface it is quite possible that stresses in the endothelial membrane perpendicular to the direction of flow may not be negligible locally. Therefore, it is possible that the membrane is in fact subjected to a biaxial stress field. The open channel probability is then most properly correlated with a strain energy density based upon a two-dimensional stress field. This would be a worthwhile future enhancement to the model.

The magnitude of calcium influx is influenced by the load-bearing contribution of the cytoskeleton. This contribution is embodied in the parameter varepsilon, which has not been experimentally determined. It is hypothesized that many direct mechanical effects are initiated by deformation of the membrane. The deformation of the membrane depends upon its load-bearing capacity relative to submembranous components. This suggests that differences in cytoskeletal and plasmalemmal structure could account for heterogeneities in the calcium response and other responses among different cell types. Establishing the value of varepsilon would help validate this hypothesis of mechanotransduction.

The simulations raise the question of what is the magnitude of shear stress-induced calcium influx. Schwarz et al. (3) reported a calcium current of 80 picoamps in cells exposed to a shear stress of 1 N/m2 using voltage clamp techniques. This corresponds to a calcium influx of 138 µM/sec, assuming a cytosolic volume of 3 picoliters. Use of such a large influx in conjunction with the published capacities of the plasma membrane efflux mechanisms and cytosolic buffers completely removes the biphasic nature of the transient and yields a plateau level of calcium several orders of magnitude larger than seen in vitro. The current measured by a voltage-clamp procedure may be nonphysiologic, in that the pipette represents a large additional conductance of calcium in parallel to plasmalemmal efflux mechanisms. It is also possible that this current is composed of other ions in addition to calcium. Thus, the magnitude of shear stress-induced calcium influx is not well established and requires further experiments.

The model predicts an elevated plateau in response to shear stress. Helmlinger et al. (24) (whose data are used as the benchmark experiment in this work) indicate that the calcium level returns to the initial baseline when the cells are stimulated with both thrombin and shear stress. In the case of shear stress only, however, they found an elevated plateau in response to shear stress. It is unclear why the plateau level should vary between these two experiments. Several mechanisms that could explain a return to baseline are not included in the model. These include deactivation of thrombin receptors and calcium entry channels. This suggests that further experiments are needed to establish both the behavior and mechanism for the calcium level at asymptotically long times.

The model predicts the modulation of the calcium transient in the physiologically relevant range. Thus, the findings of this study may have profound implications in vivo. Transport of biologically active molecules and membrane permeability depend upon the hemodynamic flow pattern. Events that alter the flow pattern, such as clot lysis and stenosis formation, can thereby influence the initiation and progression of disease. The model can better approximate the in vivo circulation with some enhancements. These enhancements include a pulsatile velocity profile and the use of a geometry similar to a blood vessel.

This model is intended as a first step in the mathematical description of flow effects upon intracellular calcium. Because calcium as a second messenger regulates many intracellular activities, the model can be an important tool in understanding the bridge between shear stress and the multitude of cellular processes it influences.


Dagger   To whom reprint requests should be addressed.


We are grateful to Dr. Gabriel Helmlinger for the experiments in which the response of the calcium transients in HUVECs to thrombin was established. T.F.W. was supported by National Science Foundation Grant BCS-911761 and National Institutes of Health Grant HL52218. R.C.B. is an Established Investigator of the American Heart Association (AHA) and his efforts were also supported by a National Grant-in-Aid from the AHA.


HUVEC, human umbilical vein endothelial cell.

rating: 4.00 from 1 votes | updated on: 4 Dec 2007 | views: 9084 |

Rate article:

excellent! bad…