Bayesian estimation of physiological parameters governing a dynamic two‐compartment model of exhaled nitric oxide

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€ ogman 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 Pietropaoli et al. 1999;H€ ogman 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 .
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 ,ĉ i is a numerical approximation of c i , andc i 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 (secondsminutes) time scales.
• Daw NO = 2prlp. 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 2p/r in equation 1.
• J 0 aw NO = 2prlpc 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 2pc w /r in equation 1.

Methods
Our proposed airway model is a variant of the two-compartment approach , the primary distinction being we make no assumptions regarding the flow rate. In its simplest form, the twocompartment 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 ; 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 . 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): @ @t cðz; tÞ ¼ À vðtÞ @ @z cðz; tÞ þ d @ 2 @z 2 cðz; tÞ þ 2p r c w À cðz; tÞ ½ 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 solutionŝ c 0 ;ĉ 1 ; . . .;ĉ n , 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 During exhalation, air is expelled from the alveolar compartment, entering the airway compartment through the right boundary of the model cylinder. To account for dead space in the sampling instrument, the airway cylinder is extended beyond the mouth by the corresponding volume. 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 functionvðtÞ 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, d Fe NO =ĉ i . In addition to the estimated flow rate function vðtÞ, the approximate solutions also depend on the airway parameter values. In these examples, CA NO = 2, J 0 aw NO =800, and Daw NO = 5 (values identical to those used in a previous simulation study [Citation Eckel et al. 2014]).
CombiningvðtÞ with the parameters CA NO , J 0 aw NO , and Daw NO enables the use of numerical integration to calculate the sequence of approximationsĉ i . 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ĉ i is calculated at 100s-1000s of time steps t i , resulting in a squence of approximations fĉ i g n i¼0 . 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 functionsvðtÞ 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 0 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 valuẽ c i . To account for this residual variation, we assume lnðc i Þ is normally distributed with mean ln (c i ) and variance r 2 = 0.1 2 , that isc i has a log-normal distribution. This is the same subject shown in Figure 3 (and  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 0 aw NO , and Daw NO , are based on the steady-state model (George et al. 2004) relating the parameters of interest to the exhaled concentration ( _ V is the volumetric flow rate, which is equivalent to the product v(t)pr 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 = ( . . . )+e (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 (. . .)+e (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)  . 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€ ogman and Meril€ ainen (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   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 0 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  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 0 aw NO and Daw NO generated by the HMA and non-LinLog 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 ð xÞ and standard error ðseð xÞÞ 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 x, 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 (6 ¼) 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 0 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 0 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 0 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) Table 1. For each combination of model parameter and estimation method, in the first two columns we report the sample mean and standard error of the 100 simulation study estimates. The third column reports the P-value resulting from a test of equality between the mean of the sampling distribution and the true value used during simulation. For these simulated subjects, the dynamic approach is the only one we conclude to be unbiased for CA NO , and for all three parameters it provides consistently greater precision than any other method. 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. 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 flowadjusted 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., _ V 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 Dz = 0.2 cm is sufficient (Shin and George 2002); however, our simulations indicate this is Table 2. Parameter estimates based on real FE NO input data. Interval estimates based on frequentist (conf) or Bayesian (cred) approaches have also been calculated when applicable. The estimated parameter values are significantly larger than the simulation study inputs, although the rank ordering of the estimates largely agrees with the rank ordering in Table 1 inadequate. Specifically, the simulated profiles illustrated in Figure 4 are based on using 1000 grid points (corresponding to a spatial resolution Dz=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 Dz ¼ 0:08 3 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 0 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. flows into the region through the lower boundary is approximately v(t)Dtpr 2 c(z,t). Similarly, the mass that flows out of the region through the upper boundary is approximately v(t)Dtpr 2 c(z + Dz,t), hence the net change in mass due to advection (flow) is approximately Àv(t)[c (z + Dz,t) À c(z,t)]pr 2 Dt 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 Dt the net mass diffusing into the region from the airway wall is approximately equal to the product p[c w À c(z,t)]2prDzDt.

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 + Dz,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 Dt the mass diffusing into the region will be approximately dc z (z + Dz)pr 2 Dt. Similarly, the mass diffusing out of the region will be approximately dc z (z,t) pr 2 Dt, hence the net change in mass due to (axial) diffusion is approximately d[c z (z + Dz,t) À c z (z,t)]pr 2 Dt.
Imposing conservation of mass requires that over the time interval (t,t + Dt), the net change inside the region must equal the net mass passing through the boundary. This implies the approximate equality [c(z,t + Dt) À c(z,t)] pr 2 Dz % À v(t)[c(z + Dz,t) À c(z,t)]pr 2 Dt + d[c z (z + Dz,t) À c z (z,t)]pr 2 Dt + p[c w À c(z,t)]2prDzDt. Dividing both sides by pr 2 DtDz then taking Dt,Dz?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 0 (z) + (2p/r) c w À(2p/r)c(z). Multiplying through by the airway volume pr 2 l, dividing through by _ Vl, and rearranging terms yields c 0 ðzÞ þ ½Daw NO =ð _ VlÞcðzÞ ¼ J 0 aw NO =ð _ VlÞ. Generically, this is a linear ODE with constant coefficients of the form c 0 (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 Figure A1. For an arbitrary airway segment of length Dz, and an arbitrary time interval of length Dt, we assume any net change in mass inside the segment must equal the net mass passing through the boundary. As illustrated here, air moves vertically through the airway during exhalation; therefore, the net change due to advection (flow) represents the difference between the mass entering through the lower boundary and exiting through the upper boundary. The airway wall is also assumed to be permeable to NO. Depending on the relative concentrations of the airway wall c w and the alveolar compartment c(z alv ,t), it will serve as either a net source or sink of NO. Typically, c w > c(z,t), yielding a net diffusion of NO from the airway wall to lumen (as pictured here). This also implies that the concentration in the lumen will increase as air moves vertically though the airway, that is, c(z + Dz,t) > c(z,t). This concentration gradient results in net NO diffusion opposite to the direction of flow; therefore, the net change in mass due to diffusion will be the difference between the mass diffusing in through the upper boundary and the mass diffusing out through the lower boundary.