Abstract
Nitric oxide (NO) produced by the vascular endothelium is an important biologic messenger that regulates vessel tone and permeability and inhibits platelet adhesion and aggregation. NO exerts its control of vessel tone by interacting with guanylyl cyclase in the vascular smooth muscle to initiate a series of reactions that lead to vessel dilation. Previous efforts to investigate this interaction by mathematical modeling of NO diffusion and reaction have been hampered by the lack of information on the production and degradation rate of NO. We use a mathematical model and previously published experimental data to estimate the rate of NO production, 6.8 × 10^{−14}μmol ⋅ μm^{−2} ⋅ s^{−1}; the NO diffusion coefficient, 3,300 μm^{2}s^{−1}; and the NO consumption rate coefficient in the vascular smooth muscle, 0.01 s^{−1} (1storder rate expression) or 0.05 μM^{−1} ⋅ s^{−1} (2ndorder rate expression). The modeling approach is discussed in detail. It provides a general framework for modeling the NO produced from the endothelium and for estimating relevant physical parameters.
 mass transfer
 kinetics
 parameter estimation
 diffusion
nitric oxide (NO) plays a diverse role as a biologic messenger, regulator, and cytotoxic agent. In the circulatory system, NO produced by the endothelium is a major factor in preventing platelet aggregation and in regulating microvascular tone and permeability (18). This control is exercised through its interaction with soluble guanylyl cyclase in the vascular tissue. However, the mechanism by which NO, a reactive free radical, is transferred from the endothelium to the lumen and vascular smooth muscle is still under debate.
Mathematical modeling can be a useful tool for investigating the simultaneous diffusion and reaction of NO. However, these models must be based on known mechanisms and employ physiologically relevant rate constants. These rate constants are not well characterized, nor is it clear how the NO production and decomposition rates affect the diffusion distance of NO in the blood, the endothelium, and the surrounding smooth muscle.
In vivo, NO diffuses from the endothelium into the luminal and abluminal regions. The NO that diffuses into the abluminal region travels into the vascular smooth muscle and binds with guanylyl cyclase. Because of its reactivity, NO may be consumed by numerous reactions before it reaches its target. In particular, NO may react with O_{2} or free radicals (19), bind with hemecontaining proteins (20), or interact with ironsulfur centers (18), as well as participate in several other reactions that result in its degradation (1). Because of the complexity of the kinetics, the interaction of NO in the smooth muscle is not well characterized and is usually ignored or treated as a firstorder reaction with a halflife (t _{½}). The value of t _{½} is usually estimated from the decomposition of NO in aqueous solution, where the reaction of NO with O_{2} is second order in NO (7,17, 27), or from measurements made from perfused tissue (9). However,t _{½} alone provides little insight when diffusion and reaction occur simultaneously.
Several previous models of NO reaction and diffusion depended on the in situ NO concentration measurements obtained by Malinski et al. (16), who used an electrochemical technique to measure the NO concentration at the endothelium and in the vascular smooth muscle of a rabbit aorta. Although the primary goal of their work was to demonstrate the utility of their porphyrinic microsensor (15), their data have had significant impact on NO diffusion modeling efforts.
The first mathematical model of NO diffusion and reaction in a biologic system was that of Lancaster (1113). He investigated whether the NO concentration in a cell was primarily influenced by the cell or the surrounding cells and examined the influences of NO sinks. He modeled the endothelium as a line of NO point sources and sinks and used an NO production rate based on the maximum NO concentration from the measurements of Malinski et al. (16). Lancaster’s initial work was closely followed by that of Wood and Garthwaite (28), who modeled neuronal NO production and diffusion in the brain. They modeled a neuron as a point source of NO and assumed that NO degraded in the brain tissue as a firstorder reaction witht _{½} of 0.5–5 s. The data of Malinski et al. were used as a guide in approximating the NO production rate. These researchers used experimental data to justify their choices for the NO production rate, but the accuracy of the estimates and the effects of other parameters were not examined. Furthermore, the finite size of the NOemitting region has not been taken into account.
The NO concentration in tissue depends on the diffusion coefficient and the degradation (reaction) rate as well as the production rate. Our purpose here is to provide a general procedure for estimating these three parameters, to determine the error in these estimates, and to investigate the relationships between the parameters. To accomplish this, we set out a general mathematical model for NO production from the endothelium and consumption in the surrounding tissue. This model is then applied to the experiment of Malinski et al. (16). The parameters and their variability are determined by applying nonlinear methods to fit the model to the data.
ESTIMATING PARAMETERS FROM EXPERIMENTAL DATA
System description.
We modeled NO production, diffusion, and reaction as a threesection system composed of the luminal region, the endothelium, and the abluminal region. NO is produced from the surfaces of the endothelium and diffuses into each of these regions. The spatial variation of NO concentration will be different in each section, reflecting the diffusion coefficient and the reaction rate expression of that section.
To obtain realistic model predictions, the size of the NOproducing region, the NO production rate, the rate constants for the NO reaction, and the diffusion coefficient are needed for each region. For the buffer in the lumen, the NO reaction rate constant and the diffusion coefficient are known. In the tissue we obtain the parameters by analyzing the experimental data of Malinski et al. (16), who stimulated a microscopic portion of rabbit aortic endothelium by injecting a bolus of bradykinin at the surface of the endothelium and measured the NO concentration electrochemically. Our model of the NO production from this experiment is depicted in Fig.1 A. In the experiment of Malinski et al., one sensor was placed adjacent to an endothelial cell and another in a vascular smooth muscle cell 100 μm away. The luminal side of the endothelium is bathed in a physiological salt solution. A microscopic portion of the endothelium was then stimulated by using a microinjector to place a bolus of bradykinin in proximity to the endothelial sensor.
The size of the microscopic region producing NO, characterized bya, is unknown and must be estimated. A single endothelial cell has a diameter of 10–20 μm, so the distance between the sensors was only 5–10 times the size of the source. This implies that the size of the producing region may be an important factor in determining the NO concentration at the sensor in the smooth muscle. In fact, we anticipate the producing region to be much larger than a single cell, because the size of the bolus of bradykinin and the effect of bradykinin diffusion are unknown, even though Malinski et al. (16) stated that a single cell was stimulated.
NO is produced from the endothelium at some rate per unit area (q˙_{NO}). The NO diffuses into the vessel wall as well as into the luminal region. In the physiological salt solution in the lumen, NO reacts with O_{2} at a known rate. The reaction rate in the vessel wall is unknown.
We approximated this microscopic source by a thin, flat “pancake”shaped ellipsoid with foci at ±a. This geometry can be described using a system of oblate spheroidal coordinates, as depicted in Fig. 2. The ellipsoidal source is suited for modeling microscale NO distribution, because it can represent a single cell or a finite group of endothelial cells. In Fig. 2 the thickness of the ellipsoid (2ε) is exaggerated for illustration, typically a/ε ≈ 100 for the size of the ellipsoidal source estimated below. As ε approaches 0, the ellipsoid appears increasingly disklike, with a approaching the diameter of the disk.
Assumptions.
The following assumptions were made to simplify the analysis.
1) NO production occurs on the cytosolic membrane of endothelial cells. One of the primary enzymes that produces NO, constitutive NO synthase (cNOS), is membrane associated (2, 10). This suggests that NO production by cNOS can be approximated by a surface reaction. This production rate fixes the flux of NO at the endothelial cell surface and the corresponding boundary conditions. Thus the endothelium is treated as two producing surfaces: one corresponding to the luminal surface and one corresponding to the abluminal surface.
2) The reaction of NO depends only on the NO concentration. The NO concentration in a stimulated endothelial cell in vitro is on the order of 1 μM (16). The concentration of O_{2} in tissue in vivo is ∼27 μM (22), and the concentration in airsaturated buffer is about eight times higher. Because these concentrations are much larger than that of NO, we assume that the NO reaction rate depends only on the NO concentration (c_{NO}). In this study, two reaction rate expressions are used: NO degradation is proportional to c_{NO} (1st order) or to
3) NO diffuses as if in a binary solution. We approximate the diffusion of NO by a binary system in which NO is a trace quantity. The diffusion coefficient (D) of NO in buffer at 37°C has been estimated as 3,300 (16) and 4,500 μm^{2}/s (5). We expect the effective diffusion coefficient in the vascular smooth muscle to be less. The diffusion coefficient of NO in lipid membrane was ∼10–40% of the value in buffer (5). Because NO is dilute, the effective diffusion coefficient is assumed to be independent of concentration.
This assumption also means that we neglect the different properties of lipids and aqueous solution, so the NO concentration will be continuous throughout all regions. In general, only the chemical potential of NO is continuous. Because NO is dilute and uncharged, NO would be expected to behave ideally in solution, so partial pressure, rather than concentration, is continuous. However, there are few data about NO solution properties, so as a first approximation we will assume NO concentration to be continuous throughout.
Governing equations and boundary conditions.
Here the general form of the equations and boundary conditions are stated. They are valid for any coordinate system. In Oblate spheroidal coordinates, these equations are written for the coordinate system that describes the microscopic NOproducing source. Although not discussed in this article, the equations of this section could be used to determine the NO concentration in other systems of physiological interest.
The concentration of a diffusing, reacting substance, such as NO, is described by the species mass balance (26). For NO, this balance can be written as
The NO distribution in each of the three regions, the lumen, endothelium, and abluminal region, is governed by Eq. 1, with the regions linked through their boundaries. Six boundary conditions and three initial conditions are needed. The initial conditions are
Two boundary conditions can be fixed far from the producing surface where the NO concentration changes slowly; therefore
At the NOproducing surface, the sum of fluxes from the singular surface is equal to the surface production rate (see appendix for more detail). When the NO consumption rate on one side of the endothelial cell (e.g., if the lumen contained a hemoglobin solution) is much greater than that on the other (abluminal) side, the NO concentration will be higher on the abluminal side, causing a concentration gradient directed toward the lumen. In this case, _{ec} = 0 at the luminal surface and _{lu} is composed ofq˙_{NO,lu}(t) and a portion ofq˙_{NO,ab}(t). This is depicted in Fig. 1 B.
Oblate spheroidal coordinates.
As discussed above, we believe that the size of the source may be important and should be determined in the analysis. For this analysis the size is taken into account by approximating a group of endothelial cells as a thin ellipsoid (Fig. 1
A), which is modeled in oblate spheroidal coordinates. This description would also be valid for one cell. In this system the ellipsoidal NOproducing surface is represented by a coordinate surface. The oblate spheroidal coordinate system is shown in Fig. 2. The coordinates have the following ranges: −∞ ≤ ξ < ∞, 0 ≤ η ≤ π/2, and 0 ≤ φ < 2π (out of the plane of the page). The expressions for the gradient and the Laplacian in this coordinate system can be determined from Happel and Brenner (8) or can be found explicitly in Moon and Spencer (21). For axisymmetric NO production, the general NO mass balance (Eq. 1
) in this coordinate system is given by
The boundary conditions, Eqs. 4
and
5, for the luminal surface ξ = ξ_{0} = sinh^{−1} ε/aare
Equations 69 can be solved for concentration as a function of ξ and η. As a partial differential equation in time and two space dimensions, the solution requires significant computational effort. We did solve this problem (the 2dimensional problem), but for most of the following discussion, we use a simplified approximation of the system. To estimate the parameters of the system, repeated solution of the reactiondiffusion equation is required. The effort required to repeatedly solve the twodimensional problem, Eqs. 69, is not justified by the limited data available.
For parameter estimation and to understand the character of the solution, we want to simplify the problem in a way that allows a onedimensional solution. There are two limiting cases of oblate spheroidal coordinates: the sphere and the disk. In the limiting case of a → 0 or a finite and sinh ξ → ∞, the spheroidal coordinate surfaces become spherical. In this case, the concentration depends only on r = a sinh ξ, so the concentration dependence on η may be neglected. In the other limiting case, as sinh ξ → 0 the surface collapses to a disk. Near the edges of the disk, the concentration depends strongly on η. The concentration data were obtained at points that were on the axis of the NOproducing region, so, depending on the size of the region (a), these measurements may be less influenced by the edges. This suggests that we investigate the effect of neglecting the concentration dependence on η. Because this dependence enters through the flux term Eq. 9, let us neglect η and replace the boundary condition Eq. 9
with the approximation
Because we have neglected the η dependence of the concentration at the endothelial surface, let us also simplify the NO mass balance by neglecting the concentration dependence in the η direction everywhere
The maximum thickness of the ellipse (which corresponds to the approximate endothelial cell thickness, 2ε) is taken as 2.5 μm. The remaining boundary conditions are given by Eq. 3
For the problem solved here, the production term is described asq˙_{NO}(t) =q˙_{NO}, where the NO production per unit area is a constant. In general, the production rate would be expected to vary with time, and for the data of Malinski et al. (16), it does appear to decrease somewhat after ∼30 s. For this analysis, we consider it constant for the first 30 s and omit the subsequent data. Malinski et al. show NO concentration over a much longer period, and it appears that the simplest case, constant production, is a reasonable first approximation.
NO reaction rates.
The reaction in the smooth muscle is represented as a composite rate summing the contributions from numerous NO decomposition mechanisms. The simplest description for this composite reaction in the luminal or abluminal region is given by the empirical rate expression
NO diffusing into the lumen reacts with O_{2} in aqueous solution. This reaction has been shown to have thirdorder kinetics (27)
Numerical technique.
The onedimensional diffusion problem, Eqs. 7, 8, and 1012, was solved numerically. We transformed each partial differential equation into a system of ordinary differential equations in time by discretizing the spatial derivative (6). A similar approach has been used to model the NO concentration profile in solution resulting from NO production by a monolayer of cultured cells (14).
We used secondorder centered differencing, because it is straightforward to implement. Other discretization methods, such as orthogonal collocation, may be more efficient for computations. The derivatives in the flux boundary conditions at the cell surfaces were discretized using forward or backward secondorder accurate differencing. Because the concentration in the endothelium and close vicinity was of primary interest, variable grid spacing was used. This allowed the points near the endothelial surfaces, where the concentration gradients are steep, to be more closely spaced.
As is common for systems derived from parabolic partial differential equations, the system of ordinary differential equations was stiff. Our solution used the routines odeint, stiff, and stifbs (23); the routines stiff and stifbs were modified to take advantage of the tridiagonal structure of the Jacobian matrix. We typically used a grid with 100–200 points. A finer grid does not materially change the result. Because numerical derivatives are computed for the parameter estimation routines and the sensitivity analysis, a highaccuracy solution is required. The error criterion for stifbs was 3 × 10^{−6}. Although the boundary condition, Eq.12, applies to an infinite domain, we solved the differenced equations over a finite grid by applying the boundary condition,Eq. 12, to the outermost points. Several maximum distances were used to ensure that the result did not vary appreciably.
The twodimensional problem was solved similarly, but a uniformly spaced grid of 51 × 10 in. (ξ, η) was used. Because the Jacobian was no longer tridiagonal, time updates were obtained by solving the resulting system using the conjugant gradient method linbcg (23).
Parameter estimation.
Our objective was to determine the values for the parameter setp = (a, q˙_{NO},k
_{n,lu}, D
_{lu},D
_{ab}), in the onedimensional problem, Eqs. 7,8, and
1012. We did this by minimizing the χ^{2} merit function.
When numerical parameter estimation routines are used, it is usually desirable to scale the parameters so that they have a value on the order of 1.0. Here the scaling was done by making the parameter set nondimensional by using the following characteristic values: length (l
_{0}) of 100 μm [distance from the endothelium to the sensor in the smooth muscle in the experiment of Malinski et al. (16)], concentration (c_{0}) of 1 μM [order of NO concentration measured by Malinski et al. (16)], and diffusion coefficient (D
_{0}) of 1,500 μm^{2}/s (typical NO diffusion coefficient in lipid, assumption 3)
The χ^{2} merit function, Eq. 16, was minimized using the LevenbergMarquardt method. This method is implemented in marqcoef (23). All points were weighted equally. The experimental error is unknown, so we arbitrarily selected a standard deviation of 5% of the maximum reading; the estimated parameter values depended only weakly on this value. The parameter error estimates were made using the formal covariance matrix of the fit, which is computed by the routine marqcoef. The derivatives of the merit function with respect to the parameters were computed numerically. Numerous simulations were performed using a range of initial parameter values.
RESULTS
Parameter estimation.
The parameters were estimated for two different kinetic rate expressions for the NO reaction in the vascular smooth muscle, first and second order, and for different parameters held constant. Incases 1 and 2, D _{lu} was fixed at 4,500 μm^{2}/s, and the parameter set (a,q˙_{NO}, k _{n,ab},D _{ab}) was fit for first and secondorder reaction (Table 1). In cases 3 and 4, D _{lu} was added to the parameter set. The value forD _{lu} was estimated to account for any enhancements to luminal mass transfer that might have occurred during the experimental measurements. Examples of fluid motion that increase NO loss to the lumen include such effects as thermal gradients (which may occur in experimental work at elevated temperature) and small disturbances from the injection of the bradykinin bolus. If such effects are present, then the value estimated forD _{lu} should be considered an “effective” luminal diffusion coefficient. In case 5, we examined the parameters that would be required to model the measured concentrations if there were no NO degradation.
There were many local minima in the χ^{2} function, and many initial values in parameter space were tried. The differences of the fit between the cases in Table 1 were statistically insignificant.
The computed NO concentration profiles from the parameters in Table 1are shown in Fig. 3, along with the experimental data of Malinski et al. (16). The concentration profiles in Fig. 3 were determined for cases 1 and 2. The concentration curves computed for other cases in Table 1 were essentially superimposed over cases 1 and 2 (data not shown). The parameters were estimated using the onedimensional model, but the production rate used in the twodimensional model, 6.8 × 10^{−14}μmol ⋅ μm^{−2} ⋅ s^{−1}, was determined by scaling the onedimensional production rate by 1.25, the difference between the one and twodimensional models, as shown inappendix . It is important to know the effect of the parameters on the model prediction, so the sensitivity of the parameters on the onedimensional model are also computed. Because the approximation used to obtain the onedimensional model primarily affects the estimate of production rate, the sensitivity estimates for the onedimensional model should be valid approximations of the twodimensional model as well.
Confidence regions.
When several parameters are estimated, there are usually interactions. For example, increasing the NO production rate, decreasing the NO reaction rate, or increasing the luminal diffusion coefficient increases the NO concentration in the smooth muscle. These interactions may be quantified by computing the formal covariance matrix of the fitted parameters (Table 2). The information from the covariance matrix can be displayed graphically as joint confidence regions (23). The confidence regions of pairs of parameters for case 1 of Table 1 are shown in Fig.4. These regions depict the interdependence between pairs of estimated parameters at the 90% confidence level, with the other parameters held constant at their computed value listed in Table 1. The curves in Fig. 4 are normalized by the fitted value from Table 1.
The value with the greatest certainty is the NO production (q˙_{NO}). Its confidence regions (Fig. 4,A–C) are narrow, and much of the regions are centered on the normalized optimal value, 1. The “tilted” orientation of the confidence regions of q˙_{NO} with a andD _{ab} (Fig. 4, A and C) suggests that these pairs of parameters are correlated. However, the regions are fairly compact. There is substantial uncertainty in the value ofk _{1,ab}, as evidenced by its extended confidence regions (Fig. 4, B, D, and E). The value ofk _{1,ab} is correlated only slightly with a, and it is nearly independent of D _{ab}.
The confidence regions for case 2, the secondorder rate expression, are similar (Fig. 4, dashed curves). Figure 4, Aand F, is more elongated, reflecting larger uncertainty in the value for a. Similarly, there is much uncertainty ink _{2,ab}. The major difference between the cases is that the value of k _{2,ab} depends quite strongly ona, whereas in Fig. 4, B, D, and E,k _{1,ab} is essentially independent.
Parameter sensitivity.
A measure of the sensitivity of the NO concentration to the parameters can be determined by computing the relative change in NO concentration caused by a small change in a parameter. This measure is called the sensitivity coefficient and is defined by
Over the time course of the experiment, the sensitivity coefficients for the NO production rate are among the largest of all the sensitivity coefficients (Figs. 5 and 6). They vary little with time or reaction order. This high sensitivity explains why the NO production rate can be estimated accurately. The reaction rate coefficients have the least influence on NO concentration for the firstorder rate expression. This result explains the large variability in the estimates ofk _{1,ab} and k _{2,ab}.
The sensitivity to the diffusion coefficients is more complex. Regardless of the reaction rate expression, the concentration in the smooth muscle is very sensitive to D _{ab} at early time. This suggests that an NO probe with a fast response and measurements of early time data would be useful in determining the value for D _{ab}. At the endothelial cell,S _{D(ab)} and S _{D(lu)} are of similar magnitude for both reaction rate expressions.
DISCUSSION
The constant NO production approximation provides a reasonable fit to the data of Malinski et al. (16), as shown in Fig. 3. Although the NO production rate probably varies with time because of the delay in cNOS activation and negativefeedback inhibition from NO (24), NO production rate has the least variability, and the value is consistent over the cases examined in Table 1. Because of the physical difficulties in stimulating a microscopic region of cells, there are several other potential sources of error. It is unlikely that the NOemitting region is actually circular or that the production is precisely uniform over the region. Any of these factors could account for the difference between the experiment and the computed concentration. Despite these factors, we find that the models presented here approximate the magnitude and trend of the data reasonably well. The NO concentration predicted by the twodimensional model tracks the data trend more closely than the onedimensional approximation, as expected if the constant flux approximation is valid.
Estimation of NO production rate.
The NO production rate is required for any modeling study of NO concentration in a biologic system. For our computations we assumed this rate to be constant; however, as more is learned of endothelial NO synthase and its regulation, it may become possible to predict the NO production as a function of time,q˙_{NO}(t). The theoreticalq˙_{NO}(t) could be tested by modeling the production of NO in a manner analogous to that done here.
The estimated NO production rate was somewhat correlated with aand D _{ab} (Fig. 4, A and C and Table2), so these parameters cannot be estimated independently from this experiment. That is, q˙_{NO} will always depend on the value selected for a and D _{ab}. Despite this, q˙_{NO }varies over a fairly narrow range.
To obtain the actual production rate, 6.8 × 10^{−14}μmol ⋅ μm^{−2} ⋅ s^{−1}, we use the difference between the two and onedimensional models (see Fig. 8) to correct the estimate of Table 1.
For the firstorder rate expression,S _{q˙NO} is relatively constant, as seen in Fig. 5, A and C, so the c_{NO,ec} and c_{NO,ab} are directly proportional to NO production rate over the course of the experiment. This might be anticipated, since the boundary conditions and the mass balance equations are linear in NO production rate. Therefore, a small change in the NO production rate would proportionally change the NO concentration at each point. The secondorder rate expression is nonlinear in NO concentration, and this is reflected by the decreased dependence of the smooth muscle NO concentration on NO production rate, as indicated in Fig. 6, A andC.
Estimation of NO degradation rate constant.
There was significant uncertainty in the values of the rate coefficients that characterize the NO consumption. We were unable to differentiate between the first and secondorder rate expressions from these data. If the reaction is second order, then it may be inferred that the interaction of NO with materials in tissue is more complicated than simple binding, an interaction that would be characterized by firstorder degradation. One possible secondorder reaction has been suggested by Lancaster (11). The NO reaction in the tissue is not primarily with O_{2}, because the estimated value fork _{2,ab} is more than an order of magnitude larger than it would be if the decomposition of NO were governed by Eq.15 (k _{2,ab} ≈ 2.0 × 10^{−3}μM^{−1} ⋅ s^{−1}).
The reaction rate coefficient influences the concentration most at longer time scales, as shown in Figs. 5, A and C, and6, A and C. Therefore, the steadystate NO concentration will be useful for estimating the reaction rate coefficient. Because the NO production began to fall at 20–30 s, only the data obtained before 30 s were used in the analysis. Had longtime steadystate data been available, the estimates for the rate coefficients would probably be more reliable. Because the reaction rate depends only on local concentration, we expect the sensitivity of the one and twodimensional models to be similar, if the production rates are such that the concentration profiles are similar.
It is common to characterize the reaction of NO byt _{½}. For the firstorder reaction (Table 1,cases 1 and 3), t _{½} = 30–115 s^{−1}. For a secondorder reaction,t _{½} is determined by t _{½} = 1/(c^{0} k _{2}), where c^{0} is the initial concentration. By use of the estimated value for the secondorder rate coefficient (Table 1, cases 2 and 4),k _{2,ab} = 0.05 μM^{−1} ⋅ s^{−1}, and by setting c_{NO} = 1 μM, t _{½} ≈ 20 s. There is additional NO loss through diffusion that is not included int _{½}, so t _{½} should be used with caution.
The possibility of complicated interactions between NO and tissue has implications in modeling the reaction and diffusion of NO. If the reaction of NO with the vascular smooth muscle is first order, thenEq. 1 is linear, and the concentration in the tissue from several sources can be summed (11, 28). This will not be the case if the rate equation is nonlinear in NO concentration, such as the secondorder kinetics discussed here or MichaelisMenten kinetics. That the contribution from different sources cannot be directly added is a complication resulting from the mathematics of the nonlinear differential equations that represent the diffusion process, rather than the physics of the diffusion process itself.
Estimation of the size of the NOproducing area.
In the experiment of Malinski et al. (16), the NO concentration was measured at the endothelium and in the vascular smooth muscle 100 μm away. Because even a few endothelial cells are of this order of magnitude, a realistic model must take into account the size of the NOproducing region. Furthermore, even if the NOproducing region could be approximated by a point source or infinite plane, the validity of these approximations cannot be determined a priori, so an estimate of size is still necessary. The finite size is taken into account by solving the diffusion equation, Eq. 1, in oblate spheroidal coordinates. The values for the characteristic size a for this system are given in Table 1. They are the same order of magnitude as the distance. For this reason, a pointsource approximation of NO production would not accurately describe the concentration profile. Assuming that NO is produced from an infinite region introduces additional error, as can be seen in the limiting behavior discussed inappendix .
The values for a were fairly consistent between the cases presented in Table 1. For the firstorder rate expression, awas weakly correlated with D _{ab} andk _{1,ab}, as shown in Fig. 4 D and Table 2. For the secondorder rate expression, there was more uncertainty ina and the value was more strongly correlated withD _{ab} and k _{2,ab}, as shown in Fig.4 D and Table 2. A further effect of reaction order is that the firstorder rate expression is more sensitive to a than the secondorder rate expression, as can be seen in Figs. 5, A andC, and 6, A and C.
The sensitivity of the concentration to a for the models is shown quantitatively by the slope of the curves in Fig. 8. Fora > 150 μm the twodimensional model is somewhat more sensitive to a than the onedimensional model, as indicated by the slightly steeper slope, although the difference diminishes asa increases. Because of the similar character of these solutions, we expect that the preceding comments about the effect ofa are also valid for the twodimensional model.
The computed size of the source is much larger than a single cell but is about what might be expected from the experiment. Malinski et al. (16) injected a 10 nM bolus of bradykinin. At 0.1 M, the bolus would have a volume of 0.1 μl, so a spherical bolus would have a radius of 280 μm. Even a smaller bolus would be dispersed by diffusion and convection from injection.
Because there are only two sensors and they are on the axis of the ellipsoid, we have no direct information about the size of the region; it must be implied from the model. If offaxis data were available (i.e., there were more sensors), the onedimensional model may be found to be too simple, in which case the twodimensional model would have to be used.
Estimation of diffusion coefficient.
Previous analyses (11, 14, 28) used the diffusion coefficient determined by Malinski et al. (16) for tissue and aqueous systems. In our analysis we estimated D _{ab} while holdingD _{lu} fixed (Table 1, cases 1 and 2) or including it with the variables estimated (Table 1, cases 3and 4). When D _{lu} was estimated, we found it to be unexpectedly high. There are two reasons whyD _{lu} may be unrealistically high. First, experimental difficulties may result from minute convection currents caused by maintaining temperature control or the action of the microinjector. This point underscores the importance of mass transfer to the lumen. Second, it may be an artifact of the parameter estimation algorithm. Because the lumen can act as a mass sink whenD _{lu} is not fixed, it can offset the effect of the other parameters. Thus the parameters may be selected from a much wider range by adjusting the mass transfer to the lumen accordingly. This is why cases 3 and 4 have more variability in the parameters. In either case, the effect of D _{lu} in reducing the χ^{2} was not very large, so cases 1and 2 are considered to be the most reliable.
NO diffuses through the vascular smooth muscle with surprising speed. Even the lowest estimate, D _{ab} = 2,000 μm^{2}/s (Table 1, case 3), is relatively high compared with the diffusion rate in lipid (5). The value for the highest estimate, D _{ab} = 3,400 μm^{2}/s, is approximately that measured in an aqueous salt solution by Malinski et al. (16).
The maximum NO concentration as measured by Malinski et al. (16) (Fig.3) shows a decrease from 1.3 to 0.85 μM over 100 μm. This rapid drop in NO concentration with distance has been taken by some researchers (25) to imply that NO must be reacting with components of the cellular membrane. Although we agree that reactions do occur, this reasoning is misleading. Depending on the NO production rate and the mass transfer to the lumen, the concentration gradient may be just the result of mass transfer from a finite NO source, as seen in Table 1,case 5.
Comparison with previous mathematical models.
This analysis differs from past analysis by accounting for the size of the NOproducing region, by considering NO production to be by surface reaction, by using the time course of the concentration data, and by allowing different properties in the tissue and luminal regions.
Lancaster (11) and Wood and Garthwaite (28) used the data of Malinski et al. (16) to estimate NO production rate. Wood and Garthwaite (28) assumed an NO point source surrounded by a 0.5μmradius sphere. They estimated the strength of the point source, 2.1 × 10^{−17}mol/s, which gave an NO concentration of 1 μM at the surface of the sphere. In terms of total NO production, our estimates (Table 1,cases 1 and 2) are 1.3 × 10^{−14} mol/s. The difference is primarily the result of different assumptions of the source: point or finite.
The phenomenological production rate constant given by Lancaster (11), 10.3 μM/s, is difficult to compare directly with the production rate obtained here. However, the equation for NO diffusion from a onedimensional point source given by Lancaster (11) is identical to that for NO diffusion from an infinite plane. Therefore, in matching this equation to the steadystate data (see Fig. 2 in Ref. 11), Lancaster effectively chose an NO production rate on the basis of NO production from an infinite planar source. We have used an infinite planar source to verify our computations and found that the production rate was similar to that predicted by the oblate spheroidal geometry, although the steady state was attained much more quickly. The similarity is a consequence of the NOproducing region being larger than the distance between NO sensors.
Laurent et al. (14) assumed NO production rates of 1 × 10^{−10} to 1.6 × 10^{−8} M/s. If it is assumed that the average thickness of the endothelium is 2.5 μm, this would be equivalent to 2.5 × 10^{−19} to 4 × 10^{−17}μmol ⋅ μm^{−2} ⋅ s^{−1}. This production rate is three to five orders of magnitude less than that computed here. Because Laurent et al. (14) were interested in the NO production from activated inducible NO synthase, rather than bradykininstimulated cNOS, using a much lower production rate is reasonable.
Previous researchers considering the NO consumption rate in the tissue have assumed it to be first order, with a t _{½}ranging from 0.5 to 5 s (k _{1,ab} = 1–0.1 s^{−1}), which is much faster than the firstorder result from Table 1. This difference is probably due to two factors. First, estimates of t _{½} in experimental systems usually do not distinguish between reaction and diffusion; therefore, the rapid diffusion of NO decreases the effective t _{½}. In Table 1, the two effects are separated. Second, the reaction rate may vary with tissue location. The cardiac tissue for which Kelm and Schrader (9) estimated t _{½} of 0.1 s is rich in myoglobin, which reacts very rapidly with NO.
Wood and Garthwaite (28) observed that the reaction rate of NO is of relatively little importance near the source and in short time scales. This is consistent with Figs. 5, A and C, and 6,A and C, which show that the sensitivity of NO concentration to the rate coefficient is higher at later times. However, NO is produced by the endothelium over a fairly long time scale, so the reaction rate may be expected to be important at later times and far from the endothelium.
As pointed out by Lancaster (13), the concentration profile of a threedimensional point source falls off much too rapidly to simulate the data of Malinski et al. (16). However, for an array of point sources, the overlapping production broadens the concentration profile (11, 13, 28). Similarly, the ellipsoidal macroscopic source presented here follows the time course of the data at the endothelium and in the smooth muscle. The similarity between multiple point sources and the continuum case is because a continuous macroscopic source can be considered a distribution of point sources, when the reaction kinetics are linear. In fact, the equation for a diskshaped point source (appendix ) is derived by integrating the point source solution over the surface of the disk, effectively combining an infinity of point sources.
Implications for further modeling and experimentation.
Although not their original purpose, Malinski et al. (16) provided information that has proven quite useful in parameter estimation. As seen in Fig. 3, there is little difference in the fit or the concentration profiles between the first and secondorder reactions. Furthermore, there is considerable uncertainty in the values of the parameters, particularly in the reaction rate coefficientk _{n,ab}. A knowledge of the tissue reaction mechanism is important for understanding NO interactions in tissue and for determining the range of NO influence. Consequently, it would be useful to have a technique that clearly shows a difference in mechanisms.
An experiment similar to that of Malinski et al. (16) could be modified to distinguish between reaction mechanisms. An additional experiment in which the rate of NO diffusion into the luminal region has been decreased, possibly by decreasing the diffusion rate by adding protein to the buffer, is required. The effect of reduced diffusion can be seen in Fig. 7. When the parameters are estimated from one experiment, D _{lu} = 4,500 μm^{2}/s (Table 1, cases 1 and 2), the profiles for first and secondorder reactions are indistinguishable. By use of the same parameters, except with D _{lu}decreased to 1,000 μm^{2}/s, the NO concentration profiles for the first and secondorder cases can be distinguished. Thus experiments using two different diffusion coefficients can be used to distinguish between reaction expressions. In either case, it is important to prevent convective luminal mass transfer. The presence of even a slight amount of convection can change the interpretation of experimental data significantly. As seen in Table 1, case 5,enhanced mass transfer to the lumen can make pure diffusion indistinguishable from diffusionreaction.
Although not discussed in detail, the general equations, Eqs.15, could be used to determine the NO concentration in other systems of physiological interest: cylindrical coordinates to describe a blood vessel or rectangular coordinates to describe cultured endothelial cells or an experiment similar to that of Malinski et al. (16) in which a large region of the endothelium was stimulated. Furthermore, modeling the NO production as a surface reaction is particularly suited to describing macroscopic NO production when the thickness of the source is much less than the distance over which diffusion occurs, such as the endothelium of a blood vessel.
Acknowledgments
This work was supported by a Whitaker Foundation Biomedical Engineering Grant.
Footnotes

