## Abstract

The fractional concentration of nitric oxide in exhaled breath (fe_{NO}) is a biomarker of airway inflammation with applications in clinical asthma management and environmental epidemiology. fe_{NO} concentration depends on the expiratory flow rate. Standard fe_{NO} is assessed at 50 mL/sec, but “extended NO analysis” uses fe_{NO} measured at multiple different flow rates to estimate parameters quantifying proximal and distal sources of NO in the lower respiratory tract. Most approaches to modeling multiple flow fe_{NO} assume the concentration of NO throughout the airway has achieved a “steady‐state.” In practice, this assumption demands that subjects maintain sustained flow rate exhalations, during which both fe_{NO} and expiratory flow rate must remain constant, and the fe_{NO} maneuver is summarized by the average fe_{NO} concentration and average flow during a small interval. In this work, we drop the steady‐state assumption in the classic two‐compartment model. Instead, we have developed a new parameter estimation approach based on measuring and adjusting for a continuously varying flow rate over the entire fe_{NO} maneuver. We have developed a Bayesian inference framework for the parameters of the partial differential equation underlying this model. Based on multiple flow fe_{NO} data from the Southern California Children's Health Study, we use observed and simulated NO concentrations to demonstrate that our approach has reasonable computation time and is consistent with existing steady‐state approaches, while our inferences consistently offer greater precision than current methods.

## Introduction

The fractional concentration of nitric oxide in exhaled breath (fe_{NO}) is a biomarker of airway inflammation with clinical (e.g., asthma) and research (e.g., environmental epidemiology) applications. Nitric oxide (NO) is produced in endothelial cells in airway tissue. The discovery (Gustafsson et al. 1991) that humans exhale measurable quantities of nitric oxide (NO) was soon followed by the discovery that, for a given subject, the concentration of NO exhaled depends strongly on the exhalation rate (Högman et al. 1997; Silkoff et al. 1997). To control for this effect, guidelines for assessment (ATS 1999; ATS/ERS 2005) and interpretation (Dweik et al. 2011) of the fractional concentration of exhaled NO (fe_{NO}) have been developed around a standardized exhalation rate of 50 mL/sec (fe_{NO,50}). A significant drawback to this approach is that this flow rate provides information on NO arising primarily from proximal airway wall sources (George 2008).

By measuring fe_{NO} at multiple flow rates, it becomes possible to partition the sources of NO into distinct anatomical subregions (Tsoukias and George 1998; Pietropaoli et al. 1999; Högman et al. 2000; Silkoff and Sylvester 2000). In the widely used two‐compartment model (George et al. 2004), the respiratory system is divided into alveolar and airway compartments, and the contribution to fe_{NO} from each estimated separately. Most approaches to modeling multiple flow fe_{NO} assume the concentration of NO throughout the airway has achieved a “steady‐state.” In practice, this assumption demands that subjects maintain sustained flow rate exhalations, during which both fe_{NO} and expiratory flow rate must remain constant. The fe_{NO} maneuver is then summarized by the average fe_{NO} concentration and average flow over a small interval during which the steady‐state assumption appears reasonable. Subjects may only be tested using multiple flow sampling protocols if they can control their exhalation rate with sufficient precision to achieve constant expiratory flows near the protocol targets. Adequate expiratory control may be difficult even for healthy adults (George et al. 2004) or impossible for other groups, such as young children (Linn et al. 2009). Even with excellent expiratory control, fe_{NO} maneuvers at high flows can often display patterns inconsistent with the steady‐state assumption (Puckett et al. 2010).

To overcome limitations of the steady‐state assumption, the goal of this paper is to drop this assumption within the context of a two‐compartment model with a cylindrical airway. Instead, we have developed a new estimation approach based on measuring, and adjusting for, a continuously varying flow rate over the entire fe_{NO} maneuver. Our methodology is similar to the “backward integration” approach introduced in Tsoukias et al. (2001) to analyze samples based on a single breath, although we make fewer simplifying assumptions and explicitly account for the effect of axial diffusion. Because we make no a priori assumptions about the flow rate, there is no inherent need for subjects to control their breathing. Our approach uses measured flow data to continuously adjust the model; therefore, we can analyze data gathered at continuously varying flow rates. This offers the potential for our approach to enable extended fe_{NO} testing with subjects unable to perform existing multiple flow protocols.

In this study, we present the dynamic two compartment model and the Bayesian inference framework we developed to estimate the parameters of the partial differential equation underlying this model. Our approach was originally motivated by mutltiple flow fe_{NO} data from the Southern California Children's Health study. Using existing data from this study, we create a simulated sample of randomly selected subjects to compare our dynamic approach with existing steady‐state approaches. We show that our inference method yields parameter estimates with both superior accuracy and precision compared with existing steady‐state methods applied to the same data. We also show examples of estimates generated using real data, and we end by discussing potential simplifications to the sampling protocol our method may enable.

## Glossary

