Department of Diagnostic Radiology, Cardiovascular Sciences Research
Group, University of Wales College of Medicine, Cardiff CF4 4XN,
United Kingdom
We have developed a mathematical model of
arterial vasomotion in which irregular rhythmic activity is generated
by the nonlinear interaction of intracellular and membrane oscillators
that depend on cyclic release of
Ca2+ from internal stores and
cyclic influx of extracellular
Ca2+, respectively. Four key
control variables were selected on the basis of the pharmacological
characteristics of histamine-induced vasomotion in rabbit ear arteries:
Ca2+ concentration in the cytosol,
Ca2+ concentration in
ryanodine-sensitive stores, cell membrane potential, and the open state
probability of Ca2+-activated
K+ channels. Although not
represented by independent dynamic variables, the model also
incorporates
Na+/Ca2+
exchange, the
Na+-K+-ATPase,
Cl
fluxes, and
Ca2+ efflux via the extrusion
ATPase. Simulations reproduce a wide spectrum of experimental
observations, including 1) the
effects of interventions that modulate the functionality of
Ca2+ stores and membrane ion
channels, 2) paradoxes such as the
apparently unpredictable dual action of
Ca2+ antagonists and low
extracellular Na+ concentration,
which can abolish vasomotion or promote the appearance of
large-amplitude oscillations, and 3)
period-doubling, quasiperiodic, and intermittent routes to chaos.
Nonlinearity is essential to explain these diverse patterns of
experimental vascular response.
nonlinear dynamics; vasomotion; calcium channels; potassium
channels; sodium/calcium exchange; calcium-adenosinetriphosphatase; sodium-potassium-adenosinetriphosphatase; chloride channels
 |
INTRODUCTION |
TISSUE PERFUSION is normally characterized by a marked
degree of spatial and temporal heterogeneity. Average local perfusion in the heart, for example, may vary by at least sixfold under resting
conditions, with the consequence that some microregions are permanently
vulnerable to hypoxia (21, 62). Superimposed time-dependent
fluctuations in perfusion that can encompass a wide range of
frequencies and amplitudes within the same vascular bed are also
apparent in myocardial elements as large as 1 cm3 and even at the segmental and
lobar levels in the pulmonary circulation (18, 57). Such temporal
variability results principally from spontaneous rhythmic fluctuations
in vascular tone and diameter, a phenomenon known as vasomotion, which
can be observed at all levels of the circulation. From a functional
point of view, the asynchronous opening and closing of resistance
arteries in different parts of the same vascular network ensures that
all tissue elements receive intermittent perfusion, and at the
microcirculatory level, vasomotion is thought to enhance fluid
filtration and lymphatic drainage (see Refs. 23 and 34 for review).
Theoretical simulations suggest that large-amplitude oscillations in
vessel diameter can increase time-averaged hydraulic conductance (55),
so that vasomotion may also represent an adaptive homeodynamic response
in pathological states such as hypoxia and hemorrhagic shock, in which
it becomes particularly pronounced (34).
Although the complex waveforms of vasomotion are often highly
irregular, there is accumulating evidence that they are not simply
random in origin but are generated by nonlinearity in the mechanisms
regulating smooth muscle tone. Indeed, fluctuations in pressure and
flow in isolated rabbit ear resistance arteries exhibit patterns that
are generic in nonlinear systems, including period-doubling,
quasiperiodic, and intermittent routes for the transition from regular
to irregular dynamics, so that vasomotion may be classified as a
deterministic or "chaotic" process (11, 14, 24-28).
Nonlinearity may confer specific biological advantages as chaotic fluid
flow accelerates diffusion-limited chemical reactions and mass-transfer
processes and could therefore enhance the delivery of oxygen and
nutrients and the removal of metabolites (23). Chaos may also allow
large changes in state to occur with minimal energy expenditure through
exploitation of the high sensitivity of nonlinear systems to
perturbation, thus enhancing the flexibility and efficiency of vascular
control (23, 26).
To gain insights into the nonlinear interactions that underlie
vasomotion at the cellular level, in the present study we have constructed a mathematical model of the ion transport systems that
participate in the regulation of vascular tone. The model simulates
fluctuations in free cytosolic
Ca2+ concentration
([Ca2+]i)
within a single myocyte, with changes in
[Ca2+]i
translated into changes in force generation via the latch state hypothesis of Hai and Murphy (29). The mathematical formulation was
assisted by the construction of a geometric phase space portrait of the
dynamics of experimental signals from isolated rabbit ear resistance
arteries by the technique of time-delayed coordinate embedding (11,
23). The complexity of this representation, which is known as an
"attractor," can be quantified as a nonintegral fractal
correlation dimension (here denoted as
DC) following nonlinear analysis
with the algorithm of Grassberger and Procaccia (22). DC provides an estimate of the
minimum number of participating dynamic control variables when rounded
to the nearest upper integer, a phase space of at least three
dimensions, i.e., DC > 2 being required to classify the behavior of a system as chaotic (22). In
general, the correlation dimension of fluctuations in pressure and flow
in rabbit arteries takes a value between 2 and 4 (14, 24-27), thus
characterizing vasomotion as a low-dimensional chaotic process that can
in theory be modeled with a "minimal" approach involving just
four independent variables. This conclusion is confirmed by
pharmacological analysis, which has revealed the participation of
"slow" intracellular and "fast" membranous
Ca2+ oscillators with periods on
the order of 1-5 min and 5-30 s, respectively (11, 25). In
the present model, two variables were therefore assigned to each
subsystem to permit independent oscillatory activity. The selection of
these variables was guided by the changes in
DC observed on administration of
specific pharmacological probes (see below).
 |
EXPERIMENTAL BACKGROUND TO THE MODEL |
Intracellular Oscillator
The intracellular oscillator was based on cyclic
Ca2+-induced
Ca2+ release (CICR) from
ryanodine-sensitive stores, with sustained oscillations having an
obligatory requirement for influx of
Ca2+ into the cytosol to replenish
the store after discharge. This flux was provided by
Ca2+ movements across the cell
membrane via voltage- and receptor-operated Ca2+ channels (VOCCs and ROCCs,
respectively),
Na+/Ca2+
exchange under conditions where it operated predominantly in influx
mode, and inositol 1,4,5-trisphosphate
(InsP3)-induced Ca2+ release (IICR) from
InsP3-sensitive stores. The major
ion fluxes incorporated in the model are illustrated schematically in
Fig. 1.

