## Abstract

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 Ca^{2+} from internal stores and cyclic influx of extracellular Ca^{2+}, respectively. Four key control variables were selected on the basis of the pharmacological characteristics of histamine-induced vasomotion in rabbit ear arteries: Ca^{2+} concentration in the cytosol, Ca^{2+} concentration in ryanodine-sensitive stores, cell membrane potential, and the open state probability of Ca^{2+}-activated K^{+} channels. Although not represented by independent dynamic variables, the model also incorporates Na^{+}/Ca^{2+}exchange, the Na^{+}-K^{+}-ATPase, Cl^{−} fluxes, and Ca^{2+} efflux via the extrusion ATPase. Simulations reproduce a wide spectrum of experimental observations, including *1*) the effects of interventions that modulate the functionality of Ca^{2+} stores and membrane ion channels, *2*) paradoxes such as the apparently unpredictable dual action of Ca^{2+} 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

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 cm^{3} 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 Ca^{2+} concentration ([Ca^{2+}]_{i}) within a single myocyte, with changes in [Ca^{2+}]_{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 D_{C}) following nonlinear analysis with the algorithm of Grassberger and Procaccia (22). D_{C} 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., D_{C} > 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 Ca^{2+} 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 D_{C} observed on administration of specific pharmacological probes (see below).

## EXPERIMENTAL BACKGROUND TO THE MODEL

### Intracellular Oscillator

The intracellular oscillator was based on cyclic Ca^{2+}-induced Ca^{2+} release (CICR) from ryanodine-sensitive stores, with sustained oscillations having an obligatory requirement for influx of Ca^{2+} into the cytosol to replenish the store after discharge. This flux was provided by Ca^{2+} movements across the cell membrane via voltage- and receptor-operated Ca^{2+} channels (VOCCs and ROCCs, respectively), Na^{+}/Ca^{2+}exchange under conditions where it operated predominantly in influx mode, and inositol 1,4,5-trisphosphate (InsP_{3})-induced Ca^{2+} release (IICR) from InsP_{3}-sensitive stores. The major ion fluxes incorporated in the model are illustrated schematically in Fig. 1.

#### Ca^{2+} release from stores.

The sarcoplasmic reticulum (SR) of vascular smooth muscle is anatomically and functionally heterogeneous and possesses Ca^{2+}- and InsP_{3}-sensitive Ca^{2+} stores (32, 38, 66). Compounds that block sequestration of Ca^{2+} by the SR Ca^{2+}-ATPase, such as cyclopiazonic acid and thapsigargin, deplete both types of stores, whereas the alkaloid ryanodine selectively impairs the functionality of a specific Ca^{2+}-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 Ca^{2+} 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 Ca^{2+} spiking (20).

Experimentally, ryanodine causes a concentration-dependent fall in the average correlation dimension of rabbit ear artery vasomotion, i.e., D_{C}, to <2, whereas SR Ca^{2+}-ATPase inhibitors do not influence D_{C} 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 [Ca^{2+}]_{i}and Ca^{2+} concentration in the ryanodine-sensitive component of the SR ([Ca^{2+}]_{SR}). Mathematically, CICR was modeled as a sigmoidal function of [Ca^{2+}]_{SR}multiplied by a sigmoidal function of [Ca^{2+}]_{i}, with sequestration by the Ca^{2+}-ATPase represented as a sigmoidal function of [Ca^{2+}]_{i}. Although not essential for the genesis of oscillations, a continuous leak of Ca^{2+} from the ryanodine-sensitive store was also included, as in other models of this type (20).