*r*and*l*are the airway radius and length (respectively), in cm.*z*_{0},*z*_{mouth}, and*z*_{alv}are the locations of the sensor, mouth, and alveolar boundary, in cm from*z*_{0}.*c*(*z*,*t*) is a solution of equation 1, indicating the NO concentration at position*z*and time*t*, in ppb.*t*_{0},*t*_{1}, … ,*t*_{n}are the*n*+1 discrete measurement times, where exhalation begins at*t*_{0}and ends*t*_{n}.*c*_{i}:=*c*(*z*_{0},*t*_{i}) is the model solution at the sensor*z*_{0}, is a numerical approximation of*c*_{i}, and is the measured concentration, all in ppb and at time*t*_{i}.*v*(*t*) is the linear flow rate, in cm/s.*d*is the diffusivity of NO in air, in cm^{2}/s.*p*is the permeability of airway wall tissue to NO, in cm/s.*c*_{w}is the concentration of NO in the airway wall tissue in ppb, which is assumed to be constant.

The parameters described in the steady‐state modeling review (George et al. 2004) can be related to the dynamic model as follows:

fe

_{NO}=*c*(*z*_{0},*t*). The concentration of NO exhaled (also denoted ce_{NO}) corresponds to the concentration at the sensor (*z*_{0}), and it is the only parameter that can be measured directly. In the dynamic model this varies with time, so it is also a function of*t*.ca

_{NO}=*c*(*z*_{alv},*t*). The gas phase alveolar concentration corresponds to the concentration at the alveolar boundary (*z*_{alv}). In principle this may vary with time, but often it is assumed to be constant on short (seconds‐minutes) time scales.Daw