View larger version (33K):
[in this window]
[in a new window]
|
Fig. 1.
Schematic of key ion fluxes that contribute to 2 oscillatory subsystems
of model. Internal stores [i.e., sarcoplasmic reticulum
(SR)] possessing a ryanodine-sensitive
Ca2+-induced
Ca2+ release (CICR) mechanism are
essential for genesis of intracellular
Ca2+ oscillations, whereas stores
possessing inositol 1,4,5-trisphosphate
(InsP3) receptors are considered
to sequester Ca2+ avidly and serve
as a large "sink/source" that is permanently full. Agonists
stimulate Ca2+ influx through
voltage- and receptor-operated
Ca2+ channels (VOCC and ROCC,
respectively) and also cause
InsP3-induced
Ca2+ release (IICR) from stores
into cytosol. Cyclic changes in membrane potential
(Vm) and
Ca2+ influx via VOCC underpin
activity of membrane oscillator, with negative feedback on
Vm provided by
transport systems that promote hyperpolarization. These include
Ca2+-activated
K+ channels,
Na+-K+-ATPase,
and Cl channels.
Ca2+-ATPase pumps
( )
mediate uptake by stores and contribute to
Ca2+ extrusion from cell.
Na+/Ca2+
exchange can promote Ca2+
extrusion or influx, depending on instantaneous value of
Vm. Coupling
between 2 oscillatory subsystems is effected through cytosolic
Ca2+ concentration.
|
|
Ca2+ release
from stores.
The sarcoplasmic reticulum (SR) of vascular smooth muscle is
anatomically and functionally heterogeneous and possesses
Ca2+- and
InsP3-sensitive
Ca2+ stores (32, 38, 66).
Compounds that block sequestration of
Ca2+ by the SR
Ca2+-ATPase, such as cyclopiazonic
acid and thapsigargin, deplete both types of stores, whereas the
alkaloid ryanodine selectively impairs the functionality of a specific
Ca2+-release channel (7, 38). In
rabbit ear arteries, ryanodine appears to target a subcomponent of the
SR that is closely coupled to the contractile machinery of the smooth
muscle cell as it suppresses the activity of an intracellular
oscillator without affecting average perfusion pressure (25, 27). By
contrast, cyclopiazonic acid and thapsigargin induce a depressor
response and initially transform chaotic behavior to highly specific
patterns of nonlinear response known as mixed-mode oscillations before
high concentrations suppress vasomotion completely (11, 27). These
observations suggest that Ca2+ is
released from distinct store subtypes in rabbit ear arteries so that,
in the formulation of the present model, CICR and IICR were assumed to
originate from separate stores, as in the Goldbeter-Dupont-Berridge model of Ca2+ spiking (20).
Experimentally, ryanodine causes a concentration-dependent fall in the
average correlation dimension of rabbit ear artery vasomotion, i.e.,
DC, to <2, whereas SR
Ca2+-ATPase inhibitors do not
influence DC in vessels exhibiting
chaotic rhythmic activity until high concentrations induce relatively simple patterns of mixed-mode behavior (11, 25, 27). The key control
variables for the intracellular oscillator were therefore taken as
[Ca2+]i
and Ca2+ concentration in the
ryanodine-sensitive component of the SR ([Ca2+]SR).
Mathematically, CICR was modeled as a sigmoidal function of
[Ca2+]SR
multiplied by a sigmoidal function of
[Ca2+]i,
with sequestration by the
Ca2+-ATPase represented as a
sigmoidal function of
[Ca2+]i.
Although not essential for the genesis of oscillations, a continuous
leak of Ca2+ from the
ryanodine-sensitive store was also included, as in other models of this
type (20).
To model the InsP3-sensitive pool,
cytosolic InsP3 levels were
considered to be nonoscillatory and to increase directly with histamine
concentration, inasmuch as experimental evidence that periodic
fluctuations in InsP3
concentration contribute to oscillations in
[Ca2+]i
under physiological conditions has been obtained only in fibroblasts (30). Indeed, in other cell types, injection of nonmetabolizable analogs of InsP3 induces
fluctuations in Ca2+
concentration, confirming that oscillatory behavior can be triggered by
constant elevated levels of InsP3
(69). Importantly, the role of histamine in the genesis of vasomotion
appears to be permissive, inasmuch as changes in the concentration
employed to induce rhythmic activity in rabbit ear arteries influence
the superficial form of the responses observed but not their dynamical
complexity quantified as DC (24).
For the purposes of a minimal model, therefore, it was assumed that the
InsP3-sensitive store was
permanently full as a consequence of avid uptake and that the
Ca2+ flux into the cytosol via the
IICR mechanism could simply be modeled with a "constant" term
that increases with agonist concentration.
Voltage-operated
Ca2+ influx.
VOCCs provide the dominant pathway for
Ca2+ influx in vascular smooth
muscle, and the L-type Ca2+
antagonist verapamil can suppress contraction in rabbit ear arteries (11, 25). The open state probability of VOCCs is a steep sigmoidal function of cell membrane potential (51), and their contribution to
vasomotion was modeled as a noninactivating whole cell conductance with
a Bolzmann-like activation term, the slope of which as a function of
voltage was taken from Nelson et al. (51). In addition to regulating
VOCC conductance by causing membrane depolarization, constrictor
agonists may also promote the opening of such channels directly by
shifting their activation curve to the left, such that open state
probability is higher at a given voltage (51).
Receptor-operated
Ca2+ influx.
A receptor-gated cationic channel that is insensitive to membrane
depolarization and to inhibitors of VOCCs, is not modulated by
[Ca2+]i
or other second messengers, and possesses a threefold higher selectivity for Ca2+ than for
Na+ has been identified in
membrane patches from rabbit ear artery myocytes (6). Agonists such as
ATP and norepinephrine thus promote an inward flux of monovalent and
divalent cations, such that ~10% of the current stimulated under
conditions of voltage clamp is carried by
Ca2+ (3, 5, 6). Because the nature
of the ROCC remains poorly defined (51), in the present model
Ca2+ movements via this pathway
were simplified as a constant term, the magnitude of which was directly
related to agonist concentration. Under certain conditions, constrictor
agonists may stimulate Ca2+ influx
into vascular smooth muscle without causing depolarization (9), and it
was therefore assumed that receptor-operated
Ca2+ influx does not modulate
membrane potential.
Ca2+
extrusion ATPase.
In rabbit ear arteries, inhibition of the membrane
Ca2+ extrusion pump with the
vanadate ion evokes a large pressor response and a small decrease in
DC (14). In different simulations,
linear or parabolic activation curves were used to approximate the
dependence of the pump kinetics on
[Ca2+]i
over the physiologically relevant range (15), and its voltage dependence was formulated as a linear function of membrane potential, the slope of which was such as to give a 40% greater rate of
Ca2+ efflux at 0 mV than at
100 mV, in accord with observations that the rate of
extracellular Na+ concentration
([Na+]o)-independent
Ca2+ efflux from vascular myocytes
increases linearly with membrane potential (16). The pump exchanges
Ca2+ for
H+, but the stoichiometry of the
exchange may vary with local pH (15, 16), and under different
experimental conditions it has been reported to be electrogenic (16) or
electroneutral (53). The possible influence of
Ca2+ extrusion via this mechanism
on membrane potential was not incorporated in the simulations presented.
Na+/Ca2+
exchange.
Reductions in
[Na+]o
inhibit
Na+/Ca2+
exchange and progressively reduce
DC for vasomotion in rabbit ear
arteries to <2 (14). The exchanger can act as an efflux (forward
mode) or an influx (reverse mode) pathway for
Ca2+ movements, with the
instantaneous direction of the
Ca2+ flux depending on the
difference between membrane potential and the reversal potential of the
exchanger, which may vary between
30 and
45 mV according
to
[Ca2+]i
(4, 54). Inasmuch as the stoichiometry of the exchange is 3 Na+:1
Ca2+, below its reversal potential
this transport system will promote [Ca2+]i
extrusion associated with depolarization due to the net movement of
charge; above the reversal potential these effects will change direction. Consequently, in mathematical simulations the signs for the
contribution of this mechanism to the rate of change of [Ca2+]i
and membrane potential will always be opposed. Modulation of exchange
activity by local changes in Na+
concentration was ignored in the formulation of the present model.
Membrane Oscillator
The membrane oscillator that contributes to rabbit ear artery
vasomotion can be suppressed by concentrations of verapamil that exert
a submaximal effect on force development and by reductions in
extracellular Ca2+, with both
interventions causing DC to
decrease to <2, even when slow oscillations attributable to the
intracellular subsystem persist (25).
Ca2+ influx mediated by cyclic
opening of VOCCs thus appears to be central to the activity of the
membrane subsystem. In some artery types, electrical and mechanical
events are tightly coupled during vasomotion, with changes in membrane
potential preceding those in tension by 40-1,000 ms (35, 48). In
others, however, bursts of depolarization may be "integrated"
over time to produce tonic increases in force (37). Contraction may
also occur in the absence of membrane depolarization through
"pharmacomechanical" coupling (9). These observations suggest
that membrane potential may be considered as a distinct dynamic
variable that cannot be directly equated with changes in force.
To generate oscillatory electrical membrane activity that contributes
to rhythmic mechanical responses, mechanisms that promote hyperpolarization must operate to oppose voltage-dependent
Ca2+ influx. Experimental evidence
suggests that a spectrum of ion channels, pumps, and exchangers may act
in concert to provide such negative feedback (14). Reductions in the
activity of Ca2+-activated
K+
(KCa) channels,
Cl
channels, the
Na+-K+-ATPase,
or
Na+/Ca2+
exchange thus attenuate the membrane component of rabbit ear artery
vasomotion, with each class of intervention separately causing
DC to fall to <2 (14). Although
it is therefore clear that multiple channels, pumps, and exchangers
contribute to the regulation of tone in these vessels, it is likely
that the entrainment of weakly contributing variables results in the
suppression of a number of degrees of freedom of the system with a
corresponding reduction in its overall dynamic complexity (44). To
construct a minimal four-dimensional model that conforms to
experimental findings, the open state probability of the
KCa channel was chosen as the
fourth independent mathematical control variable. Membrane potential
can fluctuate between 0 and
50 mV during vasomotion (48, 64), so
that potassium ions, with reversal potential less than
80 mV
(52), are the only species able to drive the membrane to sufficiently
low potentials. Other membrane transport systems that influence
DC were included in the model as
modulatory fluxes.
Ionic current through KCa channels.
Multiple KCa channel isoforms may
coexist in the same artery type (52). In a detailed study, Toro et al.
(65) reconstructed the KCa
channels present in pig coronary artery in lipid bilayers and
distinguished different subtypes on the basis of their voltage and
Ca2+ sensitivity. Frequency
histograms of prevalence demonstrated two main populations of
"maxi" KCa channels with
conductances of 245 and 295 pS, each population consisting of two
isoforms. The 245-pS subtypes had half-maximal activation potentials of
80 and 6 mV, whereas the 295-pS subtypes had half-maximal
activation potentials of
28 and
66 mV (65). Subsidiary
channels were also found at ~50 and ~400 pS and also exhibited
differences in sensitivity to voltage and
[Ca2+]i.
The voltage slope coefficients of these channels were generally similar, with the channel exhibiting greatest sensitivity to
Ca2+ half-maximally activated at
~1.2 µM Ca2+
(65). In the present minimal model, representative values of these
parameters were selected to define the activation characteristics of a
single "composite" KCa
subtype. The Hill coefficient for the
Ca2+ dependence of the activation
sigmoidal was taken as 2 to conform with experimental findings (65).
The kinetics of channel inactivation were described by a first-order
rate equation with an exponential decay constant. Membrane potential
and
[Ca2+]i
may also influence the rate of KCa
channel closure (40, 52), but this additional complexity was not
incorporated in the model. Although voltage-dependent
K+ channels are functionally
active and modulate perfusion pressure in rabbit ear arteries, they do
not contribute to the dynamical complexity of vasomotion, perhaps
reflecting their insensitivity to
[Ca2+]i,
and were therefore not included in the model (14). ATP-sensitive K+ channels influence neither tone
nor vasomotion in rabbit ear arteries under well-oxygenated conditions
(25).
Ionic current through
Na+-K+-ATPase.
Under normal physiological circumstances the
Na+-K+-ATPase
is thought to contribute to vascular smooth muscle cell potential by
hyperpolarizing the membrane by up to
10 mV (1). Inasmuch as
changes in the local concentration of
K+ and
Na+ were not included as an
integral dynamical feature of the model, the contribution of this pump
was simplified as a constant hyperpolarizing flux.
Ionic current through Cl
channels.
A depolarizing Ca2+-activated
outward Cl
flux has been
reported in rabbit ear arteries (3, 39), and reductions in
extracellular Cl
concentration and niflumic acid, an inhibitor of
Ca2+-activated
Cl
channels, suppress the
activity of the membrane oscillator involved in vasomotion in rabbit
ear arteries (14). At membrane potentials more negative than the Nernst
reversal potential (approximately
25 mV), efflux of
Cl
promotes depolarization
of vascular smooth muscle, whereas at potentials more positive than
this value, efflux of chloride ions results in hyperpolarization (39).
Modulation of channel activity by
Ca2+ concentration was not
specifically incorporated in the present minimal model, however, and
Cl
fluxes were represented
simply as a constant term multiplied by the electrochemical gradient
for Cl
.
 |