To model the InsP_{3}-sensitive pool, cytosolic InsP_{3} levels were considered to be nonoscillatory and to increase directly with histamine concentration, inasmuch as experimental evidence that periodic fluctuations in InsP_{3}concentration contribute to oscillations in [Ca^{2+}]_{i}under physiological conditions has been obtained only in fibroblasts (30). Indeed, in other cell types, injection of nonmetabolizable analogs of InsP_{3} induces fluctuations in Ca^{2+}concentration, confirming that oscillatory behavior can be triggered by constant elevated levels of InsP_{3}(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 D_{C} (24). For the purposes of a minimal model, therefore, it was assumed that the InsP_{3}-sensitive store was permanently full as a consequence of avid uptake and that the Ca^{2+} flux into the cytosol via the IICR mechanism could simply be modeled with a “constant” term that increases with agonist concentration.

#### Voltage-operated Ca^{2+} influx.

VOCCs provide the dominant pathway for Ca^{2+} influx in vascular smooth muscle, and the L-type Ca^{2+}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 Ca^{2+} influx.

A receptor-gated cationic channel that is insensitive to membrane depolarization and to inhibitors of VOCCs, is not modulated by [Ca^{2+}]_{i}or other second messengers, and possesses a threefold higher selectivity for Ca^{2+} 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 Ca^{2+} (3, 5, 6). Because the nature of the ROCC remains poorly defined (51), in the present model Ca^{2+} 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 Ca^{2+} influx into vascular smooth muscle without causing depolarization (9), and it was therefore assumed that receptor-operated Ca^{2+} influx does not modulate membrane potential.

#### Ca^{2+}extrusion ATPase.

In rabbit ear arteries, inhibition of the membrane Ca^{2+} extrusion pump with the vanadate ion evokes a large pressor response and a small decrease in D_{C} (14). In different simulations, linear or parabolic activation curves were used to approximate the dependence of the pump kinetics on [Ca^{2+}]_{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 Ca^{2+} efflux at 0 mV than at −100 mV, in accord with observations that the rate of extracellular Na^{+} concentration ([Na^{+}]_{o})-independent Ca^{2+} efflux from vascular myocytes increases linearly with membrane potential (16). The pump exchanges Ca^{2+} 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 Ca^{2+} extrusion via this mechanism on membrane potential was not incorporated in the simulations presented.

#### Na^{+}/Ca^{2+}exchange.

Reductions in [Na^{+}]_{o}inhibit Na^{+}/Ca^{2+}exchange and progressively reduce D_{C} 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 Ca^{2+} movements, with the instantaneous direction of the Ca^{2+} 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 [Ca^{2+}]_{i}(4, 54). Inasmuch as the stoichiometry of the exchange is 3 Na^{+}:1 Ca^{2+}, below its reversal potential this transport system will promote [Ca^{2+}]_{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 [Ca^{2+}]_{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 Ca^{2+}, with both interventions causing D_{C} to decrease to <2, even when slow oscillations attributable to the intracellular subsystem persist (25). Ca^{2+} 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 Ca^{2+} 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 Ca^{2+}-activated K^{+}(K_{Ca}) channels, Cl^{−} channels, the Na^{+}-K^{+}-ATPase, or Na^{+}/Ca^{2+}exchange thus attenuate the membrane component of rabbit ear artery vasomotion, with each class of intervention separately causing D_{C} 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 K_{Ca} 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 D_{C} were included in the model as modulatory fluxes.

#### Ionic current through K_{Ca} channels.

Multiple K_{Ca} channel isoforms may coexist in the same artery type (52). In a detailed study, Toro et al. (65) reconstructed the K_{Ca}channels present in pig coronary artery in lipid bilayers and distinguished different subtypes on the basis of their voltage and Ca^{2+} sensitivity. Frequency histograms of prevalence demonstrated two main populations of “maxi” K_{Ca} 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 [Ca^{2+}]_{i}. The voltage slope coefficients of these channels were generally similar, with the channel exhibiting greatest sensitivity to Ca^{2+} half-maximally activated at ∼1.2 μM Ca^{2+}(65). In the present minimal model, representative values of these parameters were selected to define the activation characteristics of a single “composite” K_{Ca}subtype. The Hill coefficient for the Ca^{2+} 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 [Ca^{2+}]_{i}may also influence the rate of K_{Ca}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 [Ca^{2+}]_{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 Ca^{2+}-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 Ca^{2+}-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 Ca^{2+} 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% O_{2}-5% CO_{2}) Holman’s buffer (composition in mM: 120 NaCl, 5 KCl, 2.5 CaCl_{2}, 1.3 NaH_{2}PO_{4}, 25 NaHCO_{3}, 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*N*
^{G}-nitro-l-arginine methyl ester was included in the perfusate to inhibit endothelial nitric oxide (NO) synthesis. Histamine,*N*
^{G}-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*([Ca^{2+}]_{i}),*y*([Ca^{2+}]_{SR}),*z* (membrane potential), and*w* (open state probability of K_{Ca} 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.

Ca^{2+} fluxes into the cytosol

Ca^{2+} uptake by stores

*Membrane oscillator. *

Relationship between ion fluxes and membrane potential

Open state probability of K_{Ca}channels
Equation 1dAll 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.

### Ca^{2+}Buffering and Volume Scaling

Conventionally, *Eq. 1c* can be written in a form that relates changes in membrane potential to ionic currents, namely
Equation 1c' where C_{m} is cell membrane capacitance. In *Eq. 1c′*,*I*
_{Na/K} represents the whole cell current for the Na^{+}-K^{+}-ATPase and each coefficient *g* represents a whole cell conductance, the associated current (*I*
_{Cl},*I*
_{Ca},*I*
_{Na/Ca}, and*I*
_{K}) of which is derived by multiplication by an activation curve and electrochemical gradient (assumed constant for Cl^{−}, Na^{+}, and K^{+}, but not for Ca^{2+}).

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 [Ca^{2+}]_{i}may thus be related to the two electrogenic membrane Ca^{2+} currents included in the model, *I*
_{Ca} and*I*
_{Na/Ca}, by expressions of the form, d[Ca^{2+}]_{i}/d*t*= −fα*I*, where f is the fraction of calcium ions that are not buffered, α = 1/*Z*V_{cyt}
*F*is a scaling factor that converts current to rate of change in [Ca^{2+}]_{i}, V_{cyt} is the volume of the cytosol,*F* is Faraday’s constant, and*Z* = 2 for*I*
_{Ca} and 1 for*I*
_{Na/Ca} (since the stoichiometry of this exchanger is 3 Na^{+}:1 Ca^{2+}). From *Eq. 1, a* and *c′*, it follows that for each of these currents*g* =*G*/fα =*Z*V_{cyt}
*FG*/f and that γ = V_{cyt}
*F*/fC_{m}. The formulation of *Eq. 1c* thus introduces a scaling factor, γ, that involves the parameters of V_{cyt}, C_{m}, and the fraction of calcium ions that remains free (note that the divalent nature of Ca^{2+}, i.e., *Z* = 2 is accounted for by the coefficient 2*G*
_{Ca} 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 (Table2). 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 Ca^{2+}-buffering proteins (36, 41). For the purposes of modeling, it therefore becomes necessary to define “effective volumes” available for Ca^{2+} movements (12) that are equivalent to 100 and 400 times the physical volumes of these compartments (
and
), respectively, whereas the ratio of their actual volumes (V_{SR}/V_{cyt}) is <1:4, because the SR occupies <20% of the cytosol. It follows that
∼
, which permits the assumption that Ca^{2+} concentration changes in the cytosol and stores resulting from intracellular movements between these two compartments are directly related according to
Equation 2a as is implicit in *Eq. 1*,*a* and*b*, if transmembrane fluxes are ignored. A more general relationship of the form
Equation 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 Ca^{2+} influx via L-type Ca^{2+} channels (11, 25). The flux through a single open VOCC has been estimated at ∼10^{5} ions/s, and there may be ∼1,000 such channels in a vascular myocyte, resulting in a flux of ∼10^{6} 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}cm^{2} 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 [Ca^{2+}]_{i}at a rate d[Ca^{2+}]_{i}/d*t*= 0.01[10^{6}/(0.2 × 1.2 × 10^{−12})][1/(6.24 × 10^{23})] ≈ 0.1 μM/s, where Avogadro’s number = 6.24 × 10^{23}mol^{−1}.

In the model the contribution of Ca^{2+} influx through VOCCs to d[Ca^{2+}]_{i}/d*t* is −*G*
_{Ca}(*z* −
)/{1 +
which reduces to 0.02*G*
_{Ca} μM/s for the fixed values*z*
_{Ca2}and *R*
_{Ca} used in the simulations at an “average” membrane potential of −40 mV (Fig. 2). In practice,*G*
_{Ca} was varied over the range 1–16 (Table 2), giving values of d[Ca^{2+}]_{i}/d*t*= 0.02–0.32 μM/s, which encompass the value expected from the above order-of-magnitude calculation. To reflect the dominance of Ca^{2+} influx via VOCCs, this flux was generally made up to four times larger (depending on the instantaneous value of membrane potential) than the combined Ca^{2+} flux entering the cytosol via InsP_{3}-induced release from stores and receptor-operated channels.

The magnitude of the Ca^{2+} 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 C_{m} of 1 μF/cm^{2} (43) and a surface area of ∼10^{−5}cm^{2}, then total C_{m} ≈ 10^{−5} μF. For a flux of*n* = 10^{6} ions/s, the current carried by Ca^{2+} is*I*
_{Ca} = 2*ne* = 3.2 × 10^{−13} C/s [where the electronic charge (*e*) = 1.6 × 10^{−19} C]. It follows that d*z*/d*t*=*I*
_{Ca}/C_{m}≈ 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 (seeresults). The coefficients chosen to reflect the magnitudes of net Ca^{2+}extrusion across the membrane by the Ca^{2+}-ATPase and Na^{+}/Ca^{2+}exchange and uptake by the SR also match relative estimates of these fluxes in vascular smooth muscle (see *Inhibition of the membrane Ca ^{2+}-ATPase*).

### Integration of Ca^{2+} Dynamics by Latch Bridge Formation

*Equation 1, a–d,* generates oscillations in [Ca^{2+}]_{i}, [Ca^{2+}]_{SR}, membrane potential, and the open state probability of K_{Ca} channels and does not simulate experimental fluctuations in arterial pressure and flow directly, although these will closely follow oscillations in [Ca^{2+}]_{i}. Force development and [Ca^{2+}]_{i}were therefore related by the cross-bridge phosphorylation and latch state model proposed by Hai and Murphy (29). In this model, elevated [Ca^{2+}]_{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*K*
_{1}–*K*
_{7}are the rate constants regulating the phosphorylation and bridge formation and can be written as the following set of differential equations
Equation 3a
Equation 3b
Equation 3c
Equation 3d
Equation 3eUnder 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*K*
_{1} and*K*
_{6}) are Ca^{2+} concentration dependent and assumed to follow the sigmoidal [Ca^{2+}]
/([
+
), which implies saturation at high [Ca^{2+}]_{i}. The values for the rate constants obtained by Hai and Murphy (29) for the swine carotid artery were*K*
_{1} =*K*
_{6} = 1.7 s^{−1},*K*
_{2} =*K*
_{5} = 0.5 s^{−1},*K*
_{3} = 0.4 s^{−1},*K*
_{4} = 0.1 s^{−1},*K*
_{7} = 0.02 s^{−1}, and*x*
_{0} = 0.6 μM. In practice, it becomes apparent that the key coefficient responsible for integrating Ca^{2+} fluctuations is*K*
_{7}, and in different simulations the value of this parameter was varied to obtain different degrees of dampening of the Ca^{2+} oscillations (Table 2). The latch equations are external to *Eq. 1, a–d*, and effectively integrate oscillations in [Ca^{2+}]_{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.2
*A* and in the time domain representation of Fig. 2
*B*. Along the slow section the dynamics are dominated by sequestration of Ca^{2+} entering the cell within the SR, with the result that [Ca^{2+}]_{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, Ca^{2+} is extruded from the cell more rapidly than it enters from the extracellular space, so that [Ca^{2+}]_{i}declines and [Ca^{2+}]_{SR}remains low. When [Ca^{2+}]_{i}becomes sufficiently low, uptake again dominates over extrusion, thus permitting sequestration of Ca^{2+}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 [Ca^{2+}]_{i}and [Ca^{2+}]_{SR}are zero (Fig. 2
*A*). Analytically, this equilibrium point is thus defined by (d*x*/d*t*)‖
= (d*y*/d*t*)‖
= 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
Equation 4a
Equation 4b where
= *x* −*x** and
= *y*− *y**. The Jacobian matrix of the system is The determinant (Det =*a*
_{11}
*a*
_{22}−*a*
_{12}
*a*
_{21}) and the trace (Tr =*a*
_{11} +*a*
_{22}) of this matrix characterize the nature of the phase space trajectories around the equilibrium point (56). Thus *1*) if Det > Tr^{2}/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 < Tr^{2}/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
/d*t*and d
/d*t*become zero. *Equation 4, a* and*b*, can then be written in the form
Equation 6a
Equation 6bFor the specific configuration of Fig.2
*A* 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
Equation 7a
Equation 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. 2
*A*. 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.2
*A* 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 Ca^{2+} movements across the membrane, so that [Ca^{2+}]_{i}+ [Ca^{2+}]_{SR}effectively remains constant. It follows that the slope of the very fast phase of the limit cycle in Fig.2
*A* 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.2
*C*. 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 K_{Ca} channel open state probability are zero. The equilibrium point (*z**,*w**) is defined by (d*z*/d*t*)‖
= (d*w*/d*t*)‖
= 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
Equation 8a
Equation 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. 2
*D*). The crucial mechanism that permits the emergence of oscillatory membrane activity is the intrinsic time-dependent inactivation of the K_{Ca} 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 K_{Ca} 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 γ*G*
_{K}, which governs the net contribution of K_{Ca} channel opening to d*z*/d*t*in *Eq. 1c* (Table 2).

The dynamics of the membrane oscillator are less robust than the dynamics generated by intracellular Ca^{2+} cycling, because the slope of the sigmoidal K_{Ca} channel activation curve and its position on the*z*-axis are highly sensitive to [Ca^{2+}]_{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., Ca^{2+}) dependence of*Eq. 1d
* at the equilibrium point can be understood by applying the transformation*z _{x}
* =

*R*

_{K}ln{β/[(

*x*+

*x*)

_{w}^{2}]} and expressing the

*w*nullcline in the form Equation 9It follows that changes in [Ca

^{2+}]

_{i}translate the inflection point

*z*

_{Ca3}of the K

_{Ca}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 [Ca^{2+}]_{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 *K*
_{7} (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*K*
_{7} (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*K*
_{7} correspond to rapid decay of the latch bridge, so that the system is capable of responding to rapid changes in [Ca^{2+}]_{i}. For low values of*K*
_{7} the latch decays slowly, allowing the system to function as an integrator of [Ca^{2+}]_{i}.

### 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.4
*A* the responses are dominated by the intracellular subsystem, inasmuch as all activity was abolished by ryanodine. In Fig. 4
*D* 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.

Low concentrations of ryanodine lock the Ca^{2+} 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 γ (γ*G*
_{Ca}, γ*G*
_{Na/Ca}, and γ*G*
_{K} 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 [Ca^{2+}]_{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 [Ca^{2+}]_{i}and [Ca^{2+}]_{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 [Ca^{2+}]_{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 (d*x*/d*t*)‖_{(x*,y*)}+ (d*y*/d*t*)‖_{(x*,y*)}= 0, so that, on adding *Eq. 1, a* and*b*, we obtain
Equation 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 Ca^{2+} by the ryanodine-sensitive store do not influence average [Ca^{2+}]_{i}.

#### Blockade of VOCCs.

In some experimental preparations the Ca^{2+} 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*G*
_{Ca}, which determines the magnitude of Ca^{2+}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 Ca^{2+} influx via VOCCs and other pathways, i.e., via ROCCs, Ca^{2+}release from the InsP_{3}-sensitive store, or reverse-mode Na^{+}/Ca^{2+}exchange (respectively represented by the coefficients*A*
_{0},*A*
_{1}, and*G*
_{Na/Ca}). When the net Ca^{2+} flux is small, blockade of VOCCs may prevent adequate replenishment of stores, so that vasomotion ceases. Conversely, if the net Ca^{2+} 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^{+}/Ca^{2+}exchanger is nearly fourfold greater and the contribution of*A*
_{0} +*A*
_{1} ∼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 [Ca^{2+}]_{SR}and lower [Ca^{2+}]_{i}with two quite distinct effects: in Fig.7
*A* 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 *G*
_{Ca} = 0 (Fig. 7
*B*). Movement of the equilibrium point upward and to the left is accompanied by an increase in the amplitude of oscillatory behavior as mean [Ca^{2+}]_{i}and force development decrease, consistent with the experimental trace shown in Fig. 6
*D*. 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*.

#### Activation by constrictor agonists.

Figure 8
*A*shows a typical experimental response in which irregular rhythmic activity results from the administration of histamine in an otherwise quiescent artery. As shown in Fig. 8
*B*for a different preparation, further increments in histamine concentration may nevertheless result in overstimulation and the suppression of oscillatory behavior at high levels of perfusion pressure. In the former case the pharmacological induction of vasomotion by histamine was simulated simply by increasing the combined contribution of*A*
_{0} +*A*
_{1} (Fig.9
*A*) and is readily interpreted by a phase plane plot of *Eq. 1, a* and *b* (Fig.9
*B*). For low values of*A*
_{0} +*A*
_{1}, there is subcritical Ca^{2+} flux into the cytosol, which cannot sustain an adequate rate of Ca^{2+} sequestration in the ryanodine-sensitive component of the SR and thus oscillatory behavior (cf. initial part of Fig. 8
*A*). This situation is denoted by α in Fig.9
*B*, which represents a stable node. Increases in *A*
_{0} +*A*
_{1} move the equilibrium point through a short segment, where it is a stable focus, into an active regimen, where it represents an unstable focus that permits oscillatory behavior (β and γ in Fig.9
*B*).

Additional increases in*A*
_{0} +*A*
_{1} cause further movement of the intersection point to the right, resulting in overstimulation and loss of oscillatory behavior at a stable focus (δ in Fig. 9
*B*). This explains how high concentrations of histamine can stop the contribution of the slow intracellular subsystem to vasomotion experimentally (cf. Fig.8
*B*). Two simulations of this scenario are presented in Fig. 10. In Fig. 10
*A* the action of histamine is simulated by increasing*A*
_{0} +*A*
_{1}, whereas in Fig. 10
*B* the coefficient*G*
_{Ca} is simultaneously increased to represent enhanced Ca^{2+} influx through VOCCs. Although both simulations possess similar characteristics, the*z-w* projection of the four-dimensional attractor of the system reveals an important difference. In Fig.10
*C*, increasing*A*
_{0} +*A*
_{1} alone does not significantly influence average membrane potential, whereas in Fig.10
*D* the additional increase in the coefficient *G*
_{Ca}induces cell depolarization, in agreement with the known experimental effects of constrictor agonists (51).

#### Inhibition of the membrane Ca^{2+}-ATPase.

Two cases that illustrate the effect of the vanadate ion, which blocks the membrane Ca^{2+} extrusion ATPase, are presented to differentiate between Ca^{2+} extrusion by the membrane Ca^{2+}-ATPase and by the Na^{+}/Ca^{2+}exchanger. Figure 11,*A* and*B*, shows experimental and simulated traces exhibiting a prominent slow component in which vanadate caused a marked rise in perfusion pressure but only minor effects on the superficial form of the signals until oscillatory activity was abolished. By contrast, in the experimental and simulated signals of Fig. 11, *D* and*E*, fast and slow oscillations coexist and vasomotion persisted when Ca^{2+}-ATPase inhibition with vanadate caused a rise in pressure, although the oscillations became simpler in form. In Fig. 11
*B*, extrusion was dominated by the Ca^{2+}-ATPase, whereas Na^{+}/Ca^{2+}exchange was the principal extrusion mechanism in Fig.11
*E* (Table 2). The phase plane analysis of Fig. 11, *C* and*F*, confirms the presence of the typical triangular *x-y* projection of the attractor attributable to CICR whether extrusion is dominated by either mechanism. The simulation of Fig.11
*E* demonstrates how the cytosolic and membrane subsystems can become entrained if Na^{+}/Ca^{2+}exchange becomes the sole mechanism for Ca^{2+} extrusion, resulting in relatively simple patterns of dynamic behavior.

In the case of rabbit vena cava, Nazer and van Breemen (49) estimated experimentally that the ratio of the Ca^{2+} fluxes exiting the smooth muscle cell via the Ca^{2+}-ATPase to the Na^{+}/Ca^{2+}exchanger is on the order of 2:1 under steady-state conditions. The values of *G*
_{Na/Ca}and *D* chosen for the simulations in Fig. 11
*E* approximately reflect this ratio, with integration of Ca^{2+}movements via the Ca^{2+}-ATPase and the Na^{+}/Ca^{2+}exchanger over a large number of oscillatory cycles for the initial stationary section of the dynamics, giving average outward fluxes of 0.65 and 0.25 μM/s, respectively. The present results suggest that the Ca^{2+}-ATPase or forward-mode Na^{+}/Ca^{2+}exchange are able to extrude Ca^{2+}from the cell at a sufficient rate to sustain the dynamics of the CICR oscillator.

#### K_{Ca} channel inhibition.

Figure12
*A*shows an original pressure trace from an experiment in which TEA, which blocks a spectrum of K_{Ca} channels, promoted constriction and suppressed membrane oscillations and decreased the amplitude of slow activity attributable to the intracellular oscillator. Model-generated oscillations in force development simulating the effects of TEA are shown in Fig.12
*B*, with a phase space interpretation in Fig. 12
*C*. Reductions in the coefficient *S* in *Eq.1d
* depress the *w*nullcline and cause the equilibrium point to move to the right, consistent with membrane depolarization, enhanced Ca^{2+} influx, and a rise in mean [Ca^{2+}]_{i}and force development, as observed experimentally. This increase in [Ca^{2+}]_{i}interactively reduces the maximum in the*z* nullcline, translates the*w* nullcline to the left, and decreases the active range over which oscillations are possible. The equilibrium point consequently moves to a region of progressively smaller-amplitude voltage oscillations, until they are ultimately abolished at a stable focus or node under conditions of relative depolarization and high levels of constriction. However, until oscillations stop, the period of the limit cycle remains almost constant. Inasmuch as K_{Ca} channel blockers do not significantly alter the frequency of the membrane subsystem experimentally (14) (Fig. 12
*A*), this justifies modeling the effects of TEA by reducing the activation-to-inactivation ratio (*S*) rather than the scaling factor (λ). Also the model predicts that large values of *S* (i.e., high K_{Ca} channel activity) may suppress oscillatory activity by converting the equilibrium point to a stable focus or node (Fig. 12
*C*). Experimentally, the change in average membrane potential that follows administration of high concentrations of K_{Ca} channel blockers such as TEA may vary between ∼10 and ∼25 mV, which is consistent with the range over which the equilibrium point shifted (∼20 mV) after reductions in*S* in Fig.12
*C* (19, 70).

#### Role of Cl^{−} channels.

Figure13
*A*shows an experimental pressure trace before and after administration of 30 μM niflumic acid, which inhibits Cl^{−} channel activity. This demonstrates selective loss of the membrane oscillator, with the slow intracellular component relatively unaffected. At the concentration of niflumic acid administered, there was a small fall in perfusion pressure. Figure 13
*B* shows model-generated tension oscillations that closely simulate these experimental findings on reducing the term*G*
_{Cl} in*Eq. 1c*. A phase space plot of*Eq. 1, c* and*d*, shows that oscillations are abolished when the equilibrium point of *Eq. 1, c* and *d*, becomes a stable focus. In rat cerebral arteries, administration of Cl^{−} channel blockers has been shown to promote hyperpolarization of 10–15 mV in vessels pressurized to induce depolarization and myogenic tone but to have no effect on membrane potential if intravascular pressure is low (50). In the simulation of Fig. 13 the value of*G*
_{Cl} was selected such that the contribution of Cl^{−} fluxes to net membrane potential was on the order of ∼5 mV.

#### Na^{+}-K^{+}-ATPase inhibition.

Figure14
*A*shows an experimental pressure trace after administration of ouabain, which inhibits the membrane Na^{+}-K^{+}-ATPase and preferentially attenuates the membrane oscillator while it also reduces the amplitude of the intracellular oscillator. The action of ouabain was simulated by reducing the term*F*
_{Na/K} in*Eq. 1c* and successfully reproduced the general form of the experimental responses (Fig.14
*B*). Phase plane analysis shows that ouabain causes the *z* nullcline to translate vertically on the ordinate, such that the equilibrium point moves to the right in association with membrane depolarization (Fig.14
*C*). This in turn enhances Ca^{2+} influx and elevates [Ca^{2+}]_{i}in association with a change in the nullcline configuration, such that the active range over which oscillations are possible is reduced; i.e., there is a translation of the *w*nullcline to the left. As a consequence of the combined nullcline movement, membrane oscillations become progressively smaller in amplitude until they are ultimately abolished at a stable focus or node. Experimentally, blockade of the Na^{+}-K^{+}-ATPase with high concentrations of ouabain may result in arterial depolarization of up to ∼10 mV (1), which is comparable to the change in membrane potential observed in the model simulation on reducing the value of the coefficient*F*
_{Na/K} selected to zero.

#### Attenuated Na^{+}/Ca^{2+}exchange.

Figures15
*A* and16
*A* show experimental pressure traces during graded reductions in [Na^{+}]_{o}and illustrate how low [Na^{+}]_{o}, which depresses Na^{+}/Ca^{2+}exchange, may abolish vasomotion or selectively attenuate membrane activity while it amplifies the contribution of the intracellular oscillator. Such apparently paradoxical behavior can be simulated by reducing *G*
_{Na/Ca}in *Eq. 1, a* and*c*, with the existence of two distinct patterns of response explained by the ability of Na^{+}/Ca^{2+}exchange to mediate net efflux or net influx of Ca^{2+} across the cell membrane, depending on whether average membrane potential is below or above the reversal potential of the exchanger. At high levels of depolarization, VOCCs and reverse-mode Na^{+}/Ca^{2+}exchange mediate Ca^{2+} influx into the cytosol, thus elevating *x** and pressure; when the cell hyperpolarizes, the direction of Na^{+}/Ca^{2+}exchange allows forward-mode efflux of Ca^{2+}.

In the simulation of Fig. 15, the contribution of voltage-dependent Ca^{2+} influx is minimized by setting *G*
_{Ca} = 0 to isolate the effects of Na^{+}/Ca^{2+}exchange. In this example, membrane potential is always below −35 mV (Fig. 15
*F*), and the exchanger operates in forward mode, thereby contributing to a relatively low [Ca^{2+}]_{i}of ∼0.3 μM (Fig. 15
*C*). Subsequent blockade of the exchanger (i.e., a decrease in*G*
_{Na/Ca}) causes a marked rise in [Ca^{2+}]_{i}due to reduced extrusion and thus enhances force development. As shown in the nullcline analysis of Fig. 15
*C*, the rise in [Ca^{2+}]_{i}may be sufficient to suppress the intracellular oscillator by converting an unstable equilibrium point to a stable focus or node. In this example, the existence of a stable focus is confirmed by the*x-y* projection of the attractor of the system (Fig. 15
*D*). Inhibition of the Na^{+}/Ca^{2+}exchanger simultaneously suppresses the fast component of vasomotion by altering the position of the *z*nullcline of the membrane subsystem in such a fashion as to convert an unstable focus to a stable focus and subsequently a stable node (Fig.15, *E* and*F*). The model predicts that the cell membrane should hyperpolarize in this scenario as a direct consequence of the 3:1 stoichiometry of Na^{+}/Ca^{2+}exchange.

Figure 16 illustrates the contrasting scenario in which Ca^{2+} influx via voltage-operated channels is large, membrane potential is generally more positive than the reversal potential for Na^{+}/Ca^{2+}exchange, and reverse-mode operation promotes Ca^{2+} influx into the cytosol and contributes to a higher average [Ca^{2+}]_{i}of ∼0.4 μM (Fig. 16
*D*). Inhibition of the Na^{+}/Ca^{2+}exchanger consequently lowers [Ca^{2+}]_{i}and induces large-amplitude oscillations (Fig.16
*C*). This is analogous to the effects of verapamil described above. At low levels of [Na^{+}]_{o}, the *x-y* phase plane projection of the underlying attractor exhibits comparatively high levels of instantaneous [Ca^{2+}]_{i}, although overall there is nevertheless a reduction in*x** and consequently net force development (Fig. 16, *B* and*D*). Such behavior can be appreciated on inspection of the experimental and corresponding simulated time series. In Fig. 16, *A* and*B*, lowering [Na^{+}]_{o}does not suppress the membrane oscillator continuously, and the contribution of chaotic *z-w* limit cycle is reflected as small-amplitude oscillations found in the troughs between peaks in force. It is thus apparent that the increase in the size of the limit cycle of the intracellular oscillator has an important modulatory effect on membrane activity. High levels of instantaneous [Ca^{2+}]_{i}move the *w* nullcline outside the range that permits oscillatory activity, thus suppressing the membrane subsystem, the trajectories of which at low levels of [Na^{+}]_{o}markedly decelerate when they enter the region of the box shown in Fig.16
*F*. This phase corresponds to the peak in force development during the vasomotion cycle; oscillations subsequently become reestablished when instantaneous [Ca^{2+}]_{i}returns to low values. Changes in*G*
_{Na/Ca} affect the frequency of the intracellular oscillator according to the associated change in average [Ca^{2+}]_{i}, with a small increase apparent in the transients of Fig. 15,*A* and*B*, and a small decrease apparent in Fig. 16, *A* and*B*.

### Universal Routes to Chaos

#### Period doubling and integral periodicity.

In the Feigenbaum period-doubling “route to chaos,” progressive increases in a single control parameter cause a steady state to bifurcate first into two alternating solutions, thereby generating*period 2* dynamics; subsequent bifurcations produce *period 4, 8*, and*16* oscillations, and so on until the period-doubling cascade ultimately leads to a region of fully developed chaos. Figure 4
*D* illustrates a period-doubling bifurcation in an artery in which ryanodine converted chaos into periodic membrane behavior, and Fig. 4,*E* and*F*, shows that this scenario can be reproduced in simulations. Figure 17,*A* and*B*, illustrates experimental and model-generated examples of *period 8*behavior, which is again consistent with period doubling. In nonlinear systems, narrow windows of exactly integral periodic behavior, which may be odd in nature, are found within the chaotic regimen, with*period 3* (Fig. 17,*C* and*D*) and *period 5* dynamics (Fig. 17, *E*and *F*) typically being the patterns most commonly observed experimentally in isolated rabbit ear arteries and in simulations with the present model.

#### Quasiperiodicity and mixed-mode behavior.

The intracellular and the membrane subsystems can generate quasiperiodic patterns that typically consist of two distinct oscillations possessing an incommensurate frequency ratio (11). Experimental and simulated examples are illustrated in Fig.18, *A*and *B*. Related patterns known as mixed-mode responses, which in their simplest form represent the frequency-locked states of coupled oscillators that resonate when the ratio of their frequencies is a rational fraction, may also be observed (11). These may be classified according to the notation*k ^{m}
*, where

*k*large excursions are followed by

*m*smaller ones. Figure 18,

*C*and

*D*, provides experimental and simulated examples consisting of a mixture of oscillations of 1

^{m}type.

#### Intermittency.

Another recognized route to chaos, which involves the gradual conversion of a periodic signal to sustained chaos under the variation of a single system parameter, is known as intermittency, in which periodic behavior is interrupted by bursts of irregular behavior of increasing duration (24, 28). Figure 18,*E* and*F*, illustrates experimental signals in the vicinity of a *period 3* window in which the periodic signal is interrupted by chaos. This pattern is typical of so-called type I intermittency, which is associated with the destabilization of an equilibrium point on the first-return map constructed from the successive maxima present in the signal (31). In phase space the trajectories of nearly periodic oscillations gradually move away from this equilibrium point and enter a chaotic region but are eventually reinjected into the neighborhood of the equilibrium point, thereby generating episodes of nearly periodic behavior of varying duration.

## DISCUSSION

Correlation dimension analysis suggests that a model with at least four independent control variables is required to reproduce the dynamics of vasomotion in isolated rabbit ear arteries. The present system of equations achieves this degree of complexity by assigning two independent control variables to each of two autonomous, but nevertheless coupled, nonlinear subsystems. Variations in the parameters that determine the magnitudes of the ion fluxes included in the model permit realistic simulations of a range of pharmacological interventions, including apparently paradoxical responses, and allow the emergence of hallmark patterns of nonlinear behavior. Such versatility could not be obtained in a linear system in which the contributions of different control mechanisms were simply summated.

### Intracellular Oscillator

The intracellular subsystem can be considered an extension of the two-pool CICR model developed by Goldbeter et al. (20) to describe [Ca^{2+}]_{i}spiking in nonexcitable cells that additionally incorporates Ca^{2+} influx via L-type channels and Na^{+}/Ca^{2+}exchange. These transport systems are modulated by membrane potential and are important determinants of excitation-contraction coupling in vascular smooth muscle. The contributions of receptor-operated Ca^{2+} influx and IICR were represented as constant fluxes (*A*
_{0} and*A*
_{1}), the magnitudes of which were assumed to parallel the concentration of activating agonist and, from a mathematical viewpoint, were interchangeable.

The model provides a systematic explanation for differences between interventions that modulate CICR and other transport systems, inasmuch as the coefficients of the terms contributing to the intracellular oscillator can be classified into two principal subgroups:*1*) those common to*Eq. 1, a* and*b* (i.e., *B, C*, and *L*), and*2*) those that appear in*Eq. 1a* but not in *Eq. 1b*(*A*
_{0},*A*
_{1},*G*
_{Ca},*G*
_{Na/Ca}, and*D*). Changes in the first set of coefficients preserve the relative position of the*x* and*y* nullclines, so that alterations in CICR (e.g., coefficients *C* and*L*) do not modulate the amplitude and frequency of the CICR limit cycle or average [Ca^{2+}]_{i}, which reflects ambient force development, until oscillations stop relatively abruptly. This closely matches the pharmacological effects of ryanodine observed experimentally. By contrast, interventions that cause changes in the second set of coefficients affect the position of the *x* but not the*y* nullcline. The hyperbolic shapes of the nullclines over the range that permits oscillatory activity then maintain an inverse relationship between average [Ca^{2+}]_{i}and the duration of the slow SR refilling phase, which is the principal determinant of the amplitude and frequency of the intracellular limit cycle (Fig. 2
*A*). This analysis suggests that interventions that increase tone will be associated with decreases in the amplitude and period of intracellular oscillations, whereas dilator stimuli that act by reducing [Ca^{2+}]_{i}will cause the opposite effect. These predictions were supported by experiments with the Ca^{2+}-channel antagonist verapamil and low [Na^{+}]_{o}, which directly influence transmembrane Ca^{2+} fluxes, and agents such as TEA and ouabain, which promote voltage-dependent Ca^{2+} influx secondary to membrane depolarization. The analysis also predicts that a small fall in net Ca^{2+} influx may allow oscillations to persist with increased amplitude, whereas after more severe reductions the system may no longer be able to sustain the CICR mechanism. This explains the existence of apparently paradoxical experimental responses to verapamil and low [Na^{+}]_{o}, which can promote the emergence of slow large-amplitude oscillations or abolish vasomotion completely. Similar arguments explain how histamine stimulates vasomotion at low concentrations, whereas high concentrations may subsequently result in the loss of rhythmic behavior. The model predicts that these experimental observations can be reproduced simply by increasing the combined flux of Ca^{2+} represented by the terms*A*
_{0} and*A*
_{1}, even though this approach to the action of histamine did not simulate the net depolarization of vascular smooth muscle that generally accompanies administration of constrictor agonists. To obtain this effect of agonists on membrane potential, it was necessary to increase the coefficient that determines voltage-dependent Ca^{2+} influx,*G*
_{Ca}, in parallel with *A*
_{0} and*A*
_{1}, although this strategy did not greatly affect the form of the oscillations obtained.

### Membrane Oscillator

The formulation of the membrane subsystem in terms of cell potential and K_{Ca} channel open state probability is mathematically analogous to the FitzHugh-Nagumo model of action potential generation, in which depolarization increases the open state probability of an opposing hyperpolarizing ion channel, the spontaneous inactivation of which is governed by simple exponential decay (59). Oscillatory behavior is thus dependent on the lag between the sequential activation and inactivation of voltage-dependent Ca^{2+} influx and reflects the concerted action of all ion transport mechanisms that influence membrane potential, i.e., the terms represented by the coefficients*F*
_{Na/K},*G*
_{Cl},*G*
_{Ca},*G*
_{Na/Ca}, and*G*
_{K} in*Eq. 1c*. Experimentally, selective inhibition of any of these transport systems results in loss of the fast component of vasomotion, and the model shows that the balance between the contributing ion fluxes is critical if oscillations are to be sustained, because the active range of the*w* nullcline is short in length and highly sensitive to [Ca^{2+}]_{i}. Table 2 thus shows that, in some simulations, oscillatory K^{+}- channel dynamics could be supported by voltage-operated Ca^{2+}influx and Na^{+}/Ca^{2+}exchange alone, with fluxes due to the Na^{+}-K^{+}-ATPase and chloride ions (*F*
_{Na/Ca} and*G*
_{Cl}) being set to zero, whereas in other simulations, Na^{+}-K^{+}-ATPase and/or Cl^{−} channel activity was necessary to sustain oscillatory behavior. Phase plane analysis of the model can also explain how membrane activity can cease when the K_{Ca} channel open state probability is low or high (Fig. 12), potentially accounting for “paradoxical” experimental observations that K_{Ca}-channel blockers promote oscillatory behavior in quiescent arteries but suppress vasomotion in active preparations (19, 27, 70).

Although a relatively wide spectrum of pharmacological interventions is able to abolish membrane activity in rabbit arteries, in general the frequency of the membrane oscillator remains at ∼0.06 Hz until its contribution is no longer detectable in experimental signals (14). This finding can also be explained by phase plane analysis of the model. In contrast to the hyperbolic shapes of the*x* and*y* nullclines of the intracellular oscillator, which generate an asymmetric limit cycle, the size of which determines frequency and amplitude, the active ranges of the*z* and*w* nullclines are effectively straight-line segments with an almost circular limit cycle, the period of which changes little during simulations even when its size varies substantially. Dynamically, this is analogous to the behavior of a simple harmonic oscillator, such as a pendulum, the limit cycle of which is also circular with a period independent of oscillation amplitude.

### Routes to Chaos and the Role of Coupling

Nonlinearity in the interaction between two coupled oscillators may allow the emergence of “universal” patterns of dynamics that have been extensively characterized in the mathematical, physical, and chemical literature and can be observed experimentally in the responses of rabbit ear arteries. These include the Feigenbaum period-doubling cascade, which is almost ubiquitous in nonlinear systems, and the intermittent route to chaos. Pomeau and Manneville described three intermittency classes, each of which corresponds uniquely to the manner in which the eigenvalues of a fixed point pass through the unit circle at a local bifurcation. In the type III scenario, which has been observed in the dynamics of vasomotion in rabbit ear arteries, intermittency arises through the destabilization of a limit cycle via a subcritical period-doubling bifurcation (28). In the present study we have provided an experimental trace demonstrating the characteristics of type I intermittency, which arises at a tangent bifurcation (31), and have shown that this can be simulated by the model.

Vasomotion in rabbit ear arteries also displays the so-called quasiperiodic route to chaos, in which the strength of the coupling between two oscillators and the ratio of their frequencies influence the nature of the dynamics of the coupled system (11). If the frequency ratio is commensurate, there is a frequency-locked resonance and the trajectories of the system trace out a closed orbit on the surface of a three-dimensional torus that is exactly periodic. By contrast, if the frequency ratio of the oscillators is irrational, the trajectories of their combined motion ultimately cover the surface of a torus completely, because no orbit ever repeats twice, and the dynamics are termed quasiperiodic. Chaotic orbits arise if the motion on the torus becomes unstable. Weak coupling between the two oscillators favors the emergence of quasiperiodicity or frequency-locked states, which in the case of vasomotion can be identified as mixed-mode responses (11). Increasing levels of coupling permit the appearance of chaos but may ultimately cause the oscillators to become entrained, resulting in relatively simple patterns of behavior.

All these possibilities were evident during simulations that manipulated the coefficients of the transmembrane ion fluxes that are common to *Eq. 1, a* and*c*, and therefore regulate the coupling between the two oscillators. Thus a transition from chaos to quasiperiodicity became evident on reducing terms that appear in both equations (e.g.,*G*
_{Ca} in Fig.6
*D* or*G*
_{Na/Ca} in Fig.16
*B*), with the waveforms corresponding to the fast and the slow components then becoming almost linearly superimposed. On the other hand, chaos may also be transformed to periodic behavior if coupling increases, e.g., if there is only a single mechanism (Na^{+}/Ca^{2+}exchange) mediating extrusion of Ca^{2+} from the cell after inhibition of Ca^{2+}-ATPase by vanadate (Fig. 11
*E*). In this situation the two oscillators become entrained, because Na^{+}/Ca^{2+}exchange is common to *Eq. 1, a* and*c*. Analogously, simulations in which efflux via the Ca^{2+}-ATPase was assumed to be electrogenic, and thus the full term containing the coefficient *D* in *Eq. 1a* included in *Eq. 1c*, resulted in patterns of vasomotion that were often simpler in form than the electroneutral case, inasmuch as increased coupling between the two subsystems then reduced the geometric complexity of the underlying attractors (not shown).

### Compartmentalized Function: One-Pool Models and the Role of Diffusion

Calcium ions entering the vascular smooth muscle cell are thought to be preferentially sequestered within the peripheral component of the SR before reaching the myofilaments [the “superficial buffer barrier” hypothesis (38, 49)]. The present assumption that Ca^{2+} uptake into an InsP_{3}-sensitive pool was sufficiently rapid to keep this store permanently full and subsequently permit efflux of Ca^{2+} ions into the cytosol (at a constant rate for a given agonist concentration) may be considered an idealization of this buffer role but is necessarily a simplification, since the open state probability of the InsP_{3}-sensitive Ca^{2+} channel is a bell-shaped function of Ca^{2+} concentration (8). A single-pool model has thus shown how steady InsP_{3} levels may induce cyclic emptying and refilling of the InsP_{3}-sensitive store if there is fast activation of the InsP_{3}receptor at low [Ca^{2+}]_{i}followed by slow inactivation at high [Ca^{2+}]_{i}(41, 42). It is also possible to construct an oscillator based on a single pool of intracellular Ca^{2+}if depletion of this compartment by agonists or SR Ca^{2+}-ATPase inhibitors promotes capacitative Ca^{2+} influx (13, 43). Indeed, SR ATPase inhibitors trigger ryanodine-sensitive oscillations in [Ca^{2+}]_{i}in lymphocytes that have been suggested to depend on the time delay inherent to depletion of stores and secondary influx (13). The Ca^{2+} sensitivity of the InsP_{3} receptor was not incorporated in the present minimal model, inasmuch as the mechanism theoretically sustains ryanodine-insensitive [Ca^{2+}]_{i}oscillations without requiring Ca^{2+} influx, which conflicts with the pharmacological characteristics of vasomotion observed in rabbit ear arteries experimentally (11, 41, 42). Capacitative Ca^{2+} influx was not included in the model, inasmuch as previous studies have failed to demonstrate evidence for its involvement in the regulation of tone in this artery type (11, 27).

The different microenvironments of the peripheral and central SR result in significant gradients in Ca^{2+}concentration between the subplasmalemmal space and the bulk of the cytosol, so that localized elevations in near-membrane Ca^{2+} concentration modulate ion transport systems in a spatially regulated fashion (38). “Sparks” of high Ca^{2+} concentration generated by spontaneous focal release from the peripheral SR may thus selectively activate K_{Ca} channels and promote hyperpolarization and relaxation (38). There is also a spatially regulated interaction between the SR and the cell membrane that allows calcium ions to be released vectorially from the peripheral SR into the subplasmalemmal space and extruded from the cell via the Na^{+}/Ca^{2+}exchanger, which is located within closely adjacent caveolar membrane domains (38). Inasmuch as the exchanger is thought to be colocalized with the Na^{+}-K^{+}-ATPase (38), high levels of Ca^{2+}concentration in the subplasmalemmal space may stimulate exchange of Ca^{2+} for Na^{+}, elevate intracellular Na^{+} concentration, and thereby increase the activity of the Na^{+}-K^{+}-ATPase, which is not saturated at normal levels of intracellular Na^{+} concentration (38). The scenario could contribute to oscillatory release of Na^{+} and K^{+} from isolated arteries during rhythmic activity (58) and promote hyperpolarization that exerts a negative feedback on Ca^{2+} influx.

Explicit representation of such sources of spatial heterogeneity and compartmentalized function would require the introduction of not only further ionic species (such as Na^{+}and K^{+}) as independent dynamic variables but also second-order partial differential equations to model diffusion of Ca^{2+} and/or second messengers within the intracellular domain. Indeed, waves of elevated [Ca^{2+}]_{i}are known to propagate through the cytoplasm of vascular smooth muscle cells after activation by agonists (33, 61). Although the present minimal model readily reproduced the characteristics of vasomotion in isolated rabbit ear arteries, it is clear that additional factors would confer greater dynamical flexibility. Preliminary simulations, incorporating a further constant term in *Eq. 1c* to simulate the depolarizing inward Na^{+} current that is thought to contribute to the activation of vascular smooth muscle by agonists (51), demonstrate a significant expansion in the range of oscillatory behavior possible and greater flexibility in the selection of the coefficients representing other membrane ion transport systems, i.e.,*G*
_{Ca},*G*
_{Na/Ca},*G*
_{Cl}, and*F*
_{Na/K}. The assumption that changes in [Ca^{2+}]_{i}and [Ca^{2+}]_{SR}are directly related is also a simplification, and preliminary simulations with a more general relationship in which the ratio of the “effective” volumes of the cytosol and SR,
/
, was different from 1 (*Eq. 2*) generated patterns of vasomotion with increased complexity as a result of the asymmetry introduced between *Eq. 1a* and *Eq. 1b*.

### Mechanical Models

During externally superimposed changes in intraluminal pressure, vasomotion exhibits an inverse nonlinear relationship between oscillation frequency and amplitude (24, 63, 68). Furthermore, changes in intraluminal flow rate, which influence NO production through a shear stress-dependent mechanism, exert pronounced effects on the patterns of vasomotion observed (26, 63). Nonlinear analysis of signals from rabbit and rat arteries has shown, however, that the correlation dimension of vasomotion, D_{C}, is insensitive to changes in pressure, flow, and NO activity (24, 26, 63), supporting the view that mechanical forces do not influence the underlying cellular mechanisms in a fundamental dynamical sense, as is perhaps to be expected. For example, active constriction in response to increased transmural pressure (the myogenic response) involves ion fluxes similar to those already incorporated in the model, including Ca^{2+} release from internal stores, Ca^{2+} influx through stretch-operated cation channels, and depolarization-induced activation of VOCCs (45, 68). Increases in transmural pressure may also promote depolarization and myogenic behavior by causing the closure of K_{Ca} channels (70). Further modeling studies should allow expansion of the present formulation of vasomotion to accommodate the effects of mechanical forces, thus providing a link with models that directly relate active and passive wall stress to [Ca^{2+}]_{i}and models based on the steady-state and rate-dependent components of the myogenic response, which are capable of generating oscillatory behavior (2, 67).

### Conclusions

Previous two- and three-dimensional models have provided insights into the mechanisms that generate oscillations in [Ca^{2+}]_{i}but cannot reproduce the experimental complexity of vasomotion in isolated rabbit ear arteries. Conceptually, the present approach represents a significant advance by coupling two essentially independent nonlinear oscillatory subsystems, which allows the response to pharmacological intervention to be understood at the cellular level. The ion fluxes included in the model are common to many cell types, so that variations in the coefficients that determine their relative magnitudes may form a paradigm to explain differences between cell types within a single nonlinear framework. In the context of the vascular system, analysis and modeling of smooth muscle control mechanisms should ultimately provide new insights into the dynamic coordination of organ perfusion.

## Footnotes

Address for reprint requests and other correspondence: T. M. Griffith, Dept. of Diagnostic Radiology, Univ. of Wales Coll. of Med., Heath Park, Cardiff CF4 4XN, UK (E-mail:Griffith{at}Cardiff.AC.UK).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

- Copyright © 1999 the American Physiological Society