table of contents (calcium influx / agonist / mechanotransduction) 
Biology Articles » Biomathematics » A mathematical model of the cytosolicfree calcium response in endothelial cells to fluid shear stress » Article
Article

[ 1 ] 
In this expression, Ca_{c} is the cytosolicfree calcium concentration, q_{s} is the net flux released from intracellular stores into the cytosol in response to stimulation by agonists, q_{in} is the rate of influx of extracellular calcium into the cytosol, q_{b} is the rate of buffering of calcium ion by soluble cytosolic proteins, and q_{out} 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, c_{s}. This is an effect of shear rate. Shear stress itself affects extracellular calcium influx, q_{in}. This is described via a model of mechanosensitive ion channels. Direct mechanical effects triggering the release of calcium from intracellular stores are not included.
The coupled convection and diffusion of agonist to the endothelium may be described by a species balance on thrombin in the perfusate.
[ 2 ] 
In this equation, v is the zcomponent 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] = c_{o}). Two other boundary conditions are required. There is no agonist consumption at cell free surfaces.
[ 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
[ 4 ] 
Here R_{b} is the number of bound receptors, R_{u} is the number of unbound receptors, and c_{s} is the concentration of agonist at the cell surface. k_{on}, k_{off}, and k_{d} are on, off, and degradation rate constants, respectively. N_{A} is Avogadro's number, and the quantity S is the en face area of the cell.
The steadystate velocity profile in a rectangular duct undergoing a step increase in pressure gradient is established almost immediately (10). Therefore, an analytical steadystate expression for the profile is used (11). Differentiation of this profile yields the shear stress as a function of position and perfusion rate.
[ 5 ] 
In this equation, the dimensionless coordinates X = x/a and Y = y/b have been introduced. The quantity = b/a is the aspect ratio of the channel and Q is the perfusion rate. is the shear stress and is viscosity of the perfusate.
Using the bulk agonist concentration, c_{0}, shear stress (), and flow chamber geometry as inputs, the species continuity balance is solved to obtain the surface concentration of agonist, c_{s}, which in turn drives the activation of cell surface receptors.
[ 6 ] 
The number of bound receptors governs the generation of inositol trisphosphate, which in turn drives the release of calcium from intracellular stores, q_{s}. Thus, the calcium transient is linked dynamically to the applied shear stress, .
Influx Pathway
Shear stressactivated calcium channels have been identified in human umbilical vein endothelial cells (HUVECs) (3) and distal renal tubule cells (4). Receptoroperated 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 receptoroperated 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. Calciumactivated calcium channels and voltagegated 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/m^{2} is five orders of magnitude less than the basal passive leak of calcium into endothelial cells. Therefore, increased influx of calcium due to shear stressinduced defects in membrane structure is neglected in this study.
The steadystate flux of an ion across the membrane is described by the GoldmanHodgkinKatz 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 f_{o} and the permeability of the membrane when all channels are open, P_{max}: i.e., P = f_{o}P_{max}.
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 . The deformation of the membrane imparts a strain energy density that depends upon the applied shear stress W(). 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 stressactivated channel. This fraction is designated as f_{e}. Thus, the fraction of channels in the open state dependsupon the strain energy density as indicated in the following equation:
[ 7 ] 
Here N is area channel density, k is the Boltzmann constant, T is the absolute temperature, and is a measure of the probability that a channel is in the open state in the noload 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:
[ 8 ] 
[ 9 ] 
The tensions in the direction of the flow and orthogonal to flow in the plane of the membrane are denoted by T_{z} and T_{x}, respectively. In Eq. 8, the tension gradient in the direction of flow is related to the load fraction borne by the membrane, (1 ). The quantity is the fraction of the stress borne by the interior of the cell, which is assumed to have a solidlike 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 wellknown constitutive relation for cell membranes in pure shear deformation under uniaxial tension is given by Evans and Skalak (19).
[ 10 ] 
The quantity _{z} 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 twodimensional membrane. Integration of that expression and Eq. 8 yields the relationship between strain energy density and the applied shear stress.
[ 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 70 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, (t, ), is used herein. The temporal aspect of membrane permeability, P(t, ), is also described empirically. The expressions for (t, ) and P(t, ) 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, c_{s}, and on the influx of extracellular calcium, q_{in}. The constants required to calculate the effect of flow on c_{s} 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/m^{2}, in endothelial cells (3). In the absence of experimental data, a prediction of the magnitude of the flowinduced calcium influx is offered. The load fraction, , 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 fixed a priori, values for f_{e}, , and P_{max} are obtained by satisfaction of the following three criteria: q_{in}( = 0) = 5.33 µM/sec, q_{in}( = 3) = 10 µM/sec, q_{in,max} = 23 µM/sec.
The value of q_{in}( = 0) is that for passive leak of calcium into the cytosol (12). The value of q_{in}( = 3) was selected to obtain a value for the plateau level of calcium corresponding to an agonistfree experiment (24). The maximum value of calcium influx, q_{in,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 CRAYEL computer.
RESULTS
The calcium influx, q_{in}, 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 = 0.90, there is a threshold value of 0.5 N/m^{2} below which the influx is negligible. In the range of shear stress from 0.5 N/m^{2} to 6 N/m^{2}, the flux is shearsensitive. For shear stresses higher than the shearsensitive range, the influx is maximal. This corresponds to the situation in which all ion channels are open.
The effect of , the fraction of the load borne by the cytoplasm, is also shown in Fig. 2. The threshold value of shear stress and the sheardependent range are altered. As more of the stress is borne by the membrane, the threshold level of influx decreases. Also, the shearsensitive range narrows.
The response of cytosolicfree calcium is shown in Fig. 3. In this calculation, the perfusate is agonist free (c_{s} = 0) and contains physiologic calcium (Ca_{ex} = 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 serumfree media was observed in bovine aortic endothelial cells (25).
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 surfaceagonist 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/m^{2}. 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/m^{2} to 3 N/m^{2}. 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/m^{2} and 3 N/m^{2}. This indicates that the agonistmediated 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/m^{2}. This modulation of the agonist response by the level of shear stress has been observed experimentally (24). In contrast to the noagonist simulation, the peak here is very prominent. Also, there is no elevated plateau in the absence of extracellular calcium.^{ }
The model prediction is compared with the data in Fig. 6. The simulation for = 3 N/m^{2} 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 q_{in}. The plateau region is examined in more detail in the Discussion.
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 agoniststimulated 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 serumsupplemented media, which may contain calciummobilizing 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/m^{2}. 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 calciumdependent 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 ratelimiting processes in comparison to intracellular calcium dynamics. Resistance to diffusion of agonist is the controlling step at shear stresses below 0.1 N/m^{2}. The calcium transient moves from transportlimited to kinetically limited regime in the range of 0.1 N/m^{2} to 0.5 N/m^{2}. 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 microtubuleassociated 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 twodimensional stress field. This would be a worthwhile future enhancement to the model.
The magnitude of calcium influx is influenced by the loadbearing contribution of the cytoskeleton. This contribution is embodied in the parameter , 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 loadbearing 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 would help validate this hypothesis of mechanotransduction.
The simulations raise the question of what is the magnitude of shear stressinduced calcium influx. Schwarz et al. (3) reported a calcium current of 80 picoamps in cells exposed to a shear stress of 1 N/m^{2} 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 voltageclamp 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 stressinduced 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.
^{} 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 BCS911761 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 GrantinAid from the AHA.
HUVEC, human umbilical vein endothelial cell.
rating: 4.00 from 1 votes  updated on: 4 Dec 2007  views: 8036 