METHODS |
Experimental
Experiments were performed with isolated ear preparations from male New
Zealand White rabbits killed by injection of pentobarbitone sodium (120 mg/kg iv), as previously described (24). First-generation vessels
(1-1.5 cm long and ~150-µm diameter) branching from the central artery were identified and perfused with oxygenated (95% O2-5%
CO2) Holman's buffer
(composition in mM: 120 NaCl, 5 KCl, 2.5 CaCl2, 1.3 NaH2PO4,
25 NaHCO3, 11 glucose, and 10 sucrose, pH 7.2-7.4) at 35°C in situ. The central ear artery
was tied off distally to divert all flow into a branch artery via a
proximal cannula. Side branches of this vessel were ligated, and its
distal end was cut to allow free outflow of perfusate. Mean flow was maintained at 0.5 ml/min in all experiments, and intrinsic fluctuations in flow and perfusion pressure due to vasomotion were monitored with a
Transonic flow probe and a pressure transducer connected via a sidearm
to the perfusion circuit. In all experiments, 50 µM
NG-nitro-L-arginine
methyl ester was included in the perfusate to inhibit endothelial
nitric oxide (NO) synthesis. Histamine,
NG-nitro-L-arginine
methyl ester, verapamil hydrochloride, ouabain, tetraethylammonium
chloride (TEA), sodium vanadate, and niflumic acid were obtained from
Sigma Chemical, and solutions were freshly prepared on the day of each
experiment. In some experiments, equimolar N-methyl-D-glucamine
was substituted for Na+ to provide
a low-Na+ buffer; the pH of the
perfusate was adjusted to 7.2-7.4 with hydrochloric acid.
Mathematical Formulation
As discussed above, the variables selected to simulate arterial
vasomotion were x
([Ca2+]i),
y
([Ca2+]SR),
z (membrane potential), and
w (open state probability of KCa channels). These variables
constitute two coupled oscillatory subsystems, the interactions of
which are described mathematically by the following set of nonlinear
differential equations.
Intracellular oscillator.
Ca2+ fluxes into the cytosol
Ca2+ uptake by stores
Membrane oscillator.
Relationship between ion fluxes and membrane
potential
Open state probability of KCa
channels
|
(1d)
|
All symbols are defined in Table 1, which
summarizes the fixed physiological values of the parameters (e.g.,
slope functions and reversal potentials) used for the simulations.
Ca2+
Buffering and Volume Scaling
Conventionally, Eq. 1c can be written
in a form that relates changes in membrane potential to ionic currents,
namely
|
(1c')
|
where Cm is cell
membrane capacitance. In Eq. 1c',
INa/K represents
the whole cell current for the
Na+-K+-ATPase
and each coefficient g represents a
whole cell conductance, the associated current
(ICl,
ICa,
INa/Ca, and
IK) of which is derived by multiplication by
an activation curve and electrochemical gradient (assumed constant
for Cl
,
Na+, and
K+, but not for
Ca2+).
Calcium ions that enter the cytosol are rapidly bound by intracellular
proteins (60). Indeed, experiments in gonadotrophs and chromaffin cells
suggest that ~1 in 100 calcium ions remains free in the cytosol (43).
Changes in
[Ca2+]i
may thus be related to the two electrogenic membrane
Ca2+ currents included in the
model, ICa and
INa/Ca, by
expressions of the form,
d[Ca2+]i/dt =
f
I, where f is the
fraction of calcium ions that are not buffered,
= 1/ZVcytF
is a scaling factor that converts current to rate of change in
[Ca2+]i,
Vcyt is the volume of the cytosol,
F is Faraday's constant, and
Z = 2 for
ICa and 1 for
INa/Ca (since the
stoichiometry of this exchanger is 3 Na+:1
Ca2+). From Eq.
1, a and c', it
follows that for each of these currents g = G/f
= ZVcytFG/f
and that
= VcytF/fCm.
The formulation of Eq. 1c thus
introduces a scaling factor,
, that involves the parameters of Vcyt,
Cm, and the fraction of calcium
ions that remains free (note that the divalent nature of
Ca2+, i.e., Z = 2 is
accounted for by the coefficient
2GCa in
Eq. 1c; for other mechanisms included in the
model, Z = 1). This "lumped" approach provides a convenient formulation for simulations, inasmuch as
the contributing parameters are known only approximately. Indeed, to
obtain realistic oscillations that matched experimental data, it was
necessary to adjust
over the broad range of 0.2-8 (Table 2). In practice,
exerts an important
influence on the relative periods of the two contributing oscillators,
which can exhibit significant variation between different experimental
preparations.
Buffering also occurs within the SR, and there is evidence from COS-7
cells that only ~1 in 400 calcium ions is free in this compartment,
compared with 1 in 100 in the cytosol, because of the presence of
specialized Ca2+-buffering
proteins (36, 41). For the purposes of modeling, it therefore becomes
necessary to define "effective volumes" available for
Ca2+ movements (12) that are
equivalent to 100 and 400 times the physical volumes of these
compartments (Vecyt and
VeSR), respectively, whereas the ratio of
their actual volumes
(VSR/Vcyt) is <1:4, because the SR occupies <20% of the cytosol. It follows that VeSR ~ Vecyt, which permits the assumption that
Ca2+ concentration changes in the
cytosol and stores resulting from intracellular movements between
these two compartments are directly related according to
|
(2a)
|
as
is implicit in Eq. 1,
a and
b, if transmembrane fluxes are
ignored. A more general relationship of the form
|
(2b)
|
could
take into account structural differences between cells.
Parameter Estimation
The maintenance of force development in rabbit ear arteries is
crucially dependent on extracellular
Ca2+ influx via L-type
Ca2+ channels (11, 25). The flux
through a single open VOCC has been estimated at
~105 ions/s, and there may be
~1,000 such channels in a vascular myocyte, resulting in a flux of
~106 ions/s if ~1% of
channels are assumed to be open at about
40 mV (51). Vascular
smooth muscle cells are approximately cylindrical in shape, with a
radius of ~2.5 µm and a length of ~60 µm, giving a surface area
of ~10
5
cm2 and a volume of ~1.2 × 10
12 liters. If it is
assumed that the cytosol occupies 20% of the total cell volume and
that ~1% of the calcium ions entering the cell escape buffering by
proteins, this flux would increase
[Ca2+]i
at a rate
d[Ca2+]i/dt = 0.01[106/(0.2 × 1.2 × 10
12)][1/(6.24 × 1023)]
0.1 µM/s, where Avogadro's number = 6.24 × 1023
mol
1.
In the model the contribution of
Ca2+ influx through VOCCs to
d[Ca2+]i/dt is
GCa(z
zCa1)/{1 + e[
(z
zCa2)/RCa]}, which reduces to
0.02GCa µM/s
for the fixed values
zCa2
and RCa used in
the simulations at an "average" membrane potential of
40
mV (Fig. 2). In practice,
GCa was varied
over the range 1-16 (Table 2), giving values of
d[Ca2+]i/dt = 0.02-0.32 µM/s, which encompass the value expected from the
above order-of-magnitude calculation. To reflect the dominance of
Ca2+ influx via VOCCs, this flux
was generally made up to four times larger (depending on the
instantaneous value of membrane potential) than the combined
Ca2+ flux entering the cytosol via
InsP3-induced release from stores and receptor-operated channels.

View larger version (35K):
[in this window]
[in a new window]
|
Fig. 2.
A: phase plane plot for intracellular
oscillator. Trajectories of free cytosolic
Ca2+ concentration
([Ca2+]i)
vs. Ca2+ concentration in SR
([Ca2+]SR;
x vs.
y) generate a limit cycle
oscillation that is traversed in direction shown by arrows. Three
principal regimens can be identified. Dashed and dotted curves,
x and
y nullclines, respectively;
oscillations occur only when their intersection ( ) is within
"active ranges" of both (highlighted).
B: oscillations in
[Ca2+]i
in time domain. In this simulation, membrane oscillator was
"switched off" to reflect dynamics of
[Ca2+]i
oscillator in isolation. Consequently, signal generated is periodic in
form. C and
D: phase space plot for membrane
oscillator and its associated temporal waveform, respectively. Dashed
line, z nullcline; dotted line,
w nullcline. , Equilibrium
point; active range of each nullcline is highlighted. Limit cycle is
traversed in direction shown by arrows. In this simulation,
intracellular oscillator was switched off. KCa,
Ca2+-activated K+ channel;
Vm, membrane potential.
|
|
The magnitude of the Ca2+ flux via
VOCCs employed in the above calculation can be used to estimate the
rate of change of membrane potential during the upstroke of the
vasomotion cycle. If it is assumed that the cell has a
Cm of 1 µF/cm2 (43) and a
surface area of ~10
5
cm2, then total
Cm
10
5 µF. For a flux of
n = 106 ions/s, the current carried by
Ca2+ is
ICa = 2ne = 3.2 × 10
13 C/s [where the
electronic charge (e) = 1.6 × 10
19 C]. It follows
that
dz/dt = ICa/Cm
32 mV/s. Given that opposing ion transport systems that promote
membrane hyperpolarization are simultaneously active, this value is
consistent with experimental measurements that range from 4 to 10 mV/s:
~4 mV/s (Fig. 1 in Ref. 64), ~5 mV/s (Fig. 2 in Ref. 17), and ~10
mV/s (Fig. 2 in Ref. 70).
A semiempirical approach was adopted to estimate the operating values
of the remaining coefficients in Eq.
1c, inasmuch as, in practice, such variables as
individual channel conductances and numbers of channel per cell are, in
general, unknown. The model thus reproduces the experimental changes in
membrane potential observed when the different ion transport systems
are inhibited pharmacologically (see
RESULTS). The coefficients chosen to
reflect the magnitudes of net Ca2+
extrusion across the membrane by the
Ca2+-ATPase and
Na+/Ca2+
exchange and uptake by the SR also match relative estimates of these
fluxes in vascular smooth muscle (see Inhibition of
the membrane
Ca2+-ATPase).
Integration of
Ca2+ Dynamics by
Latch Bridge Formation
Equation 1, a-d, generates
oscillations in
[Ca2+]i,
[Ca2+]SR,
membrane potential, and the open state probability of
KCa channels and does not simulate
experimental fluctuations in arterial pressure and flow directly,
although these will closely follow oscillations in
[Ca2+]i.
Force development and
[Ca2+]i
were therefore related by the cross-bridge phosphorylation and latch
state model proposed by Hai and Murphy (29). In this model, elevated
[Ca2+]i
is coupled to smooth muscle tone through the formation of cross bridges
between the actin and myosin filaments of the contractile machinery.
There are four possible states for myosin: free nonphosphorylated cross
bridges (M), free phosporylated cross bridges (Mp), attached phosphorylated cross bridges (AMp), and dephosphorylated latch bridges
(AM). The kinetics relating these states are as
follows
where
K1-K7
are the rate constants regulating the phosphorylation and bridge
formation and can be written as the following set of differential
equations
|
(3a)
|
|
(3b)
|
|
(3c)
|
|
(3d)
|
|
(3e)
|
Under the constraint of Eq. 3e, the
system effectively reduces to three differential equations. Force
development was calculated as the sum [AM] + [AMp], which reflects total cross-bridge formation. The
phosphorylation rates of M and AM (rate constants
K1 and
K6) are
Ca2+ concentration dependent and
assumed to follow the sigmoidal
[Ca2+]2i/([Ca2+]2i + x20), which implies
saturation at high
[Ca2+]i.
The values for the rate constants obtained by Hai and Murphy (29) for
the swine carotid artery were
K1 = K6 = 1.7 s
1,
K2 = K5 = 0.5 s
1,
K3 = 0.4 s
1,
K4 = 0.1 s
1,
K7 = 0.02 s
1, and
x0 = 0.6 µM. In
practice, it becomes apparent that the key coefficient responsible for
integrating Ca2+ fluctuations is
K7, and in
different simulations the value of this parameter was varied to obtain
different degrees of dampening of the
Ca2+ oscillations (Table 2). The
latch equations are external to Eq. 1, a-d, and effectively integrate oscillations in
[Ca2+]i
over time, thus providing the system with a short-term memory over
previous cycles. However, they will have no fundamental effect on the
dynamical complexity of the simulations.
Numerical Methods
Simulations were produced by numerical integration of the set of
Eq. 1 with the Runge-Kutta-Merson
algorithm by employing an integration step in the range of
0.0005-0.003 s. The initial conditions and coefficients employed
in the simulations are given in Table 2. Pharmacological interventions
were simulated by decreasing or increasing the appropriate
coefficient(s) linearly over time according to the ion transport system
under study. This is an approximation of the experimental situation;
pharmacological inhibitors were added to the buffer reservoir in a
graded fashion, so that the concentration perceived by the isolated
arteries will in practice increase gradually to a plateau as a
consequence of mixing within the perfusion circuit. The exact point in
time at which each simulated intervention was initiated is indicated by
an arrow on the corresponding figure with the initial and final values
of the appropriate coefficient(s) given in the figure legend.
Analytic solution of Eq. 1, a-d,
is impossible, inasmuch as the equations constitute a nonlinear
fourth-order system. To gain insights into underlying mechanisms, each
contributing oscillator was therefore first examined individually by
reducing the four-dimensional system to two two-dimensional projections
of the attractor. Phase space plots of the variable pairs
x vs.
y and
z vs.
w were thus obtained to evaluate the
essential dynamical features of each subsystem and investigate
systematically the consequences of changes in the coupling between them.
 |
RESULTS |
Stability Analysis
Equation 1, a and b.
The equations describing the intracellular subsystem generate a limit
cycle oscillation that exhibits three distinct regimens that are
clearly identifiable in the phase plane plot of Fig. 2A and in the time domain
representation of Fig. 2B. Along the slow section the dynamics are dominated by sequestration
of Ca2+ entering the cell within
the SR, with the result that
[Ca2+]i
remains almost constant. This phase is followed by a "very fast"
section, in which the SR empties abruptly into the cytosol. In the
ensuing fast section, Ca2+ is
extruded from the cell more rapidly than it enters from the extracellular space, so that
[Ca2+]i
declines and
[Ca2+]SR
remains low. When
[Ca2+]i
becomes sufficiently low, uptake again dominates over extrusion, thus
permitting sequestration of Ca2+
within the SR and allowing a new cycle to ensue.
The limit cycle is formed by the rotation of the trajectories of the
system around an unstable equilibrium point
(x*,y*) located at the intersection
of the x and
y nullclines, which represent the loci
of points where the rates of change of
[Ca2+]i
and
[Ca2+]SR
are zero (Fig. 2A). Analytically,
this equilibrium point is thus defined by
(dx/dt)|
= (dy/dt)|
= 0. Stability theory indicates that oscillatory behavior will be found
only when the eigenvalues characterizing the deformation of the phase
space in the neighborhood of the equilibrium point (x*,y*) are complex conjugates with
a positive real part. These eigenvalues can be determined by
linearization of Eq. 1, a and b, in the vicinity of the fixed point
according to
|
(4a)
|
|
(4b)
|
where
= x
x* and
= y
y*. The Jacobian matrix of the system
is
|
(5)
|
The determinant (Det = a11a22
a12a21)
and the trace (Tr = a11 + a22) of this
matrix characterize the nature of the phase space trajectories around
the equilibrium point (56). Thus 1) if Det > Tr2/4, the solution is
oscillatory, decaying when Tr < 0 (focal stability) and growing when
Tr > 0 (focal instability). Focal instability is the only
configuration that leads to sustained oscillatory activity.
2) If 0 < Det < Tr2/4, the solution is
nonoscillatory, with nodal stability for Tr < 0 and nodal instability
for Tr > 0. 3) If Det < 0, the
system always exhibits saddle-point instability; i.e., it attracts
trajectories from some directions but repels them from others.
An alternative way to present these conditions is in terms of the
slopes of the nullclines in the vicinity of the equilibrium point, where
d
/dt
and
d
/dt
become zero. Equation 4, a and
b, can then be written in the form
|
(6a)
|
|
(6b)
|
For the specific configuration of Fig.
2A in which the slopes of both
nullclines near the equilibrium point are negative, it is readily shown
that their intersection will be an unstable focus only when
|
(7a)
|
|
(7b)
|
The
stability of the fixed point can thus be assessed by inspecting the
slopes of the nullclines near their point of intersection on the phase
plane. The "active" regions satisfying the conditions of
Eq. 7, a and
b, where oscillations are possible are
illustrated in Fig. 2A. Although the
equilibrium point is unstable, the trajectories of the system do not
diverge to infinity but are limited by the overall shape of the phase
space, which creates a bounded basin for the trajectories, thus
permitting sustained oscillatory behavior. In Fig.
2A it can thus be seen how
trajectories of the system are limited by the
x and
y nullclines. At low values of
x, they are forced to follow the
x nullcline closely until just past
its peak before the cycle can be completed, whereas the
y nullcline imposes a similar limit
for the trajectories traversing the phase plane at low
y values. These two limits, along with
the very fast linear segment that corresponds to emptying of the SR
into the cytosol, confer a characteristic triangular shape to the limit cycle. During the latter phase of the dynamics, the ryanodine-sensitive store empties into the cytosol on a time scale that is much faster than
Ca2+ movements across the
membrane, so that
[Ca2+]i + [Ca2+]SR
effectively remains constant. It follows that the slope of the very
fast phase of the limit cycle in Fig.
2A is
1 for all values of
x* and
y*.
Equation 1, c and d.
A phase plane representation of the dynamics of the membrane oscillator
given by Eq. 1, c and
d, is shown in Fig.
2C. As in the case of the
intracellular oscillator, the conditions necessary for oscillatory
behavior can be obtained by considering the
z and
w nullclines, i.e., the loci of points
where the rates of change of membrane potential and
KCa channel open state probability are zero. The equilibrium point
(z*,w*)
is defined by
(dz/dt)|
= (dw/dt)|
= 0.
It is readily shown that this will be an unstable focus that permits
sustained oscillations only when the slopes µ and
of the
z and
w nullclines at their point of
intersection are positive and the conditions
|
(8a)
|
|
(8b)
|
are satisfied.
The nullclines of the membrane oscillator possess greater symmetry
around the equilibrium point than those of the intracellular subsystem,
resulting in a more uniform speed around the limit cycle and nearly
sinusoidal oscillatory behavior, which is delimited by the local
maximum and minimum of the z nullcline
(Fig. 2D). The crucial mechanism
that permits the emergence of oscillatory membrane activity is the
intrinsic time-dependent inactivation of the
KCa channel, which introduces
hysteresis into the system by delaying full opening of the channel
after depolarization. The coefficients
and
in
Eq. 1, c and
d, which scale the rate of change of
membrane potential and the rate of change of the KCa channel open state
probability, respectively, contribute principally to the speed of the
membrane oscillations, whereas the ratio of the rates of channel
activation to inactivation (S) in
Eq. 1d principally determines their
form. If
or S is reduced,
oscillations exhibit progressively smaller limit cycles, leading
ultimately to a steady state. To maintain the frequency of the membrane
oscillator within an appropriate range during the simulations, it was
found necessary to limit variations in the product
GK, which
governs the net contribution of
KCa channel opening to
dz/dt
in Eq. 1c (Table 2).
The dynamics of the membrane oscillator are less robust than the
dynamics generated by intracellular
Ca2+ cycling, because the slope of
the sigmoidal KCa channel
activation curve and its position on the
z-axis are highly sensitive to
[Ca2+]i.
The equilibrium point of the membrane nullclines may thus translate
substantially as a result of modulation by the intracellular oscillator. The x (i.e.,
Ca2+) dependence of
Eq. 1d at the equilibrium point can be
understood by applying the transformation
zx = RKln{
/[(x + xw)2]}
and expressing the w nullcline in the
form
|
(9)
|
It follows that changes in
[Ca2+]i
translate the inflection point
zCa3
of the KCa channel activation
sigmoidal horizontally, such that the nullcline moves from right to
left for increasing values of x.
Latch bridge dynamics.
Figure 3 illustrates the integrating effect
of the latch bridge mechanism on oscillations in
[Ca2+]i
in which the parameters in Eq. 1 were
selected to generate fast membrane oscillations with an amplitude
~20% of the slow intracellular oscillator. The top
traces show how a high value of the latch bridge decay
constant K7 (0.2)
selectively dampens fluctuations attributable to the membrane subsystem
with little effect on the position of the extrema of the slow
oscillations; i.e., the slow component of the force and signals remain
almost in phase. Low values of
K7 (e.g., 0.01)
cause further suppression of membrane oscillations, a phase shift to
the right in the slow component and also a small rise in average tone.
In physiological terms, high values of
K7 correspond to
rapid decay of the latch bridge, so that the system is capable of
responding to rapid changes in
[Ca2+]i.
For low values of
K7 the latch
decays slowly, allowing the system to function as an integrator of
[Ca2+]i.