_{NO}= 2*πrlp*. The total airway diffusing capacity is equivalent to the product of the airway surface area and the coefficient*p*. It also equals the product of the airway volume and the source term coefficient 2*p*/*r*in equation 1.J

^{′}aw_{NO}= 2*πrlpc*_{w}. The maximum total flux of NO in the airway is equivalent to the product of the airway surface area and the coefficients*p*,*c*_{w}. It also equals the product of the airway volume and the source term constant 2*pc*_{w}/*r*in equation 1.Caw

_{NO}=*c*_{w}.

## Methods

Our proposed airway model is a variant of the two‐compartment approach (Tsoukias and George 1998), the primary distinction being we make no assumptions regarding the flow rate. In its simplest form, the two‐compartment airway is assumed to be a cylinder with fixed dimensions (Fig. 1). Unlike the airway compartment, the dimensions of the alveolar compartment may vary. However, at any moment in time the NO concentration is assumed to be constant throughout the alveolar compartment, that is, it is “perfectly mixed.” The original description of the two‐compartment model incorporated a time varying alveolar concentration (Tsoukias and George 1998); in practice it is often assumed to be constant on short (seconds‐minutes) time scales.

The airway cylinder is lined with epithelial tissue containing NO producing cells. The tissue is assumed to be NO permeable, and some of this NO will diffuse from the tissue into the lumen. Exterior to this tissue is bronchial blood, which is assumed to serve as an infinite sink for NO (Tsoukias and George 1998). The tissue NO concentration is assumed to be constant; therefore, depending on the relative NO concentration between incoming air and the airway wall, the airway tissue serves as either an infinite source or sink for airway NO. Although the biological airway ends at the mouth, the cylinder is extended by the instrument dead space volume. In this region, the “airway wall” is assumed to be impermeable to NO; otherwise, it is modeled in the same manner as the rest of the airway. In this regard, the model is equivalent to a cylindrical model with a piecewise constant airway wall permeability.

### The governing equation

The dynamics of NO in the airway are assumed to be governed by the partial differential equation (PDE): (1)

The three quantities on the right hand side are known as the advection, diffusion, and reaction terms (respectively). In this context the last quantity is also known as a source term, as it represents another source of NO, namely, the airway wall. The contribution from the airway wall is assumed to be proportional to the difference in concentration between the wall and the lumen. More details regarding the physical assumptions underlying the result (eq. 1) can be found in Appendix A.

Our approach to estimation and inference is predicated upon repeated simulation of the underlying physical model. In this framework, “simulating” the model (eq. 1) largely consists of calculating a series of numerical solutions , where exhalation begins at *t*_{0} and ends at *t*_{n}. The method of lines (MOL) technique is applied to the PDE (eq. 1), wherein the spatial (*z*) variable is discretized using finite differences: upwind for the advective term, and centered for the diffusive (see Appendix B for details).

Replacing the spatial derivatives with finite difference approximations yields a large system of ordinary differential equations (ODEs). The time variable remains continuous, and the resultant semi‐discrete problem can be solved numerically when combined with appropriate boundary and initial conditions (discussed in Appendix B). An advantage of this approach is that “off‐the‐shelf” routines designed for arbitrary systems of ODEs can be used to perform the integration (Hundsdorfer and Verwer 2003; LeVeque 2007). Calculating the solution of equation 1 also requires specifying the velocity function *v*(*t*). In a sense, *v*(*t*) “drives” the solution, because it is the only term on the right hand side of equation 1 that varies with time.

### Southern California Children's Health Study data

We demonstrate data processing and model estimation and inference using data from the most recent cohort of the Southern California Children's Health Study (CHS), originally recruited in 2002–2003. Multiple flow fe_{NO} data were collected in March–June 2010 from 1640 children, ages 12–15, in eight of the CHS communities. The study protocol requested that participants perform nine constant flow fe_{NO} maneuvers at four target flow rates: three at 50 mL/sec, and two each at: 30 mL/sec, 100 mL/sec, and 300 mL/sec (Linn et al. 2013).

Samples were collected online using Ecomedics Analyzers (CLD88‐SP with DENOX module), which employ chemiluminescence to measure NO concentrations and ultrasound to measure flow rates. The flow rate measurements are available almost instantaneously, while there is a small delay in the NO signal. This delay is due to a combination of the time required for gas transport between the flow head and the sensor, and a sliding average filter applied to the output signal. The total delay is ∼1 sec, and during testing the analyzer automatically adjusts for this to provide synchronized flow and NO time series output. The NO and flow rate measurements were exported into raw text files at 100 samples/sec. Many of these values are redundant, however, and the effective sampling rates are 50 Hz for flow and 15 Hz for NO. While these data were synchronized at the time of collection, that is not required, and post hoc corrections could also be applied prior to analysis.

### Flow rate data preprocessing

We use a Dormand and Prince (1980) based routine to perform the time integration of equation 1. Like most automatic integration routines, this method adaptively varies the time step size based on running error estimate calculations. Sharp changes in the NO concentration require shorter time steps be taken, dramatically increasing computation time. Because the concentration is flow dependent, sharp changes in the flow rate can precipitate sharp changes in the concentration.

Since the solution is calculated at adaptively chosen times, in general it may be necessary to approximate *v*(*t*) at arbitrary values of *t*. Therefore, measurements made at a fixed sampling frequency must be used to estimate the continuous velocity input required by the integration routine. An example of typical flow data for a CHS participant, at the target flow of 50 mL/sec, is shown in Figure 2. Naively interpolating the raw time series can lead to spurious high frequency oscillations in flow rate estimates, mimicking the computational challenges introduced by sharp rate changes.

To dampen these oscillations, the flow rate data is run through a low‐pass frequency filter (Smith 2007). Because the data is analyzed after all of it has been collected, two pass (forward‐backward) filtering is employed via a fourth order Butterworth filter, with a low‐pass frequency threshold of 2.5 Hz. As Figure 2 shows, filtering the signal in this manner retains the gross features, such as the spikes at the beginning, while eliminating the rapid oscillations later on. The estimated flow rate function is defined by interpolating the filtered signal.

### Model simulation

As illustrated in Figure 1, the position of the sensor is defined to be the origin, *z*_{0} = 0. Therefore, the solution at this point over time corresponds to the model prediction of fe_{NO} measured throughout the maneuver; informally, =. In addition to the estimated flow rate function , the approximate solutions also depend on the airway parameter values. In these examples, ca_{NO} = 2, J^{′}aw_{NO}=800, and Daw_{NO} = 5 (values identical to those used in a previous simulation study [Citation Eckel et al. 2014]).

Combining with the parameters Ca_{NO}, J^{′}aw_{NO}, and Daw_{NO} enables the use of numerical integration to calculate the sequence of approximations . The solid line in Figure 3 is the same filtered flow data as shown in Figure 2. The dashed line illustrates the predicted concentration at the sensor throughout the exhalation (synchronized with the flow, so the time scale is shared). The solution is calculated at 100s‐1000s of time steps *t*_{i}, resulting in a squence of approximations . Because the integration routine adjusts the step size according to the state of the system, the precise number of steps will vary depending on the input values.

### Multiple flow study simulation

The multiple flow fe_{NO} sampling protocol employed by the CHS was designed to facilitate estimation via existing steady‐state models, and it involved participants performing nine fe_{NO} maneuvers at four target flow rates. Thus, we can validate our model by applying our own method, along with existing estimation techniques, to the simulated data, and then comparing the resultant estimates to the (known) values used during simulation.

Because an exhalation flow profile is such a complicated process, we do not attempt to recreate it via simulation; rather, we use real flow data from the CHS to generate the filtered flow functions used during simulations. Specifically, we used flow rate data from a stratified random sample of 100 CHS subjects. Sampling strata were determined by fe_{NO,50} level, with 25 subjects selected from each of the categories (in ppb): <10, 10−25, 25−50, and ≥50. The sample was also restricted to subjects that successfully completed the nine maneuvers specified in the protocol.

For every subject, each of their nine maneuvers was used to simulate theoretical NO concentrations *c*_{i} as a function of time by combining the observed flow data with the model equation 1 (still using the parameter values ca_{NO}=2, J^{′}aw_{NO}=800, and Daw_{NO}=5). The deterministic PDE model leading to the dashed line in Figure 3 is capable of accurately describing many of the qualitative features of exhaled nitric oxide (Shin and George 2002). However, there will inevitably be some deviation between the model prediction *c*_{i} and the observed value . To account for this residual variation, we assume is normally distributed with mean ln (*c*_{i}) and variance *σ*^{2} = 0.1^{2}, that is has a log‐normal distribution. This is the same subject shown in Figure 3 (and Fig. 2), so the dashed line in the top left panel of Figure 4 is identical to the corresponding line in Figure 3. The dotted line is the result of adding independent normal errors to the (log‐transformed) deterministic solution, and the other 8 panels in Figure 4 are the result of repeating this process with the remaining flow profiles for this subject.

### Steady‐state vs dynamic estimation

The most common approaches to estimating ca_{NO}, J′aw_{NO}, and Daw_{NO}, are based on the steady‐state model (George et al. 2004)
(2)relating the parameters of interest to the exhaled concentration ( is the volumetric flow rate, which is equivalent to the product *v*(*t*)*πr*^{2} of the linear flow rate and the cross‐sectional area of the airway cylinder). The model (eq. 2) is a special case of the model (eq. 1), which arises by imposing additional simplifying assumptions (detailed Appendix A).

The direct approach to estimation, introduced in Silkoff et al. (2000), employs a nonlinear regression model by incorporating an additive error term, that is fe_{NO} = ( … )+*ε* (nonLin). A more recent variant, introduced in Eckel et al. (2014), first log transforms the data before performing the nonlinear regression, that is log (fe_{NO}) = log (…)+*ε* (nonLinLog). This method can also constrain parameter estimates to be non‐negative, ensuring consistency with the physical assumptions.

Earlier approaches replaced the exponential function with a polynomial expansion, such as the linear approximations employed by Pietropaoli et al. (1999) in the model (linP), and the model (linT) described in Tsoukias et al. (1998). A linear function is a poor approximation of the exponential for values near 0, hence higher order polynomials have also been used. The algorithm described in Högman and Meriläinen (2007) (HMA) employs a third order approximation, in addition to ensuring the estimate of Daw_{NO} is positive. While the majority of two‐compartment models employ a cylindrical airway, the model described in Condorelli et al. (2007) (Condorelli) employs a trumpet airway shape, along with a positive diffusion coefficient *d* > 0.

The dynamic approach differs from existing methods by dropping the steady‐state assumption. It is based on estimating the trajectory of fe_{NO} through all phases of exhalation, and there is no requirement for, nor even notion of, fe_{NO} “plateaus.” This is achieved by continuously measuring the flow rate and adjusting for its impact on fe_{NO}. Because the estimates generated by the dynamic model are based on a broader spectrum of input values, it has the potential to calculate estimates with greater precision than existing methods. Parameter estimation in this framework is done via Monte Carlo methods, and the point estimates reported here are maximum a posteriori probability (MAP) values generated during Markov chain Monte Carlo (MCMC) sampling, see Appendix C for details (for readers unfamiliar with Bayesian methodology, in this context the reported MAP values are equivalent to maximum likelihood estimates (MLE)).

To asses the uncertainty in point estimates, some methods also produce interval estimates for the parameters of interest. The linP and linT point estimates are based on ordinary least squares (OLS), and confidence intervals can be derived from the assumptions underlying OLS. One of these assumptions is that the residuals are normally distributed. However, when applied to fe_{NO} the true residual distribution is unknown, potentially negating the validity of the intervals. The nonLin and nonLinLog methods employ non‐linear least squares to calculate point estimates, and the associated intervals are based on the asymptotic normality of the maximum likelihood estimator. An advantage of this approach is that it does not require explicit distributional assumptions; however, the results only become exact as *n*→∞. In the context of fe_{NO} sample sizes tend to be small, for example *n* = 9 in the case of CHS subjects, potentially calling into question a large sample approximation. The dynamic model estimates are based on a Bayesian approach, and the associated parameter range estimates are known as credible intervals. These intervals do not require assuming (asymptotic) normality. However, in general they cannot be explicitly calculated, and instead must be estimated with a numerical procedure such as MCMC (as we do here).

## Results

The simulated maneuver‐level data was processed according to the original study protocol, so for each of the 100 simulated subjects, nine estimates of the steady‐state plateau average fe_{NO} concentration and flow rate were generated. Estimates of the parameters ca_{NO}, Daw_{NO}, and J^{′}aw_{NO} were then calculated using two steady‐state approaches (HMA and nonLinLog), along with our novel dynamic approach. The steady‐state methods use the nine estimated plateaus as inputs, while the dynamic approach uses the entire trajectory of each of the nine maneuvers for estimation.

### Simulation study parameter estimates

The box plots in Figure 5 illustrate the distribution of point estimates generated by all methods for each parameter. The broken horizontal lines indicate the values used to generate the simulated data, hence a “correct” inference would recover these values. For all three parameters the dynamic estimates are centered around the true values, while simultaneously possessing the narrowest sample distribution of any method. The median estimates for J^{′}aw_{NO} and Daw_{NO} generated by the HMA and nonLinLog methods also appear very close to the true values, but there is significantly greater variability in the distribution of estimates. The HMA and nonLinLog methods also produced estimates of ca_{NO} that tend to have a positive bias.

To quantify the results illustrated in Figure 5, the mean and standard error for each sample are shown in Table 1. If we assume the 100 simulated subjects are independent, we can use these estimates to conduct a *t*‐test of equality between the mean of the sampling distribution for , and the known parameter values used during simulation. Because this equality defines an unbiased estimator, rejecting this hypothesis would be evidence of bias in the estimation procedure. The *P*‐*value* columns in Table 1 report the results of these tests using a two‐sided (≠) alternative. For all three parameters, these results support the hypothesis that the dynamic estimates are unbiased. Moreover, it is the only method to provide an unbiased estimate of ca_{NO}, and in every case it yields the smallest standard error.

The steady‐state assumption requires the time derivative be zero, (∂/∂*c*)*c*(*z*,*t*)=0, reducing the PDE (eq. 1) to a second order ODE. By also assuming the diffusivity of NO in air is zero (*d*=0), the model simplifies further into the first order ODE underlying the equality in equation 2. In a cylindrical airway model, we found the impact of neglecting the diffusion term to be modest, but (statistically) significant.

To demonstrate this, the calculations leading to the dynamic estimates reported in Table 1 were repeated with *d*=0. This produced sample mean estimates for ca_{NO}, J^{′}aw_{NO}, and Daw_{NO} of (2.00, 804, 5.43) (respectively), with corresponding standard errors of (4.71e‐03, 1.30, 0.091). Compared with the true values of (2,5,800), the estimate of ca_{NO} is unaffected (*P*=5.10e‐01), while there is evidence of small but significant positive biases in the mean estimates for both J^{′}aw_{NO} (*P*=3.06e‐03) and Daw_{NO} (*P*=8.00e‐06).

### Application to real CHS fe_{NO} data

The plots in Figure 5 and results in Table 1 show our estimation routine reliably recovers known parameter values employed to generate simulated data. To demonstrate our approach applied to real data, we also estimated ca_{NO}, J^{′}aw_{NO}, and Daw_{NO} based on the observed fe_{NO} time series for the CHS subject whose flow samples underlie the simulated profiles in Figure 4.

The dotted lines in Figure 6 illustrate the measured fe_{NO} profiles for this subject. While the observed NO profiles have shapes similar to the profiles generated for the simulation study, the fe_{NO} values are significantly (2‐4x) higher in the real data. The dashed lines in Figure 6 illustrate the predicted model solutions based on MAP estimates of the parameter values.1 These parameter estimates, which are shown in the first row of Table 2, were calculated using the MCMC sampler described in Appendix C. Plateau NO concentrations and flow rates were also estimated using the real data, and the same steady‐state methods as before were employed to calculate the other point estimates shown in Table 2.

While there is no definitive way to asses the accuracy of these estimates, some qualitative features are worth noting. The estimates in Table 2 are significantly (2‐4x) larger than the values used in the simulation study, which is consistent with the fact that this subject has unusually high fe_{NO}. As in the simulation study, the HMA and nonLinLog methods produce very similar estimates. Perhaps more interestingly, the rank ordering of the estimates in Table 2 nearly matches the rank ordering of estimates in Table 1.

## Discussion

We have developed a parameter estimation framework for fe_{NO} which treats the observed flow rate as an exogenous input. Multiple flow testing already necessitates measuring the flow rate to ensure protocol compliance; therefore, the flow data we require is, in principle, readily available using current sampling technology. Employing flow data gathered during the CHS, we validated our model by applying current steady‐state methods to 100 subjects in a simulated multiple flow study. We have also developed a novel approach to estimation, and adopting a flow‐adjusted model allows us to use the entire fe_{NO} trajectory for inference. Amongst the simulated subjects and an example CHS participant, this approach significantly improves the precision of flow‐independent parameter estimates.

### Relationship to existing methods

A restriction inherent in a cylindrical airway model is that the cross‐sectional area is fixed. In reality, this area increases with airway generation (Weibel 1963). An alternative airway geometry that incorporates this feature is the “trumpet” model (Shin and George 2002; Condorelli et al. 2007). In the trumpet model, the cross‐sectional area of the airway increases with airway depth; therefore, for the volumetric flow rate to be conserved throughout the airway, the linear flow rate must decrease with airway depth. A cylindrical model might thus be expected to underestimate the significance of the diffusion coefficient *d*, which is corroborated by the finding that axial diffusion plays a larger role in a trumpet versus cylindrical model (Shin and George 2002). This also suggests that the dynamic parameter estimates with *d*=0 may underestimate the true effect of neglecting axial diffusion.

Assumptions regarding *d* are typically made because of computational concerns. The same is true of the assumption that (∂/∂*c*)*c*(*z*,*t*)=0; however, there are significant practical implications to this assumption as well. Because fe_{NO} is flow dependent, subjects must perform sustained exhalations at constant flow rates (i.e., in eq. 2 must be held constant). Official guidelines suggest these plateaus last for at least 3 sec, and that total exhalations last for at least 4 sec in children < 12, and 6 sec for children > 12 and adults (ATS/ERS 2005). Most multiple flow protocols recommend 2–9+ maneuvers be performed (George et al. 2004), and in general performing more maneuvers leads to better parameter estimates. The need to perform repeated, extended exhalations at well controlled flow rates has had a significant negative impact on the potential for routine multiple flow fe_{NO} testing in a clinical setting.

For comparison purposes, the dynamic model estimates are based on the same nine profiles as the steady‐state methods; however, it is not limited to this type of data. By adjusting the model for the measured flow, there is no inherent need for subjects to execute a controlled breathing pattern. Dropping the steady‐state assumption ((∂/∂*c*)*c*(*z*,*t*)=0) potentially has significant practical implications, which is discussed further in the last section.

While the majority of techniques for parameter estimation employ constant flow maneuvers, variable flow approaches have been explored as well, and in Tsoukias et al. (2001) it was shown that all three parameters could be estimated from a single maneuver with a varying flow rate. Similar to our approach, this method uses measured flow data to adjust for the effects of a variable flow rate by performing a backward integration of the flow signal to estimate the airway residence time of an arbitrary bolus sampled at time *t*. Airway NO dynamics are also assumed to be governed by a similar PDE, although the governing equation neglects to account for axial diffusion. Because of this, the model does not precisely predict phases I and II of the profile (e.g., fig. 5 in Tsoukias et al. (2001)), which is a period that should theoretically have high sensitivity to Daw_{NO}. By accounting for axial diffusion our approach should be able to better model phases I and II, and by extension could potentially produce more accurate Daw_{NO} estimates. No off‐the‐shelf software is available to implement their approach, and future work should compare/contrast the two methods under a variable flow (or tidal) breathing protocol.

### Limitations

The most significant limitation to our approach is that it is computationally demanding. The estimation routine depends on calculating a numerical solution to equation 1 thousands of times. The single biggest determinant of computation time is the spatial (*z*) resolution required to resolve the dynamics of NO in the airway. Previous studies have found that Δ*z* = 0.2 cm is sufficient (Shin and George 2002); however, our simulations indicate this is inadequate. Specifically, the simulated profiles illustrated in Figure 4 are based on using 1000 grid points (corresponding to a spatial resolution Δ*z*=0.025 cm). We then ran our estimation routine, beginning at 100 grid points, and increasing the number of points until further increases no longer impacted the results. Our experiments show at least 300 points are required, corresponding to a step size of cm. Using a modern desktop CPU, our estimation procedure runs in 5–10 min per subject.

Our approach to parameter estimation requires detailed time series on both flow and fe_{NO}. During the CHS these data were gathered using commercially available analyzers which are capable of recording at 10s (fe_{NO}) to 100s (flow) of Hz. However, some newer devices output just a single plateau value (Cristescu et al. 2013; Maniscalco et al. 2016), which is inadequate for fitting our models.

Although the cylindrical two‐compartment model can explain many features of fe_{NO}, the airway is not a cylinder, and more sophisticated airway shapes can explain phenomena the standard model cannot. An example would be the steady downward slope observed in the last two maneuvers in Figure 6. One potential explanation for this phenomena is offered by the multi‐compartment, trumpet shaped model introduced in Suresh et al. (2008). In Shelley et al. (2010), this model was shown to be capable of explaining fe_{NO} profiles which decline continuously throughout maneuvers that satisfy official guidelines regarding flow rate stability. Combining our estimation routine with a more realistic airway model could be a natural way to improve the fit of some maneuvers over the standard two‐compartment approach.

### Future directions

While more sophisticated airway shapes can potentially explain trajectories like those in Figure 6, they have typically only been applied to constant flow rate data, and often only in a laboratory setting. The primary appeal of our methodology is that there are no inherent requirements regarding breathing behavior. This potentially opens the door to estimation based on tidal breathing patterns, essentially eliminating the need for subject cooperation during testing. While some preliminary work attempting to use tidal data as a proxy for fe_{NO,50} exists (van Mastrigt et al. 2014), accurately and reliably estimating ca_{NO}, J^{′}aw_{NO}, and Daw_{NO} based on tidal fe_{NO} would be a significant advancement.

From a theoretical perspective, tidal breathing patterns may be preferable to a multiple flow protocol. While current multiple flow protocols typically gather data at 2‐5 different rates, during tidal breathing the observed flow rate continuously ranges over a spectrum of values. Therefore, the dynamic model can potentially be used to estimate the parameters of clinical interest with equal, if not greater precision, by allowing the flow rate to vary continuously. Moreover, this would simultaneously simplify the process, expanding the pool of eligible patients.

## Conflict of Interest

None declared.

## supplemental Data

**Appendix S1.** Supplemental information. [phy213276-sup-0001-AppendixS1.zip]

## Footnotes

↵1 These plots illustrate the model solution over time at one point in space (the NO sensor); however, simulating the model involves calculating the solution over the length of the airway (from the alveolar boundary to the sensor). To illustrate the time evolution of the solution over the entire domain, for each maneuver in Figure 6 an animation illustrating the dynamics of NO throughout the airway has been made available in Appendix S1.

**Funding Information**This study was funded by the National Institute of Environmental Health Sciences, 1K22ES022987, 4R01ES023262, 4T32ES013678, 5P30ES07048.

- Manuscript Received: March 31, 2017.
- Manuscript Accepted: April 5, 2017.

- © 2017 The Authors.
*Physiological Reports*published by Wiley Periodicals, Inc. on behalf of The Physiological Society and the American Physiological Society

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

## References

## Appendix A Model description

The airway is assumed to be a cylinder with constant radius *r* and constant length *l* (one or both may be implicitly defined by specifying the airway volume *πr*^{2}*l* and/or surface area 2*πrl*). The model (eq. 1) results from imposing conservation of (NO) mass throughout this domain. That is, for any time interval (*t*,*t* + Δ*t*), we assume any net change in mass inside an arbitrary region (*z*,*z* + Δ*z*) equals the net mass passing through the regions's boundary over the same time interval.

#### Advection of NO axially through the airway

Over the time interval (*t*,*t* + Δ*t*), a cross‐section of air moving at *v*(*t*) cm/s will travel a distance Δ*z*≈*v*(*t*)Δ*t* cm. Therefore, over a time interval of length Δ*t* the mass that flows into the region through the lower boundary is approximately *v*(*t*)Δ*tπr*^{2}*c*(*z*,*t*). Similarly, the mass that flows out of the region through the upper boundary is approximately *v*(*t*)Δ*tπr*^{2}*c*(*z* + Δ*z*,*t*), hence the net change in mass due to advection (flow) is approximately −*v*(*t*)[*c*(*z* + Δ*z*,*t*) − *c*(*z*,*t*)]*πr*^{2}Δ*t*

#### The airway wall as a source of NO

Depending on the concentration, the tissue lining the airway wall can act as either a source or sink for NO. Typically, we assume the concentration of NO in the airway wall *c*_{w} exceeds the alveolar concentration *c*(*z*_{alv},*t*), implying the airway tissue serves as a net source of NO. The net transfer of NO from tissue to lumen is assumed to occur at a rate proportional to the concentration difference. Denoting the proportionality constant by *p*, the flux per unit area is *p*[*c*_{w} − *c*(*z*,*t*)].Therefore, over a time interval of length Δ*t* the net mass diffusing into the region from the airway wall is approximately equal to the product *p*[*c*_{w} − *c*(*z*,*t*)]2*πr*Δ*z*Δ*t*.

#### Diffusion of NO axially through the airway

Assuming *c*_{w} > *c*(*z*_{alv},*t*) implies that during exhalation the concentration will increase as air moves though the airway, i.e. *c*(*z* + Δ*z*,*t*) > *c*(*z*,*t*). According to Fick's first law, the diffusive flux per unit area will be proportional to the concentration gradient *c*_{z}(*z*,*t*) := (∂/∂*z*)*c*(*z*,*t*). Denoting by *d* the proportionality constant, over a time interval of length Δ*t* the mass diffusing into the region will be approximately *dc*_{z}(*z* + Δ*z*)*πr*^{2}Δ*t*. Similarly, the mass diffusing out of the region will be approximately *dc*_{z}(*z*,*t*)*πr*^{2}Δ*t*, hence the net change in mass due to (axial) diffusion is approximately *d*[*c*_{z}(*z* + Δ*z*,*t*) − *c*_{z}(*z*,*t*)]*πr*^{2}Δ*t*.

Imposing conservation of mass requires that over the time interval (*t*,*t* + Δ*t*), the net change inside the region must equal the net mass passing through the boundary. This implies the approximate equality [*c*(*z*,*t* + Δ*t*) − *c*(*z*,*t*)]*πr*^{2}Δ*z* ≈ − *v*(*t*)[*c*(*z* + Δ*z*,*t*) − *c*(*z*,*t*)]*πr*^{2}Δ*t* + *d*[*c*_{z}(*z* + Δ*z*,*t*) − *c*_{z}(*z*,*t*)]*πr*^{2}Δ*t* + *p*[*c*_{w} − *c*(*z*,*t*)]2*πr*Δ*z*Δ*t*. Dividing both sides by *πr*^{2}Δ*t*Δ*z* then taking Δ*t*,Δ*z*→0 results in the model (eq. 1).

#### eq. 2 as a special case of eq. 1

To transform eq. 1 into eq. 2, we begin by assuming the airway NO concentration has achieved a steady‐state, and thus does not vary through time. Mathematically, this is equivalent to assuming (∂/∂*t*)*c*(*z*,*t*) = 0, which also implies *c*(*z*,*t*) = *c*(*z*). Because the concentration is flow dependent, a necessary condition for this to occur is for the velocity to be constant through time, that is *v*(*t*) = *v*. As a further simplification, the diffusivity of NO in air is assumed to be negligible, i.e. *d* = 0, reducing the problem further to a first order ordinary differential equation (ODE).

By substituting (∂/∂*t*)*c*(*z*,*t*) = 0, *v*(*t*) = *v*, and *d* = 0 into equation eq. 1 it simplifies into 0 = −*vc*^{′}(*z*) + (2*p*/*r*)*c*_{w}−(2*p*/*r*)*c*(*z*). Multiplying through by the airway volume *πr*^{2}*l*, dividing through by , and rearranging terms yields . Generically, this is a linear ODE with constant coefficients of the form *c*^{′}(*z*) + [*a*/*d*]*c*(*z*) = *b*/*d*. Using standard techniques for ODEs, such as separation of variables, the general solution can be shown to be *c*(*z*) = *b*/*a* + *k* exp (−*z*[*a*/*d*]), where *k* is the constant of integration. Therefore, by making the above simplifying assumptions, the general solution of (eq. 1) can be written as .

The constant *k* can be found by imposing a boundary condition, and for consistency with Tsoukias and George (1998) the origin (*z*_{0}) is shifted from the sensor to the alveolar boundary. That is, by definition the concentration at the origin is equal to the alveolar concentration, or . Combining this condition with the general solution above allows us to solve for the constant, resulting in , and hence the particular solution

This solution is applicable along the length of the airway from the alveolar boundary to the mouth, or 0 ≤ *z* ≤*l*. Once air enters the sampling device dead space, both the airway wall concentration and permeability are assumed to be zero. This is equivalent to assuming *c*^{′}(*z*) = 0 in this region, i.e. there is no spatial variation in the concentration. This implies the concentration measured at the sensor equals the concentration measured at the mouth. Mathematically, this means , and inserting *z* = *l* and into the particular solution yields the equality (eq. 2).

## Appendix B Numerical integration

#### Spatial discretizations

To solve the governing PDE (eq. 1) numerically, the spatial (*z*) derivatives are replaced with finite difference approximations based on Taylor series expansions. For the diffusive term, a centered three term Taylor series approximation is employed, (∂^{2}/∂*z*^{2})*c*(*z*,*t*) ≈ [*c*(*z* − Δ*z*,*t*) − 2*c*(*z*,*t*) + *c*(*z* + Δ*z*,*t*)][Δ*z*]^{−2}. For the advective term, a biased 4 term Taylor series approximation is employed. The direction of the bias is determined by the direction of flow, as dictated by the sign of *v*(*t*). Specifically, the approximation is oriented with an “upwind” bias; two of the terms in the approximation are chosen on the side from which the flow originates, and only one is chosen from the opposite side.

For example, when *v*(*t*) > 0 the approximation is (∂/∂*z*)*c*(*z*,*t*) ≈ [*c*(*z* − 2Δ*z*,*t*) − 6*c*(*z* − Δ*z*,*t*) + 3*c*(*z*,*t*) + 2*c*(*z* + Δ*z*,*t*)][6Δ*z*]^{−1}. An upwind discretization is employed because centered discretizations for advection can lead to spurious oscillations in the numerical approximation. A completely one‐sided discretization can prevent oscillations; however, despite the formal order of the Taylor series, such an approximation will always have first order accuracy (Hundsdorfer and Verwer 2003). By employing a two‐sided, but biased, discretization, higher order accuracy can be achieved while minimizing the potential for oscillatory solutions.

#### Boundary and initial conditions

For the PDE (eq. 1) to have a unique solution, initial (time) and inflow boundary (space) conditions must be specified. Moreover, the incoming concentration depends on the direction of flow. During exhalation it corresponds to the alveolar concentration, a parameter to be estimated. During inhalation this concentration is typically the ambient NO level; however, during testing subjects may be provided air that has been “scrubbed” of NO.

By definition, modeling involves modeling exhalation. However, because respiration is cyclic, the terminal condition in one direction of flow becomes the initial condition for the reverse flow. This relationship means that NO measured during exhalation is determined, in part, by the terminal state of the previous inhalation. In principle, the previous inhalation depends, in turn, on the preceding exhalation, which depends on the inhalation before that, continuing ad nauseam.

In practice, higher flow rates diminish this dependence, and at relatively high rates (300+ mL/sec), the terminal airway concentration is effectively independent of the initial. Although 300 mL/sec is a relatively rapid rate for exhalation, it is a relatively slow rate for inhalation. For example, in all 900 maneuvers used for simulation this threshold was cleared every time, typically by factors of at least 2‐3x.

The implication of this phenomena is that calculating accurate estimates of the airway concentration immediately after inhalation does not require knowing the initial airway concentration when inhalation began. The solution can be calculated based on a simple initial condition (i.e. zero everywhere), and the end result will be essentially unchanged. The terminal condition will also depend on the inflow concentration; however, in the case of “scrubbed” air this can be assumed to be zero.

## Appendix C MCMC parameter estimation

We adopt a Bayesian approach to inference (Gelman et al. 2004), and our goal is to characterize the posterior distribution (generically denoted *f*(*θ*|*x*) for parameters *θ* and data *x*). Using Bayes rule, the posterior can be expressed in terms of the likelihood *f*(*x*|*θ*) and a prior distribution *f*(*θ*). For many purposes, including ours, the unnormalized posterior is sufficient, simplifying the relationship between to the posterior, likelihood, and prior to the proportionality *f*(*θ*|*x*) ∝ *f*(*x*|*θ*)*f*(*θ*).

Multiple parameterizations of the model (eq. 1) are possible; for consistency with eq. 2, we work with the vector of parameters . The data for each subject consists of , where is the measured concentration at time *t*_{i} during maneuver *j* ∈ 1,2, …, 9.

To formulate a likelihood, we assume that by fixing *θ* and solving the corresponding model equation, the model solution can be used to calculate the density of the observed data. If we further assume the observed values arise from a shared parametric conditional distribution with density function *f*, and that conditional on the model solutions *c*_{ij} the are independent, then the likelihood can be written as , where *θ* appears implicitly on the right via the model solution *c*_{ij}.

When the likelihood *f*(*x*|*θ*) is combined with a prior *f*(*θ*), the (unnormalized) posterior can be easily calculated for any particular set of parameters *θ*. To efficiently explore the posterior distribution we employ a Metropolis‐Hastings style MCMC algorithm (Robert and Casella, 2005), which generically proceeds as follows:

0. Select an initial value

*θ*and calculate the likelihood*f*(*θ*|*x*).1. Propose a new value

*θ*^{′}using a transition kernel*q*(*θ*→*θ*^{′}), and calculate the likelihood .2. Accept the proposed value with probability .

3. If the proposal is accepted set

*θ*=*θ*^{′}, then return to 1; otherwise, return to 1.

The choice of transition kernel *q* can have a significant impact on the efficiency of this type of algorithm. Finding an optimal *q* can be difficult; however, there are a number of more recent MCMC algorithms which incorporate an “adaptive” transition distribution (Roberts and Rosenthal 2009). To better account for variability in the posterior across individuals, the adaptive Metropolis algorithm of Haario et al. (2001) is employed to automatically calibrate the proposal distribution against the target. This has the dual benefit of both increasing the efficiency of our sampler, while also simplifying the user experience by largely automating the choice of transition kernel.