Address for reprint requests: J. C. Liao, Dept. of Chemical Engineering, University of California, 405 Hilgard Ave., Los Angeles, CA 900951592.

Present address of M. W. Vaughn: Dept. of Chemical Engineering, University of California, Los Angeles, CA 90095.
 Copyright © 1998 the American Physiological Society
Appendix
Boundary conditions at the endothelial surface.
Production and mass transfer of a substance from the singular surfaces bounding the endothelium are described by the surface mass balance (sometimes called the jump mass balance) (26). The NO consumption in the luminal region is different from that in the abluminal region (and possibly the endothelial cell itself ). Therefore, the NO gradient and the flux of NO in each of these regions will also differ. The jump mass balance relates the fluxes to the production rate. For NO, this balance can be written
An example of nonzero interface velocity u occurs during changes in tone of cylindrical blood microvessels as part of the integrated control of the microcirculation. As the vessel dilates and contracts, the endotheliumlumen interface moves; therefore, u≠ 0.
The mass flux vector ≡ c_{NO}
v
_{NO} (mass transported per unit time per unit area) will, typically, vary with time. For a trace quantity like NO, it can be expressed in the form of Fick’s law
Appendix
Comparison of oblate spheroidal models and other limiting cases.
It would be helpful to have an analytic solution of the diffusion equation in oblate spheroidal coordinates to assist in verifying the numerical computations. Because predicting the concentration as a function of time is central to our parameter estimation, steadystate approximations, which are much simpler, are less interesting. However, even for the simplest case, i.e., no reaction, the solution is quite complicated, so we look for limiting cases to verify the computations. If the size of the source is neglected, there are two alternative treatments: 1) the extent of the source is very large, as an infinite plane, or 2) it is very small, as a point source. These treatments have the advantage that the NO diffusion is one dimensional for constant surface flux, perpendicular to the plane or in the radial direction, as from a point source.
Because oblate spheroidal coordinates are used to approximate a diskshaped NO source, we can use this geometry to test the validity of our concentration profile. This is illustrated in Fig.8, where we plot the diffusiononly (V˙_{NO} = 0) concentration at the producing surface and 100 μm from the source for a range of source sizes a by using the homogeneous disk source solution Eq.24 and the one and twodimensional numerical models in oblate spheroidal coordinates. For each source, a production rate of 5.4 × 10^{−14}μmol ⋅ μm^{−2} ⋅ s^{−1}is used, and the concentration was computed at 12 s after production began, about when experimental data (16) show the concentration at the source reached a maximum. We see that the concentration profile for the twodimensional model for a thin ellipse is very near the theoretical disk value. As expected, the simplified onedimensional model predicts a higher concentration for smaller values of a. For the range of interest, a > 150 μm, the error between the simplified onedimensional model and the twodimensional model is a maximum of 25%. Because the difference in predicted concentration is most likely due to the underestimate of production rate (see discussion inOblate spheroidal coordinates), this deviation is used to adjust the production rate estimated by the onedimensional model. The deviation could also be viewed as an underestimate in the estimated size a of the region, as can be seen by shifting the onedimensional model curves to the right. This interdependency betweena and q˙_{NO} can be expected from covariance analysis (Fig. 4 A).
Even though the production rate is adjusted using diffusiononly (V˙_{NO} = 0) computations, the adjustment should be valid when V˙_{NO} ≠ 0. This is because the NO concentration is determined primarily by diffusion. To see this, note that the dimensionless form of the diffusion equation for the abluminal region with firstorder NO degradation is
The analytic solution for a diskshaped source on a plane can be obtained in the following manner. An ellipsoid with foci at ±a approaches a disk in the limit of vanishing thickness. The differential mass balance for a disk of radius a centered at the origin of xy plane is
Equation 24 is valid for production from an infinitely thin disk when there is no reaction. In the experimental system, NO would be lost to reaction. Additionally, there may be a higher diffusive flux into the lumen where D is higher. For a given D,Eq. 24 can be used to obtain a lowerlimit estimate ofa and q˙_{NO}. For this estimate, we do a simplified fit of the data of Malinski et al. (16) by using only the maximum concentration points at the endothelium (1.3 μM, 12 s) and in the smooth muscle 100 μm away (0.84 μM, 20.3 s). The results are as follows: a = 270 μm and q˙_{NO} = 4.8 × 10^{−14}μmol ⋅ μm^{−2} ⋅ s^{−1}.
Had a large area of the endothelium been stimulated, the infinite plane source solution (3) could have been used to estimateq˙_{NO}. With production estimated for the disk source of a = 270 μm, the concentration at 20.3 s and 100 μm from an infinite planar source is 63% higher than the disk source estimate. For a ≫ 100 μm, the solutions rapidly converge, and for a = 800 μm, they differ by only 2%.