View larger version (37K):
[in this window]
[in a new window]
|
Fig. 3.
Representative simulations showing force derived from
Ca2+ oscillations by integrating
system of Eq. 2 for different values
of K7 (i.e.,
latch bridge decay constant).
|
|
Simulations of Pharmacological Interventions
Inhibition of CICR.
Figure 4,
A and
D, illustrates the experimental
effects of ryanodine on patterns of vasomotion that exhibit different
relative contributions from the two participating oscillators. In Fig. 4A the responses are dominated by the
intracellular subsystem, inasmuch as all activity was abolished by
ryanodine. In Fig. 4D both subsystems
are active and ryanodine selectively attenuated the contribution of the
intracellular oscillator but not membrane activity. In neither artery
was there a major change in mean perfusion pressure.

View larger version (33K):
[in this window]
[in a new window]
|
Fig. 4.
A-C: experimental trace and 2 simulations illustrating suppression of intracellular oscillator by
ryanodine. D-F: experimental
trace and 2 simulations showing effects of ryanodine when initial
dynamics involved intracellular and membrane subsystems. Simulations in
B and
E were obtained by linearly decreasing
coefficient C in Eq.
1, a and b, from 5,000 to 1,000 and from 6,250 to 1,250, respectively, over time period shown;
in C and
F, leak coefficient
(L) was additionally increased from
0.025 to 0.055, again in a linear fashion. Traces in
D-F illustrate a reverse
period-doubling cascade (chaos P2 P1). Changes in
coefficients were initiated at points indicated by arrows.
|
|
Low concentrations of ryanodine lock the
Ca2+ release channel in an open
subconductance state, whereas high concentrations cause channel closure
(7). This biphasic effect was simulated
1) by linearly decreasing
coefficient C in Eq.
1, a and
b, over time (Fig. 4,
B and
E) and
2) by linearly decreasing
coefficient C and simultaneously
increasing the "leak" coefficient
L (Fig. 4, C and
F). Modification of the leak term
resulted in slightly faster suppression of the intracellular oscillator
(compare, for example, the height of the force maxima on simulating the
introduction of ryanodine in Fig. 4, B
and C), but in neither case was
there a major effect on the general form of the dynamics. The two
distinct patterns of oscillatory behavior simulated in Fig. 4 differ
principally in the value of the coefficient
(
GCa,
GNa/Ca, and
GK were identical in both cases; Table 2). A large value of
effectively corresponds to a cell in which the change in membrane potential effected by a given absolute ionic flux is associated with a relatively small change in
[Ca2+]i;
i.e., the cell is large, or there is a high degree of buffering by
intracellular proteins. Consequently, the membrane component of
vasomotion is much less apparent in Fig. 4,
A-C (
= 8.35) than in
D-F (
= 0.34).
Phase plane analysis shows that ryanodine moves the
[Ca2+]i
and
[Ca2+]SR
nullclines upward on the ordinate simultaneously whether its action is
modeled by decreasing C alone or by
decreasing C and increasing
L together (Fig.
5, A and
B). This implies the following predictions: 1) the limit cycle
remains almost unchanged in size until close to the position where the
equilibrium point becomes a stable focus and oscillatory behavior
ceases, and 2) there will be no
change in ambient force development, inasmuch as this is determined by
average
[Ca2+]i
at the intersection point of the nullclines. Both of these predictions
parallel the experimental findings presented in Fig. 4,
A and
D. The lack of effect of ryanodine on
force development in simulations is to be expected on theoretical
grounds. From the mathematical definition of the equilibrium point, it
follows that
(dx/dt)|(x*,y*) + (dy/dt)|(x*,y*) = 0, so that, on adding Eq. 1, a and
b, we obtain
|
(10)
|
This expression shows analytically that the equilibrium
value x* is independent of the
coefficients B, C, and
L, thus confirming that sequestration
and release of Ca2+ by the
ryanodine-sensitive store do not influence average
[Ca2+]i.

View larger version (18K):
[in this window]
[in a new window]
|
Fig. 5.
Phase space interpretation of effects of ryanodine. On reducing
C or reducing
C and increasing
L simultaneously
(A and
B, respectively), intersection of
x and
y nullclines ( ) moves vertically
upward until oscillatory behavior ceases relatively abruptly.
L alters slopes of nullclines at high
[Ca2+]i,
but there is almost no change in shape and size of limit cycle.
|
|
Blockade of VOCCs.
In some experimental preparations the
Ca2+ antagonist verapamil
abolishes all oscillatory activity, whereas in others the amplitude of
the slow component of vasomotion is enhanced and membrane activity is
attenuated (Fig. 6). These distinct classes
of response are associated with a fall in mean perfusion pressure and
could be successfully simulated by reducing the coefficient
GCa, which determines the magnitude of Ca2+
influx via VOCCs in Eq. 1, a and
c (Fig. 6,
B and
D). The variable outcomes in such
simulations are explained by differences in the balance between
Ca2+ influx via VOCCs and other
pathways, i.e., via ROCCs, Ca2+
release from the InsP3-sensitive
store, or reverse-mode
Na+/Ca2+
exchange (respectively represented by the coefficients
A0,
A1, and
GNa/Ca). When
the net Ca2+ flux is small,
blockade of VOCCs may prevent adequate replenishment of stores, so that
vasomotion ceases. Conversely, if the net
Ca2+ flux into the cytosol via
VOCC-independent pathways is large, intracellular oscillations may
persist even in the presence of verapamil. Indeed, the critical
difference between the simulations shown in Fig. 6,
B and
D, is that the contribution of the
Na+/Ca2+
exchanger is nearly fourfold greater and the contribution of A0 + A1 ~50%
greater in the case where oscillatory behavior is maintained. Phase
space analysis provides the dynamical explanation for this apparently
paradoxical experimental behavior (Fig. 7). As a consequence of the hyperbolic shape of the active ranges of the
x and
y nullclines, verapamil moves the
equilibrium point of Eq. 1, a, and
b, to higher
[Ca2+]SR
and lower
[Ca2+]i
with two quite distinct effects: in Fig.
7A the equilibrium point is forced
beyond the active oscillatory area, thus being transformed from an
unstable focus to a stable focus or node that attracts trajectories,
with the result that oscillations cease. Conversely, the equilibrium
point may remain within the active region, as an unstable focus, even
with GCa = 0 (Fig. 7B). Movement of the
equilibrium point upward and to the left is accompanied by an increase
in the amplitude of oscillatory behavior as mean [Ca2+]i
and force development decrease, consistent with the experimental trace
shown in Fig. 6D. The increase in the
size of the limit cycle is also consistent with the slight increment in
the period of the slow component of vasomotion observed in Fig. 6,
C and D.

View larger version (33K):
[in this window]
[in a new window]
|
Fig. 6.
A and
C: experimental traces illustrating
effects of verapamil, which may abolish vasomotion or induce
large-amplitude oscillations. B and
D: model-generated force oscillations
simulating this paradoxical action of verapamil. These were obtained by
linearly decreasing coefficient
GCa in
Eq |
|